平行平板ポアズイユ流れの熱伝達のプログラム（最終版）

! ２平板間の流れ（平板ポアズイユ流れ、壁温一定）を解くプログラム
! NS eq. U(∂U/∂X)+V(∂U/∂Y)=-dP/dX+1/Re(∂²U/∂²Y)
! 連続の式 ∂U/∂X+∂V/∂Y=0
! エネルギー方程式 U(∂T/∂X)+V(∂T/∂Y)=1/(RePr)(∂²U/∂²Y)
parameter (m=10,n=5,n1=6) ! n1=n+1
real u(0:m,0:n),v(0:m,0:n), dpdx(m)
real a(n),b(n),c(n),d(n+1),e(n),f(n+1)
real t(0:m,0:n)
real tb(m), qnu(m)
real dummy_v(n);
iter_limit = 200 ! 流れ場の反復計算の反復回数の上限
eps = 1.0e-6 ! 流れ場の収束判定

xl=10
re=100.; pr =0.7
dy=0.5/n
dy2=dy*dy
dx=xl/m

! 速度u,vの初期化
do i=0,m
    do j=0,n
        u(i,j) =0.0 ; v(i,j)=0
    end do
end do

! 流入速度(X=0)の速度
do j=1,n
    u(0,j)=1.0
end do

do i=0, m ! 温度の初期化
    do j=1,n
        t(i,j)=0
    end do
end do

do i=1,m ! 壁の温度
    t(i,0)=1
end do


do i=1,m
    x= x+dx
    do iter=1, iter_limit
    ! NS方程式の係数行列の計算と値のセット
        do j=1,n
            a(j)=-v(i,j)/(2.*dy)-1./(re*dy2)
            b(j)=u(i,j)/dx+2./(re*dy2)
            if (j.ne.n) then
                c(j)=v(i,j)/(2*dy) -1./(re*dy2)
            else
                c(j)=0.
            end if
            d(j)=u(i-1,j)/dx*u(i,j)
            e(j)=1.
            f(j)=dy
        end do
        a(n)=a(n)-1./(re*dy2)
        d(n1)=0.5
        f(n)=dy/2
        f(n1)=0.

        ! 速度と圧力を解くプログラムを呼び出す
        call calcup(a,b,c,d,e,f,n,n1)
        do j=1,n
            u(i,j)=d(j)
        end do
        err=abs(1-dpdx(i)/d(n1))
        dpdx(i)=d(n1)
       
        ! 連続の式からvを計算　（ここは改良すべき・・・）
        do j=1,n-1
!            v(i,j) = v(i,j-1)+dy/dx*(u(i-1,j)-u(i,j))
            um=0.5*(u(i-1,j) + u(i-1,j-1))
            up=0.5*(u(i,j) + u(i,j-1))
            v(i,j)=v(i,j-1)+dy/dx*(um-up) ! 少し改良・・・
        end do
        if (err.lt.eps) exit
    end do
end do


! エネルギー方程式の連立方程式の係数の計算
do i=1, m
    do j=1,n
        a(j)=-v(i,j)/(2.*dy)-1./(re*pr*dy2)
        b(j)=u(i,j)/dx+2./(re*pr*dy2)
        if (j.ne.n) then
            c(j)=v(i,j)/(2*dy) -1./(re*pr*dy2)
        else
            c(j)=0.
        end if
        d(j)= u(i,j)/dx*t(i-1,j)
    end do
    d(1)=d(1)+1./(re*pr*dy2)*t(i,0)
    a(1)=0.
    a(n)=2*a(n)
    c(n)=0.

    call tdma(a,b,c,d,n) ! 連立方程式を解く
   
    do j=1,n
        t(i,j)=d(j)
    end do
end do

! 台形公式でバルク温度（混合平均温度）の計算
    do i=1,m
        s=0.5*u(i,n)*t(i,n)
        do j =1,n-1
            s = s + u(i,j)*t(i,j)
        end do
        tb(i)=2*s*dy;
    end do

! 局所の熱伝導量、並びに局所ヌセルト数を計算
    do i=1, m
        dtdy= (4*t(i,1)-t(i,2)-3*t(i,0))/(2*dy)
        qnu(i)=-dtdy/(1-tb(i))
    end do

write(6,*) ' ***** 計算結果 ***** '
write(6,*) 'U'
do i=1, m
    write(6,100) i*dx,(u(i,j),j=1,n), dpdx(i)
end do
write(6,*)
write(6,*) 'V'
do i=1, m
    write(6,100) i*dx,(v(i,j),j=1,n)
end do

write(6,*)

write(6,*)
write(6,*) 'T'
do i=1, m
    write(6,100) i*dx,(t(i,j),j=1,n), qnu(i)
end do

100 format(f8.5,1x, 12(f8.5,1x))
110 format(a, 1x, 5(f8.5,1x),a)

end

! UとdP/dXを解くサブルーティン
subroutine calcup(a,b,c,d,e,f,n,n1)
real a(n), b(n), c(n),d(n1),e(n) ,f(n1)

! 前進消去
do i=1,n-1
    ratio=a(i+1)/b(i)
    b(i+1)=b(i+1)-ratio*c(i)
    d(i+1)=d(i+1)-ratio*d(i)
    d(n1)=d(n1)-f(i)/b(i)*d(i)
    e(i+1)=e(i+1)-ratio*e(i)
    f(i+1)=f(i+1)-f(i)/b(i)*c(i)
    f(n1)=f(n1)-f(i)/b(i)*e(i)
end do
f(n1)=f(n1)-f(n)/b(n)*e(i)
d(n1)=d(n1)-f(n)/b(n)*d(n)

! 後退代入
d(n1)=d(n1)/f(n1)
d(n)=(d(n)-e(n)*d(n1))/b(n)
do i=n-1,1,-1
    d(i)=(d(i)-c(i)*d(i+1)-e(i)*d(n1))/b(i)
end do
end

!  TDMA
subroutine tdma(a,b,c,d,n)
real a(n), b(n), c(n),d(n)

! 前進消去
do i=1,n-1
    ratio=a(i+1)/b(i)
    b(i+1)=b(i+1)-ratio*c(i)
    d(i+1)=d(i+1)-ratio*d(i)
end do

! 後退代入
d(n)=d(n)/b(n)
do i=n-1,1,-1
    d(i)=(d(i)-c(i)*d(i+1))/b(i)
end do

end

出力部は自分で作るといいにゃ。そこまでは面倒を見きれないケロよ。

ｍはX方向の分割数、nはY方向の分割数、n1=n+1だケロ。Xの最大値は１０だから、m=100、n=10くらいが丁度いいんじゃないかな。mを１０００、nを１００にしても、m×nオーダーの計算法だから、エンターキーを押した瞬間に終わると思うけれど、U、V、T合せて３０万のデータ数になるから、mとnをあまり大きくすのはやめたほうがいいと思うにゃ。

タグ：数値解析