| 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 ThermalNonEquilibriumModel | ||
| 10 | * \brief This file contains the parts of the local residual to | ||
| 11 | * calculate the heat conservation in the thermal non-equilibrium model. | ||
| 12 | */ | ||
| 13 | |||
| 14 | #ifndef DUMUX_ENERGY_NONEQUILIBRIUM_LOCAL_RESIDUAL_HH | ||
| 15 | #define DUMUX_ENERGY_NONEQUILIBRIUM_LOCAL_RESIDUAL_HH | ||
| 16 | |||
| 17 | #include <cmath> | ||
| 18 | #include <dumux/common/spline.hh> | ||
| 19 | #include <dumux/common/exceptions.hh> | ||
| 20 | #include <dumux/common/properties.hh> | ||
| 21 | #include <dumux/common/typetraits/typetraits.hh> | ||
| 22 | #include <dumux/common/numeqvector.hh> | ||
| 23 | #include <dumux/flux/referencesystemformulation.hh> | ||
| 24 | |||
| 25 | namespace Dumux { | ||
| 26 | |||
| 27 | /*! | ||
| 28 | * \ingroup ThermalNonEquilibriumModel | ||
| 29 | * \brief This file contains the parts of the local residual to | ||
| 30 | * calculate the heat conservation in the thermal non-equilibrium model. | ||
| 31 | */ | ||
| 32 | // forward declaration | ||
| 33 | template <class TypeTag, int numEnergyEqFluid> | ||
| 34 | class EnergyLocalResidualNonEquilibrium; | ||
| 35 | |||
| 36 | template<class TypeTag> | ||
| 37 | class EnergyLocalResidualNonEquilibrium<TypeTag, 1/*numEnergyEqFluid*/> | ||
| 38 | { | ||
| 39 | using Scalar = GetPropType<TypeTag, Properties::Scalar>; | ||
| 40 | using NumEqVector = Dumux::NumEqVector<GetPropType<TypeTag, Properties::PrimaryVariables>>; | ||
| 41 | using Problem = GetPropType<TypeTag, Properties::Problem>; | ||
| 42 | using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>; | ||
| 43 | using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView; | ||
| 44 | using SubControlVolume = typename FVElementGeometry::SubControlVolume; | ||
| 45 | using FluxVariables = GetPropType<TypeTag, Properties::FluxVariables>; | ||
| 46 | using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>; | ||
| 47 | using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView; | ||
| 48 | using Element = typename GridView::template Codim<0>::Entity; | ||
| 49 | using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView; | ||
| 50 | using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; | ||
| 51 | |||
| 52 | using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>; | ||
| 53 | using Indices = typename ModelTraits::Indices; | ||
| 54 | |||
| 55 | static constexpr auto numEnergyEqFluid = ModelTraits::numEnergyEqFluid(); | ||
| 56 | static constexpr auto numEnergyEqSolid = ModelTraits::numEnergyEqSolid(); | ||
| 57 | static constexpr auto energyEq0Idx = Indices::energyEq0Idx; | ||
| 58 | static constexpr auto energyEqSolidIdx = Indices::energyEqSolidIdx; | ||
| 59 | |||
| 60 | static constexpr auto numPhases = ModelTraits::numFluidPhases(); | ||
| 61 | static constexpr auto numComponents = ModelTraits::numFluidComponents(); | ||
| 62 | |||
| 63 | public: | ||
| 64 | template <typename T = void> | ||
| 65 | static void fluidPhaseStorage(NumEqVector& storage, | ||
| 66 | const SubControlVolume& scv, | ||
| 67 | const VolumeVariables& volVars, | ||
| 68 | int phaseIdx) | ||
| 69 | { | ||
| 70 | static_assert(AlwaysFalse<T>::value, "Deprecated interface that has been removed! Use new interface with additional argument problem instead. Will be entirely removed after release 3.10."); | ||
| 71 | } | ||
| 72 | |||
| 73 | //! The energy storage in the fluid phase with index phaseIdx | ||
| 74 | 5010400 | static void fluidPhaseStorage(NumEqVector& storage, | |
| 75 | const Problem&, | ||
| 76 | const SubControlVolume& scv, | ||
| 77 | const VolumeVariables& volVars, | ||
| 78 | int phaseIdx) | ||
| 79 | { | ||
| 80 | // in case we have one energy equation for more than one fluid phase, | ||
| 81 | // add up parts on the one energy equation | ||
| 82 | 5010400 | storage[energyEq0Idx] += volVars.porosity() | |
| 83 | 5010400 | * volVars.density(phaseIdx) | |
| 84 | 5010400 | * volVars.internalEnergy(phaseIdx) | |
| 85 | 5010400 | * volVars.saturation(phaseIdx); | |
| 86 | |||
| 87 | } | ||
| 88 | |||
| 89 | |||
| 90 | //! The energy storage in the solid matrix | ||
| 91 | static void solidPhaseStorage(NumEqVector& storage, | ||
| 92 | const SubControlVolume& scv, | ||
| 93 | const VolumeVariables& volVars) | ||
| 94 | { | ||
| 95 | // heat conduction for the fluid phases | ||
| 96 | 9734560 | for(int sPhaseIdx = 0; sPhaseIdx<numEnergyEqSolid; ++sPhaseIdx) | |
| 97 | { | ||
| 98 | 9734560 | storage[energyEqSolidIdx+sPhaseIdx] += volVars.temperatureSolid() | |
| 99 | 9734560 | * volVars.solidHeatCapacity() | |
| 100 | 9734560 | * volVars.solidDensity() | |
| 101 | 9734560 | * (1.0 - volVars.porosity()); | |
| 102 | } | ||
| 103 | } | ||
| 104 | |||
| 105 | /*! | ||
| 106 | * \brief The dispersive energy fluxes | ||
| 107 | * | ||
| 108 | * \param flux The flux | ||
| 109 | * \param fluxVars The flux variables. | ||
| 110 | */ | ||
| 111 | static void heatDispersionFlux(NumEqVector& flux, | ||
| 112 | FluxVariables& fluxVars) | ||
| 113 | {} | ||
| 114 | |||
| 115 | //! The advective phase energy fluxes | ||
| 116 | 3000710 | static void heatConvectionFlux(NumEqVector& flux, | |
| 117 | FluxVariables& fluxVars, | ||
| 118 | int phaseIdx) | ||
| 119 | { | ||
| 120 | 6001420 | auto upwindTerm = [phaseIdx](const auto& volVars) | |
| 121 | 3000710 | { return volVars.density(phaseIdx)*volVars.mobility(phaseIdx)*volVars.enthalpy(phaseIdx); }; | |
| 122 | |||
| 123 | //in case we have one energy equation for more than one fluid phase, add up advective parts on the one energy equation | ||
| 124 | 3000710 | flux[energyEq0Idx] += fluxVars.advectiveFlux(phaseIdx, upwindTerm); | |
| 125 | |||
| 126 | //now add the diffusive part | ||
| 127 | 3000710 | const auto diffusiveFluxes = fluxVars.molecularDiffusionFlux(phaseIdx); | |
| 128 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2339200 times.
|
3000710 | const auto& elemVolVars = fluxVars.elemVolVars(); |
| 129 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2339200 times.
|
3000710 | const auto& scvf = fluxVars.scvFace(); |
| 130 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2339200 times.
|
3000710 | const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()]; |
| 131 | 3000710 | const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()]; | |
| 132 | |||
| 133 |
2/2✓ Branch 0 taken 6001420 times.
✓ Branch 1 taken 3000710 times.
|
9002130 | for (int compIdx = 0; compIdx < numComponents; ++compIdx) |
| 134 | { | ||
| 135 | //no diffusion of the main component, this is a hack to use normal fick's law which computes both diffusions (main and component). We only add the part from the component here | ||
| 136 |
2/2✓ Branch 0 taken 3000710 times.
✓ Branch 1 taken 3000710 times.
|
6001420 | if (phaseIdx == compIdx) |
| 137 | 3000710 | continue; | |
| 138 | //we need the upwind enthalpy. Even better would be the componentEnthalpy | ||
| 139 | 3000710 | auto enthalpy = 0.0; | |
| 140 |
2/2✓ Branch 0 taken 1776390 times.
✓ Branch 1 taken 1224320 times.
|
3000710 | if (diffusiveFluxes[compIdx] > 0) |
| 141 | 1776390 | enthalpy += insideVolVars.enthalpy(phaseIdx); | |
| 142 | else | ||
| 143 | 1224320 | enthalpy += outsideVolVars.enthalpy(phaseIdx); | |
| 144 | |||
| 145 | //check for the reference system and adapt units of the diffusive flux accordingly. | ||
| 146 | if (FluxVariables::MolecularDiffusionType::referenceSystemFormulation() == ReferenceSystemFormulation::massAveraged) | ||
| 147 | 3000710 | flux[energyEq0Idx] += diffusiveFluxes[compIdx]*enthalpy; | |
| 148 | else | ||
| 149 | flux[energyEq0Idx] += diffusiveFluxes[compIdx]*FluidSystem::molarMass(compIdx)*enthalpy; | ||
| 150 | } | ||
| 151 | 3000710 | } | |
| 152 | |||
| 153 | //! The diffusive energy fluxes | ||
| 154 | 3000710 | static void heatConductionFlux(NumEqVector& flux, | |
| 155 | FluxVariables& fluxVars) | ||
| 156 | { | ||
| 157 | //in case we have one energy equation for more than one fluid phase we use an effective law in the nonequilibrium fourierslaw | ||
| 158 | 3000710 | flux[energyEq0Idx] += fluxVars.heatConductionFlux(0); | |
| 159 | //heat conduction for the solid phases | ||
| 160 |
2/2✓ Branch 0 taken 3000710 times.
✓ Branch 1 taken 3000710 times.
|
6001420 | for(int sPhaseIdx = 0; sPhaseIdx<numEnergyEqSolid; ++sPhaseIdx) |
| 161 | 3000710 | flux[energyEqSolidIdx+sPhaseIdx] += fluxVars.heatConductionFlux(numPhases + sPhaseIdx); | |
| 162 | 3000710 | } | |
| 163 | |||
| 164 | /*! | ||
| 165 | * \brief heat transfer between the phases for nonequilibrium models | ||
| 166 | * | ||
| 167 | * \param source The source which ought to be simulated | ||
| 168 | * \param element An element which contains part of the control volume | ||
| 169 | * \param fvGeometry The finite-volume geometry | ||
| 170 | * \param elemVolVars The volume variables of the current element | ||
| 171 | * \param scv The sub-control volume over which we integrate the source term | ||
| 172 | */ | ||
| 173 | 2505200 | static void computeSourceEnergy(NumEqVector& source, | |
| 174 | const Element& element, | ||
| 175 | const FVElementGeometry& fvGeometry, | ||
| 176 | const ElementVolumeVariables& elemVolVars, | ||
| 177 | const SubControlVolume &scv) | ||
| 178 | { | ||
| 179 | // specialization for 2 fluid phases | ||
| 180 | 2505200 | const auto& volVars = elemVolVars[scv]; | |
| 181 | 2505200 | const Scalar characteristicLength = volVars.characteristicLength() ; | |
| 182 | |||
| 183 | // interfacial area | ||
| 184 | // Shi & Wang, Transport in porous media (2011) | ||
| 185 | 2505200 | const Scalar as = volVars.fluidSolidInterfacialArea(); | |
| 186 | |||
| 187 | // temperature fluid is the same for both fluids | ||
| 188 | 2505200 | const Scalar TFluid = volVars.temperatureFluid(0); | |
| 189 | 2505200 | const Scalar TSolid = volVars.temperatureSolid(); | |
| 190 | |||
| 191 | Scalar solidToFluidEnergyExchange ; | ||
| 192 | |||
| 193 | 2505200 | const Scalar fluidConductivity = volVars.fluidThermalConductivity(0) ; | |
| 194 | |||
| 195 | 2505200 | const Scalar factorEnergyTransfer = volVars.factorEnergyTransfer() ; | |
| 196 | |||
| 197 | 2505200 | solidToFluidEnergyExchange = factorEnergyTransfer * (TSolid - TFluid) / characteristicLength * as * fluidConductivity; | |
| 198 | |||
| 199 | 2505200 | solidToFluidEnergyExchange *= volVars.nusseltNumber(0); | |
| 200 | |||
| 201 |
2/2✓ Branch 0 taken 5010400 times.
✓ Branch 1 taken 2505200 times.
|
7515600 | for(int energyEqIdx = 0; energyEqIdx < numEnergyEqFluid+numEnergyEqSolid; ++energyEqIdx) |
| 202 | { | ||
| 203 |
2/2✓ Branch 0 taken 2505200 times.
✓ Branch 1 taken 2505200 times.
|
5010400 | switch (energyEqIdx) |
| 204 | { | ||
| 205 | 2505200 | case 0 : | |
| 206 | 2505200 | source[energyEq0Idx + energyEqIdx] += solidToFluidEnergyExchange; | |
| 207 | 2505200 | break; | |
| 208 | 2505200 | case 1 : | |
| 209 | 2505200 | source[energyEq0Idx + energyEqIdx] -= solidToFluidEnergyExchange; | |
| 210 | 2505200 | break; | |
| 211 | default: | ||
| 212 | 5010400 | DUNE_THROW(Dune::NotImplemented, | |
| 213 | "wrong index"); | ||
| 214 | } // end switch | ||
| 215 | } // end energyEqIdx | ||
| 216 | 2505200 | } // end source | |
| 217 | }; | ||
| 218 | |||
| 219 | /*! | ||
| 220 | * \ingroup ThermalNonEquilibriumModel | ||
| 221 | * \brief TODO docme | ||
| 222 | */ | ||
| 223 | template<class TypeTag> | ||
| 224 | class EnergyLocalResidualNonEquilibrium<TypeTag, 2/*numEnergyEqFluid*/> | ||
| 225 | : public EnergyLocalResidualNonEquilibrium<TypeTag, 1/*numEnergyEqFluid*/> | ||
| 226 | { | ||
| 227 | using Scalar = GetPropType<TypeTag, Properties::Scalar>; | ||
| 228 | using NumEqVector = Dumux::NumEqVector<GetPropType<TypeTag, Properties::PrimaryVariables>>; | ||
| 229 | using Problem = GetPropType<TypeTag, Properties::Problem>; | ||
| 230 | using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>; | ||
| 231 | using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView; | ||
| 232 | using SubControlVolume = typename FVElementGeometry::SubControlVolume; | ||
| 233 | using FluxVariables = GetPropType<TypeTag, Properties::FluxVariables>; | ||
| 234 | using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>; | ||
| 235 | using SolidSystem = GetPropType<TypeTag, Properties::SolidSystem>; | ||
| 236 | using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView; | ||
| 237 | using Element = typename GridView::template Codim<0>::Entity; | ||
| 238 | using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView; | ||
| 239 | using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; | ||
| 240 | |||
| 241 | using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>; | ||
| 242 | using Indices = typename ModelTraits::Indices; | ||
| 243 | |||
| 244 | static constexpr auto numPhases = ModelTraits::numFluidPhases(); | ||
| 245 | static constexpr auto numEnergyEqFluid = ModelTraits::numEnergyEqFluid(); | ||
| 246 | static constexpr auto numEnergyEqSolid = ModelTraits::numEnergyEqSolid(); | ||
| 247 | static constexpr int energyEq0Idx = Indices::energyEq0Idx; | ||
| 248 | static constexpr int energyEqSolidIdx = Indices::energyEqSolidIdx; | ||
| 249 | static constexpr int conti0EqIdx = Indices::conti0EqIdx; | ||
| 250 | |||
| 251 | static constexpr auto numComponents = ModelTraits::numFluidComponents(); | ||
| 252 | static constexpr int phase0Idx = FluidSystem::phase0Idx; | ||
| 253 | static constexpr int phase1Idx = FluidSystem::phase1Idx; | ||
| 254 | static constexpr int sPhaseIdx = numPhases; | ||
| 255 | |||
| 256 | static constexpr bool enableChemicalNonEquilibrium = ModelTraits::enableChemicalNonEquilibrium(); | ||
| 257 | |||
| 258 | public: | ||
| 259 | template <typename T = void> | ||
| 260 | static void fluidPhaseStorage(NumEqVector& storage, | ||
| 261 | const SubControlVolume& scv, | ||
| 262 | const VolumeVariables& volVars, | ||
| 263 | int phaseIdx) | ||
| 264 | { | ||
| 265 | static_assert(AlwaysFalse<T>::value, "Deprecated interface that has been removed! Use new interface with additional argument problem instead. Will be entirely removed after release 3.10."); | ||
| 266 | } | ||
| 267 | |||
| 268 | //! The energy storage in the fluid phase with index phaseIdx | ||
| 269 | 9448320 | static void fluidPhaseStorage(NumEqVector& storage, | |
| 270 | const Problem&, | ||
| 271 | const SubControlVolume& scv, | ||
| 272 | const VolumeVariables& volVars, | ||
| 273 | int phaseIdx) | ||
| 274 | { | ||
| 275 | 9448320 | storage[energyEq0Idx+phaseIdx] += volVars.porosity() | |
| 276 | 9448320 | * volVars.density(phaseIdx) | |
| 277 | 9448320 | * volVars.internalEnergy(phaseIdx) | |
| 278 | 9448320 | * volVars.saturation(phaseIdx); | |
| 279 | |||
| 280 | } | ||
| 281 | |||
| 282 | //! The advective phase energy fluxes | ||
| 283 | 4724160 | static void heatConvectionFlux(NumEqVector& flux, | |
| 284 | FluxVariables& fluxVars, | ||
| 285 | int phaseIdx) | ||
| 286 | { | ||
| 287 | 9448320 | auto upwindTerm = [phaseIdx](const auto& volVars) | |
| 288 | 4724160 | { return volVars.density(phaseIdx)*volVars.mobility(phaseIdx)*volVars.enthalpy(phaseIdx); }; | |
| 289 | |||
| 290 | // in case we have one energy equation for more than one fluid phase, add up advective parts on the one energy equation | ||
| 291 | 4724160 | flux[energyEq0Idx+phaseIdx] += fluxVars.advectiveFlux(phaseIdx, upwindTerm); | |
| 292 | |||
| 293 | // add the diffusiv part | ||
| 294 | 4724160 | const auto diffusiveFluxes = fluxVars.molecularDiffusionFlux(phaseIdx); | |
| 295 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4724160 times.
|
4724160 | const auto& elemVolVars = fluxVars.elemVolVars(); |
| 296 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4724160 times.
|
4724160 | const auto& scvf = fluxVars.scvFace(); |
| 297 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4724160 times.
|
4724160 | const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()]; |
| 298 | 4724160 | const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()]; | |
| 299 | |||
| 300 |
2/2✓ Branch 0 taken 9448320 times.
✓ Branch 1 taken 4724160 times.
|
14172480 | for (int compIdx = 0; compIdx < numComponents; ++compIdx) |
| 301 | { | ||
| 302 | // no diffusion of the main component, this is a hack to use normal fick's law which computes both diffusions (main and component). We only add the part from the component here | ||
| 303 |
2/2✓ Branch 0 taken 4724160 times.
✓ Branch 1 taken 4724160 times.
|
9448320 | if (phaseIdx == compIdx) |
| 304 | 4724160 | continue; | |
| 305 | // we need the upwind enthalpy. Even better would be the componentEnthalpy | ||
| 306 | 4724160 | auto enthalpy = 0.0; | |
| 307 |
2/2✓ Branch 0 taken 2373901 times.
✓ Branch 1 taken 2350259 times.
|
4724160 | if (diffusiveFluxes[compIdx] > 0) |
| 308 | 2373901 | enthalpy += insideVolVars.enthalpy(phaseIdx); | |
| 309 | else | ||
| 310 | 2350259 | enthalpy += outsideVolVars.enthalpy(phaseIdx); | |
| 311 | 4724160 | flux[energyEq0Idx+phaseIdx] += diffusiveFluxes[compIdx]*FluidSystem::molarMass(compIdx)*enthalpy; | |
| 312 | } | ||
| 313 | 4724160 | } | |
| 314 | |||
| 315 | //! The diffusive energy fluxes | ||
| 316 | 2362080 | static void heatConductionFlux(NumEqVector& flux, | |
| 317 | FluxVariables& fluxVars) | ||
| 318 | { | ||
| 319 |
2/2✓ Branch 0 taken 4724160 times.
✓ Branch 1 taken 2362080 times.
|
7086240 | for(int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) |
| 320 | { | ||
| 321 | 4724160 | flux[energyEq0Idx+phaseIdx] += fluxVars.heatConductionFlux(phaseIdx); | |
| 322 | } | ||
| 323 |
2/2✓ Branch 0 taken 2362080 times.
✓ Branch 1 taken 2362080 times.
|
4724160 | for(int sPhaseIdx=0; sPhaseIdx<numEnergyEqSolid; ++sPhaseIdx) |
| 324 | { | ||
| 325 | 2362080 | flux[energyEqSolidIdx+sPhaseIdx] += fluxVars.heatConductionFlux(numPhases + sPhaseIdx); | |
| 326 | } | ||
| 327 | 2362080 | } | |
| 328 | |||
| 329 | /*! | ||
| 330 | * \brief The dispersive energy fluxes | ||
| 331 | * | ||
| 332 | * \param flux The flux | ||
| 333 | * \param fluxVars The flux variables. | ||
| 334 | */ | ||
| 335 | static void heatDispersionFlux(NumEqVector& flux, | ||
| 336 | FluxVariables& fluxVars) | ||
| 337 | {} | ||
| 338 | |||
| 339 | /*! | ||
| 340 | * \brief Calculates the source term of the equation. | ||
| 341 | * | ||
| 342 | * \param source The source term which ought to be simulated | ||
| 343 | * \param element An element which contains part of the control volume | ||
| 344 | * \param fvGeometry The finite-volume geometry | ||
| 345 | * \param elemVolVars The volume variables of the current element | ||
| 346 | * \param scv The sub-control volume over which we integrate the source term | ||
| 347 | */ | ||
| 348 | 2362080 | static void computeSourceEnergy(NumEqVector& source, | |
| 349 | const Element& element, | ||
| 350 | const FVElementGeometry& fvGeometry, | ||
| 351 | const ElementVolumeVariables& elemVolVars, | ||
| 352 | const SubControlVolume &scv) | ||
| 353 | { | ||
| 354 | // specialization for 2 fluid phases | ||
| 355 | 2362080 | const auto &volVars = elemVolVars[scv]; | |
| 356 | |||
| 357 | 2362080 | const Scalar areaWN = volVars.interfacialArea(phase0Idx, phase1Idx); | |
| 358 | 2362080 | const Scalar areaWS = volVars.interfacialArea(phase0Idx, sPhaseIdx); | |
| 359 | 2362080 | const Scalar areaNS = volVars.interfacialArea(phase1Idx, sPhaseIdx); | |
| 360 | |||
| 361 | 2362080 | const Scalar Tw = volVars.temperatureFluid(phase0Idx); | |
| 362 | 2362080 | const Scalar Tn = volVars.temperatureFluid(phase1Idx); | |
| 363 | 2362080 | const Scalar Ts = volVars.temperatureSolid(); | |
| 364 | |||
| 365 | 2362080 | const Scalar lambdaWetting = volVars.fluidThermalConductivity(phase0Idx); | |
| 366 | 2362080 | const Scalar lambdaNonwetting = volVars.fluidThermalConductivity(phase1Idx); | |
| 367 |
1/2✓ Branch 0 taken 2362080 times.
✗ Branch 1 not taken.
|
2362080 | const Scalar lambdaSolid = volVars.solidThermalConductivity(); |
| 368 | |||
| 369 |
1/2✓ Branch 0 taken 2362080 times.
✗ Branch 1 not taken.
|
2362080 | const Scalar lambdaWN = harmonicMean(lambdaWetting, lambdaNonwetting); |
| 370 |
1/2✓ Branch 0 taken 2362080 times.
✗ Branch 1 not taken.
|
2362080 | const Scalar lambdaWS = harmonicMean(lambdaWetting, lambdaSolid); |
| 371 |
1/2✓ Branch 0 taken 2362080 times.
✗ Branch 1 not taken.
|
2362080 | const Scalar lambdaNS = harmonicMean(lambdaNonwetting, lambdaSolid); |
| 372 | |||
| 373 |
1/2✓ Branch 0 taken 2362080 times.
✗ Branch 1 not taken.
|
2362080 | const Scalar characteristicLength = volVars.characteristicLength() ; |
| 374 |
1/2✓ Branch 0 taken 2362080 times.
✗ Branch 1 not taken.
|
2362080 | const Scalar factorEnergyTransfer = volVars.factorEnergyTransfer() ; |
| 375 | |||
| 376 |
1/2✓ Branch 0 taken 2362080 times.
✗ Branch 1 not taken.
|
4724160 | const Scalar nusseltWN = harmonicMean(volVars.nusseltNumber(phase0Idx), volVars.nusseltNumber(phase1Idx)); |
| 377 | 2362080 | const Scalar nusseltWS = volVars.nusseltNumber(phase0Idx); | |
| 378 | 2362080 | const Scalar nusseltNS = volVars.nusseltNumber(phase1Idx); | |
| 379 | |||
| 380 | 2362080 | const Scalar wettingToNonwettingEnergyExchange = factorEnergyTransfer * (Tw - Tn) / characteristicLength * areaWN * lambdaWN * nusseltWN ; | |
| 381 | 2362080 | const Scalar wettingToSolidEnergyExchange = factorEnergyTransfer * (Tw - Ts) / characteristicLength * areaWS * lambdaWS * nusseltWS ; | |
| 382 | 2362080 | const Scalar nonwettingToSolidEnergyExchange = factorEnergyTransfer * (Tn - Ts) / characteristicLength * areaNS * lambdaNS * nusseltNS ; | |
| 383 | |||
| 384 |
2/2✓ Branch 0 taken 7086240 times.
✓ Branch 1 taken 2362080 times.
|
9448320 | for(int phaseIdx = 0; phaseIdx < numEnergyEqFluid+numEnergyEqSolid; ++phaseIdx) |
| 385 | { | ||
| 386 |
3/3✓ Branch 0 taken 2362080 times.
✓ Branch 1 taken 2362080 times.
✓ Branch 2 taken 2362080 times.
|
7086240 | switch (phaseIdx) |
| 387 | { | ||
| 388 | 2362080 | case phase0Idx: | |
| 389 | 2362080 | source[energyEq0Idx + phaseIdx] += ( - wettingToNonwettingEnergyExchange - wettingToSolidEnergyExchange); | |
| 390 | 2362080 | break; | |
| 391 | 2362080 | case phase1Idx: | |
| 392 | 2362080 | source[energyEq0Idx + phaseIdx] += (+ wettingToNonwettingEnergyExchange - nonwettingToSolidEnergyExchange); | |
| 393 | 2362080 | break; | |
| 394 | 2362080 | case sPhaseIdx: | |
| 395 | 2362080 | source[energyEq0Idx + phaseIdx] += (+ wettingToSolidEnergyExchange + nonwettingToSolidEnergyExchange); | |
| 396 | 2362080 | break; | |
| 397 | default: | ||
| 398 | 7086240 | DUNE_THROW(Dune::NotImplemented, | |
| 399 | "wrong index"); | ||
| 400 | } // end switch | ||
| 401 | |||
| 402 | |||
| 403 | using std::isfinite; | ||
| 404 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 7086240 times.
|
7086240 | if (!isfinite(source[energyEq0Idx + phaseIdx])) |
| 405 |
0/36✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✗ Branch 40 not taken.
✗ Branch 41 not taken.
✗ Branch 43 not taken.
✗ Branch 44 not taken.
✗ Branch 46 not taken.
✗ Branch 47 not taken.
✗ Branch 49 not taken.
✗ Branch 50 not taken.
✗ Branch 53 not taken.
✗ Branch 54 not taken.
|
7086240 | DUNE_THROW(NumericalProblem, "Calculated non-finite source, " << "Tw="<< Tw << " Tn="<< Tn<< " Ts="<< Ts); |
| 406 | }// end phases | ||
| 407 | |||
| 408 | // we only need to do this for when there is more than 1 fluid phase | ||
| 409 | if (enableChemicalNonEquilibrium) | ||
| 410 | { | ||
| 411 | // Here comes the catch: We are not doing energy conservation for the whole | ||
| 412 | // system, but rather for each individual phase. | ||
| 413 | // -> Therefore the energy fluxes over each phase boundary need be | ||
| 414 | // individually accounted for. | ||
| 415 | // -> Each particle crossing a phase boundary does carry some mass and | ||
| 416 | // thus energy! | ||
| 417 | // -> Therefore, this contribution needs to be added. | ||
| 418 | // -> the particle always brings the energy of the originating phase. | ||
| 419 | // -> Energy advectivly transported into a phase = the moles of a component that go into a phase | ||
| 420 | // * molMass * enthalpy of the component in the *originating* phase | ||
| 421 | |||
| 422 | const auto& fluidState = volVars.fluidState(); | ||
| 423 | |||
| 424 |
2/2✓ Branch 0 taken 7086240 times.
✓ Branch 1 taken 2362080 times.
|
9448320 | for(int phaseIdx = 0; phaseIdx < numEnergyEqFluid+numEnergyEqSolid; ++phaseIdx) |
| 425 | { | ||
| 426 |
3/3✓ Branch 0 taken 2362080 times.
✓ Branch 1 taken 2362080 times.
✓ Branch 2 taken 2362080 times.
|
7086240 | switch (phaseIdx) |
| 427 | { | ||
| 428 | case phase0Idx: | ||
| 429 | // sum up the transferred energy by the components into the wetting phase | ||
| 430 |
2/2✓ Branch 0 taken 4724160 times.
✓ Branch 1 taken 2362080 times.
|
7086240 | for(int compIdx = 0; compIdx < numComponents; ++compIdx) |
| 431 | { | ||
| 432 | 4724160 | const unsigned int eqIdx = conti0EqIdx + compIdx + phaseIdx*numComponents; | |
| 433 | 4724160 | source[energyEq0Idx + phaseIdx] += (source[eqIdx] | |
| 434 | 4724160 | * FluidSystem::molarMass(compIdx) | |
| 435 | 4724160 | * FluidSystem::componentEnthalpy(fluidState, phase1Idx, compIdx) ); | |
| 436 | } | ||
| 437 | break; | ||
| 438 | case phase1Idx: | ||
| 439 | // sum up the transferred energy by the components into the nonwetting phase | ||
| 440 |
2/2✓ Branch 0 taken 4724160 times.
✓ Branch 1 taken 2362080 times.
|
7086240 | for(int compIdx =0; compIdx<numComponents; ++compIdx) |
| 441 | { | ||
| 442 | 4724160 | const unsigned int eqIdx = conti0EqIdx + compIdx + phaseIdx*numComponents; | |
| 443 | 4724160 | source[energyEq0Idx + phaseIdx] += (source[eqIdx] | |
| 444 | 4724160 | * FluidSystem::molarMass(compIdx) | |
| 445 | 4724160 | *FluidSystem::componentEnthalpy(fluidState, phase0Idx, compIdx)); | |
| 446 | } | ||
| 447 | break; | ||
| 448 | case sPhaseIdx: | ||
| 449 | break; // no sorption | ||
| 450 | default: | ||
| 451 | 7086240 | DUNE_THROW(Dune::NotImplemented, | |
| 452 | "wrong index"); | ||
| 453 | } // end switch | ||
| 454 | } // end phases | ||
| 455 | } // EnableChemicalNonEquilibrium | ||
| 456 | 2362080 | } // end source | |
| 457 | }; | ||
| 458 | } // end namespace Dumux | ||
| 459 | |||
| 460 | #endif | ||
| 461 |