第5-3回　オイラー方程式（Sod shock tube）[python]

$$\frac{\partial Q}{\partial t}+ \frac{\partial F}{\partial x} =0$$

$$Q=\left\{ \begin{array}{c} \rho \\ \rho u \\ e \end{array}\right\} , F=\left\{ \begin{array}{c} \rho u \\ p+\rho u^2 \\ (e+p)u \end{array}\right\} $$

今回は今までのまとめとして、Sod shock tubeと呼ばれる衝撃波管に対して数値解析を行います。この衝撃波管は厳密解が求められており、数値解の検証によく使われます。

pythonのコードは最後に書いておきます。

衝撃波管の条件は以下の通りです。

$$\left(\begin{array}{c} u_L\\rho_L\\ p_L \end{array}\right) =\left(\begin{array}{c} 0.0\\1.0\\1.0 \end{array}\right) , \left(\begin{array}{c} u_R\\rho_R\\ p_R \end{array}\right) =\left(\begin{array}{c} 0.0\\0.125\\0.1 \end{array}\right) $$

具体的な解析内容は以下の通りです。

時間ステップ数 : \(nstep=300\)
時間刻み幅 : \(dt=0.001\)
x方向に取るセル数 : \(nx0=100\)
x方向の刻み幅 : \(dx=0.01\)
初期条件 : \(Q(x)=Q_L (x \leq 0.5), Q(x)=Q_R (x>0.5)\)
境界条件 : \(Q(bound_L)=Q_L, \frac{\partial Q(bound_R)}{\partial x}=0\)
比熱比：\(\gamma=1.4\)

結果になります。
0.3秒まで0.001刻みの時間履歴になります。

続いて\(t=2s\)における解析解と数値解の比較になります。

青：解析解　　赤：数値解

sod shock tubeの解析解は計算しておらず、比較のため重ねただけです。

急襲に値の変化している場所では、なまりが大きく定性的に一致したといった結論が得られます。

境界の計算方法等、コードについて不足分を次の記事にまとめたので、そちらと合わせて確認いただけたらと思います。よろしくお願いします。

だいぶ物理現象を形にできてきている気がして、楽しかったです。

補足を第5-4回に記載しています。

by hide

*コードを見るときは，エディターに貼り付けてみた方が見やすいです。

import numpy
import os
import matplotlib.pyplot as pyplot

# --------------------------
# -- initial value        --
# --------------------------
nstep = 300                  # 時間ステップ数
nx0 = 100                    # 空間ステップ数
dt = 0.001                   # 時間刻み幅
dx = 0.01                    # 空間刻み幅

lpml = 1                     # 仮想境界セル数
nx = nx0+2*lpml              # 総空間セル数

# --定数--
gamma = 1.4    # 比熱比

# -- 出力--
dir_name = "FVS_test"        # 出力フォルダ名

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

def setup():  # 初期値入力
    global x, bol, bor, qf, Qc
    u = [0.0]*nx               # 速度
    rho = [0.0]*nx             # 密度
    p = [0.0]*nx               # 圧力
    e = [0.0]*nx               # エネルギー
    x = [0.0]*nx               # 位置

    """
    e=p/(r-1)+rho*u^2/2
    """
    for i in range(nx):
        u[i] = 0.0

        if i <= nx*0.5:
            rho[i] = 1.0
            p[i] = 1.0

        else:
            rho[i] = 0.125
            p[i] = 0.1

        e[i] = p[i]/(gamma-1)+rho[i]*(u[i]**2)/2
        x[i] = i*dx-dx/2


    bol = [0.0]*3             # 左端仮想セル
    bor = [0.0]*3             # 右端仮想セル
    for j in range(3):
        if j == 0:
            bol[j] = rho[0]
            bor[j] = rho[nx-1]
        elif j == 1:
            bol[j] = u[0]*rho[0]
            bor[j] = u[nx-1]*rho[nx-1]
        elif j == 2:
            bol[j] = e[0]
            bor[j] = e[nx-1]

    qf = [[0.0] * 3 for i in [1] * nx]  # 基本量
    Qc = [[0.0] * 3 for i in [1] * nx]  # 保存量
    for i in range(nx):
        for j in range(3):
            if j == 0:
                qf[i][j] = u[i]
                Qc[i][j] = rho[i]
            elif j == 1:
                qf[i][j] = rho[i]
                Qc[i][j] = u[i]*rho[i]
            elif j == 2:
                qf[i][j] = p[i]
                Qc[i][j] = e[i]


def cal_Q():  # 保存量を計算
    global Qc
    cal_Res()                         # 保存量計算用にまとめる

    Qc = numpy.array(Qc)              # arrayに変換
    lo_R = numpy.array(Res)           # arrayに変換

    for i in range(1, nx-1):
        Qc[i] = Qc[i]-dt/dx*lo_R[i]   # メインの計算

    Qc.tolist()                       # listに戻す
    bound()                           # 境界の計算


def bound():  # 境界の計算
    global Qc
    for i in range(3):
        Qc[0][i] = 2*bol[i]-Qc[1][i]  # 左端境界の計算
        Qc[nx-1][i] = Qc[nx-2][i]     # 右端の計算


def cal_Res():  # 境界フラックスの計算
    global Res

    Res = numpy.array([[0.0] * 3 for i in [1] * nx])
    fvs()                             # FVS法によるフラックスの作成

    for i in range(1, nx-1):
        Res[i] = Fplus[i]-Fplus[i-1]

    Res.tolist()


def fvs():  # FVS法によるフラックスの計算(セル1と2の境界をFplus[1]に格納)
    global Fplus

    Fplus = [[0.0] * 3 for i in [1] * (nx+1)]

    for i in range(0, nx-1):
        # iセルにおけるR,R^-1,Λ,|Λ|
        R, R_inv, Gam, Gam_abs = A_pm(i)
        Ap = numpy.dot((numpy.dot(R, Gam+Gam_abs)), R_inv)               # 固有値が正のものを計算

        # i+1セルにおけるR,R^-1,Λ,|Λ|
        R, R_inv, Gam, Gam_abs = A_pm(i+1)
        Am = numpy.dot((numpy.dot(R, Gam-Gam_abs)), R_inv)               # 固有値が負のものを計算

        Fplus[i] = 0.5*(numpy.dot(Ap, Qc[i]) + numpy.dot(Am, Qc[i+1]))   # フラックスを計算


def A_pm(ite):  # ヤコビアン行列の固有値もろもろ計算
    H = (Qc[ite][2]+qf[ite][2])/Qc[ite][0]  # エンタルピー
    u = qf[ite][0]
    c = numpy.sqrt((gamma-1)*(H-0.5*u**2))
    b_para = (gamma-1)/c**2
    a_para = 0.5*b_para*u**2

    R = numpy.array([[1.0, 1.0, 1.0, ], [u-c, u, u+c],
                     [H-u*c, 0.5*u**2, H+u*c]])
    R_inv = numpy.array([[0.5*(a_para+u/c), 0.5*(-b_para*u-1/c), 0.5*b_para], [
                        1-a_para, b_para*u, -b_para], [0.5*(a_para-u/c), 0.5*(-b_para*u+1/c), 0.5*b_para]])
    Gam = numpy.array([[(u-c), 0.0, 0.0], [0.0, u, 0.0], [0.0, 0.0, (u+c)]])
    Gam_abs = numpy.array(
        [[abs(u-c), 0.0, 0.0], [0.0, abs(u), 0.0], [0.0, 0.0, abs(u+c)]])

    return R, R_inv, Gam, Gam_abs


def qftoQc(qf):  # 基本量から保存量変換
    lo_Qc = [[0.0] * 3 for i in [1] * nx]
    for i in range(nx):
        for j in range(3):
            if j == 0:
                lo_Qc[i][j] = qf[i][1]
            elif j == 1:
                lo_Qc[i][j] = qf[i][1]*qf[i][0]
            elif j == 2:
                lo_Qc[i][j] = (qf[i][2]/(gamma-1) + 1.0/2.0*qf[i][1]*(qf[i][0]**2))

    return lo_Qc


def Qctoqf(Qc):  # 保存量から基本量変換
    lo_qf = [[0.0] * 3 for i in [1] * nx]
    for i in range(nx):
        for j in range(3):
            if j == 0:
                lo_qf[i][j] = Qc[i][1]/Qc[i][0]
            elif j == 1:
                lo_qf[i][j] = Qc[i][0]
            elif j == 2:
                lo_qf[i][j] = (gamma-1)*(Qc[i][2] - 1.0/2.0*Qc[i][0]*((Qc[i][1]/Qc[i][0])**2))
    return lo_qf


def output_q(odir9, f_name9, q9, x9):  # テキスト形式で出力
    outlist = ["x[m] u[m/s] rho[kg/m3] p[Pa]"]  # 出力するものの名前
    for i in range(len(q9)):
        outlist.append(str(x9[i])+" "+str(q9[i][0]) +
                       " "+str(q9[i][1])+" "+str(q9[i][2]))
    outlist = '\n'.join(outlist)

    with open(odir9+"/"+f_name9, 'wt') as f:
        f.write(outlist)


def cre_dir():  # フォルダ作成
    try:
        os.mkdir(dir_name)
    except:
        pass


# --------------------------
# -- main                 --
# --------------------------
cre_dir()  # フォルダ作成
setup()    # 初期値計算

for k in range(nstep):
    print(k)
    cal_Q()               # 保存量計算
    qf = Qctoqf(Qc)       # 保存量から基本量へ変換

    output_q(dir_name, "time"+'{:0=4}'.format(int(k))+"d-3", qf, x)  # 出力