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