Télécharger fauupt.eso

Retour à la liste

Numérotation des lignes :

  1. C FAUUPT SOURCE BECC 10/06/11 21:15:06 6690
  2. SUBROUTINE FAUUPT(nesp,nsca,idim,
  3. & gaml,rhol,pl,unl,utl,uvl,etherl,
  4. & gamr,rhor,pr,unr,utr,uvr,etherr,
  5. & yl,yr,sl,sr,v_inf,f,
  6. & cellt)
  7. C************************************************************************
  8. C
  9. C PROJET : CASTEM 2000
  10. C
  11. C NOM : FAUUPT ('FLUX for Bas Mach' and thermally perfect)
  12. C
  13. C DESCRIPTION :
  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, nsca, idim,iesp
  47. real*8 gaml, rhol, pl, unl, utl, uvl, etherl
  48. & , gamr, rhor, pr, unr, utr, uvr, etherr
  49. & , yl(*), yr(*), v_inf, f(*)
  50. & , cellt, sl(*), sr(*)
  51. C Some functions and parameters
  52. real*8 fau_m4, fau_p5, kp, ku, sigma
  53. C
  54. real*8 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. c Computation of the specific total enthalpies.
  63. c------------------------------------------------------------------
  64. htotl = ((unl*unl) + (utl*utl) + (uvl*uvl)) / 2.0d0
  65. htotl = htotl + etherl
  66. htotl = htotl + (pl / rhol)
  67. htotr = ((unr*unr) + (utr*utr) + (uvr*uvr)) / 2.0d0
  68. htotr = htotr + etherr
  69. htotr = htotr + (pr / rhor)
  70. * write(*,*) htotl, htotr
  71. c------------------------------------------------------------------
  72. c Computing reference velocity
  73. c see Eq.(70) of the Ref.3).
  74. c------------------------------------------------------------------
  75. asonl=sqrt(gaml*pl/rhol)
  76. asonr=sqrt(gamr*pr/rhor)
  77. qq=0.5D0*sqrt((unl*unl)+(utl*utl)+(uvl*uvl))
  78. c------------------------------------------------------------------
  79. if(qq .lt. v_inf) then
  80. uco_l = v_inf
  81. else
  82. uco_l = qq
  83. endif
  84. c-----------------------------
  85. if(uco_l .ge. asonl) then
  86. uco_l = asonl
  87. endif
  88. c-------------------------------------------------------------------
  89. qq=0.5D0*sqrt((unr*unr)+(utr*utr)+(uvr*uvr))
  90. if(qq .lt. v_inf) then
  91. uco_r = v_inf
  92. else
  93. uco_r = qq
  94. endif
  95. c-------------------------
  96. if(uco_r .ge. asonr) then
  97. uco_r = asonr
  98. endif
  99. c-------------------------------------------------------------------
  100. c mco0 is the equivalent of M_{\inf} in Ref 3.
  101. c-------------------------------------------------------------------
  102. mco0 = max((uco_r / asonr), (uco_l / asonl))
  103. c write(*,*) 'mco0 = ', mco0
  104. c-------------------------------------------------------------------
  105. c Computation of the speed of sound;
  106. c numerical speed of sound at the interface is taken as an average
  107. c of the speeds of sounds of the neighbouring cells
  108. c-------------------------------------------------------------------
  109. asonm=0.5d0*(asonl + asonr)
  110. c-------------------------------------------------------------------
  111. c Computing numerical Mach number
  112. c ml, mr, mbar2 (used in the term for Mu, which does not involves
  113. c the cut-off speed, see equation (70), ref 3.)
  114. c Then, even if mc0 = 1, the term involving kp acts!
  115. c-------------------------------------------------------------------
  116. ml = unl / asonm
  117. mr = unr / asonm
  118. mbar2 = ((unl*unl)+(unr*unr))/(2.0d0*asonm*asonm)
  119. * write(*,*) asonm
  120. * write(*,*) ml, mr, mbar2
  121. c-------------------------------
  122. c-------------------------------------------------------------------
  123. c Scaling function for the speed of sound
  124. c We redefine mzero (see (71), ref 3).
  125. c-------------------------------------------------------------------
  126. C
  127. mzero = max((mco0 * mco0), mbar2)
  128. mzero = min(1.0D0, mzero)
  129. mzero = mzero ** 0.5D0
  130. C
  131. fson = mzero * (2.0D0 - mzero)
  132. alpha = (3.0D0 / 16.0D0) * (-4.0D0 + (5.0D0 * fson * fson))
  133. C write(*,*) 'alpha = ', alpha
  134. C-------------------------------------------------------------------
  135. C mhalf (see (73), ref 3)
  136. C First, contribution of Mu; then contribution of ml and mr.
  137. C-------------------------------------------------------------------
  138. C
  139. rhom = 0.5D0 * (rhol + rhor)
  140. mhalf = max(1.0D0 - (sigma * mbar2), 0.0D0)
  141. mhalf = -1.0D0 * mhalf * kp * (pr - pl) /
  142. & (fson * rhom * asonm * asonm)
  143. C
  144. mhalf = mhalf + fau_m4(ml,1.0D0) + fau_m4(mr,-1.0D0)
  145. c write(*,*) 'mhalf = ', mhalf
  146. uhalf = mhalf * asonm
  147. c write(*,*) 'unl, unr ' , unl, unr
  148. c write(*,*) 'mhalf, uhalf ', mhalf, uhalf
  149. C-------------------------------------------------------------------
  150. C phalf
  151. C-------------------------------------------------------------------
  152. fp5p = fau_p5(ml,1.0D0,alpha)
  153. fp5m = fau_p5(mr,-1.0D0,alpha)
  154. phalf = -1.0D0 * ku * fp5m * fp5p * (rhol + rhor) * fson * asonm *
  155. & (unr - unl)
  156. phalf = phalf + (fp5p * pl) + (fp5m * pr)
  157. c write(*,*) 'phalf = ', phalf
  158. c---------------------------------------------------------------------
  159. c Computing numerical fluxes
  160. c---------------------------------------------------------------------
  161. if(uhalf .ge. 0.0d0) then
  162. f(1) = uhalf * rhol
  163. f(2) = (f(1) * unl) + phalf
  164. f(3) = f(1) * utl
  165. if (idim .eq. 3) then
  166. f(4) = f(1) * uvl
  167. endif
  168. f(2 + idim) = f(1) * htotl
  169. do iesp = 1, nesp, 1
  170. f(2 + idim + iesp) = f(1) * yl(iesp)
  171. enddo
  172. do iesp = 1, nsca, 1
  173. f(2 + idim + nesp + iesp) = f(1) * sl(iesp)
  174. enddo
  175. else
  176. f(1) = uhalf * rhor
  177. f(2) = (f(1) * unr) + phalf
  178. f(3) = f(1) * utr
  179. if (idim .eq. 3) then
  180. f(4) = f(1) * uvr
  181. endif
  182. f(idim + 2) = f(1) * htotr
  183. do iesp = 1, nesp, 1
  184. f(2 + idim + iesp) = f(1) * yr(iesp)
  185. enddo
  186. do iesp = 1, nsca, 1
  187. f(2 + idim + nesp + iesp) = f(1) * sr(iesp)
  188. enddo
  189. endif
  190. c----------------------------------------------------------------------
  191. cellt = max(abs(unl),abs(unr)) + asonm
  192. c write(*,*) 'cellt = ', cellt
  193. c----------------------------------------------------------------------
  194. return
  195. end
  196.  
  197.  

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