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

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

Consolidation of the culvert routines into structure classes

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