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