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