Télécharger fvlhre.eso

Retour à la liste

Numérotation des lignes :

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

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