Numérotation des lignes :

C GMAT2D    SOURCE    CHAT      05/01/13    00:18:34     5004       SUBROUTINE GMAT2D(GMT,WVECT,GA,V_INF,DIAM)CC      Voir kon171.esoCc-----------------------------------------------------c  EQUATIONS:  2D Euler Equations of Gas Dynamicscc  REFERENCE:  AIAA Journal, Sept. 1998c              "Low-Diffusion Flux-Splitting Methodsc               for Flows at All Speeds"c-----------------------------------------------------c  This subroutine provides the preconditioning matrixc  which is the product of the two matrices:c     a)  the 'original' preconditioning matrixc          given in the Reference andc     b)  the jacobian matrix which iscc                  dW/dUc        with W=(p,u,v,T) (T-temperature) andc             U=(rho,rho*u,rho*v,rho*e).c-----------------------------------------------------c  INPUT:cc        wvect   --  vector of the primitive variablesc                    (rho, u, v, p);cc        ga      --  ratio of the specific heats;cc        v_inf   --  reference velocitycc        diam    --  diameter of the cellcc  OUTPUT:cc        gmt     -- preconditioning matrix 4 by 4.c-----------------------------------------------------      IMPLICIT INTEGER(I-N)       real*8 wvect(4),gmt(4,4)       real*8 ga,gm1,a       real*8 rho,p,u,v, v_inf       real*8 rho_p,ur,theta       real*8 epsil,h,qq,m,t1,t2       real*8 diam,lambda,mref,mref2,mach,mach2,dtaum1       parameter(epsil=1.0d0)c--------------------------------------------------c  Inputc--------------------------------------------------       rho=wvect(1)       u=wvect(2)       v=wvect(3)       p=wvect(4)c--------------------------------------------------       gm1=ga-1.0d0       a=sqrt(ga*p/rho)c--------------------------------------------------       rho_p=ga/(a*a)       h=0.5d0*(u*u+v*v)+ga*p/(gm1*rho)c--------------------------------------------------c  Reference velocityc--------------------------------------------------       qq=sqrt(u*u+v*v)       if(qq .lt. (epsil*v_inf)) then          ur = epsil*v_inf       else          ur = qq       endif       if(ur .gt. a)then          ur=a       endifc----------------------------------------------------       theta=1.0d0/(ur*ur) - (1.0d0-ga)*rho_p/gac----------------------------------------------------       m=(1.0d0-ga)*(theta-rho_p)       t1=1.0d0+m*(h-qq*qq)+(theta-rho_p)*ga/rho_p       t2=h*(theta-rho_p)*((1.0d0-ga)*(h-qq*qq)+ga/rho_p)c----- We compute u+a------------       mref=ur/a       mref2=mref*mref       mach=qq/a       mach2=mach*mach       lambda=(1.0D0-mref2)**2       lambda=(lambda*mach2)+(4.0D0*mref2)       lambda=sqrt(lambda)*a       lambda=((1.0D0+mref2)*qq)+lambda       lambda=lambda*0.5D0       dtaum1=lambda/diamc--------------------------------       gmt(1,1)=t1 * dtaum1       gmt(1,2)=u*m * dtaum1       gmt(1,3)=v*m * dtaum1       gmt(1,4)=-1.0D0 * m * dtaum1c-------------------------------       gmt(2,1)=u*(t1-1.0d0) * dtaum1       gmt(2,2)=(1.0d0+u*u*m) * dtaum1       gmt(2,3)=u*v*m * dtaum1       gmt(2,4)=-1.0D0 * m*u * dtaum1c-------------------------------       gmt(3,1)=v*(t1-1.0d0) * dtaum1       gmt(3,2)=u*v*m * dtaum1       gmt(3,3)=(1.0d0+v*v*m) * dtaum1       gmt(3,4)=-1.0D0 * m*v * dtaum1c-------------------------------       gmt(4,1)=t2 * dtaum1       gmt(4,2)=u*h*m * dtaum1       gmt(4,3)=v*h*m * dtaum1       gmt(4,4)=(-m*h+1.0d0) * dtaum1c-------------------------------       return       end

© Cast3M 2003 - Tous droits réservés.
Mentions légales