２次元移流方程式

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

7章まで読んだらやってみよう！その2

今回は２次元の課題にトライ． 8章に2次元のテスト問題が書かれているけれど，いきなりナビエ・ストークス方程式を解き始めるよりは，もう少し単純な課題で２次元計算の基礎を知っておいたほうがいいと思う．

なお，本課題では別冊テキスト「流体計算で覚えるPython3」（暗黙の型宣言著）に掲載されているテスト問題を使用します．

課題

「流体計算で覚えるPython3」７章に書かれている回転円錐問題をFortranで解いてみよう．

• 支配方程式は，２次元移流方程式．（7.1節に掲載）
• 初期条件は，2次元ベル型関数．（7.5節に掲載）
• ヒント
  • 補足課題3「時間積分」で作成したプログラム一式をベースに２次元化すればOK．

解答例

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

まずはmain.f90のInitializeブロックで初期条件を定義． 移流速度c_xはy，c_yはxで決まる，という点に注意．うっかり逆にしないこと！

    Initialize : block
                 :
                c_x(i,j) = -2.d0*PI*y(j)     
                c_y(i,j) =  2.d0*PI*x(i)   
                 :
　　　call fileout_2dim(file_init, x, y, c_x, c_y, f%var) !初期条件出力
    end block Initialize

初期条件を定義できたら，いったん出力してみて確認． 図1のt/T=0のようになっていればOK．
次に，２次元移流方程式を解くためのコードを書く．ここでは重要なポイントだけ紹介．

• main.f90

df/dx, df/dyをそれぞれ計算し，fを求める．というのを繰り返し，時間積分していく．

    do n = 1, Nt
        TVDRK : block
            integer :: i, j
            call cacheTVDRK_2D(f%var)
            do step=1,3
                do j=2, Ny-1
                    call compute1stDerivative(f%dx(:,j), f%var(:,j), dx, Gx)
                end do
                do i=2, Nx-1
                    call compute1stDerivative(f%dy(i,:), f%var(i,:), dy, Gy)
                end do
                call integrateTVDRK_2D_Adv(f%var, f%dx, f%dy, dt, step, c_x, c_y)
                ! 境界でゼロを強制
                f%var(1 ,: ) = 0d0
                f%var(Nx,: ) = 0d0
                f%var(: ,1 ) = 0d0
                f%var(: ,Ny) = 0d0         
            enddo
        end block TVDRK
    end do

• module_TimeIntegration.f90

２次元移流方程式を時間方向に離散化して，次ステップのfを求める． ここではTVD Runge-Kutta法で時間積分．

    subroutine integrateTVDRK_2D_Adv(f, d_f_dx, d_f_dy, dt, step, conv_x, conv_y)
        implicit none
        real(8), intent(inout) :: f(:,:)
        real(8), intent(in   ) :: d_f_dx(:,:), d_f_dy(:,:)
        real(8), intent(in   ) :: dt
        integer, intent(in   ) :: step
        real(8), intent(in   ) :: conv_x(:,:), conv_y(:,:)
        real(8), parameter :: a(3) = [0d0, 3d0/4d0, 1d0/3d0]
        real(8), parameter :: b(3) = [1d0, 1d0/4d0, 2d0/3d0]

        integer :: Nx, Ny
        integer :: i, j

        Nx = size(f,1)
        Ny = size(f,2)

        do j=1, Ny
            do i=1, Nx
                f(i,j) = a(step)*fn2(i,j) + b(step)*(f(i,j)-conv_x(i,j)*dt*d_f_dx(i,j)-conv_y(i,j)*dt*d_f_dy(i,j))
            end do
        end do

    end subroutine integrateTVDRK_2D

ひととおりプログラムができたら，円錐が1/4周，2/4周，3/4周，1周したときの計算結果を出力し，正しく反時計回りにまわっていることを確認． 可視化ソフトは「AVえせ」というフリーウェアを使用．# あやしいソフトじゃないよ！

Fig.1. Time variation of rotating cone

格子点数やクーラン数を変えてみて，1周したときの計算結果を可視化して確認する．

まず，総格子点数を変えてみた場合．可視化結果にほとんどちがいはみられない．（格子点数の分だけベクトルが表示される）

Fig.2. Rotating cone at t/T=4/4. (CFL = 0.04)

次に，クーラン数(CFL)を変えてみた場合．こちらも，ほとんどちがいはみられない．

Fig.3. Rotating cone at t/T=4/4. (N=1600)

1周したときの領域全体のエラーを出力する．

    print *, "CFL:", CFL, "error:", sum(abs(f%var-f0%var))/(Nx*Ny) ! f%varは計算結果．f0%varは初期値．

格子点数とクーラン数を変えたときのエラーをグラフにすると，

Fig.4. Error

クーラン数よりも格子点数によってエラーが決まっている．