Equation of State
PelePhysics allows the user to use different equation of state (EOS) as the constitutive equation and close the compressible Navier-Stokes system of equations. All the routines needed to fully define an EOS are implemented through PelePhysics module. Available models include:
An ideal gas mixture model (similar to the CHEMKIN-II approach)
A simple GammaLaw model
Cubic models such as Soave-Redlich-Kwong; Peng-Robinson support was started but is currently stalled.
Examples of EOS implementation can be seen in PelePhysics/Eos. The following sections will fully describe the implementation of Soave-Redlich-Kwong, a non-ideal cubic EOS, for a general mixture of species. Integration with CEPTR, for a chemical mechanism described in a chemkin format, will also be highlighted. For an advanced user interested in implementing a new EOS this chapter should provide a good starting point.
Soave-Redlich-Kwong (SRK)
The cubic model is built on top of the ideal gas models, and is selected by specifying its name as the Eos_dir during the build (the Chemistry_Model must also be specified). Any additional parameters (e.g., attractions, repulsions, critical states) are either included in the underlying FUEGO database used to generate the source file model implementation, or else are inferred from the input model data.
SRK EOS as a function of Pressure (p), Temperature(T), and \(\tau\) (specific volume) is given by
where \(Y_k\) are species mass fractions, \(R\) is the universal gas constant, and \(b_m\) and \(a_m\) are mixture repulsion and attraction terms, respectively.
Mixing rules
For a mixture of species, the following mixing rules are used to compute \(b_m\) and \(a_m\).
where \(b_i\) and \(a_i\) for each species is defined using critical pressure and temperature.
where
where \(\omega_i\) are the accentric factors and
For chemically unstable species such as radicals, critical temperatures and pressures are not available. For species where critical properties are not available, we use the Lennard-Jones potential for that species to construct attractive and repulsive coefficients.
where \(\sigma_i\), \(\epsilon_i\) are the Lennard-Jones potential molecular diameter and well-depth, respectively, \(m_i\) the molecular mass, and \(k_b\) is Boltzmann’s constant.
In terms of implementation, a routine called MixingRuleAmBm can be found in the SRK eos implementation. The following code block shows the subroutine which receives species mass fractions and temperature as input. The outputs of this routine are \(b_m\) and \(a_m\) .
do i = 1, nspecies
Tr = T*oneOverTc(i)
amloc(i) = (1.0d0 + Fomega(i)*(1.0d0-sqrt(Tr))) *sqrtAsti(i)
bm = bm + massFrac(i)*Bi(i)
enddo
do j = 1, nspecies
do i = 1, nspecies
am = am + massFrac(i)*massFrac(j)*amloc(i)*amloc(j)
end do
end do
Thermodynamic Properties
Most of the thermodynamic properties can be calculated from the equation of state and involve derivatives of various thermodynamic quantities and of EOS parameters. In the following, some of these thermodynamic properties for SRK and the corresponding routines are presented.
Specific heat
For computing mixture specific heat at constant volume and pressure, the ideal gas contribution and the departure from the ideal gas are computed. Specific heat at constant volume can be computed using the following
For SRK EOS, the formula for \(c_v\) reduces to
where \(c_v^{id}\) is the specific heat at constant volume. Mixture specific heat at constant volume is implemented through the routine SRK_EOS_GetMixtureCv
subroutine SRK_EOS_GetMixtureCv(state)
implicit none
type (eos_t), intent(inout) :: state
real(amrex_real) :: tau, K1
state % wbar = 1.d0 / sum(state % massfrac(:) * inv_mwt(:))
call MixingRuleAmBm(state%T,state%massFrac,state%am,state%bm)
tau = 1.0d0/state%rho
! Derivative of the EOS AM w.r.t Temperature - needed for calculating enthalpy, Cp, Cv and internal energy
call Calc_dAmdT(state%T,state%massFrac,state%am,state%dAmdT)
! Second Derivative of the EOS AM w.r.t Temperature - needed for calculating enthalpy, Cp, Cv and internal energy
call Calc_d2AmdT2(state%T,state%massFrac,state%d2AmdT2)
! Ideal gas specific heat at constant volume
call ckcvbs(state%T, state % massfrac, iwrk, rwrk, state % cv)
! Real gas specific heat at constant volume
state%cv = state%cv + state%T*state%d2AmdT2* (1.0d0/state%bm)*log(1.0d0+state%bm/tau)
end subroutine SRK_EOS_GetMixtureCv
Specific heat at constant pressure is given by
where all the derivatives in the above expression for SRK EOS are given by
subroutine SRK_EOS_GetMixtureCp(state)
implicit none
type (eos_t), intent(inout) :: state
real(amrex_real) :: tau, K1
real(amrex_real) :var: : Cpig
real(amrex_real) :: eosT1Denom, eosT2Denom, eosT3Denom
real(amrex_real) :: InvEosT1Denom,InvEosT2Denom,InvEosT3Denom
real(amrex_real) :: dhmdT,dhmdtau
real(amrex_real) :: Rm
state % wbar = 1.d0 / sum(state % massfrac(:) * inv_mwt(:))
call MixingRuleAmBm(state%T,state%massFrac,state%am,state%bm)
tau = 1.0d0/state%rho
! Derivative of the EOS AM w.r.t Temperature - needed for calculating enthalpy, Cp, Cv and internal energy
call Calc_dAmdT(state%T,state%massFrac,state%dAmdT)
! Second Derivative of the EOS AM w.r.t Temperature - needed for calculating enthalpy, Cp, Cv and internal energy
call Calc_d2AmdT2(state%T,state%massFrac,state%d2AmdT2)
K1 = (1.0d0/state%bm)*log(1.0d0+state%bm/tau)
eosT1Denom = tau-state%bm
eosT2Denom = tau*(tau+state%bm)
eosT3Denom = tau+state%bm
InvEosT1Denom = 1.0d0/eosT1Denom
InvEosT2Denom = 1.0d0/eosT2Denom
InvEosT3Denom = 1.0d0/eosT3Denom
Rm = (Ru/state%wbar)
! Derivative of Pressure w.r.t to Temperature
state%dPdT = Rm*InvEosT1Denom - state%dAmdT*InvEosT2Denom
! Derivative of Pressure w.r.t to tau (specific volume)
state%dpdtau = -Rm*state%T*InvEosT1Denom*InvEosT1Denom + state%am*(2.0*tau+state%bm)*InvEosT2Denom*InvEosT2Denom
! Ideal gas specific heat at constant pressure
call ckcpbs(state % T, state % massfrac, iwrk, rwrk,Cpig)
! Derivative of enthalpy w.r.t to Temperature
dhmdT = Cpig + state%T*state%d2AmdT2*K1 - state%dAmdT*InvEosT3Denom + Rm*state%bm*InvEosT1Denom
! Derivative of enthalpy w.r.t to tau (specific volume)
dhmdtau = -(state%T*state%dAmdT - state%am)*InvEosT2Denom + state%am*InvEosT3Denom*InvEosT3Denom - &
Rm*state%T*state%bm*InvEosT1Denom*InvEosT1Denom
! Real gas specific heat at constant pressure
state%cp = dhmdT - (dhmdtau/state%dpdtau)*state%dPdT
end subroutine SRK_EOS_GetMixtureCp
Internal energy and Enthalpy
Similarly mixture internal energy for SRK EOS is given by
and mixture enthalpy \(h_m = e + p \tau\)
and the implementation can be found in the routine SRK_EOS_GetMixture_H.
Speed of Sound
The sound speed for SRK EOS is given by
Species enthalpy
For computation of kinetics and transport fluxes we will also need the species partial enthalpies and the chemical potential. The species enthalpies for SRK EOS are given by
where
Chemical potential
The chemical potentials are the derivative of the free energy with respect to composition. Here the free energy f` is given by
Then
Other primitive variable derivatives
The Godunov (FV) algorithm also needs some derivatives to express source terms in terms of primitive variables. In particular one needs
and
All of the terms needed to evaluate this quantity are known except for