| Line | Branch | Exec | Source |
|---|---|---|---|
| 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 | // SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder | ||
| 5 | // SPDX-License-Identifier: GPL-3.0-or-later | ||
| 6 | // | ||
| 7 | /*! | ||
| 8 | * \file | ||
| 9 | * \ingroup MPNCTests | ||
| 10 | * \brief Problem showcasing the capabilities of the kinetic model. | ||
| 11 | * | ||
| 12 | * The whole domain is porous medium, but the upper half has properties mimicking | ||
| 13 | * the ones of a free-flow domain. | ||
| 14 | * This way a poor man's coupling approach is accomplished: Without the | ||
| 15 | * complications of coupling, the main characteristics of a porous and a free-flow | ||
| 16 | * domain are depicted. | ||
| 17 | * | ||
| 18 | * The porous domain is bypassed with dry air. This way the equilibration | ||
| 19 | * process on top of the porous domain can be studied. | ||
| 20 | * | ||
| 21 | * \author Philipp Nuske | ||
| 22 | */ | ||
| 23 | #ifndef DUMUX_EVAPORATION_ATMOSPHERE_PROBLEM_HH | ||
| 24 | #define DUMUX_EVAPORATION_ATMOSPHERE_PROBLEM_HH | ||
| 25 | |||
| 26 | #include <dumux/common/properties.hh> | ||
| 27 | #include <dumux/common/parameters.hh> | ||
| 28 | #include <dumux/common/boundarytypes.hh> | ||
| 29 | #include <dumux/common/numeqvector.hh> | ||
| 30 | |||
| 31 | #include <dumux/porousmediumflow/problem.hh> | ||
| 32 | #include <dumux/porousmediumflow/mpnc/pressureformulation.hh> | ||
| 33 | #include <dumux/material/binarycoefficients/h2o_n2.hh> | ||
| 34 | #include <dumux/material/constraintsolvers/misciblemultiphasecomposition.hh> | ||
| 35 | |||
| 36 | namespace Dumux { | ||
| 37 | |||
| 38 | /*! | ||
| 39 | * \ingroup MPNCTests | ||
| 40 | * | ||
| 41 | * \brief Problem that simulates the coupled heat and mass transfer processes | ||
| 42 | resulting from the evaporation of liquid water from | ||
| 43 | * a porous medium sub-domain into a gas filled "quasi-freeflow" sub-domain. | ||
| 44 | */ | ||
| 45 | template <class TypeTag> | ||
| 46 | class EvaporationAtmosphereProblem: public PorousMediumFlowProblem<TypeTag> | ||
| 47 | { | ||
| 48 | using ParentType = PorousMediumFlowProblem<TypeTag>; | ||
| 49 | using Scalar = GetPropType<TypeTag, Properties::Scalar>; | ||
| 50 | using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>; | ||
| 51 | using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>; | ||
| 52 | using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>; | ||
| 53 | using NumEqVector = Dumux::NumEqVector<PrimaryVariables>; | ||
| 54 | using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView; | ||
| 55 | using ElementFluxVariablesCache = typename GetPropType<TypeTag, Properties::GridFluxVariablesCache>::LocalView; | ||
| 56 | using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView; | ||
| 57 | using SubControlVolume = typename FVElementGeometry::SubControlVolume; | ||
| 58 | using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; | ||
| 59 | using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView; | ||
| 60 | using Element = typename GridView::template Codim<0>::Entity; | ||
| 61 | using GlobalPosition = typename Element::Geometry::GlobalCoordinate; | ||
| 62 | using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>; | ||
| 63 | using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>; | ||
| 64 | using FluidState = GetPropType<TypeTag, Properties::FluidState>; | ||
| 65 | using ParameterCache = typename FluidSystem::ParameterCache; | ||
| 66 | using GridVariables = GetPropType<TypeTag, Properties::GridVariables>; | ||
| 67 | |||
| 68 | using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>; | ||
| 69 | using Indices = typename ModelTraits::Indices; | ||
| 70 | |||
| 71 | enum { dimWorld = GridView::dimensionworld }; | ||
| 72 | enum { numPhases = ModelTraits::numFluidPhases() }; | ||
| 73 | enum { numComponents = ModelTraits::numFluidComponents() }; | ||
| 74 | enum { s0Idx = Indices::s0Idx }; | ||
| 75 | enum { p0Idx = Indices::p0Idx }; | ||
| 76 | enum { conti00EqIdx = Indices::conti0EqIdx }; | ||
| 77 | enum { energyEq0Idx = Indices::energyEqIdx }; | ||
| 78 | enum { liquidPhaseIdx = FluidSystem::liquidPhaseIdx }; | ||
| 79 | enum { gasPhaseIdx = FluidSystem::gasPhaseIdx }; | ||
| 80 | enum { wCompIdx = FluidSystem::H2OIdx }; | ||
| 81 | enum { nCompIdx = FluidSystem::N2Idx }; | ||
| 82 | enum { numEnergyEqFluid = ModelTraits::numEnergyEqFluid() }; | ||
| 83 | enum { numEnergyEqSolid = ModelTraits::numEnergyEqSolid() }; | ||
| 84 | |||
| 85 | static constexpr bool enableChemicalNonEquilibrium = getPropValue<TypeTag, Properties::EnableChemicalNonEquilibrium>(); | ||
| 86 | using ConstraintSolver = MiscibleMultiPhaseComposition<Scalar, FluidSystem>; | ||
| 87 | |||
| 88 | // formulations | ||
| 89 | static constexpr auto pressureFormulation = ModelTraits::pressureFormulation(); | ||
| 90 | static constexpr auto mostWettingFirst = MpNcPressureFormulation::mostWettingFirst; | ||
| 91 | static constexpr auto leastWettingFirst = MpNcPressureFormulation::leastWettingFirst; | ||
| 92 | |||
| 93 | public: | ||
| 94 | 1 | EvaporationAtmosphereProblem(std::shared_ptr<const GridGeometry> gridGeometry) | |
| 95 |
9/26✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 1 times.
✓ Branch 15 taken 1 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 1 times.
✗ Branch 19 not taken.
✓ Branch 21 taken 1 times.
✗ Branch 22 not taken.
✓ Branch 24 taken 1 times.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
|
3 | : ParentType(gridGeometry) |
| 96 | { | ||
| 97 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | percentOfEquil_ = getParam<Scalar>("BoundaryConditions.percentOfEquil"); |
| 98 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | nTemperature_ = getParam<Scalar>("FluidSystem.nTemperature"); |
| 99 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | nPressure_ = getParam<Scalar>("FluidSystem.nPressure"); |
| 100 |
2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
|
1 | outputName_ = getParam<std::string>("Problem.Name"); |
| 101 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | TInitial_ = getParam<Scalar>("InitialConditions.TInitial"); |
| 102 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | SwPMInitial_ = getParam<Scalar>("InitialConditions.SwPMInitial"); |
| 103 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | SwFFInitial_ = getParam<Scalar>("InitialConditions.SwFFInitial"); |
| 104 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | pnInitial_ = getParam<Scalar>("InitialConditions.pnInitial"); |
| 105 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | pnInjection_ = getParam<Scalar>("InitialConditions.pnInjection"); |
| 106 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | TInject_ = getParam<Scalar>("BoundaryConditions.TInject"); |
| 107 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | massFluxInjectedPhase_ = getParam<Scalar>("BoundaryConditions.massFluxInjectedPhase"); |
| 108 | |||
| 109 | // initialize the tables of the fluid system | ||
| 110 | 2 | FluidSystem::init(TInitial_ - 15.0, 453.15, nTemperature_, // T_min, T_max, n_T | |
| 111 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | 0.75*pnInitial_, 2.25*pnInitial_, nPressure_); // p_min, p_max, n_p |
| 112 | 1 | } | |
| 113 | |||
| 114 | void setGridVariables(std::shared_ptr<GridVariables> gridVariables) | ||
| 115 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | { gridVariables_ = gridVariables; } |
| 116 | |||
| 117 | const GridVariables& gridVariables() const | ||
| 118 | 3144960 | { return *gridVariables_; } | |
| 119 | |||
| 120 | /*! | ||
| 121 | * \name Problem parameters | ||
| 122 | */ | ||
| 123 | // \{ | ||
| 124 | |||
| 125 | /*! | ||
| 126 | * \brief Returns the problem name | ||
| 127 | * | ||
| 128 | * This is used as a prefix for files generated by the simulation. | ||
| 129 | */ | ||
| 130 | const std::string name() const | ||
| 131 |
2/6✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
1 | { return outputName_ ; } |
| 132 | |||
| 133 | // \} | ||
| 134 | |||
| 135 | /*! | ||
| 136 | * \name Boundary conditions | ||
| 137 | */ | ||
| 138 | // \{ | ||
| 139 | |||
| 140 | /*! | ||
| 141 | * \brief Specifies which kind of boundary condition should be | ||
| 142 | * used for which equation on a given boundary segment. | ||
| 143 | * | ||
| 144 | * \param globalPos The global position | ||
| 145 | */ | ||
| 146 | 6536 | BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const | |
| 147 | { | ||
| 148 | 6536 | BoundaryTypes bcTypes; | |
| 149 | // Default: Neumann | ||
| 150 | 6536 | bcTypes.setAllNeumann(); | |
| 151 | |||
| 152 | // Put a dirichlet somewhere: we need this for convergence | ||
| 153 |
6/6✓ Branch 0 taken 1102 times.
✓ Branch 1 taken 1178 times.
✓ Branch 2 taken 1102 times.
✓ Branch 3 taken 1178 times.
✓ Branch 4 taken 1102 times.
✓ Branch 5 taken 1178 times.
|
6840 | if(onRightBoundary_(globalPos) && globalPos[1] > this->spatialParams().heightPM() + eps_) |
| 154 | { | ||
| 155 | bcTypes.setAllDirichlet(); | ||
| 156 | } | ||
| 157 | 6536 | return bcTypes; | |
| 158 | } | ||
| 159 | |||
| 160 | /*! | ||
| 161 | * \brief Evaluates the boundary conditions for a Dirichlet boundary segment. | ||
| 162 | * | ||
| 163 | * \param globalPos The global position | ||
| 164 | * | ||
| 165 | */ | ||
| 166 | PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const | ||
| 167 | { | ||
| 168 | 1102 | return initial_(globalPos); | |
| 169 | } | ||
| 170 | |||
| 171 | /*! | ||
| 172 | * \brief Evaluates the boundary conditions for a Neumann boundary segment. | ||
| 173 | */ | ||
| 174 | 205276 | NumEqVector neumann(const Element& element, | |
| 175 | const FVElementGeometry& fvGeometry, | ||
| 176 | const ElementVolumeVariables& elemVolVars, | ||
| 177 | const ElementFluxVariablesCache& elemFluxVarsCache, | ||
| 178 | const SubControlVolumeFace& scvf) const | ||
| 179 | { | ||
| 180 | 205276 | NumEqVector values(0.0); | |
| 181 | 615828 | const auto& globalPos = fvGeometry.scv(scvf.insideScvIdx()).dofPosition(); | |
| 182 | 205276 | const Scalar massFluxInjectedPhase = massFluxInjectedPhase_ ; | |
| 183 | |||
| 184 | ParameterCache dummyCache; | ||
| 185 | 205276 | FluidState fluidState; | |
| 186 | |||
| 187 |
2/2✓ Branch 0 taken 410552 times.
✓ Branch 1 taken 205276 times.
|
615828 | for (int phaseIdx=0; phaseIdx<numPhases; phaseIdx++) |
| 188 | { | ||
| 189 | 821104 | fluidState.setPressure(phaseIdx, pnInitial_); | |
| 190 | } | ||
| 191 | |||
| 192 | 205276 | fluidState.setTemperature(gasPhaseIdx, TInject_ ); | |
| 193 | 205276 | fluidState.setTemperature(liquidPhaseIdx, TInitial_ ); // this value is a good one, TInject does not work | |
| 194 | |||
| 195 | // This solves the system of equations defining x=x(p,T) | ||
| 196 | 205276 | ConstraintSolver::solve(fluidState, | |
| 197 | dummyCache) ; | ||
| 198 | |||
| 199 | // Now let's make the air phase less than fully saturated with water | ||
| 200 | 410552 | fluidState.setMoleFraction(gasPhaseIdx, wCompIdx, fluidState.moleFraction(gasPhaseIdx, wCompIdx)*percentOfEquil_ ) ; | |
| 201 | 410552 | fluidState.setMoleFraction(gasPhaseIdx, nCompIdx, 1.-fluidState.moleFraction(gasPhaseIdx, wCompIdx) ) ; | |
| 202 | |||
| 203 | // compute density of injection phase | ||
| 204 | 205276 | const Scalar density = FluidSystem::density(fluidState, | |
| 205 | dummyCache, | ||
| 206 | gasPhaseIdx); | ||
| 207 | 205276 | fluidState.setDensity(gasPhaseIdx, density); | |
| 208 | 205276 | const Scalar molarDensity = FluidSystem::molarDensity(fluidState, | |
| 209 | dummyCache, | ||
| 210 | gasPhaseIdx); | ||
| 211 | 205276 | fluidState.setMolarDensity(gasPhaseIdx, molarDensity); | |
| 212 | |||
| 213 |
2/2✓ Branch 0 taken 410552 times.
✓ Branch 1 taken 205276 times.
|
615828 | for(int phaseIdx=0; phaseIdx<numPhases; phaseIdx++) |
| 214 | { | ||
| 215 | 410552 | const Scalar h = FluidSystem::enthalpy(fluidState, | |
| 216 | dummyCache, | ||
| 217 | phaseIdx); | ||
| 218 | 821104 | fluidState.setEnthalpy(phaseIdx, h); | |
| 219 | } | ||
| 220 | |||
| 221 |
2/2✓ Branch 0 taken 87172 times.
✓ Branch 1 taken 118104 times.
|
205276 | const Scalar molarFlux = massFluxInjectedPhase / fluidState.averageMolarMass(gasPhaseIdx); |
| 222 | |||
| 223 | // actually setting the fluxes | ||
| 224 |
2/4✓ Branch 0 taken 87172 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 87172 times.
✗ Branch 3 not taken.
|
174344 | if (onLeftBoundary_(globalPos) && this->spatialParams().inFF_(globalPos)) |
| 225 | { | ||
| 226 | 84360 | values[conti00EqIdx + gasPhaseIdx * numComponents + wCompIdx] | |
| 227 | 84360 | = -molarFlux * fluidState.moleFraction(gasPhaseIdx, wCompIdx); | |
| 228 | 84360 | values[conti00EqIdx + gasPhaseIdx * numComponents + nCompIdx] | |
| 229 | 84360 | = -molarFlux * fluidState.moleFraction(gasPhaseIdx, nCompIdx); | |
| 230 | // energy equations are specified mass specifically | ||
| 231 | 42180 | values[energyEq0Idx + gasPhaseIdx] = - massFluxInjectedPhase | |
| 232 | 84360 | * fluidState.enthalpy(gasPhaseIdx) ; | |
| 233 | } | ||
| 234 | 205276 | return values; | |
| 235 | } | ||
| 236 | |||
| 237 | /*! | ||
| 238 | * \name Volume terms | ||
| 239 | */ | ||
| 240 | // \{ | ||
| 241 | |||
| 242 | /*! | ||
| 243 | * \brief Evaluates the initial value for a control volume. | ||
| 244 | * | ||
| 245 | * \param globalPos The global position | ||
| 246 | */ | ||
| 247 | PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const | ||
| 248 | { | ||
| 249 | 465 | return initial_(globalPos); | |
| 250 | } | ||
| 251 | |||
| 252 | /*! | ||
| 253 | * \brief Evaluates the source term for all balance equations within a given | ||
| 254 | * sub-controlvolume. | ||
| 255 | * | ||
| 256 | * \param element The finite element | ||
| 257 | * \param fvGeometry The finite volume geometry of the element | ||
| 258 | * \param elemVolVars The volume variables of the element | ||
| 259 | * \param scv The sub-control volume | ||
| 260 | * | ||
| 261 | * Positive values mean that mass is created, negative ones mean that it vanishes. | ||
| 262 | */ | ||
| 263 | ✗ | NumEqVector source(const Element &element, | |
| 264 | const FVElementGeometry& fvGeometry, | ||
| 265 | const ElementVolumeVariables& elemVolVars, | ||
| 266 | const SubControlVolume &scv) const | ||
| 267 | { | ||
| 268 | 2362080 | NumEqVector values(0.0); | |
| 269 | ✗ | return values; | |
| 270 | |||
| 271 | } | ||
| 272 | |||
| 273 | private: | ||
| 274 | // the internal method for the initial condition | ||
| 275 | 1567 | PrimaryVariables initial_(const GlobalPosition &globalPos) const | |
| 276 | { | ||
| 277 | 1567 | PrimaryVariables priVars(0.0); | |
| 278 | 1567 | const Scalar T = TInitial_; | |
| 279 | Scalar S[numPhases]; | ||
| 280 | |||
| 281 |
2/4✓ Branch 0 taken 1567 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1567 times.
✗ Branch 3 not taken.
|
3134 | if (this->spatialParams().inPM_(globalPos)){ |
| 282 | 240 | S[liquidPhaseIdx] = SwPMInitial_; | |
| 283 | 240 | S[gasPhaseIdx] = 1. - S[liquidPhaseIdx] ; | |
| 284 | } | ||
| 285 |
2/4✓ Branch 0 taken 1327 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1327 times.
✗ Branch 3 not taken.
|
2654 | else if (this->spatialParams().inFF_(globalPos)){ |
| 286 | 1327 | S[liquidPhaseIdx] = SwFFInitial_; | |
| 287 | 1327 | S[gasPhaseIdx] = 1. - S[liquidPhaseIdx] ; | |
| 288 | } | ||
| 289 | else | ||
| 290 | ✗ | DUNE_THROW(Dune::InvalidStateException, | |
| 291 | "You should not be here: x=" << globalPos[0] << " y= "<< globalPos[dimWorld-1]); | ||
| 292 | |||
| 293 |
2/2✓ Branch 0 taken 1567 times.
✓ Branch 1 taken 1567 times.
|
3134 | for (int i = 0; i < numPhases - 1; ++i) |
| 294 | 3134 | priVars[s0Idx + i] = S[i]; | |
| 295 | |||
| 296 | // capillary pressure Params | ||
| 297 | 1567 | FluidState equilibriumFluidState; | |
| 298 | |||
| 299 | //set saturation to initial values, this needs to be done in order for the fluidState to tell me pc | ||
| 300 |
2/2✓ Branch 0 taken 3134 times.
✓ Branch 1 taken 1567 times.
|
4701 | for (int phaseIdx = 0; phaseIdx < numPhases ; ++phaseIdx) |
| 301 | { | ||
| 302 | 3134 | equilibriumFluidState.setSaturation(phaseIdx, S[phaseIdx]); | |
| 303 | 6268 | equilibriumFluidState.setTemperature(phaseIdx, TInitial_ ); | |
| 304 | } | ||
| 305 | |||
| 306 | // obtain pc according to saturation | ||
| 307 | 1567 | const int wPhaseIdx = this->spatialParams().template wettingPhaseAtPos<FluidSystem>(globalPos); | |
| 308 | 3134 | const auto& fm = this->spatialParams().fluidMatrixInteractionAtPos(globalPos); | |
| 309 | 1567 | const auto capPress = fm.capillaryPressures(equilibriumFluidState, wPhaseIdx); | |
| 310 | |||
| 311 | Scalar p[numPhases]; | ||
| 312 |
2/4✓ Branch 0 taken 1567 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1567 times.
✗ Branch 3 not taken.
|
3134 | if (this->spatialParams().inPM_(globalPos)){ |
| 313 | // Use homogeneous pressure in the domain and let the newton find the pressure distribution | ||
| 314 | using std::abs; | ||
| 315 | 480 | p[liquidPhaseIdx] = pnInitial_ - abs(capPress[liquidPhaseIdx]); | |
| 316 | 720 | p[gasPhaseIdx] = p[liquidPhaseIdx] + abs(capPress[liquidPhaseIdx]); | |
| 317 | } | ||
| 318 |
2/4✓ Branch 0 taken 1327 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1327 times.
✗ Branch 3 not taken.
|
2654 | else if (this->spatialParams().inFF_(globalPos)){ |
| 319 | 1327 | p[gasPhaseIdx] = pnInitial_ ; | |
| 320 | 1327 | p[liquidPhaseIdx] = pnInitial_ ; | |
| 321 | } | ||
| 322 | else | ||
| 323 | ✗ | DUNE_THROW(Dune::InvalidStateException, "You should not be here: x=" << globalPos[0] << " y= "<< globalPos[dimWorld-1]); | |
| 324 | |||
| 325 | if(pressureFormulation == mostWettingFirst){ | ||
| 326 | // This means that the pressures are sorted from the most wetting to the least wetting-1 in the primary variables vector. | ||
| 327 | // For two phases this means that there is one pressure as primary variable: pw | ||
| 328 | priVars[p0Idx] = p[liquidPhaseIdx]; | ||
| 329 | } | ||
| 330 | else if(pressureFormulation == leastWettingFirst){ | ||
| 331 | // This means that the pressures are sorted from the least wetting to the most wetting-1 in the primary variables vector. | ||
| 332 | // For two phases this means that there is one pressure as primary variable: pn | ||
| 333 | 1567 | priVars[p0Idx] = p[gasPhaseIdx]; | |
| 334 | } | ||
| 335 | else DUNE_THROW(Dune::InvalidStateException, "EvaporationAtmosphereProblem does not support the chosen pressure formulation."); | ||
| 336 | |||
| 337 |
2/2✓ Branch 0 taken 4701 times.
✓ Branch 1 taken 1567 times.
|
6268 | for (int energyEqIdx=0; energyEqIdx< numEnergyEqFluid+numEnergyEqSolid; ++energyEqIdx) |
| 338 | 9402 | priVars[energyEq0Idx + energyEqIdx] = T; | |
| 339 | |||
| 340 |
2/2✓ Branch 0 taken 3134 times.
✓ Branch 1 taken 1567 times.
|
4701 | for (int phaseIdx=0; phaseIdx<numPhases; phaseIdx++) |
| 341 | 6268 | equilibriumFluidState.setPressure(phaseIdx, p[phaseIdx]); | |
| 342 | |||
| 343 | // This solves the system of equations defining x=x(p,T) | ||
| 344 | ParameterCache dummyCache; | ||
| 345 | 1567 | ConstraintSolver::solve(equilibriumFluidState, | |
| 346 | dummyCache) ; | ||
| 347 | |||
| 348 | 1567 | FluidState dryFluidState(equilibriumFluidState); | |
| 349 | // Now let's make the air phase less than fully saturated with vapor | ||
| 350 | 3134 | dryFluidState.setMoleFraction(gasPhaseIdx, wCompIdx, dryFluidState.moleFraction(gasPhaseIdx, wCompIdx) * percentOfEquil_ ) ; | |
| 351 | 1567 | dryFluidState.setMoleFraction(gasPhaseIdx, nCompIdx, 1.0-dryFluidState.moleFraction(gasPhaseIdx, wCompIdx) ) ; | |
| 352 | |||
| 353 | /* Difference between kinetic and MPNC: | ||
| 354 | * number of component related primVar and how they are calculated (mole fraction, fugacities, resp.) | ||
| 355 | */ | ||
| 356 | if(enableChemicalNonEquilibrium){ | ||
| 357 |
2/2✓ Branch 1 taken 3134 times.
✓ Branch 2 taken 1567 times.
|
4701 | for(int phaseIdx=0; phaseIdx < numPhases; ++ phaseIdx) |
| 358 | { | ||
| 359 |
2/2✓ Branch 0 taken 6268 times.
✓ Branch 1 taken 3134 times.
|
9402 | for(int compIdx=0; compIdx <numComponents; ++compIdx){ |
| 360 | 6268 | int offset = compIdx + phaseIdx * numComponents ; | |
| 361 | |||
| 362 |
2/4✓ Branch 0 taken 6268 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 6268 times.
✗ Branch 3 not taken.
|
12536 | if (this->spatialParams().inPM_(globalPos)){ |
| 363 | 2880 | priVars[conti00EqIdx + offset] = equilibriumFluidState.moleFraction(phaseIdx,compIdx) ; | |
| 364 | } | ||
| 365 |
2/4✓ Branch 0 taken 5308 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 5308 times.
✗ Branch 3 not taken.
|
10616 | else if (this->spatialParams().inFF_(globalPos)){ |
| 366 | 15924 | priVars[conti00EqIdx + offset] = dryFluidState.moleFraction(phaseIdx,compIdx) ; | |
| 367 | } | ||
| 368 | else | ||
| 369 | ✗ | DUNE_THROW(Dune::InvalidStateException, "You should not be here: x=" << globalPos[0] << " y= "<< globalPos[dimWorld-1]); | |
| 370 | } | ||
| 371 | } | ||
| 372 | } | ||
| 373 | else | ||
| 374 | { | ||
| 375 | // in the case I am using the "standard" mpnc model, the variables to be set are the "fugacities" | ||
| 376 | const Scalar fugH2O = FluidSystem::H2O::vaporPressure(T) ; | ||
| 377 | const Scalar fugN2 = p[gasPhaseIdx] - fugH2O ; | ||
| 378 | |||
| 379 | priVars[conti00EqIdx + FluidSystem::N2Idx] = fugN2 ; | ||
| 380 | priVars[conti00EqIdx + FluidSystem::H2OIdx] = fugH2O ; | ||
| 381 | |||
| 382 | Scalar xl[numComponents]; | ||
| 383 | Scalar beta[numComponents]; | ||
| 384 | |||
| 385 | const Scalar Henry = BinaryCoeff::H2O_N2::henry(TInitial_); | ||
| 386 | const Scalar satVapPressure = FluidSystem::H2O::vaporPressure(TInitial_); | ||
| 387 | xl[FluidSystem::H2OIdx] = x_[liquidPhaseIdx][wCompIdx]; | ||
| 388 | xl[FluidSystem::N2Idx] = x_[liquidPhaseIdx][nCompIdx]; | ||
| 389 | beta[FluidSystem::H2OIdx] = satVapPressure ; | ||
| 390 | beta[FluidSystem::N2Idx] = Henry ; | ||
| 391 | |||
| 392 | for (int i = 0; i < numComponents; ++i) | ||
| 393 | priVars[conti00EqIdx + i] = xl[i]*beta[i]; // this should be really fug0Idx but the compiler only knows one or the other | ||
| 394 | } | ||
| 395 | 1567 | return priVars; | |
| 396 | } | ||
| 397 | |||
| 398 | /*! | ||
| 399 | * \brief Returns whether the tested position is on the left boundary of the domain. | ||
| 400 | */ | ||
| 401 | bool onLeftBoundary_(const GlobalPosition & globalPos) const | ||
| 402 |
10/10✓ Branch 0 taken 87172 times.
✓ Branch 1 taken 118104 times.
✓ Branch 2 taken 87172 times.
✓ Branch 3 taken 118104 times.
✓ Branch 4 taken 87172 times.
✓ Branch 5 taken 118104 times.
✓ Branch 6 taken 87172 times.
✓ Branch 7 taken 118104 times.
✓ Branch 8 taken 87172 times.
✓ Branch 9 taken 118104 times.
|
1026380 | { return globalPos[0] < this->gridGeometry().bBoxMin()[0] + eps_; } |
| 403 | |||
| 404 | /*! | ||
| 405 | * \brief Returns whether the tested position is on the right boundary of the domain. | ||
| 406 | */ | ||
| 407 | bool onRightBoundary_(const GlobalPosition & globalPos) const | ||
| 408 |
10/10✓ Branch 0 taken 2280 times.
✓ Branch 1 taken 4256 times.
✓ Branch 2 taken 2280 times.
✓ Branch 3 taken 4256 times.
✓ Branch 4 taken 2280 times.
✓ Branch 5 taken 4256 times.
✓ Branch 6 taken 2280 times.
✓ Branch 7 taken 4256 times.
✓ Branch 8 taken 2280 times.
✓ Branch 9 taken 4256 times.
|
32680 | { return globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_; } |
| 409 | |||
| 410 | /*! | ||
| 411 | * \brief Returns whether the tested position is on the lower boundary of the domain. | ||
| 412 | */ | ||
| 413 | bool onLowerBoundary_(const GlobalPosition & globalPos) const | ||
| 414 | { return globalPos[dimWorld-1] < this->gridGeometry().bBoxMin()[dimWorld-1] + eps_; } | ||
| 415 | |||
| 416 | private: | ||
| 417 | static constexpr Scalar eps_ = 1e-6; | ||
| 418 | Scalar percentOfEquil_ ; | ||
| 419 | int nTemperature_; | ||
| 420 | int nPressure_; | ||
| 421 | std::string outputName_; | ||
| 422 | Scalar heatIntoSolid_; | ||
| 423 | Scalar TInitial_ ; | ||
| 424 | Scalar SwPMInitial_ ; | ||
| 425 | Scalar SwFFInitial_ ; | ||
| 426 | Scalar SnInitial_; | ||
| 427 | Scalar pnInitial_; | ||
| 428 | Scalar pnInjection_; | ||
| 429 | Dune::ParameterTree inputParameters_; | ||
| 430 | Scalar x_[numPhases][numComponents] ; | ||
| 431 | |||
| 432 | Scalar TInject_; | ||
| 433 | |||
| 434 | Scalar massFluxInjectedPhase_ ; | ||
| 435 | |||
| 436 | std::shared_ptr<GridVariables> gridVariables_; | ||
| 437 | |||
| 438 | public: | ||
| 439 | |||
| 440 | Dune::ParameterTree getInputParameters() const | ||
| 441 | { return inputParameters_; } | ||
| 442 | }; | ||
| 443 | |||
| 444 | } // end namespace | ||
| 445 | |||
| 446 | #endif | ||
| 447 |