Python: Brent法で露出管の限界座屈圧力に対する必要板厚を求める

続けざまだが、必要に迫られて作っている、Brent法による数値計算事例をもう一つアップしておく。「Python便利！」である。

露出型水圧鉄管の限界座屈圧力は以下の式で計算される。

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

$p_k$	鉄管の限界座屈圧力
$t$	板厚（余裕厚を除く）
$D_0'$	設計外径
$E_s$	鉄管の弾性係数 (=206,000MPa)
$\nu_s$	鉄管のポアソン比 (=0.3)

ここで、通常余裕厚は、$\epsilon=1.5 mm$ とされ、余裕厚を含む設計板厚 $t_0$ および設計外圧 $D_0'$ は以下の通り計算される。

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

ここに、$D_0$ は設計内径である。

露出型水圧鉄管の場合、抜水時などへの万が一の対応として、0.02MPaの外圧（実質は負圧）に耐えられるよう設計することが水門鉄管技術基準で定められている。これに対する必要板厚は、安全率1.5を見込んで、限界座屈圧力 $p_k=0.03 MPa$ 以上となるよう板厚を定める必要がある。

ここでは、与えられた限界座屈圧力を満足する板厚をBrent法によって計算することを考える。
限界座屈圧力算定式を変形し、以下の $f=0$ となる非線形方程式を $t$ について解けば良いことになる。

\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}

Brent法での解法に必要な２つの初期値は、プログラム中で、t1=1.0, t2=50.0 として与えている。

プログラム事例を以下に示す。

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.02   # 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() 

計算結果は以下の通り。

t + eps= 15.695548237490522
t0= 16.0

That's all.
Thank you.