1 : C FUNRE1 SOURCE BECC 09/12/07 21:15:31 6579
2 : subroutine funre1(nordpo, tmaxcv, rgas_a, acv_a,
3 : & rgas_b, acv_b, q, k0, T_ag, T_bg,
4 : & r_a, p_a, u_a, T_a,
5 : & r_b, p_b, u_b, T_b, Dsw, logan)
6 : c************************************************************************
7 : c
8 : c projet :
9 : c
10 : c name : funre1
11 : c
12 : c description : computation of the weak state behind a reactive
13 : c shock
14 : c
15 : c language : fortran 77
16 : c
17 : c author : a. beccantini den/dm2s/sfme/ltmf
18 : c
19 : c************************************************************************
20 : c
20 : c callen by : racre2
21 : c
21 : c called program : none
22 : c
23 : c************************************************************************
24 : c
25 : c
26 : c**** input
27 : c
28 : c nordpo = polynomials degree
29 : c tmaxcv = threshold temperature for cv(t)
30 : c rgas_a = r of the gas in r and rs
31 : c acv_a = \sum_{i=1,nesp} y_i*acv_a{i,j}
32 : c rgas_b = r of the gas in rss
33 : c acv_b = \sum_{i=1,nesp} y_i*acv_b{i,j}
34 : c q = realeased heat per mass unit
35 : c k0 = fundamental flame speed
36 : C T_ag,T_bg = guess temperatures for a and b
37 : C (to compute properties)
38 : C r_a = density in rs
39 : C p_a = pressure in rs
40 : c u_a = velocity in rs
41 : C T_a = temperature in rs
42 : c**** output
43 : c r_b = approximate density in rss (behind the shock)
44 : c p_b = approximate pressure in rss
45 : c u_b = approximate speed in rss
46 : C T_b = approximate temperature in rss
47 : C (NB if T_bg = T_b and T_ag = T_a, rss thus
48 : C computed is exact)
49 : C Dsw = approximate shock wave speed... it could
50 : C differs from u_a + k0 if (k0lim < k0)
51 : C Nevertheless, if logsto = .true., if k0lim
52 : C < k0 we stop!
53 : c logan = an anomaly is detected
54 : c
55 : c************************************************************************
56 : c
57 : c created the 07/12/2009
58 : c
59 : c implicit none
60 : integer nordpo
61 : real*8 tmaxcv, rgas_a, acv_a(1:(nordpo+1))
62 : & ,rgas_b, acv_b(1:(nordpo+1))
63 : real*8 q, k0
64 : real*8 T_ag, T_bg
65 : real*8 r_a, p_a, u_a, T_a
66 : real*8 r_b, p_b, u_b, T_b, Dsw
67 : real*8 gam_a, e_a, cv_a, gam_b, e_b, cv_b
68 : real*8 ttq, k02ad, h_a
69 : real*8 a, b, c, delta, delta1, k0cj2, k0new
70 : real*8 csi
71 : logical logan
72 : CC
73 : CC**** We write the input
74 : CC
75 : C write(*,*) 'We enter funre1'
76 : C
77 : C**** Computation of gam_a and gam_b
78 : C
78 : call prith1(nordpo, acv_a, Tmaxcv, T_ag, e_a, cv_a)
79 : gam_a = (cv_a + Rgas_a) / cv_a
79 : call prith1(nordpo, acv_b, Tmaxcv, T_bg, e_b, cv_b)
80 : gam_b = (cv_b + Rgas_b) / cv_b
81 : ttq = q + ((cv_b * T_bg) - e_b) -
82 : & ((cv_a * T_ag) - e_a)
83 : C h_a and h_b computed as we deal with cp gases
84 : h_a = (gam_a / (gam_a - 1.0D0)) * Rgas_a * T_a
85 : C
86 : r_a = P_a / (Rgas_a * T_a)
87 : C
88 : C**** Computation of r_b
89 : C
90 : k02ad = k0 * k0 / (Rgas_a * T_a)
91 : if (k02ad .lt. 1.0D0)then
92 : C if ( .true. ) then
93 : C if ( .false. ) then
94 : C write(*,*) 'funre1.f. Parte 1'
95 : c = k02ad
96 : b = k02ad + 1.0D0
97 : b = b * gam_b / (gam_b + 1.0D0)
98 : a = 2.0D0 * (h_a + ttq) / (Rgas_a * T_a)
99 : a = a + k02ad
100 : a = a * (gam_b - 1.0D0) / (gam_b + 1.0D0)
101 : delta = (b * b) - (a * c)
102 : if (delta .lt. 0.0D0) then
103 : C if ( .true. ) then
104 : C write(*,*) 'funre1.f. Parte 1.2'
105 : C We are in CJDF regime
106 : a = 1.0D0 / (gam_b * gam_b)
107 : b = 1.0D0 - a
108 : b = b * (h_a + ttq)
109 : b = b - (Rgas_a * T_a)
110 : c = (Rgas_a * T_a)**2
111 : delta1 = (b * b) - (a * c)
112 : if (delta1 .lt. 0.0D0) then
113 : write(*,*) 'subroutine funre1.f'
114 : write(*,*) 'We cannot compute k0cj'
115 : write(*,*) 'ANOMALY DETECTED'
116 : logan = .true.
117 : goto 9999
118 : endif
119 : k0cj2 = (b - sqrt(delta1)) / a
120 : if (k0cj2 .lt. 0.0D0) then
121 : write(*,*) 'subroutine funre1.f'
122 : write(*,*) 'k0cj^2 < 0'
123 : write(*,*) 'ANOMALY DETECTED'
124 : logan = .true.
125 : goto 9999
126 : endif
127 : k02ad = k0cj2 / (Rgas_a * T_a)
128 : b = k02ad + 1.0D0
129 : b = b * gam_b / (gam_b + 1.0D0)
130 : a = 2.0D0 * (h_a + ttq) / (Rgas_a * T_a)
131 : a = a + k02ad
132 : a = a * (gam_b - 1.0D0) / (gam_b + 1.0D0)
133 : csi = b / a
134 : k0new = k0cj2 ** 0.5D0
135 : else
136 : k0new = k0
137 : csi = b + sqrt(delta)
138 : csi = csi / a
139 : endif
140 : r_b = r_a * csi
141 : else
142 : C write(*,*) 'funre1.f. Parte 2'
143 : c = 2.0D0 * (h_a + ttq)
144 : c = c / (k0 * k0)
145 : c = c + 1.0D0
146 : c = c * (gam_b - 1.0D0) / (gam_b + 1.0D0)
147 : b = Rgas_a * T_a / (k0 * k0)
148 : b = b + 1.0D0
149 : b = b * gam_b / (gam_b + 1.0D0)
150 : delta = (b * b) - c
151 : if (delta .lt. 0.0D0) then
152 : C if ( .true. ) then
153 : C write(*,*) 'funre1.f. Parte 2.2'
154 : a = 1.0D0 / (gam_b * gam_b)
155 : b = 1.0D0 - a
156 : b = b * (h_a + ttq)
157 : b = b - (Rgas_a * T_a)
158 : c = (Rgas_a * T_a)
159 : c = c * c
160 : delta1 = (b * b) - (a * c)
161 : if (delta1 .lt. 0.0D0) then
162 : write(*,*) 'subroutine funre1.f'
163 : write(*,*) 'We cannot compute k0cj'
164 : write(*,*) 'ANOMALY DETECTED'
165 : logan = .true.
166 : goto 9999
167 : endif
168 : k0cj2 = (b - sqrt(delta1)) / a
169 : if (k0cj2 .lt. 0.0D0) then
170 : write(*,*) 'subroutine funre1.f'
171 : write(*,*) 'k0cj^2 < 0'
172 : write(*,*) 'ANOMALY DETECTED'
173 : logan = .true.
174 : goto 9999
175 : endif
176 : c = 2.0D0 * (h_a + ttq)
177 : c = c / k0cj2
178 : c = c + 1.0D0
179 : c = c * (gam_b - 1.0D0) / (gam_b + 1.0D0)
180 : b = Rgas_a * T_a / k0cj2
181 : b = b + 1.0D0
182 : b = b * gam_b / (gam_b + 1.0D0)
183 : delta1 = b * b - c
184 : C if (abs (delta1) > 1.0D-6) then
185 : if (abs (delta1) .gt. 1.0D-6) then
186 : write(*,*) 'subroutine funre1.f'
187 : write(*,*) 'CJDF'
188 : write(*,*) 'delta != 0'
189 : write(*,*) 'ANOMALY DETECTED'
190 : logan = .true.
191 : goto 9999
192 : endif
193 : csi = b
194 : k0new = k0cj2 ** 0.5D0
195 : else
196 : csi = b - sqrt (delta)
197 : k0new = k0
198 : endif
199 : r_b = r_a / csi
200 : endif
201 : csi = r_a / r_b
202 : Dsw = k0new + u_a
203 : u_b = Dsw - (k0new * csi)
204 : p_b = p_a + (r_a * k0new * k0new * (1.0d0 - csi))
205 : if (p_b .lt. 0.0D0) then
206 : write(*,*) 'subroutine funre1.f'
207 : write(*,*) 'p_b < 0'
208 : write(*,*) 'ANOMALY DETECTED'
209 : logan = .true.
210 : goto 9999
211 : endif
212 : T_b = P_b / (Rgas_b * r_b)
213 : C write(*,*) 'r_b, p_b, u_b, T_b, Dsw, k0new'
214 : C write(*,'(6E12.4)') r_b, p_b, u_b, T_b, Dsw
215 : C write(*,*) r_b, p_b, u_b, T_b, Dsw, k0new
216 : C stop
217 : CC write(*,*) 'We exit funre1'
218 : CC
219 : 9999 continue
220 : return
221 : end
222 :
© Cast3M 2003 - All rights reserved.
Disclaimer