3.4
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
cctpfa/maxwellstefanslaw.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 *****************************************************************************/
25#ifndef DUMUX_DISCRETIZATION_CC_TPFA_MAXWELL_STEFAN_LAW_HH
26#define DUMUX_DISCRETIZATION_CC_TPFA_MAXWELL_STEFAN_LAW_HH
27
28#include <dune/common/fmatrix.hh>
29#include <dune/common/float_cmp.hh>
30
33
39
40namespace Dumux {
41
42// forward declaration
43template <class TypeTag, DiscretizationMethod discMethod, ReferenceSystemFormulation referenceSystem>
44class MaxwellStefansLawImplementation;
45
50template <class TypeTag, ReferenceSystemFormulation referenceSystem>
52{
56 using FVElementGeometry = typename GridGeometry::LocalView;
57 using SubControlVolume = typename GridGeometry::SubControlVolume;
58 using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
59 using Extrusion = Extrusion_t<GridGeometry>;
60 using GridView = typename GridGeometry::GridView;
62 using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
63 using Element = typename GridView::template Codim<0>::Entity;
64 using ElementFluxVariablesCache = typename GetPropType<TypeTag, Properties::GridFluxVariablesCache>::LocalView;
67
68 static const int dim = GridView::dimension;
69 static const int dimWorld = GridView::dimensionworld;
72
73 static_assert(referenceSystem == ReferenceSystemFormulation::massAveraged, "only the mass averaged reference system is supported for the Maxwell-Stefan formulation");
74
75
76 using ComponentFluxVector = Dune::FieldVector<Scalar, numComponents>;
77 using ReducedComponentVector = Dune::FieldVector<Scalar, numComponents-1>;
78 using ReducedComponentMatrix = Dune::FieldMatrix<Scalar, numComponents-1, numComponents-1>;
79
80public:
81 // state the discretization method this implementation belongs to
83 //return the reference system
85 { return referenceSystem; }
86
88
93
99 static ComponentFluxVector flux(const Problem& problem,
100 const Element& element,
101 const FVElementGeometry& fvGeometry,
102 const ElementVolumeVariables& elemVolVars,
103 const SubControlVolumeFace& scvf,
104 const int phaseIdx,
105 const ElementFluxVariablesCache& elemFluxVarsCache)
106 {
107 //this is to calculate the maxwellStefan diffusion in a multicomponent system.
108 //see: Multicomponent Mass Transfer. R. Taylor u. R. Krishna. J. Wiley & Sons, New York 1993
109 ComponentFluxVector componentFlux(0.0);
110 ReducedComponentVector moleFracInside(0.0);
111 ReducedComponentVector moleFracOutside(0.0);
112 ReducedComponentVector reducedFlux(0.0);
113 ReducedComponentMatrix reducedDiffusionMatrixInside(0.0);
114 ReducedComponentMatrix reducedDiffusionMatrixOutside(0.0);
115
116 // get inside/outside volume variables
117 const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
118 const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
119 const auto rhoInside = insideVolVars.density(phaseIdx);
120 const auto rhoOutside = outsideVolVars.density(phaseIdx);
121 //calculate the mole fraction vectors
122 for (int compIdx = 0; compIdx < numComponents-1; compIdx++)
123 {
124 //calculate x_inside
125 const auto xInside = insideVolVars.moleFraction(phaseIdx, compIdx);
126 //calculate outside molefraction with the respective transmissibility
127 const auto xOutside = outsideVolVars.moleFraction(phaseIdx, compIdx);
128
129 moleFracInside[compIdx] = xInside;
130 moleFracOutside[compIdx] = xOutside;
131 }
132
133 //we cannot solve that if the matrix is 0 everywhere
134 if(!(Dune::FloatCmp::eq<Scalar>(insideVolVars.saturation(phaseIdx), 0) || Dune::FloatCmp::eq<Scalar>(outsideVolVars.saturation(phaseIdx), 0)))
135 {
136 const auto insideScvIdx = scvf.insideScvIdx();
137 const auto& insideScv = fvGeometry.scv(insideScvIdx);
138 const Scalar omegai = calculateOmega_(scvf,
139 insideScv,
140 insideVolVars.extrusionFactor());
141
142 //now we have to do the tpfa: J_i = J_j which leads to: tij(xi -xj) = -rho Bi^-1 omegai(x*-xi) with x* = (omegai Bi^-1 + omegaj Bj^-1)^-1 (xi omegai Bi^-1 + xj omegaj Bj^-1) with i inside and j outside
143 reducedDiffusionMatrixInside = setupMSMatrix_(problem, element, fvGeometry, insideVolVars, insideScv, phaseIdx);
144
145 //if on boundary
146 if (scvf.boundary() || scvf.numOutsideScvs() > 1)
147 {
148 moleFracOutside -= moleFracInside;
149 reducedDiffusionMatrixInside.solve(reducedFlux, moleFracOutside);
150 reducedFlux *= omegai*rhoInside;
151 }
152 else //we need outside cells as well if we are not on the boundary
153 {
154 Scalar omegaj;
155 const auto outsideScvIdx = scvf.outsideScvIdx();
156 const auto& outsideScv = fvGeometry.scv(outsideScvIdx);
157
158 reducedDiffusionMatrixOutside = setupMSMatrix_(problem, element, fvGeometry, outsideVolVars, outsideScv, phaseIdx);
159
160 if (dim == dimWorld)
161 // assume the normal vector from outside is anti parallel so we save flipping a vector
162 omegaj = -1.0*calculateOmega_(scvf,
163 outsideScv,
164 outsideVolVars.extrusionFactor());
165 else
166 omegaj = calculateOmega_(fvGeometry.flipScvf(scvf.index()),
167 outsideScv,
168 outsideVolVars.extrusionFactor());
169
170 reducedDiffusionMatrixInside.invert();
171 reducedDiffusionMatrixOutside.invert();
172 reducedDiffusionMatrixInside *= omegai*rhoInside;
173 reducedDiffusionMatrixOutside *= omegaj*rhoOutside;
174
175 //in the helpervector we store the values for x*
176 ReducedComponentVector helperVector(0.0);
177 ReducedComponentVector gradientVectori(0.0);
178 ReducedComponentVector gradientVectorj(0.0);
179
180 reducedDiffusionMatrixInside.mv(moleFracInside, gradientVectori);
181 reducedDiffusionMatrixOutside.mv(moleFracOutside, gradientVectorj);
182
183 auto gradientVectorij = (gradientVectori + gradientVectorj);
184
185 //add the two matrixes to each other
186 reducedDiffusionMatrixOutside += reducedDiffusionMatrixInside;
187
188 reducedDiffusionMatrixOutside.solve(helperVector, gradientVectorij);
189
190 //Bi^-1 omegai rho(x*-xi)
191 helperVector -=moleFracInside;
192 reducedDiffusionMatrixInside.mv(helperVector, reducedFlux);
193 }
194
195 reducedFlux *= -Extrusion::area(scvf);
196 for (int compIdx = 0; compIdx < numComponents-1; compIdx++)
197 {
198 componentFlux[compIdx] = reducedFlux[compIdx];
199 componentFlux[numComponents-1] -=reducedFlux[compIdx];
200 }
201 }
202 return componentFlux ;
203 }
204
205private:
206 static Scalar calculateOmega_(const SubControlVolumeFace& scvf,
207 const SubControlVolume &scv,
208 const Scalar extrusionFactor)
209 {
210 auto distanceVector = scvf.ipGlobal();
211 distanceVector -= scv.center();
212 distanceVector /= distanceVector.two_norm2();
213
214 Scalar omega = (distanceVector * scvf.unitOuterNormal());
215 omega *= extrusionFactor;
216
217 return omega;
218 }
219
220 static ReducedComponentMatrix setupMSMatrix_(const Problem& problem,
221 const Element& element,
222 const FVElementGeometry& fvGeometry,
223 const VolumeVariables& volVars,
224 const SubControlVolume& scv,
225 const int phaseIdx)
226 {
227 ReducedComponentMatrix reducedDiffusionMatrix(0.0);
228
229 //this is to not devide by 0 if the saturation in 0 and the effectiveDiffusionCoefficient becomes zero due to that
230 if(Dune::FloatCmp::eq<Scalar>(volVars.saturation(phaseIdx), 0))
231 return reducedDiffusionMatrix;
232
233 const auto avgMolarMass = volVars.averageMolarMass(phaseIdx);
234
235 for (int compIIdx = 0; compIIdx < numComponents-1; compIIdx++)
236 {
237 const auto xi = volVars.moleFraction(phaseIdx, compIIdx);
238
239 const auto Mn = FluidSystem::molarMass(numComponents-1);
240 Scalar tin = volVars.effectiveDiffusionCoefficient(phaseIdx, compIIdx, numComponents-1);
241
242 // set the entries of the diffusion matrix of the diagonal
243 reducedDiffusionMatrix[compIIdx][compIIdx] += xi*avgMolarMass/(tin*Mn);
244
245 // now set the rest of the entries (off-diagonal and additional entries for diagonal)
246 for (int compJIdx = 0; compJIdx < numComponents; compJIdx++)
247 {
248 // we don't want to calculate e.g. water in water diffusion
249 if (compJIdx == compIIdx)
250 continue;
251
252 const auto xj = volVars.moleFraction(phaseIdx, compJIdx);
253 const auto Mi = FluidSystem::molarMass(compIIdx);
254 const auto Mj = FluidSystem::molarMass(compJIdx);
255 Scalar tij = volVars.effectiveDiffusionCoefficient(phaseIdx, compIIdx, compJIdx);
256 reducedDiffusionMatrix[compIIdx][compIIdx] += xj*avgMolarMass/(tij*Mi);
257
258 if (compJIdx < numComponents-1)
259 reducedDiffusionMatrix[compIIdx][compJIdx] += xi*(avgMolarMass/(tin*Mn) - avgMolarMass/(tij*Mj));
260 }
261 }
262 return reducedDiffusionMatrix;
263 }
264
265};
266} // end namespace
267
268#endif
The infrastructure to retrieve run-time parameters from Dune::ParameterTrees.
The available discretization methods in Dumux.
Helper classes to compute the integration elements.
The reference frameworks and formulations available for splitting total fluxes into a advective and d...
Container storing the diffusion coefficients required by the Maxwell- Stefan diffusion law....
Classes related to flux variables caching.
DiscretizationMethod
The available discretization methods in Dumux.
Definition: method.hh:37
ReferenceSystemFormulation
The formulations available for Fick's law related to the reference system.
Definition: referencesystemformulation.hh:45
Definition: adapt.hh:29
typename Extrusion< T >::type Extrusion_t
Convenience alias for obtaining the extrusion type.
Definition: extrusion.hh:177
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property
Definition: propertysystem.hh:150
Definition: maxwellstefanslaw.hh:37
static constexpr ReferenceSystemFormulation referenceSystemFormulation()
Definition: cctpfa/maxwellstefanslaw.hh:84
static ComponentFluxVector flux(const Problem &problem, const Element &element, const FVElementGeometry &fvGeometry, const ElementVolumeVariables &elemVolVars, const SubControlVolumeFace &scvf, const int phaseIdx, const ElementFluxVariablesCache &elemFluxVarsCache)
Returns the diffusive fluxes of all components within a fluid phase across the given sub-control volu...
Definition: cctpfa/maxwellstefanslaw.hh:99
The empty filler class corresponding to EmptyCache.
Definition: fluxvariablescaching.hh:32
Definition: fluxvariablescaching.hh:67
Container storing the diffusion coefficients required by the Maxwell- Stefan diffusion law....
Definition: maxwellstefandiffusioncoefficients.hh:45
Declares all properties used in Dumux.