出口边界.txt

来自「用FORTRAN语言实现二维欧拉方程的求解」· 文本 代码 · 共 35 行

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+2,j,1)=roub              !边界单元变量赋值
   U(M+2,j,2)=U(M+2,j,1)*ub
   U(M+2,j,3)=U(M+2,j,1)*vb
   U(M+2,j,4)=pb/(gama-1)+0.5*U(M+2,j,1)*(ub**2+vb**2) 
   

   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,:)         !第二层虚单元