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

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

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