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

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

Deleting unnecessary structure test scripts

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