1 : C RACREA SOURCE BECC 09/12/07 21:16:02 6579
2 : subroutine racrea(nordpo, Tmaxcv,
3 : & Rgas_l, acv_l, Rg_rss, acv_z, Rgas_r, acv_r,
4 : & q, k0,
5 : & rho_l, p_l, u_l, T_l, et_l, gam_l,
6 : & rho_r, p_r, u_r, T_r, et_r, gam_r,
7 : & rho_ls, p_ls, u_ls, T_ls, c_ls, d_ls,
8 : & rho_z, p_rss, u_rss, T_g_z, Dsw,
9 : & rho_rs, p_rs, u_rs, T_rs, c_rs, d_rs,
10 : & logan, lognc)
11 : C
12 : c************************************************************************
13 : c
14 : c project :
15 : c
16 : c name : racrea
17 : c
18 : c description : computation of the shock-shock field-by-field
19 : c decomposition of the reactive Riemann problem for
20 : c thermally perfect gases.
21 : C
22 : C We have to distinguish between 5 different states
23 : C
24 : C l = left state;
25 : C burnt state (csi = csimax)
26 : C |
27 : C | Left GNL wave
28 : C |
29 : C ls = left star state;
30 : C connected to l by a GNL wave
31 : C (csi = csimax)
32 : C |
33 : C | Contact discontinuity
34 : C |
35 : C rss = right star star state;
36 : C connected to rs by a reactive wave
37 : C (csi = csimax)
38 : C |
39 : C | Reactive shock (weak defl.)
40 : C |
41 : C rss = right star state;
42 : C connected to r by a GNL wave
43 : C (csi = csimin)
44 : C |
45 : C | Right GNL wave
46 : C |
47 : C r = right state;
48 : C unburnt state (csi = csimin)
49 : c
50 : c language : fortran 77
51 : c
52 : c author : a. beccantini den/dm2s/sfme/ltmf
53 : c
54 : c************************************************************************
55 : c
55 : c called by : fcolre
56 : c
56 : c called program : funco1, funrea
57 : c
58 : c************************************************************************
59 : c
60 : c
61 : c**** input
62 : c
63 : c nordpo = polynomials degree of cv(t)
64 : c tmaxcv = threshold temperature for cv(t)
65 : c rgas_l = gas constant (J/kg/K)in l (and in ls)
66 : c acv_l = cv_l = \sum_{i=0,nordpo} (acv_l{i} T^i)
67 : c rgas_rss = gas constant in rss
68 : c acv_z = cv_rss = \sum_{i=0,nordpo} (acv_z{i} T^i)
69 : c rgas_r = gas constant r (and in rs)
70 : c acv_r = cv_r = \sum_{i=0,nordpo} (acv_r{i} T^i)
71 : c
72 : c q = released heat per mass unit from rs to rss
73 : c k0 = fundamental flame speed
74 : c
75 : C rho_l = density in l
76 : C p_l = pressure in
77 : c u_l = velocity in l
78 : C T_l = temperature in l (redundant with rho_l and
79 : c p_l)
80 : c et_l = (sensible) internal energy in l (redundant
81 : c with rho_l and p_l)
82 : c gam_l = specific heat ratio in l (redundant
83 : c with rho_l and p_l)
84 : c
85 : C rho_r = density in r
86 : C p_r = pressure in r
87 : c u_r = velocity in r
88 : C T_r = temperature in r (redundant with rho_r and
89 : c p_r)
90 : c et_r = (sensible) internal energy in r (redundant
91 : c with rho_r and p_r)
92 : c gam_r = specific heat ratio in r (redundant
93 : c with rho_r and p_r)
94 : c
95 : c**** output
96 : c
97 : c
98 : C rho_ls = density in ls
99 : C p_ls = pressure in ls
100 : c u_ls = velocity in ls
101 : c T_ls = temperature in ls
102 : c c_ls = sound speed in ls
103 : c
104 : C rho_z = density in rss
105 : C p_rss = pressure in rss
106 : c u_rss = velocity in rss
107 : c T_rss = temperature in rss
108 : c Dsw = reactive shock velocity in the laboratory
109 : C frame
110 : C
111 : C rho_rs = density in rs
112 : C p_rs = pressure in rs
113 : c u_rs = velocity in rs
114 : c T_rs = temperature in rs
115 : c c_rs = sound speed in rs
116 : c
117 : c logan = an anomaly is detected
118 : c lognc = the apporach has some problems.
119 : c
120 : c************************************************************************
121 : c
122 : c 07/12/2009 Created.
123 : c
124 : c implicit none
125 : integer nmax, nordpo, icon
126 : parameter (nmax=50)
127 : real*8 eps, opeps, zero
128 : parameter (eps=1.0d-3,opeps=(1.d0+eps),zero=1.0d-12)
129 : real*8 Tmaxcv, Rgas_l, Rg_rss, Rgas_r
130 : & , acv_l(0:nordpo), acv_r(0:nordpo), acv_z(0:nordpo)
131 : & , q, k0
132 : & , rho_l, p_l, u_l, T_l, et_l, gam_l
133 : & , rho_r, p_r, u_r, T_r, et_r, gam_r
134 : & , rho_ls, p_ls, u_ls, et_ls, c_ls
135 : & , rho_z, p_rss, u_rss, et_rss, Dsw
136 : & , rho_rs, p_rs, u_rs, et_rs, c_rs
137 : real*8 Tref, pref, cref, uref, c_l, c_r, u_rm
138 : real*8 u_ls1, p_ls1, T_ls, T_ls2
139 : C real*8 du_ls_ov_dT_ls, dp_ls_ov_dT_ls
140 : real*8 der3, der4
141 : real*8 u_rss1, p_rss1, T_rs, T_rs2
142 : real*8 cell1, cell2, cell3, cell4
143 : & , cell5, cell6, cell7, cell8, cell9
144 : & , cv_rs, cv_ls
145 : real*8 T_g_z, T_g_z1
146 : C real*8 du_rss_ov_dT_rs, dp_rss_ov_dT_rs
147 : real*8 der1, der2
148 : real*8 detj, dT_ls, dT_rs
149 : real*8 umx_ls, upx_ls, rc1_ls, d_ls
150 : real*8 umx_rs, upx_rs, rc1_rs, d_rs
151 : logical lognc, logan, logeq, loe_ls, loe_rs
152 : logical logdeb
153 : parameter (logdeb = .false.)
154 : CC
155 : CC**** We write the inputs
156 : CC
157 : C write(*,*) 'racrea.f'
158 : C write(*,'(I2,E12.4)') nordpo, tmaxcv
159 : C write(*,'(6E12.4)') rgas_l, (acv_l(icon), icon=0,nordpo)
160 : C write(*,'(6E12.4)') rg_rss, (acv_z(icon), icon=0,nordpo)
161 : C write(*,'(6E12.4)') rgas_r, (acv_r(icon), icon=0,nordpo)
162 : C write(*,'(2E12.4)') q, k0
163 : C write(*,'(6E12.4)') rho_l, p_l, u_l, T_l, et_l, gam_l
164 : C write(*,'(6E12.4)') rho_r, p_r, u_r, T_r, et_r, gam_r
165 : C
166 : C**** Reference scales for a compressible solver.
167 : C For a low Mach number solver, the cut off speed is better than
168 : C the sound speed
169 : C
170 : Tref = 0.5d0 * ((T_l + T_r) + abs (T_l - T_r))
171 : C WARNING: c_l and c_r are input of funco1
172 : C and are not only used to compute the reference speed
173 : c_l = (gam_l * p_l / rho_l) ** 0.5d0
174 : c_r = (gam_r * p_r / rho_r) ** 0.5d0
175 : cref = 0.5d0 * ((c_l + c_r) + abs (c_l - c_r))
176 : uref = cref
177 : C
178 : pref = 0.5d0 * ((p_l + p_r) + abs (p_l - p_r))
179 : C
180 : C**** Redefinition of the right speed u_rm.
181 : C u_rm, used in the function funco1, which is designed to evaluate
182 : C a left GNL wave and not a right one.
182 : C Note that, after calling funco1 to compute the rs state, the
183 : C speed u_rs must be multiplied by -1.
184 : C
185 : u_rm = -1.0D0 * u_r
186 : C
187 : C**** Initialization of the states ls, rs, rss
188 : C
189 : C Starting point ls = l
190 : C
191 : u_ls = u_l
192 : p_ls = p_l
193 : T_ls = T_l
194 : C
195 : C Starting point rs = r
196 : C
197 : u_rs = u_r
198 : p_rs = p_r
199 : T_rs = T_r
200 : C We need to define *_rs since used in funrea
201 : C We initialize T_g_z = T_l since rss is a burnt state, like the
202 : C left one.
203 : rho_rs = rho_r
204 : et_rs = et_r
205 : T_g_z = T_l
206 : C
207 : C Starting point rss
208 : C NB Note that rss state is here considered function of T_rs only
208 : C (and of the guess temperature for T_rss, here called T_g_z,
209 : C in the case of thermally perfect gases).
210 : C
211 : C
211 : call funrea((nordpo+1), Tmaxcv, Rg_rss, acv_z,
212 : & rho_rs, p_rs, u_rs, et_rs, q, k0, T_g_z,
213 : & rho_z, p_rss, u_rss, Dsw, et_rss, logan)
214 : C
215 : if (logan) then
216 : write(*,'(/A8)') 'racrea.f'
217 : write(*,*) 'icon =', 0
218 : C write(*,*) 'rho_l, p_l, u_l, T_l, et_l, gam_l'
219 : C write(*,'(6E16.6)') rho_l, p_l, u_l, T_l, et_l, gam_l
220 : C write(*,*) 'rho_r, p_r, u_r, T_r, et_r, gam_r'
221 : C write(*,'(6E16.6)') rho_r, p_r, u_r, T_r, et_r, gam_r
222 : write(*,*) 'ANOMALY DETECTED 1'
223 : goto 9999
224 : endif
225 : T_g_z = p_rss / (Rg_rss * rho_z)
226 : C
227 : C**** Iterative procedure to compute the intermediate states
228 : C We use a quasi-Newton approach, with numerical evaluation of the
229 : C derivatives, involving the two variables T_rs and T_ls
230 : C (temperature of the states rs and ls).
231 : C The temperature for the rss state is
232 : C evaluated at the end of each iteration and is exact in the
232 : C particular case of calorically perfect gases. For the moment we
233 : C do not involve this variable in the quasi-Newton approach.
234 : C Summarizing we have to compute
235 : C
236 : C F = (f_1, f_2)^T = 0
237 : C
238 : C i.e.
239 : C
240 : C f_1(T_ls, T_rs; T_rss) = u_ls(T_ls, T_rs; T_rss) -
241 : C u_rss(T_ls, T_rs; T_rss) = 0
242 : C f_2(T_ls, T_rs; T_rss) = p_ls(T_ls, T_rs; T_rss) -
243 : C p_rss(T_ls, T_rs; T_rss) = 0
244 : C
245 : C At each iteration we compute X = (T_ls, T_rs)^T by solving
246 : C
247 : C F(X_1;T_rss1) + F'(X_1;T_rss1) (X - X_1) = 0
248 : C
249 : C (X - X_1) = - (F'(X_1;T_rss1))^{-1} F(X_1;T_rss1)
250 : C
251 : C Then
252 : C T_rss = fun(X,T_rss1)
253 : C
254 : do icon = 1, nmax, 1
255 : C
256 : C Linearisation of the problem around ls1 = ls
257 : C rss1 = rss
258 : C rs1 = rs
259 : C We need to save the old values of u_rs and p_rs.
260 : C
261 : u_ls1 = u_ls
262 : p_ls1 = p_ls
263 : C
264 : C Computation of of the derivatives
265 : C dp_ls/dT_ls and du_ls/dT_ls
266 : C
267 : T_ls2 = T_ls * opeps
267 : call funco1((nordpo+1), Tmaxcv, Rgas_l, acv_l,
268 : & rho_l, T_l, (Rgas_l * T_l), et_l, u_l, gam_l, c_l,
269 : & T_ls2, cell1, cv_ls, rho_ls, p_ls, et_ls, u_ls,
270 : & cell2, cell3, cell4,
271 : & cell5, cell6, cell7,
272 : & logeq, cell8, cell9)
273 : C State _ls output of funco1 only used to evaluate derivatives!
274 : der3 = (u_ls - u_ls1) / (T_ls * eps)
275 : der4 = (p_ls - p_ls1) / (T_ls * eps)
276 : C write(*,*) 'T_ls , eps', T_ls, eps
277 : C write(*,*) 'der3, der4 = ', der3, der4
278 : C
279 : C We do not need save the variables _rs in _rs1.
280 : C We need save the old values of u_rss1 and p_rss1 and the
281 : C guess temperature for rss.
282 : C
283 : T_g_z1 = T_g_z
284 : u_rss1 = u_rss
285 : p_rss1 = p_rss
286 : C
287 : C Computation of of the derivatives
288 : C dp_rss/dT_rs2 and du_rss/dT_rs2
289 : C
290 : T_rs2 = T_rs * opeps
290 : call funco1((nordpo+1), Tmaxcv, Rgas_r, acv_r,
291 : & rho_r, T_r, (Rgas_r * T_r), et_r, u_rm, gam_r, c_r,
292 : & T_rs2, cell1, cv_rs, rho_rs, p_rs, et_rs, u_rs,
293 : & cell2, cell3, cell4,
294 : & cell5, cell6, cell7,
295 : & logeq, cell8, cell9)
296 : u_rs = -1.0D0 * u_rs
297 : C -1 * necessary since funco1 is designed for left GNL waves.
298 : C State _rs output of funco1, namely
299 : C rho_rs, p_rs, et_rs, u_rs
300 : C only used to evaluate derivatives as input of funrea!
300 : call funrea((nordpo+1), Tmaxcv, Rg_rss, acv_z,
301 : & rho_rs, p_rs, u_rs, et_rs, q, k0, T_g_z,
302 : & rho_z, p_rss, u_rss, Dsw, et_rss, logan)
303 : if (logan) then
304 : write(*,'(/A8)') 'racrea.f'
305 : write(*,*) 'icon =', icon
306 : write(*,*) 'ANOMALY DETECTED 2'
307 : goto 9999
308 : endif
309 : C State _rss output of funrea only used to evaluate derivatives!
310 : der1 = (u_rss - u_rss1) / (T_rs * eps)
311 : der2 = (p_rss - p_rss1) / (T_rs * eps)
312 : C write(*,*) der1, der2
313 : C
314 : C | || | | |
315 : C |der3 -der1||dT_ls| | u_ls1 - u_rss1 |
316 : C | || | = - | |
317 : C |der4 -der2||dT_rs| | p_ls1 - p_rss1 |
318 : C | || | | |
319 : C
320 : C | | | || |
321 : C |dT_ls| 1 |der2 -der1|| u_ls1 - u_rss1 |
322 : C | | = ------ | || |
323 : C |dT_rs| det(M) |der4 -der3 || p_ls1 - p_rss1 |
324 : C | | | || |
325 : C
326 : detj = (der4 * der1)
327 : & - (der3 * der2)
328 : C
329 : dT_ls = (der2 * (u_ls1 - u_rss1)) -
330 : & (der1 * (p_ls1 - p_rss1))
331 : dT_ls = dT_ls / detj
332 : dT_rs = (der4 * (u_ls1 - u_rss1)) -
333 : & (der3 * (p_ls1 - p_rss1))
334 : dT_rs = dT_rs / detj
335 : C We limited the jump
336 : if (abs (dT_ls) .gt. (0.5d0 * T_ls)) then
337 : dT_ls = sign(dT_ls, (0.5d0 * T_ls))
338 : endif
339 : if (abs (dT_rs) .gt. (0.5d0 * T_rs)) then
340 : dT_rs = sign(dT_rs, (0.5d0 * T_rs))
341 : endif
342 : C
343 : T_ls = T_ls + dT_ls
344 : T_rs = T_rs + dT_rs
345 : C
346 : C Computation of the states.
347 : C ls
347 : call funco1((nordpo+1), Tmaxcv, Rgas_l, acv_l,
348 : & rho_l, T_l, (Rgas_l * T_l), et_l, u_l, gam_l, c_l,
349 : & T_ls, cell1, cv_ls, rho_ls, p_ls, et_ls, u_ls,
350 : & cell2, umx_ls, upx_ls,
351 : & rc1_ls, cell6, cell7,
352 : & loe_ls, cell8, cell9)
353 : C rs
353 : call funco1((nordpo+1), Tmaxcv, Rgas_r, acv_r,
354 : & rho_r, T_r, (Rgas_r * T_r), et_r, u_rm, gam_r, c_r,
355 : & T_rs, cell1, cv_rs, rho_rs, p_rs, et_rs, u_rs,
356 : & cell2, umx_rs, upx_rs,
357 : & rc1_rs, cell6, cell7,
358 : & loe_rs, cell8, cell9)
359 : u_rs = -1.0D0 * u_rs
360 : C rss
360 : call funrea((nordpo+1), Tmaxcv, Rg_rss, acv_z,
361 : & rho_rs, p_rs, u_rs, et_rs, q, k0, T_g_z,
362 : & rho_z, p_rss, u_rss, Dsw, et_rss, logan)
363 : if (logan) then
364 : write(*,'(/A8)') 'racrea.f'
365 : write(*,*) 'icon =', icon
366 : write(*,*) 'ANOMALY DETECTED 3'
367 : goto 9999
368 : endif
369 : C
370 : C Computation of T_g_z
371 : C (T_rss = fun(X,T_rss1)).
372 : C
373 : T_g_z = p_rss / (Rg_rss * rho_z)
374 : C
375 : C write(*,*) icon, T_ls, T_rs, T_g_z
376 : C
377 : lognc = abs (dT_rs) .lt. (zero * tref)
378 : lognc = lognc .and. (abs (dT_ls) .lt. (zero * tref))
379 : lognc = lognc .and. (abs (T_g_z - T_g_z1) .lt.
380 : & (zero * tref))
381 : lognc = lognc .and. (abs(p_rss - p_ls) .lt. (zero * pref))
382 : lognc = lognc .and. (abs(u_rss - u_ls) .lt. (zero * uref))
383 : if (lognc) goto 9998
384 : enddo
385 : 9998 continue
386 : C write(*,*) icon
387 : C
388 : C**** Sound speed and shock speed
389 : C
390 : c_ls = ((Rgas_l + cv_ls) * p_ls) / (cv_ls * rho_ls)
391 : c_ls = c_ls ** 0.5d0
392 : c_rs = ((Rgas_r + cv_rs) * p_rs) / (cv_rs * rho_rs)
393 : c_rs = c_rs ** 0.5d0
394 : lognc = .not. lognc
395 : if (T_rs .le. (1.00001D0 * T_r)) then
396 : C In case of rarefaction wave we take d_rs = d_r
397 : d_rs = u_r + c_r
398 : else
399 : d_rs = u_r + sqrt(rc1_rs / umx_rs / upx_rs)
400 : endif
401 : if (T_ls .le. (1.00001D0 * T_l)) then
402 : C In case of rarefaction wave we take d_ls = d_l
403 : d_ls = u_l - c_l
404 : else
405 : d_ls = u_l - sqrt(rc1_ls / umx_ls / upx_ls)
406 : endif
407 : if (lognc) then
408 : write(*,'(/A8)') 'racrea.f'
409 : write(*,*) 'Quasi-Newton method does not converge'
410 : endif
411 : if (logdeb) then
412 : write(*,*) 'Intersection '
413 : write(*,*) rho_l, p_l, u_l, T_l
414 : write(*,*) rho_ls, p_ls, u_ls, T_ls, d_ls, loe_rs
415 : write(*,*) rho_z, p_rss, u_rss, T_g_z
416 : write(*,*) rho_rs, p_rs, u_rs, T_rs, d_rs, loe_rs
417 : write(*,*) rho_r, p_r, u_r, T_r
418 : endif
419 : C
420 : 9999 continue
421 : return
422 : end
423 :
© Cast3M 2003 - All rights reserved.
Disclaimer