Télécharger fausm4.eso

Retour à la liste

Numérotation des lignes :

  1. C FAUSM4 SOURCE CHAT 05/01/12 23:56:50 5004
  2. C FAUSM3 SOURCE BECC 01/07/27 21:16:20 4170
  3. SUBROUTINE FAUSM4(NESP,
  4. & GAMG,ROG,PG,UNG,UTG,UVG,
  5. & GAMD,ROD,PD,UND,UTD,UVD,
  6. & YG,YD,V_INF,F,
  7. & CELLT)
  8. C************************************************************************
  9. C
  10. C PROJET : CASTEM 2000
  11. C
  12. C NOM : FAUSM4
  13. C
  14. C DESCRIPTION : Voir KONJA2
  15. C
  16. C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
  17. C
  18. C AUTEUR : S. KUDRIAKOV, DM2S/SFME/LTMF
  19. C
  20. C************************************************************************
  21. C
  22. c----------------------------------------------------------------------
  23. c GENERAL DESCRIPTION:
  24. c This subroutine provides the numerical flux function
  25. c defined at the cell interface; this flux is given in
  26. c the NORMAL DIRECTION (nvect)
  27. c The low-mach number corrections are made for the flux functions
  28. c
  29. c EQUATIONS: 3D Euler equations of gas dynamics
  30. c
  31. c
  32. c REFERENCE: 1) JCP, 129, 364-382 (1996)
  33. c " A Sequel to AUSM: AUSM+ ";
  34. c M.S.Liou
  35. c 2) AIAA Journal, Sept. 1998
  36. c "Low-Diffusion Flux-Splitting Methods for Flows at All Speeds"
  37. c J.R.Edwards and M.S.Liou
  38. c----------------------------------------------------------------------
  39. c INPUT:
  40. c
  41. c NESP -- number of species in the multispecies gas
  42. c
  43. c alpha -- parameter of the AUSM+ scheme in the Pressure function;
  44. c ( -3/4 <= alpha <= 3/16 ) (imposed as a parameter)
  45. c
  46. c beta -- parameter of the AUSM+ scheme in the Mach function;
  47. c ( -1/16 <= beta <= 1/2 ) (imposed as a parameter)
  48. c
  49. c (gamg,rhog,pg,ung,utg,uvg) -- vector of the primitive variables
  50. c at the left cell;
  51. c
  52. c (gamd,rhod,pd,und,utd,uvd) -- vector of the primitive variables
  53. c at the right cell;
  54. c
  55. c yg -- vector of the mass fractionc of the species
  56. c at the left cell;
  57. c
  58. c yd -- vector of the mass fractions of the species
  59. c at the right cell;
  60. c
  61. c v_inf -- parameter for choosing the reference velocity
  62. c when the magnitude of the physical velocity
  63. c is close to zero
  64. c----------------------------------------------------------------------
  65. c
  66. c OUTPUT:
  67. c f -- numerical flux-function in the NORMAL DIRECTION
  68. c----------------------------------------------------------------------
  69. c Variable 'cans' (set after the descriptions of all the variables)
  70. c can be used for switching off the low-Mach number additions
  71. c by simply assigning 'cans=0.0d0'
  72. c----------------------------------------------------------------------
  73. IMPLICIT INTEGER(I-N)
  74. integer nesp,i
  75. real*8 gamg,rog,pg,ung,utg,uvg
  76. real*8 gamd,rod,pd,und,utd,uvd
  77. real*8 alpha,beta,upr_l,upr_r
  78. real*8 gm1l,gm1r,f(*)
  79. real*8 un_l, un_r, ut_l, ut_r
  80. real*8 uv_l, uv_r
  81. real*8 ml,mr,Mplus,Mmin,mmid
  82. real*8 mpl_m, mmin_m,am
  83. real*8 rold_l,pold_l,eold_l
  84. real*8 rold_r,pold_r,eold_r
  85. real*8 Pplus,Pmin,pmid
  86. real*8 hr_l,hr_r,top,bot
  87. real*8 br1,br2
  88. real*8 aleft, arigh
  89. real*8 epsil,qq,amw,Mmin1,Mplus1
  90. real*8 fmid,mlw,mrw,termp
  91. real*8 ur_r,ur_l,urm,mhalf,mhalfr
  92. real*8 canc,cellt,v_inf,rum
  93. real*8 yg(*),yd(*)
  94. parameter(alpha = 0.1875D0, beta = 0.125D0)
  95. parameter(epsil = 1.0d0)
  96. canc=1.0d0
  97. upr_l=0.0d0
  98. upr_r=0.0d0
  99. c-------------------------------------------------------------
  100. gm1l=gamg-1.0d0
  101. gm1r=gamd-1.0d0
  102. c----------------------------
  103. rold_l=rog
  104. pold_l=pg
  105. c----------------------------
  106. rold_r=rod
  107. pold_r=pd
  108. c------------------------------------------------------------------
  109. c Computation of the specific total energy on the left and right.
  110. c------------------------------------------------------------------
  111. eold_l=(ung*ung+utg*utg+uvg*uvg)/2.0d0
  112. eold_l=eold_l+pold_l/(gm1l*rold_l)
  113. eold_r=(und*und+utd*utd+uvd*uvd)/2.0d0
  114. eold_r=eold_r+pold_r/(gm1r*rold_r)
  115. c------------------------------------------------------------------
  116. c Computing reference velocity and its derivatives
  117. c see Eq.(2) of the Ref.2).
  118. c------------------------------------------------------------------
  119. aleft=sqrt(gamg*pold_l/rold_l)
  120. arigh=sqrt(gamd*pold_r/rold_r)
  121. qq=sqrt(ung*ung+utg*utg+uvg*uvg)
  122. if(qq .lt. (epsil*v_inf)) then
  123. ur_l = epsil*v_inf
  124. else
  125. ur_l=qq
  126. endif
  127. c-------------------------------------------------------------------
  128. if(ur_l .ge. aleft) then
  129. ur_l=aleft
  130. endif
  131. c-------------------------------------------------------------------
  132. if(ur_l .lt. upr_l) ur_l=upr_l
  133. c-------------------------------------------------------------------
  134. qq=sqrt(und*und+utd*utd+uvd*uvd)
  135. if(qq .lt. (epsil*v_inf)) then
  136. ur_r = epsil*v_inf
  137. else
  138. ur_r=qq
  139. endif
  140. c------------------------------------------------------------------
  141. if(ur_r .ge. arigh) then
  142. ur_r=arigh
  143. endif
  144. c------------------------------------------------------------------
  145. if(ur_r .lt. upr_r) ur_r=upr_r
  146. c-------------------------------------------------------------------
  147. c Reference velocity at the interface is taken as an average
  148. c of the reference velocities of the neighbouring cells
  149. c-------------------------------------------------------------------
  150. urm=0.5d0*(ur_l+ur_r)
  151. c-------------------------------------------------------------------
  152. c Computation of the speed of sound;
  153. c numerical speed of sound at the interface is taken as an average
  154. c of the speeds of sounds of the neighbouring cells
  155. c-------------------------------------------------------------------
  156. am=0.5d0*(aleft+arigh)
  157. c-------------------------------------------------------------------
  158. c Computing numerical Mach number; see p.370, under (A1).
  159. c-------------------------------------------------------------------
  160. un_l=ung
  161. un_r=und
  162. c--------------
  163. ut_l=utg
  164. ut_r=utd
  165. c--------------
  166. uv_l=uvg
  167. uv_r=uvd
  168. c-------------------------------------------------------------------
  169. ml=un_l/am
  170. mr=un_r/am
  171. mhalf=0.5d0*(un_l+un_r)/am
  172. c-------------------------------
  173. mhalfr=urm/am
  174. c-------------------------------------------------------------------
  175. c Scaling function for the speed of sound;see Eq.(32) of the Ref. 2)
  176. c-------------------------------------------------------------------
  177. top=(1.0d0-mhalfr*mhalfr)*(1.0d0-mhalfr*mhalfr)
  178. top=top*mhalf*mhalf+4.0d0*mhalfr*mhalfr
  179. bot=1.0d0+mhalfr*mhalfr
  180. if(abs(canc-0.0d0).lt.0.000001d0) then
  181. fmid=1.0d0
  182. else
  183. fmid=sqrt(top)/bot
  184. endif
  185. c------------------------------------------------------------------
  186. c 'New' speed of sound 'amw' defined as a product of the scaling
  187. c function 'fmid' and the 'Old' speed of sound 'am'; see (31) of Ref.2)
  188. c------------------------------------------------------------------
  189. amw=fmid*am
  190. mlw=un_l/amw
  191. mrw=un_r/amw
  192. c--------------------------
  193. am=amw
  194. c-------------------------------------------------------------------
  195. c Redefinition of the numerical mach numbers
  196. c See Eqs.(33) and (34) of the Ref. 2)
  197. c-------------------------------------------------------------------
  198. if(abs(canc-0.0d0).lt.0.000001d0) then
  199. top=2.0d0
  200. bot=0.0d0
  201. else
  202. top=1.0d0+mhalfr*mhalfr
  203. bot=1.0d0-mhalfr*mhalfr
  204. endif
  205. ml=0.5d0*(top*mlw+bot*mrw)
  206. mr=0.5d0*(top*mrw+bot*mlw)
  207. c-------------------------------------------------------------------
  208. c Mplus and Mmin are calligraphic lettes M+ and M- from the paper,
  209. c see (19a) and (19b), p.367. of the Ref.1)
  210. c-------------------------------------------------------------------
  211. if(abs(ml) .ge. 1.0d0) then
  212. Mplus=(ml+abs(ml))/2.0d0
  213. else
  214. Mplus=(ml+1.0d0)*(ml+1.0d0)/4.0d0
  215. Mplus=Mplus+beta*(ml*ml-1.0d0)*(ml*ml-1.0d0)
  216. endif
  217. Mplus1=(ml+abs(ml))/2.0d0
  218. c-------------------------------------------------------------------
  219. if(abs(mr) .ge. 1.0d0) then
  220. Mmin=(mr-abs(mr))/2.0d0
  221. else
  222. Mmin=-(mr-1.0d0)*(mr-1.0d0)/4.0d0
  223. Mmin=Mmin-beta*(mr*mr-1.0d0)*(mr*mr-1.0d0)
  224. endif
  225. Mmin1=(mr-abs(mr))/2.0d0
  226. c-----------------------------------------------------------
  227. c mmid is m_{1/2} (notation as in the paper),
  228. c see Eq.(13), p.366 of the Ref.1)
  229. c-----------------------------------------------------------
  230. mmid=Mplus+Mmin
  231. c----------------------------------------------------------------
  232. c computing the main convective variables
  233. c mpl_m is m^{+}_{1/2} (paper's notation) and
  234. c mmin_m is m^{-}_{1/2} (paper's notation),
  235. c see Eq. (A2) on p.370 of the Ref.1)
  236. c----------------------------------------------------------------
  237. termp=(Mmin1-Mmin+Mplus-Mplus1)*(1.0d0/(mhalfr*mhalfr)-1.0d0)
  238. termp=termp*(pold_l-pold_r)/(pold_l/rold_l+pold_r/rold_r)
  239. c-------------------------------------------------------------
  240. if(mmid .ge. 0.0d0) then
  241. mpl_m = mmid
  242. else
  243. mpl_m = 0.0d0
  244. endif
  245. c------------------------------------------------------------------
  246. if(mmid .le. 0.0d0) then
  247. mmin_m = mmid
  248. else
  249. mmin_m = 0.0d0
  250. endif
  251. c---------------------------------------------------------------
  252. c Computing the calligraphic P+ and P- with their derivatives,
  253. c see (21a) & (21b) on p.368 of the Ref.1)
  254. c---------------------------------------------------------------
  255. if(ml .ge. 1.0d0) then
  256. Pplus = 1.0d0
  257. else
  258. if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
  259. Pplus=(ml+1.0d0)*(ml+1.0d0)*(2.0d0-ml)/4.0d0
  260. Pplus=Pplus+alpha*ml*(ml*ml-1.0d0)*(ml*ml-1.0d0)
  261. else
  262. Pplus = 0.0d0
  263. endif
  264. endif
  265. c---------------------------------------------------------------
  266. if(mr .ge. 1.0d0) then
  267. Pmin = 0.0d0
  268. else
  269. if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
  270. Pmin=(mr-1.0d0)*(mr-1.0d0)*(2.0d0+mr)/4.0d0
  271. Pmin=Pmin-alpha*mr*(mr*mr-1.0d0)*(mr*mr-1.0d0)
  272. else
  273. Pmin = 1.0d0
  274. endif
  275. endif
  276. c-------------------------------------------------------------------
  277. c computing pmid - p_{1/2}, see Eq.(20b), p.367 of the Ref.1)
  278. c-------------------------------------------------------------------
  279. pmid=Pplus*pold_l + Pmin*pold_r
  280. c-------------------------------------------------------------------
  281. rum=am*(mpl_m*rold_l+mmin_m*rold_r)+canc*am*termp
  282. c-------------------------------------------------------------------
  283. c Computing numerical fluxes
  284. c-------------------------------------------------------------------
  285. f(1)=rum
  286. c-------------------------------------------------------------------
  287. if(rum .ge. 0.0d0) then
  288. br1=rum*un_l
  289. else
  290. br1=rum*un_r
  291. endif
  292. f(2)=br1+pmid
  293. c-------------------------------------------------------------
  294. if(rum .ge. 0.0d0) then
  295. br2=rum*ut_l
  296. else
  297. br2=rum*ut_r
  298. endif
  299. f(3)=br2
  300. c-------------------------------------------------------------
  301. if(rum .ge. 0.0d0) then
  302. br2=rum*uv_l
  303. else
  304. br2=rum*uv_r
  305. endif
  306. f(4)=br2
  307. c-------------------------------------------------------------
  308. hr_l=(ung*ung+utg*utg+uvg*uvg)/2.0d0+gamg*pold_l/gm1l/rold_l
  309. hr_r=(und*und+utd*utd+uvd*uvd)/2.0d0+gamd*pold_r/gm1r/rold_r
  310. if(rum .ge. 0.0d0) then
  311. f(5)=rum*hr_l
  312. else
  313. f(5)=rum*hr_r
  314. endif
  315. c---------------------------------------------------------------------
  316. do 777 i=1,nesp
  317. if(rum .ge. 0.0d0) then
  318. br1=rum*yg(i)
  319. else
  320. br1=rum*yd(i)
  321. endif
  322. f(5+i)=br1
  323. 777 continue
  324. c----------------------------------------------------------------------
  325. cellt=1.0d0/(0.5d0*abs(un_l+un_r)+am)
  326. c----------------------------------------------------------------------
  327. return
  328. end
  329.  
  330.  
  331.  
  332.  
  333.  
  334.  
  335.  
  336.  
  337.  
  338.  
  339.  
  340.  
  341.  
  342.  
  343.  
  344.  
  345.  
  346.  
  347.  
  348.  
  349.  
  350.  
  351.  
  352.  
  353.  
  354.  
  355.  
  356.  
  357.  
  358.  
  359.  
  360.  
  361.  
  362.  
  363.  
  364.  
  365.  
  366.  
  367.  
  368.  
  369.  
  370.  
  371.  
  372.  
  373.  
  374.  
  375.  
  376.  
  377.  
  378.  
  379.  
  380.  
  381.  

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