C FCOLRE SOURCE BECC 09/12/07 21:15:17 6579 subroutine fcolre(ndim, nesp, nLHS, nordpo, & reacoe, aprop, Runiv, & Tmaxcv, & wvec_l, & wvec_r, & wwork, y, acv_l, acv_r, acv_z, acv_w , & flu, f_lag, cellt, & logan) C c************************************************************************ c c project : c c name : fcolre c c description : Reactive Euler equations for a mixture of reactive C thermally perfect gases. c c Shock-shock fds + entropy fix, implemented in a C "Discrete Equation Method" fashion. c c Voir: c 1) Beccantini, report ??? 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, racrea, flufun C C called by : comres 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 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 idim = 1, the total number of variables C is 4 + 1 + nesp = 5 + nesp. c c wwork, y, acv_l, C acv_r, acv_z, 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_z(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 07/12/2009 created c c************************************************************************ c c n.b.: all variables are declared c C implicit none integer ndim, nesp, nLHS, nordpo integer icomp real*8 zero, rind_l parameter (zero=1.0d-10) real*8 reacoe(nesp), aprop(0:nordpo+2,nesp), Runiv, Tmaxcv & , xst, wvec_l(4+ndim+nesp), wvec_r(4+ndim+nesp) & , flu(4+ndim+nesp) & , f_lag(4+ndim+nesp), cellt real*8 wwork(4+ndim+nesp) & , acv_l(0:nordpo), acv_r(0:nordpo), acv_z(0:nordpo) & , acv_w (0:nordpo) $ , y(1:nesp) real*8 csi0, csi1, q C Left real*8 gam_l, Rgas_l, cv_l, hf_l, rho_l, u_l(1:3), p_l, T_l $ , et_l, csi_l, un_l C Left star... same Rgas, acv, hf and csi as left real*8 rho_ls, p_ls, T_ls, d_ls, un_ls C Right star star real*8 ga_rss, Rg_rss, cv_rss, hf_rss, cs_rss, rh_rss $ , p_rss, et_rss, T_rss, un_rss C Right star... same Rgas, acv, hf and csi as right real*8 rho_rs, p_rs, d_rs, T_rs, un_rs C Right real*8 gam_r, Rgas_r, cv_r, hf_r, rho_r, u_r(1:3), p_r, T_R & , et_r & , csi_r, un_r real*8 c_l, c_r, c_ls, c_rs real*8 k0, D_resh real*8 vit_l, vit_ls, vit_rs, vit_r real*8 coefy, omcoey real*8 flurho, fluru, fluret, trekin C logical logan, lognc, logdeb parameter (logdeb = .false.) C C We have to distinguish between 5 different states C C l = burnt state (csi = csimax) C | Left GNL wave C ls = connected to l by a GNL wave (csi = csimax) C | Contact discontinuty C rsw = connected to rs by a reactive wave (csi = csimax) C | Reactive shock (weak deflagration wave) C rs = connected to r by a GNL wave (csi = csimin) C | Right GNL wave C r = unburnt state (csi = csimin) C C We define the left state such that csi_l > csi_r C csi0 = wvec_l(ndim+3) csi1 = wvec_r(ndim+3) if (csi0 .gt. csi1) then k0 = wvec_r(4+ndim+nesp) rind_l = 1.0D0 C C The right star star state is the burnt state relatively to the C right and right stare states. It has the same value of csi as C the left state. C Note that if this state has not necessary the same mass C fractions as the left state. C We compute the physical properties of the right star state, C namely C ga_rss, Rg_rss, cv_rss, acv_z, hf_rss. C The other physical values are of course meaningless! C do icomp = 1, (nesp + 4 + ndim), 1 wwork(icomp) = wvec_r(icomp) enddo wwork(3+ndim) = wvec_l(3+ndim) C call prithe(ndim, nesp, nLHS, nordpo, aprop, reacoe, Runiv, $ Tmaxcv,wwork,ga_rss, Rg_rss, cv_rss, acv_z, hf_rss, $ rh_rss,u_r, p_rss, T_rss ,y, et_rss, cs_rss, logan) if (logan) then write(*,'(/A8)') 'fcolre.f' write(*,*) 'ANOMALY DETECTED 3' goto 9999 endif C C The left state is considered as the burnt state. C We compute the physical properties of the left state, namely C gam_l, Rgas_l, cv_l, acv_l, hf_l. C We also compute the physical variable on the left, namely C rho_l, u_l, p_l, T_L, et_l and also csi_r 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)') 'fcolre.f' write(*,*) 'ANOMALY DETECTED 1' goto 9999 endif C C The right state is considered as the unburnt state. C We compute the physical properties of the right state, namely C gam_r, Rgas_r, cv_r, acv_r, hf_r. C We also compute the physical variable on the left, namely C rho_r, u_r, p_r, T_r, et_r, csi_r 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)') 'fcolre.f' write(*,*) 'ANOMALY DETECTED 2' goto 9999 endif else rind_l = -1.0D0 k0 = wvec_l(4+ndim+nesp) C do icomp = 1, (nesp + 4 + ndim), 1 wwork(icomp)=wvec_l(icomp) enddo wwork(3+ndim)=wvec_r(3+ndim) C call prithe(ndim, nesp, nLHS, nordpo, aprop, reacoe, Runiv, $ Tmaxcv,wwork,ga_rss, Rg_rss, cv_rss, acv_z, hf_rss, $ rh_rss,u_l, p_rss, T_rss, y, et_rss, cs_rss, logan) if (logan) then write(*,'(/A8)') 'fcolre.f' write(*,*) 'ANOMALY DETECTED 6' goto 9999 endif C call prithe(ndim, nesp, nLHS, nordpo, aprop, reacoe, Runiv, $ Tmaxcv,wvec_r,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)') 'fcolre.f' write(*,*) 'ANOMALY DETECTED 4' goto 9999 endif C un_l=-1*u_l(1) C call prithe(ndim, nesp, nLHS, nordpo, aprop, reacoe, Runiv, $ Tmaxcv,wvec_l,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)') 'fcolre.f' write(*,*) 'ANOMALY DETECTED 5' goto 9999 endif C un_r=-1*u_r(1) endif un_r = rind_l * u_r(1) un_l = rind_l * u_l(1) C C**** Test C if (logdeb ) then write(*,'(A23)') 'rho_r, p_r, t_r, u_r' write(*,'(4E16.6)') rho_r, p_r, t_r, un_r write(*,'(A23)') 'rho_l, p_l, t_l, u_l' write(*,'(4E16.6)') rho_l, p_l, t_l, un_l endif C C The reaction should be exothermic C q = hf_r - hf_rss if (q .lt. (-1.0D0 * zero * max( abs(hf_rss), abs(hf_r) )) ) then write(*,'(/A8)') 'fcolre.f' write(*,*) 'Endothermic reaction ??? ' write(*,*) 'q =', q logan = .true. goto 9999 endif q = 0.5d0 * (q + abs(q)) C C Computation of the intermediate states C C Vacuum (TO DO!!!) C C Intermediate states in the case of no-vacuum C C write(*,*) k0, q C write(*,*) rho_l, p_l, un_l, T_l, et_l, gam_l C write(*,*) rho_r, p_r, un_r, T_r, et_r, gam_r C write(*,*) rh_rss, p_rss, un_l, T_rss, et_rss, ga_rss C write(*,*) csi0, csi1 C call racrea(nordpo, Tmaxcv, & Rgas_l, acv_l, Rg_rss, acv_z, Rgas_r, acv_r, & q, k0, & rho_l, p_l, un_l, T_l, et_l, gam_l, & rho_r, p_r, un_r, T_r, et_r, gam_r, & rho_ls, p_ls, un_ls, T_ls, c_ls, d_ls, & rh_rss, p_rss, un_rss, T_rss, D_resh, & rho_rs, p_rs, un_rs, T_rs, c_rs, d_rs, & logan, lognc) if (lognc) then C C Maybe this is due to unphysical problems involved in this C approach C do icomp = 1, (4 + ndim + nesp), 1 flu(icomp) = 0.0d0 enddo C do icomp = 1, (4 + ndim + nesp), 1 f_lag(icomp) = 0.0d0 enddo goto 9999 endif C if (logan) then write(*,'(/A8)') 'fcolre.f' write(*,*) 'rho_l, p_l, un_l, T_l, et_l, gam_l' write(*,'(6E16.6)') rho_l, p_l, un_l, T_l, et_l, gam_l write(*,*) 'rho_r, p_r, un_r, T_r, et_r, gam_r' write(*,'(6E16.6)') rho_r, p_r, un_r, T_r, et_r, gam_r write(*,*) ' icomp, wvec_l, wvec_r' do icomp = 1, (4 + nesp + ndim) write(*,'(I6,2E16.6)') icomp, wvec_l(icomp), wvec_r(icomp) enddo write(*,*) 'ANOMALY DETECTED' goto 9999 endif C if (logdeb) then write(*,*) rho_l, p_l, un_l, T_l write(*,*) rho_ls, p_ls, un_ls, T_ls write(*,*) rh_rss, p_rss, un_rss, T_rss write(*,*) rho_rs, p_rs, un_rs, T_rs write(*,*) rho_r, p_r, un_r, T_r C stop endif C C**** Wave speed C C Speed of the GNL-waves (rarefaction or shock) C C C L LS C ****** ******* C * or * C * LS L* C ******* ******* C C On the rarefaction we lineairize over the variables C rho, T, w C c_l = gam_l * p_l / rho_l c_l = c_l ** 0.5d0 c_r = gam_r * p_r / rho_r c_r = c_r ** 0.5d0 C vit_l = min(d_ls, (un_l - c_l)) vit_l = 0.5d0*((d_ls + (un_l - c_l)) & - abs( d_ls - (un_l - c_l) )) C vit_ls = max(d_ls, (un_ls - c_ls)) vit_ls = 0.5d0*((d_ls + (un_ls - c_ls)) & + abs( d_ls - (un_ls - c_ls) )) C vit_ls exactly => vit_l vit_ls = max(vit_ls, vit_l) C vit_rs = min(d_rs, (u_rs + c_rs)) vit_rs = 0.5d0*((d_rs + (un_rs + c_rs)) & - abs( d_rs - (un_rs + c_rs) )) C vit_r = max(d_r, (un_r + c_r)) vit_r = 0.5d0*((d_rs + (un_r + c_r)) & + abs( d_rs - (un_r + c_r) )) C vit_rs exactly =< vit_r vit_rs = min(vit_rs, vit_r) C C write(*,*) vit_l, vit_ls, un_ls, un_rss, D_resh, vit_rs, vit_r C stop C C**** Maximum eigenvalue C cellt = abs(vit_r) cellt = 0.5d0 * (abs(vit_l) + cellt + abs ( abs(vit_l) - cellt)) C C**** Eulerian flux (1D) C xst = 0.0d0 call flufun(nordpo, & Tmaxcv, acv_w , & xst, & vit_l, vit_ls, un_ls, D_resh, vit_rs, vit_r, & Rgas_l, acv_l, rho_l, un_l, T_l, hf_l, & rho_ls, T_ls, & Rg_rss, acv_z, rh_rss, un_rss, T_rss, hf_rss, & rho_rs, un_rs, T_rs, & Rgas_r, acv_r, rho_r, un_r, T_r, hf_r, & flurho, fluru, fluret) C In order to obtain the flux in the true left-right frame, I have C to multiply the speed * rind_l flurho = flurho * rind_l fluret = fluret * rind_l flu(1) = flurho flu(2) = fluru flu(2+ndim) = fluret flu(3+ndim) = 0.0d0 C Note that the Larrouturou theorem does not hold C Namely un_ls (= un_rss) has not the same sign as flurho C write (*,*) un_ls, un_rss, flurho C un_ls > 0 => coefy = 1 coefy = 0.5d0 * (1.0d0 + sign(1.0d0, (rind_l * un_ls))) 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 Kinetic energy flux trekin = 0.0D0 do icomp = 2, ndim, 1 trekin = trekin + (((coefy * wvec_l(icomp + 1)) + & (omcoey * wvec_r(icomp + 1))) ** 2) enddo trekin = 0.5D0 * trekin flu(2+ndim) = flu(2+ndim) + (flurho * trekin) C C**** Lagrangian flux over the reactive shock C xst = D_resh call flufun(nordpo, & Tmaxcv, acv_w , & xst, & vit_l, vit_ls, un_ls, D_resh, vit_rs, vit_r, & Rgas_l, acv_l, rho_l, un_l, T_l, hf_l, & rho_ls, T_ls, & Rg_rss, acv_z, rh_rss, un_rss, T_rss, hf_rss, & rho_rs, un_rs, T_rs, & Rgas_r, acv_r, rho_r, un_r, T_r, hf_r, & flurho, fluru, fluret) C flurho = flurho * rind_l fluret = fluret * rind_l f_lag(1) = flurho f_lag(2) = fluru f_lag(2+ndim) = fluret f_lag(3+ndim) = -1.0D0 * rind_l * D_resh C C Since we suppose that left state is the burnt one, C D_resh > un_ls, i.e. over (x/t) = D_rest we have the C right passive scalar ! C indl = 1 => burnt state in wvec_l => passive scalar in r C indl = -1 => burnt state in wvec_r => passive scalar in l C Note that the result does not change as q-> 0 C (virtual combustion wave). C coefy = 0.5d0 * (1.0d0 + sign(1.0d0, (1.0D0 * rind_l))) omcoey = 1.0d0 - coefy C do icomp = 3, (ndim + 1), 1 f_lag(icomp) = flurho * ((coefy * wvec_r(icomp)) + & (omcoey * wvec_l(icomp))) enddo do icomp = (ndim + 4), (ndim + nesp + 4), 1 f_lag(icomp) = flurho * ((coefy * wvec_r(icomp)) + & (omcoey * wvec_l(icomp))) enddo C Kinetic energy flux trekin = 0.0D0 do icomp = 2, ndim, 1 trekin = trekin + (((coefy * wvec_r(icomp + 1)) + & (omcoey * wvec_l(icomp + 1))) ** 2) enddo trekin = 0.5D0 * trekin C write(*,*) trekin f_lag(2+ndim) = f_lag(2+ndim) + (flurho * trekin) C 9999 continue return end