Line |
Branch |
Exec |
Source |
1 |
|
|
// |
2 |
|
|
// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder |
3 |
|
|
// SPDX-License-Identifier: GPL-3.0-or-later |
4 |
|
|
// |
5 |
|
|
/*! |
6 |
|
|
* \file |
7 |
|
|
* \ingroup MPNCTests |
8 |
|
|
* \brief Problem where hot, pure liquid water is injected from the left hand side into an initially |
9 |
|
|
* isothermal domain. |
10 |
|
|
* |
11 |
|
|
* The water is fully evaporated by a strong heat source. |
12 |
|
|
* A local thermal non-equilibrium model is used: i.e. two different (fluid, solid) |
13 |
|
|
* temperatures are primary variables. |
14 |
|
|
* |
15 |
|
|
* \author Philipp Nuske |
16 |
|
|
*/ |
17 |
|
|
|
18 |
|
|
#ifndef DUMUX_COMBUSTION_PROBLEM_ONE_COMPONENT_HH |
19 |
|
|
#define DUMUX_COMBUSTION_PROBLEM_ONE_COMPONENT_HH |
20 |
|
|
|
21 |
|
|
#include <dumux/common/properties.hh> |
22 |
|
|
#include <dumux/common/parameters.hh> |
23 |
|
|
#include <dumux/common/boundarytypes.hh> |
24 |
|
|
#include <dumux/common/numeqvector.hh> |
25 |
|
|
|
26 |
|
|
#include <dumux/porousmediumflow/problem.hh> |
27 |
|
|
#include <dumux/porousmediumflow/mpnc/initialconditionhelper.hh> |
28 |
|
|
|
29 |
|
|
namespace Dumux { |
30 |
|
|
|
31 |
|
|
/*! |
32 |
|
|
* \ingroup MPNCTests |
33 |
|
|
* \brief Problem where water is injected from the left hand side into a porous media filled domain, |
34 |
|
|
* which is supplied with energy from the right hand side to evaporate the water. |
35 |
|
|
*/ |
36 |
|
|
template<class TypeTag> |
37 |
|
|
class CombustionProblemOneComponent: public PorousMediumFlowProblem<TypeTag> |
38 |
|
|
{ |
39 |
|
|
using ParentType = PorousMediumFlowProblem<TypeTag>; |
40 |
|
|
using Scalar = GetPropType<TypeTag, Properties::Scalar>; |
41 |
|
|
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>; |
42 |
|
|
using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>; |
43 |
|
|
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>; |
44 |
|
|
using NumEqVector = Dumux::NumEqVector<PrimaryVariables>; |
45 |
|
|
|
46 |
|
|
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>; |
47 |
|
|
using ElementVolumeVariables = typename GridVariables::GridVolumeVariables::LocalView; |
48 |
|
|
using ElementFluxVariablesCache = typename GridVariables::GridFluxVariablesCache::LocalView; |
49 |
|
|
|
50 |
|
|
using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView; |
51 |
|
|
using SubControlVolume = typename FVElementGeometry::SubControlVolume; |
52 |
|
|
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; |
53 |
|
|
using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView; |
54 |
|
|
using Element = typename GridView::template Codim<0>::Entity; |
55 |
|
|
using GlobalPosition = typename Element::Geometry::GlobalCoordinate; |
56 |
|
|
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>; |
57 |
|
|
using FluidState = GetPropType<TypeTag, Properties::FluidState>; |
58 |
|
|
using ParameterCache = typename FluidSystem::ParameterCache; |
59 |
|
|
|
60 |
|
|
using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>; |
61 |
|
|
using Indices = typename ModelTraits::Indices; |
62 |
|
|
|
63 |
|
|
static constexpr int dimWorld = GridView::dimensionworld; |
64 |
|
|
static constexpr int numComponents = ModelTraits::numFluidComponents(); |
65 |
|
|
static constexpr int s0Idx = Indices::s0Idx; |
66 |
|
|
static constexpr int p0Idx = Indices::p0Idx; |
67 |
|
|
static constexpr int conti00EqIdx = Indices::conti0EqIdx; |
68 |
|
|
static constexpr int energyEq0Idx = Indices::energyEqIdx; |
69 |
|
|
static constexpr int numEnergyEqFluid = ModelTraits::numEnergyEqFluid(); |
70 |
|
|
static constexpr int numEnergyEqSolid = ModelTraits::numEnergyEqSolid(); |
71 |
|
|
static constexpr int energyEqSolidIdx = energyEq0Idx + numEnergyEqFluid + numEnergyEqSolid - 1; |
72 |
|
|
static constexpr int wPhaseIdx = FluidSystem::wPhaseIdx; |
73 |
|
|
static constexpr int nPhaseIdx = FluidSystem::nPhaseIdx; |
74 |
|
|
static constexpr int wCompIdx = FluidSystem::H2OIdx; |
75 |
|
|
static constexpr int nCompIdx = FluidSystem::N2Idx; |
76 |
|
|
|
77 |
|
|
static constexpr auto numPhases = ModelTraits::numFluidPhases(); |
78 |
|
|
|
79 |
|
|
public: |
80 |
|
✗ |
CombustionProblemOneComponent(std::shared_ptr<const GridGeometry> gridGeometry) |
81 |
|
✗ |
: ParentType(gridGeometry) |
82 |
|
|
{ |
83 |
|
✗ |
outputName_ = getParam<std::string>("Problem.Name"); |
84 |
|
✗ |
nRestart_ = getParam<Scalar>("Constants.nRestart"); |
85 |
|
✗ |
TInitial_ = getParam<Scalar>("InitialConditions.TInitial"); |
86 |
|
✗ |
TRight_ = getParam<Scalar>("InitialConditions.TRight"); |
87 |
|
✗ |
pnInitial_ = getParam<Scalar>("InitialConditions.pnInitial"); |
88 |
|
✗ |
TBoundary_ = getParam<Scalar>("BoundaryConditions.TBoundary"); |
89 |
|
✗ |
SwBoundary_ = getParam<Scalar>("BoundaryConditions.SwBoundary"); |
90 |
|
✗ |
SwOneComponentSys_= getParam<Scalar>("BoundaryConditions.SwOneComponentSys"); |
91 |
|
✗ |
massFluxInjectedPhase_ = getParam<Scalar>("BoundaryConditions.massFluxInjectedPhase"); |
92 |
|
✗ |
heatFluxFromRight_ = getParam<Scalar>("BoundaryConditions.heatFluxFromRight"); |
93 |
|
✗ |
coldTime_ =getParam<Scalar>("BoundaryConditions.coldTime"); |
94 |
|
✗ |
time_ = 0.0; |
95 |
|
✗ |
} |
96 |
|
|
|
97 |
|
✗ |
void setTime(Scalar time) |
98 |
|
✗ |
{ time_ = time; } |
99 |
|
|
|
100 |
|
|
void setGridVariables(std::shared_ptr<GridVariables> gridVariables) |
101 |
|
✗ |
{ gridVariables_ = gridVariables; } |
102 |
|
|
|
103 |
|
|
const GridVariables& gridVariables() const |
104 |
|
✗ |
{ return *gridVariables_; } |
105 |
|
|
|
106 |
|
|
/*! |
107 |
|
|
* \name Problem Params |
108 |
|
|
*/ |
109 |
|
|
// \{ |
110 |
|
|
/*! |
111 |
|
|
* \brief The problem name. |
112 |
|
|
* |
113 |
|
|
* This is used as a prefix for files generated by the simulation. |
114 |
|
|
*/ |
115 |
|
|
const std::string name() const |
116 |
|
✗ |
{ return outputName_;} |
117 |
|
|
|
118 |
|
|
// \} |
119 |
|
|
|
120 |
|
|
/*! |
121 |
|
|
* \brief Evaluates the source term for all balance equations within a given |
122 |
|
|
* sub-control volume. |
123 |
|
|
* |
124 |
|
|
* \param element The finite element |
125 |
|
|
* \param fvGeometry The finite volume geometry of the element |
126 |
|
|
* \param elemVolVars The volume variables of the element |
127 |
|
|
* \param scv The local index of the sub-control volume |
128 |
|
|
* |
129 |
|
|
* Positive values mean that mass is created, negative ones mean that it vanishes. |
130 |
|
|
*/ |
131 |
|
✗ |
NumEqVector source(const Element &element, |
132 |
|
|
const FVElementGeometry& fvGeometry, |
133 |
|
|
const ElementVolumeVariables& elemVolVars, |
134 |
|
|
const SubControlVolume &scv) const |
135 |
|
|
{ |
136 |
|
✗ |
NumEqVector values(0.0); |
137 |
|
|
|
138 |
|
✗ |
const auto& globalPos = scv.dofPosition(); |
139 |
|
|
|
140 |
|
✗ |
const Scalar volume = scv.volume(); |
141 |
|
✗ |
const Scalar numScv = fvGeometry.numScv(); // box: numSCV, cc:1 |
142 |
|
|
|
143 |
|
✗ |
if (time_ > coldTime_ ) |
144 |
|
|
{ |
145 |
|
✗ |
if (onRightBoundaryPorousMedium_(globalPos)) |
146 |
|
|
{ |
147 |
|
|
// Testing the location of a vertex, but function is called for each associated scv. Compensate for that |
148 |
|
✗ |
values[energyEqSolidIdx] = heatFluxFromRight_ / volume / numScv; |
149 |
|
|
} |
150 |
|
|
} |
151 |
|
✗ |
return values; |
152 |
|
|
} |
153 |
|
|
|
154 |
|
|
/*! |
155 |
|
|
* \brief Specifies which kind of boundary condition should be |
156 |
|
|
* used for which equation on a given boundary segment. |
157 |
|
|
* |
158 |
|
|
* \param globalPos The position for which the bc type should be evaluated |
159 |
|
|
*/ |
160 |
|
✗ |
BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const |
161 |
|
|
{ |
162 |
|
✗ |
BoundaryTypes bcTypes; |
163 |
|
|
// Default: Neumann |
164 |
|
✗ |
bcTypes.setAllNeumann(); |
165 |
|
|
|
166 |
|
✗ |
if(onRightBoundary_(globalPos) ) { |
167 |
|
|
bcTypes.setAllDirichlet(); |
168 |
|
|
} |
169 |
|
✗ |
return bcTypes; |
170 |
|
|
} |
171 |
|
|
|
172 |
|
|
/*! |
173 |
|
|
* \brief Evaluates the boundary conditions for a Dirichlet boundary segment. |
174 |
|
|
* |
175 |
|
|
* \param globalPos The global position |
176 |
|
|
* |
177 |
|
|
* For this method, the \a values parameter stores primary variables. |
178 |
|
|
*/ |
179 |
|
|
PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const |
180 |
|
|
{ |
181 |
|
✗ |
return initial_(globalPos); |
182 |
|
|
} |
183 |
|
|
|
184 |
|
|
/*! |
185 |
|
|
* \brief Evaluates the boundary conditions for a Neumann boundary segment. |
186 |
|
|
* |
187 |
|
|
* \param element The finite element |
188 |
|
|
* \param fvGeometry The finite-volume geometry in the box scheme |
189 |
|
|
* \param elemVolVars The volume variables of the element |
190 |
|
|
* \param elemFluxVarsCache Flux variables caches for all faces in stencil |
191 |
|
|
* \param scvf The sub-control volume face |
192 |
|
|
* |
193 |
|
|
* For this method, the \a values parameter stores the mass flux |
194 |
|
|
* in normal direction of each phase. Negative values mean influx. |
195 |
|
|
*/ |
196 |
|
✗ |
NumEqVector neumann(const Element &element, |
197 |
|
|
const FVElementGeometry& fvGeometry, |
198 |
|
|
const ElementVolumeVariables& elemVolVars, |
199 |
|
|
const ElementFluxVariablesCache& elemFluxVarsCache, |
200 |
|
|
const SubControlVolumeFace& scvf) const |
201 |
|
|
{ |
202 |
|
✗ |
NumEqVector values(0.0); |
203 |
|
|
|
204 |
|
✗ |
const auto& globalPos = fvGeometry.scv(scvf.insideScvIdx()).dofPosition(); |
205 |
|
✗ |
const auto& scvIdx = scvf.insideScvIdx(); |
206 |
|
✗ |
const Scalar massFluxInjectedPhase = massFluxInjectedPhase_; |
207 |
|
|
|
208 |
|
✗ |
FluidState fluidState; |
209 |
|
|
|
210 |
|
✗ |
const Scalar pn = elemVolVars[scvIdx].pressure(nPhaseIdx); |
211 |
|
✗ |
const Scalar pw = elemVolVars[scvIdx].pressure(wPhaseIdx); |
212 |
|
|
|
213 |
|
✗ |
fluidState.setPressure(nPhaseIdx, pn); |
214 |
|
✗ |
fluidState.setPressure(wPhaseIdx, pw); |
215 |
|
|
|
216 |
|
✗ |
fluidState.setTemperature(TBoundary_); |
217 |
|
|
ParameterCache dummyCache; |
218 |
|
✗ |
fluidState.setMoleFraction(wPhaseIdx, nCompIdx, 0.0); |
219 |
|
✗ |
fluidState.setMoleFraction(wPhaseIdx, wCompIdx, 1.0); |
220 |
|
|
// compute density of injection phase |
221 |
|
✗ |
const Scalar density = FluidSystem::density(fluidState, |
222 |
|
|
dummyCache, |
223 |
|
|
wPhaseIdx); |
224 |
|
✗ |
fluidState.setDensity(wPhaseIdx, density); |
225 |
|
✗ |
const Scalar molarDensity = FluidSystem::molarDensity(fluidState, |
226 |
|
|
dummyCache, |
227 |
|
|
wPhaseIdx); |
228 |
|
✗ |
fluidState.setMolarDensity(wPhaseIdx, molarDensity); |
229 |
|
|
|
230 |
|
✗ |
for(int phaseIdx=0; phaseIdx<numPhases; phaseIdx++) { |
231 |
|
✗ |
const Scalar h = FluidSystem::enthalpy(fluidState, |
232 |
|
|
dummyCache, |
233 |
|
|
phaseIdx); |
234 |
|
✗ |
fluidState.setEnthalpy(phaseIdx, h); |
235 |
|
|
} |
236 |
|
|
|
237 |
|
✗ |
const Scalar molarFlux = massFluxInjectedPhase / fluidState.averageMolarMass(wPhaseIdx); |
238 |
|
|
|
239 |
|
✗ |
if (onLeftBoundary_(globalPos)) |
240 |
|
|
{ |
241 |
|
✗ |
values[conti00EqIdx + wCompIdx] = - molarFlux * fluidState.moleFraction(wPhaseIdx, wCompIdx); |
242 |
|
✗ |
values[conti00EqIdx + nCompIdx] = - molarFlux * fluidState.moleFraction(wPhaseIdx, nCompIdx); |
243 |
|
✗ |
values[energyEq0Idx] = - massFluxInjectedPhase * fluidState.enthalpy(wPhaseIdx); |
244 |
|
|
} |
245 |
|
✗ |
return values; |
246 |
|
|
} |
247 |
|
|
|
248 |
|
|
/*! |
249 |
|
|
* \brief Evaluates the initial value for a control volume. |
250 |
|
|
* |
251 |
|
|
* For this method, the \a values parameter stores primary |
252 |
|
|
* variables. |
253 |
|
|
* |
254 |
|
|
* \param globalPos The global position |
255 |
|
|
*/ |
256 |
|
|
PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const |
257 |
|
|
{ |
258 |
|
✗ |
return initial_(globalPos); |
259 |
|
|
} |
260 |
|
|
|
261 |
|
|
private: |
262 |
|
|
// the internal method for the initial condition |
263 |
|
✗ |
PrimaryVariables initial_(const GlobalPosition &globalPos) const |
264 |
|
|
{ |
265 |
|
✗ |
const Scalar curPos = globalPos[0]; |
266 |
|
✗ |
const Scalar slope = (SwBoundary_ - SwOneComponentSys_)/this->spatialParams().lengthPM(); |
267 |
|
|
Scalar S[numPhases]; |
268 |
|
✗ |
const Scalar thisSaturation = SwOneComponentSys_ + curPos * slope; |
269 |
|
|
|
270 |
|
✗ |
S[wPhaseIdx] = SwBoundary_; |
271 |
|
✗ |
if (inPM_(globalPos) ) |
272 |
|
✗ |
S[wPhaseIdx] = thisSaturation; |
273 |
|
✗ |
S[nPhaseIdx] = 1. - S[wPhaseIdx]; |
274 |
|
|
|
275 |
|
✗ |
Scalar thisTemperature = TInitial_; |
276 |
|
✗ |
if(onRightBoundary_(globalPos)) |
277 |
|
✗ |
thisTemperature = TRight_; |
278 |
|
|
|
279 |
|
✗ |
FluidState fs; |
280 |
|
✗ |
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) |
281 |
|
|
{ |
282 |
|
✗ |
fs.setSaturation(phaseIdx, S[phaseIdx]); |
283 |
|
✗ |
fs.setTemperature(thisTemperature); |
284 |
|
|
} |
285 |
|
|
|
286 |
|
|
//obtain pc according to saturation |
287 |
|
✗ |
const int wettingPhaseIdx = this->spatialParams().template wettingPhaseAtPos<FluidSystem>(globalPos); |
288 |
|
✗ |
const auto& fm = this->spatialParams().fluidMatrixInteractionAtPos(globalPos); |
289 |
|
✗ |
const auto capPress = fm.capillaryPressures(fs, wettingPhaseIdx); |
290 |
|
|
Scalar p[numPhases]; |
291 |
|
|
|
292 |
|
|
using std::abs; |
293 |
|
✗ |
p[wPhaseIdx] = pnInitial_ - abs(capPress[wPhaseIdx]); |
294 |
|
✗ |
p[nPhaseIdx] = p[wPhaseIdx] + abs(capPress[wPhaseIdx]); |
295 |
|
✗ |
for (int phaseIdx = 0; phaseIdx < numPhases; phaseIdx++) |
296 |
|
✗ |
fs.setPressure(phaseIdx, p[phaseIdx]); |
297 |
|
|
|
298 |
|
✗ |
fs.setMoleFraction(wPhaseIdx, wCompIdx, 1.0); |
299 |
|
✗ |
fs.setMoleFraction(wPhaseIdx, nCompIdx, 0.0); |
300 |
|
✗ |
fs.setMoleFraction(nPhaseIdx, wCompIdx, 1.0); |
301 |
|
✗ |
fs.setMoleFraction(nPhaseIdx, nCompIdx, 0.0); |
302 |
|
|
|
303 |
|
✗ |
const int refPhaseIdx = inPM_(globalPos) ? wPhaseIdx : nPhaseIdx; |
304 |
|
|
|
305 |
|
|
MPNCInitialConditionHelper<PrimaryVariables, FluidSystem, ModelTraits> helper; |
306 |
|
✗ |
helper.solve(fs, MPNCInitialConditions::NotAllPhasesPresent{.refPhaseIdx = refPhaseIdx}); |
307 |
|
✗ |
auto priVars = helper.getPrimaryVariables(fs); |
308 |
|
|
|
309 |
|
|
// additionally set the temperature for thermal non-equilibrium for the phases |
310 |
|
✗ |
priVars[energyEq0Idx] = thisTemperature; |
311 |
|
✗ |
priVars[energyEqSolidIdx] = thisTemperature; |
312 |
|
✗ |
return priVars; |
313 |
|
|
} |
314 |
|
|
|
315 |
|
|
/*! |
316 |
|
|
* \brief Returns whether the tested position is on the left boundary of the domain. |
317 |
|
|
*/ |
318 |
|
|
bool onLeftBoundary_(const GlobalPosition & globalPos) const |
319 |
|
✗ |
{ return globalPos[0] < this->gridGeometry().bBoxMin()[0] + eps_;} |
320 |
|
|
|
321 |
|
|
/*! |
322 |
|
|
* \brief Returns whether the tested position is on the right boundary of the domain. |
323 |
|
|
*/ |
324 |
|
|
bool onRightBoundary_(const GlobalPosition & globalPos) const |
325 |
|
✗ |
{ return globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_;} |
326 |
|
|
|
327 |
|
|
/*! |
328 |
|
|
* \brief Returns whether the tested position is in a specific region (right) in the domain |
329 |
|
|
* \todo this needs to be more sophisticated in order to allow for meshes with nodes not directly on the boundary |
330 |
|
|
*/ |
331 |
|
|
bool onRightBoundaryPorousMedium_(const GlobalPosition & globalPos) const |
332 |
|
|
{ |
333 |
|
|
using std::abs; |
334 |
|
✗ |
return (abs(globalPos[0] - (this->spatialParams().lengthPM())) < eps_ ); |
335 |
|
|
} |
336 |
|
|
|
337 |
|
|
/*! |
338 |
|
|
* \brief Returns whether the tested position is in the porous medium. |
339 |
|
|
*/ |
340 |
|
|
bool inPM_(const GlobalPosition & globalPos) const |
341 |
|
✗ |
{ return !this->spatialParams().inOutFlow(globalPos); } |
342 |
|
|
|
343 |
|
|
private: |
344 |
|
|
static constexpr Scalar eps_ = 1e-6; |
345 |
|
|
int nTemperature_; |
346 |
|
|
int nPressure_; |
347 |
|
|
std::string outputName_; |
348 |
|
|
int nRestart_; |
349 |
|
|
Scalar TInitial_; |
350 |
|
|
Scalar TRight_; |
351 |
|
|
|
352 |
|
|
Scalar pnInitial_; |
353 |
|
|
|
354 |
|
|
Dune::ParameterTree inputParameters_; |
355 |
|
|
Scalar x_[numPhases][numComponents]; |
356 |
|
|
|
357 |
|
|
Scalar TBoundary_; |
358 |
|
|
Scalar SwBoundary_; |
359 |
|
|
Scalar SwOneComponentSys_; |
360 |
|
|
|
361 |
|
|
Scalar massFluxInjectedPhase_; |
362 |
|
|
Scalar heatFluxFromRight_; |
363 |
|
|
Scalar lengthPM_; |
364 |
|
|
Scalar coldTime_; |
365 |
|
|
|
366 |
|
|
Scalar time_; |
367 |
|
|
std::shared_ptr<GridVariables> gridVariables_; |
368 |
|
|
}; |
369 |
|
|
|
370 |
|
|
} // end namespace Dumux |
371 |
|
|
|
372 |
|
|
#endif |
373 |
|
|
|