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

Last change on this file since 7998 was 7998, checked in by habili, 14 years ago

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