quantityreconstruction.hh

Go to the documentation of this file.

// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
//
// SPDX-FileCopyrightText: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
// SPDX-License-Identifier: GPL-3.0-or-later
//
 
#ifndef DUMUX_TWOPVE_QUANTITY_RECONSTRUCTION_HH
#define DUMUX_TWOPVE_QUANTITY_RECONSTRUCTION_HH
 
#include <algorithm>
#include <array>
#include <cmath>
 
#include <dune/common/exceptions.hh>
 
#include <dumux/common/exceptions.hh>
#include <dumux/nonlinear/findscalarroot.hh>
 
namespace Dumux {
 
namespace TwoPVE {

template<typename Scalar>
struct GasPlumeDistancesData
{
    Scalar current;
    Scalar minimum;
};

template<typename Scalar>
struct PhaseDensitiesData
{
    Scalar wetting;
    Scalar nonwetting;
};

template<typename Scalar>
struct PhaseViscositiesData
{
    Scalar wetting;
    Scalar nonwetting;
};

template<typename Scalar>
struct ResidualSaturationsData
{
    Scalar wetting;
    Scalar nonwetting;
};

template<typename Scalar>
struct BrooksCoreyParametersData
{
    Scalar lambda;
    Scalar entryPressure;
};
 
} // end namespace TwoPVE

template<class Scalar, class FluidSystem>
class TwoPVEQuantityReconstruction
{
    static constexpr int wettingPhaseIdx = FluidSystem::phase0Idx;
    static constexpr int nonwettingPhaseIdx = FluidSystem::phase1Idx;
    static constexpr int numPhases = FluidSystem::numPhases;
 
    using GasPlumeDistances = TwoPVE::GasPlumeDistancesData<Scalar>;
    using PhaseDensities = TwoPVE::PhaseDensitiesData<Scalar>;
    using PhaseViscosities = TwoPVE::PhaseViscositiesData<Scalar>;
    using ResidualSaturations = TwoPVE::ResidualSaturationsData<Scalar>;
    using BrooksCoreyParameters = TwoPVE::BrooksCoreyParametersData<Scalar>;
 
public:

    Scalar computeGasPlumeDist(const PhaseDensities& densities,
                               const ResidualSaturations& residualSaturations,
                               Scalar gravityNorm,
                               Scalar domainHeight,
                               Scalar saturationWCoarse,
                               Scalar minGasPlumeDist,
                               const BrooksCoreyParameters& brooksCoreyParameters) const
    {
        const Scalar densityW = densities.wetting;
        const Scalar densityNw = densities.nonwetting;
        const Scalar swr = residualSaturations.wetting;
        const Scalar snr = residualSaturations.nonwetting;
 
        if (gravityNorm == 0.0)
            DUNE_THROW(Dune::InvalidStateException, "The two-phase vertical-equilibrium model requires nonzero gravity because its reconstruction assumes gravity-driven vertical segregation");
 
        if (densityW <= densityNw)
            DUNE_THROW(Dune::InvalidStateException, "The two-phase VE reconstruction requires the wetting phase to be denser than the nonwetting phase. Got rho_w=" << densityW << " and rho_n=" << densityNw);
 
        // residual of the water balance: water content of the column minus the coarse-level water content,
        // monotonically increasing in the gas plume distance
        using std::max; using std::min;
        const auto massConservation = [&](const Scalar gasPlumeDist)
        {
            const Scalar bottomOfMobileGas = max(gasPlumeDist, 0.0);
            const Scalar waterBelowMinimum = max(min(gasPlumeDist, minGasPlumeDist), 0.0);
            const Scalar waterTrappedRegion = (1.0-snr) * max(bottomOfMobileGas - max(minGasPlumeDist, 0.0), 0.0);
            const Scalar waterAbovePlume = integrateSaturationWAbovePlume_(bottomOfMobileGas, domainHeight, gasPlumeDist, densities, residualSaturations, gravityNorm, brooksCoreyParameters);
            return waterBelowMinimum + waterTrappedRegion + waterAbovePlume - saturationWCoarse * domainHeight;
        };
 
        if (massConservation(domainHeight) <= 0.0)
            return domainHeight;
 
        // the water content approaches the residual water content as the gas plume distance decreases below the column
        if (saturationWCoarse <= swr)
            DUNE_THROW(NumericalProblem, "The coarse-level wetting-phase saturation " << saturationWCoarse << " does not exceed the residual saturation " << swr);
 
        // move the lower end of the bracket downwards until the residual is no longer positive
        const Scalar fringeHeight = brooksCoreyParameters.entryPressure/((densityW - densityNw)*gravityNorm);
        Scalar lowerBound = 0.0;
        Scalar residualAtLowerBound = massConservation(lowerBound);
        for (Scalar expansion = fringeHeight; residualAtLowerBound > 0.0; expansion *= 2.0)
        {
            lowerBound -= expansion;
            residualAtLowerBound = massConservation(lowerBound);
        }
 
        if (residualAtLowerBound == 0.0)
            return lowerBound;
 
        return findScalarRootBrent(lowerBound, domainHeight, massConservation);
    }

    Scalar computeCapillaryPressureCoarse(Scalar gasPlumeDist,
                                          const PhaseDensities& densities,
                                          Scalar gravityNorm,
                                          Scalar entryPressureBC) const
    {
        return gravityNorm * gasPlumeDist * (densities.nonwetting - densities.wetting) + entryPressureBC;
    }

    std::array<Scalar, numPhases> reconstructPressure(Scalar gasPlumeDist,
                                                      const PhaseDensities& densities,
                                                      Scalar gravityNorm,
                                                      Scalar heightAboveBottom,
                                                      Scalar pressureWCoarse,
                                                      Scalar entryPressureBC) const
    {
        const Scalar referenceDensityW = densities.wetting;
        const Scalar referenceDensityNw = densities.nonwetting;
 
        // the nonwetting-phase pressure below the gas plume exceeds the wetting-phase pressure by the entry pressure (capillary fringe model)
        std::array<Scalar, numPhases> reconstructedPressures;
        reconstructedPressures[wettingPhaseIdx] = pressureWCoarse - referenceDensityW * gravityNorm * heightAboveBottom;
        if (heightAboveBottom <= gasPlumeDist)
            reconstructedPressures[nonwettingPhaseIdx] = reconstructedPressures[wettingPhaseIdx] + entryPressureBC;
        else
            reconstructedPressures[nonwettingPhaseIdx] = pressureWCoarse - referenceDensityW * gravityNorm * gasPlumeDist - referenceDensityNw * gravityNorm * (heightAboveBottom - gasPlumeDist) + entryPressureBC;
 
        return reconstructedPressures;
    }

    Scalar reconstructSaturation(const GasPlumeDistances& gasPlumeDistances,
                                 const PhaseDensities& densities,
                                 const ResidualSaturations& residualSaturations,
                                 Scalar gravityNorm,
                                 Scalar heightAboveBottom,
                                 Scalar deltaZ,
                                 const BrooksCoreyParameters& brooksCoreyParameters) const
    {
        const Scalar gasPlumeDist = gasPlumeDistances.current;
        const Scalar minGasPlumeDist = gasPlumeDistances.minimum;
 
        const Scalar lowerBound = heightAboveBottom - deltaZ/2.0;
        const Scalar upperBound = heightAboveBottom + deltaZ/2.0;
        const Scalar cellHeight = upperBound - lowerBound;
 
        const Scalar snr = residualSaturations.nonwetting;
        const Scalar saturationWBelowPlume = 1.0;
        const Scalar saturationWTrappedRegion = 1.0 - snr;
        const auto integrateSaturationWAbovePlume = [&](const Scalar lower, const Scalar upper)
        {
            return integrateSaturationWAbovePlume_(lower, upper, gasPlumeDist, densities, residualSaturations, gravityNorm, brooksCoreyParameters);
        };
 
        Scalar saturationIntegral = 0.0;
 
        // differentiate between six cases:
        //     1) fine-level element is completely under the gas plume,
        //     2) fine-level element is completely above the gas plume,
        //     3) fine-level element is completely in trapped gas region
        //     4) fine-level element intersects only historical minimum plume distance
        //     5) fine-level element intersects only the current gas plume distance
        //     6) fine-level element intersects both plume distances
 
        if (upperBound <= minGasPlumeDist) // 1)
            saturationIntegral = saturationWBelowPlume*cellHeight;
        else if (lowerBound >= gasPlumeDist) // 2)
            saturationIntegral = integrateSaturationWAbovePlume(lowerBound, upperBound);
        else if (lowerBound >= minGasPlumeDist && upperBound <= gasPlumeDist) // 3)
            saturationIntegral = saturationWTrappedRegion*cellHeight;
        else if (lowerBound < minGasPlumeDist && upperBound <= gasPlumeDist) // 4)
        {
            const Scalar integralBelowMinimum = saturationWBelowPlume * (minGasPlumeDist - lowerBound);
            const Scalar integralTrappedRegion = saturationWTrappedRegion * (upperBound - minGasPlumeDist);
            saturationIntegral = integralBelowMinimum + integralTrappedRegion;
        }
        else if (lowerBound >= minGasPlumeDist && upperBound > gasPlumeDist) // 5)
        {
            const Scalar integralTrappedRegion = saturationWTrappedRegion * (gasPlumeDist - lowerBound);
            const Scalar integralAbovePlume = integrateSaturationWAbovePlume(gasPlumeDist, upperBound);
            saturationIntegral = integralTrappedRegion + integralAbovePlume;
        }
        else // 6)
        {
            const Scalar integralBelowMinimum = saturationWBelowPlume * (minGasPlumeDist - lowerBound);
            const Scalar integralTrappedRegion = saturationWTrappedRegion * (gasPlumeDist - minGasPlumeDist);
            const Scalar integralAbovePlume = integrateSaturationWAbovePlume(gasPlumeDist, upperBound);
            saturationIntegral = integralBelowMinimum + integralTrappedRegion + integralAbovePlume;
        }
 
        return saturationIntegral/cellHeight;
    }

    Scalar reconstructCapillaryPressure(Scalar gasPlumeDist,
                                        const PhaseDensities& densities,
                                        Scalar gravityNorm,
                                        Scalar heightAboveBottom,
                                        Scalar entryPressureBC) const
    {
        const Scalar referenceDensityW = densities.wetting;
        const Scalar referenceDensityNw = densities.nonwetting;
 
        if (heightAboveBottom <= gasPlumeDist)
            return entryPressureBC;
        else
            return referenceDensityW * gravityNorm * (heightAboveBottom - gasPlumeDist) + entryPressureBC - referenceDensityNw * gravityNorm * (heightAboveBottom - gasPlumeDist);
    }

    std::array<Scalar, numPhases> reconstructMobilitiesFine(const GasPlumeDistances& gasPlumeDistances,
                                                            const PhaseDensities& densities,
                                                            const PhaseViscosities& viscosities,
                                                            Scalar gravityNorm,
                                                            Scalar heightAboveBottom,
                                                            Scalar deltaZ,
                                                            const BrooksCoreyParameters& brooksCoreyParameters) const
    {
        const Scalar gasPlumeDist = gasPlumeDistances.current;
        const Scalar lowerBound = heightAboveBottom - deltaZ/2.0;
        const Scalar upperBound = heightAboveBottom + deltaZ/2.0;
 
        const Scalar relPermWBelowPlume = 1.0;
        const Scalar relPermNwBelowPlume = 0.0;
 
        // Brooks-Corey relative permeabilities in terms of the effective saturation Se = u^(-lambda):
        // krw = Se^(2/lambda + 3) and krn = (1 - Se)^2 (1 - Se^(2/lambda + 1)), expanded into powers of u
        const Scalar lambdaBC = brooksCoreyParameters.lambda;
        const auto integratePower = [&](const Scalar exponent, const Scalar lower, const Scalar upper)
        {
            return integratePowerAbovePlume_(exponent, lower, upper, gasPlumeDist, densities, gravityNorm, brooksCoreyParameters.entryPressure);
        };
        const auto integrateRelPermWAbovePlume = [&](const Scalar lower, const Scalar upper)
        {
            return integratePower(2.0 + 3.0*lambdaBC, lower, upper);
        };
        const auto integrateRelPermNwAbovePlume = [&](const Scalar lower, const Scalar upper)
        {
            const Scalar integral = (upper - lower)
                                    - 2.0*integratePower(lambdaBC, lower, upper)
                                    + integratePower(2.0*lambdaBC, lower, upper)
                                    - integratePower(2.0 + lambdaBC, lower, upper)
                                    + 2.0*integratePower(2.0 + 2.0*lambdaBC, lower, upper)
                                    - integratePower(2.0 + 3.0*lambdaBC, lower, upper);
 
            // close to the gas plume distance, where krn vanishes, the terms cancel up to round-off errors of either sign
            using std::max;
            return max(integral, 0.0);
        };
 
        // the relative permeabilities are piecewise defined, with a kink or jump at the gas plume distance
        const auto integrateRelativePermeability = [&](const auto& integralAbovePlume, const Scalar relPermBelowPlume)
        {
            const Scalar cellHeight = upperBound - lowerBound;
 
            // differentiate between three cases:
            //     1) fine-level element is completely under the gas plume,
            //     2) fine-level element is completely above the gas plume,
            //     3) fine-level element intersects with gas plume
            //     for zp_min<z<zp, there might be snr present and the rel perm for gas is treated as 0, thus the check z<zp is enough
            if (upperBound <= gasPlumeDist)
                return relPermBelowPlume * cellHeight;
            else if (lowerBound >= gasPlumeDist)
                return integralAbovePlume(lowerBound, upperBound);
            else
                return relPermBelowPlume * (gasPlumeDist - lowerBound) + integralAbovePlume(gasPlumeDist, upperBound);
        };
 
        // cell averages of the relative permeabilities divided by the viscosities
        std::array<Scalar, numPhases> mobilities;
        mobilities[wettingPhaseIdx] = integrateRelativePermeability(integrateRelPermWAbovePlume, relPermWBelowPlume)/(upperBound - lowerBound)/viscosities.wetting;
        mobilities[nonwettingPhaseIdx] = integrateRelativePermeability(integrateRelPermNwAbovePlume, relPermNwBelowPlume)/(upperBound - lowerBound)/viscosities.nonwetting;
        return mobilities;
    }
 
private:
    Scalar integratePowerAbovePlume_(Scalar exponent,
                                     Scalar lowerBound,
                                     Scalar upperBound,
                                     Scalar gasPlumeDist,
                                     const PhaseDensities& densities,
                                     Scalar gravityNorm,
                                     Scalar entryPressureBC) const
    {
        using std::expm1; using std::log1p; using std::pow;
        const Scalar lengthScale = entryPressureBC/((densities.wetting - densities.nonwetting)*gravityNorm);
        const Scalar uLower = 1.0 + (lowerBound - gasPlumeDist)/lengthScale;
        const Scalar logRatio = log1p((upperBound - lowerBound)/(lengthScale*uLower));
        const Scalar k = 1.0 - exponent;
        if (k == 0.0)
            return lengthScale*logRatio;
 
        // expm1 keeps the antiderivative accurate for exponents close to one, where it approaches the logarithm
        return lengthScale*pow(uLower, k)*expm1(k*logRatio)/k;
    }

    Scalar integrateSaturationWAbovePlume_(Scalar lowerBound,
                                           Scalar upperBound,
                                           Scalar gasPlumeDist,
                                           const PhaseDensities& densities,
                                           const ResidualSaturations& residualSaturations,
                                           Scalar gravityNorm,
                                           const BrooksCoreyParameters& brooksCoreyParameters) const
    {
        const Scalar swr = residualSaturations.wetting;
        const Scalar snr = residualSaturations.nonwetting;
        const Scalar integralEffectiveSaturation = integratePowerAbovePlume_(brooksCoreyParameters.lambda, lowerBound, upperBound, gasPlumeDist, densities, gravityNorm, brooksCoreyParameters.entryPressure);
        return swr*(upperBound - lowerBound) + (1.0 - swr - snr)*integralEffectiveSaturation;
    }
};
 
} // end namespace Dumux
 
#endif