source: trunk/anuga_core/source/anuga/structures/boyd_box_routine.py @ 7969

Last change on this file since 7969 was 7969, checked in by steve, 14 years ago

Playing with boyd function
File size: 6.1 KB
Line 
1#! /usr/bin/python
2
3# To change this template, choose Tools | Templates
4# and open the template in the editor.
5
6__author__="steve"
7__date__ ="$23/08/2010 5:18:51 PM$"
8
9
10
11def boyd_box(height, width, flow_width, inflow_specific_energy):
12    """Boyd's generalisation of the US department of transportation culvert methods
13
14        WARNING THIS IS A SIMPLISTIC APPROACH and OUTLET VELOCITIES ARE LIMITED TO EITHER
15        FULL PIPE OR CRITICAL DEPTH ONLY
16        For Supercritical flow this is UNDERESTIMATING the Outlet Velocity
17        The obtain the CORRECT velocity requires an iteration of Depth to Establish
18        the Normal Depth of flow in the pipe.
19
20        It is proposed to provide this in a seperate routine called
21    boyd_generalised_culvert_model_complex
22
23    The Boyd Method is based on methods described by the following:
24        1.
25        US Dept. Transportation Federal Highway Administration (1965)
26        Hydraulic Chart for Selection of Highway Culverts.
27        Hydraulic Engineering Circular No. 5 US Government Printing
28        2.
29        US Dept. Transportation Federal Highway Administration (1972)
30        Capacity charts for the Hydraulic design of highway culverts.
31        Hydraulic Engineering Circular No. 10 US Government Printing
32    These documents provide around 60 charts for various configurations of culverts and inlets.
33
34        Note these documents have been superceded by:
35        2005 Hydraulic Design of Highway Culverts, Hydraulic Design Series No. 5 (HDS-5),
36        Which combines culvert design information previously contained in Hydraulic Engineering Circulars
37        (HEC) No. 5, No. 10, and No. 13 with hydrologic, storage routing, and special culvert design information.
38        HEC-5 provides 20 Charts
39        HEC-10 Provides an additional 36 Charts
40        HEC-13 Discusses the Design of improved more efficient inlets
41        HDS-5 Provides 60 sets of Charts
42
43        In 1985 Professor Michael Boyd Published "Head-Discharge Relations for Culverts", and in
44        1987 published "Generalised Head Discharge Equations for Culverts".
45        These papers reviewed the previous work by the US DOT and provided a simplistic approach for 3 configurations.
46
47        It may be possible to extend the same approach for additional charts in the original work, but to date this has not been done.
48        The additional charts cover a range of culvert shapes and inlet configurations
49
50    """
51
52    # Calculate flows for inflow control
53
54    Q_inflow_unsubmerged = 0.540*g**0.5*width*inflow_specific_energy**1.50 # Flow based on inflow Ctrl inflow Unsubmerged
55    Q_inflow_submerged = 0.702*g**0.5*width*height**0.89*inflow_specific_energy**0.61  # Flow based on inflow Ctrl inflow Submerged
56
57    if log_filename is not None:
58        s = 'Q_inflow_unsubmerged = %.6f, Q_inflow_submerged = %.6f' %(Q_inflow_unsubmerged, Q_inflow_submerged)
59        log_to_file(log_filename, s)
60
61    # FIXME(Ole): Are these functions really for inflow control?
62    if Q_inflow_unsubmerged < Q_inflow_submerged:
63        Q = Q_inflow_unsubmerged
64        dcrit = (Q**2/g/width**2)**0.333333
65        if dcrit > height:
66            dcrit = height
67        flow_area = width*dcrit
68        outflow_culvert_depth = dcrit
69        case = 'inflow unsubmerged Box Acts as Weir'
70    else:
71        Q = Q_inflow_submerged
72        flow_area = width*height
73        outflow_culvert_depth = height
74        case = 'inflow submerged Box Acts as Orifice'
75
76    dcrit = (Q**2/g/width**2)**0.333333
77
78    outflow_culvert_depth = dcrit
79    if outflow_culvert_depth > height:
80        outflow_culvert_depth = height  # Once again the pipe is flowing full not partfull
81        flow_area = width*height  # Cross sectional area of flow in the culvert
82        perimeter = 2*(width+height)
83        case = 'inflow CTRL outflow unsubmerged PIPE PART FULL'
84    else:
85        flow_area = width * outflow_culvert_depth
86        perimeter = width+2*outflow_culvert_depth
87        case = 'inflow CTRL Culvert is open channel flow we will for now assume critical depth'
88
89    if delta_total_energy < inflow_specific_energy:
90        # Calculate flows for outflow control
91
92        # Determine the depth at the outflow relative to the depth of flow in the Culvert
93        if outflow_depth > height:        # The outflow is Submerged
94            outflow_culvert_depth=height
95            flow_area=width*height       # Cross sectional area of flow in the culvert
96            perimeter=2.0*(width+height)
97            case = 'outflow submerged'
98        else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
99            dcrit = (Q**2/g/width**2)**0.333333
100            outflow_culvert_depth=dcrit   # For purpose of calculation assume the outflow depth = Critical Depth
101            if outflow_culvert_depth > height:
102                outflow_culvert_depth=height
103                flow_area=width*height
104                perimeter=2.0*(width+height)
105                case = 'outflow is Flowing Full'
106            else:
107                flow_area=width*outflow_culvert_depth
108                perimeter=(width+2.0*outflow_culvert_depth)
109                case = 'outflow is open channel flow'
110
111        hyd_rad = flow_area/perimeter
112
113        if log_filename is not None:
114            s = 'hydraulic radius at outflow = %f' % hyd_rad
115            log_to_file(log_filename, s)
116
117        # outflow control velocity using tail water
118        culvert_velocity = sqrt(delta_total_energy/((sum_loss/2/g)+(manning**2*culvert_length)/hyd_rad**1.33333))
119        Q_outflow_tailwater = flow_area * culvert_velocity
120
121        if log_filename is not None:
122            s = 'Q_outflow_tailwater = %.6f' % Q_outflow_tailwater
123            log_to_file(log_filename, s)
124        Q = min(Q, Q_outflow_tailwater)
125
126    return Q
127
128
129if __name__ == "__main__":
130
131
132    g=9.81
133    culvert_slope=0.1  # Downward
134
135    inlet_depth=2.0
136    outlet_depth=0.0
137
138    inlet_velocity=0.0,
139    outlet_velocity=0.0,
140
141    culvert_length=4.0
142    culvert_width=1.2
143    culvert_height=0.75
144
145    culvert_type='box'
146    manning=0.013
147    sum_loss=0.0
148
149    inlet_specific_energy=inlet_depth #+0.5*v**2/g
150    z_in = 0.0
151    z_out = -culvert_length*culvert_slope/100
152    E_in = z_in+inlet_depth # +
153    E_out = z_out+outlet_depth # +
154    delta_total_energy = E_in-E_out
155
156    Q = boyd_box(culvert_height, culvert_width, culvert_width, inlet_specific_energy)
157
158    print 'Q ',Q
Note: See TracBrowser for help on using the repository browser.