damyarou

python, GMT などのプログラム

Hydraulic jump type energy dissipator


\begin{equation}
v_1=0.9\cdot \sqrt{2\cdot g\cdot H}
\end{equation}

\begin{equation}
d_1=\cfrac{Q}{b\cdot v_1}
\end{equation}

\begin{equation}
F_1=\cfrac{v_1}{\sqrt{g\cdot d_1}}
\end{equation}

\begin{equation}
\cfrac{d_2}{d_1}=\cfrac{1}{2}\cdot \left(\sqrt{1+8\cdot F_1{}^2}-1\right)
\end{equation}

\begin{equation}
\cfrac{W}{d_1}=\cfrac{(1+2\cdot F_1{}^2)\cdot \sqrt{1+8\cdot F_1{}^2}-1-5\cdot F_1{}^2}
{1+4\cdot F_1{}^2-\sqrt{1+8\cdot F_1{}^2}}
-\left(\cfrac{\sqrt{g}}{C}\cdot F_1\right)^{2/3}
\end{equation}

\begin{equation}
L=4.5\cdot d_2
\end{equation}

\begin{align}
& Q   & & \text{discharge}\\
& b   & & \text{width of apron}\\
& H   & & \text{height from apron level to flood water level of reservoir}\\
& F_1 & & \text{Froude number at starting point of jump}\\
& v_1 & & \text{velocity at starting point of jump}\\
& d_1 & & \text{water depth at staring point of jump}\\
& d_2 & & \text{water depth at end point of jump}\\
& W   & & \text{height of end sill}\\
& C   & & \text{discharge coefficient of end sill}\\
& L   & & \text{length of energy dissipator}\\
& g   & & \text{gravity acceleration}
\end{align}
import numpy as np

def main():
    g=9.8
    q=42.0
    h0=554.1-321.0
    v1=26.5
    b=4.0

    d1=np.round(q/b/v1,decimals=3)
    fr=np.round(v1/np.sqrt(g*d1),decimals=3)
    d1d2=np.round(1/2*(np.sqrt(1+8*fr**2)-1),decimals=3)
    d2=np.round(d1*d1d2,decimals=3)
    print('v1=',v1)
    print('d1=',d1)
    print('Fr=',fr)
    print('d1d2=',d1d2)
    print('d2=',d2)

    c=1.704
    c1=(1+2*fr**2)*np.sqrt(1+8*fr**2)-1-5*fr**2
    c2=1+4*fr**2-np.sqrt(1+8*fr**2)
    cc=np.round(c1/c2-(np.sqrt(g)/c*fr)**(2/3),decimals=3)
    ww=np.round(d1*cc,decimals=3)
    print('cc=',cc)
    print('W=',ww)
    print('L=',4.5*d2)


#==============
# Execution
#==============
if __name__ == '__main__': main()

Calculation results are shown below.

v1= 26.5
d1= 0.396
Fr= 13.452
d1d2= 18.531
d2= 7.338
cc= 10.31
W= 4.083
L= 33.021

Thank you.

Normal depth of steep channel


\begin{gather}
Q=A\cdot v \\
v=\cfrac{1}{n}\cdot R^{2/3}\cdot (\sin\theta)^{1/2} \\
A=b\cdot h \qquad R=\cfrac{b\cdot h}{b+2\cdot h}
\end{gather}
Q(known) discharge
b(known) width of rectangular open channe
n(known) Manning's roughness coefficient
\theta(known) gradient of channel invert
h(unknown) normal depth

The critical depth of steep channel is calculated as following equation.


\begin{equation}
h_c=\left(\cfrac{Q^2}{g\cdot b^2\cdot \cos\theta}\right)^{1/3}
\end{equation}
# normal depth and critical depth of rectangular cross section
import numpy as np
from scipy import optimize


def cal_hc(q,b,cs):
    # critical depth
    g=9.8
    hc=(q**2/g/b**2/cs)**(1/3)
    return hc


def func(h,q,b,n,sn):
    f=q-b*h/n*(b*h/(b+2*h))**(2/3)*sn**(1/2)    
    return f    


def main():
    q=42.0  # discharge
    b=4.0   # channel width
    n=0.014 # Manning's roughness coefficient
    theta=np.radians(37.0)
    sn=np.sin(theta)
    cs=np.cos(theta)
    
    h1=0.0
    h2=10.0
    hh=optimize.brentq(func,h1,h2,args=(q,b,n,sn))
    v=q/b/hh
    
    print('hn=',hh) # normal depth
    print('v=',v)
    
    hc=cal_hc(q,b,cs)
    print('hc=',hc) # critical depth
    
#==============
# Execution
#==============
if __name__ == '__main__': main()

Calculation results are shown below.

hn= 0.3962334595851377
v= 26.499528866122652
hc= 2.415097316241126

Thank you.

Soil pressure coefficient considering earthquake by Mononobe-Okabe

Active soil pressure coefficient considering earthquake by Mononobe-Okabe can be obtained by following equation. (USACE EM 1110-2-2100 Appendix G)


\begin{equation}
P_a=K_a\cdot (1-k_v)\cdot \left\{\cfrac{1}{2}\cdot \gamma_s\cdot d^2+q\cdot d\right\}
\end{equation}

\begin{equation}
K_a=\cfrac{\cos^2 (\phi-\alpha-\theta)}{
\cos\theta\cdot \cos^2\alpha\cdot \cos(\alpha+\delta+\theta)\cdot 
\left\{
1+\sqrt{\cfrac{\sin(\phi+\delta)\cdot \sin(\phi-\beta-\theta)}{\cos(\alpha+\beta+\theta)\cdot \cos(\alpha-\beta)}}
\right\}^2
}
\end{equation}

\begin{equation}
\theta=\tan^{-1}\left(\cfrac{k_h}{1-k_v}\right)
\end{equation}

\begin{align}
& P_a      & & \text{active soil pressure}           & \quad & K_a    & & \text{active poil pressure coefficient} \\
& \gamma_s & & \text{unit weight of soil}            & \quad & d      & & \text{soil depth} \\
& q        & & \text{surcharge load}                 & \quad & \phi   & & \text{internal friction angle of soil} \\
& \alpha   & & \text{gradient of wall surface}       & \quad & \beta  & & \text{gradient of ground surface} \\
& \delta   & & \text{friction angle of wall surface} & \quad & \theta & & \text{seismic inertia angle} \\
& k_h      & & \text{horizontal seismic coefficient} & \quad & k_v    & & \text{vertical seismic coefficient}
\end{align}
import numpy as np

kh=0.077 # horizontal seismic coefficient
#kh=0 # horizontal seismic coefficient
kv=0 # vertical seismic coefficient
phi=np.radians(30.0) # internal friction angle of soil
alpha=0 # gradient of wall surface
beta=0 # gradient of ground surface
delta=0 # friction angle of wall surface
theta=np.arctan(kh/(1-kv)) # seismic inertia angle

c1=np.cos(phi-alpha-theta)
c2=np.cos(theta) * (np.cos(alpha))**2 * np.cos(alpha+delta+theta)
c3=1 + np.sqrt(np.sin(phi+delta) * np.sin(phi-beta-theta) / np.cos(alpha+beta+theta) / np.cos(alpha-beta))
Ka=c1**2/c2/c3**2

print(Ka)

Calculation results are shown below.

Ka=0.381 # earthquake (kh=0.077, kv=0)
Ka=0.333 # normal (kh=0, kv=0)

That's all. Thak you.

Python: calculation of required plate thickness of exposed type penstock against buckling pressure by Brent method

A critical buckling pressure of exposed type penstocks can be obtained by following equation.


\begin{equation}
p_k=\cfrac{2\cdot E_s}{1-\nu_s{}^2}\cdot \left(\cfrac{t}{D_0'}\right)^3
\end{equation}
p_kcritical buckling pressure of penstock
tplate thickness excluding corrosion allowance
D_0'design external diameter
E_selastic modulus of penstock (=206,000MPa)
\nu_sPoisson's modulus of penstock (=0.3)

usually, a corrosion allowance is set to \epsilon=1.5 mm, and the design plate thickness including corrosion allowance t_0 and the design external diameter D_0' can be calculated as follows.


\begin{equation}
t_0=t+\epsilon \qquad
D_0'=D_0+t_0
\end{equation}

Where, D_0 means the design internal diameter.

In case of exposed type penstocks, it is required to withstand the buckling pressure of 0.03 MPa which includes a safety factor of 1.5, and the plate thickness with this capacity shall be defined as a lower limit of plate thickness.

The equation to be solved by Brent method is shown below. The equation has a shape of f=0.


\begin{equation} 
f=p_k-\cfrac{2\cdot E_s}{1-\nu_s{}^2}\cdot \left(\cfrac{t}{D_0+t+\epsilon}\right)^3=0
\end{equation}

Required two initial values for Brent method are given as t1=1.0, t2=50.0 in the program.

Program code is shown below.

import numpy as np
from scipy import optimize


def func(t,d0,eps,pk):
    Es=206000 # elastic modulus of steel
    po=0.3    # Poisson's ratio of steel
    f=pk-2*Es/(1-po**2)*(t/(d0+t+eps))**3
    return f


def main():
    d0=4000.0 # internal diameter of penstock
    eps=1.5   # corrosion alloowance
    pk=0.03   # critical buckling pressure
    t1=1.0    # initial value for Brent method
    t2=50.0   # initial value for Brent method
    tt=optimize.brentq(func,t1,t2,args=(d0,eps,pk))
    t0=np.ceil(tt+eps) # required plate thickness
    print('t + eps=',tt+eps)
    print('t0=',t0)
   

#==============
# Execution
#==============
if __name__ == '__main__': main() 

The calculation results are shown below.

t + eps= 17.758192911648965
t0= 18.0

That's all. Thank you.

Non-uniform flow analysis by Python: calculation of water surface frofile of settling basin

A program to calculate water surface profile of settling basin which has an open channel with varied width and varied invert level. Since the condition is subcritical flow, the calculation will be carried out from downstream section to upstream section.

Basic equation of non-uniform flow analysis is shown below.


\begin{equation*}
\left(\cfrac{Q^2}{2 g A_2{}^2}+h_2+z_2\right) 
- \left(\cfrac{Q^2}{2 g A_1{}^2}+h_1+z_1\right)
= \cfrac{1}{2}\left(\cfrac{n_2{}^2 Q^2}{R_1{}^{4/3} A_1{}^2} + \cfrac{n_1{}^2 Q^2}{R_2{}^{4/3} A_2{}^2}\right)\cdot(x_2-x_1)
\end{equation*}


QDischarge (constant value)
x_2, z_2, A_2, R_2, n_2, h_2Properties of upstream section (distance, invert level, flow section area, hydraulic radius, Manning's roughness coefficient, water depth)
x_1, z_1, A_1, R_1, n_1, h_1Properties of downstream section (distance, invert level, flow section area, hydraulic radius, Manning's roughness coefficient, water depth)

Where, subscripu '1' is for downstream section and subscript '2' is for upstream section.

The equations for calculations of section area A and hydraulic radius R for rectangular cross section are shown below.


\begin{equation*}
A=b\cdot h \qquad R=\cfrac{b\cdot h}{b+2\cdot h}
\end{equation*}

where, b means channel width, h means water depth.

Since this case is a condition under subcritical flow, the conditions of downstream with subscript of '1' are all known values, and calculation is carried out to obtain the water depth of h_2 using the know conditions.

The section properties are defined as a function of x in def sec(x,h). Since the calculation model has many varied section, the part of definition of section properties (def sec(x,h)) is long. However, the core part for bi-section method is simple.

Two required initial values for bi-section method are defined as ha=3.0 and hb=7.0 in the function of def bisection(h1,x1,x2,q). These initial values should be changed appropriately, depending on the problem to be solved.

A function def sec(x,h) has a fuction to output section properties (z, A, R, n) from input data of distance (x) and water depth (h). When this function def sec(x,h) is changed, this program can be applied to a non-uniform flow analysis of open channel with various cross section.

A program is shown below.

# Non-Uniform Flow Analysis (Subcritical flow)
import numpy as np


def sec(x, h):
    # definition of section property
    # x  : distance
    # h  : water depth
    # zz : invert level
    # aa : secion area
    # rr : hydraulic radius
    # nn : Manning's roughness coefficient
    n0=0.014 # roughness coefficient (normal value)
    zz,aa,rr,nn=0,0,0,0
    if 0.0 <= x < 11.0:
        ds=11
        nn=n0
        z1=562.2; b1=4.0
        z2=560.5; b2=26.0
        zz=z1-(z1-z2)/ds*x
        bb=b1+(b2-b1)/ds*x
        aa=bb*h
        rr=bb*h/(bb+2*h)
    if 11.0 <= x < 19.0:
        ds=8.0
        nn=n0
        z1=560.5; b1=12
        z2=560.5; b2=12
        zz=z1-(z1-z2)/ds*(x-11)
        bb=b1+(b2-b1)/ds*(x-11)       
        ah=bb*h
        rh=bb*h/(bb+2*h)
        aa=ah*2
        rr=rh*2
    if 19.0 <= x < 47.0:
        ds=29.0
        nn=n0
        z1=560.5     ; b1=12.0
        z2=z1+ds*0.02; b2=12.0
        zz=z1-(z1-z2)/ds*(x-19)
        bb=b1+(b2-b1)/ds*(x-19)       
        ah=bb*h
        rh=bb*h/(bb+2*h)
        aa=ah*2
        rr=rh*2
    if 47.0 <= x < 62.0:
        ds=15.0
        nn=n0
        z1=561.06    ; b1=12.0
        z2=z1+ds*0.02; b2=6.0
        zz=z1-(z1-z2)/ds*(x-47)
        bb=b1+(b2-b1)/ds*(x-47)       
        ah=bb*h
        rh=bb*h/(bb+2*h)
        aa=ah*2
        rr=rh*2
    if 62.0 <= x < 69.0:
        ds=7.0
        nn=n0
        z1=561.36    ; b1=6.0
        z2=z1+ds*0.02; b2=6.0
        zz=z1-(z1-z2)/ds*(x-62)
        bb=b1+(b2-b1)/ds*(x-62)       
        ah=bb*h
        rh=bb*h/(bb+2*h)
        aa=ah*2
        rr=rh*2
    if 69.0 <= x < 69.5:
        ds=0.5
        nn=n0
        z1=561.5; b1=6.0
        z2=562.0; b2=6.0
        zz=z1-(z1-z2)/ds*(x-69)
        bb=b1+(b2-b1)/ds*(x-69)       
        ah=bb*h
        rh=bb*h/(bb+2*h)
        aa=ah*2
        rr=rh*2
    if 69.5 <= x < 73.5:
        ds=4.0
        nn=n0
        z1=562.0; b1=6.0
        z2=562.0; b2=6.0
        zz=z1-(z1-z2)/ds*(x-69.5)
        bb=b1+(b2-b1)/ds*(x-69.5)       
        ah=bb*h
        rh=bb*h/(bb+2*h)
        aa=ah*2
        rr=rh*2
    if 73.5 <= x < 74.0:
        ds=0.5
        nn=n0
        z1=562.0; b1=6.0
        z2=561.5; b2=6.5
        zz=z1-(z1-z2)/ds*(x-73.5)
        bb=b1+(b2-b1)/ds*(x-73.5)       
        ah=bb*h
        rh=bb*h/(bb+2*h)
        aa=ah*2
        rr=rh*2
    if 74.0 <= x < 77.0:
        ds=3.0
        nn=n0
        z1=561.5; b1=6.5
        z2=561.5; b2=6.5
        zz=z1-(z1-z2)/ds*(x-74.0)
        bb=b1+(b2-b1)/ds*(x-74.0)       
        ah=bb*h
        rh=bb*h/(bb+2*h)
        aa=ah*2
        rr=rh*2
    if 77.0 <= x < 96.625:
        ds=19.625
        nn=n0
        z1=561.5; b1=6.5
        z2=563.0; b2=6.5
        zz=z1-(z1-z2)/ds*(x-77.0)
        bb=b1+(b2-b1)/ds*(x-77.0)       
        ah=bb*h
        rh=bb*h/(bb+2*h)
        aa=ah*2
        rr=rh*2
    if 96.625 <= x <= 102.125:
        ds=5.5
        nn=n0
        z1=563.0; b1=6.5
        z2=563.0; b2=6.5
        zz=z1-(z1-z2)/ds*(x-96.625)
        bb=b1+(b2-b1)/ds*(x-96.625)       
        ah=bb*h
        rh=bb*h/(bb+2*h)
        aa=ah*2
        rr=rh*2
    return zz,aa,rr,nn


def func(h2,h1,x1,x2,q):
    g=9.8
    z1,a1,r1,n1=sec(x1,h1)
    z2,a2,r2,n2=sec(x2,h2)
    e2=q**2/2/g/a2**2+h2+z2
    e1=q**2/2/g/a1**2+h1+z1
    e3=0.5*(n1**2*q**2/r1**(4/3)/a1**2 + n2**2*q**2/r2**(4/3)/a2**2)*(x2-x1)
    return e2-e1-e3


def bisection(h1,x1,x2,q):
    ha=3.0 # lower initial value for bisection method
    hb=7.0 # upper initial value for bisection method
    for k in range(100):
        hi=0.5*(ha+hb)
        fa=func(ha,h1,x1,x2,q)
        fb=func(hb,h1,x1,x2,q)
        fi=func(hi,h1,x1,x2,q)
        if fa*fi<0: hb=hi
        if fb*fi<0: ha=hi
        #print(fa,fi,fb)
        if np.abs(hb-ha)<1e-10: break
    return hi


def main():
    g=9.8    # gravity acceleration
    q=42.0   # discharge
    # starting point (sub-critical flow)
    h1=3.9 # water depth at starting point
    x1=0.0   # origin of distance coordinate
    z1,a1,r1,n1=sec(x1,h1) # section property
    v1=q/a1  # flow velocity
    print('{0:>8s}{1:>8s}{2:>8s}{3:>8s}{4:>10s}{5:>8s}'.format('x','z','h','z+h','Energy','v'))
    print('{0:8.3f}{1:8.3f}{2:8.3f}{3:8.3f}{4:10.5f}{5:8.3f}'.format(x1,z1,h1,z1+h1,z1+h1+v1**2/2/g,v1))

    # calculation point
    #xp=np.arange(1,103,1)
    xp=np.array([11,19,47,62,69,69.5,73.5,74,77,96.625,102.125])

    # water level calculation to upstream direction
    for x2 in xp:
        h2=bisection(h1,x1,x2,q)
        z2,a2,r2,n2=sec(x2,h2)
        v2=q/a2
        print('{0:8.3f}{1:8.3f}{2:8.3f}{3:8.3f}{4:10.5f}{5:8.3f}'.format(x2,z2,h2,z2+h2,z2+h2+v2**2/2/g,v2))
        x1=x2 # distance
        h1=h2 # water depth


#==============
# Execution
#==============
if __name__ == '__main__': main()        

Calculated results are shown below. A water level drop is larger than expected.

       x       z       h     z+h    Energy       v
   0.000 562.200   3.900 566.100 566.46982   2.692
  11.000 560.500   5.971 566.471 566.47522   0.293
  19.000 560.500   5.971 566.471 566.47523   0.293
  47.000 561.060   5.410 566.470 566.47528   0.323
  62.000 561.360   5.091 566.451 566.47541   0.687
  69.000 561.500   4.950 566.450 566.47553   0.707
  69.500 562.000   4.444 566.444 566.47554   0.788
  73.500 562.000   4.444 566.444 566.47562   0.788
  74.000 561.500   4.954 566.454 566.47563   0.652
  77.000 561.500   4.954 566.454 566.47567   0.652
  96.625 563.000   3.431 566.431 566.47615   0.942
 102.125 563.000   3.431 566.431 566.47634   0.942

That's all. Thank you.

Python : Calculation of normal depth of open channel with rectangular cross section by 'scipy.optimize.brentq'

Calculation program of normal depth of open channel with rectangular cross sectionby Brent's method is introduced in this post. Basic equations for normal depth calculation of open channel with rectangular cross section are shown below.


\begin{gather}
Q=A\cdot v \\
v=\cfrac{1}{n}\cdot R^{2/3}\cdot i^{1/2} \\
R=\cfrac{b\cdot h}{b+2\cdot h}
\end{gather}
Q(known) discharge
b(known) width of rectangular open channe
n(known) Manning's roughness coefficient
i(known) gradient of channel invert
h(unknown) normal depth

Following equation is solved for h.


\begin{equation}
f=Q-\cfrac{b\cdot h}{n}\cdot \left(\cfrac{b\cdot h}{b+2\cdot h}\right)^{2/3}\cdot i^{1/2}=0
\end{equation}

The critical depth calculation of h_c is extra.

# normal depth and critical depth of rectangular cross section
import numpy as np
from scipy import optimize


def cal_hc(q,b):
    # critical depth
    g=9.8
    hc=(q**2/g/b**2)**(1/3)
    return hc


def func(h,q,b,n,i):
    f=q-b*h/n*(b*h/(b+2*h))**(2/3)*i**(1/2)    
    return f    


def main():
    q=42.0  # discharge
    b=4.0   # channel width
    n=0.014 # Manning's roughness coefficient
    i=0.001 # invert gradient
    
    h1=0.0
    h2=10.0
    hh=optimize.brentq(func,h1,h2,args=(q,b,n,i))

    print('hn=',hh) # normal depth
    
    hc=cal_hc(q,b)
    print('hc=',hc) # critical depth
    
#==============
# Execution
#==============
if __name__ == '__main__': main()

Calculation results are shown below.

hn= 3.866645305835682
hc= 2.2407023732785825

That's all. Thank you.