source: trunk/anuga_core/source/anuga/structures/boyd_pipe_operator.py @ 8049

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

Added in Rudy's amended boyd_box_operator

  • Property svn:executable set to *
File size: 12.4 KB
Line 
1import anuga
2import math
3import types
4
5class Boyd_pipe_operator(anuga.Structure_operator):
6    """Culvert flow - transfer water from one location to another via a circular pipe culvert.
7    Sets up the geometry of problem
8   
9    This is the base class for culverts. Inherit from this class (and overwrite
10    compute_discharge method for specific subclasses)
11   
12    Input: Two points, pipe_size (diameter),
13    mannings_rougness,
14    """ 
15
16    def __init__(self,
17                 domain,
18                 end_point0, 
19                 end_point1,
20                 losses,
21                 diameter=None,
22                 apron=None,
23                 manning=0.013,
24                 enquiry_gap=0.2,
25                 use_momentum_jet=True,
26                 use_velocity_head=True,
27                 description=None,
28                 label=None,
29                 structure_type='boyd_pipe',
30                 logging=False,
31                 verbose=False):
32
33
34                     
35        anuga.Structure_operator.__init__(self,
36                                          domain,
37                                          end_point0, 
38                                          end_point1,
39                                          width=diameter,
40                                          height=None,
41                                          apron=apron,
42                                          manning=manning,
43                                          enquiry_gap=enquiry_gap,                                                       
44                                          description=description,
45                                          label=label,
46                                          structure_type=structure_type,
47                                          logging=logging,
48                                          verbose=verbose)           
49
50 
51        if type(losses) == types.DictType:
52            self.sum_loss = sum(losses.values())
53        elif type(losses) == types.ListType:
54            self.sum_loss = sum(losses)
55        else:
56            self.sum_loss = losses
57       
58        self.use_momentum_jet = use_momentum_jet
59        self.use_velocity_head = use_velocity_head
60       
61        self.culvert_length = self.get_culvert_length()
62        self.culvert_diameter = self.get_culvert_diameter()
63
64        self.max_velocity = 10.0
65
66        self.inlets = self.get_inlets()
67
68        # Stats
69       
70        self.discharge = 0.0
71        self.velocity = 0.0
72       
73        self.case = 'N/A'
74       
75   
76    def discharge_routine(self):
77
78        local_debug ='false'
79       
80        #import pdb
81        #pdb.set_trace()
82       
83        if self.inflow.get_enquiry_height() > 0.01: #this value was 0.01: Remember this needs to be compared to the Invert Lvl
84            if local_debug =='true':
85                anuga.log.critical('Specific E & Deltat Tot E = %s, %s'
86                             % (str(self.inflow.get_enquiry_specific_energy()),
87                                str(self.delta_total_energy)))
88                anuga.log.critical('culvert type = %s' % str(culvert_type))
89            # Water has risen above inlet
90
91
92            msg = 'Specific energy at inlet is negative'
93            assert self.inflow.get_enquiry_specific_energy() >= 0.0, msg
94
95            if self.use_velocity_head :
96                self.driving_energy = self.inflow.get_enquiry_specific_energy()
97            else:
98                self.driving_energy = self.inflow.get_enquiry_height()
99                """
100        For a circular pipe the Boyd method reviews 3 conditions
101        1. Whether the Pipe Inlet is Unsubmerged (acting as weir flow into the inlet)
102        2. Whether the Pipe Inlet is Fully Submerged (acting as an Orifice)
103        3. Whether the energy loss in the pipe results in the Pipe being controlled by Channel Flow of the Pipe
104
105        For these conditions we also would like to assess the pipe flow characteristics as it leaves the pipe
106        """
107
108        diameter = self.culvert_diameter
109
110        local_debug ='false'
111        if self.inflow.get_average_height() > 0.01: #this should test against invert
112            if local_debug =='true':
113                anuga.log.critical('Specific E & Deltat Tot E = %s, %s'
114                             % (str(self.inflow.get_average_specific_energy()),
115                                str(self.delta_total_energy)))
116                anuga.log.critical('culvert type = %s' % str(culvert_type))
117            # Water has risen above inlet
118
119
120            msg = 'Specific energy at inlet is negative'
121            assert self.inflow.get_average_specific_energy() >= 0.0, msg
122
123            # Calculate flows for inlet control for circular pipe
124            Q_inlet_unsubmerged = 0.421*anuga.g**0.5*diameter**0.87*self.inflow.get_average_specific_energy()**1.63 # Inlet Ctrl Inlet Unsubmerged
125            Q_inlet_submerged = 0.530*anuga.g**0.5*diameter**1.87*self.inflow.get_average_specific_energy()**0.63   # Inlet Ctrl Inlet Submerged
126            # Note for to SUBMERGED TO OCCUR self.inflow.get_average_specific_energy() should be > 1.2 x diameter.... Should Check !!!
127
128
129            Q = min(Q_inlet_unsubmerged, Q_inlet_submerged)
130
131            # THE LOWEST Value will Control Calcs From here
132            # Calculate Critical Depth Based on the Adopted Flow as an Estimate
133            dcrit1 = diameter/1.26*(Q/anuga.g**0.5*diameter**2.5)**(1/3.75)
134            dcrit2 = diameter/0.95*(Q/anuga.g**0.5*diameter**2.5)**(1/1.95)
135            # From Boyd Paper ESTIMATE of Dcrit has 2 criteria as
136            if dcrit1/diameter  > 0.85:
137                outlet_culvert_depth = dcrit2
138            else:
139                outlet_culvert_depth = dcrit1
140            #outlet_culvert_depth = min(outlet_culvert_depth, diameter)
141            # Now determine Hydraulic Radius Parameters Area & Wetted Perimeter
142            if outlet_culvert_depth >= diameter:
143                outlet_culvert_depth = diameter  # Once again the pipe is flowing full not partfull
144                flow_area = (diameter/2)**2 * math.pi  # Cross sectional area of flow in the culvert
145                perimeter = diameter * math.pi
146                flow_width= diameter
147                self.case = 'Inlet CTRL Outlet submerged Circular PIPE FULL'
148                if local_debug == 'true':
149                    anuga.log.critical('Inlet CTRL Outlet submerged Circular '
150                                 'PIPE FULL')
151            else:
152                #alpha = anuga.acos(1 - outlet_culvert_depth/diameter)    # Where did this Come From ????/
153                alpha = anuga.acos(1-2*outlet_culvert_depth/diameter)*2
154                #flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))        # Pipe is Running Partly Full at the INLET   WHRE did this Come From ?????
155                flow_area = diameter**2/8*(alpha - math.sin(alpha))   # Equation from  GIECK 5th Ed. Pg. B3
156                flow_width= diameter*math.sin(alpha/2.0)
157                perimeter = alpha*diameter/2.0
158                self.case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'
159                if local_debug =='true':
160                    anuga.log.critical('INLET CTRL Culvert is open channel flow '
161                                 'we will for now assume critical depth')
162                    anuga.log.critical('Q Outlet Depth and ALPHA = %s, %s, %s'
163                                 % (str(Q), str(outlet_culvert_depth),
164                                    str(alpha)))
165            if self.delta_total_energy < self.inflow.get_average_specific_energy():  #  OUTLET CONTROL !!!!
166                # Calculate flows for outlet control
167
168                # Determine the depth at the outlet relative to the depth of flow in the Culvert
169                if self.outflow.get_average_height() > diameter:       # Outlet is submerged Assume the end of the Pipe is flowing FULL
170                    outlet_culvert_depth=diameter
171                    flow_area = (diameter/2)**2 * math.pi  # Cross sectional area of flow in the culvert
172                    perimeter = diameter * math.pi
173                    flow_width= diameter
174                    self.case = 'Outlet submerged'
175                    if local_debug =='true':
176                        anuga.log.critical('Outlet submerged')
177                else:   # Culvert running PART FULL for PART OF ITS LENGTH   Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
178                    # IF  self.outflow.get_average_height() < diameter
179                    dcrit1 = diameter/1.26*(Q/anuga.g**0.5*diameter**2.5)**(1/3.75)
180                    dcrit2 = diameter/0.95*(Q/anuga.g**0.5*diameter**2.5)**(1/1.95)
181                    if dcrit1/diameter >0.85:
182                        outlet_culvert_depth= dcrit2
183                    else:
184                        outlet_culvert_depth = dcrit1
185                    if outlet_culvert_depth > diameter:
186                        outlet_culvert_depth = diameter  # Once again the pipe is flowing full not partfull
187                        flow_area = (diameter/2)**2 * math.pi  # Cross sectional area of flow in the culvert
188                        perimeter = diameter * math.pi
189                        flow_width= diameter
190                        self.case = 'Outlet unsubmerged PIPE FULL'
191                        if local_debug =='true':
192                            anuga.log.critical('Outlet unsubmerged PIPE FULL')
193                    else:
194                        alpha = anuga.acos(1-2*outlet_culvert_depth/diameter)*2
195                        flow_area = diameter**2/8*(alpha - math.sin(alpha))   # Equation from  GIECK 5th Ed. Pg. B3
196                        flow_width= diameter*math.sin(alpha/2.0)
197                        perimeter = alpha*diameter/2.0
198                        self.case = 'Outlet is open channel flow we will for now assume critical depth'
199                        if local_debug == 'true':
200                            anuga.log.critical('Q Outlet Depth and ALPHA = %s, %s, %s'
201                                         % (str(Q), str(outlet_culvert_depth),
202                                            str(alpha)))
203                            anuga.log.critical('Outlet is open channel flow we '
204                                         'will for now assume critical depth')
205            if local_debug == 'true':
206                anuga.log.critical('FLOW AREA = %s' % str(flow_area))
207                anuga.log.critical('PERIMETER = %s' % str(perimeter))
208                anuga.log.critical('Q Interim = %s' % str(Q))
209            hyd_rad = flow_area/perimeter
210
211
212
213            # Outlet control velocity using tail water
214            if local_debug =='true':
215                anuga.log.critical('GOT IT ALL CALCULATING Velocity')
216                anuga.log.critical('HydRad = %s' % str(hyd_rad))
217            # Calculate Pipe Culvert Outlet Control Velocity.... May need initial Estimate First ??
218           
219            culvert_velocity = math.sqrt(self.delta_total_energy/((self.sum_loss/2/anuga.g)+(self.manning**2*self.culvert_length)/hyd_rad**1.33333))
220            Q_outlet_tailwater = flow_area * culvert_velocity
221           
222           
223            if local_debug =='true':
224                anuga.log.critical('VELOCITY = %s' % str(culvert_velocity))
225                anuga.log.critical('Outlet Ctrl Q = %s' % str(Q_outlet_tailwater))
226
227            Q = min(Q, Q_outlet_tailwater)
228            if local_debug =='true':
229                anuga.log.critical('%s,%.3f,%.3f'
230                             % ('dcrit 1 , dcit2 =',dcrit1,dcrit2))
231                anuga.log.critical('%s,%.3f,%.3f,%.3f'
232                             % ('Q and Velocity and Depth=', Q,
233                                culvert_velocity, outlet_culvert_depth))
234
235            culv_froude=math.sqrt(Q**2*flow_width/(anuga.g*flow_area**3))
236            if local_debug =='true':
237                anuga.log.critical('FLOW AREA = %s' % str(flow_area))
238                anuga.log.critical('PERIMETER = %s' % str(perimeter))
239                anuga.log.critical('Q final = %s' % str(Q))
240                anuga.log.critical('FROUDE = %s' % str(culv_froude))
241
242            # Determine momentum at the outlet
243            barrel_velocity = Q/(flow_area + anuga.velocity_protection/flow_area)
244
245        else: # self.inflow.get_average_height() < 0.01:
246            Q = barrel_velocity = outlet_culvert_depth = 0.0
247
248        # Temporary flow limit
249        if barrel_velocity > self.max_velocity:
250            barrel_velocity = self.max_velocity
251            Q = flow_area * barrel_velocity
252
253        return Q, barrel_velocity, outlet_culvert_depth
254
255       
256       
257       
Note: See TracBrowser for help on using the repository browser.