GCC Code Coverage Report


Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/test/porousmediumflow/mpnc/thermalnonequilibrium/problem.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 0 98 0.0%
Functions: 0 6 0.0%
Branches: 0 130 0.0%

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