source: trunk/anuga_core/source/anuga/structures/boyd_box_operator.py @ 7995

Last change on this file since 7995 was 7995, checked in by steve, 12 years ago

Added stats

File size: 11.2 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
5
6import structure_operator
7
8class Boyd_box_operator(structure_operator.Structure_operator):
9    """Culvert flow - transfer water from one rectangular box to another.
10    Sets up the geometry of problem
11   
12    This is the base class for culverts. Inherit from this class (and overwrite
13    compute_discharge method for specific subclasses)
14   
15    Input: Two points, pipe_size (either diameter or width, height),
16    mannings_rougness,
17    """ 
18
19    def __init__(self,
20                 domain,
21                 end_point0, 
22                 end_point1,
23                 width,
24                 height=None,
25                 apron=None,
26                 manning=0.013,
27                 enquiry_gap=0.2,
28                 use_momentum_jet=True,
29                 use_velocity_head=True,
30                 verbose=False):
31                     
32        structure_operator.Structure_operator.__init__(self,
33                                                       domain,
34                                                       end_point0, 
35                                                       end_point1,
36                                                       width,
37                                                       height,
38                                                       apron,
39                                                       manning,
40                                                       enquiry_gap,
41                                                       verbose)           
42       
43        self.use_momentum_jet = use_momentum_jet
44        self.use_velocity_head = use_velocity_head
45       
46        self.culvert_length = self.get_culvert_length()
47        self.culvert_width = self.get_culvert_width()
48        self.culvert_height = self.get_culvert_height()
49
50
51        self.sum_loss = 0.0
52        self.max_velocity = 10.0
53        self.log_filename = None
54
55        self.inlets = self.get_inlets()
56
57
58        # Stats
59       
60        self.discharge = 0.0
61        self.velocity = 0.0
62       
63       
64    def __call__(self):
65       
66        timestep = self.domain.get_timestep()
67       
68        self.__determine_inflow_outflow()
69       
70        Q, barrel_speed, outlet_depth = self.__discharge_routine()
71
72        #inflow  = self.routine.get_inflow()
73        #outflow = self.routine.get_outflow()
74
75        old_inflow_height = self.inflow.get_average_height()
76        old_inflow_xmom = self.inflow.get_average_xmom()
77        old_inflow_ymom = self.inflow.get_average_ymom()
78           
79        if old_inflow_height > 0.0 :
80                Qstar = Q/old_inflow_height
81        else:
82                Qstar = 0.0
83
84        factor = 1.0/(1.0 + Qstar*timestep/self.inflow.get_area())
85
86        new_inflow_height = old_inflow_height*factor
87        new_inflow_xmom = old_inflow_xmom*factor
88        new_inflow_ymom = old_inflow_ymom*factor
89           
90
91        self.inflow.set_heights(new_inflow_height)
92
93        #inflow.set_xmoms(Q/inflow.get_area())
94        #inflow.set_ymoms(0.0)
95
96
97        self.inflow.set_xmoms(new_inflow_xmom)
98        self.inflow.set_ymoms(new_inflow_ymom)
99
100
101        loss = (old_inflow_height - new_inflow_height)*self.inflow.get_area()
102
103           
104        # set outflow
105        if old_inflow_height > 0.0 :
106                timestep_star = timestep*new_inflow_height/old_inflow_height
107        else:
108            timestep_star = 0.0
109
110           
111        outflow_extra_height = Q*timestep_star/self.outflow.get_area()
112        outflow_direction = - self.outflow.outward_culvert_vector
113        outflow_extra_momentum = outflow_extra_height*barrel_speed*outflow_direction
114           
115
116        gain = outflow_extra_height*self.outflow.get_area()
117           
118        #print Q, Q*timestep, barrel_speed, outlet_depth, Qstar, factor, timestep_star
119        #print '  ', loss, gain
120
121        # Stats
122        self.discharge  = outflow_extra_height*self.outflow.get_area()/timestep
123        self.velocity = self.discharge/outlet_depth/self.width
124
125        new_outflow_height = self.outflow.get_average_height() + outflow_extra_height
126
127        if self.use_momentum_jet :
128            # FIXME (SR) Review momentum to account for possible hydraulic jumps at outlet
129            #new_outflow_xmom = outflow.get_average_xmom() + outflow_extra_momentum[0]
130            #new_outflow_ymom = outflow.get_average_ymom() + outflow_extra_momentum[1]
131
132            new_outflow_xmom = barrel_speed*new_outflow_height*outflow_direction[0]
133            new_outflow_ymom = barrel_speed*new_outflow_height*outflow_direction[1]
134
135        else:
136            #new_outflow_xmom = outflow.get_average_xmom()
137            #new_outflow_ymom = outflow.get_average_ymom()
138
139            new_outflow_xmom = 0.0
140            new_outflow_ymom = 0.0
141
142
143        self.outflow.set_heights(new_outflow_height)
144        self.outflow.set_xmoms(new_outflow_xmom)
145        self.outflow.set_ymoms(new_outflow_ymom)
146
147
148    def __determine_inflow_outflow(self):
149        # Determine flow direction based on total energy difference
150
151        if self.use_velocity_head:
152            self.delta_total_energy = self.inlets[0].get_enquiry_total_energy() - self.inlets[1].get_enquiry_total_energy()
153        else:
154            self.delta_total_energy = self.inlets[0].get_enquiry_stage() - self.inlets[1].get_enquiry_stage()
155
156
157        self.inflow  = self.inlets[0]
158        self.outflow = self.inlets[1]
159       
160
161        if self.delta_total_energy < 0:
162            self.inflow  = self.inlets[1]
163            self.outflow = self.inlets[0]
164            self.delta_total_energy = -self.delta_total_energy
165
166   
167    def __discharge_routine(self):
168
169        local_debug ='false'
170       
171        if self.inflow.get_enquiry_height() > 0.01: #this value was 0.01:
172            if local_debug =='true':
173                log.critical('Specific E & Deltat Tot E = %s, %s'
174                             % (str(self.inflow.get_enquiry_specific_energy()),
175                                str(self.delta_total_energy)))
176                log.critical('culvert type = %s' % str(culvert_type))
177            # Water has risen above inlet
178
179            if self.log_filename is not None:
180                s = 'Specific energy  = %f m' % self.inflow.get_enquiry_specific_energy()
181                log_to_file(self.log_filename, s)
182
183            msg = 'Specific energy at inlet is negative'
184            assert self.inflow.get_enquiry_specific_energy() >= 0.0, msg
185
186            if self.use_velocity_head :
187                driving_energy = self.inflow.get_enquiry_specific_energy()
188            else:
189                driving_energy = self.inflow.get_enquiry_height()
190
191            height = self.culvert_height
192            width = self.culvert_width
193            flow_width = self.culvert_width
194
195            Q_inlet_unsubmerged = 0.540*g**0.5*width*driving_energy**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
196            Q_inlet_submerged = 0.702*g**0.5*width*height**0.89*driving_energy**0.61  # Flow based on Inlet Ctrl Inlet Submerged
197
198            # FIXME(Ole): Are these functions really for inlet control?
199            if Q_inlet_unsubmerged < Q_inlet_submerged:
200                Q = Q_inlet_unsubmerged
201                dcrit = (Q**2/g/width**2)**0.333333
202                if dcrit > height:
203                    dcrit = height
204                flow_area = width*dcrit
205                outlet_culvert_depth = dcrit
206                case = 'Inlet unsubmerged Box Acts as Weir'
207            else:
208                Q = Q_inlet_submerged
209                flow_area = width*height
210                outlet_culvert_depth = height
211                case = 'Inlet submerged Box Acts as Orifice'
212
213            dcrit = (Q**2/g/width**2)**0.333333
214
215            outlet_culvert_depth = dcrit
216            if outlet_culvert_depth > height:
217                outlet_culvert_depth = height  # Once again the pipe is flowing full not partfull
218                flow_area = width*height  # Cross sectional area of flow in the culvert
219                perimeter = 2*(width+height)
220                case = 'Inlet CTRL Outlet unsubmerged PIPE PART FULL'
221            else:
222                flow_area = width * outlet_culvert_depth
223                perimeter = width+2*outlet_culvert_depth
224                case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'
225
226            if self.delta_total_energy < driving_energy:
227                # Calculate flows for outlet control
228
229                # Determine the depth at the outlet relative to the depth of flow in the Culvert
230                if self.outflow.get_enquiry_height() > height:        # The Outlet is Submerged
231                    outlet_culvert_depth=height
232                    flow_area=width*height       # Cross sectional area of flow in the culvert
233                    perimeter=2.0*(width+height)
234                    case = 'Outlet submerged'
235                else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
236                    dcrit = (Q**2/g/width**2)**0.333333
237                    outlet_culvert_depth=dcrit   # For purpose of calculation assume the outlet depth = Critical Depth
238                    if outlet_culvert_depth > height:
239                        outlet_culvert_depth=height
240                        flow_area=width*height
241                        perimeter=2.0*(width+height)
242                        case = 'Outlet is Flowing Full'
243                    else:
244                        flow_area=width*outlet_culvert_depth
245                        perimeter=(width+2.0*outlet_culvert_depth)
246                        case = 'Outlet is open channel flow'
247
248                hyd_rad = flow_area/perimeter
249
250                if self.log_filename is not None:
251                    s = 'hydraulic radius at outlet = %f' % hyd_rad
252                    log_to_file(self.log_filename, s)
253
254                # Outlet control velocity using tail water
255                culvert_velocity = math.sqrt(self.delta_total_energy/((self.sum_loss/2/g)+(self.manning**2*self.culvert_length)/hyd_rad**1.33333))
256                Q_outlet_tailwater = flow_area * culvert_velocity
257
258                if self.log_filename is not None:
259                    s = 'Q_outlet_tailwater = %.6f' % Q_outlet_tailwater
260                    log_to_file(self.log_filename, s)
261                Q = min(Q, Q_outlet_tailwater)
262            else:
263                pass
264                #FIXME(Ole): What about inlet control?
265
266            culv_froude=math.sqrt(Q**2*flow_width/(g*flow_area**3))
267            if local_debug =='true':
268                log.critical('FLOW AREA = %s' % str(flow_area))
269                log.critical('PERIMETER = %s' % str(perimeter))
270                log.critical('Q final = %s' % str(Q))
271                log.critical('FROUDE = %s' % str(culv_froude))
272
273            # Determine momentum at the outlet
274            barrel_velocity = Q/(flow_area + velocity_protection/flow_area)
275
276        # END CODE BLOCK for DEPTH  > Required depth for CULVERT Flow
277
278        else: # self.inflow.get_enquiry_height() < 0.01:
279            Q = barrel_velocity = outlet_culvert_depth = 0.0
280
281        # Temporary flow limit
282        if barrel_velocity > self.max_velocity:
283            barrel_velocity = self.max_velocity
284            Q = flow_area * barrel_velocity
285
286        return Q, barrel_velocity, outlet_culvert_depth
287       
288       
289       
Note: See TracBrowser for help on using the repository browser.