3.1-git
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
reaction.hh
Go to the documentation of this file.
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 * See the file COPYING for full copying permissions. *
5 * *
6 * This program is free software: you can redistribute it and/or modify *
7 * it under the terms of the GNU General Public License as published by *
8 * the Free Software Foundation, either version 3 of the License, or *
9 * (at your option) any later version. *
10 * *
11 * This program is distributed in the hope that it will be useful, *
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
14 * GNU General Public License for more details. *
15 * *
16 * You should have received a copy of the GNU General Public License *
17 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
18 *****************************************************************************/
27#ifndef DUMUX_THERMOCHEM_REACTION_HH
28#define DUMUX_THERMOCHEM_REACTION_HH
29
30namespace Dumux {
31
38class ThermoChemReaction {
39
40public:
44 template<class VolumeVariables>
45 typename VolumeVariables::PrimaryVariables::value_type
46 thermoChemReaction(const VolumeVariables &volVars) const
47 {
48 using FluidSystem = typename VolumeVariables::FluidSystem;
49 using SolidSystem = typename VolumeVariables::SolidSystem;
50
51 static constexpr auto H2OIdx = FluidSystem::compIdx(FluidSystem::MultiPhaseFluidSystem::H2OIdx);
52 static constexpr int cPhaseIdx = SolidSystem::comp0Idx;
53 static constexpr int hPhaseIdx = SolidSystem::comp1Idx;
54
55 using Scalar = typename VolumeVariables::PrimaryVariables::value_type;
56
57 // calculate the equilibrium temperature Teq
58 Scalar T= volVars.temperature();
59 Scalar Teq = 0;
60
61 Scalar moleFractionVapor = 1e-3;
62
63 if(volVars.moleFraction(0, H2OIdx) > 1e-3)
64 moleFractionVapor = volVars.moleFraction(0, H2OIdx);
65
66 if(volVars.moleFraction(0, H2OIdx) >= 1.0) moleFractionVapor = 1;
67
68 Scalar pV = volVars.pressure(0) *moleFractionVapor;
69 Scalar vaporPressure = pV*1.0e-5;
70 Scalar pFactor = log(vaporPressure);
71
72 Teq = -12845;
73 Teq /= (pFactor - 16.508); //the equilibrium temperature
74
75 Scalar realSolidDensityAverage = (volVars.solidVolumeFraction(hPhaseIdx)*volVars.solidComponentDensity(hPhaseIdx)
76 + volVars.solidVolumeFraction(cPhaseIdx)*volVars.solidComponentDensity(cPhaseIdx))
77 / (volVars.solidVolumeFraction(hPhaseIdx)
78 + volVars.solidVolumeFraction(cPhaseIdx));
79
80 if(realSolidDensityAverage <= volVars.solidComponentDensity(cPhaseIdx))
81 {
82 realSolidDensityAverage = volVars.solidComponentDensity(cPhaseIdx);
83 }
84
85 if(realSolidDensityAverage >= volVars.solidComponentDensity(hPhaseIdx))
86 {
87 realSolidDensityAverage = volVars.solidComponentDensity(hPhaseIdx);
88 }
89
90 Scalar qMass = 0.0;
91
92 // discharge or hydration
93 if (T < Teq){
94 Scalar dXH_dt1 = 0.0;
95 Scalar dXH_dt2 = 0.0;
96
97 Scalar xH = (realSolidDensityAverage-volVars.solidComponentDensity(cPhaseIdx))/(volVars.solidComponentDensity(hPhaseIdx)- volVars.solidComponentDensity(cPhaseIdx));
98
99 if(xH < 1.0e-5) {xH = 1.0e-5; }
100 if(xH >=1 ) {xH = 1 - 1e-5; }
101
102 Scalar R = 8.314 ; // [J/molK]
103 Scalar peq = 1e5*exp( (-12845)/T + 16.508);
104
105 if(peq >= pV) {peq=899954;}
106
107 Scalar dXH_dt = 0;
108
109 // reaction kinetics for T-Teq > 50 K
110 if(Teq -T > 50.25){
111
112 Scalar A =exp(-8.9486e4/(R*T));
113 Scalar B = pow(((pV/peq)-1),0.83);
114 Scalar D = 1-xH;
115 Scalar C = pow((-log(D)),0.666);
116
117 dXH_dt = 1.3945e4*A *B*3*D*C;
118
119 }
120
121 // reaction kinetics for T-Teq < 50 K
122 if(Teq -T < 49.75){
123
124 Scalar E = exp(5.3332e4/T);
125 Scalar F = pow((pV*1e-5),6);
126 Scalar G = (1-xH);
127
128 dXH_dt = 1.004e-34 *E * F * G;
129
130 }
131
132 // linearization of the point of discontinuity
133 if(Teq-T <=50.25 && Teq-T >=49.75){
134
135 Scalar op = ((Teq-T)- 49.75)*2;
136 Scalar A =exp(-8.9486e4/(R*T));
137 Scalar B = pow(((pV/peq)-1),0.83);
138 Scalar D = 1-xH;
139 Scalar C = pow((-log(D)),0.666);
140 dXH_dt1 = 1.3945e4*A *B*3*D*C;
141
142 Scalar E = exp(5.3332e4/T);
143 Scalar F = pow((pV*1e-5),6);
144 Scalar G = (1-xH);
145 dXH_dt2 = 1.004e-34 *E * F * G;
146
147 dXH_dt = dXH_dt1*op + dXH_dt2*(1-op);
148 }
149
150 // no reaction at equilibrium
151 if(Teq -T <= 0)
152 dXH_dt = 0;
153 Scalar deltaRhoS = volVars.solidComponentDensity(hPhaseIdx) - volVars.solidComponentDensity(cPhaseIdx);
154 qMass = dXH_dt*deltaRhoS;
155 }
156
157 return qMass;
158 }
159
160
164 template<class VolumeVariables>
165 typename VolumeVariables::PrimaryVariables::value_type
166 thermoChemReactionSimple(const VolumeVariables &volVars) const
167 {
168 using FluidSystem = typename VolumeVariables::FluidSystem;
169 using SolidSystem = typename VolumeVariables::SolidSystem;
170
171 static constexpr auto H2OIdx = FluidSystem::compIdx(FluidSystem::MultiPhaseFluidSystem::H2OIdx);
172 static constexpr int cPhaseIdx = SolidSystem::comp0Idx;
173 static constexpr int hPhaseIdx = SolidSystem::comp1Idx;
174
175 using Scalar = typename VolumeVariables::PrimaryVariables::value_type;
176
177 // calculate the equilibrium temperature Teq
178 Scalar T= volVars.temperature();
179 Scalar Teq = 0;
180
181 Scalar moleFractionVapor = 1e-3;
182
183 if(volVars.moleFraction(0, H2OIdx) > 1e-3)
184 moleFractionVapor = volVars.moleFraction(0, H2OIdx);
185
186 if(volVars.moleFraction(0, H2OIdx) >= 1.0) moleFractionVapor = 1;
187
188 Scalar pV = volVars.pressure(0) *moleFractionVapor;
189 Scalar vaporPressure = pV*1.0e-5;
190 Scalar pFactor = log(vaporPressure);
191
192 Teq = -12845;
193 Teq /= (pFactor - 16.508); //the equilibrium temperature
194
195
196 Scalar realSolidDensityAverage = (volVars.solidVolumeFraction(hPhaseIdx)*volVars.solidComponentDensity(hPhaseIdx)
197 + volVars.solidVolumeFraction(cPhaseIdx)*volVars.solidComponentDensity(cPhaseIdx))
198 / (volVars.solidVolumeFraction(hPhaseIdx)
199 + volVars.solidVolumeFraction(cPhaseIdx));
200
201 if(realSolidDensityAverage <= volVars.solidComponentDensity(cPhaseIdx))
202 {
203 realSolidDensityAverage = volVars.solidComponentDensity(cPhaseIdx);
204 }
205
206 if(realSolidDensityAverage >= volVars.solidComponentDensity(hPhaseIdx))
207 {
208 realSolidDensityAverage = volVars.solidComponentDensity(hPhaseIdx);
209 }
210
211 Scalar qMass = 0.0;
212
213 // discharge or hydration
214 if (T < Teq){
215 Scalar massFracH2O_fPhase = volVars.massFraction(0, H2OIdx);
216 Scalar krh = 0.2;
217
218 Scalar rHydration = - massFracH2O_fPhase* (volVars.solidComponentDensity(hPhaseIdx)- realSolidDensityAverage)
219 * krh * (T-Teq)/ Teq;
220
221 qMass = rHydration;
222 }
223
224 // charge or hydration
225 else if(T > Teq){
226
227 Scalar krd = 0.05;
228
229 Scalar rDehydration = -(volVars.solidComponentDensity(cPhaseIdx)- realSolidDensityAverage)
230 * krd * (Teq-T)/ Teq;
231
232 qMass = rDehydration;
233 }
234
235 if(Teq -T == 0) qMass = 0;
236
237 return qMass;
238 }
239
240};
241
242} // end namespace Dumux
243
244#endif
Dumux
make the local view function available whenever we use the grid geometry
Definition: adapt.hh:29
Dumux::ThermoChemReaction
Class for the evaluation of the reaction rate of Calciumoxide to Halciumhydroxide.
Definition: reaction.hh:38
Dumux::ThermoChemReaction::thermoChemReactionSimple
VolumeVariables::PrimaryVariables::value_type thermoChemReactionSimple(const VolumeVariables &volVars) const
Evaluates the simple chemical reaction kinetics (see Nagel et al. 2014)
Definition: reaction.hh:166
Dumux::ThermoChemReaction::thermoChemReaction
VolumeVariables::PrimaryVariables::value_type thermoChemReaction(const VolumeVariables &volVars) const
Evaluates the reaction kinetics (see Nagel et al. 2014 ).
Definition: reaction.hh:46
Generated by doxygen 1.9.3