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

Last change on this file since 5586 was 5586, checked in by ole, 14 years ago

Rating curves for culverts

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