1 : C FCOLRE SOURCE BECC 09/12/07 21:15:17 6579
2 : subroutine fcolre(ndim, nesp, nLHS, nordpo,
3 : & reacoe, aprop, Runiv,
4 : & Tmaxcv,
5 : & wvec_l,
6 : & wvec_r,
7 : & wwork, y, acv_l, acv_r, acv_z, acv_w ,
8 : & flu, f_lag, cellt,
9 : & logan)
10 : C
11 : c************************************************************************
12 : c
13 : c project :
14 : c
15 : c name : fcolre
16 : c
17 : c description : Reactive Euler equations for a mixture of reactive
18 : C thermally perfect gases.
19 : c
20 : c Shock-shock fds + entropy fix, implemented in a
21 : C "Discrete Equation Method" fashion.
22 : c
23 : c Voir:
24 : c 1) Beccantini, report ???
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
34 : c called sub : prithe, racrea, flufun
35 : C
35 : C called by : comres
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 wvec_l = left primitive variables
62 : c
63 : c wvec_r = right primitive variables
64 : C
65 : C NB Primitive variables =
66 : C
67 : C w(1) = rho
68 : C w(2) = ux
69 : C w(3) = uy
70 : C w(4) = uz
71 : C w(2+ndim) = p
72 : C w(3+ndim) = csi
73 : C w(3+ndim+1) = yf(1)
74 : C ...
75 : C w(3+ndim+nLHS) = yf(nLHS)
76 : C w(3+ndim+nLHS+1) = yi(nLHS+1)
77 : C ...
78 : C w(3+ndim+nesp-1) = yi(nesp-1)
79 : C w(3+ndim+nesp) = dy1t = yi(1) - yf(1)
80 : C w(4+ndim+nesp) = k0
81 : C
82 : C For idim = 1, the total number of variables
83 : C is 4 + 1 + nesp = 5 + nesp.
84 : c
85 : c wwork, y, acv_l,
86 : C acv_r, acv_z,
87 : C acv_w = temporary work vectors
88 : C wwork(4+ndim+nesp), y(1:nesp),
89 : C acv_l(0:nordpo),
90 : C acv_r(0:nordpo),
91 : C acv_z(0:nordpo),
92 : C acv_w (0:nordpo)
93 : C
94 : c**** output
95 : c
96 : c flu = Eulerian interfacial flux in (n,t1,t2), i.e.
97 : c rho*un mass
98 : C rho*un*un + p momentum along n
99 : C rho*un*ut1 momentum along t1
100 : * rho*un*ut2 momentum along t2
101 : C rho*un*ht energy
102 : c 0.0 alpha
103 : C rho*un*y_{1,f} mass * y_{1,f}
104 : C ...
105 : C rho*un*y_{nLHS+1,i} mass * y_{nLHS+1,i}
106 : c ...
107 : C rho*un*dy_{1,t} mass * dy_{1,t}
108 : C rho*un*k0 mass * k0
109 : C
110 : C NB. flu(1:(4+ndim+nesp))
111 : c
112 : C f_lag = Lagrangian interfacial flux on the reactive
113 : C discontinuty in (n,t1,t2), i.e.
114 : c rho*(un - D_resh) mass
115 : C rho*un*un + p - (D_resh rho un) momentum along n
116 : C rho*(un - D_resh)*ut1 momentum along t1
117 : C rho*(un - D_resh)*ut2 momentum along t2
118 : C rho*un*ht - (D_resh rho et) energy
119 : c - D_resh alpha
120 : c rho*(un - D_resh)*y_{1,f} mass * y_{1,f}
121 : C ...
122 : C
123 : C NB. f_lag(1:(4+ndim+nesp))
124 : C
125 : c
126 : c cellt = stability condition, i.e.
127 : c
128 : c dt/dx < cellt (dimension 1/velocity)
129 : c
130 : c logan = si true, anomaly detected
131 : c
132 : c************************************************************************
133 : c
134 : C 07/12/2009 created
135 : c
136 : c************************************************************************
137 : c
138 : c n.b.: all variables are declared
139 : c
140 : C implicit none
141 : integer ndim, nesp, nLHS, nordpo
142 : integer icomp
143 : real*8 zero, rind_l
144 : parameter (zero=1.0d-10)
145 : real*8 reacoe(nesp), aprop(0:nordpo+2,nesp), Runiv, Tmaxcv
146 : & , xst, wvec_l(4+ndim+nesp), wvec_r(4+ndim+nesp)
147 : & , flu(4+ndim+nesp)
148 : & , f_lag(4+ndim+nesp), cellt
149 : real*8 wwork(4+ndim+nesp)
150 : & , acv_l(0:nordpo), acv_r(0:nordpo), acv_z(0:nordpo)
151 : & , acv_w (0:nordpo)
152 : $ , y(1:nesp)
153 : real*8 csi0, csi1, q
154 : C Left
155 : real*8 gam_l, Rgas_l, cv_l, hf_l, rho_l, u_l(1:3), p_l, T_l
156 : $ , et_l, csi_l, un_l
157 : C Left star... same Rgas, acv, hf and csi as left
158 : real*8 rho_ls, p_ls, T_ls, d_ls, un_ls
159 : C Right star star
160 : real*8 ga_rss, Rg_rss, cv_rss, hf_rss, cs_rss, rh_rss
161 : $ , p_rss, et_rss, T_rss, un_rss
162 : C Right star... same Rgas, acv, hf and csi as right
163 : real*8 rho_rs, p_rs, d_rs, T_rs, un_rs
164 : C Right
165 : real*8 gam_r, Rgas_r, cv_r, hf_r, rho_r, u_r(1:3), p_r, T_R
166 : & , et_r
167 : & , csi_r, un_r
168 : real*8 c_l, c_r, c_ls, c_rs
169 : real*8 k0, D_resh
170 : real*8 vit_l, vit_ls, vit_rs, vit_r
171 : real*8 coefy, omcoey
172 : real*8 flurho, fluru, fluret, trekin
173 : C
174 : logical logan, lognc, logdeb
175 : parameter (logdeb = .false.)
176 : C
177 : C We have to distinguish between 5 different states
178 : C
179 : C l = burnt state (csi = csimax)
180 : C | Left GNL wave
181 : C ls = connected to l by a GNL wave (csi = csimax)
182 : C | Contact discontinuty
183 : C rsw = connected to rs by a reactive wave (csi = csimax)
184 : C | Reactive shock (weak deflagration wave)
185 : C rs = connected to r by a GNL wave (csi = csimin)
186 : C | Right GNL wave
187 : C r = unburnt state (csi = csimin)
188 : C
189 : C We define the left state such that csi_l > csi_r
190 : C
191 : csi0 = wvec_l(ndim+3)
192 : csi1 = wvec_r(ndim+3)
193 : if (csi0 .gt. csi1) then
194 : k0 = wvec_r(4+ndim+nesp)
195 : rind_l = 1.0D0
196 : C
197 : C The right star star state is the burnt state relatively to the
198 : C right and right stare states. It has the same value of csi as
199 : C the left state.
200 : C Note that if this state has not necessary the same mass
201 : C fractions as the left state.
202 : C We compute the physical properties of the right star state,
203 : C namely
204 : C ga_rss, Rg_rss, cv_rss, acv_z, hf_rss.
205 : C The other physical values are of course meaningless!
206 : C
207 : do icomp = 1, (nesp + 4 + ndim), 1
208 : wwork(icomp) = wvec_r(icomp)
209 : enddo
210 : wwork(3+ndim) = wvec_l(3+ndim)
211 : C
211 : call prithe(ndim, nesp, nLHS, nordpo, aprop, reacoe, Runiv,
212 : $ Tmaxcv,wwork,ga_rss, Rg_rss, cv_rss, acv_z, hf_rss,
213 : $ rh_rss,u_r, p_rss, T_rss ,y, et_rss, cs_rss, logan)
214 : if (logan) then
215 : write(*,'(/A8)') 'fcolre.f'
216 : write(*,*) 'ANOMALY DETECTED 3'
217 : goto 9999
218 : endif
219 : C
220 : C The left state is considered as the burnt state.
221 : C We compute the physical properties of the left state, namely
222 : C gam_l, Rgas_l, cv_l, acv_l, hf_l.
223 : C We also compute the physical variable on the left, namely
224 : C rho_l, u_l, p_l, T_L, et_l and also csi_r
225 : C
225 : call prithe(ndim, nesp, nLHS, nordpo, aprop, reacoe, Runiv,
226 : $ Tmaxcv,wvec_l,gam_l, Rgas_l, cv_l, acv_l, hf_l, rho_l,
227 : $ u_l, p_l, T_l,y, et_l, csi_l, logan)
228 : if (logan) then
229 : write(*,'(/A8)') 'fcolre.f'
230 : write(*,*) 'ANOMALY DETECTED 1'
231 : goto 9999
232 : endif
233 : C
234 : C The right state is considered as the unburnt state.
235 : C We compute the physical properties of the right state, namely
236 : C gam_r, Rgas_r, cv_r, acv_r, hf_r.
237 : C We also compute the physical variable on the left, namely
238 : C rho_r, u_r, p_r, T_r, et_r, csi_r
239 : C
239 : call prithe(ndim, nesp, nLHS, nordpo, aprop, reacoe, Runiv,
240 : $ Tmaxcv,wvec_r,gam_r, Rgas_r, cv_r, acv_r, hf_r, rho_r,
241 : $ u_r, p_r, T_r,y, et_r, csi_r, logan)
242 : if (logan) then
243 : write(*,'(/A8)') 'fcolre.f'
244 : write(*,*) 'ANOMALY DETECTED 2'
245 : goto 9999
246 : endif
247 : else
248 : rind_l = -1.0D0
249 : k0 = wvec_l(4+ndim+nesp)
250 : C
251 : do icomp = 1, (nesp + 4 + ndim), 1
252 : wwork(icomp)=wvec_l(icomp)
253 : enddo
254 : wwork(3+ndim)=wvec_r(3+ndim)
255 : C
255 : call prithe(ndim, nesp, nLHS, nordpo, aprop, reacoe, Runiv,
256 : $ Tmaxcv,wwork,ga_rss, Rg_rss, cv_rss, acv_z, hf_rss,
257 : $ rh_rss,u_l, p_rss, T_rss, y, et_rss, cs_rss, logan)
258 : if (logan) then
259 : write(*,'(/A8)') 'fcolre.f'
260 : write(*,*) 'ANOMALY DETECTED 6'
261 : goto 9999
262 : endif
263 : C
263 : call prithe(ndim, nesp, nLHS, nordpo, aprop, reacoe, Runiv,
264 : $ Tmaxcv,wvec_r,gam_l, Rgas_l, cv_l, acv_l, hf_l, rho_l,
265 : $ u_l, p_l, T_l,y, et_l, csi_l, logan)
266 : if (logan) then
267 : write(*,'(/A8)') 'fcolre.f'
268 : write(*,*) 'ANOMALY DETECTED 4'
269 : goto 9999
270 : endif
271 : C un_l=-1*u_l(1)
272 : C
272 : call prithe(ndim, nesp, nLHS, nordpo, aprop, reacoe, Runiv,
273 : $ Tmaxcv,wvec_l,gam_r, Rgas_r, cv_r, acv_r, hf_r, rho_r,
274 : $ u_r, p_r, T_r,y, et_r, csi_r, logan)
275 : if (logan) then
276 : write(*,'(/A8)') 'fcolre.f'
277 : write(*,*) 'ANOMALY DETECTED 5'
278 : goto 9999
279 : endif
280 : C un_r=-1*u_r(1)
281 : endif
282 : un_r = rind_l * u_r(1)
283 : un_l = rind_l * u_l(1)
284 : C
285 : C**** Test
286 : C
287 : if (logdeb ) then
288 : write(*,'(A23)') 'rho_r, p_r, t_r, u_r'
289 : write(*,'(4E16.6)') rho_r, p_r, t_r, un_r
290 : write(*,'(A23)') 'rho_l, p_l, t_l, u_l'
291 : write(*,'(4E16.6)') rho_l, p_l, t_l, un_l
292 : endif
293 : C
294 : C The reaction should be exothermic
295 : C
296 : q = hf_r - hf_rss
297 : if (q .lt. (-1.0D0 * zero * max( abs(hf_rss), abs(hf_r) )) ) then
298 : write(*,'(/A8)') 'fcolre.f'
299 : write(*,*) 'Endothermic reaction ??? '
300 : write(*,*) 'q =', q
301 : logan = .true.
302 : goto 9999
303 : endif
304 : q = 0.5d0 * (q + abs(q))
305 : C
306 : C Computation of the intermediate states
307 : C
308 : C Vacuum (TO DO!!!)
309 : C
310 : C Intermediate states in the case of no-vacuum
311 : C
312 : C write(*,*) k0, q
313 : C write(*,*) rho_l, p_l, un_l, T_l, et_l, gam_l
314 : C write(*,*) rho_r, p_r, un_r, T_r, et_r, gam_r
315 : C write(*,*) rh_rss, p_rss, un_l, T_rss, et_rss, ga_rss
316 : C write(*,*) csi0, csi1
317 : C
317 : call racrea(nordpo, Tmaxcv,
318 : & Rgas_l, acv_l, Rg_rss, acv_z, Rgas_r, acv_r,
319 : & q, k0,
320 : & rho_l, p_l, un_l, T_l, et_l, gam_l,
321 : & rho_r, p_r, un_r, T_r, et_r, gam_r,
322 : & rho_ls, p_ls, un_ls, T_ls, c_ls, d_ls,
323 : & rh_rss, p_rss, un_rss, T_rss, D_resh,
324 : & rho_rs, p_rs, un_rs, T_rs, c_rs, d_rs,
325 : & logan, lognc)
326 : if (lognc) then
327 : C
328 : C Maybe this is due to unphysical problems involved in this
329 : C approach
330 : C
331 : do icomp = 1, (4 + ndim + nesp), 1
332 : flu(icomp) = 0.0d0
333 : enddo
334 : C
335 : do icomp = 1, (4 + ndim + nesp), 1
336 : f_lag(icomp) = 0.0d0
337 : enddo
338 : goto 9999
339 : endif
340 : C
341 : if (logan) then
342 : write(*,'(/A8)') 'fcolre.f'
343 : write(*,*) 'rho_l, p_l, un_l, T_l, et_l, gam_l'
344 : write(*,'(6E16.6)') rho_l, p_l, un_l, T_l, et_l, gam_l
345 : write(*,*) 'rho_r, p_r, un_r, T_r, et_r, gam_r'
346 : write(*,'(6E16.6)') rho_r, p_r, un_r, T_r, et_r, gam_r
347 : write(*,*) ' icomp, wvec_l, wvec_r'
348 : do icomp = 1, (4 + nesp + ndim)
349 : write(*,'(I6,2E16.6)') icomp, wvec_l(icomp), wvec_r(icomp)
350 : enddo
351 : write(*,*) 'ANOMALY DETECTED'
352 : goto 9999
353 : endif
354 : C
355 : if (logdeb) then
356 : write(*,*) rho_l, p_l, un_l, T_l
357 : write(*,*) rho_ls, p_ls, un_ls, T_ls
358 : write(*,*) rh_rss, p_rss, un_rss, T_rss
359 : write(*,*) rho_rs, p_rs, un_rs, T_rs
360 : write(*,*) rho_r, p_r, un_r, T_r
361 : C stop
362 : endif
363 : C
364 : C**** Wave speed
365 : C
366 : C Speed of the GNL-waves (rarefaction or shock)
367 : C
368 : C
369 : C L LS
370 : C ****** *******
371 : C * or *
372 : C * LS L*
373 : C ******* *******
374 : C
375 : C On the rarefaction we lineairize over the variables
376 : C rho, T, w
377 : C
378 : c_l = gam_l * p_l / rho_l
379 : c_l = c_l ** 0.5d0
380 : c_r = gam_r * p_r / rho_r
381 : c_r = c_r ** 0.5d0
382 : C vit_l = min(d_ls, (un_l - c_l))
383 : vit_l = 0.5d0*((d_ls + (un_l - c_l))
384 : & - abs( d_ls - (un_l - c_l) ))
385 : C vit_ls = max(d_ls, (un_ls - c_ls))
386 : vit_ls = 0.5d0*((d_ls + (un_ls - c_ls))
387 : & + abs( d_ls - (un_ls - c_ls) ))
388 : C vit_ls exactly => vit_l
389 : vit_ls = max(vit_ls, vit_l)
390 : C vit_rs = min(d_rs, (u_rs + c_rs))
391 : vit_rs = 0.5d0*((d_rs + (un_rs + c_rs))
392 : & - abs( d_rs - (un_rs + c_rs) ))
393 : C vit_r = max(d_r, (un_r + c_r))
394 : vit_r = 0.5d0*((d_rs + (un_r + c_r))
395 : & + abs( d_rs - (un_r + c_r) ))
396 : C vit_rs exactly =< vit_r
397 : vit_rs = min(vit_rs, vit_r)
398 : C
399 : C write(*,*) vit_l, vit_ls, un_ls, un_rss, D_resh, vit_rs, vit_r
400 : C stop
401 : C
402 : C**** Maximum eigenvalue
403 : C
404 : cellt = abs(vit_r)
405 : cellt = 0.5d0 * (abs(vit_l) + cellt + abs ( abs(vit_l) - cellt))
406 : C
407 : C**** Eulerian flux (1D)
408 : C
409 : xst = 0.0d0
409 : call flufun(nordpo,
410 : & Tmaxcv, acv_w ,
411 : & xst,
412 : & vit_l, vit_ls, un_ls, D_resh, vit_rs, vit_r,
413 : & Rgas_l, acv_l, rho_l, un_l, T_l, hf_l,
414 : & rho_ls, T_ls,
415 : & Rg_rss, acv_z, rh_rss, un_rss, T_rss, hf_rss,
416 : & rho_rs, un_rs, T_rs,
417 : & Rgas_r, acv_r, rho_r, un_r, T_r, hf_r,
418 : & flurho, fluru, fluret)
419 : C In order to obtain the flux in the true left-right frame, I have
420 : C to multiply the speed * rind_l
421 : flurho = flurho * rind_l
422 : fluret = fluret * rind_l
423 : flu(1) = flurho
424 : flu(2) = fluru
425 : flu(2+ndim) = fluret
426 : flu(3+ndim) = 0.0d0
427 : C Note that the Larrouturou theorem does not hold
428 : C Namely un_ls (= un_rss) has not the same sign as flurho
429 : C write (*,*) un_ls, un_rss, flurho
430 : C un_ls > 0 => coefy = 1
431 : coefy = 0.5d0 * (1.0d0 + sign(1.0d0, (rind_l * un_ls)))
432 : omcoey = 1.0d0 - coefy
433 : C Tangential speeds (ndim .ge. 2)
434 : do icomp = 3, (ndim + 1), 1
435 : flu(icomp) = flurho * ((coefy * wvec_l(icomp)) +
436 : & (omcoey * wvec_r(icomp)))
437 : enddo
438 : do icomp = (ndim + 4), (ndim + nesp + 4), 1
439 : flu(icomp) = flurho * ((coefy * wvec_l(icomp)) +
440 : & (omcoey * wvec_r(icomp)))
441 : enddo
442 : C Kinetic energy flux
443 : trekin = 0.0D0
444 : do icomp = 2, ndim, 1
445 : trekin = trekin + (((coefy * wvec_l(icomp + 1)) +
446 : & (omcoey * wvec_r(icomp + 1))) ** 2)
447 : enddo
448 : trekin = 0.5D0 * trekin
449 : flu(2+ndim) = flu(2+ndim) + (flurho * trekin)
450 : C
451 : C**** Lagrangian flux over the reactive shock
452 : C
453 : xst = D_resh
453 : call flufun(nordpo,
454 : & Tmaxcv, acv_w ,
455 : & xst,
456 : & vit_l, vit_ls, un_ls, D_resh, vit_rs, vit_r,
457 : & Rgas_l, acv_l, rho_l, un_l, T_l, hf_l,
458 : & rho_ls, T_ls,
459 : & Rg_rss, acv_z, rh_rss, un_rss, T_rss, hf_rss,
460 : & rho_rs, un_rs, T_rs,
461 : & Rgas_r, acv_r, rho_r, un_r, T_r, hf_r,
462 : & flurho, fluru, fluret)
463 : C
464 : flurho = flurho * rind_l
465 : fluret = fluret * rind_l
466 : f_lag(1) = flurho
467 : f_lag(2) = fluru
468 : f_lag(2+ndim) = fluret
469 : f_lag(3+ndim) = -1.0D0 * rind_l * D_resh
470 : C
471 : C Since we suppose that left state is the burnt one,
472 : C D_resh > un_ls, i.e. over (x/t) = D_rest we have the
473 : C right passive scalar !
474 : C indl = 1 => burnt state in wvec_l => passive scalar in r
475 : C indl = -1 => burnt state in wvec_r => passive scalar in l
476 : C Note that the result does not change as q-> 0
477 : C (virtual combustion wave).
478 : C
479 : coefy = 0.5d0 * (1.0d0 + sign(1.0d0, (1.0D0 * rind_l)))
480 : omcoey = 1.0d0 - coefy
481 : C
482 : do icomp = 3, (ndim + 1), 1
483 : f_lag(icomp) = flurho * ((coefy * wvec_r(icomp)) +
484 : & (omcoey * wvec_l(icomp)))
485 : enddo
486 : do icomp = (ndim + 4), (ndim + nesp + 4), 1
487 : f_lag(icomp) = flurho * ((coefy * wvec_r(icomp)) +
488 : & (omcoey * wvec_l(icomp)))
489 : enddo
490 : C Kinetic energy flux
491 : trekin = 0.0D0
492 : do icomp = 2, ndim, 1
493 : trekin = trekin + (((coefy * wvec_r(icomp + 1)) +
494 : & (omcoey * wvec_l(icomp + 1))) ** 2)
495 : enddo
496 : trekin = 0.5D0 * trekin
497 : C write(*,*) trekin
498 : f_lag(2+ndim) = f_lag(2+ndim) + (flurho * trekin)
499 : C
500 : 9999 continue
501 : return
502 : end
503 :
© Cast3M 2003 - All rights reserved.
Disclaimer