modi.f90

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

      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)
   D=sqrt(ror/rol)
   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_E =(EL+ER*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

   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

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

   lamdmax_x(i,j)=abs(kx*ave_vx+ky*ave_vy)+c*coeff

   call Matrix(T_,LAMD,T,Uch,UR-UL)
   F_epslion(i,j,:)=0.5*(FL+FR)-0.5*Uch
   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

      call Matrix_vector(P,V(i,j+1,:),fore)
      call Matrix_vector(P,V(i,j,:),mid)
      call Matrix_vector(P,V(i,j-1,:),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,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)) 

      call Matrix_vector(T_,Si(i,j+1,:),dUch)
      UR=U(i,j+1,:)-0.5*dUch
      call Matrix_vector(T_,Si(i,j,:),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)
   D=sqrt(ror/rol)
   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))) !sqrt((gama-1)*(ave_H-0.5*(ave_vx**2+ave_vy**2))) 
 
   qk=kx*ave_vx+ky*ave_vy
   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

   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

   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)
   lamdmax_y(i,j)=abs(kx*ave_vx+ky*ave_vy)+c*coeff

   call Matrix(T_,LAMD,T,Uch,UR-UL)

   F_eta(i,j,:)=0.5*(FL+FR)-0.5*Uch
!   write(*,*) 'flux=',sx*ave_vx+sy*ave_vy-c*coeff2
   end do
end do
end subroutine MUSCL


subroutine Matrix(T_,LAMD,T,Uch,dU)
use IMSL

implicit none 
real:: T_(dim,dim),LAMD(dim,dim),T(dim,dim),Kc(dim,dim)
real:: Uch(dim),dU(dim)
 
 Kc=(T_.x.LAMD).x.T
 Uch=Kc.x.dU
! write(*,*) 'Kc=',Kc
end subroutine Matrix


subroutine Matrix_vector(P,U,Uch)
use IMSL

implicit none 
real:: P(dim,dim),U(dim),Uch(dim)
 
 Uch=P.x.U
end subroutine Matrix_vector
    
!**************计算进出口质量流量*******************!

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),oldu(M+4,N+4,4)
  real::ro,max_resden_ro,max_resden_u,max_resden_v,max_resden_E
  double precision vx,vy,p,Ma,rou,c
  real::dt
  real             max_old,max_new,ro_old,ro_new,mass_in,mass_out,max_resden
 call Mesh()
call initialize()          !物理区域参数初始化
call cellArea()
call boundary_condition()

call cellArea()
do time=1,step1
   oldU=U
   write(*,*) 'step=',time
   call MUSCL()


!   write(*,*) 'fx=',F_eta(3,3,1)
    do i=3,M+2
	   do j=3,N+2
           dt=1.0*cfl/(lamdmax_x(i,j)+lamdmax_y(i,j))
		   U(i,j,:)=U(i,j,:)-dt/area(i,j)*(F_epslion(i,j,:)-F_epslion(i-1,j,:)+F_eta(i,j,:)-F_eta(i,j-1,:))
		   !write(*,*) 'dt=',dt
	   enddo
	enddo

	 call boundary_condition()
     do i=3,M+2
	    do j=3,N+2
		  !ro_old=oldU(i,j,1)
		  !if(ro_old>max_old)then
		   ! max_old=ro_old
		  !end if
		  !ro_new=U(i,j,1)
		  !if(ro_new>max_new)then
		   ! max_new=ro_new
		  !end if
		  ro=U(i,j,1)
		  Residual(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)=(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)=(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)=(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)then
		  max_resden=Residual(2)
		  end if
		enddo
	 enddo
!	 write(41,*) step,abs(max_old-max_new)/max_old
	 write(*,*) time,max_resden
     call massflow()
  
!	 write(*,*) 'in=',mass_in
	! write(40,*) step,mass_in !abs(mass_in-mass_out)/mass_in
	 
enddo
  !往文件里写数据
  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


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