流体解析のあれこれ

流体解析をこよなく愛する有限体積法の信者によるブログ

等間隔格子における二階の微分に対する四次中心差分法の導出

有限体積法

多くの場合，二階の微分の近似には前回導出した二次中心差分法で十分であるが，さらに高精度な近似も可能である。 下図において格子 i周りの四次の項まで考慮するテイラー展開を用いて等間隔格子における二階の微分に対する二次中心差分法を導出する。

一次元の等間隔格子の概念図

f_{i+2} = f_i + \dfrac{2 \Delta}{1!} \dfrac{\partial f}{\partial x} + \dfrac{\left( 2 \Delta \right)^2}{2!} \dfrac{\partial^2 f}{\partial x^2} + \dfrac{\left( 2 \Delta \right)^3}{3!} \dfrac{\partial^3 f}{\partial x^3} + \dfrac{\left( 2 \Delta \right)^4}{4!} \dfrac{\partial^4 f}{\partial x^4}+ \mathcal{O} \left( \Delta^5 \right)

f_{i+1} = f_i + \dfrac{\Delta}{1!} \dfrac{\partial f}{\partial x} + \dfrac{\Delta^2}{2!} \dfrac{\partial^2 f}{\partial x^2} + \dfrac{\Delta^3}{3!} \dfrac{\partial^3 f}{\partial x^3} + \dfrac{\Delta^4}{4!} \dfrac{\partial^4 f}{\partial x^4} + \mathcal{O} \left( \Delta^5 \right)

f_{i-1} = f_i - \dfrac{\Delta}{1!} \dfrac{\partial f}{\partial x} + \dfrac{\Delta^2}{2!} \dfrac{\partial^2 f}{\partial x^2} - \dfrac{\Delta^3}{3!} \dfrac{\partial^3 f}{\partial x^3} + \dfrac{\Delta^4}{4!} \dfrac{\partial^4 f}{\partial x^4} + \mathcal{O} \left( \Delta^5 \right)

f_{i-2} = f_i - \dfrac{2 \Delta}{1!} \dfrac{\partial f}{\partial x} + \dfrac{\left( 2 \Delta \right)^2}{2!} \dfrac{\partial^2 f}{\partial x^2} - \dfrac{\left( 2 \Delta \right)^3}{3!} \dfrac{\partial^3 f}{\partial x^3} + \dfrac{\left( 2 \Delta \right)^4}{4!} \dfrac{\partial^4 f}{\partial x^4} + \mathcal{O} \left( \Delta^5 \right)

これらの式から一階，三階と四階の微分を消去し，整理すると次式が求められる。

\dfrac{\partial^2 \phi}{\partial x^2} = \dfrac{- \phi_{i+2} + 16 \phi_{i+1} - 30 \phi_i + 16 \phi_{i-1} - \phi_{i-2}}{12 \Delta^2}

以上より，差分差分法として二階の微分をコンピュータが計算することができる四則演算のみで近似することができた。 なお，等間隔格子を用いる場合，下図のセル界面における勾配を用いる有限体積法として二階の微分を記述することも可能である。

一次元のセル界面を含む等間隔格子の概念図

\dfrac{\partial^2 \phi}{\partial x^2} = \dfrac{- \phi_{i+2} + 16 \phi_{i+1} - 30 \phi_i + 16 \phi_{i-1} - \phi_{i-2}}{12 \Delta^2} = \dfrac{1}{\Delta} \left( \dfrac{- \phi_{i+2} + 15 \phi_{i+1} - 15 \phi_i + \phi_{i-1}}{12 \Delta} - \dfrac{- \phi_{i+1} + 15 \phi_i - 15 \phi_{i-1} + \phi_{i-2}}{12 \Delta} \right)

等間隔格子における二階の微分に対する二次中心差分法の導出

有限体積法

Navier-Stokes式中の粘性項，エネルギー保存式中の熱伝導の項など，二階の微分は多くの物理現象を表現するために用いられる。 下図において格子 i周りの二次の項まで考慮するテイラー展開を用いて等間隔格子における二階の微分に対する二次中心差分法を導出する。

一次元の等間隔格子の概念図

\phi_{i+1} = \phi_i + \dfrac{\Delta}{1!} \dfrac{\partial \phi}{\partial x} + \dfrac{\Delta^2}{2!} \dfrac{\partial^2 \phi}{\partial x^2} + \mathcal{O} \left( \Delta^3 \right)

\phi_{i-1} = \phi_i - \dfrac{\Delta}{1!} \dfrac{\partial \phi}{\partial x} + \dfrac{\Delta^2}{2!} \dfrac{\partial^2 \phi}{\partial x^2} + \mathcal{O} \left( \Delta^3 \right)

これらの式から一階の微分を消去し，整理すると次式が求められる。

\dfrac{\partial^2 \phi}{\partial x^2} = \dfrac{\phi_{i+1} - 2 \phi_i + \phi_{i-1}}{\Delta^2}

以上より，差分差分法として二階の微分をコンピュータが計算することができる四則演算のみで近似することができた。 なお，等間隔格子を用いる場合，下図のセル界面における勾配を用いる有限体積法として二階の微分を記述することも可能である。

\dfrac{\partial^2 \phi}{\partial x^2} = \dfrac{\phi_{i+1} - 2 \phi_i + \phi_{i-1}}{\Delta^2} = \dfrac{1}{\Delta} \left( \dfrac{\phi_{i+1} - \phi_i}{\Delta} - \dfrac{\phi_i - \phi_{i-1}}{\Delta} \right)

一次元のセル界面を含む等間隔格子の概念図

次回はさらに高次な離散化手法である等間隔格子における二階微分に対する四次中心差分法を導出する。

非圧縮性流体解析における圧力について

解法

いわゆる非圧縮性流体解析では，流体は圧縮しないと仮定する。別の言い方をすると，流体の密度は圧力によって変化しないことになる。流体の密度が圧力によって変化しないとすると，連続(質量保存)の式とNavier-Stokes式から密度をとおした圧力の影響はなくなり，Navier-Stokes (運動量保存)式に圧力勾配の形のみで圧力の影響が残る。そのため，非圧縮性流体解析では，圧力の絶対値は必要なく相対値のみを求めればよい。圧力の解は無数に存在することになるため，圧力を一意に定めるためには恣意的に圧力の基準点を設定する必要がある。流出境界に圧力指定境界条件を用いる場合には，流出境界面を圧力の基準点とする場合が多く，流出境界に自由流出境界条件を用いる場合には解析領域内のどこかに圧力の基準点を設定する場合が多い。非圧縮性流体解析にもかかわらず，密度を圧力の関数として評価する場合，流入境界面を圧力の基準点としない限り，流入境界面の圧力が変化することで密度が変化し，質量流量が変化するだけでなく，質量の保存性に影響を及ぼすことがあるため注意が必要である。

常微分方程式の数値解法(6)

常微分方程式

SciPyのodeintを用いて以下のOrdinary Differential Equations (ODE)を解き，文献[1]の数値解を再現する．

\displaystyle d \alpha / d \tau = s \left( \eta - \eta \alpha + \alpha -q \alpha^2 \right)
\displaystyle d \eta / d \tau = s^{-1} \left( - \eta - \eta \alpha + f \rho \right)
\displaystyle d \rho / d \tau = w \left( \alpha - \rho \right)

ここで， s = 77.27， q = 8.375 \times 10^{-6}， f = 1.0， w = 0.1610であり，初期値は \alpha = 4.0， \eta = 1.1， \rho = 4.0とした．

[1] Field, R. J. and Noyes, R. M., Oscillations in chemical systems. IV. Limit cycle behavior in a model of a real chemical reaction, J. Chem. Phys., 60(5), 1877-1884 (1974)

Field and Noyes (1974)中のFIG. 1の再現

#
# Field, R. J. and Noyes, R. M.,
# Oscillations in chemical systems. IV.
# Limit cycle behavior in a model of a real chemical reaction,
# J. Chem. Phys., 60(5), 1877-1884 (1974)
#
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# right-hand side
def rhs(y, x, s, q, f, w):
    return [s*(y[1] - y[1]*y[0] + y[0] - q*y[0]**2),
            (- y[1] - y[0]*y[1] + f*y[2])/s,
            w*(y[0] - y[2])]

# parameters
s = 77.27
q = 8.375e-6
f = 1.0
w = 0.1610

# initial condition
y0 = [4.0, 1.1, 4.0]

# output interval
x = np.arange(0.0, 330.0, 0.1)

# solve ode
y = odeint(rhs, y0, x, args=(s, q, f, w))

# output
figs = plt.figure()
fig1 = figs.add_subplot(3, 1, 1)
fig2 = figs.add_subplot(3, 1, 2)
fig3 = figs.add_subplot(3, 1, 3)
fig1.set_yscale('log')
fig1.set_ylabel(r'$\alpha$')
fig1.plot(x, y[:,0])
fig2.set_yscale('log')
fig2.set_ylabel(r'$\eta$')
fig2.plot(x, y[:,1])
fig3.set_yscale('log')
fig3.set_xlabel(r'$\tau$')
fig3.set_ylabel(r'$\rho$')
fig3.plot(x, y[:,2])
figs.tight_layout()
#plt.show()
plt.savefig('Field1974.png')

常微分方程式の数値解法(5)

常微分方程式

SciPyのodeintを用いて以下のOrdinary Differential Equations (ODE)を解き，文献[1,2]のChemical Akzo Nobel problemの数値解を再現する．

[1] Amat, S. et al., On a Variational Method for Stiff Differential Equations Arising from CHemistry Kinetics, Mathematics, 7, 459 (2019) [2] Kafash, B. et al., Chemical Akzo Nobel Problem; Mathematical Description and Its Solution Via Model Maker Software, MATCH Commum. Math. Comput. Chem., 71, 279-286 (2014)

文献[1]の再現

#
# Chemical Akzo Nobel Problem in
#
# [1] Amat, S. et al.,
#     On a Variational Method for Stiff Differential Equations
#     Arising from CHemistry Kinetics,
#     Mathematics, 7, 459 (2019)
#
# [2] Kafash, B. et al.,
#     Chemical Akzo Nobel Problem; Mathematical Description
#     and Its Solution Via Model Maker Software,
#     MATCH Commum. Math. Comput. Chem., 71, 279-286 (2014)
#
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# right-hand side
def rhs(y, x, k0, k1, k2, K, k3, klA, Ks, p, H):

    r0 = k0 * y[0]**4 * np.sqrt(y[1])
    r1 = k1 * y[2] * y[3]
    r2 = k1 * y[0] * y[4] / K
    r3 = k2 * y[0] * y[3]**2
    y5 = Ks * y[0] * y[3]
    r4 = k3 * y5**2 * np.sqrt(y[1])
    Fin = klA * (p / H - y[1])

    yp0 = - 2.0 * r0 + r1 - r2 - r3
    yp1 = - 0.5 * r0 - r3 - 0.5e0 * r4 + Fin
    yp2 = r0 - r1 + r2
    yp3 = - r1 + r2 - 2.0 * r3
    yp4 = r1 - r2 + r4
    return [yp0, yp1, yp2, yp3, yp4]

# parameters
k0 = 18.7
k1 = 0.58
k2 = 0.09
K = 34.4
k3 = 0.42
klA = 3.3
Ks = 115.83
p = 0.9
H = 737.0

# initial condition
# 0: FLB, 1: CO2, 2: FLBT, 3: ZHU, 4: ZLA, 5: FLB.ZHU
FLB     = 0.444
CO2     = 0.00123
FLBT    = 0.0
ZHU     = 0.007
ZLA     = 0.0
FLB_ZHU = 0.367
ZHU     = FLB_ZHU/(Ks*FLB)
y0      = [FLB, CO2, FLBT, ZHU, ZLA]

# output interval
x = np.arange(0.0, 200.0, 0.1)

# solve ode
y = odeint(rhs, y0, x, args = (k0, k1, k2, K, k3, klA, Ks, p, H))

# output
plt.xlabel('t')
plt.ylabel('y')
plt.plot(x, y[:,0],           label='FLB')
plt.plot(x, y[:,1]*1e2,       label='CO2 (\u00d7$10^2$)')
plt.plot(x, y[:,2],           label='FLBT')
plt.plot(x, y[:,3]*1e1,       label='ZHU (\u00d7$10^1$)')
plt.plot(x, y[:,4]*1e1,       label='ZLA (\u00d7$10^1$)')
plt.plot(x, Ks*y[:,0]*y[:,3], label='FLB.ZHU')
plt.legend()
#plt.show()
plt.savefig('ChemicalAkzoNobelProblem.png')

常微分方程式の数値解法(4)

常微分方程式

SciPyのodeintを用いて以下のOrdinary Differential Equations (ODE)を解き，文献[1]のProblem HIRESの数値解を再現する．

\displaystyle y_1^{\prime} = - 1.7 y_1 + 0.43 y_2 + 8.32 y_3 + 0.0007
\displaystyle y_2^{\prime} = 1.7 y_1 - 8.75 y_2
\displaystyle y_3^{\prime} = - 10.03 y_3 + 0.43 y_4 + 0.035 y_5
\displaystyle y_4^{\prime} = 8.32 y_2 + 1.71 y_3 - 1.12 y_4
\displaystyle y_5^{\prime} = - 1.745 y_5 + 0.43 y_6 + 0.43 y_7
\displaystyle y_6^{\prime} = - 280 y_6 y_8 + 0.69 y_4 + 1.71 y_5 - 0.43 y_6 + 0.69 y_7
\displaystyle y_7^{\prime} = 280 y_6 y_8 - 1.81 y_7
\displaystyle y_8^{\prime} = - 280 y_6 y_8 + 1.81 y_7

[1] Amat, S. et al., On a Variational Method for Stiff Differential Equations Arising from Chemistry Kinetics, Mathematics, 7, 459 (2019)

Amat, S. (2019)のFigure 4の再現

#
# Problem HIRES
#
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# right-hand side
def rhs(y, x):
        return[- 1.71*y[0] + 0.43*y[1] + 8.32*y[2] + 0.0007,
                 1.71*y[0] - 8.75*y[1],
               - 10.03*y[2] + 0.43*y[3] + 0.035*y[4],
                 8.32*y[1] + 1.71*y[2] - 1.12*y[3],
               - 1.745*y[4] + 0.43*y[5] + 0.43*y[6],
               - 280.0*y[5]*y[7] + 0.69*y[3] + 1.71*y[4] - 0.43*y[5] + 0.69*y[6],
                 280.0*y[5]*y[7] - 1.81*y[6],
               - 280.0*y[5]*y[7] + 1.81*y[6]]

# initial condition
y0 = [1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0057]

# output interval
x = np.arange(0.0,400.0,0.01)

# solve ode
y = odeint(rhs, y0, x)

# outout
figs = plt.figure(figsize=(6.4,4.8*2))
fig1 = figs.add_subplot(4, 2, 1)
fig2 = figs.add_subplot(4, 2, 2)
fig3 = figs.add_subplot(4, 2, 3)
fig4 = figs.add_subplot(4, 2, 4)
fig5 = figs.add_subplot(4, 2, 5)
fig6 = figs.add_subplot(4, 2, 6)
fig7 = figs.add_subplot(4, 2, 7)
fig8 = figs.add_subplot(4, 2, 8)
fig1.set_xlim(0.0,5.0)
fig2.set_xlim(0.0,5.0)
fig3.set_xlim(0.0,5.0)
fig4.set_xlim(0.0,5.0)
fig5.set_xlim(0.0,400.0)
fig6.set_xlim(0.0,400.0)
fig7.set_xlim(0.0,5.0)
fig8.set_xlim(0.0,5.0)
fig1.set_ylabel('y1')
fig2.set_ylabel('y2')
fig3.set_ylabel('y3')
fig4.set_ylabel('y5')
fig5.set_ylabel('y5')
fig6.set_ylabel('y6')
fig7.set_ylabel('y7')
fig8.set_ylabel('y8')
fig7.set_xlabel('t')
fig8.set_xlabel('t')
fig1.plot(x, y[:,0])
fig2.plot(x, y[:,1])
fig3.plot(x, y[:,2])
fig4.plot(x, y[:,3])
fig5.plot(x, y[:,4])
fig6.plot(x, y[:,5])
fig7.plot(x, y[:,6])
fig8.plot(x, y[:,7])
figs.tight_layout()
#plt.show()
figs.savefig('tmp.png')

常微分方程式の数値解法(3)

常微分方程式

SciPyのodeintを用いて以下の厳密解のあるOrdinary Differential Equations (ODE)を解き，文献[1]の数値解を再現する．

\displaystyle u_1^{\prime} = 9 u_1 + 24 u_2 + 5 \cos t - \dfrac{1}{3} \sin t
\displaystyle u_2^{\prime} = - 24 u_1 - 51 u_2 - 9 \cos t + \dfrac{1}{3} \sin t

なお，文献[1]の u_2^{\prime}の右辺第3項の -95 \cos tは -9 \cos tの間違いである．これら常微分方程式の厳密解は以下に示すとおりである．

\displaystyle u_1 \left( t \right) = 2 e^{-3 t} - e^{-39 t} + \dfrac{1}{3} \cos t
\displaystyle u_2 \left( t \right) = - e^{-3 t} + 2 e^{-39 t} - \dfrac{1}{3} \cos t
\displaystyle

[1] Thohura, S. and Rahman, A., Numerical Approach for Solving Stiff Differential Equations: A Comparative Study, GLOBAL JOURNAl OF SCIENCE FRONTIER RESEARCH MATHEMATICS AND DECISION SCIENCES, 13 (6) (2013)

文献[1]のFig. 2の再現

#
# Thohura, S. and Rahman, A.,
# Numerical Approach for Solving Stiff Differential Equations:
# A Comparative Study,
# GLOBAL JOURNAl OF SCIENCE FRONTIER RESEARCH
# MATHEMATICS AND DECISION SCIENCES, 13 (6) (2013)
#
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# right-hand side
def rhs(y, x):
    return [ 9.0*y[0] + 24.0*y[1] + 5.0*np.cos(x) - np.sin(x)/3.0,
           -24.0*y[0] - 51.0*y[1] - 9.0*np.cos(x) + np.sin(x)/3.0]

# exact solution
def exact(x):
    return [2.0*np.exp(-3.0*x) -     np.exp(-39.0*x) + np.cos(x)/3.0,
           -    np.exp(-3.0*x) + 2.0*np.exp(-39.0*x) - np.cos(x)/3.0]

# initial condition
y0 = [4.0/3.0, 2.0/3.0]

# output interval
x = np.arange(0.0, 1.0, 0.01)

# solve ode
y = odeint(rhs, y0, x)

#y_e = [[0 for i in range(100)] for j in range(2)]
y_e = np.empty((100, 2))

# exact solution
i = 0
for t in x:
    y_e[i,:] = exact(t)
    i += 1

# output
plt.xlabel('t')
plt.ylabel('u')
plt.plot(x, y_e[:,0],marker='o', markerfacecolor='none',color='blue')
plt.plot(x, y_e[:,1],marker='o', markerfacecolor='none',color='orange')
plt.plot(x, y[:,0], label='u1', color='blue')
plt.plot(x, y[:,1], label='u2', color='orange')
plt.legend()
#plt.show()
plt.savefig('Thohura2013-Fig2.png')