25#ifndef DUMUX_DISCRETIZATION_CC_TPFA_MAXWELL_STEFAN_LAW_HH
26#define DUMUX_DISCRETIZATION_CC_TPFA_MAXWELL_STEFAN_LAW_HH
28#include <dune/common/fmatrix.hh>
29#include <dune/common/float_cmp.hh>
43template <
class TypeTag, DiscretizationMethod discMethod, ReferenceSystemFormulation referenceSystem>
44class MaxwellStefansLawImplementation;
50template <
class TypeTag, ReferenceSystemFormulation referenceSystem>
56 using FVElementGeometry =
typename GridGeometry::LocalView;
57 using SubControlVolume =
typename GridGeometry::SubControlVolume;
58 using SubControlVolumeFace =
typename GridGeometry::SubControlVolumeFace;
60 using GridView =
typename GridGeometry::GridView;
63 using Element =
typename GridView::template Codim<0>::Entity;
68 static const int dim = GridView::dimension;
69 static const int dimWorld = GridView::dimensionworld;
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>;
85 {
return referenceSystem; }
99 static ComponentFluxVector
flux(
const Problem& problem,
100 const Element& element,
101 const FVElementGeometry& fvGeometry,
102 const ElementVolumeVariables& elemVolVars,
103 const SubControlVolumeFace& scvf,
105 const ElementFluxVariablesCache& elemFluxVarsCache)
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);
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);
122 for (
int compIdx = 0; compIdx < numComponents-1; compIdx++)
125 const auto xInside = insideVolVars.moleFraction(phaseIdx, compIdx);
127 const auto xOutside = outsideVolVars.moleFraction(phaseIdx, compIdx);
129 moleFracInside[compIdx] = xInside;
130 moleFracOutside[compIdx] = xOutside;
134 if(!(Dune::FloatCmp::eq<Scalar>(insideVolVars.saturation(phaseIdx), 0) || Dune::FloatCmp::eq<Scalar>(outsideVolVars.saturation(phaseIdx), 0)))
136 const auto insideScvIdx = scvf.insideScvIdx();
137 const auto& insideScv = fvGeometry.scv(insideScvIdx);
138 const Scalar omegai = calculateOmega_(scvf,
140 insideVolVars.extrusionFactor());
143 reducedDiffusionMatrixInside = setupMSMatrix_(problem, element, fvGeometry, insideVolVars, insideScv, phaseIdx);
146 if (scvf.boundary() || scvf.numOutsideScvs() > 1)
148 moleFracOutside -= moleFracInside;
149 reducedDiffusionMatrixInside.solve(reducedFlux, moleFracOutside);
150 reducedFlux *= omegai*rhoInside;
155 const auto outsideScvIdx = scvf.outsideScvIdx();
156 const auto& outsideScv = fvGeometry.scv(outsideScvIdx);
158 reducedDiffusionMatrixOutside = setupMSMatrix_(problem, element, fvGeometry, outsideVolVars, outsideScv, phaseIdx);
162 omegaj = -1.0*calculateOmega_(scvf,
164 outsideVolVars.extrusionFactor());
166 omegaj = calculateOmega_(fvGeometry.flipScvf(scvf.index()),
168 outsideVolVars.extrusionFactor());
170 reducedDiffusionMatrixInside.invert();
171 reducedDiffusionMatrixOutside.invert();
172 reducedDiffusionMatrixInside *= omegai*rhoInside;
173 reducedDiffusionMatrixOutside *= omegaj*rhoOutside;
176 ReducedComponentVector helperVector(0.0);
177 ReducedComponentVector gradientVectori(0.0);
178 ReducedComponentVector gradientVectorj(0.0);
180 reducedDiffusionMatrixInside.mv(moleFracInside, gradientVectori);
181 reducedDiffusionMatrixOutside.mv(moleFracOutside, gradientVectorj);
183 auto gradientVectorij = (gradientVectori + gradientVectorj);
186 reducedDiffusionMatrixOutside += reducedDiffusionMatrixInside;
188 reducedDiffusionMatrixOutside.solve(helperVector, gradientVectorij);
191 helperVector -=moleFracInside;
192 reducedDiffusionMatrixInside.mv(helperVector, reducedFlux);
195 reducedFlux *= -Extrusion::area(scvf);
196 for (
int compIdx = 0; compIdx < numComponents-1; compIdx++)
198 componentFlux[compIdx] = reducedFlux[compIdx];
199 componentFlux[numComponents-1] -=reducedFlux[compIdx];
202 return componentFlux ;
206 static Scalar calculateOmega_(
const SubControlVolumeFace& scvf,
207 const SubControlVolume &scv,
208 const Scalar extrusionFactor)
210 auto distanceVector = scvf.ipGlobal();
211 distanceVector -= scv.center();
212 distanceVector /= distanceVector.two_norm2();
214 Scalar omega = (distanceVector * scvf.unitOuterNormal());
215 omega *= extrusionFactor;
220 static ReducedComponentMatrix setupMSMatrix_(
const Problem& problem,
221 const Element& element,
222 const FVElementGeometry& fvGeometry,
223 const VolumeVariables& volVars,
224 const SubControlVolume& scv,
227 ReducedComponentMatrix reducedDiffusionMatrix(0.0);
230 if(Dune::FloatCmp::eq<Scalar>(volVars.saturation(phaseIdx), 0))
231 return reducedDiffusionMatrix;
233 const auto avgMolarMass = volVars.averageMolarMass(phaseIdx);
235 for (
int compIIdx = 0; compIIdx < numComponents-1; compIIdx++)
237 const auto xi = volVars.moleFraction(phaseIdx, compIIdx);
239 const auto Mn = FluidSystem::molarMass(numComponents-1);
240 Scalar tin = volVars.effectiveDiffusionCoefficient(phaseIdx, compIIdx, numComponents-1);
243 reducedDiffusionMatrix[compIIdx][compIIdx] += xi*avgMolarMass/(tin*Mn);
246 for (
int compJIdx = 0; compJIdx < numComponents; compJIdx++)
249 if (compJIdx == compIIdx)
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);
258 if (compJIdx < numComponents-1)
259 reducedDiffusionMatrix[compIIdx][compJIdx] += xi*(avgMolarMass/(tin*Mn) - avgMolarMass/(tij*Mj));
262 return reducedDiffusionMatrix;
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
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.