# source:anuga_core/source/anuga/culvert_flows/culvert_class.py@5435

Last change on this file since 5435 was 5435, checked in by ole, 15 years ago

Moved culvert class code from development into ANUGA

File size: 14.9 KB
Line
1from anuga.shallow_water.shallow_water_domain import Inflow, General_forcing
2from anuga.culvert_flows.culvert_polygons import create_culvert_polygons
3from anuga.utilities.system_tools import log_to_file
4
5class Culvert_flow:
6    """Culvert flow - transfer water from one hole to another
7
8    Using Momentum as Calculated by Culvert Flow !!
9    Could be Several Methods Investigated to do This !!!
10
11    2008_May_08
12    To Ole:
13    OK so here we need to get the Polygon Creating code to create polygons for the culvert Based on
14    the two input Points (X0,Y0) and (X1,Y1) - need to be passed to create polygon
15
16    The two centers are now passed on to create_polygon.
17
18
19    Input: Two points, pipe_size (either diameter or width, height), mannings_rougness,
20    inlet/outlet energy_loss_coefficients, internal_bend_coefficent,
21    top-down_blockage_factor and bottom_up_blockage_factor
22
23
24    And the Delta H enquiry should be change from Openings in line 412 to the enquiry Polygons infront
25    of the culvert
26    At the moment this script uses only Depth, later we can change it to Energy...
27
28    Once we have Delta H can calculate a Flow Rate and from Flow Rate an Outlet Velocity
29    The Outlet Velocity x Outlet Depth = Momentum to be applied at the Outlet...
30
31    """
32
33    def __init__(self,
34                 domain,
35                 label=None,
36                 description=None,
37                 end_point0=None,
38                 end_point1=None,
39                 width=None,
40                 height=None,
41                 diameter=None,
42                 manning=None,          # Mannings Roughness for Culvert
43                 invert_level0=None,    # Invert level if not the same as the Elevation on the Domain
44                 invert_level1=None,    # Invert level if not the same as the Elevation on the Domain
45                 loss_exit=None,
46                 loss_entry=None,
47                 loss_bend=None,
48                 loss_special=None,
49                 blockage_topdwn=None,
50                 blockage_bottup=None,
51                 culvert_routine=None,
52                 verbose=False):
53
54        from Numeric import sqrt, sum
55
56        # Input check
57        if diameter is not None:
58            self.culvert_type = 'circle'
59            self.diameter = diameter
60            if height is not None or width is not None:
61                msg = 'Either diameter or width&height must be specified, but not both.'
62                raise Exception, msg
63        else:
64            if height is not None:
65                if width is None:
66                    self.culvert_type = 'square'
67                    width = height
68                else:
69                    self.culvert_type = 'rectangle'
70            elif width is not None:
71                if height is None:
72                    self.culvert_type = 'square'
73                    height = width
74            else:
75                msg = 'Either diameter or width&height must be specified.'
76                raise Exception, msg
77
78            if height == width:
79                self.culvert_type = 'square'
80
81            self.height = height
82            self.width = width
83
84
85        assert self.culvert_type in ['circle', 'square', 'rectangle']
86
87        # Set defaults
88        if manning is None: manning = 0.012   # Set a Default Mannings Roughness for Pipe
89        if loss_exit is None: loss_exit = 1.00
90        if loss_entry is None: loss_entry = 0.50
91        if loss_bend is None: loss_bend=0.00
92        if loss_special is None: loss_special=0.00
93        if blockage_topdwn is None: blockage_topdwn=0.00
94        if blockage_bottup is None: blockage_bottup=0.00
95        if culvert_routine is None: culvert_routine=boyd_generalised_culvert_model
96        if label is None: label = 'culvert_flow_' + id(self)
97
98        # Open log file for storing some specific results...
99        self.log_filename = label + '.log'
100        self.label = label
101
102
103        # Create the fundamental culvert polygons from POLYGON
104        P = create_culvert_polygons(end_point0,
105                                    end_point1,
106                                    width=width,
107                                    height=height)
108
109        if verbose is True:
110            pass
111            #plot_polygons([[end_point0, end_point1],
112            #               P['exchange_polygon0'],
113            #               P['exchange_polygon1'],
114            #               P['enquiry_polygon0'],
115            #               P['enquiry_polygon1']],
116            #              figname='culvert_polygon_output')
117
118        self.openings = []
119        self.openings.append(Inflow(domain,
120                                    polygon=P['exchange_polygon0']))
121
122        self.openings.append(Inflow(domain,
123                                    polygon=P['exchange_polygon1']))
124
125
126        # Assume two openings for now: Referred to as 0 and 1
127        assert len(self.openings) == 2
128
129        # Store basic geometry
130        self.end_points = [end_point0, end_point1]
131        self.invert_levels = [invert_level0, invert_level1]
132        self.enquiry_polygons = [P['enquiry_polygon0'], P['enquiry_polygon1']]
133        self.vector = P['vector']
134        self.length = P['length']; assert self.length > 0.0
135        self.verbose = verbose
136        self.last_time = 0.0
137        self.temp_keep_delta_t = 0.0
138
139
140        # Store hydraulic parameters
141        self.manning = manning
142        self.loss_exit = loss_exit
143        self.loss_entry = loss_entry
144        self.loss_bend = loss_bend
145        self.loss_special = loss_special
146        self.sum_loss = loss_exit + loss_entry + loss_bend + loss_special
147        self.blockage_topdwn = blockage_topdwn
148        self.blockage_bottup = blockage_bottup
149
150        # Store culvert routine
151        self.culvert_routine = culvert_routine
152
153
154        # Create objects to update momentum (a bit crude at this stage)
155
156
157        xmom0 = General_forcing(domain, 'xmomentum',
158                                polygon=P['exchange_polygon0'])
159
160        xmom1 = General_forcing(domain, 'xmomentum',
161                                polygon=P['exchange_polygon1'])
162
163        ymom0 = General_forcing(domain, 'ymomentum',
164                                polygon=P['exchange_polygon0'])
165
166        ymom1 = General_forcing(domain, 'ymomentum',
167                                polygon=P['exchange_polygon1'])
168
169        self.opening_momentum = [ [xmom0, ymom0], [xmom1, ymom1] ]
170
171
172        # Print something
173        s = 'Culvert Effective Length = %.2f m' %(self.length)
174        log_to_file(self.log_filename, s)
175
176        s = 'Culvert Direction is %s\n' %str(self.vector)
177        log_to_file(self.log_filename, s)
178
179    def __call__(self, domain):
180        from anuga.utilities.numerical_tools import mean
181        from anuga.utilities.polygon import inside_polygon
182        from anuga.config import g, epsilon
183        from Numeric import take, sqrt
184        from anuga.config import velocity_protection
185
186
187        log_filename = self.log_filename
188
189        # Time stuff
190        time = domain.get_time()
191        delta_t = time-self.last_time
192        s = '\nTime = %.2f, delta_t = %f' %(time, delta_t)
193        log_to_file(log_filename, s)
194
195        msg = 'Time did not advance'
196        if time > 0.0: assert delta_t > 0.0, msg
197
198
199        # Get average water depths at each opening
200        openings = self.openings   # There are two Opening  and 
201        for i, opening in enumerate(openings):
202            stage = domain.quantities['stage'].get_values(location='centroids',
203                                                          indices=opening.exchange_indices)
204            elevation = domain.quantities['elevation'].get_values(location='centroids',
205                                                                  indices=opening.exchange_indices)
206
207            # Indices corresponding to energy enquiry field for this opening
208            coordinates = domain.get_centroid_coordinates() # Get all centroid points (x,y)
209            enquiry_indices = inside_polygon(coordinates, self.enquiry_polygons[i])
210
211
212            # Get model values for points in enquiry polygon for this opening
213            dq = domain.quantities
214            stage = dq['stage'].get_values(location='centroids', indices=enquiry_indices)
215            xmomentum = dq['xmomentum'].get_values(location='centroids', indices=enquiry_indices)
216            ymomentum = dq['ymomentum'].get_values(location='centroids', indices=enquiry_indices)
217            elevation = dq['elevation'].get_values(location='centroids', indices=enquiry_indices)
218            depth = stage - elevation
219
220            # Compute mean values of selected quantitites in the enquiry area in front of the culvert
221            # Epsilon handles a dry cell case
222            ux = xmomentum/(depth+velocity_protection/depth)   # Velocity (x-direction)
223            uy = ymomentum/(depth+velocity_protection/depth)   # Velocity (y-direction)
224            v = mean(sqrt(ux**2+uy**2))      # Average velocity
225            w = mean(stage)                  # Average stage
226
227            # Store values at enquiry field
228            opening.velocity = v
229
230
231            # Compute mean values of selected quantitites in the exchange area in front of the culvert
232            # Stage and velocity comes from enquiry area and elevation from exchange area
233
234            # Use invert level instead of elevation if specified
235            invert_level = self.invert_levels[i]
236            if invert_level is not None:
237                z = invert_level
238            else:
239                elevation = dq['elevation'].get_values(location='centroids', indices=opening.exchange_indices)
240                z = mean(elevation)                   # Average elevation
241
242            # Estimated depth above the culvert inlet
243            d = w - z
244
245            if d < 0.0:
246                # This is possible since w and z are taken at different locations
247                #msg = 'D < 0.0: %f' %d
248                #raise msg
249                d = 0.0
250
251            # Ratio of depth to culvert height.
252            # If ratio > 1 then culvert is running full
253            ratio = d/self.height
254            opening.ratio = ratio
255
256            # Average measures of energy in front of this opening
257            Es = d + 0.5*v**2/#  Specific energy in exchange area
258            Et = w + 0.5*v**2/#  Total energy in the enquiry area
259            opening.total_energy = Et
260            opening.specific_energy = Es
261
262            # Store current average stage and depth with each opening object
263            opening.depth = d
264            opening.stage = w
265            opening.elevation = z
266
267
268        #################  End of the FOR loop ################################################
269
270
271        # We now need to deal with each opening individually
272
273        # Determine flow direction based on total energy difference
274        delta_Et = openings.total_energy - openings.total_energy
275
276        if delta_Et > 0:
277            #print 'Flow U/S ---> D/S'
278            inlet=openings
279            outlet=openings
280
281            inlet.momentum = self.opening_momentum
282            outlet.momentum = self.opening_momentum
283
284        else:
285            #print 'Flow D/S ---> U/S'
286            inlet=openings
287            outlet=openings
288
289            inlet.momentum = self.opening_momentum
290            outlet.momentum = self.opening_momentum
291
292            delta_Et = -delta_Et
293
294        msg = 'Total energy difference is negative'
295        assert delta_Et > 0.0, msg
296
297        delta_h = inlet.stage - outlet.stage
298        delta_z = inlet.elevation - outlet.elevation
299        culvert_slope = (delta_z/self.length)
300
301        if culvert_slope < 0.0:
303            # Flow will be purely controlled by uphill outlet face
304            print 'WARNING: Flow is running uphill. Watch Out!', inlet.elevation, outlet.elevation
305
306
307        s = 'Delta total energy = %.3f' %(delta_Et)
308        log_to_file(log_filename, s)
309
310
311        Q, barrel_velocity, culvert_outlet_depth = self.culvert_routine(self, inlet, outlet, delta_Et, g)
312        #####################################################
313        barrel_momentum = barrel_velocity*culvert_outlet_depth
314
315        s = 'Barrel velocity = %f' %barrel_velocity
316        log_to_file(log_filename, s)
317
318        # Compute momentum vector at outlet
319        outlet_mom_x, outlet_mom_y = self.vector * barrel_momentum
320
321        s = 'Directional momentum = (%f, %f)' %(outlet_mom_x, outlet_mom_y)
322        log_to_file(log_filename, s)
323
324        delta_t = time - self.last_time
325        if delta_t > 0.0:
326            xmomentum_rate = outlet_mom_x - outlet.momentum.value
327            xmomentum_rate /= delta_t
328
329            ymomentum_rate = outlet_mom_y - outlet.momentum.value
330            ymomentum_rate /= delta_t
331
332            s = 'X Y MOM_RATE = (%f, %f) ' %(xmomentum_rate, ymomentum_rate)
333            log_to_file(log_filename, s)
334        else:
335            xmomentum_rate = ymomentum_rate = 0.0
336
337
338        # Set momentum rates for outlet jet
339        outlet.momentum.rate = xmomentum_rate
340        outlet.momentum.rate = ymomentum_rate
341
342        # Remember this value for next step (IMPORTANT)
343        outlet.momentum.value = outlet_mom_x
344        outlet.momentum.value = outlet_mom_y
345
346        if int(domain.time*100) % 100 == 0:
347            s = 'T=%.5f, Culvert Discharge = %.3f f'\
348                %(time, Q)
349            s += ' Depth= %0.3f  Momentum = (%0.3f, %0.3f)'\
350                 %(culvert_outlet_depth, outlet_mom_x,outlet_mom_y)
351            s += ' Momentum rate: (%.4f, %.4f)'\
352                 %(xmomentum_rate, ymomentum_rate)
353            s+='Outlet Vel= %.3f'\
354                %(barrel_velocity)
355            log_to_file(log_filename, s)
356
357
358
359
360
361        # Execute flow term for each opening
362        # This is where Inflow objects are evaluated and update the domain
363        for opening in self.openings:
364            opening(domain)
365
366        # Execute momentum terms
367        # This is where Inflow objects are evaluated and update the domain
368        outlet.momentum(domain)
369        outlet.momentum(domain)
370
371        # Store value of time
372        self.last_time = time
373
374
Note: See TracBrowser for help on using the repository browser.