1 : C FUNCJ SOURCE BECC 09/12/07 21:15:28 6579
2 : subroutine funcj(nordpo, Tmaxcv, acv_r, Rgas_r, acv_b, Rgas_b,
3 : & k0, q,
4 : & T_ag, T_bg,
5 : & p_r, T_r, u_r,
6 : & logdef,
7 : & r_acj, T_acj, p_acj, u_acj, d_acj,
8 : & r_bcj, T_bcj, p_bcj, u_bcj, d_bcj,
9 : & k0cjdt, logan)
10 : C
11 : C INPUT
12 : C
13 : C nordpo = order of polynomial for cp and cv (see also
14 : C Tmaxcv)
15 : C
16 : C Tmaxcv = maximum temperature for cv polynomial expansion
17 : C cv(T) = cv(Tmaxcv) if T > Tmaxcv
18 : C
19 : C acv_r = to compute cv and ether for r;
20 : C vector such that
21 : C cv = \sum_{i=1,nordpo+1} acv(i) T^{i-1}
22 : C ether = \sum_{i=1,nordpo+1} acv(i) T^{i} / (i)
23 : C
24 : C Rgas_r = gas constant for r
25 : C
26 : C acv_b = to compute cv and ether for b (rss) ;
27 : C
28 : C Rgas_b = gas constant for (rss)
29 : C
30 : C k0 = fundamental flame speed
31 : C
32 : C q = released heat rs->rss
33 : C
34 : C T_ag, T_bg = guess temperature for a (rs) and b (rss)
35 : C used to evaluate gamma and linearisation of
36 : C enthalpies
37 : C
38 : C p_r, T_r, u_r = pressure, temperature and velocities for r
39 : C
40 : C
41 : C
42 : C OUTPUT
43 : C
44 : C logdef = logical. If .true., CJDF is admissible.
45 : C Otherwise, it is not and the output is CJDT
46 : C state
47 : C
48 : C d_bcj = reactive shock speed
49 : C
50 : C d_acj = non-reactive shock speed in CJDF case
51 : C reactive shock speed in CJDT case
52 : C
53 : C r_acj, T_acj,
54 : C p_acj, u_acj = unburnt state behind the CJDF or CJDT
55 : C
56 : C r_bcj, T_bcj,
57 : C p_bcj, u_bcj = unburnt state behind the CJDF or CJDT
58 : C
59 : C logan = if .true., an anomaly has been detected
60 : C
61 : c implicit none
62 : integer nordpo
63 : real*8 Tmaxcv, acv_r(1:(nordpo+1)), Rgas_r
64 : & , acv_b(1:(nordpo+1)), Rgas_b, k0, q
65 : & , T_ag, T_bg, p_r, T_r, u_r
66 : & , r_acj, T_acj, p_acj, u_acj, d_acj
67 : & , r_bcj, T_bcj, p_bcj, u_bcj, d_bcj
68 : real*8 zero
69 : parameter(zero = 0.0D0)
70 : real*8 et_a, cv_a, et_b, cv_b, gam_a, gam_b
71 : & , qti, qtiti, a, b, c
72 : & , ht_acj, m
73 : & , cv_r, et_r, k0cjdt, kcjdt
74 : & , r_r
75 : logical logan, logdef, logdeb
76 : parameter(logdeb = .false.)
77 : if (logdeb) then
78 : write(*,*) 'Subroutine funcj.f'
79 : endif
80 : C
81 : logan = .false.
82 : C
83 : C Properties (gam_a, gam_b)
84 : C
84 : call prith1(nordpo, acv_r, Tmaxcv, T_ag, et_a, cv_a)
85 : gam_a = (cv_a + Rgas_r) / cv_a
85 : call prith1(nordpo, acv_b, Tmaxcv, T_bg, et_b, cv_b)
86 : gam_b = (cv_b + Rgas_b) / cv_b
87 : C
88 : C**** Correction of q as it is zero
89 : C
90 : q = max(q, (zero * Rgas_r * T_r))
91 : qti = q + ((cv_b * T_bg) - et_b)
92 : C qtiti used in case of CJDF
93 : qtiti = qti - ((cv_a * T_ag) - et_a)
94 : qtiti = max(qtiti, 0.0D0)
95 : if (logdeb )then
96 : write(*,*) 'gam_a, gam_b', gam_a, gam_b
97 : write(*,*) 'q, qti, qtiti'
98 : write(*,*) q, qti, qtiti
99 : endif
100 : C
101 : C**** Computation of the CJDT
102 : C It should be done here, in order to know K0,CJDT
103 : C
103 : call prith1(nordpo, acv_r, Tmaxcv, T_r, et_r, cv_r)
104 : a = 1.0D0 / gam_b
105 : a = a * a
106 : c = Rgas_r * T_r
107 : c = c * c
108 : b = 1.0D0 - a
109 : b = b * ((et_r + (Rgas_r * T_r)) + qti)
110 : b = b - (Rgas_r * T_r)
111 : kcjdt = (b*b) - (a*c)
112 : if ((kcjdt .lt. 0.0D0))then
113 : kcjdt = 0.0D0
114 : if (logdeb) then
115 : write(*,*) 'subroutine funcj'
116 : write(*,*) 'Computation of kcjdt'
117 : write(*,*) 'Delta < 0', kcjdt
118 : write(*,*) 'Delta = ', kcjdt
119 : write(*,*) 'b*b = ', b*b
120 : write(*,*) 'sqrt(|Delta|)/b = ', (sqrt(abs(kcjdt))/b)
121 : endif
122 : C logan = .true.
123 : C goto 9999
124 : endif
125 : kcjdt = b + sqrt(kcjdt)
126 : kcjdt = kcjdt / a
127 : if (kcjdt .lt. 0.0D0)then
128 : write(*,*) 'subroutine funcj'
129 : write(*,*) 'Computation of kcjdt'
130 : write(*,*) 'kcjdt*kcjdt < 0'
131 : logan = .true.
132 : goto 9999
133 : endif
134 : kcjdt = sqrt(kcjdt)
135 : d_bcj = kcjdt + u_r
136 : C
137 : C**** Computation of T_bcj (case of cjdt detonation)
138 : C
139 : r_r = P_r / (Rgas_r * T_r)
140 : r_bcj = (gam_b - 1.0D0) / (gam_b + 1.0D0)
141 : r_bcj = r_bcj *
142 : & ((2.0D0 * ((et_r + (Rgas_r * T_r)) + qti) / (kcjdt*kcjdt))
143 : & + 1.0D0)
144 : if (r_bcj .lt. 0.0D0) then
145 : write(*,*) 'subroutine funcj'
146 : write(*,*) 'Computation of r_bcj'
147 : write(*,*) 'r_bcj < 0'
148 : logan = .true.
149 : goto 9999
150 : endif
151 : r_bcj = sqrt(r_bcj)
152 : r_bcj = r_r / r_bcj
153 : m = r_r * kcjdt
154 : p_bcj = p_r -
155 : & ((m * m) * ((1.0D0 / r_bcj) - (1.0D0 / r_r)))
156 : if (p_bcj .lt. 0.0D0) then
157 : write(*,*) 'subroutine funcj'
158 : write(*,*) 'Computation of p_bcj'
159 : write(*,*) 'p_bcj < 0'
160 : logan = .true.
161 : goto 9999
162 : endif
163 : T_bcj = p_bcj / (Rgas_b * r_bcj)
164 : u_bcj = d_bcj - (m / r_bcj)
165 : C
166 : C******* Computation of VN state
167 : C
168 : b = gam_a / (gam_a + 1.0D0)
169 : b = b * (1.0D0 + (P_r / (r_r * kcjdt * kcjdt)))
170 : c = 2.0D0 * ((et_r + (Rgas_r * T_r)) - (et_a - (cv_a * T_ag)))
171 : c = c / (kcjdt * kcjdt)
172 : c = c + 1.0D0
173 : c = c * (gam_a - 1.0D0) / (gam_a + 1.0D0)
174 : r_acj = b*b - c
175 : if (r_acj .lt. 0.0D0) then
176 : r_acj = 0.0D0
177 : if (logdeb) then
178 : write(*,*) 'subroutine funcj'
179 : write(*,*) 'Computation of r_acj'
180 : write(*,*) 'delta < 0', r_acj
181 : endif
182 : C logan = .true.
183 : C goto 9999
184 : endif
185 : r_acj = b - sqrt( r_acj )
186 : r_acj = r_R / r_acj
187 : m = kcjdt * r_R
188 : P_acj = (1.0D0 / r_acj) - (1.0D0 / r_r)
189 : P_acj = P_R - (m**2 * P_acj)
190 : if (P_acj .lt. 0.0D0) then
191 : write(*,*) 'subroutine funcj'
192 : write(*,*) 'Computation of P_acj'
193 : write(*,*) 'P_acj < 0'
194 : logan = .true.
195 : goto 9999
196 : endif
197 : T_acj = P_acj / (Rgas_R * r_acj)
198 : u_acj = D_BCJ - (m / r_acj)
199 : k0cjdt = d_bcj - u_acj
200 : CC
201 : CC******* Test. We write the values thus determined
202 : CC
203 : C write(*,*) 'T_acj, p_acj, r_acj, u_acj '
204 : C write(*,*) T_acj, p_acj, r_acj, u_acj
205 : C write(*,*) 'T_bcj, p_bcj, r_bcj, u_bcj, d_bcj'
206 : C write(*,*) T_bcj, p_bcj, r_bcj, u_bcj, d_bcj
207 : C stop
208 : C
209 : if (k0 .ge. k0cjdt) then
210 : logdef = .false.
211 : d_acj = d_bcj
212 : else
213 : C
214 : C CJDF
215 : C
216 : logdef = .true.
217 : a = 1.0D0 / gam_b
218 : a = a * a
219 : b = 1.0D0 - a
220 : b = b * (gam_a / (gam_a - 1.0D0))
221 : b = b - 1.0D0
222 : b = b * (k0 * k0)
223 : c = 2.0D0 * (k0 * k0)
224 : c = c * (1.0D0 - a)
225 : c = c * qtiti
226 : c = c - ((k0**4) * a)
227 : T_acj = b * b + c
228 : if (T_acj .lt. 0.0D0) then
229 : T_acj = 0.0D0
230 : if (logdeb) then
231 : write(*,*) 'subroutine funcj'
232 : write(*,*) 'computation of T_acj'
233 : write(*,*) 'Delta < 0', T_acj
234 : endif
235 : C logan = .true.
236 : C goto 9999
237 : endif
238 : T_acj = sqrt(t_acj)
239 : T_acj = T_acj + b
240 : T_acj = T_acj / Rgas_r
241 : C
242 : C CJDF state in case of right-travelling shock
243 : C
243 : call funsho(nordpo, Tmaxcv, acv_r, Rgas_r,
244 : & P_r, T_r, u_r, .true.,
245 : & T_acj, p_acj, r_acj, u_acj, d_acj)
246 : C
247 : r_bcj = (gam_b - 1.0D0) / (gam_b + 1.0D0)
248 : ht_acj = gam_a / (gam_a - 1.0D0) * Rgas_r * T_acj
249 : r_bcj = r_bcj *
250 : & ((2.0D0 * (ht_acj + qtiti) / (k0*k0))
251 : & + 1.0D0)
252 : if (r_bcj .lt. 0.0D0) then
253 : write(*,*) 'subroutine funcj'
254 : write(*,*) 'computation of r_bcj'
255 : write(*,*) 'Delta < 0'
256 : logan = .true.
257 : goto 9999
258 : endif
259 : r_bcj = sqrt(r_bcj)
260 : r_bcj = r_acj / r_bcj
261 : m = r_acj * k0
262 : p_bcj = p_acj -
263 : & ((m * m) * ((1.0D0 / r_bcj) - (1.0D0 / r_acj)))
264 : if (p_bcj .lt. 0.0D0) then
265 : write(*,*) 'subroutine funcj'
266 : write(*,*) 'Computation of P_bcj'
267 : write(*,*) 'p_bcj < 0'
268 : logan = .true.
269 : goto 9999
270 : endif
271 : T_bcj = p_bcj / (Rgas_b * r_bcj)
272 : d_bcj = k0 + u_acj
273 : u_bcj = d_bcj - (m/r_bcj)
274 : CC
275 : CC**** Test. We write the values thus determined
276 : CC
277 : C write(*,*) 'T_acj, p_acj, r_acj, u_acj, d_acj '
278 : C write(*,*) T_acj, p_acj, r_acj, u_acj, d_acj
279 : C write(*,*) 'T_bcj, p_bcj, r_bcj, u_bcj, d_bcj'
280 : C write(*,*) T_bcj, p_bcj, r_bcj, u_bcj, d_bcj
281 : C stop
282 : C
283 : if (d_bcj .gt. d_acj) then
284 : if (logdeb) then
285 : write(*,*) 'subroutine funcj'
286 : write(*,*) 'Deflagration regime'
287 : write(*,*) 'd_bcj, d_acj ', d_bcj, d_acj
288 : write(*,*)
289 : $ 'Reactive shock faster than the non-reactive one'
290 : endif
291 : C logan = .true.
292 : C goto 9999
293 : endif
294 : endif
295 : C write(*,*) 'k0cjdt, d_acj, d_bcj ', k0cjdt, d_acj, d_bcj
296 : C
297 : 9999 continue
298 : return
299 : end
300 : C
301 :
© Cast3M 2003 - All rights reserved.
Disclaimer