1 : C PRITHE SOURCE BECC 10/09/10 21:15:01 6472
2 : subroutine prithe(ndim,nesp,nLHS,nordpo,aprop,reacoe,Runiv,Tmaxcv
3 : $ ,w,gamma,Rgas,cvtot,acvtot,hftot,rho,u,p,T,y,ether,csi,logan)
4 : C
5 : C This subroutine creates some thermodynamical variables and
6 : C properties starting from the primitive variables.
7 : C
8 : C
9 : C
10 : C INPUT
11 : C *****
12 : C
13 : C ndim = dimension (1,2,3)
14 : C
15 : C nesp = number of species in the Euler equations
16 : C
17 : C nLHS = number of species in the LHS of the chemical reaction
18 : C nLHS < nesp
19 : C
20 : C nordpo = order of polynomial for cp and cv (see also Tmaxcv)
21 : C
22 : C aprop = gas properties: ,reacoe,Runiv,
23 : C aprop(0:nordpo+2,iesp) :
24 : C (coefcv(i),i=0,nordpo,1), Wi, hf0i)
25 : C
26 : C reacoe(1:nesp) = coefficient of each species in the global chemical
27 : C reaction
28 : C
29 : C Runiv = universal constant of the gas
30 : C
31 : C Tmaxcv = maximum temperature for cv polynomial expansion
32 : C cv(T) = cv(Tmaxcv) if T > Tmaxcv
33 : C
34 : C w = primitive variables
35 : C w(1) = rho
36 : C w(2) = ux
37 : C w(3) = uy
38 : C w(4) = uz
39 : C w(2+ndim) = p
40 : C w(3+ndim) = csi
41 : C w(3+ndim+1) = yf(1)
42 : C ...
43 : C w(3+ndim+nLHS) = yf(nLHS)
44 : C w(3+ndim+nLHS+1) = yi(nLHS+1)
45 : C ...
46 : C w(3+ndim+nesp-1) = yi(nesp-1)
47 : C w(3+ndim+nesp) = dy1t = yi(1) - yf(1)
48 : C w(4+ndim+nesp) = k0
49 : C
50 : C
51 : C
52 : C OUTPUT
53 : C ******
54 : C
55 : C gamma = specific heat ratio
56 : C
57 : C Rgas = gas constant
58 : C
59 : C cvtot = cv
60 : C
61 : C acvtot = coefficients such that for the mixture
62 : C cv = \sum_{i=1,nordpo+1} acvtot(i) T^{i-1}
63 : C ether = \sum_{i=1,nordpo+1} acvtot(i) T^{i} / (i)
64 : C
65 : C hftot = \sum_i Y_i hf0i = formation energy of the mixture at 0K
66 : C
67 : C rho = density
68 : C
69 : C u(1:ndim) = velocity
70 : C
71 : C p = pressure
72 : C
73 : C T = temperature
74 : C
75 : C y = mass fractions (1:nesp), \sum y(i) = 1
76 : C
77 : C ether = sensible energy at 0K
78 : C
79 : C csi = progress variable
80 : C
81 : C
82 : C ERROR VARIABLE
83 : C **************
84 : C
85 : C logan = if it becomes .false. an anomaly has been detected.
86 : C
87 : C************************************************************************
88 : C
89 : C HISTORY AND MODIFICATIONS
90 : C
91 : c 07/12/2009 Created.
92 : C
93 : C 30/07/2010 We check that species mass fractions are not lower
94 : C than yneg > 0 (and not zero).
95 : C
96 : C************************************************************************
97 : C
98 : C implicit none
99 : integer ndim, nesp, nLHS, nordpo
100 : & , iesp, icoef, idim
101 : real*8 aprop(0:(nordpo+2),nesp), reacoe(nesp), Runiv, Tmaxcv
102 : & , w(1:(4+ndim+nesp))
103 : & , gamma, Rgas, cvtot, acvtot(1:(nordpo+1)), hftot
104 : & , rho, u(1:ndim), p, T, y(1:nesp), ether
105 : & , csi, funt, dcv, coef, coef1, T1, yneg
106 : parameter (yneg = -1.0D-4)
107 : logical logan
108 : C
109 : C Initialisation of some variable
110 : C
111 : y(nesp)=1.0d0
112 : hftot = 0.0d0
113 : do icoef = 1 ,(nordpo+1), 1
114 : acvtot(icoef) = 0.0d0
115 : enddo
116 : Rgas = 0.0d0
117 : C
118 : C Density, velociy, pressure
119 : C
120 : rho = w(1)
121 : do idim = 1, ndim, 1
122 : u(idim) = w(1+idim)
123 : enddo
124 : p = w(2+ndim)
125 : C
126 : C Csi and mass fractions
127 : C
128 : csi = w(3+ndim)
129 : if ((csi . lt. 0.0d0) .or. (csi .gt. 1.0d0)) then
130 : write(*,*) 'prithe.f'
131 : write(*,*) 'csi = ', csi
132 : write(*,*) 'Negative csi.'
133 : logan = .true.
134 : goto 9999
135 : endif
136 : C coef = dy1t / (w_1 * a_1)
137 : coef = w(3+ndim+nesp) / (aprop(nordpo+1,1) * reacoe(1))
138 : do iesp = 1, nLHS, 1
139 : coef1 = aprop(nordpo+1,iesp) * reacoe(iesp)
140 : if (coef1 .le. 0.0D0) then
141 : write(*,*) 'prithe.f'
142 : write(*,*) 'Incompatibility in reacoe variable for ', iesp
143 : write(*,*) 'reacoe should be > 0. reacoe =', reacoe(iesp)
144 : logan=.true.
145 : goto 9999
146 : endif
147 : C y_j = y_{j,f} + ((1 - csi) * (y_{j,i} - y_{j,f}))
148 : y(iesp) = w(3+ndim+iesp) + ((1.0d0 - csi) * coef1 * coef)
149 : if(y(iesp) .lt. yneg)then
150 : write(*,*) 'prithe.f'
151 : write(*,*) 'Negative mass fraction, species ', iesp
152 : write(*,*) y(iesp)
153 : logan=.true.
154 : goto 9999
155 : endif
156 : y(nesp) = y(nesp) - y(iesp)
157 : enddo
158 : C
159 : do iesp = (nLHS + 1), (nesp - 1), 1
160 : coef1 = aprop(nordpo+1,iesp) * reacoe(iesp)
161 : if (coef1 .gt. 0.0D0) then
162 : write(*,*) 'prithe.f'
163 : write(*,*) 'Incompatibility in reacoe variable for ', iesp
164 : write(*,*) 'reacoe should be <= 0. reacoe =', reacoe(iesp)
165 : logan=.true.
166 : goto 9999
167 : endif
168 : C y_j = y_{j,i} + csi * (y_{j,f} - y_{j,i})
169 : y(iesp) = w(3+ndim+iesp) - (csi * coef1 * coef)
170 : if(y(iesp) .lt. yneg)then
171 : write(*,*) 'prithe.f'
172 : write(*,*) 'Negative mass fraction, species ', iesp
173 : write(*,*) y(iesp)
174 : logan=.true.
175 : goto 9999
176 : endif
177 : y(nesp) = y(nesp) - y(iesp)
178 : enddo
179 : C
180 : if(y(nesp) .lt. yneg)then
181 : write(*,*) 'prithe.f'
182 : write(*,*) 'Negative mass fraction, species ', nesp
183 : write(*,*) y(nesp)
184 : logan=.true.
185 : goto 9999
186 : endif
187 : C
188 : C Rgas, hftot, Temperature and acvtot
189 : C In order to avoid one more loop, this could have be done before
190 : C
191 : do iesp=1,nesp,1
192 : C Rgas = y(iesp) / W_iesp
193 : Rgas = Rgas + (y(iesp)/aprop(nordpo+1,iesp))
194 : C hftot = hftot + (y_i * hf0i)
195 : hftot = hftot + (y(iesp) * aprop(nordpo+2,iesp))
196 : do icoef = 0, nordpo, 1
197 : acvtot(icoef+1) = acvtot(icoef+1) +
198 : & (y(iesp) * aprop(icoef,iesp))
199 : enddo
200 : enddo
201 : Rgas = Rgas * Runiv
202 : T = p / (Rgas * rho)
203 : C
204 : C gamma, ether
205 : C
206 : C T1 = min (T, Tmaxcv)
207 : T1 = 0.5d0 * ((T + Tmaxcv) - abs (T - Tmaxcv))
208 : C
209 : funt = 1.0D0
210 : cvtot = acvtot(1)
211 : ether = cvtot * T1
212 : do icoef = 2, nordpo+1, 1
213 : funt = funt * T1
214 : dcv = acvtot(icoef) * funt
215 : cvtot = cvtot + dcv
216 : ether = ether + (dcv * T1 / icoef)
217 : enddo
218 : gamma = (Rgas + cvtot) / cvtot
219 : C
220 : C if ( T > Tmaxcv ) => T1 = (T - Tmaxcv)
221 : C else T1 = 0
222 : C i.e.
223 : C T1 = max ((T - Tmaxcv),0)
224 : C
225 : T1 = 0.5d0 * ((T - Tmaxcv) + abs (T - Tmaxcv))
226 : ether = ether + (cvtot * T1)
227 : C
228 : 9999 continue
229 : return
230 : end
231 : C
232 :
233 :
© Cast3M 2003 - All rights reserved.
Disclaimer