main.f90

利用该程序

module Euler

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

 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=65000   !出口参数

 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

    call cellArea()
	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 !*area(i,j)
     U(i,j,2)=rou_in*u_in !*area(i,j)
     U(i,j,3)=0
     U(i,j,4)=(p_in/(gama-1.0)+0.5*rou_in*u_in**2) !*area(i,j)
    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()
   
   real::rou,vx,vy,p,kx,ky,sx,sy,coeff,pp,pb,vx_p,vy_p,vx_w,vy_w
   real::temp,V,dx,dy,rouw,pw,c
   real::ub,vb,pl,ul,vl,up,vp,roup,roub,roul,u_in,v_in
   real:: c_in
   integer:: i,j

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

   
   call cellArea() 


do i=3,M+2  
!#################上壁面边界条件#################!
 
   rou=U(i,N+2,1)                    !流场最外层紧靠边界的单元格心处参数,同上面的L              
   vx=U(i,N+2,2)/U(i,N+2,1)
   vy=U(i,N+2,3)/U(i,N+2,1)
   p=(gama-1.0)*(U(i,N+2,4)-0.5*rou*(vx**2+vy**2))
   c=sqrt(gama*p/rou)
   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)
   
   pw=p ! +rou*c*(kx/coeff*vx+ky/coeff*vy)  !W 代表无滑移边界单元参数，过渡，实际不存在，参考文献David L.Whitfield
   rouw=rou !+(pw-p)/(c*c)
   vx_w=vx-kx/coeff*(kx/coeff*vx+ky/coeff*vy)
   vy_w=vy-ky/coeff*(kx/coeff*vx+ky/coeff*vy)

   pp=2*pw-p      !P 代表无滑移边界外的虚单元参数，参考文献David L.Whitfield
   roup=2*rouw-rou
   vx_p=2*vx_w-vx 
   vy_p=2*vy_w-vy
   
   U(i,N+3,1)=roup                           
   U(i,N+3,2)=roup*vx_p
   U(i,N+3,3)=roup*vy_p
   U(i,N+3,4)=pp/(gama-1.0)+0.5*roup*(vx_p**2+vy_p**2)   

   U(i,N+4,:)=U(i,N+3,:)                       !直接外推至第二层虚单元

!#################下壁面边界条件#################!

   rou=U(i,3,1)
   vx=U(i,3,2)/U(i,3,1)
   vy=U(i,3,3)/U(i,3,1)
   p=(gama-1.0)*(U(i,3,4)-0.5*rou*(vx**2+vy**2))
   c=sqrt(gama*p/rou)
   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)
   
   pw=p !-rou*c*(kx/coeff*vx+ky/coeff*vy)  !W 代表无滑移边界单元参数，过渡，实际不存在，参考文献David L.Whitfield
   rouw=rou !+(pw-p)/(c*c)
   vx_w=vx-kx/coeff*(kx/coeff*vx+ky/coeff*vy)
   vy_w=vy-ky/coeff*(kx/coeff*vx+ky/coeff*vy)

   pp=2*pw-p      !P 代表无滑移边界外的虚单元参数，参考文献David L.Whitfield
   roup=2*rouw-rou
   vx_p=2*vx_w-vx 
   vy_p=2*vy_w-vy

   U(i,2,1)=roup
   U(i,2,2)=roup*vx_p 
   U(i,2,3)=roup*vy_p
   U(i,2,4)=pp/(gama-1.0)+0.5*roup*(vx_p**2+vy_p**2)!press_p
   U(i,1,:)=U(i,2,:)
 end do
!#################亚音流入口#################!

do j=3,N+2
   ul=U(3,j,2)/U(3,j,1)         ! L 代表流场之内紧靠边界的单元格心处参数，参考文献David L.Whitfield
   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))
   roul=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)

   pb=0.5*(p_in+pl+roul*c*(kx/coeff*(u_in-ul)+ky/coeff*(v_in-vl)))
   roub=rou_in+(pb-p_in)/c/c        ! b 代表流场最外层边界单元格心处参数，用于中间过渡，实际不存在
   ub=u_in+kx/coeff*(p_in-pb)/roul/c
   vb=v_in+ky/coeff*(p_in-pb)/roul/c

   pp=2*pb-pl                     ! p 代表第一层虚单元格心处参数
   roup=2*roub-roul
   up=2*ub-ul
   vp=2*vb-vl
   
   U(2,j,1)=roup              !第一层虚单元
   U(2,j,2)=U(2,j,1)*up
   U(2,j,3)=U(2,j,1)*vp
   U(2,j,4)=pp/(gama-1)+0.5*U(2,j,1)*(up**2+vp**2)

   U(1,j,:)=U(2,j,:)         !第二层虚单元
end do

!#################亚音流出口#################!
do j=3,N+2
   !采用直接外推的方法
   !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,:)
   !另一种方法
   ul=U(M+2,j,2)/U(M+2,j,1)         ! L=a 代表流场之内紧靠边界的单元格心处参数，参考文献David L.Whitfield
   vl=U(M+2,j,3)/U(M+2,j,1)
   pl=(gama-1.0)*(U(M+2,j,4)-0.5*(ul**2+vl**2)*U(M+2,j,1))
   roul=U(M+2,j,1)
   c=sqrt(gama*(gama-1)*(U(M+2,j,4)/U(M+2,j,1)-0.5*(ul**2+vl**2)))

   dx=x(M+2,j+1)-x(M+2,j) !
   dy=y(M+2,j+1)-y(M+2,j) !
   sx=dy
   sy=-dx
   kx=sx/area(M+2,j)
   ky=sy/area(M+2,j)
   coeff=sqrt(kx*kx+ky*ky)
    
   pb=p_out
   roub=roul+(pb-pl)/c/c
   ub=ul+kx/coeff*(pl-pb)/roul/c
   vb=vl+ky/coeff*(pl-pb)/roul/c 

   pp=2*pb-pl                     ! p 代表第一层虚单元格心处参数
   roup=2*roub-roul
   up=2*ub-ul
   vp=2*vb-vl

   U(M+3,j,1)=roup             !第一层虚单元
   U(M+3,j,2)=U(M+3,j,1)*up
   U(M+3,j,3)=U(M+3,j,1)*vp
   U(M+3,j,4)=pp/(gama-1)+0.5*U(M+3,j,1)*(up**2+vp**2)
   U(M+4,j,:)=U(M+3,j,:)         !第二层虚单元
 end do

end subroutine boundary_condition



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

subroutine MUSCL()


implicit none
real::V(M+4,N+4,4),UL(4),UR(4),derta
real::Si(M+4,N+4,4),FL(4),FR(4)
real::T(4,4),T_(4,4),P(4,4)
real::LAMD(4,4)=0
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,c,D,ave_V
real::dx,dy,kx,ky,sx,sy,qk,coeff,coeff2,q(4),for(4),aft(4),mid(4),Uch(4),dUch(4)
real::sx1,sx2,sy1,sy2,kx1,kx2,ky1,ky2,coeff3,coeff4

integer i,j,k

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))


!*******************通量F计算*******************************!

!****************用特征型变量求插值因子****************!

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
      mid=matmul(P,V(i,j,:))
      for=matmul(P,V(i+1,j,:))
   	  aft=matmul(P,V(i-1,j,:))   
	  for=for-mid
	  aft =mid-aft
      q =2*for*aft/(for*for+aft*aft+epsilon)
      Si(i,j,:)=0.5*q*(for+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)) 
      dUch=matmul(T_,Si(i,j,:))	
	  UR=U(i+1,j,:)-0.5*dUch      										   !所以此处需还原为守恒型变量，Uch=T_*U=P_*Up
      UL=U(i,j,:)+0.5*dUch                         !界面左侧变量

	  ror   =UR(1)                                 !界面右侧参数
      vxr   =UR(2)/UR(1)
      vyr   =UR(3)/UR(1)
      ER    =UR(4)/UR(1)
      pressR=ror*(gama-1)*(ER-0.5*(vxr**2+vyr**2))
      HR    =gama*ER-0.5*(gama-1)*(vxr**2+vyr**2)

	  rol   =UL(1)                                 !界面左侧参数
      vxl   =UL(2)/UL(1)
      vyl   =UL(3)/UL(1)
      EL    =UL(4)/UL(1)
      pressL=rol*(gama-1)*(EL-0.5*(vxl**2+vyl**2))
      HL    =gama*EL-0.5*(gama-1)*(vxl**2+vyl**2)
	   
      FR(1)=UR(1)*(sx*vxr+sy*vyr)                   !界面右侧通量
      FR(2)=UR(2)*(sx*vxr+sy*vyr)+sx*pressR
      FR(3)=UR(3)*(sx*vxr+sy*vyr)+sy*pressR
      FR(4)=UR(1)*(sx*vxr+sy*vyr)*HR