3.3.0
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
staggered/freeflow/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 *****************************************************************************/
24#ifndef DUMUX_DISCRETIZATION_STAGGERED_MAXWELL_STEFAN_LAW_HH
25#define DUMUX_DISCRETIZATION_STAGGERED_MAXWELL_STEFAN_LAW_HH
26
27#include <dune/common/fmatrix.hh>
28
33
37
38namespace Dumux {
39
40// forward declaration
41template <class TypeTag, DiscretizationMethod discMethod, ReferenceSystemFormulation referenceSystem>
42class MaxwellStefansLawImplementation;
43
48template <class TypeTag, ReferenceSystemFormulation referenceSystem>
50{
54 using FVElementGeometry = typename GridGeometry::LocalView;
55 using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
56 using Extrusion = Extrusion_t<GridGeometry>;
57 using GridView = typename GridGeometry::GridView;
58 using Element = typename GridView::template Codim<0>::Entity;
60 using CellCenterPrimaryVariables = GetPropType<TypeTag, Properties::CellCenterPrimaryVariables>;
63
65
66 static const int numComponents = ModelTraits::numFluidComponents();
67 static const int numPhases = ModelTraits::numFluidPhases();
68 static constexpr bool useMoles = getPropValue<TypeTag, Properties::UseMoles>();
69
70 using ReducedComponentVector = Dune::FieldVector<Scalar, numComponents-1>;
71 using ReducedComponentMatrix = Dune::FieldMatrix<Scalar, numComponents-1, numComponents-1>;
72
73 static_assert(ModelTraits::numFluidPhases() == 1, "Only one phase allowed supported!");
74
75 static_assert(referenceSystem == ReferenceSystemFormulation::massAveraged, "only the mass averaged reference system is supported for the Maxwell-Stefan formulation");
76
77public:
78 // state the discretization method this implementation belongs to
80 //return the reference system
82 { return referenceSystem; }
83
88
90
96 template<class ElementVolumeVariables>
97 static CellCenterPrimaryVariables flux(const Problem& problem,
98 const Element& element,
99 const FVElementGeometry& fvGeometry,
100 const ElementVolumeVariables& elemVolVars,
101 const SubControlVolumeFace& scvf)
102 {
103 //this is to calculate the maxwellStefan diffusion in a multicomponent system.
104 //see: Multicomponent Mass Transfer. R. Taylor u. R. Krishna. J. Wiley & Sons, New York 1993
105 CellCenterPrimaryVariables componentFlux(0.0);
106 ReducedComponentVector moleFracInside(0.0);
107 ReducedComponentVector moleFracOutside(0.0);
108 ReducedComponentVector reducedFlux(0.0);
109 ReducedComponentMatrix reducedDiffusionMatrixInside(0.0);
110 ReducedComponentMatrix reducedDiffusionMatrixOutside(0.0);
111
112 // get inside/outside volume variables
113 const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
114 const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
115 const auto rhoInside = insideVolVars.density();
116 const auto rhoOutside = outsideVolVars.density();
117
118 //to implement outflow boundaries correctly we need to loop over all components but the main component as only for the transported ones we implement the outflow boundary. diffusion then is 0.
119 for(int compIdx = 0; compIdx < numComponents; ++compIdx)
120 {
121 if(compIdx == FluidSystem::getMainComponent(0))
122 continue;
123
124 // get equation index
125 const auto eqIdx = Indices::conti0EqIdx + compIdx;
126 if(scvf.boundary())
127 {
128 const auto bcTypes = problem.boundaryTypes(element, scvf);
129 if((bcTypes.isOutflow(eqIdx)) || (bcTypes.isSymmetry()))
130 return componentFlux;
131 }
132 }
133
134 //calculate the mole fraction vectors
135 for (int compIdx = 0; compIdx < numComponents-1; compIdx++)
136 {
137 //calculate x_inside
138 const auto xInside = insideVolVars.moleFraction(compIdx);
139 //calculate outside molefraction with the respective transmissibility
140 const auto xOutside = outsideVolVars.moleFraction(compIdx);
141
142 moleFracInside[compIdx] = xInside;
143 moleFracOutside[compIdx] = xOutside;
144 }
145
146 //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
147 reducedDiffusionMatrixInside = setupMSMatrix_(problem, fvGeometry, insideVolVars, scvf);
148
149 reducedDiffusionMatrixOutside = setupMSMatrix_(problem, fvGeometry, outsideVolVars, scvf);
150
151 const auto insideScvIdx = scvf.insideScvIdx();
152 const auto& insideScv = fvGeometry.scv(insideScvIdx);
153 const auto outsideScvIdx = scvf.outsideScvIdx();
154 const auto& outsideScv = fvGeometry.scv(outsideScvIdx);
155
156 const Scalar insideDistance = (insideScv.dofPosition() - scvf.ipGlobal()).two_norm();
157
158 const Scalar omegai = calculateOmega_(insideDistance, insideVolVars.extrusionFactor());
159
160 //if on boundary
161 if (scvf.boundary())
162 {
163 moleFracOutside -= moleFracInside;
164 reducedDiffusionMatrixInside.solve(reducedFlux, moleFracOutside);
165 reducedFlux *= omegai*rhoInside;
166 }
167 else
168 {
169 const Scalar outsideDistance = (outsideScv.dofPosition() - scvf.ipGlobal()).two_norm();
170 const Scalar omegaj = calculateOmega_(outsideDistance, outsideVolVars.extrusionFactor());
171
172 reducedDiffusionMatrixInside.invert();
173 reducedDiffusionMatrixOutside.invert();
174 reducedDiffusionMatrixInside *= omegai*rhoInside;
175 reducedDiffusionMatrixOutside *= omegaj*rhoOutside;
176
177 //in the helpervector we store the values for x*
178 ReducedComponentVector helperVector(0.0);
179 ReducedComponentVector gradientVectori(0.0);
180 ReducedComponentVector gradientVectorj(0.0);
181
182 reducedDiffusionMatrixInside.mv(moleFracInside, gradientVectori);
183 reducedDiffusionMatrixOutside.mv(moleFracOutside, gradientVectorj);
184
185 auto gradientVectorij = (gradientVectori + gradientVectorj);
186
187 //add the two matrixes to each other
188 reducedDiffusionMatrixOutside += reducedDiffusionMatrixInside;
189
190 reducedDiffusionMatrixOutside.solve(helperVector, gradientVectorij);
191
192 //Bi^-1 omegai rho(x*-xi)
193 helperVector -=moleFracInside;
194 reducedDiffusionMatrixInside.mv(helperVector, reducedFlux);
195 }
196
197 reducedFlux *= -Extrusion::area(scvf);
198
199 for (int compIdx = 0; compIdx < numComponents-1; compIdx++)
200 {
201 componentFlux[compIdx] = reducedFlux[compIdx];
202 componentFlux[numComponents-1] -= reducedFlux[compIdx];
203 }
204
205 return componentFlux ;
206 }
207
208private:
209 static Scalar calculateOmega_(const Scalar distance,
210 const Scalar extrusionFactor)
211 {
212 Scalar omega = 1/distance;
213 omega *= extrusionFactor;
214
215 return omega;
216 }
217
218 static ReducedComponentMatrix setupMSMatrix_(const Problem& problem,
219 const FVElementGeometry& fvGeometry,
220 const VolumeVariables& volVars,
221 const SubControlVolumeFace& scvf)
222 {
223 ReducedComponentMatrix reducedDiffusionMatrix(0.0);
224
225 for (int compIIdx = 0; compIIdx < numComponents-1; compIIdx++)
226 {
227 const auto xi = volVars.moleFraction(compIIdx);
228 const auto avgMolarMass = volVars.averageMolarMass(0);
229 const auto Mn = FluidSystem::molarMass(numComponents-1);
230 const Scalar tin = getEffectiveDiffusionCoefficient_(volVars, compIIdx, numComponents-1);
231
232 // set the entries of the diffusion matrix of the diagonal
233 reducedDiffusionMatrix[compIIdx][compIIdx] += xi*avgMolarMass/(tin*Mn);
234
235 for (int compJIdx = 0; compJIdx < numComponents; compJIdx++)
236 {
237 // we don't want to calculate e.g. water in water diffusion
238 if (compJIdx == compIIdx)
239 continue;
240
241 const auto xj = volVars.moleFraction(compJIdx);
242 const auto Mi = FluidSystem::molarMass(compIIdx);
243 const auto Mj = FluidSystem::molarMass(compJIdx);
244 const Scalar tij = getEffectiveDiffusionCoefficient_(volVars, compIIdx, compJIdx);
245 reducedDiffusionMatrix[compIIdx][compIIdx] += xj*avgMolarMass/(tij*Mi);
246 if (compJIdx < numComponents-1)
247 reducedDiffusionMatrix[compIIdx][compJIdx] += xi*(avgMolarMass/(tin*Mn) - avgMolarMass/(tij*Mj));
248 }
249 }
250 return reducedDiffusionMatrix;
251 }
252
253 static Scalar getEffectiveDiffusionCoefficient_(const VolumeVariables& volVars, const int phaseIdx, const int compIdx)
254 {
255 return volVars.effectiveDiffusionCoefficient(phaseIdx, FluidSystem::getMainComponent(phaseIdx), compIdx);
256 }
257};
258} // end namespace
259
260#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...
Classes related to flux variables caching.
Container storing the diffusion coefficients required by the Maxwell- Stefan diffusion law....
ctype distance(const Dune::FieldVector< ctype, dimWorld > &a, const Dune::FieldVector< ctype, dimWorld > &b)
Compute the shortest distance between two points.
Definition: geometry/distance.hh:138
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 (equivalent to old macro GET_PROP_TYPE(....
Definition: propertysystem.hh:149
Definition: maxwellstefanslaw.hh:37
The empty filler class corresponding to EmptyCache.
Definition: fluxvariablescaching.hh:32
Definition: fluxvariablescaching.hh:64
Container storing the diffusion coefficients required by the Maxwell- Stefan diffusion law....
Definition: maxwellstefandiffusioncoefficients.hh:45
static CellCenterPrimaryVariables flux(const Problem &problem, const Element &element, const FVElementGeometry &fvGeometry, const ElementVolumeVariables &elemVolVars, const SubControlVolumeFace &scvf)
Returns the diffusive fluxes of all components within a fluid phase across the given sub-control volu...
Definition: staggered/freeflow/maxwellstefanslaw.hh:97
static constexpr ReferenceSystemFormulation referenceSystemFormulation()
Definition: staggered/freeflow/maxwellstefanslaw.hh:81
Declares all properties used in Dumux.