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