3.3.0
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
immiscibleflash.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 *****************************************************************************/
25#ifndef DUMUX_IMMISCIBLE_FLASH_HH
26#define DUMUX_IMMISCIBLE_FLASH_HH
27
28#include <dune/common/fvector.hh>
29#include <dune/common/fmatrix.hh>
30
33
34namespace Dumux {
35
63template <class Scalar, class FluidSystem>
65{
66 static constexpr int numPhases = FluidSystem::numPhases;
67 static constexpr int numComponents = FluidSystem::numComponents;
68 static_assert(numPhases == numComponents,
69 "Immiscibility assumes that the number of phases"
70 " is equal to the number of components");
71
72 using ParameterCache = typename FluidSystem::ParameterCache;
73
74 static constexpr int numEq = numPhases;
75
76 using Matrix = Dune::FieldMatrix<Scalar, numEq, numEq>;
77 using Vector = Dune::FieldVector<Scalar, numEq>;
78
79public:
80 using ComponentVector = Dune::FieldVector<Scalar, numComponents>;
81
86 explicit ImmiscibleFlash(int wettingPhaseIdx = 0)
87 : wettingPhaseIdx_(wettingPhaseIdx)
88 {}
89
96 template <class FluidState>
97 void guessInitial(FluidState &fluidState,
98 ParameterCache &paramCache,
99 const ComponentVector &globalMolarities)
100 {
101 // the sum of all molarities
102 Scalar sumMoles = 0;
103 for (int compIdx = 0; compIdx < numComponents; ++compIdx)
104 sumMoles += globalMolarities[compIdx];
105
106 for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
107 // pressure. use atmospheric pressure as initial guess
108 fluidState.setPressure(phaseIdx, 2e5);
109
110 // saturation. assume all fluids to be equally distributed
111 fluidState.setSaturation(phaseIdx, 1.0/numPhases);
112 }
113 }
114
125 template <class MaterialLaw, class FluidState>
126 void solve(FluidState& fluidState,
127 ParameterCache& paramCache,
128 const MaterialLaw& material,
129 const ComponentVector& globalMolarities)
130 {
132 // Check if all fluid phases are incompressible
134 bool allIncompressible = true;
135 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
136 if (FluidSystem::isCompressible(phaseIdx)) {
137 allIncompressible = false;
138 break;
139 }
140 }
141
142 if (allIncompressible) {
143 paramCache.updateAll(fluidState);
144 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
145 Scalar rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
146 Scalar rhoMolar = FluidSystem::molarDensity(fluidState, paramCache, phaseIdx);
147 fluidState.setDensity(phaseIdx, rho);
148 fluidState.setMolarDensity(phaseIdx, rhoMolar);
149
150 Scalar saturation =
151 globalMolarities[/*compIdx=*/phaseIdx]
152 / fluidState.molarDensity(phaseIdx);
153 fluidState.setSaturation(phaseIdx, saturation);
154 }
155 }
156
157// /////////////////////////
158 // Newton method
160
161 // Jacobian matrix
162 Matrix J;
163 // solution, i.e. phase composition
164 Vector deltaX;
165 // right hand side
166 Vector b;
167
168 completeFluidState_(fluidState, paramCache, material);
169
170 const int nMax = 50; // <- maximum number of newton iterations
171 for (int nIdx = 0; nIdx < nMax; ++nIdx) {
172 // calculate Jacobian matrix and right hand side
173 linearize_(J, b, fluidState, paramCache, material, globalMolarities);
174
175 // Solve J*x = b
176 deltaX = 0;
177
178 try { J.solve(deltaX, b); }
179 catch (Dune::FMatrixError& e)
180 {
181 /*
182 std::cout << "error: " << e << "\n";
183 std::cout << "b: " << b << "\n";
184 */
185
186 throw NumericalProblem(e.what());
187 }
188
189 /*
190 printFluidState_(fluidState);
191 std::cout << "J:\n";
192 for (int i = 0; i < numEq; ++i) {
193 for (int j = 0; j < numEq; ++j) {
194 std::ostringstream os;
195 os << J[i][j];
196
197 std::string s(os.str());
198 do {
199 s += " ";
200 } while (s.size() < 20);
201 std::cout << s;
202 }
203 std::cout << "\n";
204 }
205 std::cout << "residual: " << b << "\n";
206 std::cout << "deltaX: " << deltaX << "\n";
207 std::cout << "---------------\n";
208 */
209
210 // update the fluid quantities.
211 Scalar relError = update_(fluidState, paramCache, material, deltaX);
212
213 if (relError < 1e-9)
214 return;
215 }
216
217 /*
218 printFluidState_(fluidState);
219 std::cout << "globalMolarities: ";
220 for (int compIdx = 0; compIdx < numComponents; ++ compIdx)
221 std::cout << globalMolarities[compIdx] << " ";
222 std::cout << "\n";
223 */
224
225 DUNE_THROW(NumericalProblem,
226 "Flash calculation failed."
227 " {c_alpha^kappa} = {" << globalMolarities << "}, T = "
228 << fluidState.temperature(/*phaseIdx=*/0));
229 }
230
231
232protected:
233 template <class FluidState>
234 void printFluidState_(const FluidState &fs)
235 {
236 std::cout << "saturations: ";
237 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
238 std::cout << fs.saturation(phaseIdx) << " ";
239 std::cout << "\n";
240
241 std::cout << "pressures: ";
242 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
243 std::cout << fs.pressure(phaseIdx) << " ";
244 std::cout << "\n";
245
246 std::cout << "densities: ";
247 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
248 std::cout << fs.density(phaseIdx) << " ";
249 std::cout << "\n";
250
251 std::cout << "global component molarities: ";
252 for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
253 Scalar sum = 0;
254 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
255 sum += fs.saturation(phaseIdx)*fs.molarity(phaseIdx, compIdx);
256 }
257 std::cout << sum << " ";
258 }
259 std::cout << "\n";
260 }
261
262 template <class MaterialLaw, class FluidState>
263 void linearize_(Matrix &J,
264 Vector &b,
265 FluidState &fluidState,
266 ParameterCache &paramCache,
267 const MaterialLaw& material,
268 const ComponentVector &globalMolarities)
269 {
270 FluidState origFluidState(fluidState);
271 ParameterCache origParamCache(paramCache);
272
273 Vector tmp;
274
275 // reset jacobian
276 J = 0;
277
278 calculateDefect_(b, fluidState, fluidState, globalMolarities);
279
280 // assemble jacobian matrix
281 for (int pvIdx = 0; pvIdx < numEq; ++ pvIdx) {
283 // approximately calculate partial derivatives of the
284 // fugacity defect of all components in regard to the mole
285 // fraction of the i-th component. This is done via
286 // forward differences
287
288 // deviate the mole fraction of the i-th component
289 Scalar x_i = getQuantity_(fluidState, pvIdx);
290 const Scalar eps = 1e-10/quantityWeight_(fluidState, pvIdx);
291 setQuantity_<MaterialLaw>(fluidState, paramCache, material, pvIdx, x_i + eps);
292
293 // compute derivative of the defect
294 calculateDefect_(tmp, origFluidState, fluidState, globalMolarities);
295 tmp -= b;
296 tmp /= eps;
297 // store derivative in jacobian matrix
298 for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
299 J[eqIdx][pvIdx] = tmp[eqIdx];
300
301 // fluid state and parameter cache to their original values
302 fluidState = origFluidState;
303 paramCache = origParamCache;
304
305 // end forward differences
307 }
308 }
309
310 template <class FluidState>
311 void calculateDefect_(Vector &b,
312 const FluidState &fluidStateEval,
313 const FluidState &fluidState,
314 const ComponentVector &globalMolarities)
315 {
316 // global molarities are given
317 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
318 b[phaseIdx] =
319 fluidState.saturation(phaseIdx)
320 * fluidState.molarity(phaseIdx, /*compIdx=*/phaseIdx);
321
322 b[phaseIdx] -= globalMolarities[/*compIdx=*/phaseIdx];
323 }
324 }
325
326 template <class MaterialLaw, class FluidState>
327 Scalar update_(FluidState &fluidState,
328 ParameterCache &paramCache,
329 const MaterialLaw& material,
330 const Vector &deltaX)
331 {
332 Scalar relError = 0;
333 for (int pvIdx = 0; pvIdx < numEq; ++ pvIdx) {
334 Scalar tmp = getQuantity_(fluidState, pvIdx);
335 Scalar delta = deltaX[pvIdx];
336 using std::max;
337 using std::abs;
338 using std::min;
339 relError = max(relError, abs(delta)*quantityWeight_(fluidState, pvIdx));
340
341 if (isSaturationIdx_(pvIdx)) {
342 // dampen to at most 20% change in saturation per
343 // iteration
344 delta = min(0.2, max(-0.2, delta));
345 }
346 else if (isPressureIdx_(pvIdx)) {
347 // dampen to at most 30% change in pressure per
348 // iteration
349 delta = min(0.30*fluidState.pressure(0), max(-0.30*fluidState.pressure(0), delta));
350 }
351
352 setQuantityRaw_(fluidState, pvIdx, tmp - delta);
353 }
354
355 completeFluidState_<MaterialLaw>(fluidState, paramCache, material);
356
357 return relError;
358 }
359
360 template <class MaterialLaw, class FluidState>
361 void completeFluidState_(FluidState &fluidState,
362 ParameterCache &paramCache,
363 const MaterialLaw& material)
364 {
365 // calculate the saturation of the last phase as a function of
366 // the other saturations
367 Scalar sumSat = 0.0;
368 for (int phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx)
369 sumSat += fluidState.saturation(phaseIdx);
370
371 if (sumSat > 1.0) {
372 // make sure that the last saturation does not become
373 // negative
374 for (int phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx)
375 {
376 Scalar S = fluidState.saturation(phaseIdx);
377 fluidState.setSaturation(phaseIdx, S/sumSat);
378 }
379 sumSat = 1;
380 }
381
382 // set the last saturation
383 fluidState.setSaturation(/*phaseIdx=*/numPhases - 1, 1.0 - sumSat);
384
385 // update the pressures using the material law (saturations
386 // and first pressure are already set because it is implicitly
387 // solved for.)
388 const auto pc = material.capillaryPressures(fluidState, wettingPhaseIdx_);
389 for (int phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx)
390 fluidState.setPressure(phaseIdx,
391 fluidState.pressure(0)
392 + (pc[phaseIdx] - pc[0]));
393
394 // update the parameter cache
395 paramCache.updateAll(fluidState, /*except=*/ParameterCache::Temperature|ParameterCache::Composition);
396
397 // update all densities
398 for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
399 Scalar rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
400 Scalar rhoMolar = FluidSystem::molarDensity(fluidState, paramCache, phaseIdx);
401 fluidState.setDensity(phaseIdx, rho);
402 fluidState.setMolarDensity(phaseIdx, rhoMolar);
403 }
404 }
405
406 bool isPressureIdx_(int pvIdx)
407 { return pvIdx == 0; }
408
409 bool isSaturationIdx_(int pvIdx)
410 { return 1 <= pvIdx && pvIdx < numPhases; }
411
412 // retrieves a quantity from the fluid state
413 template <class FluidState>
414 Scalar getQuantity_(const FluidState &fs, int pvIdx)
415 {
416 assert(pvIdx < numEq);
417
418 // first pressure
419 if (pvIdx < 1) {
420 int phaseIdx = 0;
421 return fs.pressure(phaseIdx);
422 }
423 // saturations
424 else {
425 assert(pvIdx < numPhases);
426 int phaseIdx = pvIdx - 1;
427 return fs.saturation(phaseIdx);
428 }
429 }
430
431 // set a quantity in the fluid state
432 template <class MaterialLaw, class FluidState>
433 void setQuantity_(FluidState &fs,
434 ParameterCache &paramCache,
435 const MaterialLaw& material,
436 int pvIdx,
437 Scalar value)
438 {
439 assert(0 <= pvIdx && pvIdx < numEq);
440
441 if (pvIdx < 1) {
442 // -> first pressure
443 Scalar delta = value - fs.pressure(0);
444
445 // set all pressures. here we assume that the capillary
446 // pressure does not depend on absolute pressure.
447 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
448 fs.setPressure(phaseIdx, fs.pressure(phaseIdx) + delta);
449 paramCache.updateAllPressures(fs);
450
451 // update all densities
452 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
453 Scalar rho = FluidSystem::density(fs, paramCache, phaseIdx);
454 Scalar rhoMolar = FluidSystem::molarDensity(fs, paramCache, phaseIdx);
455 fs.setDensity(phaseIdx, rho);
456 fs.setMolarDensity(phaseIdx, rhoMolar);
457 }
458 }
459 else {
460 assert(pvIdx < numPhases);
461
462 // -> first M-1 saturations
463 Scalar delta = value - fs.saturation(/*phaseIdx=*/pvIdx - 1);
464 fs.setSaturation(/*phaseIdx=*/pvIdx - 1, value);
465
466 // set last saturation (-> minus the change of the
467 // saturation of the other phases)
468 fs.setSaturation(/*phaseIdx=*/numPhases - 1,
469 fs.saturation(numPhases - 1) - delta);
470
471 // update all fluid pressures using the capillary pressure
472 // law
473 const auto pc = material.capillaryPressures(fs, wettingPhaseIdx_);
474 for (int phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx)
475 fs.setPressure(phaseIdx,
476 fs.pressure(0)
477 + (pc[phaseIdx] - pc[0]));
478 paramCache.updateAllPressures(fs);
479
480 // update all densities
481 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
482 Scalar rho = FluidSystem::density(fs, paramCache, phaseIdx);
483 Scalar rhoMolar = FluidSystem::molarDensity(fs, paramCache, phaseIdx);
484 fs.setDensity(phaseIdx, rho);
485 fs.setMolarDensity(phaseIdx, rhoMolar);
486
487 }
488 }
489 }
490
491 // set a quantity in the fluid state
492 template <class FluidState>
493 void setQuantityRaw_(FluidState &fs, int pvIdx, Scalar value)
494 {
495 assert(pvIdx < numEq);
496
497 // first pressure
498 if (pvIdx < 1) {
499 int phaseIdx = 0;
500 fs.setPressure(phaseIdx, value);
501 }
502 // saturations
503 else {
504 assert(pvIdx < numPhases);
505 int phaseIdx = pvIdx - 1;
506
507 // make sure that the first M-1 saturations does not get
508 // negative
509 using std::max;
510 value = max(0.0, value);
511 fs.setSaturation(phaseIdx, value);
512 }
513 }
514
515 template <class FluidState>
516 Scalar quantityWeight_(const FluidState &fs, int pvIdx)
517 {
518 // first pressure
519 if (pvIdx < 1)
520 return 1.0/1e5;
521 // first M - 1 saturations
522 else {
523 assert(pvIdx < numPhases);
524 return 1.0;
525 }
526 }
527
528private:
529 int wettingPhaseIdx_ = 0;
530
531};
532
533} // end namespace Dumux
534
535#endif
Some exceptions thrown in DuMux
Makes the capillary pressure-saturation relations available under the M-phase API for material laws.
Definition: adapt.hh:29
std::string saturation(int phaseIdx) noexcept
I/O name of saturation for multiphase systems.
Definition: name.hh:43
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
Exception thrown if a fixable numerical problem occurs.
Definition: exceptions.hh:39
Determines the pressures and saturations of all fluid phases given the total mass of all components.
Definition: immiscibleflash.hh:65
void guessInitial(FluidState &fluidState, ParameterCache &paramCache, const ComponentVector &globalMolarities)
Guess initial values for all quantities.
Definition: immiscibleflash.hh:97
void printFluidState_(const FluidState &fs)
Definition: immiscibleflash.hh:234
void solve(FluidState &fluidState, ParameterCache &paramCache, const MaterialLaw &material, const ComponentVector &globalMolarities)
Calculates the chemical equilibrium from the component fugacities in a phase.
Definition: immiscibleflash.hh:126
bool isPressureIdx_(int pvIdx)
Definition: immiscibleflash.hh:406
void setQuantityRaw_(FluidState &fs, int pvIdx, Scalar value)
Definition: immiscibleflash.hh:493
void calculateDefect_(Vector &b, const FluidState &fluidStateEval, const FluidState &fluidState, const ComponentVector &globalMolarities)
Definition: immiscibleflash.hh:311
Scalar getQuantity_(const FluidState &fs, int pvIdx)
Definition: immiscibleflash.hh:414
void completeFluidState_(FluidState &fluidState, ParameterCache &paramCache, const MaterialLaw &material)
Definition: immiscibleflash.hh:361
void linearize_(Matrix &J, Vector &b, FluidState &fluidState, ParameterCache &paramCache, const MaterialLaw &material, const ComponentVector &globalMolarities)
Definition: immiscibleflash.hh:263
Scalar update_(FluidState &fluidState, ParameterCache &paramCache, const MaterialLaw &material, const Vector &deltaX)
Definition: immiscibleflash.hh:327
void setQuantity_(FluidState &fs, ParameterCache &paramCache, const MaterialLaw &material, int pvIdx, Scalar value)
Definition: immiscibleflash.hh:433
bool isSaturationIdx_(int pvIdx)
Definition: immiscibleflash.hh:409
Scalar quantityWeight_(const FluidState &fs, int pvIdx)
Definition: immiscibleflash.hh:516
ImmiscibleFlash(int wettingPhaseIdx=0)
Contruct a new flash.
Definition: immiscibleflash.hh:86
Dune::FieldVector< Scalar, numComponents > ComponentVector
Definition: immiscibleflash.hh:80