35 enum { numComponents = FluidSystem::numComponents };
40 template <
class Flu
idState,
class Params,
class LhsEval =
typename Flu
idState::Scalar>
41 static LhsEval computeFugacityCoefficient(
const FluidState& fs,
51 LhsEval Vm = params.molarVolume(phaseIdx);
52 LhsEval T = fs.temperature(phaseIdx);
53 LhsEval p = fs.pressure(phaseIdx);
54 LhsEval A = params.A(phaseIdx);
55 LhsEval B = params.B(phaseIdx);
56 LhsEval Bi = params.Bi(phaseIdx, compIdx);
57 LhsEval m1 = params.m1(phaseIdx);
58 LhsEval m2 = params.m2(phaseIdx);
61 LhsEval Bi_B = Bi / B;
65 LhsEval Z = p * Vm / RT;
69 for (
unsigned compJIdx = 0; compJIdx < numComponents; ++compJIdx) {
70 A_s += params.aCache(phaseIdx, compIdx, compJIdx) * fs.moleFraction(phaseIdx, compJIdx);
80 alpha = -log(Z - B) + Bi_B * (Z - 1);
81 beta = log((Z + m2 * B) / (Z + m1 * B)) * A / ((m1 - m2) * B);
82 gamma = (2 / A) * A_s - Bi_B;
83 ln_phi = alpha + (beta * gamma);
85 fugCoeff = exp(ln_phi);
93 fugCoeff = min(1e10, fugCoeff);
98 fugCoeff = max(1e-10, fugCoeff);
104 template <
class Flu
idState,
class Params>
105 static typename FluidState::Scalar computeMolarVolume(
const FluidState& fs,
110 Valgrind::CheckDefined(fs.temperature(phaseIdx));
111 Valgrind::CheckDefined(fs.pressure(phaseIdx));
113 using Evaluation =
typename FluidState::Scalar;
116 const Evaluation& T = fs.temperature(phaseIdx);
117 const Evaluation& p = fs.pressure(phaseIdx);
118 const Evaluation& A = params.A(phaseIdx);
119 const Evaluation& B = params.B(phaseIdx);
120 const Evaluation& m1 = params.m1(phaseIdx);
121 const Evaluation& m2 = params.m2(phaseIdx);
124 const Evaluation a1 = 1.0;
125 const Evaluation a2 = (m1 + m2 - 1) * B - 1;
126 const Evaluation a3 = A + m1 * m2 * B * B - (m1 + m2) * B * (B + 1);
127 const Evaluation a4 = -A * B - m1 * m2 * B * B * (B + 1);
128 Valgrind::CheckDefined(a1);
129 Valgrind::CheckDefined(a2);
130 Valgrind::CheckDefined(a3);
131 Valgrind::CheckDefined(a4);
135 Valgrind::SetUndefined(Vm);
136 Evaluation Z[3] = {0.0, 0.0, 0.0};
140 const Evaluation RT_p = R * T / p;
146 Vm = max(1e-7, Z[2] * RT_p);
148 Vm = max(1e-7, Z[0] * RT_p);
150 else if (numSol == 1) {
154 Vm = max(1e-7, Z[0] * RT_p);
157 Valgrind::CheckDefined(Vm);
158 assert(std::isfinite(scalarValue(Vm)));
unsigned cubicRoots(SolContainer *sol, Scalar a, Scalar b, Scalar c, Scalar d)
Invert a cubic polynomial analytically.
Definition PolynomialUtils.hpp:331