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

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

Finished culvert flow work with a good working solution

File size: 15.8 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'
97        label += '_' + str(id(self)) 
98       
99        # Open log file for storing some specific results...
100        self.log_filename = label + '.log' 
101        self.label = label
102
103        # Print something
104        log_to_file(self.log_filename, self.label)       
105        log_to_file(self.log_filename, description)
106        log_to_file(self.log_filename, self.culvert_type)       
107
108
109        # Create the fundamental culvert polygons from POLYGON
110        if self.culvert_type == 'circle':
111            # Redefine width and height for use with create_culvert_polygons
112            width = height = diameter
113       
114        P = create_culvert_polygons(end_point0,
115                                    end_point1,
116                                    width=width,   
117                                    height=height)
118       
119        if verbose is True:
120            pass
121            #plot_polygons([[end_point0, end_point1],
122            #               P['exchange_polygon0'],
123            #               P['exchange_polygon1'],
124            #               P['enquiry_polygon0'],
125            #               P['enquiry_polygon1']],
126            #              figname='culvert_polygon_output')
127       
128        self.openings = []
129        self.openings.append(Inflow(domain,
130                                    polygon=P['exchange_polygon0']))
131
132        self.openings.append(Inflow(domain,
133                                    polygon=P['exchange_polygon1']))                                   
134
135
136        # Assume two openings for now: Referred to as 0 and 1
137        assert len(self.openings) == 2
138       
139        # Store basic geometry
140        self.end_points = [end_point0, end_point1]
141        self.invert_levels = [invert_level0, invert_level1]               
142        self.enquiry_polygons = [P['enquiry_polygon0'], P['enquiry_polygon1']]
143        self.vector = P['vector']
144        self.length = P['length']; assert self.length > 0.0
145        self.verbose = verbose
146        self.last_time = 0.0       
147       
148
149        # Store hydraulic parameters
150        self.manning = manning
151        self.loss_exit = loss_exit
152        self.loss_entry = loss_entry
153        self.loss_bend = loss_bend
154        self.loss_special = loss_special
155        self.sum_loss = loss_exit + loss_entry + loss_bend + loss_special
156        self.blockage_topdwn = blockage_topdwn
157        self.blockage_bottup = blockage_bottup
158
159        # Store culvert routine
160        self.culvert_routine = culvert_routine
161
162       
163        # Create objects to update momentum (a bit crude at this stage)
164
165       
166        xmom0 = General_forcing(domain, 'xmomentum',
167                                polygon=P['exchange_polygon0'])
168
169        xmom1 = General_forcing(domain, 'xmomentum',
170                                polygon=P['exchange_polygon1'])
171
172        ymom0 = General_forcing(domain, 'ymomentum',
173                                polygon=P['exchange_polygon0'])
174
175        ymom1 = General_forcing(domain, 'ymomentum',
176                                polygon=P['exchange_polygon1'])
177
178        self.opening_momentum = [ [xmom0, ymom0], [xmom1, ymom1] ]
179       
180
181        # Print something
182        s = 'Culvert Effective Length = %.2f m' %(self.length)
183        log_to_file(self.log_filename, s)
184   
185        s = 'Culvert Direction is %s\n' %str(self.vector)
186        log_to_file(self.log_filename, s)       
187       
188    def __call__(self, domain):
189        from anuga.utilities.numerical_tools import mean
190        from anuga.utilities.polygon import inside_polygon
191        from anuga.config import g, epsilon
192        from Numeric import take, sqrt
193        from anuga.config import velocity_protection       
194
195
196        log_filename = self.log_filename
197         
198        # Time stuff
199        time = domain.get_time()
200        delta_t = time-self.last_time
201        s = '\nTime = %.2f, delta_t = %f' %(time, delta_t)
202        log_to_file(log_filename, s)
203       
204        msg = 'Time did not advance'
205        if time > 0.0: assert delta_t > 0.0, msg
206
207
208        # Get average water depths at each opening       
209        openings = self.openings   # There are two Opening [0] and [1]
210        for i, opening in enumerate(openings):
211            stage = domain.quantities['stage'].get_values(location='centroids',
212                                                          indices=opening.exchange_indices)
213            elevation = domain.quantities['elevation'].get_values(location='centroids',
214                                                                  indices=opening.exchange_indices)
215
216            # Indices corresponding to energy enquiry field for this opening
217            coordinates = domain.get_centroid_coordinates() # Get all centroid points (x,y)
218            enquiry_indices = inside_polygon(coordinates, self.enquiry_polygons[i]) 
219           
220           
221            # Get model values for points in enquiry polygon for this opening
222            dq = domain.quantities
223            stage = dq['stage'].get_values(location='centroids', indices=enquiry_indices)
224            xmomentum = dq['xmomentum'].get_values(location='centroids', indices=enquiry_indices)
225            ymomentum = dq['ymomentum'].get_values(location='centroids', indices=enquiry_indices)                       
226            elevation = dq['elevation'].get_values(location='centroids', indices=enquiry_indices)
227            depth = stage - elevation
228
229            # Compute mean values of selected quantitites in the enquiry area in front of the culvert
230            # Epsilon handles a dry cell case
231            ux = xmomentum/(depth+velocity_protection/depth)   # Velocity (x-direction)
232            uy = ymomentum/(depth+velocity_protection/depth)   # Velocity (y-direction)
233            v = mean(sqrt(ux**2+uy**2))      # Average velocity
234            w = mean(stage)                  # Average stage
235
236            # Store values at enquiry field
237            opening.velocity = v
238
239
240            # Compute mean values of selected quantitites in the exchange area in front of the culvert
241            # Stage and velocity comes from enquiry area and elevation from exchange area
242           
243            # Use invert level instead of elevation if specified
244            invert_level = self.invert_levels[i]
245            if invert_level is not None:
246                z = invert_level
247            else:
248                elevation = dq['elevation'].get_values(location='centroids', indices=opening.exchange_indices)
249                z = mean(elevation)                   # Average elevation
250
251            # Estimated depth above the culvert inlet
252            d = w - z  # Used for calculations involving depth
253            if d < 0.0:
254                # This is possible since w and z are taken at different locations
255                #msg = 'D < 0.0: %f' %d
256                #raise msg
257                d = 0.0
258           
259
260
261            # Depth at exchange area used to trigger calculations
262            stage = dq['stage'].get_values(location='centroids', indices=enquiry_indices)
263            elevation = dq['elevation'].get_values(location='centroids', indices=enquiry_indices)
264            depth = stage - elevation
265            d_trigger = mean(depth)
266
267
268
269            # Ratio of depth to culvert height.
270            # If ratio > 1 then culvert is running full
271            if self.culvert_type == 'circle':
272                ratio = d/self.diameter
273            else:   
274                ratio = d/self.height 
275            opening.ratio = ratio
276               
277            # Average measures of energy in front of this opening
278            Es = d + 0.5*v**2/#  Specific energy in exchange area
279            Et = w + 0.5*v**2/#  Total energy in the enquiry area
280            opening.total_energy = Et
281            opening.specific_energy = Es           
282           
283            # Store current average stage and depth with each opening object
284            opening.depth = d
285            opening.depth_trigger = d_trigger           
286            opening.stage = w
287            opening.elevation = z
288           
289
290        #################  End of the FOR loop ################################################
291
292           
293        # We now need to deal with each opening individually
294
295        # Determine flow direction based on total energy difference
296        delta_Et = openings[0].total_energy - openings[1].total_energy
297
298        if delta_Et > 0:
299            #print 'Flow U/S ---> D/S'
300            inlet=openings[0]
301            outlet=openings[1]
302
303            inlet.momentum = self.opening_momentum[0]
304            outlet.momentum = self.opening_momentum[1]
305
306        else:
307            #print 'Flow D/S ---> U/S'
308            inlet=openings[1]
309            outlet=openings[0]
310
311            inlet.momentum = self.opening_momentum[1]
312            outlet.momentum = self.opening_momentum[0]
313           
314            delta_Et = -delta_Et
315
316        msg = 'Total energy difference is negative'
317        assert delta_Et > 0.0, msg
318
319        delta_h = inlet.stage - outlet.stage
320        delta_z = inlet.elevation - outlet.elevation
321        culvert_slope = (delta_z/self.length)
322
323        if culvert_slope < 0.0:
324            # Adverse gradient - flow is running uphill
325            # Flow will be purely controlled by uphill outlet face
326            print 'WARNING: Flow is running uphill. Watch Out!', inlet.elevation, outlet.elevation
327
328
329        s = 'Delta total energy = %.3f' %(delta_Et)
330        log_to_file(log_filename, s)
331
332       
333        Q, barrel_velocity, culvert_outlet_depth = self.culvert_routine(self, inlet, outlet, delta_Et, g)
334        #####################################################
335        barrel_momentum = barrel_velocity*culvert_outlet_depth
336               
337        s = 'Barrel velocity = %f' %barrel_velocity
338        log_to_file(log_filename, s)
339
340        # Compute momentum vector at outlet
341        outlet_mom_x, outlet_mom_y = self.vector * barrel_momentum
342           
343        s = 'Directional momentum = (%f, %f)' %(outlet_mom_x, outlet_mom_y)
344        log_to_file(log_filename, s)
345
346        delta_t = time - self.last_time
347        if delta_t > 0.0:
348            xmomentum_rate = outlet_mom_x - outlet.momentum[0].value
349            xmomentum_rate /= delta_t
350               
351            ymomentum_rate = outlet_mom_y - outlet.momentum[1].value
352            ymomentum_rate /= delta_t
353                   
354            s = 'X Y MOM_RATE = (%f, %f) ' %(xmomentum_rate, ymomentum_rate)
355            log_to_file(log_filename, s)                   
356        else:
357            xmomentum_rate = ymomentum_rate = 0.0
358
359
360        # Set momentum rates for outlet jet
361        outlet.momentum[0].rate = xmomentum_rate
362        outlet.momentum[1].rate = ymomentum_rate
363
364        # Remember this value for next step (IMPORTANT)
365        outlet.momentum[0].value = outlet_mom_x
366        outlet.momentum[1].value = outlet_mom_y                   
367
368        if int(domain.time*100) % 100 == 0:
369            s = 'T=%.5f, Culvert Discharge = %.3f f'\
370                %(time, Q)
371            s += ' Depth= %0.3f  Momentum = (%0.3f, %0.3f)'\
372                 %(culvert_outlet_depth, outlet_mom_x,outlet_mom_y)
373            s += ' Momentum rate: (%.4f, %.4f)'\
374                 %(xmomentum_rate, ymomentum_rate)                   
375            s+='Outlet Vel= %.3f'\
376                %(barrel_velocity)
377            log_to_file(log_filename, s)
378           
379               
380
381
382       
383        # Execute flow term for each opening
384        # This is where Inflow objects are evaluated and update the domain
385        for opening in self.openings:
386            opening(domain)
387
388        # Execute momentum terms
389        # This is where Inflow objects are evaluated and update the domain
390        outlet.momentum[0](domain)
391        outlet.momentum[1](domain)       
392           
393        # Store value of time   
394        self.last_time = time
395
396
Note: See TracBrowser for help on using the repository browser.