1 : C FUNREA SOURCE BECC 09/12/07 21:15:32 6579
2 : subroutine funrea(nordp1, tmaxcv, rg_rss, acv_z,
3 : & rho_rs, p_rs, u_rs, et_rs, q, k0, t_gues,
4 : & rho_z, p_rss, u_rss, Dsw, et_rss, logan)
5 : c************************************************************************
6 : c
7 : c projet :
8 : c
9 : c name : funrea
10 : c
11 : c description : computation of the weak state behind a reactive
12 : c shock
13 : c
14 : c language : fortran 77
15 : c
16 : c author : a. beccantini den/dm2s/sfme/ltmf
17 : c
18 : c************************************************************************
19 : c
19 : c callen by : fcolre
20 : c
20 : c called program : none
21 : c
22 : c************************************************************************
23 : c
24 : c
25 : c**** input
26 : c
27 : c nordp1 = polynomials degree + 1
28 : c tmaxcv = threshold temperature for cv(t)
29 : c rg_rss = r of the gas in rss
30 : c acv_z = \sum_{i=1,nesp} y_i*acv_z{i,j}
31 : C rho_rs = density in rs
32 : C p_rs = pressure in rs
33 : c u_rs = velocity in rs
34 : c et_rs = (sensible) internal energy in rs
35 : c q = realeased heat per mass unit
36 : c k0 = fundamental flame speed
37 : c t_gues = guess temperature in rss around which we
38 : c linearised the internal energy (enthalpy)
39 : c**** output
40 : c rho_z = approximate density in rss (behind the shock)
41 : c p_rss = approximate pressure in rss
42 : c u_rss = approximate speed in rss
43 : C Dsw = approximate shock wave speed... it could
44 : C differs from u_rs + k0 if (k0lim < k0)
45 : c cv_rss = approximate cv in rss
46 : c et_rss = approximate (sensible) internal energy in rs
47 : c logan = an anomaly is detected
48 : c
49 : c************************************************************************
50 : c
51 : c created the 07/12/2009
52 : C
53 : C implicit none
54 : integer nordp1
55 : integer i1
56 : real*8 tmaxcv, rg_rss, acv_z(1:nordp1)
57 : real*8 rho_rs, p_rs, u_rs, et_rs, q, k0, t_gues
58 : real*8 rho_z, p_rss, u_rss, cv_rss
59 : real*8 et_rss, funt, funcv, t2, ga_rss
60 : real*8 acoef, bcoef, ccoef, acoef1, bcoef1, ccoef1, k0lim2, k0lim
61 : & ,Dsw, ratio
62 : logical logan
63 : CC
64 : CC**** We write the input
65 : CC
66 : C write(*,*) 'We enter funrea'
67 : C write(*,'(I2,E12.4)') nordp1, tmaxcv
68 : C write(*,'(6E12.4)') rg_rss, (acv_z(icon), icon=1,nordp1)
69 : C write(*,'(7E12.4)') rho_rs, p_rs, u_rs, et_rs, q, k0, t_gues
70 : C
71 : C**** Computation of ga_rss, ether_ss
72 : C
73 : C t2 = min (t_gues, Tmaxcv)
74 : t2 = 0.5d0 * ((t_gues + tmaxcv) - abs (t_gues - tmaxcv))
75 : cv_rss = acv_z(1)
76 : et_rss = t2 * cv_rss
77 : funt = 1.0d0
78 : do i1 = 2, nordp1, 1
79 : funt = funt * t2
80 : funcv = acv_z(i1) * funt
81 : cv_rss = cv_rss + funcv
82 : et_rss = et_rss + (funcv * t2/ i1)
83 : enddo
84 : c t2 = max ((t_gues - tmaxcv),0)
85 : t2 = 0.5d0 * ((t_gues - tmaxcv) + abs (t_gues - tmaxcv))
86 : et_rss = et_rss + (cv_rss * t2)
87 : ga_rss = (cv_rss + Rg_rss) / cv_rss
88 : C
89 : C**** Computation of rho_z
90 : C
91 : k0lim2 = k0 * k0
92 : bcoef = ga_rss / (ga_rss + 1.0d0)
93 : bcoef = bcoef * (1.0d0 + (p_rs / (rho_rs * k0lim2)))
94 : C
95 : ccoef = (p_rs / rho_rs) + et_rs + q
96 : ccoef = ccoef - (et_rss - (cv_rss * t_gues))
97 : ccoef = (ccoef * 2.0D0 / k0lim2) + 1.0d0
98 : ccoef = ccoef * (ga_rss - 1.0d0) / (ga_rss + 1.0d0)
99 : acoef = (bcoef * bcoef) - ccoef
100 : if (acoef .lt. 0.0d0) then
101 : C
102 : C******* k0 is too big.
103 : C Computation of K0lim2 = k0lim * k0lim corresponding to
104 : C CJ deflagration
105 : C
106 : acoef1 = (ga_rss * (p_rs / rho_rs)) / (ga_rss + 1.0d0)
107 : acoef1 = acoef1 * acoef1
108 : bcoef1 = ga_rss / (ga_rss + 1.0d0)
109 : bcoef1 = bcoef1 * bcoef1 * (p_rs / rho_rs)
110 : ccoef1 = (p_rs / rho_rs) + et_rs + q
111 : ccoef1 = ccoef1 - (et_rss - (cv_rss * t_gues))
112 : ccoef1 = ccoef1 * (ga_rss - 1.0d0) / (ga_rss + 1.0d0)
113 : bcoef1 = ccoef1 - bcoef1
114 : ccoef1 = 1.0d0 / ((ga_rss + 1.0d0) ** 2)
115 : k0lim2 = (bcoef1 * bcoef1) - (acoef1 * ccoef1)
116 : if ( (k0lim2 .lt. 0.0d0) ) then
117 : C
118 : C Discrimant zero... quite bizarre for exothermic reaction...
119 : C it can happen but for numerical reasons...
120 : C in that case, the best is to force it to zero!!!
121 : C
122 : logan = .true.
123 : write(*,'(/A8)') 'funrea.f'
124 : write(*,*) 'Negative discrimant for k0lim'
125 : write(*,*) 'delta ', k0lim2
126 : write(*,*) 'ANOMALY DETECTED'
127 : goto 9999
128 : endif
129 : k0lim2 = k0lim2 ** 0.5d0
130 : C Before this line, k0lim2 is just the square root of the discrimant
131 : k0lim2 = (bcoef1 + k0lim2) / acoef1
132 : k0lim2 = 1.0d0 / k0lim2
133 : bcoef1 = ga_rss / (ga_rss + 1.0d0)
134 : bcoef1 = bcoef1 * (1.0d0 + (p_rs / (rho_rs * k0lim2)))
135 : ratio = bcoef1
136 : else
137 : acoef1 = acoef ** 0.5d0
138 : ratio = bcoef - acoef1
139 : endif
140 : C
141 : rho_z = rho_rs / ratio
142 : if ( (rho_z .lt. 0.0d0) ) then
143 : C
144 : C Negative density
145 : C
146 : logan = .true.
147 : write(*,'(/A8)') 'funrea.f'
148 : C write(*,*) 'rho_rs, p_rs, u_rs, et_rs, q, k0, t_gues'
149 : C write(*,'(7E16.6)') rho_rs, p_rs, u_rs, et_rs, q, k0, t_gues
150 : C write(*,*) 'rho_z, p_rss, u_rss, Dsw, et_rss'
151 : C write(*,'(5E16.6)') rho_z, p_rss, u_rss, Dsw, et_rss
152 : write(*,*) 'Negative density'
153 : write(*,*) 'rho_z', rho_z
154 : write(*,*) 'ANOMALY DETECTED'
155 : goto 9999
156 : endif
157 : k0lim = k0lim2 ** 0.5d0
158 : p_rss = p_rs + (rho_rs * k0lim2 * (1.0d0 - ratio))
159 : Dsw = k0lim + u_rs
160 : u_rss = u_rs + ((p_rss - p_rs) / (k0lim * rho_rs))
161 : CC
162 : C write(*,*) 'rho_z, p_rss, u_rss, Dsw, et_rss'
163 : C write(*,'(5E12.4)') rho_z, p_rss, u_rss, Dsw, et_rss
164 : C write(*,*) 'We exit funrea'
165 : C
166 : 9999 continue
167 : return
168 : end
169 :
© Cast3M 2003 - All rights reserved.
Disclaimer