source: anuga_core/source/anuga/culvert_flows/culvert_flow.py @ 5434

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

Added culvert framework

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