1 : C RACNRE SOURCE BECC 09/12/07 21:15:59 6579
2 : subroutine racnre(nordpo, Tmaxcv,
3 : & Rgas_l, acv_l, Rgas_r, acv_r,
4 : & r_l, p_l, u_l, T_l, gam_l,
5 : & r_r, p_r, u_r, T_r, gam_r,
6 : & d_l, d_r,
7 : & r_b, p_b, u_b, T_b, d_b,
8 : & r_a, p_a, u_a, T_a, d_a,
9 : & logan, lognc)
10 : C
11 : c************************************************************************
12 : c
13 : c project :
14 : c
15 : c name : racnre
16 : c
17 : c description : computation of the shock-shock field-by-field
18 : c decomposition of the non-reactive Riemann problem
19 : C for thermally perfect gases.
20 : C
21 : C We have to distinguish between 4 different states
22 : C
23 : C l = left state
24 : C |
25 : C | Left GNL wave
26 : C |
27 : C b = left star state;
28 : C connected to l by a GNL wave
29 : C |
30 : C | Contact discontinuity
31 : C |
32 : C a = right star state;
33 : C connected to r by a GNL wave
34 : C |
35 : C | Right GNL wave
36 : C |
37 : C r = right state;
38 : C unburnt state
39 : c
40 : c language : fortran 77
41 : c
42 : c author : a. beccantini den/dm2s/sfme/ltmf
43 : c
44 : c************************************************************************
45 : c
45 : c called by : ...
46 : c
46 : c called program : ...
47 : c
48 : c************************************************************************
49 : c
50 : c
51 : c**** input
52 : c
53 : c nordpo = polynomials degree of cv(t)
54 : c tmaxcv = threshold temperature for cv(t)
55 : c rgas_l = gas constant (J/kg/K)in l (and in ls)
56 : c acv_l = cv_l = \sum_{i=0,nordpo} (acv_l{i} T^i)
57 : c rgas_r = gas constant r (and in rs)
58 : c acv_r = cv_r = \sum_{i=0,nordpo} (acv_r{i} T^i)
59 : c
60 : C r_l = density in l
61 : C p_l = pressure in
62 : c u_l = velocity in l
63 : C T_l = temperature in l (redundant with r_l and
64 : c p_l)
65 : c gam_l = specific heat ratio in l (redundant
66 : c with r_l and p_l)
67 : c
68 : C r_r = density in r
69 : C p_r = pressure in r
70 : c u_r = velocity in r
71 : C T_r = temperature in r (redundant with r_r and
72 : c p_r)
73 : c gam_r = specific heat ratio in r (redundant
74 : c with r_r and p_r)
75 : c
76 : c**** output
77 : C
78 : C d_r, d_l = speed of states l and r
79 : c
80 : c
81 : C r_b = density in b
82 : C p_b = pressure in b
83 : c u_b = velocity in b
84 : c T_b = temperature in b
85 : c d_b = shock speed of the left shock l->b
86 : c
87 : C r_a = density in a
88 : C p_a = pressure in a
89 : c u_a = velocity in a
90 : c T_a = temperature in a
91 : c d_a = shock speed of the
92 : C right shock a->r
93 : c
94 : c logan = an anomaly is detected
95 : c lognc = the approach has some problems
96 : C
97 : c
98 : c************************************************************************
99 : c
100 : c 07/12/2009 Created.
101 : c
102 : c implicit none
103 : integer nordpo
104 : integer nmax
105 : parameter (nmax = 20)
106 : integer iter, icon
107 : real*8 eps
108 : parameter(eps = 1.0D-6)
109 : real*8 Tmaxcv, Rgas_l, Rgas_r
110 : & , acv_l(0:nordpo), acv_r(0:nordpo)
111 : & , r_l, p_l, u_l, T_l, gam_l
112 : & , r_r, p_r, u_r, T_r, gam_r
113 : & , d_l, d_r
114 : & , r_b, p_b, u_b, T_b, d_b
115 : & , r_a, p_a, u_a, T_a, d_a
116 : C
117 : real*8 Tref, pref, Dtmin, coefder
118 : real*8 T_ag, T_bg
119 : real*8 P_ag, r_ag, u_ag, d_ag
120 : real*8 T_agp, P_agp, r_agp, u_agp, d_agp
121 : real*8 dpdua
122 : real*8 P_bg, r_bg, u_bg, d_bg
123 : real*8 T_bgp, P_bgp, r_bgp, u_bgp, d_bgp
124 : real*8 dpdub
125 : real*8 a, b
126 : real*8 erro
127 : logical lognc, logan, logdeb
128 : parameter(logdeb = .false.)
129 : CC
130 : CC**** We write the inputs
131 : CC
132 : if( logdeb ) then
133 : write(*,*) ' '
134 : write(*,*) 'RP'
135 : write(*,*) 'Left state'
136 : write(*,*) 'R, cv'
137 : write(*,'(6E12.4)') rgas_l, (acv_l(icon), icon=0,nordpo)
138 : write(*,*) 'rho, p, u, T, gamma'
139 : write(*,'(6E12.4)') r_l, p_l, u_l, T_l, gam_l
140 : write(*,*) 'Right state'
141 : write(*,*) 'R, cv'
142 : write(*,'(6E12.4)') rgas_r, (acv_r(icon), icon=0,nordpo)
143 : write(*,*) 'rho, p, u, T, gamma'
144 : write(*,'(6E12.4)') r_r, p_r, u_r, T_r, gam_r
145 : endif
146 : C stop
147 : C
148 : C**** Reference scales for a compressible solver.
149 : C For a low Mach number solver, the cut off speed is better than
150 : C the sound speed
151 : C
152 : Tref = 0.5d0 * ((T_l + T_r) + abs (T_l - T_r))
153 : Dtmin = 1.0D-3 * Tref
154 : coefder = 1.001D0
155 : pref = 0.5d0 * ((p_l + p_r) + abs (p_l - p_r))
156 : C
157 : C Initialisation of temperatures in a, b, c, d
158 : C
159 : T_a = T_r
160 : T_b = T_l
161 : C
162 : do iter = 1, nmax
163 : C
164 : C Guess temperatures at this iteration
165 : C
166 : T_ag = T_a
167 : T_bg = T_b
168 : C
169 : C******* Right state
170 : C
170 : call funsho( nordpo, Tmaxcv, acv_r, Rgas_r,
171 : & P_r, T_r, u_r, .true.,
172 : & T_ag, P_ag, r_ag, u_ag, d_ag)
173 : C
174 : T_agp = max((T_ag + Dtmin), (T_ag * coefder))
174 : call funsho( nordpo, Tmaxcv, acv_r, Rgas_r,
175 : & P_r, T_r, u_r, .true.,
176 : & T_agp, P_agp, r_agp, u_agp, d_agp)
177 : C
178 : dpdua = (P_agp - P_ag) / (u_agp - u_ag)
179 : if (dpdua .lt. 0.0d0) then
180 : write(*,*) 'dpdua > 0 , ', dpdua
181 : logan = .true.
182 : write(*,*) 'subroutine racnre.f'
183 : write(*,*) 'an error has been detected'
184 : goto 9999
185 : endif
186 : C
187 : C******* Left state
188 : C
188 : call funsho( nordpo, Tmaxcv, acv_l, Rgas_l,
189 : & P_l, T_l, u_l, .false.,
190 : & T_bg, P_bg, r_bg, u_bg, d_bg)
191 : C
192 : T_bgp = max((T_bg + Dtmin), (T_bg * coefder))
192 : call funsho( nordpo, Tmaxcv, acv_l, Rgas_l,
193 : & P_l, T_l, u_l, .false.,
194 : & T_bgp, P_bgp, r_bgp, u_bgp, d_bgp)
195 : C
196 : dpdub = (P_bgp - P_bg) / (u_bgp - u_bg)
197 : if (dpdub .gt. 0.0d0) then
198 : write(*,*) 'dpdub > 0 , ', dpdub
199 : logan = .true.
200 : write(*,*) 'subroutine racnre.f'
201 : write(*,*) 'an error has been detected'
202 : goto 9999
203 : endif
204 : C
205 : C******* Intersection
206 : C
207 : a = dpdua - dpdub
208 : b = ((dpdua * u_ag) - (dpdub * u_bg)) + (p_bg - p_ag)
209 : u_b = b / a
210 : p_b = p_bg + (dpdub * (u_b - u_bg))
211 : p_a = p_ag + (dpdua * (u_b - u_ag))
212 : if ( abs (p_a - p_b) .gt. (eps * pref)) then
213 : write(*,*) 'pa - pb too big, ', (p_a - p_b)
214 : logan = .true.
215 : write(*,*) 'subroutine racnre.f'
216 : write(*,*) 'an error has been detected'
217 : goto 9999
218 : endif
219 : C p_a, p_b should be positive
220 : p_a = max (p_a, (eps * p_ag))
221 : p_b = max (p_b, (eps * p_bg))
222 : C
223 : C******* Computation of temperature T_b, T_a
224 : C
225 : T_b = T_bg + ((p_b - p_bg) * (T_bgp - T_bg) / (p_bgp - p_bg))
226 : T_b = max (T_b, (eps * T_bg))
227 : C
228 : T_a = T_ag + ((p_a - p_ag) * (T_agp - T_ag) / (p_agp - p_ag))
229 : T_a = max (T_a, (eps * T_ag))
230 : C
231 : C******* TEST
232 : C
233 : erro = abs (T_a - T_ag) + abs (T_b - T_bg)
234 : erro = erro / Tref
235 : if (logdeb) then
236 : write(*,*) 'Intersection'
237 : write(*,*) 'T_ag, T_bg = ', T_ag, T_bg
238 : write(*,*) 'dpdua, dpdub = ', dpdua, dpdub
239 : write(*,*) 'p_a, p_b = ', p_a, p_b
240 : write(*,*) 'T_a, T_b = ', T_a, T_b
241 : WRITE(*,*) 'erro =', erro
242 : endif
243 : if (erro .le. 1.0D-14) then
244 : lognc = .false.
245 : goto 9998
246 : endif
247 : enddo
248 : if (logdeb ) then
249 : write(*,*) 'Warning. Convergence non-reached'
250 : write(*,*) 'erro = ', erro
251 : endif
252 : lognc = .true.
253 : 9998 continue
253 : call funsho( nordpo, Tmaxcv, acv_r, Rgas_r,
254 : & P_r, T_r, u_r, .true.,
255 : & T_a, P_a, r_a, u_a, d_a)
255 : call funsho( nordpo, Tmaxcv, acv_l, Rgas_l,
256 : & P_l, T_l, u_l, .false.,
257 : & T_b, P_b, r_b, u_b, d_b)
258 : d_l = gam_l * P_l / r_l
259 : d_l = u_l - sqrt(d_l)
260 : d_l = min(d_l, d_b)
261 : d_r = gam_r * P_r / r_r
262 : d_r = u_r + sqrt(d_r)
263 : d_r = max(d_r, d_a)
264 : C
265 : if (logdeb) then
266 : write(*,*) 'Final result'
267 : write(*,*) 'iter = ', iter
268 : write(*,*) 'r. rho, p, u, T, d'
269 : write(*,*) r_r, p_r, u_r, T_r, d_r
270 : write(*,*) 'a. rho, p, u, T, d'
271 : write(*,*) r_a, p_a, u_a, T_a, d_a
272 : write(*,*) 'b. rho, p, u, T, d'
273 : write(*,*) r_b, p_b, u_b, T_b, d_b
274 : write(*,*) 'l. rho, p, u, T, d'
275 : write(*,*) r_l, p_l, u_l, T_l, d_l
276 : C write(*,*) 'acv'
277 : C do iter = 0, nordpo, 1
278 : C write(*,*) acv_b(iter), acv_l(iter)
279 : C enddo
280 : endif
281 : 9999 continue
282 : return
283 : end
284 :
© Cast3M 2003 - All rights reserved.
Disclaimer