2次元圧縮性流体計算（無反射境界条件）

テキスト「空力音の直接数値計算を支える技術」（暗黙の型宣言著）の補足課題7

全部読んだらやってみよう！その2

課題

補足課題6で作成した圧縮性流体計算プログラムに、無反射境界条件を適用する．

ヒント
- 特性波の方向に注意．入ってくる特性波を0にして打ち消す．

解答例

プログラム全体は載せられないので，ポイントだけ解説します．

main.f90

        TVDRK : block
                !   :
                !  略
                !   :
                !--- 時間積分の前にNSCBCを適用
                ! NSCBC x-dir
                bd = ([1, Nx])

                do j = 1, Ny
                    ! x=0側の境界．
                    !入ってくる特性波（つまり流速が正の特性波）を0にして打ち消す．
                    L11(1,j) = max((u(1,j) - C0),0.d0)*(d_p_dx(1,j) - dens(1,j)*C0*d_u_dx(1,j))   
                    ! L11の特性波の流速は(u - C0)

                    L12(1,j) = max( u(1,j)      ,0.d0)*((C0**2)*d_dens_dx(1,j) - d_p_dx(1,j))        
                    ! L12の特性波の流速はu

                    L13(1,j) = max( u(1,j)      ,0.d0)*d_v_dx(1,j)                                                   
                    ! L13の特性波の流速はu

                    L15(1,j) = max((u(1,j) + C0),0.d0)*(d_p_dx(1,j) + dens(1,j)*C0*d_u_dx(1,j)) 
                    ! L15の特性波の流速は(u + C0)

                   ! x=Lx側の境界．
                   ! 入ってくる特性波（つまり流速が負の特性波）を0にして打ち消す．
                    L11(2,j) = min((u(Nx,j) - C0),0.d0)*(d_p_dx(Nx,j) - dens(Nx,j)*C0*d_u_dx(Nx,j))    
                    L12(2,j) = min( u(Nx,j)      ,0.d0)*((C0**2)*d_dens_dx(Nx,j) - d_p_dx(Nx,j))
                    L13(2,j) = min( u(Nx,j)      ,0.d0)*d_v_dx(Nx,j)
                    L15(2,j) = min((u(Nx,j) + C0),0.d0)*(d_p_dx(Nx,j) + dens(Nx,j)*C0*d_u_dx(Nx,j))
                end do
                
                do nb = 1, 2
                    d11 = (L12(nb,:) + 0.5d0*(L15(nb,:) + L11(nb,:)))/C0**2
                    d12 = (L15(nb,:) + L11(nb,:))*0.5d0
                    d13 = (L15(nb,:) - L11(nb,:))*0.5d0/(dens(bd(nb),:)*C0)
                    d14 =  L13(nb,:)

                   ! NSCBCを適用した支配方程式(テキストの式(6.1)~(6.5))と，
                   ! もともとの支配方程式(テキストの式(2.1)~(2.5))を比較し，
                   ! 境界の値にNSCBCを適用する．
                    d_dens_u_dx(bd(nb),:)   = d11(:)
                    d_dens_u_u_dx(bd(nb),:) = u(bd(nb),:)*d11(:) + dens(bd(nb),:)*d13(:)
                    d_p_dx(bd(nb), :)  = 0.d0
                    d_visyx_dx(bd(nb), :) = 0.d0
                    d_dens_u_v_dx(bd(nb),:) = v(bd(nb),:)*d11(:) + dens(bd(nb),:)*d14(:)
                    d_dens_u_TE_dx(bd(nb),:)= 0.5d0*(u(bd(nb),:)**2 + v(bd(nb),:)**2)*d11(:) + d12(:)/(Gamma-1.d0) &
                                            + dens(bd(nb),:)*u(bd(nb),:)*d13(:) + dens(bd(nb),:)*v(bd(nb),:)*d14(:)
                    d_u_p_dx(bd(nb),:) = 0.d0
                    d_qx_dx(bd(nb),:)  = 0.d0
                end do

                ! NSCBC y-dir
                !   :
                !  略
                !   :
                
                !--- Time Integration
                !   :
                !  略
                !   :                
        end block TVDRK