微分方程式の数値計算 ウィルソンθ法

ウィルソンθ法概要

　ウィルソンθ法は線形加速度法の変形版である。線形加速度法は、運動方程式を$(t+Δt)$時点の解を計算しているが、ウィルソンθ法ではそれより先の$(t+θΔt)$時点(ただし、$θ>1$)に適用しているという違いがある。(線形加速度法についてはこちら参照)

\begin{align}
u(t+\theta Δt) &= u(t)+ (\thetaΔt)\dot{u}(t)+\frac{(\thetaΔt)^2}{2!}\ddot{u}(t)+\frac{(\thetaΔt)^3}{3!}\frac{\ddot{u}(t+\thetaΔt)-\ddot{u}(t)}{\thetaΔt}\\
\dot{u}(t+\thetaΔt) &= \dot{u}(t)+\thetaΔt\frac{\ddot{u}(t)+\ddot{u}(t+\thetaΔt)}{2}\\
\end{align}

$(t+θΔt)$時点で満たすべき振動の運動方程式は以下のようになる。質量$m$、減衰$c$、剛性$k$とする。

\begin{align}
m\ddot{u}(t+\thetaΔt)+c\dot{u}(t+\thetaΔt)+ku(t+\thetaΔt) = f(t+\thetaΔt)
\end{align}

線形加速度法と同じように、上式から$\ddot{u}(t+\thetaΔt)$を求めることができる。
$\ddot{u}(t+Δt)$ を求めるには、$\ddot{u}(t)$と$\ddot{u}(t+\thetaΔt)$の線形補間で$\ddot{u}(t+Δt)$を求める。

\begin{align}
\ddot{u}(t+Δt) &= \ddot{u(t)}+\frac{\ddot{u}(t+\thetaΔt) - \ddot{u}(t)}{\thetaΔt}×Δt\\
&= \left( 1-\frac{1}{\theta}\right)\ddot{u(t)}+\frac{1}{\theta}\ddot{u}(t+\thetaΔt)
\end{align}

得られた$\ddot{u}(t+Δt)$を用いて、線形加速度法同様に以下式により$u(t+Δt)、\dot{u}(t+Δt)$を求める。

\begin{align}
u(t+Δt) &= u(t)+ Δt\dot{u}(t)+\frac{(Δt)^2}{2!}\ddot{u}(t)+\frac{(Δt)^3}{3!}\frac{\ddot{u}(t+Δt)-\ddot{u}(t)}{Δt}\\
\dot{u}(t+Δt) &= \dot{u}(t)+Δt\frac{\ddot{u}(t)+\ddot{u}(t+Δt)}{2}\\
\end{align}

ウィルソンθ法の直接計算法

以下の流れで解を求める。
①$u(t+\thetaΔt)、\dot{u}(t+\thetaΔt)$を運動方程式に代入する。
②$\ddot{u}(t+\thetaΔt)$について解く。
③$\ddot{u}(t+Δt)$ を求める。
④得られた$\ddot{u}(t+Δt)$を$u(t+Δt)、\dot{u}(t+Δt)$の式に代入し解を求める。

具体的計算(1自由度の場合)

$u(t+\thetaΔt)、\dot{u}(t+\thetaΔt)$を運動方程式に代入すると以下のようになる。

\begin{align}
&m\ddot{u}(t+\thetaΔt)+c\left(\dot{u}(t)+\thetaΔt\frac{\ddot{u}(t)+\ddot{u}(t+\thetaΔt)}{2}\right) \\
&+k\left(u(t)+ (\thetaΔt)\dot{u}(t)+\frac{(\thetaΔt)^2}{2!}\ddot{u}(t)+\frac{(\thetaΔt)^3}{6}\frac{\ddot{u}(t+\thetaΔt)-\ddot{u}(t)}{\thetaΔt}\right) = f(t+\thetaΔt) \\
\end{align}

少し整理すると以下のようになる。

\begin{align}
&m\ddot{u}(t+\thetaΔt)+c\left(\dot{u}(t)+\thetaΔt\frac{\ddot{u}(t)+\ddot{u}(t+\thetaΔt)}{2}\right) \\
&+k\left(u(t)+ (\thetaΔt)\dot{u}(t)+\frac{(\thetaΔt)^2}{3}\ddot{u}(t)+\frac{(\thetaΔt)^2}{6}\ddot{u}(t+\thetaΔt)\right) = f(t+\thetaΔt)
\end{align}

$\ddot{u}(t+\thetaΔt)$について解くと以下のようになる。

\therefore \ddot{u}(t+\thetaΔt)=\frac{f(t+\thetaΔt)-c\left(\dot{u}(t)+\frac{(\thetaΔt)}{2}\ddot{u}(t)\right) - k\left(u(t)+ (\thetaΔt)\dot{u}(t)+\frac{(\thetaΔt)^2}{3}\ddot{u}(t)\right)}{m+\frac{(\thetaΔt)}{2}c+\frac{(\thetaΔt)^2}{6}k}

得られた$\ddot{u}(t+\thetaΔt)$から、$\ddot{u}(t+Δt)$を求める。

\begin{align}
\ddot{u}(t+Δt) &= \left( 1-\frac{1}{\theta}\right)\ddot{u(t)}+\frac{1}{\theta}\ddot{u}(t+\thetaΔt)
\end{align}

得られた$\ddot{u}(t+Δt)$を、以下の$u(t+Δt)、\dot{u}(t+Δt)$の式に代入し解を求める。

\begin{align}
\dot{u}(t+Δt) &= \dot{u}(t)+Δt\frac{\ddot{u}(t)+\ddot{u}(t+Δt)}{2}\\
u(t+Δt) &= u(t)+ Δt\dot{u}(t)+\frac{(Δt)^2}{3}\ddot{u}(t)+\frac{(Δt)^2}{6}\ddot{u}(t+Δt)\\
\end{align}

具体的計算(多自由度の場合)

$u(t+Δt)、\dot{u}(t+Δt)$、$\ddot{u}(t+Δt)$が多自由度の場合でも計算可能である。その場合は行列演算になることに注意して以下のようになる。(質量マトリクス$M$、減衰マトリクス$C$、剛性マトリクス$K$とする。)

\begin{align}
\vec{\ddot{u}}(t+\thetaΔt)&= \left(M+\frac{(\thetaΔt)}{2}C+\frac{(\thetaΔt)^2}{6}K\right)^{-1} \left(\vec{f}(t+\thetaΔt)-C\left(\vec{\dot{u}}(t)+\frac{(\thetaΔt)}{2}\vec{\ddot{u}}(t)\right) - K\left(\vec{u}(t)+ (\thetaΔt)\vec{\dot{u}}(t)+\frac{(\thetaΔt)^2}{3}\vec{\ddot{u}}(t)\right)\right)\\
\vec{\ddot{u}}(t+Δt) &= \left( 1-\frac{1}{\theta}\right)\vec{\ddot{u}}(t)+\frac{1}{\theta}\vec{\ddot{u}}(t+\thetaΔt)\\

\vec{\dot{u}}(t+Δt) &= \vec{\dot{u}}(t)+Δt\frac{\vec{\ddot{u}}(t)+\vec{\ddot{u}}(t+Δt)}{2}\\
\vec{u}(t+Δt) &= \vec{u}(t)+ Δt\vec{\dot{u}}(t)+\frac{(Δt)^2}{3}\vec{\ddot{u}}(t)+\frac{(Δt)^2}{6}\vec{\ddot{u}}(t+Δt)\\
\end{align}

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
}

cWilsonParams = {
    "THETA": 1.4
}


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 solverWilsonTheta(m,c,k,xn,vn,an,fnt,dt,theta):
    tdt = theta*dt
    ant_num = fnt - c*(vn + tdt*an/2.0) - k*(xn + tdt*vn + (tdt**2)*an/3.0)
    ant_den = m + tdt*c/2.0 + (tdt**2)*k/6.0
    ant = ant_num/ant_den
    coeff = 1.0/theta
    an1 = (1.0 - coeff)*an + coeff*ant
    vn1 = vn + dt*(an + an1)/2.0
    xn1 = xn + dt*vn + (dt**2)*an/3.0 + (dt**2)*an1/6.0
    return(an1, vn1, xn1)


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_wt = np.zeros(num_step)
    v_wt = np.zeros(num_step)
    x_wt = 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"]
    th = b = cWilsonParams["THETA"]
    for i in range(t.size):
        if (i < 1):
            # 初期値
            [a_nm[i], v_nm[i], x_nm[i]] = [a0, v0, x0]
            [a_wt[i], v_wt[i], x_wt[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)
            # ウィルソンΘ法による計算
            a_wt[i],v_wt[i], x_wt[i] = solverWilsonTheta(m,c,k,
                                                         x_wt[i-1],
                                                         v_wt[i-1],
                                                         a_wt[i-1],
                                                         f,dt,th)

    # グラフ確認
    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_wt, 'r--', lw=1.0, label='WilsonΘ method')
    plt.xlim(0, 3.5)
    plt.ylim(-0.04, 0.08)
    plt.xlabel('Time (s)')
    plt.ylabel('Displacement x(m)')
    plt.legend()
    plt.grid()
    plt.show()
    # plt.savefig("./1dofVibSim_WhilsonTheta.png")

解析解は関数「calcAnalyticalSolution」で計算している。
参考で線形加速度法(ニューマークβ法)との比較も行った。計算は「solverNewmarkBeta」関数にて行っているが、関数の説明についてはこちらを参照。

関数 solverWilsonTheta の簡単な説明

    tdt = theta*dt
    ant_num = fnt - c*(vn + tdt*an/2.0) - k*(xn + tdt*vn + (tdt**2)*an/3.0)
    ant_den = m + tdt*c/2.0 + (tdt**2)*k/6.0
    ant = ant_num/ant_den
    coeff = 1.0/theta
    an1 = (1.0 - coeff)*an + coeff*ant
    vn1 = vn + dt*(an + an1)/2.0
    xn1 = xn + dt*vn + (dt**2)*an/3.0 + (dt**2)*an1/6.0

「solverNewmarkBeta」関数の処理とおおよそ同じであり、$(t+θΔt)$時点の加速度$\ddot{u}(t+θΔt)$ は変数ant にて計算している。変数antとanから$(t+Δt)$時点の加速度$\ddot{u}(t+Δt)$ を変数an1にて計算している。その後の計算は「solverNewmarkBeta」関数の処理と同じである。

プログラムの実行結果

実行結果を以下に示す。

解析解を黒実線、線形加速度法(ニューマークβ法)による数値計算結果は緑の実線、ウィルソンθ法による数値計算結果を赤点線で示している。解析解と概ね一致している。ウィルソンθ法はニューマークβ法に対して先読みした時刻まで加速度が一定変化するという計算を行っているため、加速度変化量が大きい時刻0.1sや0.5sあたりでは、ウィルソンθ法の方がxの変化を大きく見積もっている傾向がみられる。

参考文献

戸川隼人著.「有限要素法による振動解析」サイエンス社.1975年