第4-3回　有限体積法の基礎(Burgers方程式)[python]

$$\frac{\partial u}{\partial t}+ u\frac{\partial u}{\partial x} =0$$

今までのまとめとして、上記のバーガース方程式(Burgers eq.)を解いていきます。pythonのコードは最後に張っておきます。

この方程式は一次元Navier-Stokes方程式の圧力項を無視したものになります。

有限体積法のため、上記の式を保存型に変形します。

$$\frac{\partial u}{\partial t}+ \frac{\partial}{\partial x}\left(\frac{u^2}{2}\right) =0$$

保存量 \(Q=u\),フラックス \(F=\frac{u^2}{2}\)です。

有限体積法による離散化は第4-1回で示した通りですので、次のようになります。

$$ \frac{Q_{i}^{n+1}-Q_{i}^n}{\Delta t}+\frac{ ( \tilde{f} _{i+1/2}-
\tilde{f} _{i-1/2}) }{\Delta x}$$

さらに\(Q=u\)と第4-2回で説明したフラックスを当てはめると次の式が求まります。

$$u^{n+1}_i=u^n_i- \frac{ \Delta t}{ \Delta x } ( \tilde{f} _{i+1/2}-
\tilde{f} _{i-1/2}) $$

$$=\frac{1}{4} [(u_{i+1}^2-u_{i-1}^2)-| u_{i+1}+u_{i}|(u_{i+1}-u_{i})+| u_{i}+u_{i-1}|(u_{i}-u_{i-1})] $$

この式を用いて、次のように初期値を決め計算を行います。

時間ステップ数 : \(nstep=500\)
時間刻み幅 : \(dt=0.005\)
x方向に取るセル数 : \(nx=100\)
x方向の刻み幅 : \(dx=0.01\)
初期条件 : \(u(x)=1.0 (x \leq 1.0), u(x)=0.0 (x>0.1)\)
境界条件 : \(u(0)=1, u(L)=0\)

解析結果になります。

移流していく様子が捉えられました。一次精度ですので、垂直の状態で移動はしてませんが、今後の足掛かりには十分な結果が得られました。

コードは今後のことも考え、編集しやすいように関数を定義しました。
今後はクラスなども使うかもしれませんが、よろしくお願いします。

by hide

import math
import numpy
import matplotlib.pyplot as pyplot
import matplotlib.animation as animation

# -----------------------
# -- initial condition --
# -----------------------
nstep = 500               # 時間ステップ数
dt = 0.005                # 時間刻み幅
nx = 100                  # 空間分割数
dx = 0.01                 # 空間刻み幅

# -----------------------
# -- function          --
# -----------------------

def setup():  # 初期条件を計算
    global u, x           # グローバル変数の定義

    x = [0.0]*nx          # 位置xの代入
    for i in range(nx):
        x[i] = (dx*(i+0.5))

    u = [0.0]*nx          # 速度uの代入
    for i in range(nx):
        if x[i] <= 0.1:
            u[i] = 1.0

        else:
            u[i] = 0.0

def cal_u():  # 時間ステップを進める計算
    global u
    new_u = [0.0]*nx     # 新しいuを用意

    for j in range(1, nx-1):
        f_j_p12 = (u[j+1]**2+u[j]**2)/4-abs(u[j+1]+u[j]) * \
            (u[j+1]-u[j])/4       
        f_j_m12 = (u[j]**2+u[j-1]**2)/4-abs(u[j]+u[j-1]) * \
            (u[j]-u[j-1])/4         

        new_u[j] = (u[j]-dt/dx*(f_j_p12-f_j_m12)) # メインの計算
    u = new_u[:]         # uの上書き

def bound():  # 境界条件の代入
    global u
    u[0] = 1.0
    u[nx-1] = 0.0

# -----------------------
# -- preparation       --
# -----------------------
setup()                  # 初期値の計算

fig = pyplot.figure()    # 図の準備
img = pyplot.plot(x, u)
ims = []
ims.append(img)

# -----------------------
# -- main              --
# -----------------------
for i in range(nstep):
    cal_u()              # uの計算
    bound()              # 境界の計算

    img = pyplot.plot(x, u)  # 図の作成
    ims.append(img)

pyplot.xlabel("length L [m]")
pyplot.ylabel("velocity [m/s]")

ani = animation.ArtistAnimation(fig, ims, interval=10)
pyplot.show()