Changeset 5437 for anuga_core

Timestamp:

Jun 25, 2008, 7:10:21 PM (16 years ago)

Author:

ole

Message:

Finished culvert flow work with a good working solution

Location:

anuga_core/source/anuga/culvert_flows

Files:

• culvert_class.py (modified) (6 diffs)
• culvert_routines.py (modified) (8 diffs)

anuga_core/source/anuga/culvert_flows/culvert_class.py

@@ -94 +94 @@
         if blockage_bottup is None: blockage_bottup=0.00
         if culvert_routine is None: culvert_routine=boyd_generalised_culvert_model
-        if label is None: label = 'culvert_flow_' + id(self) 
+        if label is None: label = 'culvert_flow'
+        label += '_' + str(id(self)) 
         
         # Open log file for storing some specific results...
@@ -100 +101 @@
         self.label = label
 
-
-        # Create the fundamental culvert polygons from POLYGON 
+        # Print something
+        log_to_file(self.log_filename, self.label)        
+        log_to_file(self.log_filename, description)
+        log_to_file(self.log_filename, self.culvert_type)        
+
+
+        # Create the fundamental culvert polygons from POLYGON
+        if self.culvert_type == 'circle':
+            # Redefine width and height for use with create_culvert_polygons
+            width = height = diameter
+        
         P = create_culvert_polygons(end_point0,
                                     end_point1,
@@ -135 +145 @@
         self.verbose = verbose
         self.last_time = 0.0        
-        self.temp_keep_delta_t = 0.0
         
 
@@ -239 +248 @@
                 elevation = dq['elevation'].get_values(location='centroids', indices=opening.exchange_indices)
                 z = mean(elevation)                   # Average elevation
-                
+
             # Estimated depth above the culvert inlet
-            d = w - z
-
+            d = w - z  # Used for calculations involving depth
             if d < 0.0:
                 # This is possible since w and z are taken at different locations
@@ -249 +257 @@
                 d = 0.0
             
+
+
+            # Depth at exchange area used to trigger calculations
+            stage = dq['stage'].get_values(location='centroids', indices=enquiry_indices)
+            elevation = dq['elevation'].get_values(location='centroids', indices=enquiry_indices)
+            depth = stage - elevation
+            d_trigger = mean(depth)
+
+
+
             # Ratio of depth to culvert height.
             # If ratio > 1 then culvert is running full
-            ratio = d/self.height  
+            if self.culvert_type == 'circle':
+                ratio = d/self.diameter
+            else:    
+                ratio = d/self.height  
             opening.ratio = ratio
                 
@@ -260 +281 @@
             opening.specific_energy = Es            
             
-            # Store current average stage and depth with each opening object 
+            # Store current average stage and depth with each opening object
             opening.depth = d
+            opening.depth_trigger = d_trigger            
             opening.stage = w
             opening.elevation = z

anuga_core/source/anuga/culvert_flows/culvert_routines.py

@@ -58 +58 @@
     sum_loss = culvert.sum_loss
     length = culvert.length
-    
-    if inlet.depth >= 0.01:
+
+    if inlet.depth_trigger >= 0.01 and inlet.depth >= 0.01:
         # Calculate driving energy
         E = inlet.specific_energy
@@ -91 +91 @@
                 outlet_culvert_depth = inlet.depth
                 width = diameter*sin(alpha)
-                perimeter = alpha*diameter                
+                #perimeter = alpha*diameter                
                 case = 'Inlet unsubmerged'
             else:    
@@ -98 +98 @@
                 outlet_culvert_depth = diameter
                 width = diameter
-                perimeter = diameter
+                #perimeter = diameter
                 case = 'Inlet submerged'                    
 
@@ -129 +129 @@
                     case = 'Outlet is open channel flow'
 
+                hyd_rad = flow_area/perimeter
+                s = 'hydraulic radius at outlet = %f' %hyd_rad
+                log_to_file(log_filename, s)
+
+                # Outlet control velocity using tail water
+                culvert_velocity = sqrt(delta_Et/((sum_loss/2*g)+(manning**2*length)/hyd_rad**1.33333)) 
+                Q_outlet_tailwater = flow_area * culvert_velocity
+
+                s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater
+                log_to_file(log_filename, s)
+                Q = min(Q, Q_outlet_tailwater)
+                    
+
 
         else:
@@ -148 +161 @@
                 flow_area = width*inlet.depth
                 outlet_culvert_depth = inlet.depth
-                perimeter=(width+2.0*inlet.depth)                
+                #perimeter=(width+2.0*inlet.depth)                
                 case = 'Inlet unsubmerged'
             else:    
@@ -154 +167 @@
                 flow_area = width*height
                 outlet_culvert_depth = height
-                perimeter=2.0*(width+height)                
+                #perimeter=2.0*(width+height)                
                 case = 'Inlet submerged'                    
 
@@ -176 +189 @@
                     perimeter=(width+2.0*outlet.depth)
                     case = 'Outlet is open channel flow'
-                    
-
-
-
-        # Common code for rectangular and circular culvert types    
-        hyd_rad = flow_area/perimeter
-        s = 'hydraulic radius at outlet = %f' %hyd_rad
-        log_to_file(log_filename, s)
-
-        # Outlet control velocity using tail water
-        culvert_velocity = sqrt(delta_Et/((sum_loss/2*g)+(manning**2*length)/hyd_rad**1.33333)) 
-        Q_outlet_tailwater = flow_area * culvert_velocity
-
-        s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater
-        log_to_file(log_filename, s)
-        Q = min(Q, Q_outlet_tailwater)
-
+
+                hyd_rad = flow_area/perimeter
+                s = 'hydraulic radius at outlet = %f' %hyd_rad
+                log_to_file(log_filename, s)
+
+                # Outlet control velocity using tail water
+                culvert_velocity = sqrt(delta_Et/((sum_loss/2*g)+(manning**2*length)/hyd_rad**1.33333)) 
+                Q_outlet_tailwater = flow_area * culvert_velocity
+
+                s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater
+                log_to_file(log_filename, s)
+                Q = min(Q, Q_outlet_tailwater)
+
+
+        # Common code for circle and square geometries
         log_to_file(log_filename, 'Case: "%s"' %case)
-    
         flow_rate_control=Q
 
@@ -209 +219 @@
         # Determine momentum at the outlet 
         barrel_velocity = Q/(flow_area + velocity_protection/flow_area)
+
+
     else: #inlet.depth < 0.01:
         Q = barrel_velocity = outlet_culvert_depth = 0.0