source: trunk/anuga_core/source/anuga/structures/boyd_box_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.

File size: 9.4 KB
RevLine 
[8001]1import anuga
[7993]2import math
[7997]3import types
[7993]4
[8001]5class Boyd_box_operator(anuga.Structure_operator):
[7993]6    """Culvert flow - transfer water from one rectangular box to another.
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 (either diameter or width, height),
13    mannings_rougness,
[8018]14    """
[7993]15
[8018]16
[7993]17    def __init__(self,
18                 domain,
19                 end_point0, 
20                 end_point1,
[7997]21                 losses,
[7993]22                 width,
23                 height=None,
24                 apron=None,
25                 manning=0.013,
26                 enquiry_gap=0.2,
27                 use_momentum_jet=True,
28                 use_velocity_head=True,
[7996]29                 description=None,
[8018]30                 label=None,
31                 structure_type='boyd_box',
32                 logging=False,
[7993]33                 verbose=False):
[8018]34
[7993]35                     
[8001]36        anuga.Structure_operator.__init__(self,
37                                          domain,
38                                          end_point0, 
39                                          end_point1,
40                                          width,
41                                          height,
42                                          apron,
43                                          manning,
44                                          enquiry_gap,                                                       
45                                          description,
[8018]46                                          label,
47                                          structure_type,
48                                          logging,
49                                          verbose)     
[7993]50       
[7997]51       
52        if type(losses) == types.DictType:
53            self.sum_loss = sum(losses.values())
54        elif type(losses) == types.ListType:
[8001]55            self.sum_loss = sum(losses)
[7997]56        else:
[8001]57            self.sum_loss = losses
[7997]58       
[7993]59        self.use_momentum_jet = use_momentum_jet
60        self.use_velocity_head = use_velocity_head
61       
62        self.culvert_length = self.get_culvert_length()
63        self.culvert_width = self.get_culvert_width()
64        self.culvert_height = self.get_culvert_height()
65
66        self.max_velocity = 10.0
[7995]67
[7993]68        self.inlets = self.get_inlets()
[7995]69
70
71        # Stats
[7993]72       
[7995]73        self.discharge = 0.0
74        self.velocity = 0.0
[7993]75
76   
[8008]77    def discharge_routine(self):
[7993]78
79        local_debug ='false'
80       
81        if self.inflow.get_enquiry_height() > 0.01: #this value was 0.01:
82            if local_debug =='true':
[8001]83                anuga.log.critical('Specific E & Deltat Tot E = %s, %s'
[7993]84                             % (str(self.inflow.get_enquiry_specific_energy()),
85                                str(self.delta_total_energy)))
[8001]86                anuga.log.critical('culvert type = %s' % str(culvert_type))
[7993]87            # Water has risen above inlet
88
89            if self.log_filename is not None:
90                s = 'Specific energy  = %f m' % self.inflow.get_enquiry_specific_energy()
91                log_to_file(self.log_filename, s)
92
93            msg = 'Specific energy at inlet is negative'
94            assert self.inflow.get_enquiry_specific_energy() >= 0.0, msg
95
96            if self.use_velocity_head :
[7996]97                self.driving_energy = self.inflow.get_enquiry_specific_energy()
[7993]98            else:
[7996]99                self.driving_energy = self.inflow.get_enquiry_height()
[7993]100
101            height = self.culvert_height
102            width = self.culvert_width
103            flow_width = self.culvert_width
[7996]104            # intially assume the culvert flow is controlled by the inlet
105            # check unsubmerged and submerged condition and use Min Q
106            # but ensure the correct flow area and wetted perimeter are used
[8001]107            Q_inlet_unsubmerged = 0.544*anuga.g**0.5*width*self.driving_energy**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
108            Q_inlet_submerged = 0.702*anuga.g**0.5*width*height**0.89*self.driving_energy**0.61  # Flow based on Inlet Ctrl Inlet Submerged
[7993]109
110            # FIXME(Ole): Are these functions really for inlet control?
111            if Q_inlet_unsubmerged < Q_inlet_submerged:
112                Q = Q_inlet_unsubmerged
[8001]113                dcrit = (Q**2/anuga.g/width**2)**0.333333
[7993]114                if dcrit > height:
115                    dcrit = height
[7996]116                    flow_area = width*dcrit
117                    perimeter= 2.0*(width+dcrit)
118                else: # dcrit < height
119                    flow_area = width*dcrit
120                    perimeter= 2.0*dcrit+width
[7993]121                outlet_culvert_depth = dcrit
122                case = 'Inlet unsubmerged Box Acts as Weir'
[7996]123            else: # Inlet Submerged but check internal culvert flow depth
[7993]124                Q = Q_inlet_submerged
[8001]125                dcrit = (Q**2/anuga.g/width**2)**0.333333
[7996]126                if dcrit > height:
127                    dcrit = height
128                    flow_area = width*dcrit
129                    perimeter= 2.0*(width+dcrit)
130                else: # dcrit < height
131                    flow_area = width*dcrit
132                    perimeter= 2.0*dcrit+width
133                outlet_culvert_depth = dcrit
[7993]134                case = 'Inlet submerged Box Acts as Orifice'
135
[8001]136            dcrit = (Q**2/anuga.g/width**2)**0.333333
[7996]137            # May not need this .... check if same is done above
[7993]138            outlet_culvert_depth = dcrit
139            if outlet_culvert_depth > height:
140                outlet_culvert_depth = height  # Once again the pipe is flowing full not partfull
141                flow_area = width*height  # Cross sectional area of flow in the culvert
142                perimeter = 2*(width+height)
143                case = 'Inlet CTRL Outlet unsubmerged PIPE PART FULL'
144            else:
145                flow_area = width * outlet_culvert_depth
146                perimeter = width+2*outlet_culvert_depth
147                case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'
[7996]148            # Initial Estimate of Flow for Outlet Control using energy slope
149            #( may need to include Culvert Bed Slope Comparison)
150            hyd_rad = flow_area/perimeter
[8001]151            culvert_velocity = math.sqrt(self.delta_total_energy/((self.sum_loss/2/anuga.g)+(self.manning**2*self.culvert_length)/hyd_rad**1.33333))
[7996]152            Q_outlet_tailwater = flow_area * culvert_velocity
153           
154           
155            if self.delta_total_energy < self.driving_energy:
[7993]156                # Calculate flows for outlet control
157
158                # Determine the depth at the outlet relative to the depth of flow in the Culvert
159                if self.outflow.get_enquiry_height() > height:        # The Outlet is Submerged
160                    outlet_culvert_depth=height
161                    flow_area=width*height       # Cross sectional area of flow in the culvert
162                    perimeter=2.0*(width+height)
163                    case = 'Outlet submerged'
164                else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
[8001]165                    dcrit = (Q**2/anuga.g/width**2)**0.333333
[7993]166                    outlet_culvert_depth=dcrit   # For purpose of calculation assume the outlet depth = Critical Depth
167                    if outlet_culvert_depth > height:
168                        outlet_culvert_depth=height
169                        flow_area=width*height
170                        perimeter=2.0*(width+height)
171                        case = 'Outlet is Flowing Full'
172                    else:
173                        flow_area=width*outlet_culvert_depth
174                        perimeter=(width+2.0*outlet_culvert_depth)
175                        case = 'Outlet is open channel flow'
176
177                hyd_rad = flow_area/perimeter
178
179                if self.log_filename is not None:
180                    s = 'hydraulic radius at outlet = %f' % hyd_rad
181                    log_to_file(self.log_filename, s)
182
[7996]183                # Final Outlet control velocity using tail water
[8001]184                culvert_velocity = math.sqrt(self.delta_total_energy/((self.sum_loss/2/anuga.g)+(self.manning**2*self.culvert_length)/hyd_rad**1.33333))
[7993]185                Q_outlet_tailwater = flow_area * culvert_velocity
186
187                if self.log_filename is not None:
188                    s = 'Q_outlet_tailwater = %.6f' % Q_outlet_tailwater
189                    log_to_file(self.log_filename, s)
190                Q = min(Q, Q_outlet_tailwater)
191            else:
192                pass
193                #FIXME(Ole): What about inlet control?
194
[8001]195            culv_froude=math.sqrt(Q**2*flow_width/(anuga.g*flow_area**3))
[7993]196            if local_debug =='true':
[8001]197                anuga.log.critical('FLOW AREA = %s' % str(flow_area))
198                anuga.log.critical('PERIMETER = %s' % str(perimeter))
199                anuga.log.critical('Q final = %s' % str(Q))
200                anuga.log.critical('FROUDE = %s' % str(culv_froude))
[7993]201
202            # Determine momentum at the outlet
[8001]203            barrel_velocity = Q/(flow_area + anuga.velocity_protection/flow_area)
[7993]204
205        # END CODE BLOCK for DEPTH  > Required depth for CULVERT Flow
206
207        else: # self.inflow.get_enquiry_height() < 0.01:
208            Q = barrel_velocity = outlet_culvert_depth = 0.0
209
210        # Temporary flow limit
211        if barrel_velocity > self.max_velocity:
212            barrel_velocity = self.max_velocity
213            Q = flow_area * barrel_velocity
214
215        return Q, barrel_velocity, outlet_culvert_depth
216       
217       
218       
Note: See TracBrowser for help on using the repository browser.