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

Last change on this file since 8002 was 8002, checked in by habili, 12 years ago

using "import anuga"

  • Property svn:executable set to *
File size: 16.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                 verbose=False):
29                     
30        anuga.Structure_operator.__init__(self,
31                                          domain,
32                                          end_point0, 
33                                          end_point1,
34                                          width=diameter,
35                                          height=None,
36                                          apron=apron,
37                                          manning=manning,
38                                          enquiry_gap=enquiry_gap,                                                       
39                                          description=description,
40                                          verbose=verbose)           
41
42 
43        if type(losses) == types.DictType:
44            self.sum_loss = sum(losses.values())
45        elif type(losses) == types.ListType:
46            self.sum_loss = sum(losses)
47        else:
48            self.sum_loss = losses
49       
50        self.use_momentum_jet = use_momentum_jet
51        self.use_velocity_head = use_velocity_head
52       
53        self.culvert_length = self.get_culvert_length()
54        self.culvert_diameter = self.get_culvert_diameter()
55
56        self.max_velocity = 10.0
57        self.log_filename = None
58
59        self.inlets = self.get_inlets()
60
61
62        # Stats
63       
64        self.discharge = 0.0
65        self.velocity = 0.0
66       
67       
68    def __call__(self):
69       
70        timestep = self.domain.get_timestep()
71       
72        self.__determine_inflow_outflow()
73       
74        Q, barrel_speed, outlet_depth = self.__discharge_routine()
75
76        #inflow  = self.routine.get_inflow()
77        #outflow = self.routine.get_outflow()
78
79        old_inflow_height = self.inflow.get_average_height()
80        old_inflow_xmom = self.inflow.get_average_xmom()
81        old_inflow_ymom = self.inflow.get_average_ymom()
82           
83        if old_inflow_height > 0.0 :
84                Qstar = Q/old_inflow_height
85        else:
86                Qstar = 0.0
87
88        factor = 1.0/(1.0 + Qstar*timestep/self.inflow.get_area())
89
90        new_inflow_height = old_inflow_height*factor
91        new_inflow_xmom = old_inflow_xmom*factor
92        new_inflow_ymom = old_inflow_ymom*factor
93           
94
95        self.inflow.set_heights(new_inflow_height)
96
97        #inflow.set_xmoms(Q/inflow.get_area())
98        #inflow.set_ymoms(0.0)
99
100
101        self.inflow.set_xmoms(new_inflow_xmom)
102        self.inflow.set_ymoms(new_inflow_ymom)
103
104
105        loss = (old_inflow_height - new_inflow_height)*self.inflow.get_area()
106
107           
108        # set outflow
109        if old_inflow_height > 0.0 :
110                timestep_star = timestep*new_inflow_height/old_inflow_height
111        else:
112            timestep_star = 0.0
113
114           
115        outflow_extra_height = Q*timestep_star/self.outflow.get_area()
116        outflow_direction = - self.outflow.outward_culvert_vector
117        outflow_extra_momentum = outflow_extra_height*barrel_speed*outflow_direction
118           
119
120        gain = outflow_extra_height*self.outflow.get_area()
121           
122        #print Q, Q*timestep, barrel_speed, outlet_depth, Qstar, factor, timestep_star
123        #print '  ', loss, gain
124
125        # Stats
126        self.discharge  = Q#outflow_extra_height*self.outflow.get_area()/timestep
127        self.velocity = barrel_speed#self.discharge/outlet_depth/self.width
128
129        new_outflow_height = self.outflow.get_average_height() + outflow_extra_height
130
131        if self.use_momentum_jet :
132            # FIXME (SR) Review momentum to account for possible hydraulic jumps at outlet
133            #new_outflow_xmom = outflow.get_average_xmom() + outflow_extra_momentum[0]
134            #new_outflow_ymom = outflow.get_average_ymom() + outflow_extra_momentum[1]
135
136            new_outflow_xmom = barrel_speed*new_outflow_height*outflow_direction[0]
137            new_outflow_ymom = barrel_speed*new_outflow_height*outflow_direction[1]
138
139        else:
140            #new_outflow_xmom = outflow.get_average_xmom()
141            #new_outflow_ymom = outflow.get_average_ymom()
142
143            new_outflow_xmom = 0.0
144            new_outflow_ymom = 0.0
145
146
147        self.outflow.set_heights(new_outflow_height)
148        self.outflow.set_xmoms(new_outflow_xmom)
149        self.outflow.set_ymoms(new_outflow_ymom)
150
151
152    def __determine_inflow_outflow(self):
153        # Determine flow direction based on total energy difference
154
155        if self.use_velocity_head:
156            self.delta_total_energy = self.inlets[0].get_enquiry_total_energy() - self.inlets[1].get_enquiry_total_energy()
157        else:
158            self.delta_total_energy = self.inlets[0].get_enquiry_stage() - self.inlets[1].get_enquiry_stage()
159
160
161        self.inflow  = self.inlets[0]
162        self.outflow = self.inlets[1]
163       
164
165        if self.delta_total_energy < 0:
166            self.inflow  = self.inlets[1]
167            self.outflow = self.inlets[0]
168            self.delta_total_energy = -self.delta_total_energy
169
170   
171    def __discharge_routine(self):
172
173        local_debug ='false'
174       
175        if self.inflow.get_enquiry_height() > 0.01: #this value was 0.01: Remember this needs to be compared to the Invert Lvl
176            if local_debug =='true':
177                anuga.log.critical('Specific E & Deltat Tot E = %s, %s'
178                             % (str(self.inflow.get_enquiry_specific_energy()),
179                                str(self.delta_total_energy)))
180                anuga.log.critical('culvert type = %s' % str(culvert_type))
181            # Water has risen above inlet
182
183            if self.log_filename is not None:
184                s = 'Specific energy  = %f m' % self.inflow.get_enquiry_specific_energy()
185                log_to_file(self.log_filename, s)
186
187            msg = 'Specific energy at inlet is negative'
188            assert self.inflow.get_enquiry_specific_energy() >= 0.0, msg
189
190            if self.use_velocity_head :
191                self.driving_energy = self.inflow.get_enquiry_specific_energy()
192            else:
193                self.driving_energy = self.inflow.get_enquiry_height()
194                """
195        For a circular pipe the Boyd method reviews 3 conditions
196        1. Whether the Pipe Inlet is Unsubmerged (acting as weir flow into the inlet)
197        2. Whether the Pipe Inlet is Fully Submerged (acting as an Orifice)
198        3. Whether the energy loss in the pipe results in the Pipe being controlled by Channel Flow of the Pipe
199
200        For these conditions we also would like to assess the pipe flow characteristics as it leaves the pipe
201        """
202
203        diameter = self.culvert_diameter
204
205        local_debug ='false'
206        if self.inflow.get_average_height() > 0.01: #this should test against invert
207            if local_debug =='true':
208                anuga.log.critical('Specific E & Deltat Tot E = %s, %s'
209                             % (str(self.inflow.get_average_specific_energy()),
210                                str(self.delta_total_energy)))
211                anuga.log.critical('culvert type = %s' % str(culvert_type))
212            # Water has risen above inlet
213
214            if self.log_filename is not None:
215                s = 'Specific energy  = %f m' % self.inflow.get_average_specific_energy()
216                log_to_file(self.log_filename, s)
217
218            msg = 'Specific energy at inlet is negative'
219            assert self.inflow.get_average_specific_energy() >= 0.0, msg
220
221            # Calculate flows for inlet control for circular pipe
222            Q_inlet_unsubmerged = 0.421*anuga.g**0.5*diameter**0.87*self.inflow.get_average_specific_energy()**1.63 # Inlet Ctrl Inlet Unsubmerged
223            Q_inlet_submerged = 0.530*anuga.g**0.5*diameter**1.87*self.inflow.get_average_specific_energy()**0.63   # Inlet Ctrl Inlet Submerged
224            # Note for to SUBMERGED TO OCCUR self.inflow.get_average_specific_energy() should be > 1.2 x diameter.... Should Check !!!
225
226            if self.log_filename is not None:
227                s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' % (Q_inlet_unsubmerged, Q_inlet_submerged)
228                log_to_file(self.log_filename, s)
229            Q = min(Q_inlet_unsubmerged, Q_inlet_submerged)
230
231            # THE LOWEST Value will Control Calcs From here
232            # Calculate Critical Depth Based on the Adopted Flow as an Estimate
233            dcrit1 = diameter/1.26*(Q/anuga.g**0.5*diameter**2.5)**(1/3.75)
234            dcrit2 = diameter/0.95*(Q/anuga.g**0.5*diameter**2.5)**(1/1.95)
235            # From Boyd Paper ESTIMATE of Dcrit has 2 criteria as
236            if dcrit1/diameter  > 0.85:
237                outlet_culvert_depth = dcrit2
238            else:
239                outlet_culvert_depth = dcrit1
240            #outlet_culvert_depth = min(outlet_culvert_depth, diameter)
241            # Now determine Hydraulic Radius Parameters Area & Wetted Perimeter
242            if outlet_culvert_depth >= diameter:
243                outlet_culvert_depth = diameter  # Once again the pipe is flowing full not partfull
244                flow_area = (diameter/2)**2 * math.pi  # Cross sectional area of flow in the culvert
245                perimeter = diameter * math.pi
246                flow_width= diameter
247                case = 'Inlet CTRL Outlet submerged Circular PIPE FULL'
248                if local_debug == 'true':
249                    anuga.log.critical('Inlet CTRL Outlet submerged Circular '
250                                 'PIPE FULL')
251            else:
252                #alpha = anuga.acos(1 - outlet_culvert_depth/diameter)    # Where did this Come From ????/
253                alpha = anuga.acos(1-2*outlet_culvert_depth/diameter)*2
254                #flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))        # Pipe is Running Partly Full at the INLET   WHRE did this Come From ?????
255                flow_area = diameter**2/8*(alpha - math.sin(alpha))   # Equation from  GIECK 5th Ed. Pg. B3
256                flow_width= diameter*math.sin(alpha/2.0)
257                perimeter = alpha*diameter/2.0
258                case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'
259                if local_debug =='true':
260                    anuga.log.critical('INLET CTRL Culvert is open channel flow '
261                                 'we will for now assume critical depth')
262                    anuga.log.critical('Q Outlet Depth and ALPHA = %s, %s, %s'
263                                 % (str(Q), str(outlet_culvert_depth),
264                                    str(alpha)))
265            if self.delta_total_energy < self.inflow.get_average_specific_energy():  #  OUTLET CONTROL !!!!
266                # Calculate flows for outlet control
267
268                # Determine the depth at the outlet relative to the depth of flow in the Culvert
269                if self.outflow.get_average_height() > diameter:       # Outlet is submerged Assume the end of the Pipe is flowing FULL
270                    outlet_culvert_depth=diameter
271                    flow_area = (diameter/2)**2 * math.pi  # Cross sectional area of flow in the culvert
272                    perimeter = diameter * math.pi
273                    flow_width= diameter
274                    case = 'Outlet submerged'
275                    if local_debug =='true':
276                        anuga.log.critical('Outlet submerged')
277                else:   # Culvert running PART FULL for PART OF ITS LENGTH   Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
278                    # IF  self.outflow.get_average_height() < diameter
279                    dcrit1 = diameter/1.26*(Q/anuga.g**0.5*diameter**2.5)**(1/3.75)
280                    dcrit2 = diameter/0.95*(Q/anuga.g**0.5*diameter**2.5)**(1/1.95)
281                    if dcrit1/diameter >0.85:
282                        outlet_culvert_depth= dcrit2
283                    else:
284                        outlet_culvert_depth = dcrit1
285                    if outlet_culvert_depth > diameter:
286                        outlet_culvert_depth = diameter  # Once again the pipe is flowing full not partfull
287                        flow_area = (diameter/2)**2 * math.pi  # Cross sectional area of flow in the culvert
288                        perimeter = diameter * math.pi
289                        flow_width= diameter
290                        case = 'Outlet unsubmerged PIPE FULL'
291                        if local_debug =='true':
292                            anuga.log.critical('Outlet unsubmerged PIPE FULL')
293                    else:
294                        alpha = anuga.acos(1-2*outlet_culvert_depth/diameter)*2
295                        flow_area = diameter**2/8*(alpha - math.sin(alpha))   # Equation from  GIECK 5th Ed. Pg. B3
296                        flow_width= diameter*math.sin(alpha/2.0)
297                        perimeter = alpha*diameter/2.0
298                        case = 'Outlet is open channel flow we will for now assume critical depth'
299                        if local_debug == 'true':
300                            anuga.log.critical('Q Outlet Depth and ALPHA = %s, %s, %s'
301                                         % (str(Q), str(outlet_culvert_depth),
302                                            str(alpha)))
303                            anuga.log.critical('Outlet is open channel flow we '
304                                         'will for now assume critical depth')
305            if local_debug == 'true':
306                anuga.log.critical('FLOW AREA = %s' % str(flow_area))
307                anuga.log.critical('PERIMETER = %s' % str(perimeter))
308                anuga.log.critical('Q Interim = %s' % str(Q))
309            hyd_rad = flow_area/perimeter
310
311            if self.log_filename is not None:
312                s = 'hydraulic radius at outlet = %f' %hyd_rad
313                log_to_file(self.log_filename, s)
314
315            # Outlet control velocity using tail water
316            if local_debug =='true':
317                anuga.log.critical('GOT IT ALL CALCULATING Velocity')
318                anuga.log.critical('HydRad = %s' % str(hyd_rad))
319            # Calculate Pipe Culvert Outlet Control Velocity.... May need initial Estimate First ??
320           
321            culvert_velocity = math.sqrt(self.delta_total_energy/((self.sum_loss/2/anuga.g)+(self.manning**2*self.culvert_length)/hyd_rad**1.33333))
322            Q_outlet_tailwater = flow_area * culvert_velocity
323           
324           
325            if local_debug =='true':
326                anuga.log.critical('VELOCITY = %s' % str(culvert_velocity))
327                anuga.log.critical('Outlet Ctrl Q = %s' % str(Q_outlet_tailwater))
328            if self.log_filename is not None:
329                s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater
330                log_to_file(self.log_filename, s)
331            Q = min(Q, Q_outlet_tailwater)
332            if local_debug =='true':
333                anuga.log.critical('%s,%.3f,%.3f'
334                             % ('dcrit 1 , dcit2 =',dcrit1,dcrit2))
335                anuga.log.critical('%s,%.3f,%.3f,%.3f'
336                             % ('Q and Velocity and Depth=', Q,
337                                culvert_velocity, outlet_culvert_depth))
338
339            culv_froude=math.sqrt(Q**2*flow_width/(anuga.g*flow_area**3))
340            if local_debug =='true':
341                anuga.log.critical('FLOW AREA = %s' % str(flow_area))
342                anuga.log.critical('PERIMETER = %s' % str(perimeter))
343                anuga.log.critical('Q final = %s' % str(Q))
344                anuga.log.critical('FROUDE = %s' % str(culv_froude))
345
346            # Determine momentum at the outlet
347            barrel_velocity = Q/(flow_area + anuga.velocity_protection/flow_area)
348
349        else: # self.inflow.get_average_height() < 0.01:
350            Q = barrel_velocity = outlet_culvert_depth = 0.0
351
352        # Temporary flow limit
353        if barrel_velocity > self.max_velocity:
354            barrel_velocity = self.max_velocity
355            Q = flow_area * barrel_velocity
356
357        return Q, barrel_velocity, outlet_culvert_depth
358
359       
360       
361       
Note: See TracBrowser for help on using the repository browser.