3.2-git
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
porousmediumflow/3pwateroil/volumevariables.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 *****************************************************************************/
26#ifndef DUMUX_3P2CNI_VOLUME_VARIABLES_HH
27#define DUMUX_3P2CNI_VOLUME_VARIABLES_HH
28
29#include <vector>
30#include <iostream>
31
32#include <dumux/common/math.hh>
35
39
42
44
45namespace Dumux {
46
52template <class Traits>
55, public EnergyVolumeVariables<Traits, ThreePWaterOilVolumeVariables<Traits> >
56{
59 using Scalar = typename Traits::PrimaryVariables::value_type;
60 using ModelTraits = typename Traits::ModelTraits;
61 using FS = typename Traits::FluidSystem;
62 static constexpr int numFluidComps = ParentType::numFluidComponents();
63
64 enum {
66
67 wCompIdx = FS::wCompIdx,
68 nCompIdx = FS::nCompIdx,
69
70 wPhaseIdx = FS::wPhaseIdx,
71 gPhaseIdx = FS::gPhaseIdx,
72 nPhaseIdx = FS::nPhaseIdx,
73
74 switch1Idx = ModelTraits::Indices::switch1Idx,
75 switch2Idx = ModelTraits::Indices::switch2Idx,
76 pressureIdx = ModelTraits::Indices::pressureIdx
77 };
78
79 // present phases
80 enum {
81 threePhases = ModelTraits::Indices::threePhases,
82 wPhaseOnly = ModelTraits::Indices::wPhaseOnly,
83 gnPhaseOnly = ModelTraits::Indices::gnPhaseOnly,
84 wnPhaseOnly = ModelTraits::Indices::wnPhaseOnly,
85 gPhaseOnly = ModelTraits::Indices::gPhaseOnly,
86 wgPhaseOnly = ModelTraits::Indices::wgPhaseOnly
87 };
88
89 using EffDiffModel = typename Traits::EffectiveDiffusivityModel;
90 using DiffusionCoefficients = typename Traits::DiffusionType::DiffusionCoefficientsContainer;
91
92public:
94 using FluidState = typename Traits::FluidState;
96 using FluidSystem = typename Traits::FluidSystem;
98 using Indices = typename ModelTraits::Indices;
100 using SolidState = typename Traits::SolidState;
102 using SolidSystem = typename Traits::SolidSystem;
106 static constexpr bool onlyGasPhaseCanDisappear()
107 { return Traits::ModelTraits::onlyGasPhaseCanDisappear(); }
108
117 template<class ElemSol, class Problem, class Element, class Scv>
118 void update(const ElemSol &elemSol,
119 const Problem &problem,
120 const Element &element,
121 const Scv& scv)
122 {
123 ParentType::update(elemSol, problem, element, scv);
124 const auto& priVars = elemSol[scv.localDofIndex()];
125 const auto phasePresence = priVars.state();
126
127 // capillary pressure parameters
128 using MaterialLaw = typename Problem::SpatialParams::MaterialLaw;
129 const auto& materialParams = problem.spatialParams().materialLawParams(element, scv, elemSol);
130
132 {
133 /* first the saturations */
134 if (phasePresence == threePhases)
135 {
136 sw_ = priVars[switch1Idx];
137 sn_ = priVars[switch2Idx];
138 sg_ = 1. - sw_ - sn_;
139 }
140 else if (phasePresence == wPhaseOnly)
141 {
142 sw_ = 1.;
143 sn_ = 0.;
144 sg_ = 0.;
145 }
146 else if (phasePresence == gnPhaseOnly)
147 {
148 sw_ = 0.;
149 sn_ = priVars[switch1Idx];
150 sg_ = 1. - sn_;
151 }
152 else if (phasePresence == wnPhaseOnly)
153 {
154 sn_ = priVars[switch2Idx];
155 sw_ = 1. - sn_;
156 sg_ = 0.;
157 }
158 else if (phasePresence == gPhaseOnly)
159 {
160 sw_ = 0.;
161 sn_ = 0.;
162 sg_ = 1.;
163 }
164 else if (phasePresence == wgPhaseOnly)
165 {
166 sw_ = priVars[switch1Idx];
167 sn_ = 0.;
168 sg_ = 1. - sw_;
169 }
170 else DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
172
173 fluidState_.setSaturation(wPhaseIdx, sw_);
174 fluidState_.setSaturation(gPhaseIdx, sg_);
175 fluidState_.setSaturation(nPhaseIdx, sn_);
176
177 /* now the pressures */
178 if (phasePresence == threePhases || phasePresence == gnPhaseOnly || phasePresence == gPhaseOnly || phasePresence == wgPhaseOnly)
179 {
180 pg_ = priVars[pressureIdx];
181
182 // calculate capillary pressures
183 Scalar pcgw = MaterialLaw::pcgw(materialParams, sw_);
184 Scalar pcnw = MaterialLaw::pcnw(materialParams, sw_);
185 Scalar pcgn = MaterialLaw::pcgn(materialParams, sw_ + sn_);
186
187 Scalar pcAlpha = MaterialLaw::pcAlpha(materialParams, sn_);
188 Scalar pcNW1 = 0.0; // TODO: this should be possible to assign in the problem file
189
190 pn_ = pg_- pcAlpha * pcgn - (1.-pcAlpha)*(pcgw - pcNW1);
191 pw_ = pn_ - pcAlpha * pcnw - (1.-pcAlpha)*pcNW1;
192 }
193 else if (phasePresence == wPhaseOnly || phasePresence == wnPhaseOnly)
194 {
195 pw_ = priVars[pressureIdx];
196
197 // calculate capillary pressures
198 Scalar pcgw = MaterialLaw::pcgw(materialParams, sw_);
199 Scalar pcnw = MaterialLaw::pcnw(materialParams, sw_);
200 Scalar pcgn = MaterialLaw::pcgn(materialParams, sw_ + sn_);
201
202 Scalar pcAlpha = MaterialLaw::pcAlpha(materialParams, sn_);
203 Scalar pcNW1 = 0.0; // TODO: this should be possible to assign in the problem file
204
205 pn_ = pw_ + pcAlpha * pcnw + (1.-pcAlpha)*pcNW1;
206 pg_ = pn_ + pcAlpha * pcgn + (1.-pcAlpha)*(pcgw - pcNW1);
207 }
208 else DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
210
211 fluidState_.setPressure(wPhaseIdx, pw_);
212 fluidState_.setPressure(gPhaseIdx, pg_);
213 fluidState_.setPressure(nPhaseIdx, pn_);
214
215 /* now the temperature */
216 if (phasePresence == wPhaseOnly || phasePresence == wnPhaseOnly || phasePresence == gPhaseOnly)
217 {
218 temp_ = priVars[switch1Idx];
219 }
220 else if (phasePresence == threePhases)
221 {
222 // temp from inverse pwsat and pnsat which have to sum up to pg
223 Scalar temp = FluidSystem::inverseVaporPressureCurve(fluidState_, gPhaseIdx, wCompIdx); // initial guess
224 for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
225 {
226 fluidState_.setTemperature(phaseIdx, temp);
227 }
228 solidState_.setTemperature(temp);
229 Scalar defect = pg_ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx)
230 - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);
231
232 using std::abs;
233 while(abs(defect) > 0.01) // simply a small number chosen ...
234 {
235 Scalar deltaT = 1.e-8 * temp;
236 for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
237 {
238 fluidState_.setTemperature(phaseIdx, temp+deltaT);
239 }
240 solidState_.setTemperature(temp+deltaT);
241 Scalar fUp = pg_ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx)
242 - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);
243
244 for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
245 {
246 fluidState_.setTemperature(phaseIdx, temp-deltaT);
247 }
248 solidState_.setTemperature(temp-deltaT);
249 Scalar fDown = pg_ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx)
250 - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);
251
252 temp = temp - defect * 2. * deltaT / (fUp - fDown);
253
254 for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
255 {
256 fluidState_.setTemperature(phaseIdx, temp);
257 }
258 solidState_.setTemperature(temp);
259 defect = pg_ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx)
260 - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);
261 }
262 temp_ = temp;
263 }
264 else if (phasePresence == wgPhaseOnly)
265 {
266 // temp from inverse pwsat
267 temp_ = FluidSystem::inverseVaporPressureCurve(fluidState_, gPhaseIdx, wCompIdx);
268 }
269 else if (phasePresence == gnPhaseOnly)
270 {
271 // temp from inverse pnsat
272 temp_ = FluidSystem::inverseVaporPressureCurve(fluidState_, gPhaseIdx, nCompIdx);
273 }
274 else DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
276
277 for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
278 {
279 fluidState_.setTemperature(phaseIdx, temp_);
280 }
281 solidState_.setTemperature(temp_);
282
283 // now comes the tricky part: calculate phase composition
284 if (phasePresence == threePhases) {
285
286 // all phases are present, phase compositions are a
287 // result of the the gas <-> liquid equilibrium.
288 Scalar partPressH2O = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx);
289 Scalar partPressNAPL = pg_ - partPressH2O;
290
291 Scalar xgn = partPressNAPL/pg_;
292 Scalar xgw = partPressH2O/pg_;
293
294 // Henry
295 Scalar xwn = partPressNAPL / FluidSystem::henryCoefficient(fluidState_, wPhaseIdx,nCompIdx);
296 Scalar xww = 1.-xwn;
297
298 // Not yet filled with real numbers for the NAPL phase
299 Scalar xnw = partPressH2O / FluidSystem::henryCoefficient(fluidState_, nPhaseIdx,wCompIdx);
300 Scalar xnn = 1.-xnw;
301
302 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
303 fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
304 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
305 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
306 fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
307 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
308
309 Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
310 Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
311 Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
312 Scalar rhoWMolar = FluidSystem::molarDensity(fluidState_, wPhaseIdx);
313 Scalar rhoGMolar = FluidSystem::molarDensity(fluidState_, gPhaseIdx);
314 Scalar rhoNMolar = FluidSystem::molarDensity(fluidState_, nPhaseIdx);
315
316 fluidState_.setDensity(wPhaseIdx, rhoW);
317 fluidState_.setDensity(gPhaseIdx, rhoG);
318 fluidState_.setDensity(nPhaseIdx, rhoN);
319 fluidState_.setMolarDensity(wPhaseIdx, rhoWMolar);
320 fluidState_.setMolarDensity(gPhaseIdx, rhoGMolar);
321 fluidState_.setMolarDensity(nPhaseIdx, rhoNMolar);
322 }
323 else if (phasePresence == wPhaseOnly) {
324 // only the water phase is present, water phase composition is
325 // stored explicitly.
326
327 // extract mole fractions in the water phase
328 Scalar xwn = priVars[switch2Idx];
329 Scalar xww = 1 - xwn;
330
331 // write water mole fractions in the fluid state
332 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
333 fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
334
335 // note that the gas phase is actually not existing!
336 // thus, this is used as phase switch criterion
337 Scalar xgn = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx) / pg_;
338 Scalar xgw = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx) / pg_;
339
340
341 // note that the NAPL phase is actually not existing!
342 // thus, this is used as phase switch criterion
343 // maybe solubility would be better than this approach via Henry
344 Scalar xnn = xwn * FluidSystem::henryCoefficient(fluidState_, wPhaseIdx,nCompIdx) / (xgn * pg_);
345 Scalar xnw = xgw*pg_ / FluidSystem::henryCoefficient(fluidState_, nPhaseIdx,wCompIdx);
346
347 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
348 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
349 fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
350 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
351
352 Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
353 Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
354 Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
355 Scalar rhoWMolar = FluidSystem::molarDensity(fluidState_, wPhaseIdx);
356 Scalar rhoGMolar = FluidSystem::molarDensity(fluidState_, gPhaseIdx);
357 Scalar rhoNMolar = FluidSystem::molarDensity(fluidState_, nPhaseIdx);
358
359 fluidState_.setDensity(wPhaseIdx, rhoW);
360 fluidState_.setDensity(gPhaseIdx, rhoG);
361 fluidState_.setDensity(nPhaseIdx, rhoN);
362 fluidState_.setMolarDensity(wPhaseIdx, rhoWMolar);
363 fluidState_.setMolarDensity(gPhaseIdx, rhoGMolar);
364 fluidState_.setMolarDensity(nPhaseIdx, rhoNMolar);
365 }
366 else if (phasePresence == gnPhaseOnly) {
367
368 // only gas and NAPL phases are present
369
370 Scalar xnw = priVars[switch2Idx];
371 Scalar xnn = 1.-xnw;
372 Scalar xgn = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx) / pg_;
373 Scalar xgw = 1.-xgn;
374
375 // note that the water phase is actually not present
376 // the values are used as switching criteria
377 Scalar xww = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx) / pg_;
378
379 // write mole fractions in the fluid state
380 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
381 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
382 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
383 fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
384 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
385
386 Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
387 Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
388 Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
389 Scalar rhoWMolar = FluidSystem::molarDensity(fluidState_, wPhaseIdx);
390 Scalar rhoGMolar = FluidSystem::molarDensity(fluidState_, gPhaseIdx);
391 Scalar rhoNMolar = FluidSystem::molarDensity(fluidState_, nPhaseIdx);
392
393 fluidState_.setDensity(wPhaseIdx, rhoW);
394 fluidState_.setDensity(gPhaseIdx, rhoG);
395 fluidState_.setDensity(nPhaseIdx, rhoN);
396 fluidState_.setMolarDensity(wPhaseIdx, rhoWMolar);
397 fluidState_.setMolarDensity(gPhaseIdx, rhoGMolar);
398 fluidState_.setMolarDensity(nPhaseIdx, rhoNMolar);
399
400 }
401 else if (phasePresence == wnPhaseOnly) {
402 // water and NAPL are present, phase compositions are a
403 // mole fractions of non-existing gas phase are used as switching criteria
404 Scalar partPressH2O = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx);
405 Scalar partPressNAPL = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);
406
407 Scalar xgn = partPressNAPL/pg_;
408 Scalar xgw = partPressH2O/pg_;
409
410 // Henry
411 Scalar xwn = partPressNAPL / FluidSystem::henryCoefficient(fluidState_, wPhaseIdx,nCompIdx);
412 Scalar xww = 1.-xwn;
413
414 Scalar xnw = partPressH2O / FluidSystem::henryCoefficient(fluidState_, nPhaseIdx,wCompIdx);
415 Scalar xnn = 1.-xnw;
416
417 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
418 fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
419 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
420 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
421 fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
422 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
423
424 Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
425 Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
426 Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
427 Scalar rhoWMolar = FluidSystem::molarDensity(fluidState_, wPhaseIdx);
428 Scalar rhoGMolar = FluidSystem::molarDensity(fluidState_, gPhaseIdx);
429 Scalar rhoNMolar = FluidSystem::molarDensity(fluidState_, nPhaseIdx);
430
431 fluidState_.setDensity(wPhaseIdx, rhoW);
432 fluidState_.setDensity(gPhaseIdx, rhoG);
433 fluidState_.setDensity(nPhaseIdx, rhoN);
434 fluidState_.setMolarDensity(wPhaseIdx, rhoWMolar);
435 fluidState_.setMolarDensity(gPhaseIdx, rhoGMolar);
436 fluidState_.setMolarDensity(nPhaseIdx, rhoNMolar);
437 }
438 else if (phasePresence == gPhaseOnly) {
439 // only the gas phase is present, gas phase composition is
440 // stored explicitly here below.
441
442 const Scalar xgn = priVars[switch2Idx];
443 Scalar xgw = 1 - xgn;
444
445 // write mole fractions in the fluid state
446 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
447 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
448
449 // note that the water and NAPL phase is actually not present
450 // the values are used as switching criteria
451 Scalar xww = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx) / pg_;
452 Scalar xnn = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx) / pg_;
453
454 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
455 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
456
457 Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
458 Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
459 Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
460 Scalar rhoWMolar = FluidSystem::molarDensity(fluidState_, wPhaseIdx);
461 Scalar rhoGMolar = FluidSystem::molarDensity(fluidState_, gPhaseIdx);
462 Scalar rhoNMolar = FluidSystem::molarDensity(fluidState_, nPhaseIdx);
463
464 fluidState_.setDensity(wPhaseIdx, rhoW);
465 fluidState_.setDensity(gPhaseIdx, rhoG);
466 fluidState_.setDensity(nPhaseIdx, rhoN);
467 fluidState_.setMolarDensity(wPhaseIdx, rhoWMolar);
468 fluidState_.setMolarDensity(gPhaseIdx, rhoGMolar);
469 fluidState_.setMolarDensity(nPhaseIdx, rhoNMolar);
470 }
471 else if (phasePresence == wgPhaseOnly) {
472 // only water and gas phases are present
473 const Scalar xgn = priVars[switch2Idx];
474 Scalar xgw = 1 - xgn;
475
476 // write mole fractions in the fluid state
477 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
478 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
479
480
481 Scalar xwn = xgn*pg_/FluidSystem::henryCoefficient(fluidState_, wPhaseIdx,nCompIdx);
482 Scalar xww = 1.-xwn;
483
484 // write mole fractions in the fluid state
485 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
486 fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
487
488 // note that the NAPL phase is actually not existing!
489 // thus, this is used as phase switch criterion
490 Scalar xnn = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx) / pg_;
491
492 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
493
494 Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
495 Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
496 Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
497 Scalar rhoWMolar = FluidSystem::molarDensity(fluidState_, wPhaseIdx);
498 Scalar rhoGMolar = FluidSystem::molarDensity(fluidState_, gPhaseIdx);
499 Scalar rhoNMolar = FluidSystem::molarDensity(fluidState_, nPhaseIdx);
500
501 fluidState_.setDensity(wPhaseIdx, rhoW);
502 fluidState_.setDensity(gPhaseIdx, rhoG);
503 fluidState_.setDensity(nPhaseIdx, rhoN);
504 fluidState_.setMolarDensity(wPhaseIdx, rhoWMolar);
505 fluidState_.setMolarDensity(gPhaseIdx, rhoGMolar);
506 fluidState_.setMolarDensity(nPhaseIdx, rhoNMolar);
507 }
508 else
509 assert(false); // unhandled phase state
510 } // end of if(!UseSimpleModel), i.e. the more complex version with six phase states
511 else // use the simpler model with only two phase states
512 {
513 /* first the saturations */
514 if (phasePresence == threePhases)
515 {
516 sw_ = priVars[switch1Idx];
517 sn_ = priVars[switch2Idx];
518 sg_ = 1. - sw_ - sn_;
519 }
520 else if (phasePresence == wnPhaseOnly)
521 {
522 sn_ = priVars[switch2Idx];
523 sw_ = 1. - sn_;
524 sg_ = 0.;
525 }
526 else DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
528
529 fluidState_.setSaturation(wPhaseIdx, sw_);
530 fluidState_.setSaturation(gPhaseIdx, sg_);
531 fluidState_.setSaturation(nPhaseIdx, sn_);
532
533 /* now the pressures */
534 if (phasePresence == threePhases)
535 {
536 pg_ = priVars[pressureIdx];
537
538 // calculate capillary pressures
539 Scalar pcgw = MaterialLaw::pcgw(materialParams, sw_);
540 Scalar pcnw = MaterialLaw::pcnw(materialParams, sw_);
541 Scalar pcgn = MaterialLaw::pcgn(materialParams, sw_ + sn_);
542
543 Scalar pcAlpha = MaterialLaw::pcAlpha(materialParams, sn_);
544 Scalar pcNW1 = 0.0; // TODO: this should be possible to assign in the problem file
545
546 pn_ = pg_- pcAlpha * pcgn - (1.-pcAlpha)*(pcgw - pcNW1);
547 pw_ = pn_ - pcAlpha * pcnw - (1.-pcAlpha)*pcNW1;
548 }
549 else if (phasePresence == wnPhaseOnly)
550 {
551 pw_ = priVars[pressureIdx];
552
553 // calculate capillary pressures
554 Scalar pcgw = MaterialLaw::pcgw(materialParams, sw_);
555 Scalar pcnw = MaterialLaw::pcnw(materialParams, sw_);
556 Scalar pcgn = MaterialLaw::pcgn(materialParams, sw_ + sn_);
557
558 Scalar pcAlpha = MaterialLaw::pcAlpha(materialParams, sn_);
559 Scalar pcNW1 = 0.0; // TODO: this should be possible to assign in the problem file
560
561 pn_ = pw_ + pcAlpha * pcnw + (1.-pcAlpha)*pcNW1;
562 pg_ = pn_ + pcAlpha * pcgn + (1.-pcAlpha)*(pcgw - pcNW1);
563 }
564 else DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
566
567 fluidState_.setPressure(wPhaseIdx, pw_);
568 fluidState_.setPressure(gPhaseIdx, pg_);
569 fluidState_.setPressure(nPhaseIdx, pn_);
570
571 /* now the temperature */
572 if (phasePresence == wnPhaseOnly)
573 {
574 temp_ = priVars[switch1Idx];
575 }
576 else if (phasePresence == threePhases)
577 {
578 if(sn_<=1.e-10) // this threshold values is chosen arbitrarily as a small number
579 {
580 Scalar tempOnlyWater = FluidSystem::inverseVaporPressureCurve(fluidState_, gPhaseIdx, wCompIdx);
581 temp_ = tempOnlyWater;
582 }
583 if(sw_<=1.e-10) // this threshold values is chosen arbitrarily as a small number
584 {
585 Scalar tempOnlyNAPL = FluidSystem::inverseVaporPressureCurve(fluidState_, gPhaseIdx, nCompIdx);
586 temp_ = tempOnlyNAPL;
587 }
588 else
589 {
590 // temp from inverse pwsat and pnsat which have to sum up to pg
591 Scalar tempOnlyNAPL = FluidSystem::inverseVaporPressureCurve(fluidState_, gPhaseIdx, nCompIdx);
592 Scalar tempOnlyWater = FluidSystem::inverseVaporPressureCurve(fluidState_, gPhaseIdx, wCompIdx);
593 for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
594 {
595 fluidState_.setTemperature(phaseIdx, tempOnlyWater);
596 }
597 solidState_.setTemperature(tempOnlyWater);
598 Scalar defect = pg_ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx)
599 - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);
600
601 Scalar temp = tempOnlyWater; // initial guess
602 int counter = 0;
603 using std::abs;
604 while(abs(defect) > 0.01) // simply a small number chosen ...
605 {
606 Scalar deltaT = 1.e-6; // fixed number, but T should always be in the order of a few hundred Kelvin
607 for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
608 {
609 fluidState_.setTemperature(phaseIdx, temp+deltaT);
610 }
611 solidState_.setTemperature(temp+deltaT);
612 Scalar fUp = pg_ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx)
613 - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);
614
615 for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
616 {
617 fluidState_.setTemperature(phaseIdx, temp-deltaT);
618 }
619 solidState_.setTemperature(temp-deltaT);
620 Scalar fDown = pg_ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx)
621 - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);
622
623 temp = temp - defect * 2. * deltaT / (fUp - fDown);
624
625 for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
626 {
627 fluidState_.setTemperature(phaseIdx, temp);
628 }
629 solidState_.setTemperature(temp);
630 defect = pg_ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx)
631 - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);
632 counter +=1;
633 if (counter>10) break;
634 }
635 if ((sw_>1.e-10)&&(sw_<0.01))
636 temp = temp + (sw_ - 1.e-10) * (temp - tempOnlyNAPL) / (0.01 - 1.e-10);
637 if ((sn_>1.e-10)&&(sn_<0.01))
638 temp = temp + (sn_ - 1.e-10) * (temp - tempOnlyWater) / (0.01 - 1.e-10);
639 temp_ = temp;
640 }
641 }
642 else DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
644
645 for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
646 {
647 fluidState_.setTemperature(phaseIdx, temp_);
648 }
649 solidState_.setTemperature(temp_);
650
651 // now comes the tricky part: calculate phase composition
652 if (phasePresence == threePhases) {
653
654 // all phases are present, phase compositions are a
655 // result of the the gas <-> liquid equilibrium.
656 Scalar partPressH2O = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx);
657 Scalar partPressNAPL = pg_ - partPressH2O;
658 // regularized evaporation for small liquid phase saturations
659 // avoids negative saturations of liquid phases
660 if (sw_<0.02) partPressH2O *= sw_/0.02;
661 if (partPressH2O < 0.) partPressH2O = 0;
662 if (sn_<0.02) partPressNAPL *= sn_ / 0.02;
663 if (partPressNAPL < 0.) partPressNAPL = 0;
664
665 Scalar xgn = partPressNAPL/pg_;
666 Scalar xgw = partPressH2O/pg_;
667
668 // Immiscible liquid phases, mole fractions are just dummy values
669 Scalar xwn = 0;
670 Scalar xww = 1.-xwn;
671
672 // Not yet filled with real numbers for the NAPL phase
673 Scalar xnw = 0;
674 Scalar xnn = 1.-xnw;
675
676 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
677 fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
678 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
679 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
680 fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
681 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
682
683 Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
684 Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
685 Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
686 Scalar rhoWMolar = FluidSystem::molarDensity(fluidState_, wPhaseIdx);
687 Scalar rhoGMolar = FluidSystem::molarDensity(fluidState_, gPhaseIdx);
688 Scalar rhoNMolar = FluidSystem::molarDensity(fluidState_, nPhaseIdx);
689
690 fluidState_.setDensity(wPhaseIdx, rhoW);
691 fluidState_.setDensity(gPhaseIdx, rhoG);
692 fluidState_.setDensity(nPhaseIdx, rhoN);
693 fluidState_.setMolarDensity(wPhaseIdx, rhoWMolar);
694 fluidState_.setMolarDensity(gPhaseIdx, rhoGMolar);
695 fluidState_.setMolarDensity(nPhaseIdx, rhoNMolar);
696 }
697 else if (phasePresence == wnPhaseOnly) {
698 // mole fractions of non-existing gas phase are used as switching criteria
699 Scalar partPressH2O = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx);
700 Scalar partPressNAPL = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);
701
702 Scalar xgn = partPressNAPL/pg_;
703 Scalar xgw = partPressH2O/pg_;
704
705 // immiscible liquid phases, mole fractions are just dummy values
706 Scalar xwn = 0;
707 Scalar xww = 1.-xwn;
708
709 Scalar xnw = 0;
710 Scalar xnn = 1.-xnw;
711
712 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
713 fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
714 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
715 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
716 fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
717 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
718
719 Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
720 Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
721 Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
722 Scalar rhoWMolar = FluidSystem::molarDensity(fluidState_, wPhaseIdx);
723 Scalar rhoGMolar = FluidSystem::molarDensity(fluidState_, gPhaseIdx);
724 Scalar rhoNMolar = FluidSystem::molarDensity(fluidState_, nPhaseIdx);
725
726 fluidState_.setDensity(wPhaseIdx, rhoW);
727 fluidState_.setDensity(gPhaseIdx, rhoG);
728 fluidState_.setDensity(nPhaseIdx, rhoN);
729 fluidState_.setMolarDensity(wPhaseIdx, rhoWMolar);
730 fluidState_.setMolarDensity(gPhaseIdx, rhoGMolar);
731 fluidState_.setMolarDensity(nPhaseIdx, rhoNMolar);
732 }
733 else DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
734 }
735
736 for (int phaseIdx = 0; phaseIdx < numPs; ++phaseIdx) {
737 // Mobilities
738 const Scalar mu =
740 phaseIdx);
741 fluidState_.setViscosity(phaseIdx,mu);
742
743 Scalar kr;
744 kr = MaterialLaw::kr(materialParams, phaseIdx,
745 fluidState_.saturation(wPhaseIdx),
746 fluidState_.saturation(nPhaseIdx),
747 fluidState_.saturation(gPhaseIdx));
748 mobility_[phaseIdx] = kr / mu;
749 Valgrind::CheckDefined(mobility_[phaseIdx]);
750 }
751
752 // material dependent parameters for NAPL adsorption
753 bulkDensTimesAdsorpCoeff_ = MaterialLaw::bulkDensTimesAdsorpCoeff(materialParams);
754
755 // porosity
756 updateSolidVolumeFractions(elemSol, problem, element, scv, solidState_, numFluidComps);
757 EnergyVolVars::updateSolidEnergyParams(elemSol, problem, element, scv, solidState_);
758
759 auto getEffectiveDiffusionCoefficient = [&](int phaseIdx, int compIIdx, int compJIdx)
760 {
761 return EffDiffModel::effectiveDiffusionCoefficient(*this, phaseIdx, compIIdx, compJIdx);
762 };
763
764 effectiveDiffCoeff_.update(getEffectiveDiffusionCoefficient);
765
766 // permeability
767 permeability_ = problem.spatialParams().permeability(element, scv, elemSol);
768 Valgrind::CheckDefined(permeability_);
769
770 fluidState_.setTemperature(temp_);
771 // the enthalpies (internal energies are directly calculated in the fluidstate
772 for (int phaseIdx = 0; phaseIdx < numPs; ++phaseIdx)
773 {
774 Scalar h = FluidSystem::enthalpy(fluidState_, phaseIdx);
775 fluidState_.setEnthalpy(phaseIdx, h);
776 }
777
778 EnergyVolVars::updateEffectiveThermalConductivity();
779 }
780
784 const FluidState &fluidState() const
785 { return fluidState_; }
786
790 const SolidState &solidState() const
791 { return solidState_; }
792
798 Scalar averageMolarMass(int phaseIdx) const
799 { return fluidState_.averageMolarMass(phaseIdx); }
800
807 Scalar saturation(const int phaseIdx) const
808 { return fluidState_.saturation(phaseIdx); }
809
817 Scalar massFraction(const int phaseIdx, const int compIdx) const
818 { return fluidState_.massFraction(phaseIdx, compIdx); }
819
827 Scalar moleFraction(const int phaseIdx, const int compIdx) const
828 { return fluidState_.moleFraction(phaseIdx, compIdx); }
829
836 Scalar density(const int phaseIdx) const
837 { return fluidState_.density(phaseIdx); }
838
845 Scalar molarDensity(const int phaseIdx) const
846 { return fluidState_.molarDensity(phaseIdx); }
847
854 Scalar pressure(const int phaseIdx) const
855 { return fluidState_.pressure(phaseIdx); }
856
864 Scalar temperature() const
865 { return fluidState_.temperature(/*phaseIdx=*/0); }
866
873 Scalar mobility(const int phaseIdx) const
874 { return mobility_[phaseIdx]; }
875
876 Scalar viscosity(const int phaseIdx) const
877 { return fluidState_.viscosity(phaseIdx); }
878
882 Scalar capillaryPressure() const
883 { return fluidState_.capillaryPressure(); }
884
888 Scalar porosity() const
889 { return solidState_.porosity(); }
890
894 Scalar permeability() const
895 { return permeability_; }
896
900 [[deprecated("Will be removed after release 3.2. Use diffusionCoefficient(phaseIdx, compIIdx, compJIdx)!")]]
901 Scalar diffusionCoefficient(int phaseIdx, int compIdx) const
902 { return diffusionCoefficient(phaseIdx, FluidSystem::getMainComponent(phaseIdx), compIdx); }
903
904 /*
905 * \brief Returns the binary diffusion coefficients for a phase in \f$[m^2/s]\f$.
906 */
907 Scalar diffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const
908 {
909 if (phaseIdx != nPhaseIdx)
910 return FluidSystem::diffusionCoefficient(fluidState_, phaseIdx);
911 else
912 return 1.e-10;
913 }
914
918 Scalar effectiveDiffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const
919 { return effectiveDiffCoeff_(phaseIdx, compIIdx, compJIdx); }
920
925 { return bulkDensTimesAdsorpCoeff_; }
926
933 Scalar internalEnergy(int phaseIdx) const
934 { return fluidState_.internalEnergy(phaseIdx); };
935
942 Scalar enthalpy(int phaseIdx) const
943 { return fluidState_.enthalpy(phaseIdx); };
944
945protected:
948
949private:
950 Scalar sw_, sg_, sn_, pg_, pw_, pn_, temp_;
951
952 Scalar moleFrac_[numPs][numFluidComps];
953 Scalar massFrac_[numPs][numFluidComps];
954
955 Scalar permeability_;
956 Scalar mobility_[numPs];
957 Scalar bulkDensTimesAdsorpCoeff_;
958
960 DiffusionCoefficients effectiveDiffCoeff_;
961
962};
963} // end namespace Dumux
964
965#endif
Some exceptions thrown in DuMux
Define some often used mathematical functions.
Some templates to wrap the valgrind macros.
Represents all relevant thermodynamic quantities of a multi-phase, multi-component fluid system assum...
Update the solid volume fractions (inert and reacitve) and set them in the solidstate.
A central place for various physical constants occuring in some equations.
void updateSolidVolumeFractions(const ElemSol &elemSol, const Problem &problem, const Element &element, const Scv &scv, SolidState &solidState, const int solidVolFracOffset)
update the solid volume fractions (inert and reacitve) and set them in the solidstate
Definition: updatesolidvolumefractions.hh:36
bool CheckDefined(const T &value)
Make valgrind complain if the object occupied by an object is undefined.
Definition: valgrind.hh:72
Definition: adapt.hh:29
std::string phasePresence() noexcept
I/O name of phase presence.
Definition: name.hh:147
std::string viscosity(int phaseIdx) noexcept
I/O name of viscosity for multiphase systems.
Definition: name.hh:74
std::string molarDensity(int phaseIdx) noexcept
I/O name of molar density for multiphase systems.
Definition: name.hh:83
std::string density(int phaseIdx) noexcept
I/O name of density for multiphase systems.
Definition: name.hh:65
The primary variable switch controlling the phase presence state variable.
Definition: 3pwateroil/primaryvariableswitch.hh:39
Contains the quantities which are are constant within a finite volume in the three-phase,...
Definition: porousmediumflow/3pwateroil/volumevariables.hh:56
const FluidState & fluidState() const
Returns the phase state for the control-volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:784
Scalar diffusionCoefficient(int phaseIdx, int compIdx) const
Returns the diffusion coefficient.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:901
SolidState solidState_
Definition: porousmediumflow/3pwateroil/volumevariables.hh:947
Scalar diffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const
Definition: porousmediumflow/3pwateroil/volumevariables.hh:907
Scalar molarDensity(const int phaseIdx) const
Returns the molar density of a given phase within the control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:845
Scalar massFraction(const int phaseIdx, const int compIdx) const
Returns the mass fraction of a given component in a given phase within the control volume in .
Definition: porousmediumflow/3pwateroil/volumevariables.hh:817
typename Traits::FluidState FluidState
The type of the object returned by the fluidState() method.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:94
Scalar permeability() const
Returns the permeability within the control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:894
typename Traits::FluidSystem FluidSystem
The type of the fluid system.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:96
Scalar viscosity(const int phaseIdx) const
Definition: porousmediumflow/3pwateroil/volumevariables.hh:876
Scalar bulkDensTimesAdsorpCoeff() const
Returns the adsorption information.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:924
typename Traits::SolidState SolidState
Export type of solid state.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:100
Scalar porosity() const
Returns the average porosity within the control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:888
Scalar enthalpy(int phaseIdx) const
Returns the total enthalpy of a phase in the sub-control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:942
Scalar mobility(const int phaseIdx) const
Returns the effective mobility of a given phase within the control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:873
Scalar capillaryPressure() const
Returns the effective capillary pressure within the control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:882
Scalar moleFraction(const int phaseIdx, const int compIdx) const
Returns the mole fraction of a given component in a given phase within the control volume in .
Definition: porousmediumflow/3pwateroil/volumevariables.hh:827
Scalar saturation(const int phaseIdx) const
Returns the effective saturation of a given phase within the control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:807
Scalar density(const int phaseIdx) const
Returns the mass density of a given phase within the control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:836
typename Traits::SolidSystem SolidSystem
Export type of solid system.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:102
void update(const ElemSol &elemSol, const Problem &problem, const Element &element, const Scv &scv)
Updates all quantities for a given control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:118
Scalar effectiveDiffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const
Returns the effective diffusion coefficients for a phase in .
Definition: porousmediumflow/3pwateroil/volumevariables.hh:918
static constexpr bool onlyGasPhaseCanDisappear()
State if only the gas phase is allowed to disappear.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:106
Scalar temperature() const
Returns temperature inside the sub-control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:864
Scalar internalEnergy(int phaseIdx) const
Returns the total internal energy of a phase in the sub-control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:933
typename ModelTraits::Indices Indices
Export the indices.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:98
const SolidState & solidState() const
Returns the phase state for the control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:790
Scalar averageMolarMass(int phaseIdx) const
Returns the average molar mass of the fluid phase.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:798
FluidState fluidState_
Definition: porousmediumflow/3pwateroil/volumevariables.hh:943
Scalar pressure(const int phaseIdx) const
Returns the effective pressure of a given phase within the control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:854
Definition: porousmediumflow/nonisothermal/volumevariables.hh:75
The isothermal base class.
Definition: porousmediumflow/volumevariables.hh:40
static constexpr int numFluidComponents()
Return number of components considered by the model.
Definition: porousmediumflow/volumevariables.hh:52
const PrimaryVariables & priVars() const
Returns the vector of primary variables.
Definition: porousmediumflow/volumevariables.hh:76
static constexpr int numFluidPhases()
Return number of phases considered by the model.
Definition: porousmediumflow/volumevariables.hh:50
void update(const ElemSol &elemSol, const Problem &problem, const Element &element, const Scv &scv)
Updates all quantities for a given control volume.
Definition: porousmediumflow/volumevariables.hh:64
Base class for the model specific class which provides access to all volume averaged quantities.
Base class for the model specific class which provides access to all volume averaged quantities.
The primary variable switch for the extended Richards model.