GCC Code Coverage Report

Line	Branch	Exec	Source
1			// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2			//
3			// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
4			// SPDX-License-Identifier: GPL-3.0-or-later
5			//
6			/*!
7			* \file
8			* \ingroup MPNCTests
9			* \brief This file contains the parts of the local residual to
10			* calculate the heat conservation in the thermal non-equilibrium model.
11			*/
12			#ifndef DUMUX_ENERGY_COMBUSTION_LOCAL_RESIDUAL_HH
13			#define DUMUX_ENERGY_COMBUSTION_LOCAL_RESIDUAL_HH
14
15			#include <cmath>
16			#include <dune/common/math.hh>
17			#include <dumux/common/spline.hh>
18			#include <dumux/common/exceptions.hh>
19			#include <dumux/common/properties.hh>
20			#include <dumux/common/numeqvector.hh>
21
22			#include <dumux/porousmediumflow/nonequilibrium/thermal/localresidual.hh>
23
24			namespace Dumux {
25
26			/*!
27			* \ingroup MPNCTests
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			template<class TypeTag>
32			class CombustionEnergyLocalResidual
33			: public EnergyLocalResidualNonEquilibrium<TypeTag, 1/*numEnergyEqFluid*/>
34			{
35			using Scalar = GetPropType<TypeTag, Properties::Scalar>;
36			using NumEqVector = Dumux::NumEqVector<GetPropType<TypeTag, Properties::PrimaryVariables>>;
37			using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
38			using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
39			using SubControlVolume = typename FVElementGeometry::SubControlVolume;
40			using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
41			using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
42			using Element = typename GridView::template Codim<0>::Entity;
43			using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
44			using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
45
46			using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
47			using Indices = typename ModelTraits::Indices;
48
49			static constexpr auto numEnergyEqFluid = ModelTraits::numEnergyEqFluid();
50			static constexpr auto numEnergyEqSolid = ModelTraits::numEnergyEqSolid();
51			static constexpr auto energyEq0Idx = Indices::energyEq0Idx;
52
53			public:
54			/*!
55			* \brief Calculates the source term of the equation
56			*
57			* \param source The source term
58			* \param element An element which contains part of the control volume
59			* \param fvGeometry The finite-volume geometry
60			* \param elemVolVars The volume variables of the current element
61			* \param scv The sub-control volume over which we integrate the source term
62			*/
63		✗	static void computeSourceEnergy(NumEqVector& source,
64			const Element& element,
65			const FVElementGeometry& fvGeometry,
66			const ElementVolumeVariables& elemVolVars,
67			const SubControlVolume &scv)
68			{
69			//specialization for 2 fluid phases
70		✗	const auto& volVars = elemVolVars[scv];
71		✗	const auto& fs = volVars.fluidState() ;
72		✗	const Scalar characteristicLength = volVars.characteristicLength() ;
73
74			//interfacial area
75			// Shi & Wang, Transport in porous media (2011)
76		✗	const Scalar as = volVars.fluidSolidInterfacialArea();
77
78			//temperature fluid is the same for both fluids
79		✗	const Scalar TFluid = volVars.temperatureFluid(0);
80		✗	const Scalar TSolid = volVars.temperatureSolid();
81
82		✗	const Scalar satW = fs.saturation(0) ;
83		✗	const Scalar satN = fs.saturation(1) ;
84
85		✗	const Scalar eps = 1e-6 ;
86			Scalar solidToFluidEnergyExchange ;
87
88			Scalar fluidConductivity ;
89		✗	if (satW < 1.0 - eps)
90			fluidConductivity = volVars.fluidThermalConductivity(1) ;
91		✗	else if (satW >= 1.0 - eps)
92			fluidConductivity = volVars.fluidThermalConductivity(0) ;
93			else
94		✗	DUNE_THROW(Dune::NotImplemented, "wrong range");
95
96		✗	const Scalar factorEnergyTransfer = volVars.factorEnergyTransfer() ;
97
98		✗	solidToFluidEnergyExchange = factorEnergyTransfer * (TSolid - TFluid) / characteristicLength * as * fluidConductivity ;
99		✗	const Scalar epsRegul = 1e-3 ;
100
101		✗	if (satW < (0 + eps) )
102			{
103		✗	solidToFluidEnergyExchange *= volVars.nusseltNumber(1) ;
104			}
105		✗	else if ( (satW >= 0 + eps) and (satW < 1.0-eps) )
106			{
107		✗	solidToFluidEnergyExchange *= (volVars.nusseltNumber(1) * satN );
108			Scalar qBoil ;
109		✗	if (satW<=epsRegul)
110			{// regularize
111			typedef Dumux::Spline<Scalar> Spline;
112		✗	Spline sp(0.0, epsRegul, // x1, x2
113			QBoilFunc(volVars, 0.0), QBoilFunc(volVars, epsRegul), // y1, y2
114			0.0,dQBoil_dSw(volVars, epsRegul)); // m1, m2
115
116		✗	qBoil = sp.eval(satW);
117			}
118
119		✗	else if (satW>= (1.0-epsRegul) )
120			{ // regularize
121			typedef Dumux::Spline<Scalar> Spline;
122		✗	Spline sp(1.0-epsRegul, 1.0, // x1, x2
123			QBoilFunc(volVars, 1.0-epsRegul), 0.0, // y1, y2
124			dQBoil_dSw(volVars, 1.0-epsRegul), 0.0 ); // m1, m2
125
126		✗	qBoil = sp.eval(satW) ;
127			}
128			else
129			{
130		✗	qBoil = QBoilFunc(volVars, satW) ;
131			}
132
133		✗	solidToFluidEnergyExchange += qBoil;
134			}
135		✗	else if (satW >= 1.0-eps)
136			{
137		✗	solidToFluidEnergyExchange *= volVars.nusseltNumber(0) ;
138			}
139			else
140		✗	DUNE_THROW(Dune::NotImplemented, "wrong range");
141
142			using std::isfinite;
143		✗	if (!isfinite(solidToFluidEnergyExchange))
144		✗	DUNE_THROW(NumericalProblem, "Calculated non-finite source, " << "TFluid="<< TFluid << " TSolid="<< TSolid);
145
146		✗	for(int energyEqIdx =0; energyEqIdx<numEnergyEqFluid+numEnergyEqSolid; ++energyEqIdx)
147			{
148		✗	switch (energyEqIdx)
149			{
150		✗	case 0 :
151		✗	source[energyEq0Idx + energyEqIdx] += solidToFluidEnergyExchange;
152		✗	break;
153		✗	case 1 :
154		✗	source[energyEq0Idx + energyEqIdx] -= solidToFluidEnergyExchange;
155		✗	break;
156		✗	default:
157		✗	DUNE_THROW(Dune::NotImplemented,
158			"wrong index");
159			} // end switch
160			}// end energyEqIdx
161		✗	}// end source
162
163
164			/*! \brief Calculates the energy transfer during boiling, i.e. latent heat
165			*
166			* \param volVars The volume variables
167			* \param satW The wetting phase saturation. Not taken from volVars, because we regularize.
168			*/
169		✗	static Scalar QBoilFunc(const VolumeVariables & volVars,
170			const Scalar satW)
171			{
172			// using saturation as input (instead of from volVars)
173			// in order to make regularization (evaluation at different points) easier
174		✗	const auto& fs = volVars.fluidState();
175		✗	const Scalar g( 9.81 );
176		✗	const Scalar gamma(0.0589);
177		✗	const Scalar TSolid = volVars.temperatureSolid();
178			using std::pow;
179			using Dune::power;
180		✗	const Scalar as = volVars.fluidSolidInterfacialArea();
181		✗	const Scalar mul = fs.viscosity(0);
182		✗	const Scalar deltahv = fs.enthalpy(1) - fs.enthalpy(0);
183		✗	const Scalar deltaRho = fs.density(0) - fs.density(1);
184		✗	const Scalar firstBracket = pow(g * deltaRho / gamma, 0.5);
185		✗	const Scalar cp = FluidSystem::heatCapacity(fs, 0);
186			// This use of Tsat is only justified if the fluid is always boiling (Tsat equals boiling conditions)
187			// If a different state is to be simulated, please use the actual fluid temperature instead.
188		✗	const Scalar Tsat = FluidSystem::vaporTemperature(fs, 1 ) ;
189		✗	const Scalar deltaT = TSolid - Tsat;
190		✗	const Scalar secondBracket = power( (cp *deltaT / (0.006 * deltahv) ) , 3);
191		✗	const Scalar Prl = volVars.prandtlNumber(0);
192		✗	const Scalar thirdBracket = pow( 1/Prl , (1.7/0.33));
193		✗	const Scalar QBoil = satW * as * mul * deltahv * firstBracket * secondBracket * thirdBracket;
194		✗	return QBoil;
195			}
196
197			/*! \brief Calculates the derivative of the energy transfer function during boiling. Needed for regularization.
198			*
199			* \param volVars The volume variables
200			* \param satW The wetting phase saturation. Not taken from volVars, because we regularize.
201			*/
202			static Scalar dQBoil_dSw(const VolumeVariables & volVars,
203			const Scalar satW)
204			{
205			// on the fly derive w.r.t. Sw.
206			// Only linearly depending on it (directly)
207		✗	return (QBoilFunc(volVars, satW) / satW ) ;
208			}
209			};
210			} // end namespace Dumux
211
212			#endif
213