source: anuga_core/source/anuga/culvert_flows/culvert_class.py @ 5585

Last change on this file since 5585 was 5585, checked in by ole, 16 years ago

Multibarrel culvert functionality.

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