工学系大学院生のブログ

2人の大学院生による雑記ブログ

第6-5回 高次精度一次元オイラー方程式[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\} $$


第5回との比較を行いたいので、同様の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\)
MUSCL法,3次精度:\(k=1/3\)
圧縮パラメータ:\(b=(3-k)/(1-k)\)

また、MUSCL法は基本量\(q\)を用いて内挿を行い、その後保存量\(Q\)に変換し計算を進めています。



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



続きまして\(t=2.0s\)における解析解との比較になります。
青:解析解  赤:数値解



第5回の結果を下記にのせておきます。
青:解析解  赤:数値解



同じ点数を用いているにも関わらず、急峻な変化も捕らえることができました。計算時間は格子を増やすよりもかからないため、建設的です。


コードが長く、複雑になってきましたが、よろしかったら参考までにお使いください。


ここまで読んでくださり、ありがとうございます。

by hide

*コードはエディターにコピーしてから読んでいただくと見やすいです。

import numpy
import os

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

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

# -- slope limiter --

k_muscl=1/3                        # muscl精度
b_muscl=(3-k_muscl)/(1-k_muscl)

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

# -- 出力--
dir_name="FVS_muscl" # 出力フォルダ名
out_name_front="time"                 # 出力ファイル名(先頭)
out_name_back="d-3"
# --------------------------
# -- 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()               # 右辺Rの計算

    Qc=numpy.array(Qc)       # Qcのarray化
    lo_R=numpy.array(Res)    # 右辺Rのarray化

    Qc_np=numpy.array(Qc[:]) # Qc^nの保存

    for i in range(1,nx-1):           # Qc^(n+1)の計算
        Qc[i]=Qc[i]-dt/dx*lo_R[i] 
    cal_Res()                         # R^(n+1)の計算
    lo_R=numpy.array(Res)             # R^(n+1)のarray化
    for i in range(1,nx-1):  # nとn+1の平均の値を用いて,二次精度のn+1の計算
        Qc[i]=1/2*(Qc_np[i]+Qc[i]-dt/dx*lo_R[i])

    Qc.tolist()
    
    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)]
    muscl()

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

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

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


def A_pm(lQc,lqf):  # ヤコビアン行列の固有値もろもろ計算
    H = (lQc[2]+lqf[2])/lQc[0]  # エンタルピー
    u = lqf[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 muscl():
    global qf,qfL,qfR,QcL,QcR
    # 1と2の間を1に収納

    qf=Qctoqf(Qc)

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

    for i in range(1,nx-2):
        for j in range(3):
            dplus_j=qf[i+1][j]-qf[i][j]
            dminus_j=qf[i][j]-qf[i-1][j]
            dplus_jp=qf[i+2][j]-qf[i+1][j]
            dminus_jp=qf[i+1][j]-qf[i][j]
            
            qfL[i][j]=qf[i][j]+1/4*((1-k_muscl)*minmod(dminus_j,dplus_j,b_muscl)+(1+k_muscl)*minmod(dplus_j,dminus_j,b_muscl))
            qfR[i][j]=qf[i+1][j]-1/4*((1-k_muscl)*minmod(dplus_jp,dminus_jp,b_muscl)+(1+k_muscl)*minmod(dminus_jp,dplus_jp,b_muscl))

    # 境界内側用
    for j in range(3):
        dplus_jp=qf[2][j]-qf[1][j]
        dminus_jp=qf[1][j]-qf[0][j]
        qfR[0][j]=qf[1][j]-1/4*((1-k_muscl)*minmod(dplus_jp,dminus_jp,b_muscl)+(1+k_muscl)*minmod(dminus_jp,dplus_jp,b_muscl))

        dplus_j=qf[nx-1][j]-qf[nx-2][j]
        dminus_j=qf[nx-2][j]-qf[nx-3][j]
        qfL[nx-2][j]=qf[nx-2][j]+1/4*((1-k_muscl)*minmod(dminus_j,dplus_j,b_muscl)+(1+k_muscl)*minmod(dplus_j,dminus_j,b_muscl))
        
            
    QcL=qftoQc(qfL)
    QcR=qftoQc(qfR)

    # 境界外側用(境界は風上)
    qfL[0]=qf[0][:]
    QcL[0]=Qc[0][:]
    qfR[nx-2]=qf[nx-1][:]
    QcR[nx-2]=Qc[nx-1][:]
    

def minmod(x,y,b):
    ans=numpy.sign(x)*max(0,min(abs(x),numpy.sign(x)*y*b))
        
    return ans

    
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(f_name):  # テキスト形式で出力
    outlist=["x[m] u[m/s] rho[kg/m3] p[Pa]"]  # 出力するものの名前
    for i in range(len(qf)):
        outlist.append(str(x[i])+" "+str(qf[i][0])+" "+str(qf[i][1])+" "+str(qf[i][2]))
    outlist='\n'.join(outlist)

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

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

# --------------------------
# -- preparetion          --
# --------------------------
cre_dir()
setup()
# --------------------------
# -- main                 --
# --------------------------

for k in range(nstep):
    
    print(k)
    cal_Q()
    qf=Qctoqf(Qc)
        
    output_q(out_name_front+'{:0=4}'.format(int(k*dt*1000))+out_name_back)

このエントリーをはてなブックマークに追加

コメント

  • どのような境界条件を課していますか?どのようにそれを実装していますか?

  • ご質問ありがとうございます。

    左端では初期条件u=0.0,rho=1.0,p=1.0を境界条件としています。
    実装時には、まず外側に一つセル(仮想セル)を作ります。この仮想セルをセル[0]、一つ内側の計算セルをセル[1]とします。これらのセル境界で上記の条件を維持したいので
    (セル[0]+セル[1])/2=境界条件 ⇔セル[0]=2*境界条件-セル[1] (コード109行目)

    右端では流出条件として勾配0の条件を与えています。
    実装時には、同様に仮想セルを用意します。この仮想セルをセル[n-1]、一つ内側の計算セルをセル[n-2]とします。勾配0なので
    (セル[n-1]-セル[n-2])/⊿x =0 ⇔ セル[n-1]=セル[n-2]

    ホームページトップ「このブログについて」
    ⇒「第5-4回 オイラー方程式(コードの補足)[python]
    ⇒①仮想セルについて
    でも説明していますので、よろしければご覧ください。

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です

トラックバックURL: https://teru-hide.com/cfd-python-shock-tube-high-resolution/trackback/
にほんブログ村 英語ブログへ
にほんブログ村