modi.f90

用FORTRAN语言实现二维欧拉方程的求解

module Euler

!!!!!!!!!!!!!!!!!!!!!module中的变量及常量!!!!!!!!!!!!!!!!!!!!!!!!
 implicit none
 !private x,y,U,F,G,JJ,R
 integer,parameter::M=100  !x方向等分数
 integer,parameter::N=40  !y方向等分数
 integer,parameter::fileid1=10 !文件1代码＝10 
 real,parameter::CFL=0.5
 real,parameter::gama=1.4
 real,parameter::R=287.06
 real,parameter::epsilon=0.0000001
 integer,parameter::step1=3000 !时间步
 integer,parameter::step2=5    !龙格库塔推进步
 integer,parameter::dim=4

 real,parameter::p_in=101325   !来流参数
 real,parameter::rou_in=1.293
 real,parameter::t_in=273.15
 real,parameter::M_in=0.35
 real,parameter::p_out=70000   !出口参数

 real::x(M+5,N+5)     
 real::y(M+5,N+5)
 real::U(M+4,N+4,4)
 real::area(M+4,N+4)
 real::F_epslion(M+4,N+4,4),F_eta(M+4,N+4,4)  !界面处通量
 real::lamdmax_x(M+4,N+4)                 !用于求当地时间步长
 real::lamdmax_y(M+4,N+4)

 common U,x,y,area,F_epslion,F_eta,lamdmax_x,lamdmax_y



!!!!!!!!!!!!!!!!!!module中的函数!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
contains



!@@@@@@@@@@@@@@@@@@@@@@@@@@@网格生成子程序@@@@@@@@@@@@@@@@@@@@@@@@@@@!
subroutine mesh()
   integer i,j 
   real::fi_x1,fi_x2,fi_y1,fi_y2
   real::kecy(M+5)=0.0
   real::yita(N+5)=0.0


   !利用剪切变换kecy=(x-xl)/(xr-xl)和yita=(y-yd)/(yu-yd)
   !可将物理平面映射到计算平面，且kecy和yita从0到1均分
   do i=3,M+3
     kecy(i)=real(i-3)/real(M)
   end do

   do j=3,N+3
     yita(j)=real(j-3)/real(N)
   end do

    do i=3,M+3
   !下边界x坐标
     x(i,3)=3.1*real(i-3)/real(M)  !x从0到3.1
   !下边界y坐标
     y(i,3)=-0.4-0.2*cos(2*x(i,3))
   !上边界x坐标
     x(i,N+3)=x(i,3)                !x从0到3.1
   !上边界y坐标
     y(i,N+3)=0.4+0.2*cos(2*x(i,N+3))
   end do

    do j=3,N+3
   !左边界x坐标
     x(3,j)=0
   !左边界y坐标
     y(3,j)=(-0.4-0.2*cos(2*x(3,3)))+(0.8+0.4*cos(2*x(3,N+3)))*real(j-3)/real(N)
   !右边界x坐标
     x(M+3,j)=3.1
   !右边界y坐标
     y(M+3,j)=y(3,j)
   end do

   do i=4,M+2
      fi_x1=1-kecy(i)
	  fi_x2=kecy(i)
	  do j=4,N+2
	    fi_y1=1-yita(j)
		fi_y2=yita(j)		
        !用超限插值法计算内部坐标值
		x(i,j)=fi_y1*x(i,3)+fi_y2*x(i,N+3)+fi_x1*(x(3,j)-fi_y1*x(3,3)-fi_y2*x(3,N+3))+fi_x2*(x(M+3,j)-fi_y1*x(M+3,3)-fi_y2*x(M+3,N+3))
        y(i,j)=fi_y1*y(i,3)+fi_y2*y(i,N+3)+fi_x1*(y(3,j)-fi_y1*y(3,3)-fi_y2*y(3,N+3))+fi_x2*(y(M+3,j)-fi_y1*y(M+3,3)-fi_y2*y(M+3,N+3))
      end do
    end do

end subroutine mesh





!@@@@@@@@@@@@@@@@@@@@@@@@@@@物理区域流场初始化子程序@@@@@@@@@@@@@@@@@@@@@@@@@@@!
subroutine initialize()
     
	 
!以进口参数初始化
 real:: c_in,u_in
 integer:: i,j

	c_in=sqrt(gama*R*t_in)
	u_in=M_in*c_in

  
!各个单元格心
  do i=1,M+4
    do j=1,N+4
     U(i,j,1)=rou_in
     U(i,j,2)=rou_in*u_in
     U(i,j,3)=0
     U(i,j,4)=p_in/(gama-1.0)+0.5*rou_in*u_in**2
    end do
  end do

 end subroutine  initialize


!@@@@@@@@@@@@@@@@@@@@@@@@@@@单元面积计算子程序@@@@@@@@@@@@@@@@@@@@@@@@@@@@@!
subroutine cellArea()
   
   implicit none
   real:: dxAC,dyAC,dxBD,dyBD
   integer i,j
   
call mesh()

do j=3,N+3
   x(2,j)=x(3,j)+(x(3,j)-x(4,j))
   y(2,j)=y(3,j)+(y(3,j)-y(4,j))
   x(1,j)=x(2,j)+(x(2,j)-x(3,j))
   y(1,j)=y(2,j)+(y(2,j)-y(3,j))

   x(M+4,j)=x(M+3,j)+(x(M+3,j)-x(M+2,j))
   y(M+4,j)=y(M+3,j)+(y(M+3,j)-y(M+2,j))
   x(M+5,j)=x(M+4,j)+(x(M+4,j)-x(M+3,j))
   y(M+5,j)=y(M+4,j)+(y(M+4,j)-y(M+3,j))
end do

do i=1,M+5
   x(i,2)=x(i,3)+(x(i,3)-x(i,4))
   y(i,2)=y(i,3)+(y(i,3)-y(i,4))
   x(i,1)=x(i,2)+(x(i,2)-x(i,3))
   y(i,1)=y(i,2)+(y(i,2)-y(i,3))

   x(i,N+4)=x(i,N+3)+(x(i,N+3)-x(i,N+2))
   y(i,N+4)=y(i,N+3)+(y(i,N+3)-y(i,N+2))
   x(i,N+5)=x(i,N+4)+(x(i,N+4)-x(i,N+3))
   y(i,N+5)=y(i,N+4)+(y(i,N+4)-y(i,N+3))
end do

do i=1,M+4
   do j=1,N+4
	dxAC=(x(i,j+1)-x(i+1,j))
	dyAC=(y(i,j+1)-y(i+1,j))
	dxBD=(x(i,j)-x(i+1,j+1))
	dyBD=(y(i,j)-y(i+1,j+1)) 

	area(i,j)=0.5*(dxAC*dyBD-dxBD*dyAC)
   end do
end do
end subroutine cellArea




!@@@@@@@@@@@@@@@@@@@@@@@@@@@边界条件处理子程序@@@@@@@@@@@@@@@@@@@@@@@@@@@@@!
subroutine boundary_condition()
   
  implicit none 

real:: ro,vx,vy,press,kx,ky,sx,sy,coeff,ro_p,press_p,vx_p,vy_p,vx_w,vy_w
real:: temp,V,dx,dy,ro_w,press_w,c
real:: ub,vb,p_b,pl,ul,vl,p_p,up,vp,rop,rob,rol,u_in,v_in
integer i,j 


	temp=(gama-1)/2/gama
!	M=sqrt(2*(pow(pt/p_in,(gama-1)/gama)-1)/(gama-1))
	c=sqrt(gama*R*T_in)
	V=M_in*c
	u_in=V
	v_in=0

    call cellArea()
!-----------------进口、出口边界条件---------------------------
do j=3,N+2
   ul=U(3,j,2)/U(3,j,1)
   vl=U(3,j,3)/U(3,j,1)
   pl=(gama-1.0)*(U(3,j,4)-0.5*(ul**2+vl**2)*U(3,j,1))
   rol=U(3,j,1)
   c=sqrt(gama*(gama-1)*(U(3,j,4)/U(3,j,1)-0.5*(ul**2+vl**2)))

   dx=x(3,j+1)-x(3,j) !
   dy=y(3,j+1)-y(3,j) !
   sx=dy
   sy=-dx
   kx=sx/area(3,j)
   ky=sy/area(3,j)
   coeff=sqrt(kx*kx+ky*ky)
   p_b=0.5*(p_in+pl+rol*c*(kx/coeff*(u_in-ul)+ky/coeff*(v_in-vl)))
   rob=rou_in+(p_b-p_in)/c/c
   ub=u_in+kx/coeff*(p_in-p_b)/rol/c
   vb=v_in+ky/coeff*(p_in-p_b)/rol/c

   p_p=2*p_b-pl
   rop=2*rob-rol
   up =2*ub-ul
   vp =2*vb-vl
   U(2,j,1)=2*rob-rol
   U(2,j,2)=U(2,j,1)*(2*ub-ul)
   U(2,j,3)=U(2,j,1)*(2*vb-vl)
   U(2,j,4)=p_p/(gama-1)+0.5*U(2,j,1)*(up**2+vp**2)
   U(1,j,:)=U(2,j,:)

   U(M+3,j,:)=U(M+2,j,:)
   U(M+3,j,4)=p_out/(gama-1.0)+0.5*(U(M+3,j,2)**2+U(M+3,j,3)**2)/U(M+3,j,1)
   U(M+4,j,:)=U(M+3,j,:)
!   write(*,*) U(NX+3,j,2)
end do
!----------------------------------------------------------------------
!-----------------------壁面边界条件-----------------------------------
do i=3,M+2
   !--------------downwall---------------------------------------------
   ro=U(i,3,1)
   vx=U(i,3,2)/U(i,3,1)
   vy=U(i,3,3)/U(i,3,1)
   press=(gama-1.0)*(U(i,3,4)-0.5*ro*(vx**2+vy**2))
   c=sqrt(gama*press/ro)
   dx=x(i+1,3)-x(i,3) !
   dy=y(i+1,3)-y(i,3) !
   sx=dy
   sy=-dx
   kx=sx/area(i,3)
   ky=sy/area(i,3)
   coeff=sqrt(kx*kx+ky*ky)
   
   press_w=press-ro*c*(kx/coeff*vx+ky/coeff*vy)
   ro_w=ro+(press_w-press)/(c*c)
   vx_w=vx-kx/coeff*(kx/coeff*vx+ky/coeff*vy)
   vy_w=vy-ky/coeff*(kx/coeff*vx+ky/coeff*vy)

   press_p=press-2*ro*c*(kx/coeff*vx+ky/coeff*vy)
   ro_p=ro+2*(press_w-press)/(c*c)
   vx_p=vx-2*kx/coeff*(kx/coeff*vx+ky/coeff*vy) !vx_w
   vy_p=vy-2*ky/coeff*(kx/coeff*vx+ky/coeff*vy) !vy_w 

   U(i,2,1)=ro
   U(i,2,2)=ro*vx_p !ro_p
   U(i,2,3)=ro*vy_p
   U(i,2,4)=press/(gama-1.0)+0.5*ro*(vx_p**2+vy_p**2)!press_p
   U(i,1,:)=U(i,2,:)
!   write(*,*) 'p=',press
   !---------------------------------------------------------------------
   !------------------------upwall---------------------------------------
   ro=U(i,N+2,1)
   vx=U(i,N+2,2)/U(i,N+2,1)
   vy=U(i,N+2,3)/U(i,N+2,1)
   press=(gama-1.0)*(U(i,N+2,4)-0.5*ro*(vx**2+vy**2))
   c=sqrt(gama*press/ro)
   dx=x(i,N+3)-x(i+1,N+3) !
   dy=y(i,N+3)-y(i+1,N+3) !
   sx=dy  !sx=len*cos(Ncita)
   sy=-dx !sy=len*sin(Ncita)
   kx=sx/area(i,N+2)
   ky=sy/area(i,N+2)
   coeff=sqrt(kx*kx+ky*ky)
   
   press_w=press+ro*c*(kx/coeff*vx+ky/coeff*vy)
   ro_w=ro+(press_w-press)/(c*c)
   vx_w=vx-kx/coeff*(kx/coeff*vx+ky/coeff*vy)
   vy_w=vy-ky/coeff*(kx/coeff*vx+ky/coeff*vy)

   press_p=press+2*ro*c*(kx/coeff*vx+ky/coeff*vy)
   ro_p=ro+2*(press_w-press)/(c*c)
   vx_p=vx-2*kx/coeff*(kx/coeff*vx+ky/coeff*vy) !vx_w
   vy_p=vy-2*ky/coeff*(kx/coeff*vx+ky/coeff*vy) !vy_w

   U(i,N+3,1)=ro
   U(i,N+3,2)=ro*vx_p
   U(i,N+3,3)=ro*vy_p
   U(i,N+3,4)=press/(gama-1.0)+0.5*ro*(vx_p**2+vy_p**2)
   U(i,N+4,:)=U(i,N+3,:)
!   write(*,*) 'vb=',kx*vx_w+ky*vy_w
!   write(*,*) 'ro2=',ro_p
!------------------------------------------------------------------------
end do
end subroutine boundary_condition





!@@@@@@@@@@@@@@@@@@@@@@@@@@@MUSCL插值及迎风计算!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@!

subroutine MUSCL()


implicit none
real:: V(M+4,N+4,dim),UL(dim),UR(dim)
real:: FL(dim),FR(dim)
real:: T(dim,dim),T_(dim,dim),P(dim,dim),LAMD(dim,dim),Si(M+4,N+4,4)
real:: rol,ror,vxl,vyl,vxr,vyr,ER,EL,HL,HR,pressL,pressR,ro,vx,vy,E,H,Vel
real:: ave_vx,ave_vy,ave_H,ave_E,c,D,ave_V
real::dx,dy,kx,ky,sx,sy,qk,coeff,coeff2,q(dim),fore(dim),aft(dim),mid(dim),Uch(dim),dUch(dim)

integer i,j

call cellArea()

!*******守恒型变量转化为原始变量***************!
   
   V(:,:,1)=U(:,:,1)
   V(:,:,2)=U(:,:,2)/U(:,:,1)
   V(:,:,3)=U(:,:,3)/U(:,:,1)
   V(:,:,4)=(gama-1)*(U(:,:,4)-0.5*V(:,:,1)*(V(:,:,2)**2+V(:,:,3)**2))

!**********************************************!

do i=2,M+3
   do j=2,N+3
      ro=U(i,j,1)
	  vx=U(i,j,2)/U(i,j,1)
	  vy=U(i,j,3)/U(i,j,1)
      E =U(i,j,4)/U(i,j,1)
      H    =gama*E-0.5*(gama-1)*(vx**2+vy**2)
      c=sqrt((gama-1)*(H-0.5*(vx**2+vy**2)))
      dx=x(i+1,j+1)-x(i+1,j) 
      dy=y(i+1,j+1)-y(i+1,j) 
      sx=dy  !
      sy=-dx !sy=-dx
      kx=sx/area(i,j)
      ky=sy/area(i,j)
      coeff=sqrt(kx*kx+ky*ky)

      P(1,1)=-(c*c)/(gama-1)
      P(1,2)=0
      P(1,3)=0
      P(1,4)=1/(gama-1)
      P(2,1)=0
      P(2,2)=-ky*ro
      P(2,3)=kx*ro
      P(2,4)=0
      P(3,1)=0
      P(3,2)=kx*ro
      P(3,3)=ky*ro
      P(3,4)=-coeff/c
      P(4,1)=0
      P(4,2)=kx*ro
      P(4,3)=ky*ro
      P(4,4)=coeff/c

      call Matrix_vector(P,V(i+1,j,:),fore)
      call Matrix_vector(P,V(i,j,:),mid)
      call Matrix_vector(P,V(i-1,j,:),aft)

	  fore=fore-mid
	  aft =mid-aft
      q   =2*fore*aft/(fore*fore+aft*aft+epsilon)
      Si(i,j,:)=0.5*q*(fore+aft)   !用特征型变量求插值因子


   end do
end do
do i=2,M+2
   do j=2,N+2
      ro=U(i,j,1)
	  vx=U(i,j,2)/U(i,j,1)
	  vy=U(i,j,3)/U(i,j,1)
      E =U(i,j,4)/U(i,j,1)
      H    =gama*E-0.5*(gama-1)*(vx**2+vy**2)
      c=sqrt((gama-1)*(H-0.5*(vx**2+vy**2)))
      dx=x(i+1,j+1)-x(i+1,j) 
      dy=y(i+1,j+1)-y(i+1,j) 
      sx=dy  !
      sy=-dx !sy=-dx
      kx=sx/area(i,j)
      ky=sy/area(i,j)
      coeff=sqrt(kx*kx+ky*ky)
      coeff2=sqrt(sx*sx+sy*sy)
	  Vel=vx*vx+vy*vy

      T_(1,1)=-(gama-1)/(c*c)
      T_(1,2)=0
      T_(1,3)=-0.5/c/coeff
      T_(1,4)=0.5/c/coeff
      T_(2,1)=-(gama-1)*vx/(c*c)
      T_(2,2)=-ky/(coeff**2)
      T_(2,3)=0.5*kx/(coeff**2)-0.5*vx/c/coeff
      T_(2,4)=0.5*kx/(coeff**2)+0.5*vx/c/coeff
      T_(3,1)=-(gama-1)*vy/(c*c)
      T_(3,2)=kx/(coeff**2)
      T_(3,3)=0.5*ky/(coeff**2)-0.5*vy/c/coeff
      T_(3,4)=0.5*ky/(coeff**2)+0.5*vy/c/coeff
      T_(4,1)=-0.5*(gama-1)/(c*c)*Vel
      T_(4,2)=(kx*vy-ky*vx)/(coeff**2)
      T_(4,3)=0.5/c/coeff*(c*qk/coeff-0.5*Vel-c*c/(gama-1))
      T_(4,4)=0.5/c/coeff*(c*qk/coeff+0.5*Vel+c*c/(gama-1)) 

      call Matrix_vector(T_,Si(i+1,j,:),dUch)
      UR=U(i+1,j,:)-0.5*dUch
      call Matrix_vector(T_,Si(i,j,:),dUch)