Télécharger fbecre.eso

Retour à la liste

Numérotation des lignes :

  1. C FBECRE SOURCE BECC 09/12/07 21:15:16 6579
  2. subroutine fbecre(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 project :
  15. c
  16. c name : fbecre
  17. C reaction shock in WDF, CJDF, CJDT, SDT
  18. c
  19. c description : Reactive Euler equations for a mixture of reactive
  20. C thermally perfect gases.
  21. c
  22. c Shock-shock fds + entropy fix, implemented in a
  23. C "Discrete Equation Method" fashion.
  24. c
  25. c Voir:
  26. c 1) Beccantini, report ???
  27. C 2) Metayer, Massoni, Saurel, "Modelling
  28. C evaporations fronts with reactive Riemann
  29. C solvers.", JCP 205, 2005
  30. c
  31. c language : fortran 77
  32. c
  33. c author : a. beccantini den/dm2s/sfme/ltmf
  34. c
  35. c************************************************************************
  36. c
  37. c called sub : prithe, racre2, flufun
  38. C
  39. C called by : comres
  40. c
  41. c************************************************************************
  42. c
  43. c**** input
  44. C
  45. C ndim = dimension of the domain
  46. c
  47. c nesp = species involved in the Euler equations
  48. C
  49. C nLHS = species involved in the LHS of the global
  50. C reaction
  51. c
  52. c nordp0 = polynomial order of the cv
  53. c
  54. c reacoe = coefficient involved in the global reaction
  55. c (positive in the LHS, negative in the RHS...
  56. c the first one should be positive)
  57. c
  58. c aprop = gas properties (cv, molar mass, formation
  59. C energy)
  60. c
  61. C Runiv = universal gas constant
  62. c
  63. C Tmaxcv = Tmax for cv expansion
  64. C
  65. C cnna, ctna, ct1na = scalar products
  66. C (\vec{n},\vec{n}_{\rm alpha})
  67. C (\vec{t},\vec{n}_{\rm alpha})
  68. C (\vec{t1},\vec{n}_{\rm alpha})
  69. C HP. We suppose that \vec{n}_{\rm alpha} is
  70. C always from the burnt to the unburnt state
  71. C
  72. c wvec_l = left primitive variables
  73. c
  74. c wvec_r = right primitive variables
  75. C
  76. C NB Primitive variables =
  77. C
  78. C w(1) = rho
  79. C w(2) = ux
  80. C w(3) = uy
  81. C w(4) = uz
  82. C w(2+ndim) = p
  83. C w(3+ndim) = csi
  84. C w(3+ndim+1) = yf(1)
  85. C ...
  86. C w(3+ndim+nLHS) = yf(nLHS)
  87. C w(3+ndim+nLHS+1) = yi(nLHS+1)
  88. C ...
  89. C w(3+ndim+nesp-1) = yi(nesp-1)
  90. C w(3+ndim+nesp) = dy1t = yi(1) - yf(1)
  91. C w(4+ndim+nesp) = k0
  92. C
  93. C For ndim = 1, the total number of variables
  94. C is 4 + 1 + nesp = 5 + nesp.
  95. c
  96. c wwork, y, acv_l,
  97. C acv_r, acv_b,
  98. C acv_w = temporary work vectors
  99. C wwork(4+ndim+nesp), y(1:nesp),
  100. C acv_l(0:nordpo),
  101. C acv_r(0:nordpo),
  102. C acv_b(0:nordpo),
  103. C acv_w (0:nordpo)
  104. C
  105. c**** output
  106. c
  107. c flu = Eulerian interfacial flux in (n,t1,t2), i.e.
  108. c rho*un mass
  109. C rho*un*un + p momentum along n
  110. C rho*un*ut1 momentum along t1
  111. * rho*un*ut2 momentum along t2
  112. C rho*un*ht energy
  113. c 0.0 alpha
  114. C rho*un*y_{1,f} mass * y_{1,f}
  115. C ...
  116. C rho*un*y_{nLHS+1,i} mass * y_{nLHS+1,i}
  117. c ...
  118. C rho*un*dy_{1,t} mass * dy_{1,t}
  119. C rho*un*k0 mass * k0
  120. C
  121. C NB. flu(1:(4+ndim+nesp))
  122. c
  123. C f_lag = Lagrangian interfacial flux on the reactive
  124. C discontinuty in (n,t1,t2), i.e.
  125. c rho*(un - D_resh) mass
  126. C rho*un*un + p - (D_resh rho un) momentum along n
  127. C rho*(un - D_resh)*ut1 momentum along t1
  128. C rho*(un - D_resh)*ut2 momentum along t2
  129. C rho*un*ht - (D_resh rho et) energy
  130. c - D_resh alpha
  131. c rho*(un - D_resh)*y_{1,f} mass * y_{1,f}
  132. C ...
  133. C
  134. C NB. f_lag(1:(4+ndim+nesp))
  135. C
  136. c
  137. c cellt = stability condition, i.e.
  138. c
  139. c dt/dx < cellt (dimension 1/velocity)
  140. c
  141. c logan = si true, anomaly detected
  142. c
  143. c************************************************************************
  144. c
  145. C 07/12/2009 created
  146. C
  147. c************************************************************************
  148. c
  149. c n.b.: all variables are declared
  150. c
  151. c implicit none
  152. integer ndim, nesp, nLHS, nordpo, iter
  153. integer icomp
  154. real*8 zero
  155. parameter (zero=1.0d-10)
  156. real*8 rind_l
  157. real*8 reacoe(nesp), aprop(0:(nordpo+2),nesp), Runiv, Tmaxcv
  158. & , wvec_l(4+ndim+nesp), wvec_r(4+ndim+nesp)
  159. & , flu(4+ndim+nesp)
  160. & , f_lag(4+ndim+nesp), cellt
  161. real*8 wwork(4+ndim+nesp)
  162. & , acv_l(0:nordpo), acv_r(0:nordpo), acv_b(0:nordpo)
  163. & , acv_w(0:nordpo)
  164. $ , y(1:nesp)
  165. real*8 csi0, csi1, k0
  166. real*8 scnna, cnna, ctna, ct1na
  167. real*8 gam_b, Rgas_b, cv_b, hf_b, rho_b, u_b(1:3), un_b, p_b, T_b
  168. & , et_b, csi_b
  169. & , gam_l, Rgas_l, cv_l, hf_l, rho_l, u_l(1:3), p_l, T_l
  170. & , et_l, csi_l
  171. & , gam_r, Rgas_r, cv_r, hf_r, rho_r, u_r(1:3), p_r, T_r
  172. & , et_r, csi_r
  173. & , un_r, un_l, q
  174. real*8 u_l1, u_r1
  175. real*8 d_l, d_r
  176. real*8 rho_d, p_d, un_d, T_d, d_d
  177. & , rho_c, p_c, un_c, T_c, d_c
  178. & , d_b
  179. & , rho_a, p_a, un_a, T_a, d_a
  180. real*8 xst, flurho, fluru, fluret, coefy, omcoey
  181. & , trekin, D_resh, k0new
  182. real*8 coefd, d_o
  183. & , k0m
  184. & , d_lm, d_rm
  185. & , rho_dm, p_dm, un_dm, T_dm, d_dm
  186. & , rho_cm, p_cm, un_cm, T_cm, d_cm
  187. & , rho_bm, p_bm, un_bm, T_bm, d_bm
  188. & , rho_am, p_am, un_am, T_am, d_am
  189. & , k0n
  190. & , d_ln, d_rn
  191. & , rho_dn, p_dn, un_dn, T_dn, d_dn
  192. & , rho_cn, p_cn, un_cn, T_cn, d_cn
  193. & , rho_bn, p_bn, un_bn, T_bn, d_bn
  194. & , rho_an, p_an, un_an, T_an, d_an
  195. & , k0o
  196. & , d_lo, d_ro
  197. & , rho_do, p_do, un_do, T_do, d_do
  198. & , rho_co, p_co, un_co, T_co, d_co
  199. & , rho_bo, p_bo, un_bo, T_bo, d_bo
  200. & , rho_ao, p_ao, un_ao, T_ao, d_ao
  201. C
  202. logical logan
  203. logical lognc
  204. logical logdeb, logk0n
  205. C logdeb = debugging ?
  206. C logk0n = .true. : instead of k0, we take k0 * abs(cnna)
  207. C .false. : instead of k0, we take k0 * abs(cnna)
  208. parameter (logdeb = .false., logk0n = .true.)
  209. C
  210. C HP $\vec{n}_{\rm \alpha}$ is always from the burnt to the unburnt
  211. C region
  212. C
  213. scnna = sign(1.0D0, cnna)
  214. C
  215. C In the reactive case, we have to distingish between 6 different
  216. C states
  217. C
  218. C l = burnt state (csi = csimax)
  219. C | Left GNL wave
  220. C ls = connected to l by a GNL wave (csi = csimax)
  221. C = d
  222. C | Contact discontinuty
  223. C rsss = connected to rss by a rarefaction wave (true
  224. C in CJDF and CJDT regime, = rss in WDF and SDT)
  225. C (csi = csimax)
  226. C = c
  227. C | Rarefaction wave
  228. C rss = connected to rs by a reactive wave (csi = csimax)
  229. C = b
  230. C | Reactive shock (WDF, CJDF, CJDT, SDT)
  231. C rs = connected to r by a GNL wave (csi = csimin)
  232. C in SDT, state behind the non-reactive shock travelling
  233. C at the same speed than the reactive shock
  234. C = a
  235. C | Right GNL wave
  236. C r = unburnt state (csi = csimin)
  237. C
  238. C
  239. C In the non-reactive case, we have to distingish between 4 different
  240. C states:
  241. C
  242. C l = burnt/unburnt state
  243. C | Left GNL wave
  244. C ls = connected to l by a GNL wave
  245. C = b
  246. C | Contact discontinuty
  247. C rs = connected to r by a GNL wave
  248. C = a
  249. C | Right GNL wave
  250. C r = burnt/unburnt state
  251. C
  252. csi0 = wvec_l(ndim+3)
  253. csi1 = wvec_r(ndim+3)
  254. C
  255. C csi0 = 0 or 1
  256. C csi1 = 0 or 1
  257. C
  258. if (abs(csi0 - csi1) .lt. 1.0D-1)then
  259. C if (.false.) then
  260. C
  261. C************************************************************************
  262. C******* NON REACTIVE CASE **********************************************
  263. C************************************************************************
  264. C
  265. call prithe(ndim, nesp, nLHS, nordpo, aprop, reacoe, Runiv,
  266. $ Tmaxcv,wvec_l,gam_l, Rgas_l, cv_l, acv_l, hf_l, rho_l,
  267. $ u_l, p_l, T_l,y, et_l, csi_l, logan)
  268. if (logan) then
  269. write(*,'(/A8)') 'fbecre.f'
  270. write(*,*) 'ANOMALY DETECTED 01'
  271. goto 9999
  272. endif
  273. C
  274. call prithe(ndim, nesp, nLHS, nordpo, aprop, reacoe, Runiv,
  275. $ Tmaxcv,wvec_r,gam_r, Rgas_r, cv_r, acv_r, hf_r, rho_r,
  276. $ u_r, p_r, T_r,y, et_r, csi_r, logan)
  277. if (logan) then
  278. write(*,'(/A8)') 'fcolre.f'
  279. write(*,*) 'ANOMALY DETECTED 02'
  280. goto 9999
  281. endif
  282. C
  283. un_l = u_l(1)
  284. un_r = u_r(1)
  285. C
  286. call racnre(nordpo, Tmaxcv,
  287. & Rgas_l, acv_l, Rgas_r, acv_r,
  288. & rho_l, p_l, un_l, T_l, gam_l,
  289. & rho_r, p_r, un_r, T_r, gam_r,
  290. & d_l, d_r,
  291. & rho_b, p_b, un_b, T_b, d_b,
  292. & rho_a, p_a, un_a, T_a, d_a,
  293. & logan, lognc)
  294. C
  295. cellt = max (abs(d_l), abs(d_r))
  296. C
  297. C******* States
  298. C
  299. C l
  300. C ***
  301. C |* r
  302. C | * b *********
  303. C | ******* *|
  304. C | | * a * |
  305. C | | ************* |
  306. C | | | | |
  307. C | | | | |
  308. C d_l d_b un_b=un_a d_a d_r
  309. C
  310. C
  311. xst = 0.0D0
  312. call flufu2(nordpo,
  313. & Tmaxcv, Rgas_r, acv_r,
  314. & Rgas_l, acv_l, acv_w,
  315. & xst,
  316. & d_l, d_b, d_a, d_r,
  317. & rho_l, T_l, un_l, hf_l,
  318. & rho_b, T_b, un_b,
  319. & rho_a, T_a, un_a,
  320. & rho_r, T_r, un_r, hf_r,
  321. & flurho, fluru, fluret)
  322. C
  323. flu(1) = flurho
  324. flu(2) = fluru
  325. flu(2+ndim) = fluret
  326. flu(3+ndim) = 0.0d0
  327. C The Larrouturou theorem does hold
  328. C Namely un_a (= un_b) has the same sign as flurho
  329. C write (*,*) un_d, un_c, flurho
  330. C un_ls > 0 => coefy = 1
  331. coefy = 0.5d0 * (1.0d0 + sign(1.0d0, flurho))
  332. omcoey = 1.0d0 - coefy
  333. C Tangential speeds (ndim .ge. 2)
  334. do icomp = 3, (ndim + 1), 1
  335. flu(icomp) = flurho * ((coefy * wvec_l(icomp)) +
  336. & (omcoey * wvec_r(icomp)))
  337. enddo
  338. do icomp = (ndim + 4), (ndim + nesp + 4), 1
  339. flu(icomp) = flurho * ((coefy * wvec_l(icomp)) +
  340. & (omcoey * wvec_r(icomp)))
  341. enddo
  342. C Kinetic energy flux
  343. trekin = 0.0D0
  344. do icomp = 2, ndim, 1
  345. trekin = trekin + (((coefy * wvec_l(icomp + 1)) +
  346. & (omcoey * wvec_r(icomp + 1))) ** 2)
  347. enddo
  348. trekin = 0.5D0 * trekin
  349. flu(2+ndim) = flu(2+ndim) + (flurho * trekin)
  350. C
  351. C******* f_lag is here not defined !!!
  352. C
  353. else
  354. C
  355. C************************************************************************
  356. C******* REACTIVE CASE **************************************************
  357. C************************************************************************
  358. C
  359. C We define the left state such that csi_l > csi_r
  360. C
  361. if (csi0 .gt. csi1) then
  362. k0 = wvec_r(4+ndim+nesp)
  363. rind_l = 1.0D0
  364. C
  365. C The rss (b) state is the burnt state relatively to the
  366. C right and right star (a) state. It has the same value of csi as
  367. C the left state.
  368. C Note that this state has not necessary the same mass
  369. C fractions as the left state.
  370. C We compute the physical properties of the rss (b) state
  371. C namely
  372. C gam_b, Rgas_b, cv_b, acv_b, hf_b.
  373. C The other physical values are of course meaningless!
  374. C
  375. do icomp = 1, (nesp + 4 + ndim), 1
  376. wwork(icomp) = wvec_r(icomp)
  377. enddo
  378. wwork(3+ndim) = wvec_l(3+ndim)
  379. C
  380. call prithe(ndim, nesp, nLHS, nordpo, aprop, reacoe, Runiv,
  381. $ Tmaxcv,wwork,gam_b, Rgas_b, cv_b, acv_b, hf_b,
  382. $ rho_b, u_b, p_b, T_b ,y, et_b, csi_b, logan)
  383. un_b = u_b(1)
  384. if (logan) then
  385. write(*,'(/A8)') 'fbecre.f'
  386. write(*,*) 'ANOMALY DETECTED 1'
  387. goto 9999
  388. endif
  389. C
  390. C The left state is considered as the burnt state.
  391. C We compute the physical properties of the left state, namely
  392. C gam_l, Rgas_l, cv_l, acv_l, hf_l.
  393. C We also compute the physical variable on the left, namely
  394. C rho_l, u_l, p_l, T_L, et_l and also csi_r
  395. C
  396. call prithe(ndim, nesp, nLHS, nordpo, aprop, reacoe, Runiv,
  397. $ Tmaxcv,wvec_l,gam_l, Rgas_l, cv_l, acv_l, hf_l, rho_l,
  398. $ u_l, p_l, T_l,y, et_l, csi_l, logan)
  399. if (logan) then
  400. write(*,'(/A8)') 'fbecre.f'
  401. write(*,*) 'ANOMALY DETECTED 2'
  402. goto 9999
  403. endif
  404. C
  405. C We compute the normal speed in the frame
  406. C (\vec{n}_{\rm alpha}, \vec{t}_{\rm alpha},
  407. C \vec{t1}_{\rm alpha})
  408. C
  409. u_l1 = u_l(1) * cnna * scnna
  410. if (ndim .ge. 2) then
  411. u_l1 = u_l1 + (u_l(2) * ctna * scnna)
  412. endif
  413. if (ndim .eq. 3) then
  414. u_l1 = u_l1 + (u_l(3) * ct1na * scnna)
  415. endif
  416. C
  417. C The right state is considered as the unburnt state.
  418. C We compute the physical properties of the right state, namely
  419. C gam_r, Rgas_r, cv_r, acv_r, hf_r.
  420. C We also compute the physical variable on the left, namely
  421. C rho_r, u_r, p_r, T_r, et_r, csi_r
  422. C
  423. call prithe(ndim, nesp, nLHS, nordpo, aprop, reacoe, Runiv,
  424. $ Tmaxcv,wvec_r,gam_r, Rgas_r, cv_r, acv_r, hf_r, rho_r,
  425. $ u_r, p_r, T_r,y, et_r, csi_r, logan)
  426. if (logan) then
  427. write(*,'(/A8)') 'fcolre.f'
  428. write(*,*) 'ANOMALY DETECTED 3'
  429. goto 9999
  430. endif
  431. u_r1 = u_r(1) * cnna * scnna
  432. if (ndim .ge. 2) then
  433. u_r1 = u_r1 + (u_r(2) * ctna * scnna)
  434. endif
  435. if (ndim .eq. 3) then
  436. u_r1 = u_r1 + (u_r(3) * ct1na * scnna)
  437. endif
  438. else
  439. rind_l = -1.0D0
  440. k0 = wvec_l(4+ndim+nesp)
  441. C
  442. do icomp = 1, (nesp + 4 + ndim), 1
  443. wwork(icomp)=wvec_l(icomp)
  444. enddo
  445. wwork(3+ndim)=wvec_r(3+ndim)
  446. C
  447. call prithe(ndim, nesp, nLHS, nordpo, aprop, reacoe, Runiv,
  448. $ Tmaxcv,wwork,gam_b, Rgas_b, cv_b, acv_b, hf_b,
  449. $ rho_b, u_b, p_b, T_b, y, et_b, csi_b, logan)
  450. un_b = u_b(1)
  451. if (logan) then
  452. write(*,'(/A8)') 'fcolre.f'
  453. write(*,*) 'ANOMALY DETECTED 4'
  454. goto 9999
  455. endif
  456. C
  457. call prithe(ndim, nesp, nLHS, nordpo, aprop, reacoe, Runiv,
  458. $ Tmaxcv,wvec_r,gam_l, Rgas_l, cv_l, acv_l, hf_l, rho_l,
  459. $ u_l, p_l, T_l,y, et_l, csi_l, logan)
  460. if (logan) then
  461. write(*,'(/A8)') 'fcolre.f'
  462. write(*,*) 'ANOMALY DETECTED 5'
  463. goto 9999
  464. endif
  465. u_l1 = u_l(1) * cnna * scnna
  466. if (ndim .ge. 2) then
  467. u_l1 = u_l1 + (u_l(2) * ctna * scnna)
  468. endif
  469. if (ndim .eq. 3) then
  470. u_l1 = u_l1 + (u_l(3) * ct1na * scnna)
  471. endif
  472. C
  473. call prithe(ndim, nesp, nLHS, nordpo, aprop, reacoe, Runiv,
  474. $ Tmaxcv,wvec_l,gam_r, Rgas_r, cv_r, acv_r, hf_r, rho_r,
  475. $ u_r, p_r, T_r,y, et_r, csi_r, logan)
  476. if (logan) then
  477. write(*,'(/A8)') 'fcolre.f'
  478. write(*,*) 'ANOMALY DETECTED 6'
  479. goto 9999
  480. endif
  481. u_r1 = u_r(1) * cnna * scnna
  482. if (ndim .ge. 2) then
  483. u_r1 = u_r1 + (u_r(2) * ctna * scnna)
  484. endif
  485. if (ndim .eq. 3) then
  486. u_r1 = u_r1 + (u_r(3) * ct1na * scnna)
  487. endif
  488. endif
  489. un_r = rind_l * u_r(1)
  490. un_l = rind_l * u_l(1)
  491. u_r1 = rind_l * u_r1
  492. u_l1 = rind_l * u_l1
  493. C
  494. C**** Test
  495. C
  496. if (logdeb) then
  497. write(*,'(A23)') 'rho_r, p_r, t_r, u_r, u_r1'
  498. write(*,'(4E16.6)') rho_r, p_r, t_r, un_r, u_r1
  499. write(*,'(A23)') 'rho_l, p_l, t_l, u_l, u_l1'
  500. write(*,'(4E16.6)') rho_l, p_l, t_l, un_l, u_l1
  501. C stop
  502. endif
  503. C
  504. C
  505. C The reaction should be exothermic
  506. C
  507. q = hf_r - hf_b
  508. if (q .lt. (-1.0D0 * zero * max( abs(hf_b), abs(hf_r) )) ) then
  509. write(*,'(/A8)') 'fcolre.f'
  510. write(*,*) 'Endothermic reaction ??? '
  511. write(*,*) 'q =', q
  512. logan = .true.
  513. goto 9999
  514. endif
  515. q = 0.5d0 * (q + abs(q))
  516. C
  517. C Computation of the speed of the wave
  518. C It can deal with vacuum
  519. C
  520. call racre2(nordpo, Tmaxcv,
  521. & Rgas_l, acv_l, Rgas_b, acv_b, Rgas_r, acv_r,
  522. & q, k0,
  523. & rho_l, p_l, u_l1, T_l, gam_l,
  524. & rho_r, p_r, u_r1, T_r, gam_r,
  525. & d_l, d_r,
  526. & rho_d, p_d, un_d, T_d, d_d,
  527. & rho_c, p_c, un_c, T_c, d_c,
  528. & rho_b, p_b, un_b, T_b, d_b,
  529. & rho_a, p_a, un_a, T_a, d_a,
  530. & logan, lognc)
  531. C
  532. if (ndim .eq. 1)then
  533. C
  534. C***************************
  535. C********** 1D *************
  536. C***************************
  537. C
  538. cellt = max (abs(d_l), abs(d_r))
  539. D_resh = d_b
  540. elseif (logk0n) then
  541. C
  542. C******************************
  543. C********** K0 n \cdot n_a ****
  544. C******************************
  545. C
  546. C
  547. C Intermediate state
  548. C It can deal with vacuum
  549. C
  550. k0new = max ((abs ((d_b - un_a) * cnna)), (k0*1.0D-3))
  551. C k0new = k0
  552. call racre2(nordpo, Tmaxcv,
  553. & Rgas_l, acv_l, Rgas_b, acv_b, Rgas_r, acv_r,
  554. & q, k0new,
  555. & rho_l, p_l, un_l, T_l, gam_l,
  556. & rho_r, p_r, un_r, T_r, gam_r,
  557. & d_l, d_r,
  558. & rho_d, p_d, un_d, T_d, d_d,
  559. & rho_c, p_c, un_c, T_c, d_c,
  560. & rho_b, p_b, un_b, T_b, d_b,
  561. & rho_a, p_a, un_a, T_a, d_a,
  562. & logan, lognc)
  563. C
  564. D_resh = d_b
  565. C
  566. if (logan) then
  567. write(*,'(/A8)') 'fbecre.f'
  568. write(*,*) 'ANOMALY DETECTED 7'
  569. goto 9999
  570. endif
  571. C
  572. cellt = max (abs(d_l), abs(d_r))
  573. C
  574. else
  575. C
  576. C******************************
  577. C********** D n \cdot n_a *****
  578. C******************************
  579. C
  580. C
  581. C coefd = rind_l * cnna
  582. C write(*,*) 'd_b, k0 = ', d_b, k0
  583. coefd = max( abs (cnna), 1.0D-3)
  584. d_o = d_b * coefd
  585. C
  586. C********** Linear interpolation
  587. C
  588. k0m = 1.0D-3 * k0
  589. call racre2(nordpo, Tmaxcv,
  590. & Rgas_l, acv_l, Rgas_b, acv_b, Rgas_r, acv_r,
  591. & q, k0m,
  592. & rho_l, p_l, un_l, T_l, gam_l,
  593. & rho_r, p_r, un_r, T_r, gam_r,
  594. & d_lm, d_rm,
  595. & rho_dm, p_dm, un_dm, T_dm, d_dm,
  596. & rho_cm, p_cm, un_cm, T_cm, d_cm,
  597. & rho_bm, p_bm, un_bm, T_bm, d_bm,
  598. & rho_am, p_am, un_am, T_am, d_am,
  599. & logan, lognc)
  600. C
  601. if (logan) then
  602. write(*,'(/A8)') 'fbecre.f'
  603. write(*,*) 'ANOMALY DETECTED 8'
  604. goto 9999
  605. endif
  606. k0m = abs (d_bm - un_am)
  607. C write(*,*) 'd_bm, k0m = ', d_bm, k0m
  608. C
  609. k0n = 1.0D3 * k0
  610. call racre2(nordpo, Tmaxcv,
  611. & Rgas_l, acv_l, Rgas_b, acv_b, Rgas_r, acv_r,
  612. & q, k0n,
  613. & rho_l, p_l, un_l, T_l, gam_l,
  614. & rho_r, p_r, un_r, T_r, gam_r,
  615. & d_ln, d_rn,
  616. & rho_dn, p_dn, un_dn, T_dn, d_dn,
  617. & rho_cn, p_cn, un_cn, T_cn, d_cn,
  618. & rho_bn, p_bn, un_bn, T_bn, d_bn,
  619. & rho_an, p_an, un_an, T_an, d_an,
  620. & logan, lognc)
  621. C
  622. if (logan) then
  623. write(*,'(/A8)') 'fbecre.f'
  624. write(*,*) 'ANOMALY DETECTED 9'
  625. goto 9999
  626. endif
  627. k0n = abs (d_bn - un_an)
  628. C write(*,*) 'd_bn, k0n = ', d_bn, k0n
  629. C
  630. if (d_o .le. d_bm) then
  631. d_l = d_lm
  632. d_r = d_rm
  633. C
  634. rho_d = rho_dm
  635. p_d = p_dm
  636. un_d = un_dm
  637. T_d = T_dm
  638. d_d = d_dm
  639. C
  640. rho_c = rho_cm
  641. p_c = p_cm
  642. un_c = un_cm
  643. T_c = T_cm
  644. d_c = d_cm
  645. C
  646. rho_b = rho_bm
  647. p_b = p_bm
  648. un_b = un_bm
  649. T_b = T_bm
  650. d_b = d_bm
  651. C
  652. rho_a = rho_am
  653. p_a = p_am
  654. un_a = un_am
  655. T_a = T_am
  656. d_a = d_am
  657. C
  658. C write(*,*) 'Sono piccolino come il Masakino !'
  659. C write(*,*) 'coefd , k0 =', coefd, k0m
  660. C write(*,*) 'do =', d_o
  661. C write(*,*) 'dn =', d_bn
  662. C write(*,*) 'dm =', d_bm
  663. C write(*,*) 'cnna =', cnna
  664. C
  665. C k0o = 1.0D-3 * k0
  666. C do iter = 1, 10
  667. C call racre2(nordpo, Tmaxcv,
  668. C & Rgas_l, acv_l, Rgas_b, acv_b, Rgas_r, acv_r,
  669. C & q, k0o,
  670. C & rho_l, p_l, un_l, T_l, gam_l,
  671. C & rho_r, p_r, un_r, T_r, gam_r,
  672. C & d_lo, d_ro,
  673. C & rho_do, p_do, un_do, T_do, d_do,
  674. C & rho_co, p_co, un_co, T_co, d_co,
  675. C & rho_bo, p_bo, un_bo, T_bo, d_bo,
  676. C & rho_ao, p_ao, un_ao, T_ao, d_ao,
  677. C & logan, lognc)
  678. C
  679. C if (logan) then
  680. C write(*,'(/A8)') 'fbecre.f'
  681. C write(*,*) 'ANOMALY DETECTED 10'
  682. C goto 9999
  683. C endif
  684. C
  685. C k0o = d_bo - un_ao
  686. C write(*,*) k0o
  687. C
  688. C write(*,*) 'd_b0, k0o = ', d_bo, k0o
  689. C k0o = k0o + 0.1D0 * (k0n - k0m)
  690. C enddo
  691. C
  692. C stop
  693. C
  694. elseif (d_o .ge. d_bn) then
  695. d_l = d_ln
  696. d_r = d_rn
  697. C
  698. rho_d = rho_dn
  699. p_d = p_dn
  700. un_d = un_dn
  701. T_d = T_dn
  702. d_d = d_dn
  703. C
  704. rho_c = rho_cn
  705. p_c = p_cn
  706. un_c = un_cn
  707. T_c = T_cn
  708. d_c = d_cn
  709. C
  710. rho_b = rho_bn
  711. p_b = p_bn
  712. un_b = un_bn
  713. T_b = T_bn
  714. d_b = d_bn
  715. C
  716. rho_a = rho_an
  717. p_a = p_an
  718. un_a = un_an
  719. T_a = T_an
  720. d_a = d_an
  721. C
  722. C write(*,*) 'Sono piu grande del Masakino !'
  723. C write(*,*) 'coefd , k0 =', coefd, k0n
  724. C write(*,*) 'do =', d_o
  725. C write(*,*) 'dn =', d_bn
  726. C write(*,*) 'dm =', d_bm
  727. C write(*,*) 'cnna =', cnna
  728. CC
  729. C k0o = 1.0D-3 * k0
  730. C do iter = 1, 11
  731. C call racre2(nordpo, Tmaxcv,
  732. C & Rgas_l, acv_l, Rgas_b, acv_b, Rgas_r, acv_r,
  733. C & q, k0o,
  734. C & rho_l, p_l, un_l, T_l, gam_l,
  735. C & rho_r, p_r, un_r, T_r, gam_r,
  736. C & d_lo, d_ro,
  737. C & rho_do, p_do, un_do, T_do, d_do,
  738. C & rho_co, p_co, un_co, T_co, d_co,
  739. C & rho_bo, p_bo, un_bo, T_bo, d_bo,
  740. C & rho_ao, p_ao, un_ao, T_ao, d_ao,
  741. C & logan, lognc)
  742. C
  743. C if (logan) then
  744. C write(*,'(/A8)') 'fbecre.f'
  745. C write(*,*) 'ANOMALY DETECTED 10'
  746. C goto 9999
  747. C endif
  748. C
  749. C k0o = d_bo - un_ao
  750. C write(*,*) k0o
  751. C
  752. C write(*,*) 'd_b0, k0o = ', d_bo, k0o
  753. C k0o = k0o + 0.1D0 * (k0n - k0m)
  754. C enddo
  755. CC
  756. C stop
  757. C
  758. else
  759. k0o = ((d_o - d_bm) * (k0n - k0m) / (d_bn - d_bm)) + k0m
  760. call racre2(nordpo, Tmaxcv,
  761. & Rgas_l, acv_l, Rgas_b, acv_b, Rgas_r, acv_r,
  762. & q, k0o,
  763. & rho_l, p_l, un_l, T_l, gam_l,
  764. & rho_r, p_r, un_r, T_r, gam_r,
  765. & d_lo, d_ro,
  766. & rho_do, p_do, un_do, T_do, d_do,
  767. & rho_co, p_co, un_co, T_co, d_co,
  768. & rho_bo, p_bo, un_bo, T_bo, d_bo,
  769. & rho_ao, p_ao, un_ao, T_ao, d_ao,
  770. & logan, lognc)
  771. C
  772. if (logan) then
  773. write(*,'(/A8)') 'fbecre.f'
  774. write(*,*) 'ANOMALY DETECTED 10'
  775. goto 9999
  776. endif
  777. C
  778. k0o = abs (d_bo - un_ao)
  779. C
  780. C write(*,*) 'd_bo, k0o = ', d_bo, k0o
  781. C
  782. do iter = 1, 10
  783. if (d_bo .le. d_o) then
  784. k0m = k0o
  785. d_bm = d_bm
  786. else
  787. k0n = k0o
  788. d_bn = d_bo
  789. endif
  790. C
  791. k0o = ((d_o - d_bm) * (k0n - k0m) / (d_bn - d_bm)) +
  792. $ k0m
  793. call racre2(nordpo, Tmaxcv,
  794. & Rgas_l, acv_l, Rgas_b, acv_b, Rgas_r, acv_r,
  795. & q, k0o,
  796. & rho_l, p_l, un_l, T_l, gam_l,
  797. & rho_r, p_r, un_r, T_r, gam_r,
  798. & d_lo, d_ro,
  799. & rho_do, p_do, un_do, T_do, d_do,
  800. & rho_co, p_co, un_co, T_co, d_co,
  801. & rho_bo, p_bo, un_bo, T_bo, d_bo,
  802. & rho_ao, p_ao, un_ao, T_ao, d_ao,
  803. & logan, lognc)
  804. C
  805. if (logan) then
  806. write(*,'(/A8)') 'fbecre.f'
  807. write(*,*) 'ANOMALY DETECTED 10'
  808. goto 9999
  809. endif
  810. C
  811. k0o = d_bo - un_ao
  812. C write(*,*) k0o
  813. C
  814. C write(*,*) 'd_b0, k0o = ', d_bo, k0o
  815. C
  816. d_l = d_lo
  817. d_r = d_ro
  818. C
  819. rho_d = rho_do
  820. p_d = p_do
  821. un_d = un_do
  822. T_d = T_do
  823. d_d = d_do
  824. C
  825. rho_c = rho_co
  826. p_c = p_co
  827. un_c = un_co
  828. T_c = T_co
  829. d_c = d_co
  830. C
  831. rho_b = rho_bo
  832. p_b = p_bo
  833. un_b = un_bo
  834. T_b = T_bo
  835. d_b = d_bo
  836. C
  837. rho_a = rho_ao
  838. p_a = p_ao
  839. un_a = un_ao
  840. T_a = T_ao
  841. d_a = d_ao
  842. C
  843. enddo
  844. C
  845. endif
  846. D_resh = d_b
  847. cellt = max (abs(d_l), abs(d_r))
  848. C
  849. endif
  850. C
  851. C******* States
  852. C
  853. C l
  854. C *** b
  855. C |* * r
  856. C | * d ** *********
  857. C | ******* * * *|
  858. C | | * c * * * |
  859. C | | ************ * a * |
  860. C | | | | ******* |
  861. C | | | | | | |
  862. C d_l d_d un_d=un_c d_c d_b d_a d_r
  863. C
  864. C
  865. C
  866. C******* Eulerian flux
  867. C
  868. xst = 0.0D0
  869. call flufu1(nordpo,
  870. & Tmaxcv, Rgas_r, acv_r, Rgas_b, acv_b,
  871. & Rgas_l, acv_l, acv_w,
  872. & xst,
  873. & d_l, d_d, d_c, d_b, d_a, d_r,
  874. & rho_l, T_l, un_l, hf_l,
  875. & rho_d, T_d, un_d,
  876. & rho_c, T_c, un_c,
  877. & rho_b, T_b, un_b, hf_b,
  878. & rho_a, T_a, un_a,
  879. & rho_r, T_r, un_r, hf_r,
  880. & flurho, fluru, fluret)
  881. C
  882. C In order to obtain the flux in the true left-right frame, I have
  883. C to multiply the speed * rind_l
  884. flurho = flurho * rind_l
  885. fluret = fluret * rind_l
  886. flu(1) = flurho
  887. flu(2) = fluru
  888. flu(2+ndim) = fluret
  889. flu(3+ndim) = 0.0d0
  890. C Note that the Larrouturou theorem does not hold
  891. C Namely un_d (= un_c) has not the same sign as flurho
  892. C write (*,*) un_d, un_c, flurho
  893. C un_ls > 0 => coefy = 1
  894. coefy = 0.5d0 * (1.0d0 + sign(1.0d0, (rind_l * un_d)))
  895. omcoey = 1.0d0 - coefy
  896. C Tangential speeds (ndim .ge. 2)
  897. do icomp = 3, (ndim + 1), 1
  898. flu(icomp) = flurho * ((coefy * wvec_l(icomp)) +
  899. & (omcoey * wvec_r(icomp)))
  900. enddo
  901. do icomp = (ndim + 4), (ndim + nesp + 4), 1
  902. flu(icomp) = flurho * ((coefy * wvec_l(icomp)) +
  903. & (omcoey * wvec_r(icomp)))
  904. enddo
  905. C Kinetic energy flux
  906. trekin = 0.0D0
  907. do icomp = 2, ndim, 1
  908. trekin = trekin + (((coefy * wvec_l(icomp + 1)) +
  909. & (omcoey * wvec_r(icomp + 1))) ** 2)
  910. enddo
  911. trekin = 0.5D0 * trekin
  912. flu(2+ndim) = flu(2+ndim) + (flurho * trekin)
  913. C
  914. C******* Lagrangian flux over the reactive shock
  915. C
  916. xst = D_resh
  917. call flufu1(nordpo,
  918. & Tmaxcv, Rgas_r, acv_r, Rgas_b, acv_b,
  919. & Rgas_l, acv_l, acv_w,
  920. & xst,
  921. & d_l, d_d, d_c, d_b, d_a, d_r,
  922. & rho_l, T_l, un_l, hf_l,
  923. & rho_d, T_d, un_d,
  924. & rho_c, T_c, un_c,
  925. & rho_b, T_b, un_b, hf_b,
  926. & rho_a, T_a, un_a,
  927. & rho_r, T_r, un_r, hf_r,
  928. & flurho, fluru, fluret)
  929. C
  930. flurho = flurho * rind_l
  931. fluret = fluret * rind_l
  932. f_lag(1) = flurho
  933. f_lag(2) = fluru
  934. f_lag(2+ndim) = fluret
  935. f_lag(3+ndim) = -1.0D0 * rind_l * D_resh
  936. C
  937. C Since we suppose that left state is the burnt one,
  938. C D_resh > un_d, i.e. over (x/t) = D_resh we have the
  939. C right passive scalar !
  940. C rind_l = 1 => burnt state in wvec_l => passive scalar in r
  941. C rind_l = -1 => burnt state in wvec_r => passive scalar in l
  942. C Note that the result does not change as q -> 0
  943. C (virtual combustion wave).
  944. C
  945. coefy = 0.5d0 * (1.0d0 + sign(1.0d0, (1.0D0 * rind_l)))
  946. omcoey = 1.0d0 - coefy
  947. C
  948. do icomp = 3, (ndim + 1), 1
  949. f_lag(icomp) = flurho * ((coefy * wvec_r(icomp)) +
  950. & (omcoey * wvec_l(icomp)))
  951. enddo
  952. do icomp = (ndim + 4), (ndim + nesp + 4), 1
  953. f_lag(icomp) = flurho * ((coefy * wvec_r(icomp)) +
  954. & (omcoey * wvec_l(icomp)))
  955. enddo
  956. C Kinetic energy flux
  957. trekin = 0.0D0
  958. do icomp = 2, ndim, 1
  959. trekin = trekin + (((coefy * wvec_r(icomp + 1)) +
  960. & (omcoey * wvec_l(icomp + 1))) ** 2)
  961. enddo
  962. trekin = 0.5D0 * trekin
  963. C write(*,*) trekin
  964. f_lag(2+ndim) = f_lag(2+ndim) + (flurho * trekin)
  965. endif
  966. C
  967. 9999 continue
  968. return
  969. end
  970.  
  971.  

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