1 : C FUNDT SOURCE BECC 09/12/07 21:15:30 6579
2 : subroutine fundt(nordpo, Tmaxcv, acv_r, Rgas_r,
3 : & acv_b, Rgas_b, q, T_ag,
4 : & P_r, T_r, u_r,
5 : & T_b, P_b, rho_b, u_b, d, k0sdt,
6 : & T_a, P_a, rho_a, u_a,
7 : & logan)
8 : C
9 : C INPUT
10 : C
11 : C nordpo = order of polynomial for cp and cv (see also
12 : C Tmaxcv)
13 : C
14 : C Tmaxcv = maximum temperature for cv polynomial expansion
15 : C cv(T) = cv(Tmaxcv) if T > Tmaxcv
16 : C
17 : C acv_r = to compute cv and ether for r (and rs);
18 : C vector such that
19 : C cv = \sum_{i=1,nordpo+1} acv(i) T^{i-1}
20 : C ether = \sum_{i=1,nordpo+1} acv(i) T^{i} / (i)
21 : C
22 : C Rgas_r = gas constant for r and (rs)
23 : C
24 : C acv_b = to compute cv and ether for rss;
25 : C vector such that
26 : C cv = \sum_{i=1,nordpo+1} acv(i) T^{i-1}
27 : C ether = \sum_{i=1,nordpo+1} acv(i) T^{i} / (i)
28 : C
29 : C Rgas_b = gas constant for rss
30 : C
31 : C q = released chemical energy
32 : C
33 : C T_ag = guess temperature for a (= rs, state ahead the
34 : C reactive shock and behind the non-reactive one
35 : C in the ZND model)
36 : C
37 : C T_b = temperature for rss
38 : C
39 : C
40 : C OUTPUT
41 : C
42 : C p_b, u_b, rho_b = pressure, temperature, velocity and density for
43 : C rss (burnt state behind the SDT)
44 : C
45 : C d = detonation speed
46 : C
47 : C k0sdt = fundamental speed of a WDF going from a to b
48 : C
49 : C p_a, T_a, u_a,
50 : C rho_a = pressure, temperature, velocity and density for
51 : C rs (unburnt state ahead the SDT)
52 : C
53 : c implicit none
54 : integer nordpo
55 : real*8 eps
56 : parameter (eps = 1.0D-12)
57 : real*8 Tmaxcv, acv_r(1:(nordpo+1)), Rgas_r
58 : & , acv_b(1:(nordpo+1)), Rgas_b, q, T_ag
59 : & , rho_r, P_r, T_r, u_r
60 : & , T_b, P_b, rho_b, u_b, d
61 : & , T_a, P_a, rho_a, u_a
62 : real*8 e_r, cv_r, h_r, cv_b, e_b, h_b, b, c, csi, c_b, b2, m, ksdt
63 : & , cv_ag, e_ag, gam_ag, k0sdt, qref
64 : logical logan, logdeb
65 : parameter (logdeb = .false.)
66 : C
66 : call prith1(nordpo, acv_r, Tmaxcv, T_r, e_r, cv_r)
67 : h_r = e_r + (Rgas_r * T_r)
67 : call prith1(nordpo, acv_b, Tmaxcv, T_b, e_b, cv_b)
68 : h_b = e_b + (Rgas_b * T_b)
69 : rho_r = P_r / (Rgas_r * T_r)
70 : C
71 : qref = Rgas_r * T_r
72 : if ((q .lt. (eps * qref)) .and.
73 : & (abs(T_b - T_r) .lt. (eps * T_r)) .and.
74 : & (abs(Rgas_r - Rgas_b) .lt. (eps * Rgas_r))) then
75 : C
76 : C******* Constant state
77 : C
78 : rho_b = rho_r
79 : T_b = T_r
80 : P_b = P_r
81 : u_b = u_r
82 : k0sdt = (Rgas_r + cv_r) / cv_r
83 : k0sdt = k0sdt * Rgas_r * T_r
84 : k0sdt = sqrt(k0sdt)
85 : d = k0sdt + u_r
86 : C
87 : T_a = T_r
88 : P_a = P_R
89 : rho_a = rho_r
90 : u_a = u_r
91 : else
92 : b = h_b - (h_r + q)
93 : b = b + (0.5D0 * ((Rgas_r * T_r) - (Rgas_b * T_b)))
94 : c = Rgas_r * T_r
95 : c_b = Rgas_b * T_b
96 : b2 = b * b
97 : csi = b2 + (c * c_b)
98 : if (csi .lt. 0.0D0) then
99 : write(*,*) 'fundt.f'
100 : write(*,*) 'delta = ', csi
101 : write(*,*) 'We take delta = 0'
102 : csi = 0.0D0
103 : endif
104 : csi = sqrt(csi) + b
105 : if (csi .lt. 0.0D0) then
106 : write(*,*) 'fundt.f'
107 : write(*,*) 'negative density'
108 : logan = .true.
109 : goto 9999
110 : endif
111 : csi = csi / c_b
112 : rho_b = rho_r * csi
113 : P_b = Rgas_b * rho_b * T_b
114 : C
115 : m = (P_b - P_r)
116 : m = m / ((1.0D0 / rho_r) - (1.0D0 / rho_b))
117 : if (m .lt. 0.0D0) then
118 : write(*,*) 'm < 0', m
119 : m = 0.0D0
120 : endif
121 : m = sqrt( m )
122 : ksdt = m / rho_r
123 : d = u_r + ksdt
124 : u_b = d - (m / rho_b)
125 : C
126 : C State a (ahead the reactive shock and behind the non-reactive
127 : C shock in the ZND model)
128 : C
128 : call prith1(nordpo, acv_r, Tmaxcv, T_ag, e_ag, cv_ag)
129 : c = 2.0D0 * (h_r - (e_ag - (cv_ag * T_ag)))
130 : c = c / (ksdt * ksdt)
131 : c = c + 1.0D0
132 : gam_ag = (cv_ag + Rgas_r) / cv_ag
133 : c = ((gam_ag - 1.0D0) / (gam_ag + 1.0D0)) * c
134 : b = Rgas_r * T_r / (ksdt * ksdt)
135 : b = b + 1.0D0
136 : b = (b * gam_ag) / (gam_ag + 1.0D0)
137 : C
138 : csi = b * b - c
139 : if (csi .lt. 0.0D0) then
140 : if (logdeb) then
141 : write(*,*) 'fundt.f'
142 : write(*,*) 'delta < 0 . ', csi
143 : write(*,*) 'Anomaly detected'
144 : endif
145 : csi = 0.0D0
146 : C logan = .true.
147 : C goto 9999
148 : endif
149 : csi = sqrt (csi)
150 : csi = b - csi
151 : if (csi .le. 0.0D0) then
152 : write(*,*) 'fundt.f'
153 : write(*,*) 'Anomaly detected'
154 : write(*,*) 'csi = ', csi
155 : logan = .true.
156 : goto 9999
157 : endif
158 : rho_a = rho_r / csi
159 : P_a = P_r + ((m * m) * ((1.0D0 / rho_r) - (1.0D0 / rho_a)))
160 : T_a = P_a / (Rgas_r * rho_a)
161 : u_a = d - m / rho_a
162 : k0sdt = d - u_a
163 : C write(*,*) 'ksdt, k0sdt ', ksdt, k0sdt
164 : endif
165 : C
166 : 9999 continue
167 : return
168 : end
169 :
© Cast3M 2003 - All rights reserved.
Disclaimer