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