Télécharger fau_up.eso

Retour à la liste

Numérotation des lignes :

  1. C FAU_UP SOURCE BECC 10/05/05 21:15:07 6674
  2. SUBROUTINE FAU_UP(nesp,
  3. & gaml,rhol,pl,unl,utl,
  4. & gamr,rhor,pr,unr,utr,
  5. & yl,yr,v_inf,f,
  6. & cellt)
  7. C************************************************************************
  8. C
  9. C PROJET : CASTEM 2000
  10. C
  11. C NOM : FAUSM2 ('FLUX for Bas Mach')
  12. C
  13. C DESCRIPTION : Voir
  14. C
  15. C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
  16. C
  17. C AUTEUR : A. Beccantini, DM2S/SFME/LTMF
  18. C
  19. C************************************************************************
  20. C
  21. c----------------------------------------------------------------------
  22. c GENERAL DESCRIPTION:
  23. c This subroutine provides the numerical flux function
  24. c defined at the cell interface; this flux is given in
  25. c the NORMAL DIRECTION (nvect)
  26. c The low-mach number corrections are made for the flux functions
  27. c
  28. c EQUATIONS: 2D Euler equations of gas dynamics
  29. c
  30. c
  31. c REFERENCE: 1) JCP, 129, 364-382 (1996)
  32. c " A Sequel to AUSM: AUSM+ ";
  33. c M.S.Liou
  34. c 2) AIAA Journal, Sept. 1998
  35. c "Low-Diffusion Flux-Splitting Methods for Flows at
  36. C All Speeds"
  37. c J.R.Edwards and M.S.Liou
  38. C 3) JCP, 214, 137-170 (2006)
  39. C A Sequel to AUSM, Part 2 AUSM+-up for all speeds
  40. c M.S.Liou
  41. C
  42. c----------------------------------------------------------------------
  43. c INPUT:
  44. c
  45. implicit integer(i-n)
  46. integer nesp
  47. real*8 gaml, rhol, pl, unl, utl
  48. & , gamr, rhor, pr, unr, utr
  49. & , yl(*), yr(*), v_inf, f(*)
  50. & , cellt
  51. C Some functions and parameters
  52. real*8 fau_m4, fau_p5, kp, ku, sigma
  53. C
  54. real*8 gm1l, gm1r, htotl, htotr, asonl, asonr, qq
  55. & , uco_l, uco_r, mco0, asonm, ml, mr, mbar2
  56. & , mzero, fson, alpha, rhom, mhalf, uhalf
  57. & , fp5p, fp5m, phalf
  58.  
  59. PARAMETER(kp=0.25d0, ku=0.75d0, sigma=1.0D0)
  60. C PARAMETER(kp=0.25d0, ku=0.0d0, sigma=1.0D0)
  61. c-------------------------------------------------------------
  62. gm1l=gaml-1.0d0
  63. gm1r=gamr-1.0d0
  64. c------------------------------------------------------------------
  65. c Computation of the specific total enthalpies.
  66. c------------------------------------------------------------------
  67. htotl = ((unl*unl) + (utl*utl)) / 2.0d0
  68. htotl = htotl + (pl / (gm1l*rhol))
  69. htotl = htotl + (pl / rhol)
  70. htotr = ((unr*unr) + (utr*utr)) / 2.0d0
  71. htotr = htotr + (pr / (gm1r*rhor))
  72. htotr = htotr + (pr / rhor)
  73. * write(*,*) htotl, htotr
  74. c------------------------------------------------------------------
  75. c Computing reference velocity
  76. c see Eq.(70) of the Ref.3).
  77. c------------------------------------------------------------------
  78. asonl=sqrt(gaml*pl/rhol)
  79. asonr=sqrt(gamr*pr/rhor)
  80. qq=sqrt((unl*unl)+(utl*utl))
  81. c------------------------------------------------------------------
  82. if(qq .lt. v_inf) then
  83. uco_l = v_inf
  84. else
  85. uco_l = qq
  86. endif
  87. c-----------------------------
  88. if(uco_l .ge. asonl) then
  89. uco_l = asonl
  90. endif
  91. c-------------------------------------------------------------------
  92. qq=sqrt((unr*unr)+(utr*utr))
  93. if(qq .lt. v_inf) then
  94. uco_r = v_inf
  95. else
  96. uco_r = qq
  97. endif
  98. c-------------------------
  99. if(uco_r .ge. asonr) then
  100. uco_r = asonr
  101. endif
  102. c-------------------------------------------------------------------
  103. c mco0 is the equivalent of M_{\inf} in Ref 3.
  104. c-------------------------------------------------------------------
  105. mco0 = max((uco_r / asonr), (uco_l / asonl))
  106. c write(*,*) 'mco0 = ', mco0
  107. c-------------------------------------------------------------------
  108. c Computation of the speed of sound;
  109. c numerical speed of sound at the interface is taken as an average
  110. c of the speeds of sounds of the neighbouring cells
  111. c-------------------------------------------------------------------
  112. asonm=0.5d0*(asonl + asonr)
  113. c-------------------------------------------------------------------
  114. c Computing numerical Mach number
  115. c ml, mr, mbar2 (used in the term for Mu, which does not involves
  116. c the cut-off speed, see equation (70), ref 3.)
  117. c Then, even if mc0 = 1, the term involving kp acts!
  118. c-------------------------------------------------------------------
  119. ml = unl / asonm
  120. mr = unr / asonm
  121. mbar2 = ((unl*unl)+(unr*unr))/(2.0d0*asonm*asonm)
  122. * write(*,*) asonm
  123. * write(*,*) ml, mr, mbar2
  124. c-------------------------------
  125. c-------------------------------------------------------------------
  126. c Scaling function for the speed of sound
  127. c We redefine mzero (see (71), ref 3).
  128. c-------------------------------------------------------------------
  129. C
  130. mzero = max((mco0 * mco0), mbar2)
  131. mzero = min(1.0D0, mzero)
  132. mzero = mzero ** 0.5D0
  133. C
  134. fson = mzero * (2.0D0 - mzero)
  135. alpha = (3.0D0 / 16.0D0) * (-4.0D0 + (5.0D0 * fson * fson))
  136. C write(*,*) 'alpha = ', alpha
  137. C-------------------------------------------------------------------
  138. C mhalf (see (73), ref 3)
  139. C First, contribution of Mu; then contribution of ml and mr.
  140. C-------------------------------------------------------------------
  141. C
  142. rhom = 0.5D0 * (rhol + rhor)
  143. mhalf = max(1.0D0 - (sigma * mbar2), 0.0D0)
  144. mhalf = -1.0D0 * mhalf * kp * (pr - pl) /
  145. & (fson * rhom * asonm * asonm)
  146. C
  147. mhalf = mhalf + fau_m4(ml,1.0D0) + fau_m4(mr,-1.0D0)
  148. c write(*,*) 'mhalf = ', mhalf
  149. uhalf = mhalf * asonm
  150. c write(*,*) 'unl, unr ' , unl, unr
  151. c write(*,*) 'mhalf, uhalf ', mhalf, uhalf
  152. C-------------------------------------------------------------------
  153. C phalf
  154. C-------------------------------------------------------------------
  155. fp5p = fau_p5(ml,1.0D0,alpha)
  156. fp5m = fau_p5(mr,-1.0D0,alpha)
  157. phalf = -1.0D0 * ku * fp5m * fp5p * (rhol + rhor) * fson * asonm *
  158. & (unr - unl)
  159. phalf = phalf + (fp5p * pl) + (fp5m * pr)
  160. c write(*,*) 'phalf = ', phalf
  161. c---------------------------------------------------------------------
  162. c Computing numerical fluxes
  163. c---------------------------------------------------------------------
  164. if(uhalf .ge. 0.0d0) then
  165. f(1) = uhalf * rhol
  166. f(2) = (f(1) * unl) + phalf
  167. f(3) = f(1) * utl
  168. f(4) = f(1) * htotl
  169. do iesp = 1, nesp, 1
  170. f(4 + iesp) = f(1) * yl(iesp)
  171. enddo
  172. else
  173. f(1) = uhalf * rhor
  174. f(2) = (f(1) * unr) + phalf
  175. f(3) = f(1) * utr
  176. f(4) = f(1) * htotr
  177. do iesp = 1, nesp, 1
  178. f(4 + iesp) = f(1) * yr(iesp)
  179. enddo
  180. endif
  181. c----------------------------------------------------------------------
  182. cellt = 1.0d0 / ((0.5d0*abs(unl+unr))+asonm)
  183. c write(*,*) 'cellt = ', cellt
  184. c----------------------------------------------------------------------
  185. return
  186. end
  187.  
  188.  

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