3.5-git
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
misciblemultiphasecomposition.hh
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*****************************************************************************
4 * See the file COPYING for full copying permissions. *
5 * *
6 * This program is free software: you can redistribute it and/or modify *
7 * it under the terms of the GNU General Public License as published by *
8 * the Free Software Foundation, either version 3 of the License, or *
9 * (at your option) any later version. *
10 * *
11 * This program is distributed in the hope that it will be useful, *
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
14 * GNU General Public License for more details. *
15 * *
16 * You should have received a copy of the GNU General Public License *
17 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
18 *****************************************************************************/
26#ifndef DUMUX_MISCIBLE_MULTIPHASE_COMPOSITION_HH
27#define DUMUX_MISCIBLE_MULTIPHASE_COMPOSITION_HH
28
29#include <dune/common/fvector.hh>
30#include <dune/common/fmatrix.hh>
31
33
34namespace Dumux {
57template <class Scalar, class FluidSystem>
59{
60 static constexpr int numPhases = FluidSystem::numPhases;
61 static constexpr int numComponents = FluidSystem::numComponents;
62 static const int numMajorComponents = FluidSystem::numPhases;
63
64public:
80 template <class FluidState, class ParameterCache>
81 static void solve(FluidState &fluidState,
82 ParameterCache &paramCache,
83 int knownPhaseIdx = 0)
84 {
85#ifndef NDEBUG
86 // currently this solver can only handle fluid systems which
87 // assume ideal mixtures of all fluids. TODO: relax this
88 // (requires solving a non-linear system of equations, i.e. using
89 // Newton method.)
90 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
91 assert(FluidSystem::isIdealMixture(phaseIdx));
92
93 }
94#endif
95
96 //get the known mole fractions from the fluidState
97 //in a 2pnc system the n>2 mole fractions are primary variables and are already set in the fluidstate
98 Dune::FieldVector<Scalar, numComponents-numMajorComponents> xKnown(0.0);
99 for (int knownCompIdx = 0; knownCompIdx < numComponents-numMajorComponents; ++knownCompIdx)
100 {
101 xKnown[knownCompIdx] = fluidState.moleFraction(knownPhaseIdx, knownCompIdx + numMajorComponents);
102 }
103
104 // compute all fugacity coefficients
105 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
106 paramCache.updatePhase(fluidState, phaseIdx);
107
108 // since we assume ideal mixtures, the fugacity
109 // coefficients of the components cannot depend on
110 // composition, i.e. the parameters in the cache are valid
111 for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
112 Scalar fugCoeff = FluidSystem::fugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx);
113 fluidState.setFugacityCoefficient(phaseIdx, compIdx, fugCoeff);
114 }
115 }
116
117
118 // create the linear system of equations which defines the
119 // mole fractions
120 Dune::FieldMatrix<Scalar, numComponents*numPhases, numComponents*numPhases> M(0.0);
121 Dune::FieldVector<Scalar, numComponents*numPhases> x(0.0);
122 Dune::FieldVector<Scalar, numComponents*numPhases> b(0.0);
123
124 // assemble the equations expressing the assumption that the
125 // sum of all mole fractions in each phase must be 1
126 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
127 int rowIdx = numComponents*(numPhases - 1) + phaseIdx;
128 b[rowIdx] = 1.0;
129
130 for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
131 int colIdx = phaseIdx*numComponents + compIdx;
132
133 M[rowIdx][colIdx] = 1.0;
134 }
135 }
136
137 // set the additional equations for the numComponents-numMajorComponents
138 // this is only relevant if numphases != numcomponents e.g. in a 2pnc system
139 // Components, of which the mole fractions are known, set to molefraction(knownCompIdx)=xKnown
140 for(int knownCompIdx = 0; knownCompIdx < numComponents-numMajorComponents; ++knownCompIdx)
141 {
142 int rowIdx = numComponents + numPhases + knownCompIdx;
143 int colIdx = knownPhaseIdx*numComponents + knownCompIdx + numMajorComponents;
144 M[rowIdx][colIdx] = 1.0;
145 b[rowIdx] = xKnown[knownCompIdx];
146 }
147
148 // assemble the equations expressing the fact that the
149 // fugacities of each component are equal in all phases
150 for (int compIdx = 0; compIdx < numComponents; ++compIdx)
151 {
152 int col1Idx = compIdx;
153 const auto entryPhase0 = fluidState.fugacityCoefficient(0, compIdx)*fluidState.pressure(0);
154
155 for (unsigned int phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx)
156 {
157 int rowIdx = (phaseIdx - 1)*numComponents + compIdx;
158 int col2Idx = phaseIdx*numComponents + compIdx;
159 M[rowIdx][col1Idx] = entryPhase0;
160 M[rowIdx][col2Idx] = -fluidState.fugacityCoefficient(phaseIdx, compIdx)*fluidState.pressure(phaseIdx);
161 }
162 }
163
164 // preconditioning of M to reduce condition number
165 for (int compIdx = 0; compIdx < numComponents; compIdx++)
166 {
167 // Multiply row of main component (Raoult's Law) with 10e-5 (order of magn. of pressure)
168 if (compIdx < numMajorComponents)
169 M[compIdx] *= 10e-5;
170
171 // Multiply row of sec. components (Henry's Law) with 10e-9 (order of magn. of Henry constant)
172 else
173 M[compIdx] *= 10e-9;
174 }
175
176 // solve for all mole fractions
177 try { M.solve(x, b); }
178 catch (const Dune::FMatrixError& e) {
179 DUNE_THROW(NumericalProblem,
180 "MiscibleMultiPhaseComposition: Failed to solve the linear equation system for the phase composition.");
181 }
182
183 // set all mole fractions and the the additional quantities in
184 // the fluid state
185 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
186 for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
187 int rowIdx = phaseIdx*numComponents + compIdx;
188 fluidState.setMoleFraction(phaseIdx, compIdx, x[rowIdx]);
189 }
190 paramCache.updateComposition(fluidState, phaseIdx);
191
192 Scalar value = FluidSystem::density(fluidState, paramCache, phaseIdx);
193 fluidState.setDensity(phaseIdx, value);
194
195 value = FluidSystem::molarDensity(fluidState, paramCache, phaseIdx);
196 fluidState.setMolarDensity(phaseIdx, value);
197 }
198 }
199};
200
201} // end namespace Dumux
202
203#endif
Some exceptions thrown in DuMux
Definition: adapt.hh:29
std::string molarDensity(int phaseIdx) noexcept
I/O name of molar density for multiphase systems.
Definition: name.hh:83
std::string density(int phaseIdx) noexcept
I/O name of density for multiphase systems.
Definition: name.hh:65
Exception thrown if a fixable numerical problem occurs.
Definition: exceptions.hh:39
Computes the composition of all phases of a N-phase, N-component fluid system assuming that all N pha...
Definition: misciblemultiphasecomposition.hh:59
static void solve(FluidState &fluidState, ParameterCache &paramCache, int knownPhaseIdx=0)
Computes the composition of all phases of a N-phase, N-component fluid system assuming that all N pha...
Definition: misciblemultiphasecomposition.hh:81