フーボルト法概要
フーボルト法は、航空機の振動解析のためにフーボルトが考案した方法である。本手法は、時刻$(t+Δt)$時点の位置$u(t+Δt)$を、現在位置と過去の位置の多項式近似により求める。下図のように、4点 $u(t+Δt)$、$u(t)$、$u(t-Δt)$、$u(t-2Δt)$ を通るラグランジュ補間関数により、位置$u$を3次の多項式近似として算出する。
位置$u$が3次式となるため、速度$\dot{u}$は2次式、加速度$\ddot{u}$は1次式となり、加速度を1次式として扱う線形加速度法(ニューマークβ法)やウィルソンθ法と共通した面があるが、多項式近似していることにより連続性が保たれるという点が本手法の特徴である。
ラグランジュ補間による位置uの多項式近似
上図4点を通る多項式$U(\tau)$をラグランジュ補間による公式を適用すると以下のようになる。
\begin{align}
U(\tau) &= \frac{\left(\tau-(t-Δt)\right)(\tau-t)\left(\tau-(t+Δt)\right)}
{\left((t-2Δt)-(t-Δt)\right)\left((t-2Δt)-t\right)\left((t-2Δt)-(t+Δt)\right)}u(t-2Δt) \\
&+ \frac{\left(\tau-(t-2Δt)\right)(\tau-t)\left(\tau-(t+Δt)\right)}
{\left((t-Δt)-(t-2Δt)\right)\left((t-Δt)-t\right)\left((t-Δt)-(t+Δt)\right)}u(t-Δt) \\
&+ \frac{\left(\tau-(t-2Δt)\right)(\tau-(t-Δt))\left(\tau-(t+Δt)\right)}
{\left(t-(t-2Δt)\right)\left(t-(t-Δt)\right)\left(t-(t+Δt)\right)}u(t) \\
&+ \frac{\left(\tau-(t-2Δt)\right)(\tau-(t-Δt))\left(\tau-t\right)}
{\left((t+Δt)-(t-2Δt)\right)\left((t+Δt)-(t-Δt)\right)\left((t+Δt)-t\right)}u(t+Δt) \\
\end{align}
例えば、$\tau=t+Δt$ のとき、
U(t+Δt)=u(t+Δt)
を表す。整理すると以下式となる。
\begin{align}
U(\tau) &= \frac{\left(\tau-t+Δt\right)(\tau-t)\left(\tau-t-Δt\right)}
{\left(-Δt\right)\left(-2Δt\right)\left(-3Δt\right)}u(t-2Δt) \\
&+ \frac{\left(\tau-t+2Δt\right)(\tau-t)\left(\tau-t-Δt\right)}
{Δt\left(-Δt\right)\left(-2Δt\right)}u(t-Δt) \\
&+ \frac{\left(\tau-t+2Δt\right)(\tau-t+Δt)\left(\tau-t-Δt\right)}
{2Δt・Δt\left(-Δt\right)}u(t) \\
&+ \frac{\left(\tau-t+2Δt\right)(\tau-t+Δt)\left(\tau-t\right)}
{3Δt・2Δt・Δt}u(t+Δt) \\
\end{align}
変数変換による位置u式の整理
ここで、変数変換
s=\frac{\tau-t}{Δt}
と置いて、$U(\tau)$を$U(s)$に変換する。
\begin{align}
U(s) &= \frac{\left(sΔt+Δt\right)sΔt\left(sΔt-Δt\right)}
{\left(-Δt\right)\left(-2Δt\right)\left(-3Δt\right)}u(t-2Δt) \\
&+ \frac{\left(sΔt+2Δt\right)sΔt\left(sΔt-Δt\right)}
{Δt\left(-Δt\right)\left(-2Δt\right)}u(t-Δt) \\
&+ \frac{\left(sΔt+2Δt\right)(sΔt+Δt)\left(sΔt-Δt\right)}
{2Δt・Δt\left(-Δt\right)}u(t) \\
&+ \frac{\left(sΔt+2Δt\right)(sΔt+Δt)sΔt}
{3Δt・2Δt・Δt}u(t+Δt) \\
&= \frac{(s+1)s(s-1)}{-6}u(t-2Δt)+\frac{(s+2)s(s-1)}{2}u(t-Δt) \\
&+ \frac{(s+2)(s+1)(s-1)}{-2}u(t)+\frac{(s+2)(s+1)s}{6}u(t+Δt) \\
&= \frac{s^3-s}{-6}u(t-2Δt)+\frac{s^3+s^2-2s}{2}u(t-Δt)+\frac{s^3+2s^2-s-2}{-2}u(t)+\frac{s^3+3s^2+2s}{6}u(t+Δt) \\
\end{align}
例えば、$s=1$ のとき、
U(1)=U(\tau=t+Δt)=u(t+Δt)
であることが確認できる。
位置uの微分による速度と加速度の算出
次に、$U(s)$を微分することによって、$\dot{u}(t+Δt)$、$\ddot{u}(t+Δt)$を求めていく。ここで、変数変換した$s$の$\tau$に関する微分
\frac{ds}{d\tau}=\frac{1}{Δt}
であることに注意して、$U(s)$を2階微分すると以下の式になる。
\begin{align}
\dot{U}(s) &= \frac{(3s^2-1)\frac{1}{Δt}}{-6}u(t-2Δt)+\frac{(3s^2+2s-2)\frac{1}{Δt}}{2}u(t-Δt) \\
&+\frac{(3s^2+4s-1)\frac{1}{Δt}}{-2}u(t)+\frac{(3s^2+6s+2)\frac{1}{Δt}}{6}u(t+Δt) \\
\ddot{U}(s) &= \frac{6s\frac{1}{Δt^2}}{-6}u(t-2Δt)+\frac{(6s+2)\frac{1}{Δt^2}}{2}u(t-Δt) \\
&+\frac{(6s+4)\frac{1}{Δt^2}}{-2}u(t)+\frac{(6s+6)\frac{1}{Δt^2}}{6}u(t+Δt) \\
\end{align}
よって、上式に$s=1$を代入することにより、$\dot{u}(t+Δt)$、$\ddot{u}(t+Δt)$は以下の式となる。
\begin{align}
\dot{U}(1) &= \dot{U}(\tau=t+Δt) = \dot{u}(t+Δt) \\
&= \frac{1}{Δt}\left( -\frac{1}{3}u(t-2Δt)+\frac{3}{2}u(t-Δt)-3u(t)+\frac{11}{6}u(t+Δt) \right) \\
\ddot{U}(1) &= \ddot{U}(\tau=t+Δt) = \ddot{u}(t+Δt) \\
&= \frac{1}{Δt^2}\left( -u(t-2Δt)+4u(t-Δt)-5u(t)+2u(t+Δt) \right) \\
\end{align}
満たすべき運動方程式
$(t+Δt)$時点で満たすべき振動の運動方程式は以下のようになる。質量$m$、減衰$c$、剛性$k$とする。
\begin{align}
m\ddot{u}(t+Δt)+c\dot{u}(t+Δt)+ku(t+Δt) = f(t+Δt)
\end{align}
上記までに、$\dot{u}(t+Δt)$、$\ddot{u}(t+Δt)$を位置$u$で表現できているため、上式に代入して位置$u$を代数計算によりすることができる。
フーボルト法の直接計算法
以下の流れで解を求める。
①$\dot{u}(t+Δt)、\ddot{u}(t+Δt)$を運動方程式に代入する。
②$u(t+Δt)$について解く。
具体的計算(1自由度の場合)
$\dot{u}(t+Δt)、\ddot{u}(t+Δt)$を運動方程式に代入すると以下のようになる。
\begin{align}
&m\frac{1}{Δt^2}\left(-u(t-2Δt)+4u(t-Δt)-5u(t)+2u(t+Δt) \right) \\
&+c\frac{1}{Δt}\left( -\frac{1}{3}u(t-2Δt)+\frac{3}{2}u(t-Δt)-3u(t)+\frac{11}{6}u(t+Δt) \right) \\
&+ku(t+Δt) \\
&= f(t+Δt) \\
\end{align}
少し整理すると以下のようになる。
\begin{align}
&m\left(-\frac{1}{2}u(t-2Δt)+2u(t-Δt)-\frac{5}{2}u(t)+u(t+Δt) \right) \\
&+cΔt\left( -\frac{1}{6}u(t-2Δt)+\frac{3}{4}u(t-Δt)-\frac{3}{2}u(t)+\frac{11}{12}u(t+Δt) \right) \\
&+\frac{Δt^2}{2}ku(t+Δt) \\
&= \frac{Δt^2}{2}f(t+Δt) \\
\end{align}
$u(t+Δt)$について解くと以下のようになる。
\begin{align}
\left( m+\frac{11}{12}cΔt+\frac{Δt^2}{2}k \right)u(t+Δt) &= \frac{Δt^2}{2}f(t+Δt) \\
&- \left( -\frac{5}{2}m-\frac{3}{2}cΔt \right)u(t) \\
&- \left( 2m+\frac{3}{4}cΔt \right)u(t-Δt) \\
&- \left( -\frac{1}{2}m-\frac{1}{6}cΔt \right)u(t-2Δt) \\
\end{align}
\therefore u(t+Δt) = \frac{\frac{Δt^2}{2}f(t+Δt) + \left(\frac{5}{2}m+\frac{3}{2}cΔt \right)u(t) -\left(2m+\frac{3}{4}cΔt\right)u(t-Δt) +\left(\frac{1}{2}m+\frac{1}{6}cΔt \right)u(t-2Δt)}{\left( m+\frac{11}{12}cΔt+\frac{Δt^2}{2}k \right)}
Pythonによるシミュレーション確認
Pythonによるフーボルト法の動作を計算シミュレーションで確認してみる。なお、解析解と比較できるよう簡単な1次元の振動問題を取り上げる。
(こちらのオイラー法やルンゲクッタ法で取り上げた例題と同様の問題設定である。)
問題設定
以下のような1自由度の振動モデルを考える。$m$を質点の質量、$c$を減衰、$k$をバネ定数とし、$x$が質点座標を表す。$\dot{x}$は質点速度、$\ddot{x}$は質点の加速度を表し、初期変位$x_0$と初速度$v_0$を与える。
上式の運動方程式は以下となる。
m\ddot{x}+c\dot{x}+kx = 0
解析解
解析解は以下のようになる。(導出過程はこちらオイラー法やルンゲクッタ法を参照)
\begin{align}
x &= Ae^{-\zeta\omega _nt}cos(\omega _dt - \phi) \\
A &= \sqrt{x_0^2 + \Bigl(\frac{\zeta\omega _nx_0+v_0}{\omega _d} \Bigr)^2} \\
\phi &= tan^{-1}\Bigl(\frac{\zeta\omega _nx_0+v_0}{x_0 \omega _d}\Bigr)
\end{align}
作成したPythonプログラム
import numpy as np
import matplotlib.pyplot as plt
import sympy as sp
cSimSettings = {
"SIM_START": 0.0,
"SIM_END": 5.0,
"SIM_STEP": 0.01,
}
cMdlParams = {
"M": 5.0,
"C": 16.0,
"K": 320.0,
"X0": 0.05, # 初期変位
"V0": 0.4, # 初期速度
"A0": 0.0,
"F": 0.0,
}
cNewmarkParams = {
"BETA": 1.0/6.0
}
def calcAnalyticalSolution(t):
# モデルパラメータ設定
m = cMdlParams["M"]
c = cMdlParams["C"]
k = cMdlParams["K"]
v0 = cMdlParams["V0"]
x0 = cMdlParams["X0"]
# 振動パラメータ計算
wn = np.sqrt(k/m) # 不減衰固有角振動数
zeta = c/(2.0*np.sqrt(m*k)) # 減衰比
wd = np.sqrt(1-zeta**2)*wn # 減衰固有角振動数
# 解の計算
Amp = np.sqrt((x0**2 + ((zeta*wn*x0 + v0)/wd)**2))
phi = np.arctan((zeta*wn*x0 + v0)/(x0*wd))
x = Amp*np.exp(-zeta*wn*t)*np.cos(wd*t - phi)
return x
def solverNewmarkBeta(m,c,k,xn,vn,an,fn1,dt,beta):
an1_num = fn1 - c*(vn + dt*an/2.0) - k*(xn + dt*vn + (1.0/2.0 - beta)*(dt**2)*an)
an1_den = m + dt*c/2.0 + beta*(dt**2)*k
an1 = an1_num/an1_den
vn1 = vn + dt*(an + an1)/2.0
xn1 = xn + dt*vn + (dt**2)*an/2.0 + beta*(dt**2)*(an1 - an)
return(an1, vn1, xn1)
def solverHoubolt(m,c,k,xn_2,xn_1,xn,fn1,dt):
xn1_num1 = 1.0/2.0*(dt**2)*fn1
xn1_num2 = (5.0/2.0*m + 3.0/2.0*dt*c)*xn
xn1_num3 = (-2.0*m - 3.0/4.0*dt*c)*xn_1
xn1_num4 = (1.0/2.0*m + 1.0/6.0*dt*c)*xn_2
xn1_den = m + 11.0/12.0*dt*c + 1.0/2.0*(dt**2)*k
xn1 = (xn1_num1 + xn1_num2 + xn1_num3 + xn1_num4)/xn1_den
return(xn1)
def derFunc(t, y):
# モデルパラメータ設定
m = cMdlParams["M"]
c = cMdlParams["C"]
k = cMdlParams["K"]
# 状態方程式の立式
u = [y[0], y[1]]
A = [[-c/m, -k/m], [1, 0]]
dy = np.dot(A, u)
return dy
def solverEuler(f, tn, dt, yn):
yn1 = yn + dt*f(tn, yn)
return yn1
if __name__ == "__main__":
# 変数準備(固定ステップで計算するためあらかじめ用意しておく)
sim_start = cSimSettings["SIM_START"]
sim_end = cSimSettings["SIM_END"]
sim_step = cSimSettings["SIM_STEP"]
num_step = int((sim_end-sim_start)/sim_step) + 1
t = np.linspace(sim_start, sim_end, num_step)
# 数値計算の格納結果
a_nm = np.zeros(num_step)
v_nm = np.zeros(num_step)
x_nm = np.zeros(num_step)
#a_hb = np.zeros(num_step)
v_hb = np.zeros(num_step)
x_hb = np.zeros(num_step)
# 解析解の計算
num_step = int((sim_end-sim_start)/0.01) + 1
ta = np.linspace(sim_start, sim_end, num_step)
x = calcAnalyticalSolution(ta)
# 数値計算による解
a0 = cMdlParams["A0"]
v0 = cMdlParams["V0"]
x0 = cMdlParams["X0"]
m = cMdlParams["M"]
c = cMdlParams["C"]
k = cMdlParams["K"]
f = cMdlParams["F"]
dt = sim_step
b = cNewmarkParams["BETA"]
for i in range(t.size):
if (i < 1):
# 初期値
[a_nm[i], v_nm[i], x_nm[i]] = [a0, v0, x0]
[_, v_hb[i], x_hb[i]] = [a0, v0, x0]
else:
# ニューマークβ法による計算
a_nm[i],v_nm[i], x_nm[i] = solverNewmarkBeta(m,c,k,
x_nm[i-1],
v_nm[i-1],
a_nm[i-1],
f,dt,b)
# フーボルト法による計算
if (i < 3):
# 前進オイラー法
[v_hb[i], x_hb[i]] = solverEuler(derFunc,
t[i-1],
sim_step,
[v_hb[i-1], x_hb[i-1]])
else:
x_hb[i] = solverHoubolt(m,c,k,
x_hb[i-3],
x_hb[i-2],
x_hb[i-1],
f,dt)
# グラフ確認
plt.plot(ta, x, 'k', lw=3.0, label='Analytical Solution')
plt.plot(t, x_nm, 'g', lw=2.0, label='Newmarkβ method')
plt.plot(t, x_hb, 'r--', lw=1.0, label='Houbolt method')
plt.xlim(0, 3.5)
plt.ylim(-0.04, 0.07)
plt.xlabel('Time (s)')
plt.ylabel('Displacement x(m)')
plt.legend()
plt.grid()
plt.show()
# plt.savefig("./1dofVibSim_Houbolt.png")
解析解は関数「calcAnalyticalSolution」で計算している。
参考で線形加速度法(ニューマークβ法)との比較も行った。計算は「solverNewmarkBeta」関数にて行っているが、関数の説明についてはこちらを参照。
関数 solverHoubolt の簡単な説明
def solverHoubolt(m,c,k,xn_2,xn_1,xn,fn1,dt):
xn1_num1 = 1.0/2.0*(dt**2)*fn1
xn1_num2 = (5.0/2.0*m + 3.0/2.0*dt*c)*xn
xn1_num3 = (-2.0*m - 3.0/4.0*dt*c)*xn_1
xn1_num4 = (1.0/2.0*m + 1.0/6.0*dt*c)*xn_2
xn1_den = m + 11.0/12.0*dt*c + 1.0/2.0*(dt**2)*k
xn1 = (xn1_num1 + xn1_num2 + xn1_num3 + xn1_num4)/xn1_den
return(xn1)
$(t+Δt)$時点の位置$u(t+Δt)$ を変数xn1 にて計算している。関数引数により、計算に必要な$f(t+Δt)$、$u(t)$、$u(t-Δt)$、$u(t-2Δt)$時点の位置をそれぞれ変数fn1、xn、xn_1、xn_2 により受け取って計算に使用している。今回は$(t+Δt)$時点の位置$u(t+Δt)$のみの計算を実装し、速度と加速度の計算は行っていない。
計算を進める上での注意点
フーボルト法による計算式を見てもらうとわかる通り、$u(t-Δt)$、$u(t-2Δt)$ 時点の位置の値を使用するため、計算開始時点の $u(t-Δt)$、$u(t-2Δt)$ の工夫が必要である。簡単な取り扱いとしては、初期条件の値を使用する。または、$u(t-2Δt)$ の値が求まるまではオイラー法やルンゲクッタ法の数値計算を使って計算ステップを進めることが考えられる。
今回は、$u(t-2Δt)$ の値が求まるまではオイラー法により計算ステップを進めることとした。プログラムでは以下の処理が該当箇所である。
# フーボルト法による計算
if (i < 3):
# 前進オイラー法
[v_hb[i], x_hb[i]] = solverEuler(derFunc,
t[i-1],
sim_step,
[v_hb[i-1], x_hb[i-1]])
else:
x_hb[i] = solverHoubolt(m,c,k,
x_hb[i-3],
x_hb[i-2],
x_hb[i-1],
f,dt)
プログラムの実行結果
実行結果を以下に示す。
解析解を黒実線、線形加速度法(ニューマークβ法)による数値計算結果は緑の実線、フーボルト法による数値計算結果を赤点線で示している。解析解と概ね一致、また両手法も概ね一致している。
参考文献
戸川隼人著.「有限要素法による振動解析」サイエンス社.1975年