159 static void solve(FluidState& fluidState,
160 ParameterCache& paramCache,
163 unsigned numAuxConstraints,
165 bool setInternalEnergy)
167 static_assert(std::is_same<typename FluidState::Scalar, Evaluation>::value,
168 "The scalar type of the fluid state must be 'Evaluation'");
175 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
176 assert(FluidSystem::isIdealMixture(phaseIdx));
181 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
182 paramCache.updatePhase(fluidState, phaseIdx);
187 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
188 Evaluation fugCoeff = decay<Evaluation>(
189 FluidSystem::fugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx));
190 fluidState.setFugacityCoefficient(phaseIdx, compIdx, fugCoeff);
196 static const int numEq = numComponents*numPhases;
197 Dune::FieldMatrix<Evaluation, numEq, numEq> M(0.0);
198 Dune::FieldVector<Evaluation, numEq> x(0.0);
199 Dune::FieldVector<Evaluation, numEq> b(0.0);
203 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
204 const Evaluation& entryCol1 =
205 fluidState.fugacityCoefficient(0, compIdx)
206 *fluidState.pressure(0);
207 unsigned col1Idx = compIdx;
209 for (
unsigned phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx) {
210 unsigned rowIdx = (phaseIdx - 1)*numComponents + compIdx;
211 unsigned col2Idx = phaseIdx*numComponents + compIdx;
213 const Evaluation& entryCol2 =
214 fluidState.fugacityCoefficient(phaseIdx, compIdx)
215 *fluidState.pressure(phaseIdx);
217 M[rowIdx][col1Idx] = entryCol1;
218 M[rowIdx][col2Idx] = -entryCol2;
225 unsigned presentPhases = 0;
226 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
227 if (!(phasePresence& (1 << phaseIdx)))
230 unsigned rowIdx = numComponents*(numPhases - 1) + presentPhases;
234 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
235 unsigned colIdx = phaseIdx*numComponents + compIdx;
237 M[rowIdx][colIdx] = 1.0;
241 assert(presentPhases + numAuxConstraints == numComponents);
244 for (
unsigned auxEqIdx = 0; auxEqIdx < numAuxConstraints; ++auxEqIdx) {
245 unsigned rowIdx = numComponents*(numPhases - 1) + presentPhases + auxEqIdx;
246 b[rowIdx] = auxConstraints[auxEqIdx].
value();
248 unsigned colIdx = auxConstraints[auxEqIdx].
phaseIdx()*numComponents + auxConstraints[auxEqIdx].
compIdx();
249 M[rowIdx][colIdx] = 1.0;
256 catch (
const Dune::FMatrixError& e) {
257 std::ostringstream oss;
258 oss <<
"Numerical problem in MiscibleMultiPhaseComposition::solve(): " << e.what() <<
"; M="<<M;
268 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
269 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
270 unsigned rowIdx = phaseIdx*numComponents + compIdx;
271 fluidState.setMoleFraction(phaseIdx, compIdx, x[rowIdx]);
273 paramCache.updateComposition(fluidState, phaseIdx);
275 const Evaluation& rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
276 fluidState.setDensity(phaseIdx, rho);
279 const Evaluation& mu = FluidSystem::viscosity(fluidState, paramCache, phaseIdx);
280 fluidState.setViscosity(phaseIdx, mu);
283 if (setInternalEnergy) {
284 const Evaluation& h = FluidSystem::enthalpy(fluidState, paramCache, phaseIdx);
285 fluidState.setEnthalpy(phaseIdx, h);