1 : C FUNSDT SOURCE BECC 09/12/07 21:15:34 6579
2 : subroutine funsdt(nordpo, Tmaxcv, acv_r, Rgas_r,
3 : & acv_b, Rgas_b, q, T_ag, T_bg,
4 : & P_r, T_r, u_r,
5 : & T_a, P_a, r_a, u_a, d, k0sdt,
6 : & T_b, P_b, r_b, u_b,
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_bg = guess temperature for b (= rss, state behind the
38 : C reactive shock in the ZND model)
39 : C
40 : C T_a = temperature for rs
41 : C
42 : C
43 : C OUTPUT
44 : C
45 : C p_b, T_b, u_b,
46 : C r_b = pressure, temperature, velocity and density for
47 : C rss (burnt state behind the SDT)
48 : C
49 : C d = detonation speed
50 : C
51 : C k0sdt = fundamental speed of a WDF going from a to b
52 : C
53 : C p_a, u_a,
54 : C r_a = pressure, velocity and density for
55 : C rs (unburnt state ahead the SDT)
56 : C
57 : c implicit none
58 : integer nordpo
59 : real*8 Tmaxcv, acv_r(1:(nordpo+1)), Rgas_r
60 : & , acv_b(1:(nordpo+1)), Rgas_b, q, T_ag, T_bg
61 : & , r_r, P_r, T_r, u_r
62 : & , T_a, P_a, r_a, u_a, d
63 : & , T_b, P_b, r_b, u_b
64 : real*8 k0sdt, e_a, cv_a, h_a, gam_a
65 : & , gam_b, e_b, cv_b
66 : & , c, b, csi, m
67 : logical logan
68 : C
69 : C**** Unburnt state
70 : C
71 : r_r = P_r / (Rgas_r * T_r)
71 : call funsho( nordpo, Tmaxcv, acv_r, Rgas_r,
72 : & P_r, T_r, u_r, .true.,
73 : & T_a, P_a, r_a, u_a, d)
74 : C
75 : C**** Now we search for the burnt state which travels with the
76 : C same speed
77 : C
78 : k0sdt = d - u_a
78 : call prith1(nordpo, acv_r, Tmaxcv, T_ag, e_a, cv_a)
79 : gam_a = (cv_a + Rgas_r) / 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 : C
82 : h_a = (gam_a / (gam_a - 1.0D0)) * Rgas_r * T_a
83 : c = h_a + q - (e_b - (cv_b * T_bg)) +
84 : & (e_a - (cv_a * T_ag))
85 : c = 2.0D0 * c
86 : c = c / (k0sdt * k0sdt)
87 : c = c + 1.0D0
88 : c = c * (gam_b - 1.0D0) / (gam_b + 1.0D0)
89 : if (c .lt. 0.0D0) then
90 : write(*,*) 'funsdt.f'
91 : write(*,*) 'c < 0'
92 : write(*,*) 'Anomaly detected'
93 : logan = .true.
94 : goto 9999
95 : endif
96 : C
97 : b = Rgas_r * T_a
98 : b = b / (k0sdt * k0sdt)
99 : b = b + 1.0D0
100 : b = b * gam_b / (gam_b + 1.0D0)
101 : csi = b * b - c
102 : if (csi .lt. 0.0D0) then
103 : C Non itersection between the rayleigh line
104 : C and the Hugoniot adiabatic
105 : if (csi .lt. 1.0D-3) then
106 : C csi can be considered as zero
107 : csi = 0.0D0
108 : else
109 : write(*,*) 'csi = ', csi
110 : write(*,*) 'csi < 0'
111 : write(*,*) 'Anomaly detected'
112 : logan = .true.
113 : goto 9999
114 : endif
115 : else
116 : csi = sqrt(csi)
117 : endif
118 : C Computation of the SDT solution
119 : C For the WDT solution, csi = b + csi
120 : csi = b - csi
121 : r_b = r_a / csi
122 : m = k0sdt * r_a
123 : p_b = p_a + ((m*m) * ((1.0D0 / r_a) - (1.0D0 / r_b)))
124 : t_b = P_b / (Rgas_b * r_b)
125 : u_b = d - (m / r_b)
126 : C
127 : 9999 continue
128 : return
129 : end
130 :
© Cast3M 2003 - All rights reserved.
Disclaimer