２次元拡散方程式

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

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

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

課題

2次元拡散方程式をFortranで解いてみよう．

• 支配方程式は，２次元拡散方程式．
• 初期条件は，2次元ベル型関数．（別冊テキスト7.7節と同様）
• ヒント
  • 補足課題4「2次元移流方程式」で作成したプログラム一式をベースに，支配方程式を拡散方程式に変更すればOK．

解答例

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

拡散方程式では，空間微分が1階微分ではなく２階微分になるので，これを計算するコードに修正． さらに，別冊テキストの離散化式(7.9)をTVD Runge-Kutta法で時間積分するコードを追加．

２次元移流方程式から２次元拡散方程式への変更はこれくらい．一瞬でできてしまう．

• main.f90

    do n = 1, Nt
        TVDRK : block
            integer :: i, j
            call cacheTVDRK_2D(f%var)
            do step=1,3
                do j=2, Ny-1
                    call compute2ndDerivative(f%dx2(:,j), f%var(:,j), dx, Gx) ! x方向の２階微分を計算
                end do
                do i=2, Nx-1
                    call compute2ndDerivative(f%dy2(i,:), f%var(i,:), dy, Gy) ! y方向の２階微分を計算
                end do
                ! 境界でゼロを強制
                f%var(1 ,: ) = 0d0
                f%var(Nx,: ) = 0d0
                f%var(: ,1 ) = 0d0
                f%var(: ,Ny) = 0d0
                call integrateTVDRK_2D_DiffusionEq(f%var, f%dx2, f%dy2, dt, step, mu) ! TVDRK法で拡散方程式の時間積分するサブルーチンを呼び出す
            enddo
        end block TVDRK
    end do

ここまでの補足課題をこなしていれば，サブルーチン「compute2ndDerivative」は module_CompactScheme.f90 で既に定義してあるはず． もし定義していなければ，テキストを参考に定義してください．

• module_TimeIntegration.f90

    subroutine integrateTVDRK_2D_DiffusionEq(f, d_f_dx2, d_f_dy2, dt, step, mu)
        implicit none
        real(8), intent(inout) :: f(:,:)
        real(8), intent(in   ) :: d_f_dx2(:,:), d_f_dy2(:,:)
        real(8), intent(in   ) :: dt
        integer, intent(in   ) :: step
        real(8), intent(in   ) :: mu
        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) + mu*dt*(d_f_dx2(i,j) + d_f_dy2(i,j))) ! TVDRK法で拡散方程式の時間積分
            end do
        end do
    end subroutine integrateTVDRK_2D_DiffusionEq

拡散数0.25の計算結果は図1のとおり．時間とともに円錐が拡散している． もっとしっかりやるなら，厳密解（２次元ガウス分布）との誤差も計算するとよいでしょう．

Fig.1. Spatial distribution around cone. (CFL=0.04, N=1600)