1次元移流方程式をSUPG(Streamline Upwind/Petrov-Galerkin)法に基づく安定化有限要素法によって解いてみます．
言語はPythonを用います．使用したバージョンは Python 3.6.1です．
参考文献：
計算力学(第2版)-有限要素法の基礎（森北出版）
第3版 有限要素法による流れのシミュレーション OpenMPに基づくFortranソースコード付（丸善出版）
1次元移流方程式
1次元移流方程式は以下のような1階の偏微分方程式で表されます．
\cfrac{\partial f}{\partial t} + u\ \cfrac{\partial f}{\partial x} = 0
ここで，$f$ は移流される物理量（例えば濃度），$u$ は移流速度です．
安定化有限要素法による離散化
解析領域 $\Omega$ を$M$ 個の1次要素の要素領域 $\Omega_e$に分割すると以下の離散化方程式が得られます．
\int_\Omega \omega\ \cfrac{\partial f}{\partial t}\ d\Omega\ + \int_\Omega \omega\ u\ \cfrac{\partial f}{\partial x}\ d\Omega\ + \sum_{e = 1}^M\int_{\Omega_e} \tau_e\ u\ \cfrac{\partial \omega}{\partial x}\left(\cfrac{\partial f}{\partial t}\ + u\ \cfrac{\partial f}{\partial x}\right)\ d\Omega = 0
ここで，$\omega$ は重み関数，$\tau_e$ は安定化パラメータであり，一般的に非定常問題の場合は要素長 $h_e$ と時間刻み$\Delta t$ を用いて
\tau_e = \left(\left( \cfrac{2}{\Delta t}\right)^2 + \left(\cfrac{2|u|}{h_e}\right)^2\ \right)^{-\frac{1}{2}}
と表されます．ある要素 $e$ に着目すると
\int_{\Omega_e} \omega\ \cfrac{\partial f}{\partial t}\ d\Omega\ + \int_{\Omega_e} \omega\ u\ \cfrac{\partial f}{\partial x}\ d\Omega\ + \int_{\Omega_e} \tau_e\ u\ \cfrac{\partial \omega}{\partial x}\left(\cfrac{\partial f}{\partial t}\ + u\ \cfrac{\partial f}{\partial x}\right)\ d\Omega = 0
が得られます．2節点1次要素を用いて要素内の $\omega$, $f$ 及び $u$ をそれぞれ
\begin{align}
\omega(x) &= \cfrac{x_2^e - x}{h_e}\ \omega_1^e + \cfrac{x - x_1^e}{h_e}\ \omega_2^e\\
&=
\begin{pmatrix}
N_1^e & N_2^e
\end{pmatrix}
\begin{pmatrix}
\omega_1^e\\
\omega_2^e\\
\end{pmatrix}
= {\bf N}_e{\bf \omega}_e,
\end{align}
\begin{align}
f(x) &= \cfrac{x_2^e - x}{h_e}\ f_1^e + \cfrac{x - x_1^e}{h_e}\ f_2^e\\
&=
\begin{pmatrix}
N_1^e & N_2^e
\end{pmatrix}
\begin{pmatrix}
f_1^e\\
f_2^e\\
\end{pmatrix}
= {\bf N}_e{\bf f}_e,
\end{align}
\begin{align}
u(x) &= \cfrac{x_2^e - x}{h_e}\ u_1^e + \cfrac{x - x_1^e}{h_e}\ u_2^e\\
&=
{\begin{pmatrix}
N_1^e & N_2^e
\end{pmatrix}}
\begin{pmatrix}
u_1^e\\
u_2^e\\
\end{pmatrix}
= {\bf N}_e{\bf u}_e
\end{align}
と補間します．同様に $f$ の時間微分も
\begin{align}
\dot{f}(x) &= \cfrac{x_2^e - x}{h_e}\ \dot{f}_1^e + \cfrac{x - x_1^e}{h_e}\ \dot{f}_2^e\\
&=
{\begin{pmatrix}
N_1^e & N_2^e
\end{pmatrix}}
\begin{pmatrix}
\dot{f}_1^e\\
\dot{f}_2^e\\
\end{pmatrix}
= {\bf N}_e\dot{{\bf f}}_e
\end{align}
と補間します．これら代入すると以下の式が得られます．
\int_{\Omega_e} {\bf N}_e{\bf \omega}_e {\bf N}_e\dot{\bf f}_e\ d\Omega\
+ \int_{\Omega_e} {\bf N}_e{\bf \omega}_e\ {\bf N}_e{\bf u}_e\cfrac{\partial {\bf N}_e}{\partial x}\ {\bf f}_e\ d\Omega\
+ \int_{\Omega_e} \tau_e\ {\bf N}_e{\bf u}_e\cfrac{\partial {\bf N}_e}{\partial x}{\bf \omega}_e\left({\bf N}_e\dot{{\bf f}}_e\ + {\bf N}_e{\bf u}_e\cfrac{\partial {\bf N}_e}{\partial x}{\bf f}_e\right)\ d\Omega = 0
形状関数 $N$ の空間微分は
\begin{align}
\cfrac{\partial {\bf N}_e}{\partial x}
&=
{\begin{pmatrix}
\cfrac{\partial N_1^e}{\partial x} & \cfrac{\partial N_2^e}{\partial x}
\end{pmatrix}}
=
{\begin{pmatrix}
-\cfrac{1}{h_e} & \cfrac{1}{h_e}
\end{pmatrix}}
\end{align}
となり，要素内では一定値となります．
${\bf N}_ e\omega_e = ({\bf N}_ e \omega_e)^T = \omega_e^T{\bf N}_ e^T$ 等を用いて式変形していくと
\omega_e^T\left(\int_{\Omega_e} {\bf N}_e^T{\bf N}_e\ d\Omega\ \dot{\bf f}_e\
+ \int_{\Omega_e} {\bf N}_e^T\ {\bf N}_e{\bf u}_e\ \cfrac{\partial {\bf N}_e}{\partial x}\ d\Omega\ {\bf f}_e\
+ \int_{\Omega_e} \tau_e\ \left(\cfrac{\partial {\bf N}_e}{\partial x}\right)^T{\bf u}_e^T{\bf N}_e^T{\bf N}_e\ d\Omega\ \dot{{\bf f}}_e\
+ \int_{\Omega_e} \tau_e\ \left(\cfrac{\partial {\bf N}_e}{\partial x}\right)^T{\bf u}_e^T{\bf N}_e^T{\bf N}_e\ d\Omega\ {\bf u}_e\cfrac{\partial {\bf N}_e}{\partial x}\ {{\bf f}}_e\right)= 0
となります．重み関数は任意に選んだ定数なので，有限要素方程式は以下のようになります．
\left({\bf M}_e + {\bf M}_{\delta e}\right)\ \dot{{\bf f}}_e + \left({\bf S}_e + {\bf S}_{\delta e}\right)\ {\bf f}_e = 0
${\bf M}_ e$ は質量行列，${\bf S}_ e$ は移流行列，${\bf M}_ {\delta e}$, ${\bf S}_ {\delta e}$ はそれぞれSUPG項から生じる質量行列と移流行列です．
これらの行列は線座標を用いて求めることが可能であり，それぞれ
{\bf M}_e = \int_{\Omega_e} {\bf N}_e^T{\bf N}_e\ d\Omega
= \cfrac{h_e}{6}
\begin{pmatrix}
2 & 1\\
1 & 2
\end{pmatrix},
\begin{align}
{\bf A}_e &= \int_{\Omega_e} {\bf N}_e^T\ {\bf N}_e{\bf u}_e\cfrac{\partial {\bf N}_e}{\partial x}\ d\Omega\\
&= \int_{\Omega_e} {\bf N}_e^T\ {\bf N}_e\ d\Omega\ {\bf u}_e\cfrac{\partial {\bf N}_e}{\partial x}\\
&= \cfrac{h_e}{6}
\begin{pmatrix}
2 & 1\\
1 & 2
\end{pmatrix}
\begin{pmatrix}
u_1^e\\
u_2^e
\end{pmatrix}
\begin{pmatrix}
-\cfrac{1}{h_e} & \cfrac{1}{h_e}
\end{pmatrix}\\
&= \cfrac{1}{6}
\begin{pmatrix}
2 & 1\\
1 & 2
\end{pmatrix}
\begin{pmatrix}
-u_1^e & u_1^e\\
-u_2^e & u_2^e
\end{pmatrix}\\
&= \cfrac{1}{6}
\begin{pmatrix}
-2u_1^e - u_2^e & 2u_1^e + u_2^e\\
-u_1^e - 2u_2^e & u_1^e + 2u_2^e
\end{pmatrix},
\end{align}
\begin{align}
{\bf M}_{\delta e} &= \int_{\Omega_e} \tau_e\ \left(\cfrac{\partial {\bf N}_e}{\partial x}\right)^T{\bf u}_e^T{\bf N}_e^T{\bf N}_e\ d\Omega\\
&= \tau_e\ \left(\cfrac{\partial {\bf N}_e}{\partial x}\right)^T{\bf u}_e^T\ \int_{\Omega_e} {\bf N}_e^T{\bf N}_e\ d\Omega\\
&= \tau_e
\begin{pmatrix}
-\cfrac{1}{h_e}\\
\cfrac{1}{h_e}
\end{pmatrix}
\begin{pmatrix}
u_1^e & u_2^e
\end{pmatrix}
\cfrac{h_e}{6}
\begin{pmatrix}
2 & 1\\
1 & 2
\end{pmatrix}\\
&= \cfrac{\tau_e}{6}
\begin{pmatrix}
-2u_1^e - u_2^e & -u_1^e - 2u_2^e\\
2u_1^e + u_2^e & u_1^e + 2u_2^e
\end{pmatrix},
\end{align}
\begin{align}
{\bf A}_{\delta e} &= \int_{\Omega_e} \tau_e\ \left(\cfrac{\partial {\bf N}_e}{\partial x}\right)^T{\bf u}_e^T{\bf N}_e^T{\bf N}_e\ d\Omega\ {\bf u}_e\cfrac{\partial {\bf N}_e}{\partial x}\\
&= \tau_e\ \left(\cfrac{\partial {\bf N}_e}{\partial x}\right)^T{\bf u}_e^T \int_{\Omega_e} {\bf N}_e^T{\bf N}_e\ d\Omega\ {\bf u}_e\cfrac{\partial {\bf N}_e}{\partial x}\\
&=\tau_e
\begin{pmatrix}
-\cfrac{1}{h_e}\\
\cfrac{1}{h_e}
\end{pmatrix}
\begin{pmatrix}
u_1^e & u_2^e
\end{pmatrix}
\cfrac{h_e}{6}
\begin{pmatrix}
2 & 1\\
1 & 2
\end{pmatrix}
\begin{pmatrix}
u_1^e\\
u_2^e
\end{pmatrix}
\begin{pmatrix}
-\cfrac{1}{h_e} & \cfrac{1}{h_e}
\end{pmatrix}\\
&= \cfrac{\tau_e}{3h_e}\left({u_1^e}^2 + u_1^e u_2^e + {u_2^e}^2 \right)
\begin{pmatrix}
1 & -1\\
-1 & 1
\end{pmatrix}
\end{align}
となります．これらを使って全体の有限要素方程式
\left({\bf M} + {\bf M}_{\delta}\right)\ \dot{{\bf f}} + \left({\bf S} + {\bf S}_{\delta}\right)\ {\bf f} = 0
を組み立てていきます．時間方向の離散化にクランク・ニコルソン法を用いると
\left({\bf M} + {\bf M}_{\delta}\right)\ \cfrac{{\bf f}^{n+1} - {\bf f}^n}{\Delta t} + \left({\bf S} + {\bf S}_{\delta}\right)\ \cfrac{1}{2}\left({\bf f}^{n+1} + {\bf f}^n\right) = 0
が得られます．$n$ は時間ステップです．さらに変形すると
\left(\cfrac{1}{\Delta t}\left({\bf M} + {\bf M}_{\delta}\right) + \cfrac{1}{2} \left({\bf S} + {\bf S}_{\delta}\right)\right) {\bf f}^{n+1}
=
\left(\cfrac{1}{\Delta t}\left({\bf M} + {\bf M}_{\delta}\right) - \cfrac{1}{2} \left({\bf S} + {\bf S}_{\delta}\right)\right) {\bf f}^n
となります．この式をよく見てみると，左辺が未知量，右辺は既知量で，${\bf A}{\bf x} = {\bf b}$ という連立1次方程式の形をしていることが分かります．つまり，この連立1次方程式を解けば次の時刻の物理量が得られるという事です．
Pythonによる実装
各要素ごとに質量・移流行列を求めて全体行列に足しこんでいきます．コードを以下に示します．なお，上記の式を愚直にコーディングしただけであり，連立1次方程式を解くためにNumPyのlinalg.solveを用いています．本来は高速化のために三重対角行列アルゴリズムや反復法などを用いるべきですが，今回は簡単のためlinalg.solveを用いています．
解析領域は$0 \le x \le 2$ とし，境界条件は $f(x = 0, t) = f(x = 2, t) = 0$ というディリクレ境界条件を与えます．
初期条件は$0 \le x \le 1$ で $f(x, t = 0) = \sin(\pi x)$，それ以外で $f(x, t = 0) = 0$ とします．移流速度は全領域で一様で $u = 1.0$ を与えます．
実装やアニメーションの作成には以下のサイトを参考にしました．
【NumPy】連立方程式を解く numpy.linalg.solve
Pythonで波動方程式の数値計算 & 結果のアニメーション
1d_advection.py
import math
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
# 各種定数の設定
dt = 0.005 # 時間刻み
xmin = 0.0
xmax = 2.0
N = 200
dx = (xmax - xmin)/N # 空間刻み
U = 1.0 # 移流速度
M = 200
# 変数の初期化
f = np.array([[0.0]*(N + 1) for i in range(M + 1)])
u = np.zeros([N + 1])
# マトリクスの初期化
A = np.array([[0.0]*(N + 1) for i in range(N + 1)])
b = np.array([0.0]*(N + 1))
# 初期条件
for i in range(N + 1):
u[i] = U
if i*dx <= 1.0:
f[0, i] = math.sin(math.pi*i*dx)
else:
f[0, i] = 0.0
for ti in range(0, M):
A.fill(0.0)
b.fill(0.0)
for i in range(N):
u1 = u[i]
u2 = u[i + 1]
tau = 1.0/math.sqrt((2.0/dt)**2 + (2.0*U/dx)**2)
# 質量行列
A[i ,i ] = A[i ,i ] + 1.0/dt*2.0/6.0*dx
A[i ,i + 1] = A[i ,i + 1] + 1.0/dt*1.0/6.0*dx
A[i + 1,i ] = A[i + 1,i ] + 1.0/dt*1.0/6.0*dx
A[i + 1,i + 1] = A[i + 1,i + 1] + 1.0/dt*2.0/6.0*dx
# 質量行列(SUPG)
A[i ,i ] = A[i ,i ] + tau/dt/6.0*(-2.0*u1 - u2)
A[i ,i + 1] = A[i ,i + 1] + tau/dt/6.0*(- u1 - 2.0*u2)
A[i + 1,i ] = A[i + 1,i ] + tau/dt/6.0*( 2.0*u1 + u2)
A[i + 1,i + 1] = A[i + 1,i + 1] + tau/dt/6.0*( u1 + 2.0*u2)
# 移流行列
A[i ,i ] = A[i ,i ] + 0.5/6.0*(-2.0*u1 - u2)
A[i ,i + 1] = A[i ,i + 1] + 0.5/6.0*( 2.0*u1 + u2)
A[i + 1,i ] = A[i + 1,i ] + 0.5/6.0*(- u1 - 2.0*u2)
A[i + 1,i + 1] = A[i + 1,i + 1] + 0.5/6.0*( u1 + 2.0*u2)
# 移流行列(SUPG)
A[i ,i ] = A[i ,i ] + 0.5*tau/3.0/dx*(u1**2 + u1*u2 + u2**2)
A[i ,i + 1] = A[i ,i + 1] - 0.5*tau/3.0/dx*(u1**2 + u1*u2 + u2**2)
A[i + 1,i ] = A[i + 1,i ] - 0.5*tau/3.0/dx*(u1**2 + u1*u2 + u2**2)
A[i + 1,i + 1] = A[i + 1,i + 1] + 0.5*tau/3.0/dx*(u1**2 + u1*u2 + u2**2)
for i in range(N):
f1 = f[ti, i]
f2 = f[ti, i + 1]
u1 = u[i]
u2 = u[i + 1]
tau = 1.0/math.sqrt((2.0/dt)**2 + (2.0*U/dx)**2)
# 質量行列
b[i ] = b[i ] + 1.0/dt*dx/6.0*(2.0*f1 + f2)
b[i + 1] = b[i + 1] + 1.0/dt*dx/6.0*( f1 + 2.0*f2)
# 質量行列(SUPG)
b[i ] = b[i ] + tau/dt/6.0*((-2.0*u1 - u2)*f1 + (-u1 - 2.0*u2)*f2)
b[i + 1] = b[i + 1] + tau/dt/6.0*(( 2.0*u1 + u2)*f1 + ( u1 + 2.0*u2)*f2)
# 移流行列
b[i ] = b[i ] - 0.5/6.0*((-2.0*u1 - u2)*f1 + (2.0*u1 + u2)*f2)
b[i + 1] = b[i + 1] - 0.5/6.0*((- u1 - 2.0*u2)*f1 + ( u1 + 2.0*u2)*f2)
# 移流行列(SUPG)
b[i ] = b[i ] - 0.5*tau/dx/3.0*(u1**2 + u1*u2 + u2**2)*( f1 - f2)
b[i + 1] = b[i + 1] - 0.5*tau/dx/3.0*(u1**2 + u1*u2 + u2**2)*(-f1 + f2)
# 境界条件
for i in range(N + 1):
A[0, i] = 0.0
A[N, i] = 0.0
A[0, 0] = 1.0
A[N, N] = 1.0
b[0] = 0.0
b[N] = 0.0
f[ti + 1,] = np.linalg.solve(A, b)
x = np.linspace(xmin, xmax, N + 1) # x軸の設定
fig = plt.figure(figsize=(6,4))
ax = fig.add_subplot(1,1,1)
# アニメ更新用の関数
def update_func(i):
# 前のフレームで描画されたグラフを消去
ax.clear()
ax.plot(x, f[i, :], color='blue')
ax.scatter(x, f[i, :], color='blue')
# 軸の設定
ax.set_ylim(-0.1, 1.1)
# 軸ラベルの設定
ax.set_xlabel('x', fontsize=12)
ax.set_ylabel('f', fontsize=12)
# サブプロットタイトルの設定
ax.set_title('Time: ' + '{:.2f}'.format(dt*i))
ani = animation.FuncAnimation(fig, update_func, frames=M, interval=100, repeat=True)
# アニメーションの保存
#ani.save('1d_advection.gif', writer='imagemagick')
# 表示
plt.show()
計算結果
若干のアンダーシュートが見られますが，サイン波が移流しているのが確認できました．