平行平板ポアズイユ流れの熱伝達のプログラム(最終版) [ネコ騙し数学]
! 2平板間の流れ(平板ポアズイユ流れ、壁温一定)を解くプログラム
! 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
! 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
出力部は自分で作るといいにゃ。そこまでは面倒を見きれないケロよ。
mはX方向の分割数、nはY方向の分割数、n1=n+1だケロ。Xの最大値は10だから、m=100、n=10くらいが丁度いいんじゃないかな。mを1000、nを100にしても、m×nオーダーの計算法だから、エンターキーを押した瞬間に終わると思うけれど、U、V、T合せて30万のデータ数になるから、mとnをあまり大きくすのはやめたほうがいいと思うにゃ。
タグ:数値解析
2017-11-19 10:17
nice!(0)
コメント(0)
コメント 0