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

Last change on this file since 8018 was 8018, checked in by steve, 14 years ago

Added logging to file for fractional step operators. There needs to an member function log_timestepping_statistics() implemented for each operator.

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