Télécharger gmat2d.eso

Retour à la liste

Numérotation des lignes :

  1. C GMAT2D SOURCE CHAT 05/01/13 00:18:34 5004
  2. SUBROUTINE GMAT2D(GMT,WVECT,GA,V_INF,DIAM)
  3. C
  4. C Voir kon171.eso
  5. C
  6. c-----------------------------------------------------
  7. c EQUATIONS: 2D Euler Equations of Gas Dynamics
  8. c
  9. c REFERENCE: AIAA Journal, Sept. 1998
  10. c "Low-Diffusion Flux-Splitting Methods
  11. c for Flows at All Speeds"
  12. c-----------------------------------------------------
  13. c This subroutine provides the preconditioning matrix
  14. c which is the product of the two matrices:
  15. c a) the 'original' preconditioning matrix
  16. c given in the Reference and
  17. c b) the jacobian matrix which is
  18. c
  19. c dW/dU
  20. c with W=(p,u,v,T) (T-temperature) and
  21. c U=(rho,rho*u,rho*v,rho*e).
  22. c-----------------------------------------------------
  23. c INPUT:
  24. c
  25. c wvect -- vector of the primitive variables
  26. c (rho, u, v, p);
  27. c
  28. c ga -- ratio of the specific heats;
  29. c
  30. c v_inf -- reference velocity
  31. c
  32. c diam -- diameter of the cell
  33. c
  34. c OUTPUT:
  35. c
  36. c gmt -- preconditioning matrix 4 by 4.
  37. c-----------------------------------------------------
  38. IMPLICIT INTEGER(I-N)
  39. real*8 wvect(4),gmt(4,4)
  40. real*8 ga,gm1,a
  41. real*8 rho,p,u,v, v_inf
  42. real*8 rho_p,ur,theta
  43. real*8 epsil,h,qq,m,t1,t2
  44. real*8 diam,lambda,mref,mref2,mach,mach2,dtaum1
  45. parameter(epsil=1.0d0)
  46. c--------------------------------------------------
  47. c Input
  48. c--------------------------------------------------
  49. rho=wvect(1)
  50. u=wvect(2)
  51. v=wvect(3)
  52. p=wvect(4)
  53. c--------------------------------------------------
  54. gm1=ga-1.0d0
  55. a=sqrt(ga*p/rho)
  56. c--------------------------------------------------
  57. rho_p=ga/(a*a)
  58. h=0.5d0*(u*u+v*v)+ga*p/(gm1*rho)
  59. c--------------------------------------------------
  60. c Reference velocity
  61. c--------------------------------------------------
  62. qq=sqrt(u*u+v*v)
  63. if(qq .lt. (epsil*v_inf)) then
  64. ur = epsil*v_inf
  65. else
  66. ur = qq
  67. endif
  68. if(ur .gt. a)then
  69. ur=a
  70. endif
  71. c----------------------------------------------------
  72. theta=1.0d0/(ur*ur) - (1.0d0-ga)*rho_p/ga
  73. c----------------------------------------------------
  74. m=(1.0d0-ga)*(theta-rho_p)
  75. t1=1.0d0+m*(h-qq*qq)+(theta-rho_p)*ga/rho_p
  76. t2=h*(theta-rho_p)*((1.0d0-ga)*(h-qq*qq)+ga/rho_p)
  77. c----- We compute u+a------------
  78. mref=ur/a
  79. mref2=mref*mref
  80. mach=qq/a
  81. mach2=mach*mach
  82. lambda=(1.0D0-mref2)**2
  83. lambda=(lambda*mach2)+(4.0D0*mref2)
  84. lambda=sqrt(lambda)*a
  85. lambda=((1.0D0+mref2)*qq)+lambda
  86. lambda=lambda*0.5D0
  87. dtaum1=lambda/diam
  88. c--------------------------------
  89. gmt(1,1)=t1 * dtaum1
  90. gmt(1,2)=u*m * dtaum1
  91. gmt(1,3)=v*m * dtaum1
  92. gmt(1,4)=-1.0D0 * m * dtaum1
  93. c-------------------------------
  94. gmt(2,1)=u*(t1-1.0d0) * dtaum1
  95. gmt(2,2)=(1.0d0+u*u*m) * dtaum1
  96. gmt(2,3)=u*v*m * dtaum1
  97. gmt(2,4)=-1.0D0 * m*u * dtaum1
  98. c-------------------------------
  99. gmt(3,1)=v*(t1-1.0d0) * dtaum1
  100. gmt(3,2)=u*v*m * dtaum1
  101. gmt(3,3)=(1.0d0+v*v*m) * dtaum1
  102. gmt(3,4)=-1.0D0 * m*v * dtaum1
  103. c-------------------------------
  104. gmt(4,1)=t2 * dtaum1
  105. gmt(4,2)=u*h*m * dtaum1
  106. gmt(4,3)=v*h*m * dtaum1
  107. gmt(4,4)=(-m*h+1.0d0) * dtaum1
  108. c-------------------------------
  109. return
  110. end
  111.  
  112.  
  113.  
  114.  
  115.  
  116.  
  117.  
  118.  
  119.  

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