Télécharger fausre.eso

Retour à la liste

Numérotation des lignes :

  1. C FAUSRE SOURCE BECC 11/01/05 21:15:01 6836
  2. subroutine fausre(ndim, nesp, nLHS, nordpo,
  3. & reacoe, aprop, Runiv,
  4. & Tmaxcv, v_inf,
  5. & cnna, ctna, ct1na,
  6. & wvec_l,
  7. & wvec_r,
  8. & wwork, y, acv_l, acv_r, acv_b, acv_w ,
  9. & flu, f_lag, cellt,
  10. & logan)
  11. C
  12. c************************************************************************
  13. c
  14. c name : fausre
  15. c
  16. c description : Reactive Euler equations for a mixture of reactive
  17. C thermally perfect gases.
  18. c
  19. c AUSM-up, implemented in a
  20. C "Discrete Equation Method" fashion.
  21. c
  22. c Voir:
  23. c 1) Beccantini, Studer, IJNMF 2010
  24. C 2) Metayer, Massoni, Saurel, "Modelling
  25. C evaporations fronts with reactive Riemann
  26. C solvers.", JCP 205, 2005
  27. c
  28. c language : fortran 77
  29. c
  30. c author : a. beccantini den/dm2s/sfme/ltmf
  31. c
  32. c************************************************************************
  33. c
  34. c called sub : prithe, racre2, flufun
  35. C
  36. c************************************************************************
  37. c
  38. c**** input
  39. C
  40. C ndim = dimension of the domain
  41. c
  42. c nesp = species involved in the Euler equations
  43. C
  44. C nLHS = species involved in the LHS of the global
  45. C reaction
  46. c
  47. c nordp0 = polynomial order of the cv
  48. c
  49. c reacoe = coefficient involved in the global reaction
  50. c (positive in the LHS, negative in the RHS...
  51. c the first one should be positive)
  52. c
  53. c aprop = gas properties (cv, molar mass, formation
  54. C energy)
  55. c
  56. C Runiv = universal gas constant
  57. c
  58. C Tmaxcv = Tmax for cv expansion
  59. C
  60. C V_inf = cut-off coefficient (velocity)
  61. C
  62. C cnna, ctna, ct1na = scalar products
  63. C (\vec{n},\vec{n}_{\rm alpha})
  64. C (\vec{t},\vec{n}_{\rm alpha})
  65. C (\vec{t1},\vec{n}_{\rm alpha})
  66. C HP. We suppose that \vec{n}_{\rm alpha} is
  67. C always from the burnt to the unburnt state
  68. C
  69. c wvec_l = left primitive variables
  70. c
  71. c wvec_r = right primitive variables
  72. C
  73. C NB Primitive variables =
  74. C
  75. C w(1) = rho
  76. C w(2) = ux
  77. C w(3) = uy
  78. C w(4) = uz
  79. C w(2+ndim) = p
  80. C w(3+ndim) = csi
  81. C w(3+ndim+1) = yf(1)
  82. C ...
  83. C w(3+ndim+nLHS) = yf(nLHS)
  84. C w(3+ndim+nLHS+1) = yi(nLHS+1)
  85. C ...
  86. C w(3+ndim+nesp-1) = yi(nesp-1)
  87. C w(3+ndim+nesp) = dy1t = yi(1) - yf(1)
  88. C w(4+ndim+nesp) = k0
  89. C
  90. C For ndim = 1, the total number of variables
  91. C is 4 + 1 + nesp = 5 + nesp.
  92. c
  93. c wwork, y, acv_l,
  94. C acv_r, acv_b,
  95. C acv_w = temporary work vectors
  96. C wwork(4+ndim+nesp), y(1:nesp),
  97. C acv_l(0:nordpo),
  98. C acv_r(0:nordpo),
  99. C acv_b(0:nordpo),
  100. C acv_w (0:nordpo)
  101. C
  102. c**** output
  103. c
  104. c flu = Eulerian interfacial flux in (n,t1,t2), i.e.
  105. c rho*un mass
  106. C rho*un*un + p momentum along n
  107. C rho*un*ut1 momentum along t1
  108. * rho*un*ut2 momentum along t2
  109. C rho*un*ht energy
  110. c 0.0 alpha
  111. C rho*un*y_{1,f} mass * y_{1,f}
  112. C ...
  113. C rho*un*y_{nLHS+1,i} mass * y_{nLHS+1,i}
  114. c ...
  115. C rho*un*dy_{1,t} mass * dy_{1,t}
  116. C rho*un*k0 mass * k0
  117. C
  118. C NB. flu(1:(4+ndim+nesp))
  119. c
  120. C f_lag = Lagrangian interfacial flux on the reactive
  121. C discontinuty in (n,t1,t2), i.e.
  122. c rho*(un - D_resh) mass
  123. C rho*un*un + p - (D_resh rho un) momentum along n
  124. C rho*(un - D_resh)*ut1 momentum along t1
  125. C rho*(un - D_resh)*ut2 momentum along t2
  126. C rho*un*ht - (D_resh rho et) energy
  127. c - D_resh alpha
  128. c rho*(un - D_resh)*y_{1,f} mass * y_{1,f}
  129. C ...
  130. C
  131. C NB. f_lag(1:(4+ndim+nesp))
  132. C
  133. c
  134. c cellt = stability condition, i.e.
  135. c
  136. c dt/dx < cellt (dimension 1/velocity)
  137. c
  138. c logan = si true, anomaly detected
  139. c
  140. c************************************************************************
  141. c
  142. C 10/04/2008 created
  143. C 19/08/2008 distinction between reactive and non-reactive case
  144. C
  145. C
  146. c************************************************************************
  147. c
  148. c n.b.: all variables are declared
  149. c
  150. c implicit none
  151. implicit integer(i-n)
  152. integer ndim, nesp, nLHS, nordpo
  153. integer nesp1, nsca1
  154. integer icomp
  155. real*8 reacoe(nesp), aprop(0:(nordpo+2),nesp), Runiv, Tmaxcv
  156. & , wvec_l(4+ndim+nesp), wvec_r(4+ndim+nesp)
  157. & , flu(4+ndim+nesp)
  158. & , f_lag(4+ndim+nesp), cellt
  159. & , cnna, ctna, ct1na
  160. real*8 wwork(4+ndim+nesp)
  161. & , acv_l(0:nordpo), acv_r(0:nordpo), acv_b(0:nordpo)
  162. & , acv_w(0:nordpo)
  163. $ , y(1:nesp)
  164. real*8 gam_l, Rgas_l, cv_l, hf_l, rho_l, u_l(1:3), p_l, T_l
  165. & , et_l, csi_l
  166. & , gam_r, Rgas_r, cv_r, hf_r, rho_r, u_r(1:3), p_r, T_r
  167. & , et_r, csi_r
  168. & , un_r, un_l, ut_l, ut_r, ut1_l, ut1_r
  169. real*8 xst, flurho, coefy, omcoey
  170. & , v_inf
  171. C
  172. logical logan
  173. logical logdeb, logk0n
  174. C logdeb = debugging ?
  175. C logk0n = .true. : instead of k0, we take k0 * abs(cnna)
  176. C .false. : instead of k0, we take k0 * abs(cnna)
  177. parameter (logdeb = .false., logk0n = .true.)
  178. C
  179. C************************************************************************
  180. C******* NON REACTIVE CASE ONLY *****************************************
  181. C************************************************************************
  182. C
  183. call prithe(ndim, nesp, nLHS, nordpo, aprop, reacoe, Runiv,
  184. $ Tmaxcv,wvec_l,gam_l, Rgas_l, cv_l, acv_l, hf_l, rho_l,
  185. $ u_l, p_l, T_l,y, et_l, csi_l, logan)
  186. if (logan) then
  187. write(*,'(/A8)') 'fvlhre.f'
  188. write(*,*) 'ANOMALY DETECTED 01'
  189. goto 9999
  190. endif
  191. C
  192. call prithe(ndim, nesp, nLHS, nordpo, aprop, reacoe, Runiv,
  193. $ Tmaxcv,wvec_r,gam_r, Rgas_r, cv_r, acv_r, hf_r, rho_r,
  194. $ u_r, p_r, T_r,y, et_r, csi_r, logan)
  195. if (logan) then
  196. write(*,'(/A8)') 'fvlhre.f'
  197. write(*,*) 'ANOMALY DETECTED 02'
  198. goto 9999
  199. endif
  200. C
  201. un_l = u_l(1)
  202. un_r = u_r(1)
  203. C
  204. C d_l d_b un_b=un_a d_a d_r
  205. C
  206. xst = 0.0D0
  207. C
  208. nesp1 = 0
  209. nsca1 = 0
  210. C
  211. if (ndim .eq. 1) then
  212. ut_l = 0.0D0
  213. ut_r = 0.0D0
  214. ut1_l = 0.0D0
  215. ut1_r = 0.0D0
  216. elseif(ndim .eq. 2) then
  217. ut_l = wvec_l(3)
  218. ut_r = wvec_r(3)
  219. ut1_l = 0.0D0
  220. ut1_r = 0.0D0
  221. elseif(ndim .eq. 3) then
  222. ut_l = wvec_l(3)
  223. ut_r = wvec_r(3)
  224. ut1_l = wvec_l(4)
  225. ut1_r = wvec_r(4)
  226. endif
  227. C
  228. C We should take into account the formation energy in the flux
  229. C
  230. call fauupt(nesp1, nsca1, ndim,
  231. & gam_l, rho_l, p_l, un_l, ut_l, ut1_l, (et_l + hf_l),
  232. & gam_r, rho_r, p_r, un_r, ut_r, ut1_r, (et_r + hf_r),
  233. C & yl,yr,sl,sr,
  234. & 1.0D0, 1.0D0, 1.0D0, 1.0D0,
  235. & v_inf, flu,
  236. & cellt)
  237. C
  238. flurho = flu(1)
  239. flu(3+ndim) = 0.0d0
  240. C The Larrouturou theorem does hold
  241. coefy = 0.5d0 * (1.0d0 + sign(1.0d0, flurho))
  242. omcoey = 1.0d0 - coefy
  243. C
  244. do icomp = (ndim + 4), (ndim + nesp + 4), 1
  245. flu(icomp) = flurho * ((coefy * wvec_l(icomp)) +
  246. & (omcoey * wvec_r(icomp)))
  247. enddo
  248. C
  249. C******* f_lag is here not defined !!!
  250. C
  251. 9999 continue
  252. return
  253. end
  254.  
  255.  
  256.  

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