46 template <
class Flu
idState,
class Params,
class LhsEval =
typename Flu
idState::Scalar>
47 static LhsEval LBC(
const FluidState& fluidState,
51 const Scalar MPa_atm = 0.101325;
53 const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
54 const auto& P = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
55 const auto& Z = Opm::decay<LhsEval>(fluidState.compressFactor(phaseIdx));
57 LhsEval sumVolume = 0.0;
58 for (
unsigned compIdx = 0; compIdx < FluidSystem::numComponents; ++compIdx) {
59 const auto& x = Opm::decay<LhsEval>(fluidState.moleFraction(phaseIdx, compIdx));
60 const Scalar v_c = FluidSystem::criticalVolume(compIdx) / 1000;
64 LhsEval rho_pc = 1.0 / sumVolume;
65 LhsEval V = (R * T * Z)/P;
66 LhsEval rho = 1.0 / V;
67 LhsEval rho_r = rho / rho_pc;
69 LhsEval xsum_T_c = 0.0;
70 LhsEval xsum_Mm = 0.0;
71 LhsEval xsum_p_ca = 0.0;
72 for (
unsigned compIdx = 0; compIdx < FluidSystem::numComponents; ++compIdx) {
73 const Scalar& p_c = FluidSystem::criticalPressure(compIdx) / 1e6;
74 const Scalar& T_c = FluidSystem::criticalTemperature(compIdx);
75 const Scalar Mm = FluidSystem::molarMass(compIdx) * 1000;
76 const auto& x = Opm::decay<LhsEval>(fluidState.moleFraction(phaseIdx, compIdx));
77 Scalar p_ca = p_c / MPa_atm;
80 xsum_p_ca += x * p_ca;
82 LhsEval zeta_tot = Opm::pow(xsum_T_c / (Opm::pow(xsum_Mm,3.0) * Opm::pow(xsum_p_ca,4.0)),1./6);
86 for (
unsigned compIdx = 0; compIdx < FluidSystem::numComponents; ++compIdx) {
87 const Scalar& p_c = FluidSystem::criticalPressure(compIdx) / 1e6;
88 const Scalar& T_c = FluidSystem::criticalTemperature(compIdx);
89 const Scalar Mm = FluidSystem::molarMass(compIdx) * 1000;
90 const auto& x = Opm::decay<LhsEval>(fluidState.moleFraction(phaseIdx, compIdx));
91 Scalar p_ca = p_c / MPa_atm;
92 Scalar zeta = std::pow(T_c / (std::pow(Mm,3.0) * std::pow(p_ca,4.0)),1./6);
94 LhsEval xrM = x * std::pow(Mm,0.5);
97 mys = 34.0e-5*Opm::pow(T_r,0.94)/zeta;
99 mys = 17.78e-5*Opm::pow(4.58*T_r - 1.67, 0.625)/zeta;
106 std::vector<Scalar> LBC = {0.10230,
112 LhsEval sumLBC = 0.0;
113 for (
int i = 0; i < 5; ++i) {
114 sumLBC += Opm::pow(rho_r,i)*LBC[i];
117 return (my0 + (Opm::pow(sumLBC,4.0) - 1e-4)/zeta_tot)/1e3;
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition Exceptions.hpp:30