1 : C FVLHRE SOURCE BECC 11/01/05 21:15:05 6836
2 : subroutine fvlhre(ndim, nesp, nLHS, nordpo,
3 : & reacoe, aprop, Runiv,
4 : & Tmaxcv,
5 : & cnna, ctna, ct1na,
6 : & wvec_l,
7 : & wvec_r,
8 : & wwork, y, acv_l, acv_r, acv_b, acv_w ,
9 : & flu, f_lag, cellt,
10 : & logan)
11 : C
12 : c************************************************************************
13 : c
14 : c name : fvlhre
15 : c
16 : c description : Reactive Euler equations for a mixture of reactive
17 : C thermally perfect gases.
18 : c
19 : c VLH, implemented in the "Discrete Equation
20 : C Method" for the combustion
21 : C
22 : c
23 : c Voir:
24 : c 1) Beccantini,Studer, 2010, IJNMF
25 : C 2) Metayer, Massoni, Saurel, "Modelling
26 : C evaporations fronts with reactive Riemann
27 : C solvers.", JCP 205, 2005
28 : c
29 : c language : fortran 77
30 : c
31 : c author : a. beccantini den/dm2s/sfme/ltmf
32 : c
33 : c************************************************************************
34 : c
34 : c called sub : prithe, racre2, flufun
35 : C
36 : c************************************************************************
37 : c
38 : c**** input
39 : C
40 : C ndim = dimension of the domain
41 : c
42 : c nesp = species involved in the Euler equations
43 : C
44 : C nLHS = species involved in the LHS of the global
45 : C reaction
46 : c
47 : c nordp0 = polynomial order of the cv
48 : c
49 : c reacoe = coefficient involved in the global reaction
50 : c (positive in the LHS, negative in the RHS...
51 : c the first one should be positive)
52 : c
53 : c aprop = gas properties (cv, molar mass, formation
54 : C energy)
55 : c
56 : C Runiv = universal gas constant
57 : c
58 : C Tmaxcv = Tmax for cv expansion
59 : C
60 : C cnna, ctna, ct1na = scalar products
61 : C (\vec{n},\vec{n}_{\rm alpha})
62 : C (\vec{t},\vec{n}_{\rm alpha})
63 : C (\vec{t1},\vec{n}_{\rm alpha})
64 : C HP. We suppose that \vec{n}_{\rm alpha} is
65 : C always from the burnt to the unburnt state
66 : C
67 : c wvec_l = left primitive variables
68 : c
69 : c wvec_r = right primitive variables
70 : C
71 : C NB Primitive variables =
72 : C
73 : C w(1) = rho
74 : C w(2) = ux
75 : C w(3) = uy
76 : C w(4) = uz
77 : C w(2+ndim) = p
78 : C w(3+ndim) = csi
79 : C w(3+ndim+1) = yf(1)
80 : C ...
81 : C w(3+ndim+nLHS) = yf(nLHS)
82 : C w(3+ndim+nLHS+1) = yi(nLHS+1)
83 : C ...
84 : C w(3+ndim+nesp-1) = yi(nesp-1)
85 : C w(3+ndim+nesp) = dy1t = yi(1) - yf(1)
86 : C w(4+ndim+nesp) = k0
87 : C
88 : C For ndim = 1, the total number of variables
89 : C is 4 + 1 + nesp = 5 + nesp.
90 : c
91 : c wwork, y, acv_l,
92 : C acv_r, acv_b,
93 : C acv_w = temporary work vectors
94 : C wwork(4+ndim+nesp), y(1:nesp),
95 : C acv_l(0:nordpo),
96 : C acv_r(0:nordpo),
97 : C acv_b(0:nordpo),
98 : C acv_w (0:nordpo)
99 : C
100 : c**** output
101 : c
102 : c flu = Eulerian interfacial flux in (n,t1,t2), i.e.
103 : c rho*un mass
104 : C rho*un*un + p momentum along n
105 : C rho*un*ut1 momentum along t1
106 : * rho*un*ut2 momentum along t2
107 : C rho*un*ht energy
108 : c 0.0 alpha
109 : C rho*un*y_{1,f} mass * y_{1,f}
110 : C ...
111 : C rho*un*y_{nLHS+1,i} mass * y_{nLHS+1,i}
112 : c ...
113 : C rho*un*dy_{1,t} mass * dy_{1,t}
114 : C rho*un*k0 mass * k0
115 : C
116 : C NB. flu(1:(4+ndim+nesp))
117 : c
118 : C f_lag = Lagrangian interfacial flux on the reactive
119 : C discontinuty in (n,t1,t2), i.e.
120 : c rho*(un - D_resh) mass
121 : C rho*un*un + p - (D_resh rho un) momentum along n
122 : C rho*(un - D_resh)*ut1 momentum along t1
123 : C rho*(un - D_resh)*ut2 momentum along t2
124 : C rho*un*ht - (D_resh rho et) energy
125 : c - D_resh alpha
126 : c rho*(un - D_resh)*y_{1,f} mass * y_{1,f}
127 : C ...
128 : C
129 : C NB. f_lag(1:(4+ndim+nesp))
130 : C
131 : c
132 : c cellt = stability condition, i.e.
133 : c
134 : c dt/dx < cellt (dimension 1/velocity)
135 : c
136 : c logan = si true, anomaly detected
137 : c
138 : c************************************************************************
139 : c
140 : C 19/05/2010 created
141 : C
142 : c************************************************************************
143 : c
144 : c n.b.: all variables are declared
145 : c
146 : c implicit none
147 : implicit integer(i-n)
148 : integer ndim, nesp, nLHS, nordpo
149 : integer nesp1, nsca1
150 : integer icomp
151 : real*8 reacoe(nesp), aprop(0:(nordpo+2),nesp), Runiv, Tmaxcv
152 : & , wvec_l(4+ndim+nesp), wvec_r(4+ndim+nesp)
153 : & , flu(4+ndim+nesp)
154 : & , f_lag(4+ndim+nesp), cellt
155 : & , cnna, ctna, ct1na
156 : real*8 wwork(4+ndim+nesp)
157 : & , acv_l(0:nordpo), acv_r(0:nordpo), acv_b(0:nordpo)
158 : & , acv_w(0:nordpo)
159 : $ , y(1:nesp)
160 : real*8 gam_l, Rgas_l, cv_l, hf_l, rho_l, u_l(1:3), p_l, T_l
161 : & , et_l, csi_l
162 : & , gam_r, Rgas_r, cv_r, hf_r, rho_r, u_r(1:3), p_r, T_r
163 : & , et_r, csi_r
164 : & , un_r, un_l, ut_l, ut_r, ut1_l, ut1_r
165 : real*8 xst, flurho, fluru, fluret, coefy, omcoey
166 : C
167 : logical logan
168 : logical logdeb, logk0n
169 : C logdeb = debugging ?
170 : C logk0n = .true. : instead of k0, we take k0 * abs(cnna)
171 : C .false. : instead of k0, we take k0 * abs(cnna)
172 : parameter (logdeb = .false., logk0n = .true.)
173 : C
174 : C************************************************************************
175 : C******* NON REACTIVE CASE ONLY *****************************************
176 : C************************************************************************
177 : C
177 : call prithe(ndim, nesp, nLHS, nordpo, aprop, reacoe, Runiv,
178 : $ Tmaxcv,wvec_l,gam_l, Rgas_l, cv_l, acv_l, hf_l, rho_l,
179 : $ u_l, p_l, T_l,y, et_l, csi_l, logan)
180 : if (logan) then
181 : write(*,'(/A8)') 'fvlhre.f'
182 : write(*,*) 'ANOMALY DETECTED 01'
183 : goto 9999
184 : endif
185 : C
185 : call prithe(ndim, nesp, nLHS, nordpo, aprop, reacoe, Runiv,
186 : $ Tmaxcv,wvec_r,gam_r, Rgas_r, cv_r, acv_r, hf_r, rho_r,
187 : $ u_r, p_r, T_r,y, et_r, csi_r, logan)
188 : if (logan) then
189 : write(*,'(/A8)') 'fvlhre.f'
190 : write(*,*) 'ANOMALY DETECTED 02'
191 : goto 9999
192 : endif
193 : C
194 : un_l = u_l(1)
195 : un_r = u_r(1)
196 : C
197 : C d_l d_b un_b=un_a d_a d_r
198 : C
199 : xst = 0.0D0
200 : C
201 : nesp1 = 0
202 : nsca1 = 0
203 : C
204 : C We should take into account the formation energy in the flux
205 : C
206 : if (ndim .eq. 1) then
207 : ut_l = 0.0D0
208 : ut_r = 0.0D0
209 : ut1_l = 0.0D0
210 : ut1_r = 0.0D0
211 : elseif(ndim .eq. 2) then
212 : ut_l = wvec_l(3)
213 : ut_r = wvec_r(3)
214 : ut1_l = 0.0D0
215 : ut1_r = 0.0D0
216 : elseif(ndim .eq. 3) then
217 : ut_l = wvec_l(3)
218 : ut_r = wvec_r(3)
219 : ut1_l = wvec_l(4)
220 : ut1_r = wvec_r(4)
221 : endif
222 : C
222 : call fvlhtf (nesp1, nsca1,
223 : & gam_l, rho_l, p_l, un_l, ut_l, ut1_l, et_l, hf_l,
224 : & gam_r, rho_r, p_r, un_r, ut_r, ut1_r, et_r, hf_r,
225 : & 1.0D0, 1.0D0, 1.0D0, 1.0D0, flu, wwork,
226 : & cellt)
227 : C
228 : C flu(1) = flurho
229 : C flu(2) = fluru
230 : flurho = flu(1)
231 : fluru = flu(2)
232 : fluret = flu(4)
233 : flu(2+ndim) = fluret
234 : flu(3+ndim) = 0.0d0
235 : C The Larrouturou theorem does hold
236 : coefy = 0.5d0 * (1.0d0 + sign(1.0d0, flurho))
237 : omcoey = 1.0d0 - coefy
238 : C Tangential speeds (ndim .ge. 2)
239 : do icomp = 3, (ndim + 1), 1
240 : flu(icomp) = flurho * ((coefy * wvec_l(icomp)) +
241 : & (omcoey * wvec_r(icomp)))
242 : enddo
243 : do icomp = (ndim + 4), (ndim + nesp + 4), 1
244 : flu(icomp) = flurho * ((coefy * wvec_l(icomp)) +
245 : & (omcoey * wvec_r(icomp)))
246 : enddo
247 : C
248 : C******* f_lag is here not defined !!!
249 : C
250 : 9999 continue
251 : return
252 : end
253 :
254 :
© Cast3M 2003 - All rights reserved.
Disclaimer