source: trunk/anuga_core/source/anuga/structures/culvert_routines.py @ 7978

Last change on this file since 7978 was 7978, checked in by habili, 13 years ago

Refactoring of the culvert code.

File size: 19.9 KB
Line 
1from anuga.config import velocity_protection
2from anuga.utilities.numerical_tools import safe_acos as acos
3
4from math import pi, sqrt, sin, cos
5from anuga.config import g
6
7class Culvert_routines:
8    """Collection of culvert routines for use with Culvert_operator
9
10    This module holds various routines to determine FLOW through CULVERTS and SIMPLE BRIDGES
11
12    Usage:
13
14    NOTE:
15     Inlet control:  self.delta_total_energy > self.inflow.get_average_specific_energy()
16     Outlet control: self.delta_total_energy < self.inflow.get_average_specific_energy()
17     where total energy is (w + 0.5*v^2/g) and
18     specific energy is (h + 0.5*v^2/g)
19    """
20
21    def __init__(self, culvert, manning=0.0):
22       
23        self.inlets = culvert.inlets
24       
25        self.inflow  = self.inlets[0]
26        self.outflow = self.inlets[1]
27       
28        self.culvert_length = culvert.get_culvert_length()
29        self.culvert_width = culvert.get_culvert_width()
30        self.culvert_height = culvert.get_culvert_height()
31        self.sum_loss = 0.0
32        self.max_velocity = 10.0
33        self.manning = manning
34        self.log_filename = None
35       
36        # Determine flow direction based on total energy difference
37        self.delta_total_energy = self.inflow.get_average_total_energy() - self.outflow.get_average_total_energy()
38
39        if self.delta_total_energy < 0:
40            self.inflow  = self.inlets[1]
41            self.outflow = self.inlets[0]
42            self.delta_total_energy = -self.delta_total_energy
43
44        #delta_z = self.self.inflow.get_average_elevation() - self.self.outflow.get_average_elevation()
45        #culvert_slope = delta_z/self.culvert.get_self.culvert_length()
46
47        # Determine controlling energy (driving head) for culvert
48        #if self.self.inflow.get_average_specific_energy() > self.self.delta_total_energy:
49        #    # Outlet control
50        #    driving_head = self.delta_total_energy
51        #else:
52            # Inlet control
53        #    driving_head = self.inflow.get_average_specific_energy()
54           
55       
56    def get_inflow(self):
57       
58        return self.inflow
59       
60       
61    def get_outflow(self):
62   
63        return self.outflow
64   
65       
66    def boyd_box(self):
67
68        """Boyd's generalisation of the US department of transportation culvert methods
69       
70        WARNING THIS IS A SIMPLISTIC APPROACH and OUTLET VELOCITIES ARE LIMITED TO EITHER
71        FULL PIPE OR CRITICAL DEPTH ONLY
72        For Supercritical flow this is UNDERESTIMATING the Outlet Velocity
73        The obtain the CORRECT velocity requires an iteration of Depth to Establish
74        the Normal Depth of flow in the pipe.
75       
76        It is proposed to provide this in a seperate routine called
77        boyd_generalised_culvert_model_complex
78       
79        The Boyd Method is based on methods described by the following:
80        1.
81        US Dept. Transportation Federal Highway Administration (1965)
82        Hydraulic Chart for Selection of Highway Culverts.
83        Hydraulic Engineering Circular No. 5 US Government Printing
84        2.
85        US Dept. Transportation Federal Highway Administration (1972)
86        Capacity charts for the Hydraulic design of highway culverts.
87        Hydraulic Engineering Circular No. 10 US Government Printing
88        These documents provide around 60 charts for various configurations of culverts and inlets.
89       
90        Note these documents have been superceded by:
91        2005 Hydraulic Design of Highway Culverts, Hydraulic Design Series No. 5 (HDS-5),
92        Which combines culvert design information previously contained in Hydraulic Engineering Circulars
93        (HEC) No. 5, No. 10, and No. 13 with hydrologic, storage routing, and special culvert design information.
94        HEC-5 provides 20 Charts
95        HEC-10 Provides an additional 36 Charts
96        HEC-13 Discusses the Design of improved more efficient inlets
97        HDS-5 Provides 60 sets of Charts
98
99        In 1985 Professor Michael Boyd Published "Head-Discharge Relations for Culverts", and in
100        1987 published "Generalised Head Discharge Equations for Culverts".
101        These papers reviewed the previous work by the US DOT and provided a simplistic approach for 3 configurations.
102       
103        It may be possible to extend the same approach for additional charts in the original work, but to date this has not been done.
104        The additional charts cover a range of culvert shapes and inlet configurations
105       
106       
107        """
108
109        local_debug ='false'
110        if self.inflow.get_average_height() > 0.1: #this value was 0.01:
111            if local_debug =='true':
112                log.critical('Specific E & Deltat Tot E = %s, %s'
113                             % (str(self.inflow.get_average_specific_energy()),
114                                str(self.delta_total_energy)))
115                log.critical('culvert type = %s' % str(culvert_type))
116            # Water has risen above inlet
117           
118            if self.log_filename is not None:                           
119                s = 'Specific energy  = %f m' % self.inflow.get_average_specific_energy()
120                log_to_file(self.log_filename, s)
121           
122            msg = 'Specific energy at inlet is negative'
123            assert self.inflow.get_average_specific_energy() >= 0.0, msg
124                         
125            height = self.culvert_height
126            width = self.culvert_width
127            flow_width = self.culvert_width           
128           
129            Q_inlet_unsubmerged = 0.540*g**0.5*width*self.inflow.get_average_specific_energy()**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
130            Q_inlet_submerged = 0.702*g**0.5*width*height**0.89*self.inflow.get_average_specific_energy()**0.61  # Flow based on Inlet Ctrl Inlet Submerged
131
132            # FIXME(Ole): Are these functions really for inlet control?   
133            if Q_inlet_unsubmerged < Q_inlet_submerged:
134                Q = Q_inlet_unsubmerged
135                dcrit = (Q**2/g/width**2)**0.333333
136                if dcrit > height:
137                    dcrit = height
138                flow_area = width*dcrit
139                outlet_culvert_depth = dcrit
140                case = 'Inlet unsubmerged Box Acts as Weir'
141            else:   
142                Q = Q_inlet_submerged
143                flow_area = width*height
144                outlet_culvert_depth = height
145                case = 'Inlet submerged Box Acts as Orifice'                   
146
147            dcrit = (Q**2/g/width**2)**0.333333
148
149            outlet_culvert_depth = dcrit
150            if outlet_culvert_depth > height:
151                outlet_culvert_depth = height  # Once again the pipe is flowing full not partfull
152                flow_area = width*height  # Cross sectional area of flow in the culvert
153                perimeter = 2*(width+height)
154                case = 'Inlet CTRL Outlet unsubmerged PIPE PART FULL'
155            else:
156                flow_area = width * outlet_culvert_depth
157                perimeter = width+2*outlet_culvert_depth
158                case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'
159
160            if self.delta_total_energy < self.inflow.get_average_specific_energy():
161                # Calculate flows for outlet control
162               
163                # Determine the depth at the outlet relative to the depth of flow in the Culvert
164                if self.outflow.get_average_height() > height:        # The Outlet is Submerged
165                    outlet_culvert_depth=height
166                    flow_area=width*height       # Cross sectional area of flow in the culvert
167                    perimeter=2.0*(width+height)
168                    case = 'Outlet submerged'
169                else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
170                    dcrit = (Q**2/g/width**2)**0.333333
171                    outlet_culvert_depth=dcrit   # For purpose of calculation assume the outlet depth = Critical Depth
172                    if outlet_culvert_depth > height:
173                        outlet_culvert_depth=height
174                        flow_area=width*height
175                        perimeter=2.0*(width+height)
176                        case = 'Outlet is Flowing Full'
177                    else:
178                        flow_area=width*outlet_culvert_depth
179                        perimeter=(width+2.0*outlet_culvert_depth)
180                        case = 'Outlet is open channel flow'
181
182                hyd_rad = flow_area/perimeter
183
184                if self.log_filename is not None:                               
185                    s = 'hydraulic radius at outlet = %f' % hyd_rad
186                    log_to_file(self.log_filename, s)
187
188                # Outlet control velocity using tail water
189                culvert_velocity = sqrt(self.delta_total_energy/((self.sum_loss/2/g)+(self.manning**2*self.culvert_length)/hyd_rad**1.33333)) 
190                Q_outlet_tailwater = flow_area * culvert_velocity
191
192                if self.log_filename is not None:                           
193                    s = 'Q_outlet_tailwater = %.6f' % Q_outlet_tailwater
194                    log_to_file(self.log_filename, s)
195                Q = min(Q, Q_outlet_tailwater)
196            else:
197                pass
198                #FIXME(Ole): What about inlet control?
199
200            culv_froude=sqrt(Q**2*flow_width/(g*flow_area**3))
201            if local_debug =='true':
202                log.critical('FLOW AREA = %s' % str(flow_area))
203                log.critical('PERIMETER = %s' % str(perimeter))
204                log.critical('Q final = %s' % str(Q))
205                log.critical('FROUDE = %s' % str(culv_froude))
206     
207            # Determine momentum at the outlet
208            barrel_velocity = Q/(flow_area + velocity_protection/flow_area)
209
210        # END CODE BLOCK for DEPTH  > Required depth for CULVERT Flow 
211
212        else: # self.inflow.get_average_height() < 0.01:
213            Q = barrel_velocity = outlet_culvert_depth = 0.0
214
215        # Temporary flow limit
216        if barrel_velocity > self.max_velocity:
217            barrel_velocity = self.max_velocity
218            Q = flow_area * barrel_velocity
219           
220        return Q, barrel_velocity, outlet_culvert_depth
221
222
223    def boyd_circle(self):
224        """
225        For a circular pipe the Boyd method reviews 3 conditions
226        1. Whether the Pipe Inlet is Unsubmerged (acting as weir flow into the inlet)
227        2. Whether the Pipe Inlet is Fully Submerged (acting as an Orifice)
228        3. Whether the energy loss in the pipe results in the Pipe being controlled by Channel Flow of the Pipe
229       
230        For these conditions we also would like to assess the pipe flow characteristics as it leaves the pipe
231        """
232       
233        diameter = self.culvert_height
234       
235        local_debug ='false'
236        if self.inflow.get_average_height() > 0.1: #this value was 0.01:
237            if local_debug =='true':
238                log.critical('Specific E & Deltat Tot E = %s, %s'
239                             % (str(self.inflow.get_average_specific_energy()),
240                                str(self.delta_total_energy)))
241                log.critical('culvert type = %s' % str(culvert_type))
242            # Water has risen above inlet
243           
244            if self.log_filename is not None:                           
245                s = 'Specific energy  = %f m' % self.inflow.get_average_specific_energy()
246                log_to_file(self.log_filename, s)
247           
248            msg = 'Specific energy at inlet is negative'
249            assert self.inflow.get_average_specific_energy() >= 0.0, msg
250                         
251            # Calculate flows for inlet control           
252            Q_inlet_unsubmerged = 0.421*g**0.5*diameter**0.87*self.inflow.get_average_specific_energy()**1.63 # Inlet Ctrl Inlet Unsubmerged
253            Q_inlet_submerged = 0.530*g**0.5*diameter**1.87*self.inflow.get_average_specific_energy()**0.63   # Inlet Ctrl Inlet Submerged
254            # Note for to SUBMERGED TO OCCUR self.inflow.get_average_specific_energy() should be > 1.2 x diameter.... Should Check !!!
255
256            if self.log_filename is not None:                                       
257                s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' % (Q_inlet_unsubmerged, Q_inlet_submerged)
258                log_to_file(self.log_filename, s)
259            Q = min(Q_inlet_unsubmerged, Q_inlet_submerged)
260
261            # THE LOWEST Value will Control Calcs From here
262            # Calculate Critical Depth Based on the Adopted Flow as an Estimate
263            dcrit1 = diameter/1.26*(Q/g**0.5*diameter**2.5)**(1/3.75)
264            dcrit2 = diameter/0.95*(Q/g**0.5*diameter**2.5)**(1/1.95)
265            # From Boyd Paper ESTIMATE of Dcrit has 2 criteria as
266            if dcrit1/diameter  > 0.85:
267                outlet_culvert_depth = dcrit2
268            else:
269                outlet_culvert_depth = dcrit1
270            #outlet_culvert_depth = min(outlet_culvert_depth, diameter)
271            # Now determine Hydraulic Radius Parameters Area & Wetted Perimeter
272            if outlet_culvert_depth >= diameter:
273                outlet_culvert_depth = diameter  # Once again the pipe is flowing full not partfull
274                flow_area = (diameter/2)**2 * pi  # Cross sectional area of flow in the culvert
275                perimeter = diameter * pi
276                flow_width= diameter
277                case = 'Inlet CTRL Outlet submerged Circular PIPE FULL'
278                if local_debug == 'true':
279                    log.critical('Inlet CTRL Outlet submerged Circular '
280                                 'PIPE FULL')
281            else:
282                #alpha = acos(1 - outlet_culvert_depth/diameter)    # Where did this Come From ????/
283                alpha = acos(1-2*outlet_culvert_depth/diameter)*2
284                #flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))        # Pipe is Running Partly Full at the INLET   WHRE did this Come From ?????
285                flow_area = diameter**2/8*(alpha - sin(alpha))   # Equation from  GIECK 5th Ed. Pg. B3
286                flow_width= diameter*sin(alpha/2.0)
287                perimeter = alpha*diameter/2.0
288                case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'
289                if local_debug =='true':
290                    log.critical('INLET CTRL Culvert is open channel flow '
291                                 'we will for now assume critical depth')
292                    log.critical('Q Outlet Depth and ALPHA = %s, %s, %s'
293                                 % (str(Q), str(outlet_culvert_depth),
294                                    str(alpha)))
295            if self.delta_total_energy < self.inflow.get_average_specific_energy():  #  OUTLET CONTROL !!!!
296                # Calculate flows for outlet control
297               
298                # Determine the depth at the outlet relative to the depth of flow in the Culvert
299                if self.outflow.get_average_height() > diameter:       # Outlet is submerged Assume the end of the Pipe is flowing FULL
300                    outlet_culvert_depth=diameter
301                    flow_area = (diameter/2)**2 * pi  # Cross sectional area of flow in the culvert
302                    perimeter = diameter * pi
303                    flow_width= diameter
304                    case = 'Outlet submerged'
305                    if local_debug =='true':
306                        log.critical('Outlet submerged')
307                else:   # Culvert running PART FULL for PART OF ITS LENGTH   Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
308                    # IF  self.outflow.get_average_height() < diameter
309                    dcrit1 = diameter/1.26*(Q/g**0.5*diameter**2.5)**(1/3.75)
310                    dcrit2 = diameter/0.95*(Q/g**0.5*diameter**2.5)**(1/1.95)
311                    if dcrit1/diameter >0.85:
312                        outlet_culvert_depth= dcrit2
313                    else:
314                        outlet_culvert_depth = dcrit1
315                    if outlet_culvert_depth > diameter:
316                        outlet_culvert_depth = diameter  # Once again the pipe is flowing full not partfull
317                        flow_area = (diameter/2)**2 * pi  # Cross sectional area of flow in the culvert
318                        perimeter = diameter * pi
319                        flow_width= diameter
320                        case = 'Outlet unsubmerged PIPE FULL'
321                        if local_debug =='true':
322                            log.critical('Outlet unsubmerged PIPE FULL')
323                    else:
324                        alpha = acos(1-2*outlet_culvert_depth/diameter)*2
325                        flow_area = diameter**2/8*(alpha - sin(alpha))   # Equation from  GIECK 5th Ed. Pg. B3
326                        flow_width= diameter*sin(alpha/2.0)
327                        perimeter = alpha*diameter/2.0
328                        case = 'Outlet is open channel flow we will for now assume critical depth'
329                        if local_debug == 'true':
330                            log.critical('Q Outlet Depth and ALPHA = %s, %s, %s'
331                                         % (str(Q), str(outlet_culvert_depth),
332                                            str(alpha)))
333                            log.critical('Outlet is open channel flow we '
334                                         'will for now assume critical depth')
335            if local_debug == 'true':
336                log.critical('FLOW AREA = %s' % str(flow_area))
337                log.critical('PERIMETER = %s' % str(perimeter))
338                log.critical('Q Interim = %s' % str(Q))
339            hyd_rad = flow_area/perimeter
340
341            if self.log_filename is not None:
342                s = 'hydraulic radius at outlet = %f' %hyd_rad
343                log_to_file(self.log_filename, s)
344
345            # Outlet control velocity using tail water
346            if local_debug =='true':
347                log.critical('GOT IT ALL CALCULATING Velocity')
348                log.critical('HydRad = %s' % str(hyd_rad))
349            culvert_velocity = sqrt(self.delta_total_energy/((self.sum_loss/2/g)+(self.manning**2*self.culvert_length)/hyd_rad**1.33333)) 
350            Q_outlet_tailwater = flow_area * culvert_velocity
351            if local_debug =='true':
352                log.critical('VELOCITY = %s' % str(culvert_velocity))
353                log.critical('Outlet Ctrl Q = %s' % str(Q_outlet_tailwater))
354            if self.log_filename is not None:               
355                s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater
356                log_to_file(self.log_filename, s)
357            Q = min(Q, Q_outlet_tailwater)
358            if local_debug =='true':
359                log.critical('%s,%.3f,%.3f'
360                             % ('dcrit 1 , dcit2 =',dcrit1,dcrit2))
361                log.critical('%s,%.3f,%.3f,%.3f'
362                             % ('Q and Velocity and Depth=', Q,
363                                culvert_velocity, outlet_culvert_depth))
364                               
365            culv_froude=sqrt(Q**2*flow_width/(g*flow_area**3))
366            if local_debug =='true':
367                log.critical('FLOW AREA = %s' % str(flow_area))
368                log.critical('PERIMETER = %s' % str(perimeter))
369                log.critical('Q final = %s' % str(Q))
370                log.critical('FROUDE = %s' % str(culv_froude))
371
372            # Determine momentum at the outlet
373            barrel_velocity = Q/(flow_area + velocity_protection/flow_area)
374
375        else: # self.inflow.get_average_height() < 0.01:
376            Q = barrel_velocity = outlet_culvert_depth = 0.0
377
378        # Temporary flow limit
379        if barrel_velocity > self.max_velocity:       
380            barrel_velocity = self.max_velocity
381            Q = flow_area * barrel_velocity
382       
383        return Q, barrel_velocity, outlet_culvert_depth
384
385
Note: See TracBrowser for help on using the repository browser.