| 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 TwoPNCMinTests | ||
| 10 | * \brief Problem where brine is evaporating at the top boundary. The system is closed at the remaining boundaries. | ||
| 11 | */ | ||
| 12 | #ifndef DUMUX_SALINIZATION_PROBLEM_HH | ||
| 13 | #define DUMUX_SALINIZATION_PROBLEM_HH | ||
| 14 | |||
| 15 | #include <dumux/common/boundarytypes.hh> | ||
| 16 | #include <dumux/common/properties.hh> | ||
| 17 | #include <dumux/common/parameters.hh> | ||
| 18 | #include <dumux/common/numeqvector.hh> | ||
| 19 | |||
| 20 | #include <dumux/porousmediumflow/problem.hh> | ||
| 21 | #include <dumux/io/gnuplotinterface.hh> | ||
| 22 | |||
| 23 | namespace Dumux { | ||
| 24 | |||
| 25 | /*! | ||
| 26 | * \ingroup TwoPNCMinTests | ||
| 27 | * \brief Problem where water is evaporating at the top of a porous media filled container saturated with brine and air, which causes precipitation at the top. | ||
| 28 | * | ||
| 29 | */ | ||
| 30 | template <class TypeTag> | ||
| 31 | class SalinizationProblem : public PorousMediumFlowProblem<TypeTag> | ||
| 32 | { | ||
| 33 | using ParentType = PorousMediumFlowProblem<TypeTag>; | ||
| 34 | using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView; | ||
| 35 | using Scalar = GetPropType<TypeTag, Properties::Scalar>; | ||
| 36 | using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>; | ||
| 37 | using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>; | ||
| 38 | using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices; | ||
| 39 | using SolidSystem = GetPropType<TypeTag, Properties::SolidSystem>; | ||
| 40 | |||
| 41 | enum | ||
| 42 | { | ||
| 43 | // primary variable indices | ||
| 44 | pressureIdx = Indices::pressureIdx, | ||
| 45 | switchIdx = Indices::switchIdx, | ||
| 46 | |||
| 47 | // component indices | ||
| 48 | // TODO: using xwNaClIdx as privaridx works here, but | ||
| 49 | // looks like magic. Can this be done differently?? | ||
| 50 | xwNaClIdx = FluidSystem::NaClIdx, | ||
| 51 | precipNaClIdx = FluidSystem::numComponents, | ||
| 52 | |||
| 53 | // Indices of the components | ||
| 54 | H2OIdx = FluidSystem::H2OIdx, | ||
| 55 | NaClIdx = FluidSystem::NaClIdx, | ||
| 56 | AirIdx = FluidSystem::AirIdx, | ||
| 57 | |||
| 58 | // Indices of the phases | ||
| 59 | liquidPhaseIdx = FluidSystem::liquidPhaseIdx, | ||
| 60 | gasPhaseIdx = FluidSystem::gasPhaseIdx, | ||
| 61 | |||
| 62 | // index of the solid phase | ||
| 63 | sPhaseIdx = SolidSystem::comp0Idx, | ||
| 64 | |||
| 65 | |||
| 66 | // Index of the primary component of G and L phase | ||
| 67 | conti0EqIdx = Indices::conti0EqIdx, //water component | ||
| 68 | conti1EqIdx = Indices::conti0EqIdx + 1, //air component | ||
| 69 | precipNaClEqIdx = Indices::conti0EqIdx + FluidSystem::numComponents, | ||
| 70 | energyEqIdx = Indices::energyEqIdx, | ||
| 71 | |||
| 72 | // Phase State | ||
| 73 | bothPhases = Indices::bothPhases, | ||
| 74 | |||
| 75 | // Grid and world dimension | ||
| 76 | dim = GridView::dimension, | ||
| 77 | dimWorld = GridView::dimensionworld, | ||
| 78 | }; | ||
| 79 | |||
| 80 | using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>; | ||
| 81 | using NumEqVector = Dumux::NumEqVector<PrimaryVariables>; | ||
| 82 | using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>; | ||
| 83 | using GridVariables = GetPropType<TypeTag, Properties::GridVariables>; | ||
| 84 | using ElementFluxVariablesCache = typename GridVariables::GridFluxVariablesCache::LocalView; | ||
| 85 | using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView; | ||
| 86 | using Element = typename GridView::template Codim<0>::Entity; | ||
| 87 | using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>; | ||
| 88 | using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>; | ||
| 89 | using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView; | ||
| 90 | using SubControlVolume = typename FVElementGeometry::SubControlVolume; | ||
| 91 | using GlobalPosition = typename SubControlVolume::GlobalPosition; | ||
| 92 | using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; | ||
| 93 | using FluidState = GetPropType<TypeTag, Properties::FluidState>; | ||
| 94 | |||
| 95 | public: | ||
| 96 | 2 | SalinizationProblem(std::shared_ptr<const GridGeometry> gridGeometry) | |
| 97 |
13/40✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 2 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 2 times.
✓ Branch 15 taken 2 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 2 times.
✗ Branch 19 not taken.
✓ Branch 21 taken 2 times.
✗ Branch 22 not taken.
✓ Branch 24 taken 2 times.
✗ Branch 25 not taken.
✓ Branch 27 taken 2 times.
✗ Branch 28 not taken.
✓ Branch 30 taken 2 times.
✗ Branch 31 not taken.
✓ Branch 33 taken 2 times.
✗ Branch 34 not taken.
✓ Branch 36 taken 2 times.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✗ Branch 39 not taken.
✗ Branch 40 not taken.
✗ Branch 41 not taken.
✗ Branch 44 not taken.
✗ Branch 45 not taken.
✗ Branch 46 not taken.
✗ Branch 47 not taken.
✗ Branch 48 not taken.
✗ Branch 49 not taken.
✗ Branch 52 not taken.
✗ Branch 53 not taken.
✗ Branch 54 not taken.
✗ Branch 55 not taken.
|
6 | : ParentType(gridGeometry) |
| 98 | { | ||
| 99 | //Fluidsystem | ||
| 100 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | nTemperature_ = getParam<int>("FluidSystem.NTemperature"); |
| 101 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | nPressure_ = getParam<int>("FluidSystem.NPressure"); |
| 102 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | pressureLow_ = getParam<Scalar>("FluidSystem.PressureLow"); |
| 103 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | pressureHigh_ = getParam<Scalar>("FluidSystem.PressureHigh"); |
| 104 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | temperatureLow_ = getParam<Scalar>("FluidSystem.TemperatureLow"); |
| 105 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | temperatureHigh_ = getParam<Scalar>("FluidSystem.TemperatureHigh"); |
| 106 |
2/4✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
2 | name_ = getParam<std::string>("Problem.Name"); |
| 107 | |||
| 108 | //problem | ||
| 109 |
2/4✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
|
2 | name_ = getParam<std::string>("Problem.Name"); |
| 110 | |||
| 111 | //initial conditions | ||
| 112 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | initPressure_ = getParam<Scalar>("Problem.InitialPressure"); |
| 113 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | initGasSaturation_ = getParam<Scalar>("Problem.InitialGasSaturation"); |
| 114 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | initSalinity_ = getParam<Scalar>("Problem.InitialSalinity"); |
| 115 | |||
| 116 | 2 | unsigned int codim = GetPropType<TypeTag, Properties::GridGeometry>::discMethod == DiscretizationMethods::box ? dim : 0; | |
| 117 | |||
| 118 |
4/8✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2 times.
✗ Branch 11 not taken.
|
6 | permeability_.resize(gridGeometry->gridView().size(codim)); |
| 119 | 4 | FluidSystem::init(/*Tmin=*/temperatureLow_, | |
| 120 | /*Tmax=*/temperatureHigh_, | ||
| 121 | 2 | /*nT=*/nTemperature_, | |
| 122 | /*pmin=*/pressureLow_, | ||
| 123 | /*pmax=*/pressureHigh_, | ||
| 124 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | /*np=*/nPressure_); |
| 125 | 2 | } | |
| 126 | |||
| 127 | /*! | ||
| 128 | * \brief The current time. | ||
| 129 | */ | ||
| 130 | ✗ | void setTime( Scalar time ) | |
| 131 | { | ||
| 132 |
1/2✓ Branch 1 taken 84 times.
✗ Branch 2 not taken.
|
84 | time_ = time; |
| 133 | ✗ | } | |
| 134 | |||
| 135 | /*! | ||
| 136 | * \brief The time step size. | ||
| 137 | * | ||
| 138 | * This is used to calculate the source term. | ||
| 139 | */ | ||
| 140 | ✗ | void setTimeStepSize( Scalar timeStepSize ) | |
| 141 | { | ||
| 142 |
1/2✓ Branch 1 taken 84 times.
✗ Branch 2 not taken.
|
84 | timeStepSize_ = timeStepSize; |
| 143 | ✗ | } | |
| 144 | |||
| 145 | |||
| 146 | /*! | ||
| 147 | * \brief The problem name. | ||
| 148 | * | ||
| 149 | * This is used as a prefix for files generated by the simulation. | ||
| 150 | */ | ||
| 151 | const std::string& name() const | ||
| 152 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | { return name_; } |
| 153 | |||
| 154 | /*! | ||
| 155 | * \brief Specifies which kind of boundary condition should be | ||
| 156 | * used for which equation on a given boundary segment. | ||
| 157 | */ | ||
| 158 | ✗ | BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const | |
| 159 | { | ||
| 160 | 27622 | BoundaryTypes bcTypes; | |
| 161 | |||
| 162 | // default to Neumann | ||
| 163 |
5/8✓ Branch 0 taken 4200 times.
✓ Branch 1 taken 1204 times.
✓ Branch 2 taken 3654 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 4032 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 14532 times.
|
27622 | bcTypes.setAllNeumann(); |
| 164 | |||
| 165 | ✗ | return bcTypes; | |
| 166 | } | ||
| 167 | |||
| 168 | /*! | ||
| 169 | * \brief Evaluates the boundary conditions for a Dirichlet boundary segment. | ||
| 170 | */ | ||
| 171 | ✗ | PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const | |
| 172 | { | ||
| 173 | ✗ | PrimaryVariables priVars(0.0); | |
| 174 | ✗ | priVars.setState(bothPhases); | |
| 175 | |||
| 176 | ✗ | return priVars; | |
| 177 | } | ||
| 178 | |||
| 179 | |||
| 180 | /*! | ||
| 181 | * \brief Evaluate the boundary conditions for a neumann | ||
| 182 | * boundary segment. | ||
| 183 | * | ||
| 184 | * This is the method for the case where the Neumann condition is | ||
| 185 | * potentially solution dependent | ||
| 186 | * | ||
| 187 | * \param element The finite element | ||
| 188 | * \param fvGeometry The finite-volume geometry | ||
| 189 | * \param elemVolVars All volume variables for the element | ||
| 190 | * \param elemFluxVarsCache Flux variables caches for all faces in stencil | ||
| 191 | * \param scvf The sub control volume face | ||
| 192 | * | ||
| 193 | * Negative values mean influx. | ||
| 194 | * E.g. for the mass balance that would be the mass flux in \f$ [ kg / (m^2 \cdot s)] \f$. | ||
| 195 | */ | ||
| 196 | 118636 | NumEqVector neumann(const Element& element, | |
| 197 | const FVElementGeometry& fvGeometry, | ||
| 198 | const ElementVolumeVariables& elemVolVars, | ||
| 199 | const ElementFluxVariablesCache& elemFluxVarsCache, | ||
| 200 | const SubControlVolumeFace& scvf) const | ||
| 201 | { | ||
| 202 | 118636 | PrimaryVariables values(0.0); | |
| 203 | |||
| 204 |
2/2✓ Branch 0 taken 1 times.
✓ Branch 1 taken 102899 times.
|
118636 | const auto& globalPos = scvf.ipGlobal(); |
| 205 |
4/4✓ Branch 0 taken 1 times.
✓ Branch 1 taken 102899 times.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 102899 times.
|
237272 | const auto& volVars = elemVolVars[scvf.insideScvIdx()]; |
| 206 |
6/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 118634 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 118634 times.
✓ Branch 4 taken 2 times.
✓ Branch 5 taken 118634 times.
|
355908 | const Scalar hmax = this->gridGeometry().bBoxMax()[1]; |
| 207 | |||
| 208 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 118634 times.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
|
118636 | static const Scalar temperatureRef = getParam<Scalar>("FreeFlow.RefTemperature"); |
| 209 | |||
| 210 |
4/4✓ Branch 0 taken 16948 times.
✓ Branch 1 taken 101688 times.
✓ Branch 2 taken 16948 times.
✓ Branch 3 taken 101688 times.
|
237272 | if (globalPos[1] > hmax - eps_) |
| 211 | { | ||
| 212 | // get free-flow properties: | ||
| 213 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 16946 times.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
|
16948 | static const Scalar moleFracRefH2O = getParam<Scalar>("FreeFlow.RefMoleFracH2O"); |
| 214 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 16946 times.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
|
16948 | static const Scalar massFracRefH2O = refMoleToRefMassFrac_(moleFracRefH2O); |
| 215 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 16946 times.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
|
16948 | static const Scalar boundaryLayerThickness = getParam<Scalar>("FreeFlow.BoundaryLayerThickness"); |
| 216 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 16946 times.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
|
16948 | static const Scalar massTransferCoefficient = getParam<Scalar>("FreeFlow.MassTransferCoefficient"); |
| 217 | |||
| 218 | // get porous medium values: | ||
| 219 | 16948 | const Scalar massFracH2OInside = volVars.massFraction(gasPhaseIdx, H2OIdx); | |
| 220 | |||
| 221 | |||
| 222 | // calculate fluxes | ||
| 223 | // liquid phase | ||
| 224 | 16948 | Scalar evaporationRateMole = 0; | |
| 225 |
1/2✓ Branch 0 taken 16948 times.
✗ Branch 1 not taken.
|
16948 | if (massFracH2OInside - massFracRefH2O > 0) |
| 226 | { | ||
| 227 | 16948 | evaporationRateMole = massTransferCoefficient | |
| 228 | 16948 | * volVars.diffusionCoefficient(gasPhaseIdx, AirIdx, H2OIdx) | |
| 229 | 16948 | * (massFracH2OInside - massFracRefH2O) | |
| 230 | 16948 | / boundaryLayerThickness | |
| 231 | 33896 | * volVars.molarDensity(gasPhaseIdx); | |
| 232 | } | ||
| 233 | else | ||
| 234 | { | ||
| 235 | ✗ | evaporationRateMole = massTransferCoefficient | |
| 236 | ✗ | * volVars.diffusionCoefficient(gasPhaseIdx, AirIdx, H2OIdx) | |
| 237 | ✗ | * (massFracH2OInside - massFracRefH2O) | |
| 238 | ✗ | / boundaryLayerThickness | |
| 239 | * 1.2; | ||
| 240 | |||
| 241 | } | ||
| 242 | |||
| 243 |
1/2✓ Branch 0 taken 16948 times.
✗ Branch 1 not taken.
|
16948 | values[conti0EqIdx] = evaporationRateMole; |
| 244 | |||
| 245 | // gas phase | ||
| 246 | // gas flows in | ||
| 247 |
2/4✓ Branch 0 taken 16948 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 16948 times.
✗ Branch 3 not taken.
|
33896 | if (volVars.pressure(gasPhaseIdx) - 1e5 > 0) { |
| 248 |
1/2✓ Branch 0 taken 2248 times.
✗ Branch 1 not taken.
|
16948 | values[conti1EqIdx] = (volVars.pressure(gasPhaseIdx) - 1e5) |
| 249 |
2/4✓ Branch 0 taken 2248 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2248 times.
✗ Branch 3 not taken.
|
101688 | /(globalPos - fvGeometry.scv(scvf.insideScvIdx()).center()).two_norm() |
| 250 | 16948 | *volVars.mobility(gasPhaseIdx) | |
| 251 | 16948 | *volVars.permeability() | |
| 252 | 16948 | *volVars.molarDensity(gasPhaseIdx) | |
| 253 | 33896 | *volVars.moleFraction(gasPhaseIdx, AirIdx); | |
| 254 | } | ||
| 255 | //gas flows out | ||
| 256 | else { | ||
| 257 | ✗ | values[conti1EqIdx] = (volVars.pressure(gasPhaseIdx) - 1e5) | |
| 258 | ✗ | /(globalPos - fvGeometry.scv(scvf.insideScvIdx()).center()).two_norm() | |
| 259 | ✗ | *volVars.mobility(gasPhaseIdx) | |
| 260 | ✗ | *volVars.permeability() | |
| 261 | ✗ | *volVars.molarDensity(gasPhaseIdx) * (1-moleFracRefH2O); | |
| 262 | } | ||
| 263 | |||
| 264 | // energy fluxes | ||
| 265 | 33896 | values[energyEqIdx] = FluidSystem::componentEnthalpy(volVars.fluidState(), gasPhaseIdx, H2OIdx) * values[conti0EqIdx] * FluidSystem::molarMass(H2OIdx); | |
| 266 | |||
| 267 | 16948 | values[energyEqIdx] += FluidSystem::componentEnthalpy(volVars.fluidState(), gasPhaseIdx, AirIdx)* values[conti1EqIdx] * FluidSystem::molarMass(AirIdx); | |
| 268 | |||
| 269 | 63296 | values[energyEqIdx] += FluidSystem::thermalConductivity(elemVolVars[scvf.insideScvIdx()].fluidState(), gasPhaseIdx) * (volVars.temperature() - temperatureRef)/boundaryLayerThickness; | |
| 270 | } | ||
| 271 | 118636 | return values; | |
| 272 | } | ||
| 273 | |||
| 274 | /*! | ||
| 275 | * \brief Evaluates the initial value for a control volume. | ||
| 276 | * | ||
| 277 | * \param globalPos The global position | ||
| 278 | * | ||
| 279 | * For this method, the \a values parameter stores primary | ||
| 280 | * variables. | ||
| 281 | */ | ||
| 282 | 28 | PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const | |
| 283 | { | ||
| 284 | 28 | PrimaryVariables priVars(0.0); | |
| 285 | 28 | priVars.setState(bothPhases); | |
| 286 | 28 | Scalar density = 1000.00; //FluidSystem::density(, liquidPhaseIdx); | |
| 287 | |||
| 288 | 56 | priVars[pressureIdx] = initPressure_ - density*9.81*globalPos[dimWorld-1]; | |
| 289 | 28 | priVars[switchIdx] = initGasSaturation_; // Sg primary variable | |
| 290 | 28 | priVars[xwNaClIdx] = massToMoleFrac_(initSalinity_); // mole fraction | |
| 291 | 56 | priVars[precipNaClIdx] = 0.0; // [kg/m^3] | |
| 292 | 84 | priVars[energyEqIdx] = this->spatialParams().temperature(); // [K] | |
| 293 | |||
| 294 | 28 | return priVars; | |
| 295 | } | ||
| 296 | |||
| 297 | /*! | ||
| 298 | * \brief Evaluates the source term for all phases within a given | ||
| 299 | * sub-controlvolume. | ||
| 300 | * | ||
| 301 | * This is the method for the case where the source term is | ||
| 302 | * potentially solution dependent and requires some quantities that | ||
| 303 | * are specific to the fully-implicit method. | ||
| 304 | * | ||
| 305 | * \param element The finite element | ||
| 306 | * \param fvGeometry The finite-volume geometry | ||
| 307 | * \param elemVolVars All volume variables for the element | ||
| 308 | * \param scv The subcontrolvolume | ||
| 309 | * | ||
| 310 | * Positive values mean that the conserved quantity is created, negative ones mean that it vanishes. | ||
| 311 | * E.g. for the mass balance that would be a mass rate in \f$ [ kg / (m^3 \cdot s)] \f$. | ||
| 312 | */ | ||
| 313 | 157380 | NumEqVector source(const Element &element, | |
| 314 | const FVElementGeometry& fvGeometry, | ||
| 315 | const ElementVolumeVariables& elemVolVars, | ||
| 316 | const SubControlVolume &scv) const | ||
| 317 | { | ||
| 318 | 157380 | NumEqVector source(0.0); | |
| 319 | |||
| 320 | 157380 | const auto& volVars = elemVolVars[scv]; | |
| 321 | |||
| 322 | 157380 | Scalar moleFracNaCl_wPhase = volVars.moleFraction(liquidPhaseIdx, NaClIdx); | |
| 323 | 157380 | Scalar massFracNaCl_Max_wPhase = this->spatialParams().solubilityLimit(); | |
| 324 | 157380 | Scalar moleFracNaCl_Max_wPhase = massToMoleFrac_(massFracNaCl_Max_wPhase); | |
| 325 | 157380 | Scalar saltPorosity = this->spatialParams().minimalPorosity(element, scv); | |
| 326 | |||
| 327 | // precipitation of amount of salt which hexeeds the solubility limit | ||
| 328 | using std::abs; | ||
| 329 | 314760 | Scalar precipSalt = volVars.porosity() * volVars.molarDensity(liquidPhaseIdx) | |
| 330 |
2/2✓ Branch 0 taken 154716 times.
✓ Branch 1 taken 2664 times.
|
157380 | * volVars.saturation(liquidPhaseIdx) |
| 331 |
2/2✓ Branch 0 taken 154716 times.
✓ Branch 1 taken 2664 times.
|
157380 | * abs(moleFracNaCl_wPhase - moleFracNaCl_Max_wPhase) |
| 332 | 157380 | / timeStepSize_; | |
| 333 |
2/2✓ Branch 0 taken 154716 times.
✓ Branch 1 taken 2664 times.
|
157380 | if (moleFracNaCl_wPhase < moleFracNaCl_Max_wPhase) |
| 334 | 154716 | precipSalt *= -1; | |
| 335 | |||
| 336 | // make sure we don't dissolve more salt than previously precipitated | ||
| 337 |
2/2✓ Branch 2 taken 154716 times.
✓ Branch 3 taken 2664 times.
|
314760 | if (precipSalt*timeStepSize_ + volVars.solidVolumeFraction(sPhaseIdx)* volVars.solidComponentMolarDensity(sPhaseIdx)< 0) |
| 338 | 309432 | precipSalt = -volVars.solidVolumeFraction(sPhaseIdx)* volVars.solidComponentMolarDensity(sPhaseIdx)/timeStepSize_; | |
| 339 | |||
| 340 | // make sure there is still pore space available for precipitation | ||
| 341 |
3/8✗ Branch 0 not taken.
✓ Branch 1 taken 157380 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 157380 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 157380 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
472140 | if (volVars.solidVolumeFraction(sPhaseIdx) >= this->spatialParams().referencePorosity(element, scv) - saltPorosity && precipSalt > 0) |
| 342 | precipSalt = 0; | ||
| 343 | |||
| 344 | 314760 | source[conti0EqIdx + NaClIdx] += -precipSalt; | |
| 345 | 314760 | source[precipNaClEqIdx] += precipSalt; | |
| 346 | 157380 | return source; | |
| 347 | } | ||
| 348 | |||
| 349 | /*! | ||
| 350 | * \brief Adds additional VTK output data to the VTKWriter. | ||
| 351 | * | ||
| 352 | * Function is called by the output module on every write. | ||
| 353 | */ | ||
| 354 | |||
| 355 | const std::vector<Scalar>& getPermeability() | ||
| 356 | { | ||
| 357 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | return permeability_; |
| 358 | } | ||
| 359 | |||
| 360 | 86 | void updateVtkOutput(const SolutionVector& curSol) | |
| 361 | { | ||
| 362 |
2/4✓ Branch 1 taken 86 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 86 times.
✗ Branch 5 not taken.
|
172 | auto fvGeometry = localView(this->gridGeometry()); |
| 363 |
5/10✓ Branch 1 taken 86 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 86 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 86 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 946 times.
✗ Branch 10 not taken.
✓ Branch 13 taken 430 times.
✗ Branch 14 not taken.
|
2064 | for (const auto& element : elements(this->gridGeometry().gridView())) |
| 364 | { | ||
| 365 |
2/4✓ Branch 1 taken 430 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 430 times.
✗ Branch 5 not taken.
|
1720 | const auto elemSol = elementSolution(element, curSol, this->gridGeometry()); |
| 366 |
2/4✓ Branch 1 taken 860 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 430 times.
✗ Branch 5 not taken.
|
1290 | fvGeometry.bindElement(element); |
| 367 |
4/4✓ Branch 0 taken 2150 times.
✓ Branch 1 taken 860 times.
✓ Branch 2 taken 1720 times.
✓ Branch 3 taken 430 times.
|
5160 | for (auto&& scv : scvs(fvGeometry)) |
| 368 | { | ||
| 369 | 2150 | VolumeVariables volVars; | |
| 370 |
1/2✓ Branch 1 taken 2150 times.
✗ Branch 2 not taken.
|
2150 | volVars.update(elemSol, *this, element, scv); |
| 371 | 2150 | const auto dofIdxGlobal = scv.dofIndex(); | |
| 372 | 4300 | permeability_[dofIdxGlobal] = volVars.permeability(); | |
| 373 | } | ||
| 374 | } | ||
| 375 | 86 | } | |
| 376 | |||
| 377 | private: | ||
| 378 | |||
| 379 | /*! | ||
| 380 | * \brief Returns the molality of NaCl (mol NaCl / kg water) for a given mole fraction. | ||
| 381 | * | ||
| 382 | * \param massFracNaClInWater The mass fraction of NaCl in liquid phase water [kg NaCl / kg solution] | ||
| 383 | */ | ||
| 384 | 157408 | static Scalar massToMoleFrac_(Scalar massFracNaClInWater) | |
| 385 | { | ||
| 386 | 157408 | const Scalar molarMassWater = FluidSystem::molarMass(H2OIdx); /* molecular weight of water [kg/mol] */ | |
| 387 | 157408 | const Scalar molarMassNaCl = FluidSystem::molarMass(NaClIdx); /* molecular weight of NaCl [kg/mol] */ | |
| 388 | |||
| 389 | 157408 | const auto moleFracNaClInWater = -molarMassWater * massFracNaClInWater / ((molarMassNaCl - molarMassWater) * massFracNaClInWater - molarMassNaCl); | |
| 390 | 157408 | return moleFracNaClInWater; | |
| 391 | } | ||
| 392 | |||
| 393 | /*! | ||
| 394 | * \brief Returns the mass fraction of H2O in the gaseous phase for a given mole fraction. | ||
| 395 | * | ||
| 396 | * \param moleFracH2OInGasRef The reference mole fraction of H2O in the gaseous phase. | ||
| 397 | */ | ||
| 398 | 2 | static Scalar refMoleToRefMassFrac_(Scalar moleFracH2OInGasRef) | |
| 399 | { | ||
| 400 | 2 | const Scalar molarMassWater = FluidSystem::molarMass(H2OIdx); // in kg/mol | |
| 401 | 2 | const Scalar molarMassAir = FluidSystem::molarMass(AirIdx); // in kg/mol | |
| 402 | |||
| 403 | 2 | auto massFracH2OInGasRef = molarMassWater * moleFracH2OInGasRef / (molarMassWater * moleFracH2OInGasRef + molarMassAir * (1 - moleFracH2OInGasRef)); | |
| 404 | 2 | return massFracH2OInGasRef; | |
| 405 | } | ||
| 406 | |||
| 407 | |||
| 408 | std::string name_; | ||
| 409 | |||
| 410 | Scalar initPressure_; | ||
| 411 | Scalar initGasSaturation_; | ||
| 412 | Scalar initSalinity_; | ||
| 413 | |||
| 414 | Scalar pressureLow_, pressureHigh_; | ||
| 415 | Scalar temperatureLow_, temperatureHigh_; | ||
| 416 | int nTemperature_; | ||
| 417 | int nPressure_; | ||
| 418 | |||
| 419 | Scalar time_ = 0.0; | ||
| 420 | Scalar timeStepSize_ = 0.0; | ||
| 421 | static constexpr Scalar eps_ = 1e-6; | ||
| 422 | |||
| 423 | std::vector<Scalar> permeability_; | ||
| 424 | |||
| 425 | Dumux::GnuplotInterface<Scalar> gnuplot_; | ||
| 426 | Dumux::GnuplotInterface<Scalar> gnuplot2_; | ||
| 427 | std::vector<Scalar> x_; | ||
| 428 | std::vector<Scalar> y_; | ||
| 429 | std::vector<Scalar> y2_; | ||
| 430 | }; | ||
| 431 | |||
| 432 | } // end namespace Dumux | ||
| 433 | |||
| 434 | #endif | ||
| 435 |