main.f90

利用该程序

	  FL(1)=UL(1)*(sx*vxl+sy*vyl)                   !界面左侧通量
      FL(2)=UL(2)*(sx*vxl+sy*vyl)+sx*pressL
      FL(3)=UL(3)*(sx*vxl+sy*vyl)+sy*pressL
      FL(4)=UL(1)*(sx*vxl+sy*vyl)*HL

      D=sqrt(ror/rol)                               !roe平均界面两侧参数
      ave_vx=(vxl+vxr*D)/(1.0+D)
      ave_vy=(vyl+vyr*D)/(1.0+D)
      ave_H =(HL+HR*D)/(1.0+D)
      ave_V =ave_vx**2+ave_vy**2;
      c=sqrt((gama-1)*(ave_H-0.5*(ave_vx**2+ave_vy**2)))
      qk=kx*ave_vx+ky*ave_vy

   T(1,1)=0.5*ave_V-c*c/(gama-1)
   T(1,2)=-ave_vx
   T(1,3)=-ave_vy
   T(1,4)=1.0
   T(2,1)=-kx*ave_vy+ky*ave_vx
   T(2,2)=-ky
   T(2,3)=kx
   T(2,4)=0
   T(3,1)=-qk+0.5*ave_V*(1-gama)*coeff/c
   T(3,2)=kx+(gama-1)*ave_vx*coeff/c
   T(3,3)=ky+(gama-1)*ave_vy*coeff/c
   T(3,4)=-(gama-1)*coeff/c
   T(4,1)=-qk-0.5*ave_V*(1-gama)*coeff/c
   T(4,2)=kx-(gama-1)*ave_vx*coeff/c
   T(4,3)=ky-(gama-1)*ave_vy*coeff/c
   T(4,4)=(gama-1)*coeff/c
        
   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)*ave_vx/(c*c)
   T_(2,2)=-ky/(coeff**2)
   T_(2,3)=0.5*kx/(coeff**2)-0.5*ave_vx/c/coeff
   T_(2,4)=0.5*kx/(coeff**2)+0.5*ave_vx/c/coeff
   T_(3,1)=-(gama-1)*ave_vy/(c*c)
   T_(3,2)=kx/(coeff**2)
   T_(3,3)=0.5*ky/(coeff**2)-0.5*ave_vy/c/coeff
   T_(3,4)=0.5*ky/(coeff**2)+0.5*ave_vy/c/coeff
   T_(4,1)=-0.5*(gama-1)/(c*c)*ave_V
   T_(4,2)=(kx*ave_vy-ky*ave_vx)/(coeff**2)
   T_(4,3)=0.5/c/coeff*(c*qk/coeff-0.5*ave_V-c*c/(gama-1))
   T_(4,4)=0.5/c/coeff*(c*qk/coeff+0.5*ave_V+c*c/(gama-1)) 

   LAMD(1,1)=abs(sx*ave_vx+sy*ave_vy)
   LAMD(2,2)=abs(sx*ave_vx+sy*ave_vy)
   LAMD(3,3)=abs(sx*ave_vx+sy*ave_vy-c*coeff2)
   LAMD(4,4)=abs(sx*ave_vx+sy*ave_vy+c*coeff2)

   derta=0.15*(abs(sx*ave_vx+sy*ave_vy)+c*coeff2) !熵校正
   do k=1,4
     if (LAMD(k,k)<derta) then
	   LAMD(k,k)=(LAMD(k,k)**2+derta**2)/2/derta 
     end if
   end do
   F_epslion(i,j,:)=0.5*(FL+FR)-0.5*matmul(matmul(matmul(T_,LAMD),T),(UR-UL))

  end do
end do

!用特征型变量*******************通量G计算*******************************!
!求插值因子
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,j+1)-x(i+1,j+1) !
      dy=y(i,j+1)-y(i+1,j+1) !
      sx=dy  !
      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,j+1,:))
      aft=matmul(P,V(i,j-1,:))
	  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,j+1)-x(i+1,j+1) !
      dy=y(i,j+1)-y(i+1,j+1) !
      sx=dy  !
      sy=-dx !
      kx=sx/area(i,j)
      ky=sy/area(i,j)
      coeff=sqrt(kx*kx+ky*ky)
      coeff2=sqrt(sx*sx+sy*sy)
      qk=kx*vx+ky*vy
      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)) 
 	                                               !注意，由于S由特征型变量算得，
	  dUch=matmul(T_,Si(i,j,:))										   !所以此处需还原为守恒型变量，Uch=T_*U=P_*Up     
      UR=U(i,j+1,:)-0.5*dUch                       !界面上侧变量
      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


	  FL(1)=UL(1)*(sx*vxl+sy*vyl)                   !界面下侧通量
      FL(2)=UL(2)*(sx*vxl+sy*vyl)+sx*pressL
      FL(3)=UL(3)*(sx*vxl+sy*vyl)+sy*pressL
      FL(4)=UL(1)*(sx*vxl+sy*vyl)*HL

      D=sqrt(ror/rol)                               !roe平均界面两侧参数
      ave_vx=(vxl+vxr*D)/(1.0+D)
      ave_vy=(vyl+vyr*D)/(1.0+D)
      ave_H =(HL+HR*D)/(1.0+D)
      ave_V =ave_vx**2+ave_vy**2;
      c=sqrt((gama-1)*(ave_H-0.5*(ave_vx**2+ave_vy**2)))
      qk=kx*ave_vx+ky*ave_vy

   T(1,1)=0.5*ave_V-c*c/(gama-1)
   T(1,2)=-ave_vx
   T(1,3)=-ave_vy
   T(1,4)=1.0
   T(2,1)=-kx*ave_vy+ky*ave_vx
   T(2,2)=-ky
   T(2,3)=kx
   T(2,4)=0
   T(3,1)=-qk+0.5*ave_V*(1-gama)*coeff/c
   T(3,2)=kx+(gama-1)*ave_vx*coeff/c
   T(3,3)=ky+(gama-1)*ave_vy*coeff/c
   T(3,4)=-(gama-1)*coeff/c
   T(4,1)=-qk-0.5*ave_V*(1-gama)*coeff/c
   T(4,2)=kx-(gama-1)*ave_vx*coeff/c
   T(4,3)=ky-(gama-1)*ave_vy*coeff/c
   T(4,4)=(gama-1)*coeff/c
        
   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)*ave_vx/(c*c)
   T_(2,2)=-ky/(coeff**2)
   T_(2,3)=0.5*kx/(coeff**2)-0.5*ave_vx/c/coeff
   T_(2,4)=0.5*kx/(coeff**2)+0.5*ave_vx/c/coeff
   T_(3,1)=-(gama-1)*ave_vy/(c*c)
   T_(3,2)=kx/(coeff**2)
   T_(3,3)=0.5*ky/(coeff**2)-0.5*ave_vy/c/coeff
   T_(3,4)=0.5*ky/(coeff**2)+0.5*ave_vy/c/coeff
   T_(4,1)=-0.5*(gama-1)/(c*c)*ave_V
   T_(4,2)=(kx*ave_vy-ky*ave_vx)/(coeff**2)
   T_(4,3)=0.5/c/coeff*(c*qk/coeff-0.5*ave_V-c*c/(gama-1))
   T_(4,4)=0.5/c/coeff*(c*qk/coeff+0.5*ave_V+c*c/(gama-1)) 

   LAMD(1,1)=abs(sx*ave_vx+sy*ave_vy)
   LAMD(2,2)=abs(sx*ave_vx+sy*ave_vy)
   LAMD(3,3)=abs(sx*ave_vx+sy*ave_vy-c*coeff2)
   LAMD(4,4)=abs(sx*ave_vx+sy*ave_vy+c*coeff2)

   derta=0.15*(abs(sx*ave_vx+sy*ave_vy)+c*coeff2)          !熵校正
   do k=1,4
     if (LAMD(k,k)<derta) then
	   LAMD(k,k)=(LAMD(k,k)**2+derta**2)/2/derta 
     end if
   end do
   	  
   !lamdmax_y(i,j)=abs(kx*ave_vx+ky*ave_vy)+c*coeff !用于求CFL条件，当地时间步长
   F_eta(i,j,:)=0.5*(FL+FR)-0.5*matmul(matmul(matmul(T_,LAMD),T),(UR-UL))

  end do
end do

!*******************求普半径*****************************************

 do i=3,M+2   !求单元的lamdamax，用于求当地时间步长
   do j=3,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,j+1)-x(i,j) 
      dy=y(i,j+1)-y(i,j) 
      sx1=dy  !
      sy1=-dx !sy=-dx     
	  kx1=sx1/area(i,j)        
      ky1=sy1/area(i,j)
      coeff3=sqrt(kx1*kx1+ky1*ky1)

      dx=x(i+1,j+1)-x(i+1,j) 
      dy=y(i+1,j+1)-y(i+1,j) 
      sx2=dy  !
      sy2=-dx !sy=-dx
      kx2=sx2/area(i,j)                     
      ky2=sy2/area(i,j) 
	  coeff4=sqrt(kx2*kx2+ky2*ky2)
	  lamdmax_x(i,j)=abs(0.5*(vx*(kx1+kx2)+vy*(ky1+ky2)))+c*0.5*(coeff3+coeff4)  !用于求CFL条件，当地时间
 
      dx=x(i,j)-x(i+1,j) 
      dy=y(i,j)-y(i+1,j) 
      sx1=dy  !
      sy1=-dx !sy=-dx
      kx1=sx1/area(i,j)        
      ky1=sy1/area(i,j)
      coeff3=sqrt(kx1*kx1+ky1*ky1)

	  dx=x(i,j+1)-x(i+1,j+1) 
      dy=y(i,j+1)-y(i+1,j+1) 
      sx2=dy  !
      sy2=-dx !sy=-dx
      kx2=sx2/area(i,j)                     
      ky2=sy2/area(i,j) 
	  coeff4=sqrt(kx2*kx2+ky2*ky2)
	  lamdmax_y(i,j)=abs(0.5*(vx*(kx1+kx2)+vy*(ky1+ky2)))+c*0.5*(coeff3+coeff4)  !用于求CFL条件，当地时间	 
   end do
 end do
end subroutine MUSCL

    
!**************计算进出口质量流量*******************!

subroutine massflow()

implicit none
real::vx,vy,ro,cita,beta,Vel,dx,dy,signal_u,signal_beta
real::mass, vn
integer j
mass=0

!******进口质量******
do j=3,N+2
   ro=U(3,j,1)
   vx=U(3,j,2)/U(3,j,1)
   vy=U(3,j,3)/U(3,j,1)
   Vel=sqrt(vx*vx+vy*vy)
   cita=atan(vy/vx)
   
   dx=x(3,j+1)-x(3,j)
   dy=y(3,j+1)-y(3,j)
   beta=atan(dy/dx)
   signal_u=abs(vx)/vx
   signal_beta=abs(beta)/beta
   vn=signal_u*signal_beta*Vel*sin(beta-cita)
   mass=mass+ro*vn*sqrt(dx*dx+dy*dy)
end do
write(*,*) 'mass_in=',mass

mass=0
!******出口质量******
do j=3,N+2
   ro=U(M+2,j,1)
   vx=U(M+2,j,2)/U(M+2,j,1)
   vy=U(M+2,j,3)/U(M+2,j,1)
   Vel=sqrt(vx*vx+vy*vy)
   cita=atan(vy/vx)
   
   dx=x(M+2,j+1)-x(M+2,j)
   dy=y(M+2,j+1)-y(M+2,j)
   beta=atan(dy/dx)
   signal_u=abs(vx)/vx
   signal_beta=abs(beta)/beta
   vn=signal_u*signal_beta*Vel*sin(beta-cita)
   mass=mass+ro*vn*sqrt(dx*dx+dy*dy)
end do
write(*,*) 'mass_out=',mass

end subroutine massflow


!!!!!!!!!!!!!!!!!!!!!!!!!!迭代开始!!!!!!!!!!!!!!!!!!!!!!!!!!
  
 
subroutine upwind()

  integer::time,foot,i,j
  real::arefa(5),Residual(4)
  real::ro,max_resden_ro,max_resden_u,max_resden_v,max_resden_E
  real:: vx,vy,p,Ma,rou,c
  real:: Uold(M+4,N+4,4)
  real::dt

  arefa(1)=0.107
  arefa(2)=0.198
  arefa(3)=0.323
  arefa(4)=0.52
  arefa(5)=1.0

  max_resden_ro=0
  max_resden_u=0
  max_resden_v=0
  max_resden_E=0

  open(fileid1,file='result.dat')        !打开文件
  write(fileid1,*) 'VARIABLES="X" "Y" "M" "p" "ro" "vx" "vy"'
  write(fileid1,*) "ZONE F=POINT I=",M,"J=",N !便于tecplot识别，第一行声明,且注意I与J赋值相反

  open(fileid2,file='Resi_rou.dat ')
  open(fileid3,file='Resi_u.dat')
  open(fileid4,file='Resi_v.dat')
  open(fileid5,file='Resi_E.dat')
  write(fileid2,*) "time    Resi_rou     "
    write(fileid3,*) "time      Resi_u    "
     write(fileid4,*)"time       Resi_v     "
      write(fileid5,*)"time       Resi_E"
    


  !call mesh()
  call initialize()          !物理区域参数初始化
  call cellArea()
  call boundary_condition()

  do time=1,step1  !时间推进循环
    Uold=U

    do foot=1,step2  !龙格库塔多步循环
      call MUSCL() 
	  call boundary_condition()  
       do i=3,M+2
	     do j=3,N+2
           dt=CFL/(lamdmax_x(i,j)+lamdmax_y(i,j))
		   U(i,j,:)=U(i,j,:)-dt*arefa(foot)/area(i,j)*(F_epslion(i,j,:)-F_epslion(i-1,j,:)+F_eta(i,j,:)-F_eta(i,j-1,:))
		 end do
	   end do
	   
	end do          !龙格库塔多步循环结束
	         
!残差计算
   	do i=3,M+2
	    do j=3,N+2
		  ro=U(i,j,1)
		  Residual(1)=abs(Uold(i,j,1)-U(i,j,1)) !(F_epslion(i,j,1)-F_epslion(i-1,j,1)+F_eta(i,j,1)-F_eta(i,j-1,1))/area(i,j)
		  Residual(2)=abs(Uold(i,j,2)-U(i,j,2))/Uold(i,j,1) !(F_epslion(i,j,2)-F_epslion(i-1,j,2)+F_eta(i,j,2)-F_eta(i,j-1,2))/ro/area(i,j)
		  Residual(3)=abs(Uold(i,j,3)-U(i,j,3))/Uold(i,j,1) !(F_epslion(i,j,3)-F_epslion(i-1,j,3)+F_eta(i,j,1)-F_eta(i,j-1,3))/ro/area(i,j)
		  Residual(4)=abs(Uold(i,j,4)-U(i,j,4))/Uold(i,j,1)   !(F_epslion(i,j,4)-F_epslion(i-1,j,4)+F_eta(i,j,4)-F_eta(i,j-1,4))/ro/area(i,j)
          if(Residual(1)>=max_resden_ro)then
		     max_resden_ro=Residual(1)
		  end if
		 
          if(Residual(2)>=max_resden_u)then
		    max_resden_u=Residual(2)
		  end if

          if(Residual(3)>=max_resden_v)then
		    max_resden_v=Residual(3)
		  end if

		  if(Residual(4)>=max_resden_E)then
		     max_resden_E=Residual(4)
		  end if

        end do
	  end do
	write(*,*) "step=",time,"  Resi_rou=", max_resden_ro
    write(fileid2,*) time ,     max_resden_ro
    write(fileid3,*) time ,     max_resden_u
    write(fileid4,*) time ,      max_resden_v
    write(fileid5,*) time ,      max_resden_E

      max_resden_ro=0
	  max_resden_u=0
	  max_resden_v=0
      max_resden_E=0
    call massflow()
  end do   !时间推进结束

  !往文件里写数据
  do j=3,N+2
    do i=3,M+2
	  rou=U(i,j,1)
	  vx=U(i,j,2)/U(i,j,1)
	  vy=U(i,j,3)/U(i,j,1)
      c =sqrt(gama*(gama-1)*(U(i,j,4)/U(i,j,1)-0.5*(vx*vx+vy*vy)))
	  p =(gama-1)*(U(i,j,4)-0.5*rou*(vx**2+vy**2))
	  Ma =sqrt(vx**2+vy**2)/c
	  write(fileid1,*) x(i+1,j+1),y(i+1,j+1),Ma,p,rou,vx,vy 
	end do
  end do
 !写完数据
end subroutine upwind		   


end module





program Euler_2d  !主函数
  
  use Euler
  !call mesh()
  !call initialize()
  !call boundary_condition()
  call upwind()

 stop
end


!程序未完，网格点上参数虚插值得到