shock tube::衝撃波管

飛行ロボット大会も大きな大きな悔しさを残しながら終わったので(24日の日記参照)、shocktubeのプログラムを書くことにした。昨日の飲み会前に終わらせたいなと思っていたのに歯医者とかあって終わらず、結局今日になって書くのは終わった。結果を表示させてみたらやっぱりどう考えてもおかしい。これからバグ取り無間地獄。

1次元shocktubeをYee 2nd-order upwind TVD schemeで解いた。コードは以下の通り。計算は一瞬で終わった。ubuntuでtimeコマンドで時間測ってみたらだいたい0.3秒くらい。一回一回速度違うんだなと思った。コンパイラはg95。なんでgfortranじゃないかっていうとgfortranはVistaにインストールできなかったから。

  • base.f90
program shocktube
implicit none
integer,parameter :: kiza_x=200
!kiza_xは空間刻み数で半分の位置が初期膜境界:偶数である必要がある
real,parameter :: dt=0.000001
!dtは時間ステップ
real,parameter :: cv=1007.0,r=287.0,gamma=1.4
!cvは定圧比熱、rは空気の気体定数,gammaは比熱比
real :: dx
!dxは空間刻み距離
real :: q(3,kiza_x),e(3,kiza_x),u(kiza_x),h(kiza_x),ee(3,kiza_x-1),rho(kiza_x),p(kiza_x),t(kiza_x),ene(kiza_x),c(3,kiza_x)
!real :: t(kiza_x)
real :: right(3,3,kiza_x-1),left(3,3,kiza_x-1),b(2),alpha(3,kiza_x-1),sigma(3,kiza_x),g(3,kiza_x-1),gam(3,kiza_x),phi(3,kiza_x)
real :: psi,sig,minmod
integer :: i,j,k
!q(物理量、場所)は保存量Q(密度、運動量、内部エネルギ)
!e(物理量、場所)は流速ベクトルE(運動量、p+ρu^2,(e+p)u)
!u(場所)は場所によるx方向速度
!h(場所)は全エンタルピ
!ee(物理量、場所)は数値流束ベクトル。一つ右の値を表す
!rho(場所)は密度ρ
!p(場所)は静圧p
!ene(場所)は全エネルギe
!right()は右固有ベクトル
!left()は左固有ベクトル
!b()は固有ベクトル計算の定数
!alpha()は各特性線をまたいだ物理量の跳び
!sigma()は計算に使うやつ
!gはminmodが中に入っててαを比較
!gamはγの演算
!phiは流束ベクトル計算に使うベクトル
!c()は固有値ベクトル
!psiはΨの関数
!sigはσの関数
!minmodはminmod limiter
real :: rho_av(kiza_x-1),u_av(kiza_x-1),h_av(kiza_x-1),c_av(kiza_x-1)
!rho_av()はRoe平均密度
!u_av()はRoe平均速度
!h_av()はRoe平均全エンタルピ
!c_av()はRoe平均特性速度
dx=0.2/kiza_x
!dxは一個あたりの刻み幅
!print*, dx
!計算開始
!do i=1,kiza_x
p(1:kiza_x/2)=5.0*10**5
p(kiza_x/2+1:kiza_x)=1.0*10**5
!静圧の初期化
u(1:kiza_x)=0.0
!x方向速度の初期化
rho(1:kiza_x/2)=6.0
rho(kiza_x/2+1:kiza_x)=1.0
!密度の初期化
do i=1,kiza_x
	q(1,i)=rho(i)
	q(2,i)=rho(i)*u(i)
	t(i)=p(i)/(rho(i)*r)
	q(3,i)=rho(i)*(cv*t(i)+0.5*u(i)**2)
	ene(i)=q(3,i)
!	print*, q(1,i),q(2,i),q(3,i)
end do
!保存量の初期化:3つめの全エネルギが合っているかあやしい。密度×定圧比熱×静温なんだけど。
!eneに全エネルギを代入してる。
do k=1,100
do i=1,kiza_x
	e(1,i)=rho(i)*u(i)
	e(2,i)=p(i)+rho(i)*u(i)**2
	e(3,i)=(ene(i)+p(i))*u(i)
	
!	print*, e(1,i),e(2,i),e(3,i)
end do
!流束ベクトルの計算
do i=1,kiza_x
	h(i)=ene(i)+p(i)
end do
!全エンタルピの計算
do i=1,kiza_x-1
	rho_av(i)=sqrt(rho(i)*rho(i+1))
	u_av(i)=(sqrt(rho(i))*u(i)+sqrt(rho(i+1))*u(i+1))/(sqrt(rho(i))+sqrt(rho(i+1)))
	h_av(i)=(sqrt(rho(i))*h(i)+sqrt(rho(i+1))*h(i+1))/(sqrt(rho(i))+sqrt(rho(i+1)))
	c_av(i)=sqrt((gamma-1)*(h_av(i)-0.5*u_av(i)**2))
!print*, rho_av(i),u_av(i),h_av(i),c_av(i)
end do
!Roe平均の密度、速度、エンタルピ,特性速度の計算
do i=1,kiza_x-1
	c(1,i)=u_av(i)-c_av(i)
	c(2,i)=u_av(i)
	c(3,i)=u_av(i)+c_av(i)
end do
!固有値ベクトル
do i=1,kiza_x-1
right(1,1:3,i)=1.0
right(2,1,i)=u_av(i)-c_av(i)
right(2,2,i)=u_av(i)
right(2,3,i)=u_av(i)+c_av(i)
right(3,1,i)=h_av(i)-u_av(i)*c_av(i)
right(3,2,i)=0.5*u_av(i)**2
right(3,3,i)=h_av(i)+u_av(i)*c_av(i)
end do
!右固有ベクトルの計算
do i=1,kiza_x-1
b(1)=(u_av(i)**2)/2.0*(gamma-1)/(c_av(i)**2)
b(2)=(gamma-1)/(c_av(i)**2)
left(1,1,i)=0.5*(b(1)+u_av(i)/c_av(i))
left(1,2,i)=-0.5*(1/c_av(i)+b(2)*u_av(i))
left(1,3,i)=0.5*b(2)
left(2,1,i)=1-b(1)
left(2,2,i)=b(2)*u_av(i)
left(2,3,i)=-b(2)
left(3,1,i)=0.5*(b(1)-u_av(i)/c_av(i))
left(3,2,i)=0.5*(1/c_av(i)-b(2)*u_av(i))
left(3,3,i)=0.5*b(2)
end do
!左固有ベクトルの計算
do i=1,kiza_x-1
!alpha(3,i)=q(1:3,i+1)-q(1:3,i)
alpha(1:3,i)=matmul(left(1:3,1:3,i),(q(1:3,i+1)-q(1:3,i)))
end do
!物理量の跳びの計算。matmul()を使用
do i=1,kiza_x-1
sigma(1,i)=sig((u_av(i)-c_av(i)),psi(u_av(i)-c_av(i)),dt,dx)
sigma(2,i)=sig(u_av(i),psi(u_av(i)),dt,dx)
sigma(3,i)=sig((u_av(i)+c_av(i)),psi(u_av(i)+c_av(i)),dt,dx)
end do
!σの計算
do i=2,kiza_x-1
	do j=1,3
		g(j,i)=minmod(alpha(j,i),alpha(j,i-1))
	end do
end do
!gの計算。最初の一つめの配列は未使用
do i=2,kiza_x-2
	do j=1,3
		if(alpha(j,i) /= 0) then
			gam(j,i)=sigma(j,i)*(g(j,i+1)-g(j,i))/alpha(j,i)
		else
			gam(j,i)=sigma(j,i)*0
		end if
	end do
end do
!gam()の計算
do i=2,kiza_x-2
	do j=1,3
		phi(j,i)=sigma(j,i)*(g(j,i+1)-g(j,i))-psi(c(j,i)+gam(j,i))*alpha(j,i)
	end do
end do
!Φの計算
do i=2,kiza_x-2
	ee(1:3,i)=0.5*((e(1:3,i+1)+e(1:3,i))+matmul(right(1:3,1:3,i),phi(1:3,i)))
end do
!数値流束ee()の計算
do i=3,kiza_x-2
	q(1:3,i)=q(1:3,i)-dt/dx*(ee(1:3,i)-ee(1:3,i-1))
end do
!次ステップのQの計算
open(10,file='result.txt')
do i=1,kiza_x
	rho(i)=q(1,i)
	u(i)=q(2,i)/q(1,i)
	p(i)=(gamma-1)*(q(3,i)-0.5*q(1,i)*q(2,i))
	t(i)=p(i)/(rho(i)*r)
	write(10,*) rho(i),'	',u(i),'	',p(i),'	',t(i)
end do
close(10)
!q()から物理量を算出している
end do
end
real function psi(z)
	real :: z,delta
	delta=10
	if(abs(z) >= delta) then
		psi=abs(z)
	else
		psi=(z**2+delta**2)/(2.0*delta)
	end if
end
!Ψの関数
real function sig(z,psi,dt,dx)
	real :: z,psi,dt,dx
	sig=0.5*(psi-dt/dx*z**2)
end
!σの関数
real function minmod(x,y)
	real :: x,y
	if(abs(x)<abs(y)) then
		minmod=x
	else if(abs(x)>abs(y)) then
		minmod=y
	else
		minmod=0
	end if
end
!minmodの関数

コメント