C FVLHRE SOURCE PV 22/04/26 21:15:03 11344 subroutine fvlhre(ndim, nesp, nLHS, nordpo, & reacoe, aprop, Runiv, & Tmaxcv, & cnna, ctna, ct1na, & wvec_l, & wvec_r, & wwork, y, acv_l, acv_r, acv_b, acv_w , & flu, f_lag, cellt, & logan) C c************************************************************************ c c name : fvlhre c c description : Reactive Euler equations for a mixture of reactive C thermally perfect gases. c c VLH, implemented in the "Discrete Equation C Method" for the combustion C c c Voir: c 1) Beccantini,Studer, 2010, IJNMF C 2) Metayer, Massoni, Saurel, "Modelling C evaporations fronts with reactive Riemann C solvers.", JCP 205, 2005 c c language : fortran 77 c c author : a. beccantini den/dm2s/sfme/ltmf c c************************************************************************ c c called sub : prithe, racre2, flufun C c************************************************************************ c c**** input C C ndim = dimension of the domain c c nesp = species involved in the Euler equations C C nLHS = species involved in the LHS of the global C reaction c c nordp0 = polynomial order of the cv c c reacoe = coefficient involved in the global reaction c (positive in the LHS, negative in the RHS... c the first one should be positive) c c aprop = gas properties (cv, molar mass, formation C energy) c C Runiv = universal gas constant c C Tmaxcv = Tmax for cv expansion C C cnna, ctna, ct1na = scalar products C (\vec{n},\vec{n}_{\rm alpha}) C (\vec{t},\vec{n}_{\rm alpha}) C (\vec{t1},\vec{n}_{\rm alpha}) C HP. We suppose that \vec{n}_{\rm alpha} is C always from the burnt to the unburnt state C c wvec_l = left primitive variables c c wvec_r = right primitive variables C C NB Primitive variables = C C w(1) = rho C w(2) = ux C w(3) = uy C w(4) = uz C w(2+ndim) = p C w(3+ndim) = csi C w(3+ndim+1) = yf(1) C ... C w(3+ndim+nLHS) = yf(nLHS) C w(3+ndim+nLHS+1) = yi(nLHS+1) C ... C w(3+ndim+nesp-1) = yi(nesp-1) C w(3+ndim+nesp) = dy1t = yi(1) - yf(1) C w(4+ndim+nesp) = k0 C C For ndim = 1, the total number of variables C is 4 + 1 + nesp = 5 + nesp. c c wwork, y, acv_l, C acv_r, acv_b, C acv_w = temporary work vectors C wwork(4+ndim+nesp), y(1:nesp), C acv_l(0:nordpo), C acv_r(0:nordpo), C acv_b(0:nordpo), C acv_w (0:nordpo) C c**** output c c flu = Eulerian interfacial flux in (n,t1,t2), i.e. c rho*un mass C rho*un*un + p momentum along n C rho*un*ut1 momentum along t1 * rho*un*ut2 momentum along t2 C rho*un*ht energy c 0.0 alpha C rho*un*y_{1,f} mass * y_{1,f} C ... C rho*un*y_{nLHS+1,i} mass * y_{nLHS+1,i} c ... C rho*un*dy_{1,t} mass * dy_{1,t} C rho*un*k0 mass * k0 C C NB. flu(1:(4+ndim+nesp)) c C f_lag = Lagrangian interfacial flux on the reactive C discontinuty in (n,t1,t2), i.e. c rho*(un - D_resh) mass C rho*un*un + p - (D_resh rho un) momentum along n C rho*(un - D_resh)*ut1 momentum along t1 C rho*(un - D_resh)*ut2 momentum along t2 C rho*un*ht - (D_resh rho et) energy c - D_resh alpha c rho*(un - D_resh)*y_{1,f} mass * y_{1,f} C ... C C NB. f_lag(1:(4+ndim+nesp)) C c c cellt = stability condition, i.e. c c dt/dx < cellt (dimension 1/velocity) c c logan = si true, anomaly detected c c************************************************************************ c C 19/05/2010 created C c************************************************************************ c c n.b.: all variables are declared c c implicit none implicit integer(i-n) integer ndim, nesp, nLHS, nordpo integer nesp1, nsca1 integer icomp real*8 reacoe(nesp), aprop(0:(nordpo+2),nesp), Runiv, Tmaxcv & , wvec_l(4+ndim+nesp), wvec_r(4+ndim+nesp) & , flu(4+ndim+nesp) & , f_lag(4+ndim+nesp), cellt & , cnna, ctna, ct1na real*8 wwork(4+ndim+nesp) & , acv_l(0:nordpo), acv_r(0:nordpo), acv_b(0:nordpo) & , acv_w(0:nordpo) $ , y(1:nesp) real*8 gam_l, Rgas_l, cv_l, hf_l, rho_l, u_l(1:3), p_l, T_l & , et_l, csi_l & , gam_r, Rgas_r, cv_r, hf_r, rho_r, u_r(1:3), p_r, T_r & , et_r, csi_r & , un_r, un_l, ut_l, ut_r, ut1_l, ut1_r real*8 xst, flurho, fluru, fluret, coefy, omcoey real*8 void(1) C logical logan logical logdeb, logk0n C logdeb = debugging ? C logk0n = .true. : instead of k0, we take k0 * abs(cnna) C .false. : instead of k0, we take k0 * abs(cnna) parameter (logdeb = .false., logk0n = .true.) C void(1)=1.d0 C************************************************************************ C******* NON REACTIVE CASE ONLY ***************************************** C************************************************************************ C call prithe(ndim, nesp, nLHS, nordpo, aprop, reacoe, Runiv, $ Tmaxcv,wvec_l,gam_l, Rgas_l, cv_l, acv_l, hf_l, rho_l, $ u_l, p_l, T_l,y, et_l, csi_l, logan) if (logan) then write(*,'(/A8)') 'fvlhre.f' write(*,*) 'ANOMALY DETECTED 01' goto 9999 endif C call prithe(ndim, nesp, nLHS, nordpo, aprop, reacoe, Runiv, $ Tmaxcv,wvec_r,gam_r, Rgas_r, cv_r, acv_r, hf_r, rho_r, $ u_r, p_r, T_r,y, et_r, csi_r, logan) if (logan) then write(*,'(/A8)') 'fvlhre.f' write(*,*) 'ANOMALY DETECTED 02' goto 9999 endif C un_l = u_l(1) un_r = u_r(1) C C d_l d_b un_b=un_a d_a d_r C xst = 0.0D0 C nesp1 = 0 nsca1 = 0 C C We should take into account the formation energy in the flux C if (ndim .eq. 1) then ut_l = 0.0D0 ut_r = 0.0D0 ut1_l = 0.0D0 ut1_r = 0.0D0 elseif(ndim .eq. 2) then ut_l = wvec_l(3) ut_r = wvec_r(3) ut1_l = 0.0D0 ut1_r = 0.0D0 elseif(ndim .eq. 3) then ut_l = wvec_l(3) ut_r = wvec_r(3) ut1_l = wvec_l(4) ut1_r = wvec_r(4) endif C call fvlhtf (nesp1, nsca1, & gam_l, rho_l, p_l, un_l, ut_l, ut1_l, et_l, hf_l, & gam_r, rho_r, p_r, un_r, ut_r, ut1_r, et_r, hf_r, & void(1),void(1),void(1),void(1), flu, wwork, & cellt) C C flu(1) = flurho C flu(2) = fluru flurho = flu(1) fluru = flu(2) fluret = flu(4) flu(2+ndim) = fluret flu(3+ndim) = 0.0d0 C The Larrouturou theorem does hold coefy = 0.5d0 * (1.0d0 + sign(1.0d0, flurho)) omcoey = 1.0d0 - coefy C Tangential speeds (ndim .ge. 2) do icomp = 3, (ndim + 1), 1 flu(icomp) = flurho * ((coefy * wvec_l(icomp)) + & (omcoey * wvec_r(icomp))) enddo do icomp = (ndim + 4), (ndim + nesp + 4), 1 flu(icomp) = flurho * ((coefy * wvec_l(icomp)) + & (omcoey * wvec_r(icomp))) enddo C C******* f_lag is here not defined !!! C 9999 continue return end