Télécharger fausm2.eso

Retour à la liste

Numérotation des lignes :

  1. C FAUSM2 SOURCE CHAT 05/01/12 23:56:44 5004
  2. SUBROUTINE FAUSM2(NESP,
  3. & GAMG,ROG,PG,UNG,UTG,
  4. & GAMD,ROD,PD,UND,UTD,
  5. & YG,YD,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 KONJA2
  14. C
  15. C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
  16. C
  17. C AUTEUR : S. KUDRIAKOV, 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 All Speeds"
  36. c J.R.Edwards and M.S.Liou
  37. c----------------------------------------------------------------------
  38. c INPUT:
  39. c
  40. c alpha -- parameter of the AUSM+ scheme in the Pressure function;
  41. c ( -3/4 <= alpha <= 3/16 ) (imposed as a parameter)
  42. c
  43. c beta -- parameter of the AUSM+ scheme in the Mach function;
  44. c ( -1/16 <= beta <= 1/2 ) (imposed as a parameter)
  45. c
  46. c wvec_l -- vector of the primitive variables (rho,ux,uy,p) at the
  47. c left cell;
  48. c
  49. c wvec_r -- vector of the primitive variables (rho,ux,uy,p) at the
  50. c right cell;
  51. c
  52. c nvect -- normal vector to the interface (2 components in 2D);
  53. c
  54. c tvect -- tangential vector to the interface;
  55. c
  56. c ga -- ratio of the specific heats (assumed constant)
  57. c----------------------------------------------------------------------
  58. c
  59. c OUTPUT:
  60. c f -- numerical flux-function in the NORMAL DIRECTION
  61. c----------------------------------------------------------------------
  62. c Variable 'cans' (set after the descriptions of all the variables)
  63. c can be used for switching off the low-Mach number additions
  64. c by simply assigning 'cans=0.0d0'
  65. c----------------------------------------------------------------------
  66. IMPLICIT INTEGER(I-N)
  67. integer nesp,i
  68. real*8 gamg,rog,pg,ung,utg
  69. real*8 gamd,rod,pd,und,utd
  70. real*8 alpha,beta
  71. real*8 gm1l,gm1r,f(*)
  72. real*8 un_l, un_r, ut_l, ut_r
  73. real*8 ml,mr,Mplus,Mmin,mmid
  74. real*8 mpl_m, mmin_m,am
  75. real*8 rold_l,pold_l,eold_l
  76. real*8 rold_r,pold_r,eold_r
  77. real*8 Pplus,Pmin,pmid
  78. real*8 hr_l,hr_r,top,bot
  79. real*8 br1,br2,rum
  80. real*8 aleft, arigh,v_inf
  81. real*8 epsil,qq,amw,Mmin1,Mplus1
  82. real*8 fmid,mlw,mrw,termp
  83. real*8 ur_r,ur_l,urm,mhalf,mhalfr
  84. real*8 canc,cellt
  85. real*8 yg(*),yd(*)
  86. parameter(alpha = 0.1875D0, beta = 0.125D0)
  87. parameter(epsil = 1.0d0)
  88. canc=1.0d0
  89. c-------------------------------------------------------------
  90. gm1l=gamg-1.0d0
  91. gm1r=gamd-1.0d0
  92. c----------------------------
  93. rold_l=rog
  94. pold_l=pg
  95. c----------------------------
  96. rold_r=rod
  97. pold_r=pd
  98. c------------------------------------------------------------------
  99. c Computation of the specific total energy on the left and right.
  100. c------------------------------------------------------------------
  101. eold_l=(ung*ung+utg*utg)/2.0d0
  102. eold_l=eold_l+pold_l/(gm1l*rold_l)
  103. eold_r=(und*und+utd*utd)/2.0d0
  104. eold_r=eold_r+pold_r/(gm1r*rold_r)
  105. c------------------------------------------------------------------
  106. c Computing reference velocity and its derivatives
  107. c see Eq.(2) of the Ref.2).
  108. c------------------------------------------------------------------
  109. aleft=sqrt(gamg*pold_l/rold_l)
  110. arigh=sqrt(gamd*pold_r/rold_r)
  111. qq=sqrt(ung*ung+utg*utg)
  112. if(qq .lt. (epsil*v_inf)) then
  113. ur_l = epsil*v_inf
  114. else
  115. ur_l=qq
  116. endif
  117. c-----------------------------
  118. if(ur_l .ge. aleft) then
  119. ur_l=aleft
  120. endif
  121. c-------------------------------------------------------------------
  122. qq=sqrt(und*und+utd*utd)
  123. if(qq .lt. (epsil*v_inf)) then
  124. ur_r = epsil*v_inf
  125. else
  126. ur_r=qq
  127. endif
  128. c-------------------------
  129. if(ur_r .ge. arigh) then
  130. ur_r=arigh
  131. endif
  132. c-------------------------------------------------------------------
  133. c Reference velocity at the interface is taken as an average
  134. c of the reference velocities of the neighbouring cells
  135. c-------------------------------------------------------------------
  136. urm=0.5d0*(ur_l+ur_r)
  137. c-------------------------------------------------------------------
  138. c Computation of the speed of sound;
  139. c numerical speed of sound at the interface is taken as an average
  140. c of the speeds of sounds of the neighbouring cells
  141. c-------------------------------------------------------------------
  142. am=0.5d0*(aleft+arigh)
  143. c-------------------------------------------------------------------
  144. c Computing numerical Mach number; see p.370, under (A1).
  145. c-------------------------------------------------------------------
  146. un_l=ung
  147. un_r=und
  148. ut_l=utg
  149. ut_r=utd
  150. c-------------------------------------------------------------------
  151. ml=un_l/am
  152. mr=un_r/am
  153. mhalf=0.5d0*(un_l+un_r)/am
  154. c-------------------------------
  155. mhalfr=urm/am
  156. c-------------------------------------------------------------------
  157. c Scaling function for the speed of sound;see Eq.(32) of the Ref. 2)
  158. c-------------------------------------------------------------------
  159. top=(1.0d0-mhalfr*mhalfr)*(1.0d0-mhalfr*mhalfr)
  160. top=top*mhalf*mhalf+4.0d0*mhalfr*mhalfr
  161. bot=1.0d0+mhalfr*mhalfr
  162. if(abs(canc-0.0d0).lt.0.000001d0) then
  163. fmid=1.0d0
  164. else
  165. fmid=sqrt(top)/bot
  166. endif
  167. c------------------------------------------------------------------
  168. c 'New' speed of sound 'amw' defined as a product of the scaling
  169. c function 'fmid' and the 'Old' speed of sound 'am'; see (31) of Ref.2)
  170. c------------------------------------------------------------------
  171. amw=fmid*am
  172. mlw=un_l/amw
  173. mrw=un_r/amw
  174. c--------------------------
  175. am=amw
  176. c-------------------------------------------------------------------
  177. c Redefinition of the numerical mach numbers
  178. c See Eqs.(33) and (34) of the Ref. 2)
  179. c-------------------------------------------------------------------
  180. if(abs(canc-0.0d0).lt.0.000001d0) then
  181. top=2.0d0
  182. bot=0.0d0
  183. else
  184. top=1.0d0+mhalfr*mhalfr
  185. bot=1.0d0-mhalfr*mhalfr
  186. endif
  187. ml=0.5d0*(top*mlw+bot*mrw)
  188. mr=0.5d0*(top*mrw+bot*mlw)
  189. c-------------------------------------------------------------------
  190. c Mplus and Mmin are calligraphic lettes M+ and M- from the paper,
  191. c see (19a) and (19b), p.367. of the Ref.1)
  192. c-------------------------------------------------------------------
  193. if(abs(ml) .ge. 1.0d0) then
  194. Mplus=(ml+abs(ml))/2.0d0
  195. else
  196. Mplus=(ml+1.0d0)*(ml+1.0d0)/4.0d0
  197. Mplus=Mplus+beta*(ml*ml-1.0d0)*(ml*ml-1.0d0)
  198. endif
  199. Mplus1=(ml+abs(ml))/2.0d0
  200. c-------------------------------------------------------------------
  201. if(abs(mr) .ge. 1.0d0) then
  202. Mmin=(mr-abs(mr))/2.0d0
  203. else
  204. Mmin=-(mr-1.0d0)*(mr-1.0d0)/4.0d0
  205. Mmin=Mmin-beta*(mr*mr-1.0d0)*(mr*mr-1.0d0)
  206. endif
  207. Mmin1=(mr-abs(mr))/2.0d0
  208. c-----------------------------------------------------------
  209. c mmid is m_{1/2} (notation as in the paper),
  210. c see Eq.(13), p.366 of the Ref.1)
  211. c-----------------------------------------------------------
  212. mmid=Mplus+Mmin
  213. c----------------------------------------------------------------
  214. c computing the main convective variables
  215. c mpl_m is m^{+}_{1/2} (paper's notation) and
  216. c mmin_m is m^{-}_{1/2} (paper's notation),
  217. c see Eq. (A2) on p.370 of the Ref.1)
  218. c----------------------------------------------------------------
  219. termp=(Mmin1-Mmin+Mplus-Mplus1)*(1.0d0/(mhalfr*mhalfr)-1.0d0)
  220. termp=termp*(pold_l-pold_r)/(pold_l/rold_l+pold_r/rold_r)
  221. c-------------------------------------------------------------
  222. if(mmid .ge. 0.0d0) then
  223. mpl_m = mmid
  224. else
  225. mpl_m = 0.0d0
  226. endif
  227. c------------------------------------------------------------------
  228. if(mmid .le. 0.0d0) then
  229. mmin_m = mmid
  230. else
  231. mmin_m = 0.0d0
  232. endif
  233. c---------------------------------------------------------------
  234. c Computing the calligraphic P+ and P- with their derivatives,
  235. c see (21a) & (21b) on p.368 of the Ref.1)
  236. c---------------------------------------------------------------
  237. if(ml .ge. 1.0d0) then
  238. Pplus = 1.0d0
  239. else
  240. if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
  241. Pplus=(ml+1.0d0)*(ml+1.0d0)*(2.0d0-ml)/4.0d0
  242. Pplus=Pplus+alpha*ml*(ml*ml-1.0d0)*(ml*ml-1.0d0)
  243. else
  244. Pplus = 0.0d0
  245. endif
  246. endif
  247. c---------------------------------------------------------------
  248. if(mr .ge. 1.0d0) then
  249. Pmin = 0.0d0
  250. else
  251. if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
  252. Pmin=(mr-1.0d0)*(mr-1.0d0)*(2.0d0+mr)/4.0d0
  253. Pmin=Pmin-alpha*mr*(mr*mr-1.0d0)*(mr*mr-1.0d0)
  254. else
  255. Pmin = 1.0d0
  256. endif
  257. endif
  258. c-------------------------------------------------------------------
  259. c computing pmid - p_{1/2}, see Eq.(20b), p.367 of the Ref.1)
  260. c-------------------------------------------------------------------
  261. pmid=Pplus*pold_l + Pmin*pold_r
  262. c---------------------------------------------------------------------
  263. rum=am*(mpl_m*rold_l+mmin_m*rold_r)+canc*am*termp
  264. c---------------------------------------------------------------------
  265. c Computing numerical fluxes
  266. c---------------------------------------------------------------------
  267. f(1)=rum
  268. c---------------------------------------------------------------------
  269. if(rum .ge. 0.0d0) then
  270. br1=rum*un_l
  271. else
  272. br1=rum*un_r
  273. endif
  274. f(2)=br1+pmid
  275. c-------------------------------------------------------------
  276. if(rum .ge. 0.0d0) then
  277. br2=rum*ut_l
  278. else
  279. br2=rum*ut_r
  280. endif
  281. f(3)=br2
  282. c-------------------------------------------------------------
  283. hr_l=(ung*ung+utg*utg)/2.0d0+gamg*pold_l/gm1l/rold_l
  284. hr_r=(und*und+utd*utd)/2.0d0+gamd*pold_r/gm1r/rold_r
  285. if(rum .ge. 0.0d0) then
  286. f(4)=rum*hr_l
  287. else
  288. f(4)=rum*hr_r
  289. endif
  290. c---------------------------------------------------------------------
  291. do 777 i=1,nesp
  292. if(rum .ge. 0.0d0) then
  293. br1=rum*yg(i)
  294. else
  295. br1=rum*yd(i)
  296. endif
  297. f(4+i)=br1
  298. 777 continue
  299. c----------------------------------------------------------------------
  300. cellt=1.0d0/(0.5d0*abs(un_l+un_r)+am)
  301. c----------------------------------------------------------------------
  302. return
  303. end
  304.  
  305.  
  306.  
  307.  
  308.  
  309.  
  310.  
  311.  
  312.  
  313.  
  314.  
  315.  
  316.  
  317.  
  318.  
  319.  
  320.  
  321.  
  322.  
  323.  
  324.  
  325.  
  326.  
  327.  
  328.  
  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.  

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