2次元圧縮性流体計算（単極渦）

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

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

課題

いよいよ圧縮性流体計算！

8.2章に書かれている例題（超音速一様流で移流する単極渦まわりの流れ場）に取り組もう．

ヒント
- 補足課題5で作成したプログラム一式をベースに，支配方程式を質量，運動量，全エネルギの保存則に変更する．
- どの物理量をどういう順番で解くか，よく考えること．特に粘性応力の計算が難しい(><)
- 変数が膨大にあるので，ケアレスミスに注意！
- 絶対にミスは起こるので、どういうふうにデバッグしていくか・・を考えるのも，この問題の目的．

解答例

空間微分，時間積分はこれまでの補足課題と同様に以下の条件とします．

空間微分 6次精度コンパクト差分（境界では4次，境界1点内側では5次精度）
時間積分 TVDルンゲ・クッタ法

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

main.f90

    do n = 1, Nt
        TVDRK : block
            integer :: i, j
            dens_u = dens*u
            dens_v = dens*v
            dens_TE = dens*TE
           
           ! TVDRKで時間進行させる保存量をキャッシュに入れる
            call cacheTVDRK_2D_NS(dens, dens_u, dens_v, dens_TE)  
            do step=1,3
                do j=1, Ny
                    ! ひたすら空間微分
                    call compute1stDerivative(d_dens_u_dx(:,j), dens(:,j)*u(:,j), dx, Gx)
                    call compute1stDerivative(d_dens_u_u_dx(:,j), dens(:,j)*u(:,j)*u(:,j), dx, Gx)
                    call compute1stDerivative(d_p_dx(:,j), p(:,j), dx, Gx)
                    call compute1stDerivative(d_dens_u_v_dx(:,j), dens(:,j)*u(:,j)*v(:,j), dx, Gx)
                    call compute1stDerivative(d_dens_u_TE_dx(:,j), dens(:,j)*u(:,j)*TE(:,j), dx, Gx)
                    call compute1stDerivative(d_u_p_dx(:,j), u(:,j)*p(:,j), dx, Gx)
                    call compute1stDerivative(d_temp_dx(:,j), temp(:,j), dx, Gx)
                    kp(:,j) = mu(:,j)*Cp/Pr
                    qx(:,j) = -kp(:,j)*d_temp_dx(:,j)
                    call compute1stDerivative(d_qx_dx(:,j), qx(:,j), dx, Gx)

                    call compute1stDerivative(d_u_dx(:,j), u(:,j), dx, Gx)
                    call compute1stDerivative(d_v_dx(:,j), v(:,j), dx, Gx)
                end do
                do i=1, Nx
                    call compute1stDerivative(d_dens_v_dy(i,:), dens(i,:)*v(i,:), dy, Gy)
                   !   :
                   ! 省略
                   !  :
                end do

                ! 粘性応力は，x方向とy方向両方の微分から決まるため，ひととおりx方向とy方向の微分値を求めた後に計算．
                visxx = mu*(2.d0*d_u_dx - (2.d0/3.d0)*(d_u_dx + d_v_dy))
                visxy = mu*(d_u_dy + d_v_dx)
                visyx = mu*(d_v_dx + d_u_dy)
                visyy = mu*(2.d0*d_v_dy - (2.d0/3.d0)*(d_u_dx + d_v_dy))
     
                ! 粘性応力の微分．もっと厳密にやるなら，テキスト8.1章「2階微分の計算法」に書かれているやり方でやりましょう．
                do j=1, Ny
                    call compute1stDerivative(d_visxx_dx(:,j), visxx(:,j), dx, Gx)
                    call compute1stDerivative(d_visyx_dx(:,j), visyx(:,j), dx, Gx)
                    call compute1stDerivative(d_U_visx_dx(:,j), u(:,j)*visxx(:,j) + v(:,j)*visyx(:,j), dx, Gx)
                end do
                do i=1, Nx
                    !   :
                    ! 省略
                    !   :
                end do

               ! 時間積分
                call integrateTVDRK_2D_NS(dens, -(d_dens_u_dx + d_dens_v_dy), dt, step, 1) !OK
                call integrateTVDRK_2D_NS(dens_u, &
                                        -(d_dens_u_u_dx + d_dens_v_u_dy) - d_p_dx + d_visxx_dx + d_visxy_dy, &
                                          dt, step, 2)
                call integrateTVDRK_2D_NS(dens_v, &
                                        -(d_dens_u_v_dx + d_dens_v_v_dy) - d_p_dy + d_visyx_dx + d_visyy_dy, &
                                          dt, step, 3)
                call integrateTVDRK_2D_NS(dens_TE, &
                                        -(d_dens_u_TE_dx + d_dens_v_TE_dy + d_u_p_dx + d_v_p_dy) &
                                        + d_U_visx_dx + d_U_visy_dy - (d_qx_dx + d_qy_dy), dt, step, 4)

               ! 各物理量を計算．これを次の時刻で空間微分する．
                u = dens_u/dens
                v = dens_v/dens
                TE = dens_TE/dens
                dens_e = dens_TE - 0.5d0*(dens*u*u + dens*v*v)
                e = dens_e/dens
                p = dens_e*(Gamma-1.d0)
                temp = e/Cv
                mu = Mu0*((Temp0 + Temp1)/(temp + Temp1))*((temp/Temp0)**(3/2)) 
     
            enddo
        end block TVDRK
    end do

module_TimeIntegration.f90

    ! TVDRKで時間進行させる保存量をキャッシュに入れる
    subroutine cacheTVDRK_2D_NS(dens, dens_u, dens_v, dens_TE)
        real(8), intent(in) :: dens(:,:), dens_u(:,:), dens_v(:,:), dens_TE(:,:)
        dens2(:,:) = dens(:,:)
        dens_u2(:,:) = dens_u(:,:)
        dens_v2(:,:) = dens_v(:,:)
        dens_TE2(:,:) = dens_TE(:,:)
    end subroutine cacheTVDRK_2D_NS

   ! TVDRKでの時間積分は，こんなふうに書いてみた．
   ! けっきょくどの物理量も式①で積分することは同じなので，コードの繰り返しを避けるため，被積分配列をcase文で選択．
    subroutine integrateTVDRK_2D_NS(f, d_f_dt, dt, step, n)
        implicit none
        real(8), intent(inout) :: f(:,:)
        real(8), intent(in   ) :: d_f_dt(:,:)
        real(8), intent(in   ) :: dt
        integer, intent(in   ) :: step
        integer, intent(in   ) :: n
        real(8), parameter :: a(3) = [0d0, 3d0/4d0, 1d0/3d0]
        real(8), parameter :: b(3) = [1d0, 1d0/4d0, 2d0/3d0]

        select case (n)
        case (1)
            fn2 = dens2
        case (2)
            fn2 = dens_u2
        case (3)
            fn2 = dens_v2
        case (4)
            fn2 = dens_TE2    
        end select  

        f(:,:) = a(step)*fn2(:,:) + b(step)*(f(:,:)-dt*d_f_dt(:,:))  ! 式①

    end subroutine integrateTVDRK_2D_NS