From f96dcbf0756f0e3a38903eb3decee0cfb04c17a1 Mon Sep 17 00:00:00 2001 From: hanchuan Date: Mon, 29 Aug 2022 14:30:02 +0200 Subject: [PATCH 1/2] [pnm][2p-test] Fix injected non-wetting mass flux at inlet --- test/porenetwork/2p/main.cc | 1 + test/porenetwork/2p/problem.hh | 26 +++++++++++++++++++++----- 2 files changed, 22 insertions(+), 5 deletions(-) diff --git a/test/porenetwork/2p/main.cc b/test/porenetwork/2p/main.cc index 4a017e4338..db6a8d0e9e 100644 --- a/test/porenetwork/2p/main.cc +++ b/test/porenetwork/2p/main.cc @@ -93,6 +93,7 @@ int main(int argc, char** argv) using SolutionVector = GetPropType; SolutionVector x(leafGridView.size(GridView::dimension)); problem->applyInitialSolution(x); + problem->calculateSumInletVolume(); auto xOld = x; // the grid variables diff --git a/test/porenetwork/2p/problem.hh b/test/porenetwork/2p/problem.hh index f977e8df3d..bed2d8f50e 100644 --- a/test/porenetwork/2p/problem.hh +++ b/test/porenetwork/2p/problem.hh @@ -74,6 +74,7 @@ public: vtpOutputFrequency_ = getParam("Problem.VtpOutputFrequency"); useFixedPressureAndSaturationBoundary_ = getParam("Problem.UseFixedPressureAndSaturationBoundary", false); pc_ = getParam("Problem.CapillaryPressure"); + nonWettingMassFlux_ = getParam("Problem.NonWettingMassFlux", 5e-8); } /*! @@ -158,11 +159,9 @@ public: { PrimaryVariables values(0.0); - // If we do not want to use global phase pressure difference with fixed saturations and pressures, - // we can instead only fix the non-wetting phase pressure and allow the wetting phase saturation to changle freely - // by applying a Nitsche-type boundary condition which tries to minimize the difference between the present pn and the given value + // We fix the mass flux of non-wetting injection at inlet of pore-network if (!useFixedPressureAndSaturationBoundary_ && isInletPore_(scv)) - values[snIdx] = (elemVolVars[scv].pressure(nPhaseIdx) - (1e5 + pc_)) * 1e8; + values[snIdx] = nonWettingMassFlux_/sumInletPoresVolume_; return values; } @@ -190,9 +189,24 @@ public: //! Evaluate the initial invasion state of a pore throat bool initialInvasionState(const Element& element) const { return false; } - // \} + //! Loop over the scv in the domain to calculate the sum volume of inner inlet pores + void calculateSumInletVolume() + { + sumInletPoresVolume_ = 0.0; + for (const auto& element : elements(this->gridGeometry().gridView())) + { + auto fvGeometry = localView(this->gridGeometry()); + fvGeometry.bind(element); + for (const auto& scv : scvs(fvGeometry)) + { + if (isInletPore_(scv)) + sumInletPoresVolume_ += this->gridGeometry().poreVolume(scv.dofIndex())/this->gridGeometry().coordinationNumber(scv.dofIndex()); + } + } + } + private: bool isInletPore_(const SubControlVolume& scv) const @@ -213,6 +227,8 @@ private: int vtpOutputFrequency_; bool useFixedPressureAndSaturationBoundary_; Scalar pc_; + Scalar nonWettingMassFlux_; + Scalar sumInletPoresVolume_; }; } //end namespace Dumux -- GitLab From 0af66ecd5a9facf757c171d4aa578439b94c5a1f Mon Sep 17 00:00:00 2001 From: hanchuan Date: Tue, 30 Aug 2022 09:27:17 +0200 Subject: [PATCH 2/2] [pnm-regularization][2p-test] Total mass can be distributed by the volume ratio or transmissibility ratio --- test/porenetwork/2p/params.input | 1 + test/porenetwork/2p/problem.hh | 21 +++++++++++++++++---- 2 files changed, 18 insertions(+), 4 deletions(-) diff --git a/test/porenetwork/2p/params.input b/test/porenetwork/2p/params.input index c406b2fac9..4f994d96b7 100644 --- a/test/porenetwork/2p/params.input +++ b/test/porenetwork/2p/params.input @@ -13,6 +13,7 @@ VtpOutputFrequency = -1 # Write every n-th time step. 0 only writes a file if an EnableGravity = false CapillaryPressure = 5000 UseFixedPressureAndSaturationBoundary = false +DistributeByVolume = false [Vtk] AddVelocity = 1 diff --git a/test/porenetwork/2p/problem.hh b/test/porenetwork/2p/problem.hh index bed2d8f50e..097121f9e0 100644 --- a/test/porenetwork/2p/problem.hh +++ b/test/porenetwork/2p/problem.hh @@ -73,6 +73,7 @@ public: { vtpOutputFrequency_ = getParam("Problem.VtpOutputFrequency"); useFixedPressureAndSaturationBoundary_ = getParam("Problem.UseFixedPressureAndSaturationBoundary", false); + distributeByVolume_ = getParam("Problem.DistributeByVolume", true); pc_ = getParam("Problem.CapillaryPressure"); nonWettingMassFlux_ = getParam("Problem.NonWettingMassFlux", 5e-8); } @@ -160,10 +161,16 @@ public: PrimaryVariables values(0.0); // We fix the mass flux of non-wetting injection at inlet of pore-network + // The total inlet mass flux is distributed according to the ratio + // of pore volume at inlet to the total volumess if (!useFixedPressureAndSaturationBoundary_ && isInletPore_(scv)) - values[snIdx] = nonWettingMassFlux_/sumInletPoresVolume_; - - return values; + { + if (distributeByVolume_) + values[snIdx] = nonWettingMassFlux_ * ( scv.volume()/sumInletPoresVolume_ ); + else + values[snIdx] = nonWettingMassFlux_ * std::pow((this->gridGeometry().poreInscribedRadius(scv.dofIndex())), 4.0)/sumInletPoresVolume_; + } + return values / scv.volume(); } // \} @@ -202,7 +209,12 @@ public: for (const auto& scv : scvs(fvGeometry)) { if (isInletPore_(scv)) - sumInletPoresVolume_ += this->gridGeometry().poreVolume(scv.dofIndex())/this->gridGeometry().coordinationNumber(scv.dofIndex()); + { + if (distributeByVolume_) + sumInletPoresVolume_ += this->gridGeometry().poreVolume(scv.dofIndex())/this->gridGeometry().coordinationNumber(scv.dofIndex()); + else + sumInletPoresVolume_ += std::pow(this->gridGeometry().poreInscribedRadius(scv.dofIndex()), 4.0); + } } } } @@ -226,6 +238,7 @@ private: int vtpOutputFrequency_; bool useFixedPressureAndSaturationBoundary_; + bool distributeByVolume_; Scalar pc_; Scalar nonWettingMassFlux_; Scalar sumInletPoresVolume_; -- GitLab