diff --git a/CHANGELOG.md b/CHANGELOG.md index 6a5481107793d05ed9bfbd52acfaa336ac795286..cb39e723b157d7b1eee5ddf32734392e5e27d589 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -20,6 +20,7 @@ that implement an empty cache (i.e. nothing is cached for the flux variables). - __Pore-network model__: Pore-network model will no longer prevent non-wetting fluid flowing out by default. Throats blocking non-wetting fluid can be specified by setting runtime parameter `InvasionState.BlockNonwettingPhaseAtThroatLabel`. +- __Examples__: Extend the porenetwork_upscaling example to include non-creeping flow simulation in pore network. The example is able now to provide not only upscaled Darcy permeability but also Forchheimer permeability and coefficient (employed in Forchheimer equation). ### Immediate interface changes not allowing/requiring a deprecation period: diff --git a/examples/porenetwork_upscaling/CMakeLists.txt b/examples/porenetwork_upscaling/CMakeLists.txt index 03925f1192758bc69a9d17c42f31789561225f95..da306ea8452d9c3515534f80e5cbbef7fc0c7282 100644 --- a/examples/porenetwork_upscaling/CMakeLists.txt +++ b/examples/porenetwork_upscaling/CMakeLists.txt @@ -1,9 +1,19 @@ dune_symlink_to_source_files(FILES "params.input") -dumux_add_test(NAME example_pnm1p_permeabilityupscaling +dumux_add_test(NAME example_pnm1p_upscaling LABELS porenetwork example SOURCES main.cc CMAKE_GUARD HAVE_UMFPACK dune-foamgrid_FOUND - COMMAND ${CMAKE_CURRENT_BINARY_DIR}/example_pnm1p_permeabilityupscaling - CMD_ARGS -Problem.ReferenceData "2.95910919e-13 2.91015548e-13 2.71752264e-13" - -Problem.TestEpsilon 1e-7) + COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py + CMD_ARGS --script fuzzyData --delimiter " " + --files ${CMAKE_SOURCE_DIR}/test/references/example_porenetwork_upscaling_noncreeping_flow_X-dir_PermeabilityratioVsForchheimerNumber.dat + ${CMAKE_CURRENT_BINARY_DIR}/X-dir-PermeabilityratioVsForchheimerNumber.dat + ${CMAKE_SOURCE_DIR}/test/references/example_porenetwork_upscaling_noncreeping_flow_X-dir_InversePrmeabilityVsInertiaToViscousRatio.dat + ${CMAKE_CURRENT_BINARY_DIR}/X-dir-InversePrmeabilityVsInertiaToViscousRatio.dat + --command "${CMAKE_CURRENT_BINARY_DIR}/example_pnm1p_upscaling params.input -Problem.AssumeCreepingFlow false -Problem.Directions 0") + +dumux_add_test(NAME example_pnm1p_creeping_flow_upscaling + TARGET example_pnm1p_upscaling + CMAKE_GUARD HAVE_UMFPACK dune-foamgrid_FOUND + COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py + CMD_ARGS --command "${CMAKE_CURRENT_BINARY_DIR}/example_pnm1p_upscaling params.input -Problem.AssumeCreepingFlow true") diff --git a/examples/porenetwork_upscaling/README.md b/examples/porenetwork_upscaling/README.md index 23bd4d1dacdf85a78138d13a32a11ad66e2dfc10..92b1b7c60a5f4e8a003cb93fff0aea4f53bb9344 100644 --- a/examples/porenetwork_upscaling/README.md +++ b/examples/porenetwork_upscaling/README.md @@ -1,22 +1,20 @@ -# Determining the Darcy permeability of a pore network +# Determining the upscaled properties of a pore network __In this example, you will learn how to__ -* simulate flow on a pore network by applying a pressure gradient in a given direction -* perform an upscaling in order to determine the intrinsic single-phase permeability $`\mathbf{K}`$ [m$`^2`$] +* simulate creeping/non-creeping flow on a pore network by applying a pressure gradient in a given direction +* perform an upscaling in order to determine the flow properties of the porous medium such as: +the Darcy (intrinsic) single-phase permeability $`\mathbf{K}`$ [m$`^2`$] using the creeping flow simulation, the Forchheimer permeability $`\mathbf{K}`$ [m$`^2`$] and the Forchheimer coefficient $`\mathbf{\beta}`$ [m$`^{-1}`$] using the non-creeping flow simulation. __Result__. -As a result of the simulation of this example, you will get the intrinsic single-phase permeabilities for each spatial direction $`K_{xx}`$, $`K_{yy}`$, $`K_{zz}`$ [m$`^2`$] -as direct output on your terminal. -Additionally, the output contains the values of the auxiliary parameters: cross-sectional area, mass flux integrated over the cross-sectional area, -and the resulting Darcy velocity as in the example below for the x-direction permeability $`K_{xx}`$: +As a result of the creeping flow simulation of this example, you will get the Darcy (intrinsic) single-phase permeabilities for each spatial direction requested by the user $`K_{xx}`$, $`K_{yy}`$, $`K_{zz}`$ [m$`^2`$] as a direct output on your terminal as the following example for the x-direction: ``` - -x-direction: Area = 1.600000e-07 m^2; Massflux = 4.509338e-13 kg/s; v_Darcy = 2.818337e-09 m/s; K = 2.818337e-13 m^2 +X-direction: +-- Darcy (intrinsic) permeability = 3.326e-12 m^2 ``` @@ -30,6 +28,35 @@ Figure 1 shows the pressure distribution within the pore network for the case of +The non-creeping flow simulation additionally gives you the Forchheimer permeability and coefficient for each spatial direction as in the example below for the x-direction: + +``` +X-direction: +-- Darcy (intrinsic) permeability = 3.326e-12 m^2 +-- Forchheimer permeability = 3.196e-12 m^2 +-- Forchheimer coefficient = 8.553e+04 m^-1 +``` + +Furthermore, the ratio of apparent permeability to Darcy permeability is plotted versus the Forchheimer number for the spatial dimensions specified by the user, see Fig. 2. + +
+
+ Permeability ratio vs. Forchheimer number +
Fig.2 - Variation of apparent to Darcy permeability ratio versus Forchheimer number.
+
+
+ +In addition, in Fig. 3, the inverse of the apparent permeability versus $`\varrho v /\mu`$ is plotted to compare the data from the pore-network model simulation and the data obtained by applying the Forchheimer equation using the Forchheimer permeability and coefficient. As evident in the figure, the Forchheimer equation is not able to accurately predict the flow behavior in the low-velocity regime. +
+
+ Inverse of apparent permeability vs. rhovmu-1 +
Fig.3 - Inverse of apparent permeability versus ρv / μ. +
+
+
+ +After building the executable, use the parameter `Problem.AssumeCreepingFlow` in the `params.input` file to select the flow regime to be simulated (i.e. set it to `true` for creeping flow and to `false` for non-creeping flow simulations). Then, run the simulation with `./example_pnm1p_upscaling`. + __Table of contents__. This description is structured as follows: [[_TOC_]] @@ -39,23 +66,35 @@ __Table of contents__. This description is structured as follows: We consider a single-phase problem within a randomly generated pore network of 20x20x20 pores cubed from which some of the pore throats were [deleted randomly](https://doi.org/10.1029/2010WR010180). The inscribed pore body radii follow a truncated log-normal distribution. -To calculate the upscaled permeability, a pressure difference of 4000 Pa is applied sequentially in every direction, while all lateral sides are closed. -The resulting mass flow rate is then used to determine $`\mathbf{K}`$, as described [later](upscalinghelper.md). +To calculate the upscaled properties, $`10`$ pressure differences in the range of $`10`$ to $`10^{10}`$ Pa/m are applied sequentially in directions specified by the user by setting the parameter `Problem.Directions` in the `params.input` file, while all lateral sides are closed. +The resulting mass flow rates are then used to determine upscaled properties as described [later](upscalinghelper.md). ## Mathematical and numerical model -In this example we are using the single-phase pore-network model of DuMux, which considers a Hagen-Poiseuille-type law to relate the volume flow from on pore body to another to discrete pressure drops $`\Delta p = p_i - p_j`$ between the pore bodies. -We require mass conservation at each pore body $`i`$: +In this example, we are using the single-phase pore-network model of DuMux. We require mass conservation at each pore body $`i`$: ```math \sum_j Q_{ij} = 0, ``` -where $`Q_{ij}`$ is the discrete volume flow rate in a throat connecting pore bodies $`i`$ and $`j`$: +where $`Q_{ij}`$ is the discrete volume flow rate in a throat connecting pore bodies $`i`$ and $`j`$. In case of creeping flow, the pore network model considers a Hagen-Poiseuille-type law to relate the volume flow from on pore body to another to discrete pressure drops $`\Delta p = p_i - p_j`$ between the pore bodies using throat conductance, $`g_{ij}`$. + +```math + Q_{ij} = g_{ij} (p_i - p_j), +``` +or + +```math + (p_i - p_j) = Q_{ij} / g_{ij} . +``` + +In the simulation of non-creeping flow, to capture inertial effects in fluid flow, an extension of the Hagen-Poiseuille-type law which includes the expansion and contraction of the pore space when moving from a throat to a pore and vice versa. ```math - Q_{ij} = g_{ij} (p_i - p_j). + (p_i - p_j) = Q_{ij} / g_{ij} + (C_{exp} + C_{cont})Q^2, ``` +where $`C_{exp}`$ and $`C_{cont}`$ are the expansion and the contraction coefficient, respectively. + # Implementation & Post processing In the following, we take a closer look at the source files for this example. diff --git a/examples/porenetwork_upscaling/doc/_intro.md b/examples/porenetwork_upscaling/doc/_intro.md index 64f17fefbd118188c98cbbd845b7e8aed2931a62..12706ab5826d4fafb61bd7feb561f3c5bc3fedf3 100644 --- a/examples/porenetwork_upscaling/doc/_intro.md +++ b/examples/porenetwork_upscaling/doc/_intro.md @@ -1,20 +1,18 @@ -# Determining the Darcy permeability of a pore network +# Determining the upscaled properties of a pore network __In this example, you will learn how to__ -* simulate flow on a pore network by applying a pressure gradient in a given direction -* perform an upscaling in order to determine the intrinsic single-phase permeability $`\mathbf{K}`$ [m$`^2`$] +* simulate creeping/non-creeping flow on a pore network by applying a pressure gradient in a given direction +* perform an upscaling in order to determine the flow properties of the porous medium such as: +the Darcy (intrinsic) single-phase permeability $`\mathbf{K}`$ [m$`^2`$] using the creeping flow simulation, the Forchheimer permeability $`\mathbf{K}`$ [m$`^2`$] and the Forchheimer coefficient $`\mathbf{\beta}`$ [m$`^{-1}`$] using the non-creeping flow simulation. __Result__. -As a result of the simulation of this example, you will get the intrinsic single-phase permeabilities for each spatial direction $`K_{xx}`$, $`K_{yy}`$, $`K_{zz}`$ [m$`^2`$] -as direct output on your terminal. -Additionally, the output contains the values of the auxiliary parameters: cross-sectional area, mass flux integrated over the cross-sectional area, -and the resulting Darcy velocity as in the example below for the x-direction permeability $`K_{xx}`$: +As a result of the creeping flow simulation of this example, you will get the Darcy (intrinsic) single-phase permeabilities for each spatial direction requested by the user $`K_{xx}`$, $`K_{yy}`$, $`K_{zz}`$ [m$`^2`$] as a direct output on your terminal as the following example for the x-direction: ``` - -x-direction: Area = 1.600000e-07 m^2; Massflux = 4.509338e-13 kg/s; v_Darcy = 2.818337e-09 m/s; K = 2.818337e-13 m^2 +X-direction: +-- Darcy (intrinsic) permeability = 3.326e-12 m^2 ``` @@ -28,6 +26,35 @@ Figure 1 shows the pressure distribution within the pore network for the case of +The non-creeping flow simulation additionally gives you the Forchheimer permeability and coefficient for each spatial direction as in the example below for the x-direction: + +``` +X-direction: +-- Darcy (intrinsic) permeability = 3.326e-12 m^2 +-- Forchheimer permeability = 3.196e-12 m^2 +-- Forchheimer coefficient = 8.553e+04 m^-1 +``` + +Furthermore, the ratio of apparent permeability to Darcy permeability is plotted versus the Forchheimer number for the spatial dimensions specified by the user, see Fig. 2. + +
+
+ Permeability ratio vs. Forchheimer number +
Fig.2 - Variation of apparent to Darcy permeability ratio versus Forchheimer number.
+
+
+ +In addition, in Fig. 3, the inverse of the apparent permeability versus $`\varrho v /\mu`$ is plotted to compare the data from the pore-network model simulation and the data obtained by applying the Forchheimer equation using the Forchheimer permeability and coefficient. As evident in the figure, the Forchheimer equation is not able to accurately predict the flow behavior in the low-velocity regime. +
+
+ Inverse of apparent permeability vs. rhovmu-1 +
Fig.3 - Inverse of apparent permeability versus ρv / μ. +
+
+
+ +After building the executable, use the parameter `Problem.AssumeCreepingFlow` in the `params.input` file to select the flow regime to be simulated (i.e. set it to `true` for creeping flow and to `false` for non-creeping flow simulations). Then, run the simulation with `./example_pnm1p_upscaling`. + __Table of contents__. This description is structured as follows: [[_TOC_]] @@ -37,23 +64,35 @@ __Table of contents__. This description is structured as follows: We consider a single-phase problem within a randomly generated pore network of 20x20x20 pores cubed from which some of the pore throats were [deleted randomly](https://doi.org/10.1029/2010WR010180). The inscribed pore body radii follow a truncated log-normal distribution. -To calculate the upscaled permeability, a pressure difference of 4000 Pa is applied sequentially in every direction, while all lateral sides are closed. -The resulting mass flow rate is then used to determine $`\mathbf{K}`$, as described [later](upscalinghelper.md). +To calculate the upscaled properties, $`10`$ pressure differences in the range of $`10`$ to $`10^{10}`$ Pa/m are applied sequentially in directions specified by the user by setting the parameter `Problem.Directions` in the `params.input` file, while all lateral sides are closed. +The resulting mass flow rates are then used to determine upscaled properties as described [later](upscalinghelper.md). ## Mathematical and numerical model -In this example we are using the single-phase pore-network model of DuMux, which considers a Hagen-Poiseuille-type law to relate the volume flow from on pore body to another to discrete pressure drops $`\Delta p = p_i - p_j`$ between the pore bodies. -We require mass conservation at each pore body $`i`$: +In this example, we are using the single-phase pore-network model of DuMux. We require mass conservation at each pore body $`i`$: ```math \sum_j Q_{ij} = 0, ``` -where $`Q_{ij}`$ is the discrete volume flow rate in a throat connecting pore bodies $`i`$ and $`j`$: +where $`Q_{ij}`$ is the discrete volume flow rate in a throat connecting pore bodies $`i`$ and $`j`$. In case of creeping flow, the pore network model considers a Hagen-Poiseuille-type law to relate the volume flow from on pore body to another to discrete pressure drops $`\Delta p = p_i - p_j`$ between the pore bodies using throat conductance, $`g_{ij}`$. + +```math + Q_{ij} = g_{ij} (p_i - p_j), +``` +or + +```math + (p_i - p_j) = Q_{ij} / g_{ij} . +``` + +In the simulation of non-creeping flow, to capture inertial effects in fluid flow, an extension of the Hagen-Poiseuille-type law which includes the expansion and contraction of the pore space when moving from a throat to a pore and vice versa. ```math - Q_{ij} = g_{ij} (p_i - p_j). + (p_i - p_j) = Q_{ij} / g_{ij} + (C_{exp} + C_{cont})Q^2, ``` +where $`C_{exp}`$ and $`C_{cont}`$ are the expansion and the contraction coefficient, respectively. + # Implementation & Post processing In the following, we take a closer look at the source files for this example. diff --git a/examples/porenetwork_upscaling/doc/main.md b/examples/porenetwork_upscaling/doc/main.md index 7c85039daa388ce1d3efb57681d2ec576d0ce455..d484fe132dd19b6ea5958004d9f1a94ba9f29d8f 100644 --- a/examples/porenetwork_upscaling/doc/main.md +++ b/examples/porenetwork_upscaling/doc/main.md @@ -7,8 +7,8 @@ # Part 2: Main program flow The main program flow is implemented in file `main.cc` described below. -For each spatial direction x, y and z, flow through the network is simulated and the resulting mass flow rate -is used to determine the permeability. +For each spatial direction x, y and z, flow through the network is simulated using several pressure gradients and the resulting mass flow rates +are used to determine the upscaled properties. The code documentation is structured as follows: @@ -30,7 +30,10 @@ Pore-Network-Model to evaluate the upscaled Darcy permeability of a given networ #include +#include + #include // for floating point comparison +#include #include // for GetPropType #include // for getParam @@ -38,14 +41,14 @@ Pore-Network-Model to evaluate the upscaled Darcy permeability of a given networ #include // for ILU0BiCGSTABBackend #include // for LinearPDESolver +#include #include -#include #include -#include +#include +#include #include // for pore-network grid -#include #include // for getting the total mass flux leaving the network #include "upscalinghelper.hh" @@ -54,21 +57,18 @@ Pore-Network-Model to evaluate the upscaled Darcy permeability of a given networ -### Beginning of the main function +### The driver function -```cpp -int main(int argc, char** argv) try -{ - using namespace Dumux; +The template argument `TypeTag` determines if we run the example assuming +a creeping flow regime or not. Which regime is selected is set with the parameter +`Problem.AssumeCreepingFlow` in the input file. - // maybe initialize MPI and/or multithreading backend - Dumux::initialize(argc, argv); - - // We parse the command line arguments. - Parameters::init(argc, argv); +```cpp +namespace Dumux { - // Convenience alias for the type tag of the problem. - using TypeTag = Properties::TTag::PNMUpscaling; +template +void runExample() +{ ``` ### Create the grid and the grid geometry @@ -91,14 +91,12 @@ int main(int argc, char** argv) try ### Initialize the problem and grid variables ```cpp - using SpatialParams = GetPropType; - auto spatialParams = std::make_shared(gridGeometry); using Problem = GetPropType; - auto problem = std::make_shared(gridGeometry, spatialParams); + auto problem = std::make_shared(gridGeometry); // instantiate and initialize the discrete and exact solution vectors using SolutionVector = GetPropType; - SolutionVector x(gridGeometry->numDofs()); + SolutionVector x(gridGeometry->numDofs()); // zero-initializes // instantiate and initialize the grid variables using GridVariables = GetPropType; @@ -106,6 +104,23 @@ int main(int argc, char** argv) try gridVariables->init(x); ``` +### Instantiate the solver +We use the `NewtonSolver` class, which is instantiated on the basis +of an assembler and a linear solver. When the `solve` function of the +`NewtonSolver` is called, it uses the assembler and linear +solver classes to assemble and solve the non-linear system. + +```cpp + using Assembler = FVAssembler; + auto assembler = std::make_shared(problem, gridGeometry, gridVariables); + + using LinearSolver = UMFPackBackend; + auto linearSolver = std::make_shared(); + + using NewtonSolver = NewtonSolver; + NewtonSolver nonLinearSolver(assembler, linearSolver); +``` + ### Initialize VTK output ```cpp @@ -113,42 +128,31 @@ int main(int argc, char** argv) try using VtkWriter = PoreNetwork::VtkOutputModule, SolutionVector>; VtkWriter vtkWriter(*gridVariables, x, problem->name()); VtkOutputFields::initOutputModule(vtkWriter); - vtkWriter.addField(gridGeometry->poreVolume(), "poreVolume", Vtk::FieldType::vertex); - vtkWriter.addField(gridGeometry->throatShapeFactor(), "throatShapeFactor", Vtk::FieldType::element); - vtkWriter.addField(gridGeometry->throatCrossSectionalArea(), "throatCrossSectionalArea", Vtk::FieldType::element); ``` -### Instantiate the solver -We use the `LinearPDESolver` class, which is instantiated on the basis -of an assembler and a linear solver. When the `solve` function of the -`LinearPDESolver` is called, it uses the assembler and linear -solver classes to assemble and solve the linear system around the provided -solution and stores the result therein. +specify the field type explicitly since it may not be possible +to deduce this from the vector size in a pore network ```cpp - using Assembler = FVAssembler; - auto assembler = std::make_shared(problem, gridGeometry, gridVariables); - - using LinearSolver = UMFPackBackend; - auto linearSolver = std::make_shared(); - LinearPDESolver solver(assembler, linearSolver); - solver.setVerbose(false); // suppress output during solve() + vtkWriter.addField(gridGeometry->poreVolume(), "poreVolume", Vtk::FieldType::vertex); + vtkWriter.addField(gridGeometry->throatShapeFactor(), "throatShapeFactor", Vtk::FieldType::element); + vtkWriter.addField(gridGeometry->throatCrossSectionalArea(), "throatCrossSectionalArea", Vtk::FieldType::element); ``` ### Prepare the upscaling procedure. -Specify the directions for which the permeability shall be determined (default: x, y, z for 3D). +Set up a helper class to determine the total mass flux leaving the network ```cpp - using GridView = typename GetPropType::GridView; - const auto defaultDirections = GridView::dimensionworld == 3 ? std::vector{0, 1, 2} - : std::vector{0, 1}; - const auto directions = getParam>("Problem.Directions", defaultDirections); + const auto boundaryFlux = PoreNetwork::BoundaryFlux(*gridVariables, assembler->localResidual(), x); ``` -Set up a helper class to determine the total mass flux leaving the network +### The actual upscaling procedure +#### Instantiate the upscaling helper ```cpp - const auto boundaryFlux = PoreNetwork::BoundaryFlux(*gridVariables, assembler->localResidual(), x); + using Scalar = GetPropType; + using UpscalingHelper = UpscalingHelper; + UpscalingHelper upscalingHelper; ``` Set the side lengths used for applying the pressure gradient and calculating the REV outflow area. @@ -162,82 +166,103 @@ determine it automatically based on the network's bounding box. if (hasParam("Problem.SideLength")) return getParam("Problem.SideLength"); else - return UpscalingHelper::getSideLengths(*gridGeometry); + return upscalingHelper.getSideLengths(*gridGeometry); }(); // pass the side lengths to the problem problem->setSideLengths(sideLengths); ``` -### The actual upscaling procedure -Iterate over all directions specified before, apply the pressure gradient, calculated the mass flux -and finally determine the permeability. +Get the maximum and minimum pressure gradient and the population of sample points specified in the input file ```cpp + const Scalar minPressureGradient = getParam("Problem.MinimumPressureGradient", 1e1); + const Scalar maxPressureGradient = getParam("Problem.MaximumPressureGradient", 1e10); + + if (!(minPressureGradient < maxPressureGradient)) + DUNE_THROW(Dune::InvalidStateException, "Maximum pressure gradient must be greater than minimum pressure gradient"); + + const int numberOfSamples = getParam("Problem.NumberOfPressureGradients", 1); +``` + +Iterate over all directions specified before, apply several pressure gradient, calculated the mass flux +and finally determine the the upscaled properties. + +```cpp + const auto directions = getParam>("Problem.Directions", std::vector{0, 1, 2}); + upscalingHelper.setDirections(directions); for (int dimIdx : directions) { - // reset the solution - x = 0; - // set the direction in which the pressure gradient will be applied problem->setDirection(dimIdx); - // solve problem - solver.solve(x); + for (int i = 0; i < numberOfSamples; i++) + { + // reset the solution + x = 0; - // write the vtu file for the given direction - vtkWriter.write(dimIdx); + // set the pressure gradient to be applied + Scalar pressureGradient = maxPressureGradient * std::exp(i + 1 - numberOfSamples); + if (i == 0) + pressureGradient = std::min(minPressureGradient, pressureGradient); - // get the Scalar type - using Scalar = GetPropType; + problem->setPressureGradient(pressureGradient); - // calculate the permeability - const Scalar totalFluidMassFlux = boundaryFlux.getFlux(std::vector{problem->outletPoreLabel()})[0]; - const Scalar K = UpscalingHelper::getDarcyPermeability(*problem, totalFluidMassFlux); + // solve problem + nonLinearSolver.solve(x); - // optionally compare with reference data - static const auto referenceData = getParam>("Problem.ReferenceData", std::vector{}); - if (!referenceData.empty()) - { - static const Scalar eps = getParam("Problem.TestEpsilon"); - if (Dune::FloatCmp::ne(K, referenceData[dimIdx], eps)) - { - std::cerr << "Calculated permeability of " << K << " does not match with reference value of " << referenceData[dimIdx] << std::endl; - return 1; - } + // set the sample points + const Scalar totalFluidMassFlux = boundaryFlux.getFlux(std::vector{ problem->outletPoreLabel() })[0]; + upscalingHelper.setDataPoints(*problem, totalFluidMassFlux); } + + // write a vtu file for the given direction for the last sample + vtkWriter.write(dimIdx); } - // program end, return with 0 exit code (success) - return 0; -} + // calculate and report the upscaled properties + constexpr bool isCreepingFlow = std::is_same_v; + upscalingHelper.calculateUpscaledProperties(*problem, isCreepingFlow); + upscalingHelper.report(isCreepingFlow); + + // compare the Darcy permeability with reference data if provided in input file and report in case of inconsistency + static const auto referenceData = getParam>("Problem.ReferencePermeability", std::vector{}); + if (!referenceData.empty()) + upscalingHelper.compareWithReference(referenceData); + + // plot the results just for non-creeping flow + // creeping flow would just result in a straight line (permeability is independent of the pressure gradient) + if (!isCreepingFlow) + upscalingHelper.plot(); +}; + +} // end namespace Dumux ``` -### Exception handling -In this part of the main file we catch and print possible exceptions that could -occur during the simulation. -
Click to show error handler +### The main function +
Click to show main ```cpp - -catch (const Dumux::ParameterException &e) +int main(int argc, char** argv) { - std::cerr << std::endl << e << " ---> Abort!" << std::endl; - return 1; -} -catch (const Dune::DGFException & e) -{ - std::cerr << "DGF exception thrown (" << e << - "). Most likely, the DGF file name is wrong " - "or the DGF file is corrupted, " - "e.g. missing hash at end of file or wrong number (dimensions) of entries." - << " ---> Abort!" << std::endl; - return 2; -} -catch (const Dune::Exception &e) -{ - std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl; - return 3; + using namespace Dumux; + + // We parse the command line arguments. + Parameters::init(argc, argv); + + // Convenience alias for the type tag of the problem. + using CreepingFlowTypeTag = Properties::TTag::PNMUpscalingCreepingFlow; + using NonCreepingFlowTypeTag = Properties::TTag::PNMUpscalingNonCreepingFlow; + // // [[/codeblock]] + + // user decides whether creeping flow or non-creeping flow should be run + if (getParam("Problem.AssumeCreepingFlow", false)) + runExample(); + else + runExample(); + + // program end, return with 0 exit code (success) + return 0; } ``` diff --git a/examples/porenetwork_upscaling/doc/main_intro.md b/examples/porenetwork_upscaling/doc/main_intro.md index e793cab0321bb575a81884d3af0437a445b6a9e0..bb5d5dfd270ad5ef45fbc3e306ea778efef2e22a 100644 --- a/examples/porenetwork_upscaling/doc/main_intro.md +++ b/examples/porenetwork_upscaling/doc/main_intro.md @@ -1,8 +1,8 @@ # Part 2: Main program flow The main program flow is implemented in file `main.cc` described below. -For each spatial direction x, y and z, flow through the network is simulated and the resulting mass flow rate -is used to determine the permeability. +For each spatial direction x, y and z, flow through the network is simulated using several pressure gradients and the resulting mass flow rates +are used to determine the upscaled properties. The code documentation is structured as follows: diff --git a/examples/porenetwork_upscaling/doc/problem.md b/examples/porenetwork_upscaling/doc/problem.md index 7a5514983565fb2146b2ecdd2a5a9fb743550727..778c5e350e2312437d4b6ab3ff471134ff713589 100644 --- a/examples/porenetwork_upscaling/doc/problem.md +++ b/examples/porenetwork_upscaling/doc/problem.md @@ -35,6 +35,19 @@ type tag, which we want to modify or for which no meaningful default can be set. #include // for `TTag::PNMOneP` ``` +The class that contains a collection of single-phase flow throat transmissibilities +among them the transmisibility model to be used can be specified in AdvectionType class + +```cpp +#include +``` + +The class that provides specializations for both creeping and non-creeping advection types. + +```cpp +#include +``` + The local residual for incompressible flow is included. The one-phase flow model (included above) uses a default implementation of the local residual for single-phase flow. However, in this example we are using an @@ -63,13 +76,15 @@ The classes that define the problem and parameters used in this simulation
### `TypeTag` definition -A `TypeTag` for our simulation is defined, which inherits properties from the -single-phase flow model and the box scheme. +Two `TypeTag` for our simulation are defined, one for creeping flow and another for non-creeping flow, +which inherit properties from the single-phase pore network model. The non-creeping flow inherits +all properties from the creeping flow simulation but sets an own property for the `AdvectionType`. ```cpp namespace Dumux::Properties { namespace TTag { -struct PNMUpscaling { using InheritsFrom = std::tuple; }; +struct PNMUpscalingCreepingFlow { using InheritsFrom = std::tuple; }; +struct PNMUpscalingNonCreepingFlow { using InheritsFrom = std::tuple; }; } ``` @@ -81,43 +96,49 @@ default can be set, are specialized for our type tag `PNMUpscaling`. ```cpp // We use `dune-foamgrid`, which is especially tailored for 1D networks. template -struct Grid +struct Grid { using type = Dune::FoamGrid<1, 3>; }; // The problem class specifying initial and boundary conditions: template -struct Problem +struct Problem { using type = UpscalingProblem; }; -//! The spatial parameters to be employed. +//! The spatial parameters template -struct SpatialParams +struct SpatialParams { -private: using GridGeometry = GetPropType; using Scalar = GetPropType; public: using type = PoreNetwork::UpscalingSpatialParams; }; -//! The advection type. +//! The advection type for creeping flow template -struct AdvectionType +struct AdvectionType { -private: using Scalar = GetPropType; using TransmissibilityLaw = PoreNetwork::TransmissibilityPatzekSilin; public: using type = PoreNetwork::CreepingFlow; }; -// We use a single liquid phase consisting of a component with constant fluid properties. +//! The advection type for non-creeping flow (includes model for inertia effects) template -struct FluidSystem +struct AdvectionType { -private: using Scalar = GetPropType; + using TransmissibilityLaw = PoreNetwork::TransmissibilityPatzekSilin; public: + using type = PoreNetwork::NonCreepingFlow; +}; + +// We use a single liquid phase consisting of a component with constant fluid properties. +template +struct FluidSystem +{ + using Scalar = GetPropType; using type = FluidSystems::OnePLiquid >; }; ``` @@ -127,7 +148,7 @@ that contains functionality related to analytic differentiation. ```cpp template -struct LocalResidual +struct LocalResidual { using type = OnePIncompressibleLocalResidual; }; } // end namespace Dumux::Properties @@ -169,46 +190,60 @@ class UpscalingProblem : public PorousMediumFlowProblem
Click to show convenience aliases ```cpp - using ParentType = PorousMediumFlowProblem; - using GridGeometry = GetPropType; - using FVElementGeometry = typename GetPropType::LocalView; - using SubControlVolume = typename FVElementGeometry::SubControlVolume; using Scalar = GetPropType; + using GridGeometry = GetPropType; + using FVElementGeometry = typename GridGeometry::LocalView; + using SubControlVolume = typename GridGeometry::SubControlVolume; + using Element = typename GridGeometry::GridView::template Codim<0>::Entity; + using GlobalPosition = typename Element::Geometry::GlobalCoordinate; using PrimaryVariables = GetPropType; using Indices = typename GetPropType::Indices; using BoundaryTypes = Dumux::BoundaryTypes; - using Element = typename GridGeometry::GridView::template Codim<0>::Entity; - using GlobalPosition = typename Element::Geometry::GlobalCoordinate; ```
-### The constructor of our problem. + +#### The constructor of our problem. ```cpp public: - template - UpscalingProblem(std::shared_ptr gridGeometry, std::shared_ptr spatialParams) - : ParentType(gridGeometry, spatialParams) + UpscalingProblem(std::shared_ptr gridGeometry) + : ParentType(gridGeometry) { // the applied pressure gradient - pressureGradient_ = getParam("Problem.PressureGradient"); + pressureGradient_ = getParam("Problem.PressureGradient", 1e5); // We can either use pore labels (given in the grid file) to identify inlet and outlet pores // or use the network's bounding box to find these pores automatically. Using labels is usually much // more accurate, so this is the default here. useLabels_ = getParam("Problem.UseLabels", true); - // an epsilon value for the bounding box approach + // an epsilon value for the floating point comparisons to determine inlet/outlet pores eps_ = getParam("Problem.Epsilon", 1e-7); } ``` +#### Temperature +We need to specify a constant temperature for our isothermal problem. +Fluid properties that depend on temperature will be calculated with this value. + +```cpp + Scalar temperature() const + { return 283.15; } +``` + +Set the pressure gradient to be applied to the network + +```cpp + void setPressureGradient(Scalar pressureGradient) + { pressureGradient_ = pressureGradient; } +``` + #### Boundary conditions This function is used to define the __type of boundary conditions__ used depending on the location. -Here, we use Dirichlet boundary conditions (fixed pressures) at the inlet and outlet and Neumann -boundary conditions at all remaining boundaries. Note that the PNM only supports Neumann no-flow boundaries. -The specify a certain mass flux, we would have to use a source term on the boundary pores (which is not done in this example). +Here, we use Dirichlet boundary conditions (fixed pressures) at the inlet and outlet. Note that the PNM does not support Neumann boundaries. +To specify a certain mass flux on a boundary, we would have to use a source term on the boundary pores (which is not done in this example). ```cpp BoundaryTypes boundaryTypes(const Element &element, const SubControlVolume& scv) const @@ -218,8 +253,6 @@ The specify a certain mass flux, we would have to use a source term on the bound // fix the pressure at the inlet and outlet pores if (isInletPore_(scv)|| isOutletPore_(scv)) bcTypes.setAllDirichlet(); - else // Neumann (no-flow) for the remaining boundaries - bcTypes.setAllNeumann(); return bcTypes; } @@ -282,16 +315,23 @@ and the length of the domain at the inlet. // Return the applied pressure gradient. Scalar pressureGradient() const { return pressureGradient_; } +``` +
+ +Return the label of inlet pores assuming a previously set direction. - // Return the label of inlet pores assuming a previously set direction. +```cpp int inletPoreLabel() const { static constexpr std::array label = {1, 3, 5}; return label[direction_]; } +``` + +Return the label of outlet pores assuming a previously set direction. - // Return the label of outlet pores assuming a previously set direction. +```cpp int outletPoreLabel() const { static constexpr std::array label = {2, 4, 6}; @@ -315,25 +355,21 @@ private: else return scv.dofPosition()[direction_] > this->gridGeometry().bBoxMax()[direction_] - eps_; } +``` + +private data members - // private data members +```cpp Scalar eps_; Scalar pressureGradient_; int direction_; GlobalPosition length_; bool useLabels_; -``` - - - -```cpp - }; } // end namespace Dumux ``` -[[/codeblock]] diff --git a/examples/porenetwork_upscaling/doc/upscalinghelper.md b/examples/porenetwork_upscaling/doc/upscalinghelper.md index 8f5ec93b78060be0649270e9ed4b482f9ce1ecc7..7b671b301a78f462acc72013c1e8258c7828800c 100644 --- a/examples/porenetwork_upscaling/doc/upscalinghelper.md +++ b/examples/porenetwork_upscaling/doc/upscalinghelper.md @@ -4,27 +4,37 @@ | [:arrow_left: Back to the main documentation](../README.md) | [:arrow_left: Go back to part 2](main.md) | |---|---:| -# Part 3: Upscaling helper +The upscaling helper evaluates the pore-network simulation results for each direction $`i`$ and calculates the upscaled properties in this direction. Firstly, it evaluates the apparent velocity as: -The upscaling helper evaluates the pore-network simulation results for each direction $`i`$ and calculates the upscaled intrinsic permeability in this direction using: +```math + v_{\mathrm{Apparent},i} = \frac{q_{\mathrm{mass,tot},i} / \varrho}{A_{\mathrm{tot},i}} +``` + +where $`q_{\mathrm{mass,tot},i}`$ is the total mass flow leaving the network over the REV's boundary with area +$`A_{\mathrm{tot},i}`$ in $`i`$-direction. $`\varrho `$ is the fluid mass density. Then, we calculate upscaled permeability as: ```math - K_i = v_{\mathrm{Darcy},i} / \nabla p_i ~ \mu. + K_D = v_{\mathrm{Apparent},i} ~ \mu / \nabla p_i. ``` -$`\nabla p_i`$ is a given pressure gradient in $`i`$-direction and $`\mu`$ the fluid dynamic viscosity. - -We evaluate the Darcy velocity as +$`\nabla p_i`$ is a given pressure gradient in $`i`$-direction and $`\mu`$ the fluid dynamic viscosity. In the creeping flow simulation, the calculated permeability, $`K_D`$, is the Darcy (intrinsic) permeability of the system. +To simulate non-creeping flow, we use Forchheimer's equation to upscale the properties. ```math - v_{\mathrm{Darcy},i} = \frac{q_{\mathrm{mass,tot},i} / \varrho}{A_{\mathrm{tot},i}} + \nabla p_i = \frac{\mu}{K_f} v_{\mathrm{Apparent},i} + \varrho \beta v_{\mathrm{Apparent},i}^2, ``` +where $`K_f`$ is the Forchheimer permeability and $`\beta`$ is the Forchheimer coefficient. -where $`q_{\mathrm{mass,tot},i}`$ is the total mass flow leaving the network over the REV's boundary with area -$`A_{\mathrm{tot},i}`$ in $`i`$-direction. $`\varrho `$ is the fluid mass density. +Although it is sometimes assumed for the sake of simplicity that $`K_f = K_D`$, they are not exactly the same properties. As the velocity increases, the flow regime in a porous medium shifts from the Darcy to the Forchheimer regime. This change in the flow regime means that in addition to viscous dissipation (the first term in the Forchheimer equation), we also need to consider inertia effects (the second term in the Forchheimer equation). Moreover, the shift in the flow regime, means a different viscous dissipation than what the same porous medium experiences a Darcy flow. In other words, moving from the Darcy to the Forchheimer regime establishes a new velocity field in the porous medium which causes the difference in viscous dissipation in addition to inertia effects. The first and second terms in the Forchheimer equation are correlated. For more detail, we refer to the study conducted by [Dukhan and Minjeur (2010)](http://dx.doi.org/10.1007/s10934-010-9393-1). -The code documentation is structured as follows: +In this example the end of the creeping flow (Darcy regime) is defined as the moment when the pressure drop calculated by the Dracy (intrinsic) permeability $`\Delta p_i = v_{\mathrm{Apparent},i} \mu l/ K_D`$ becomes less than 99% of the total pressure drop [Muljadi et al.](https://doi.org/10.1016/j.advwatres.2015.05.019). +To calculate upscaled properties, we rearrange Forchehimer's equation as: -[[_TOC_]] +```math + \frac{\nabla p_i v_{\mathrm{Apparent},i}}{\mu} = \frac{1}{K_f} + \frac{\varrho v_{\mathrm{Apparent},i}}{\mu} \beta . +``` +Finding the linear regression line of $`\nabla p_i v_{\mathrm{Apparent},i}/\mu `$ versus $`\varrho v_{\mathrm{Apparent},i}/\mu `$ and using the intercept and the slope of the regression line, we can respectively calculate the Forchheimer permeability and coefficient. It should be noted that the calculation of the Forchheimer permeability can be affected by the pressure range applied to the porous medium as well as the number of sample points that are used in the regression process. We compute the Darcy (intrinsic) permeability as the maximum permeability of the sample of data of the system which happens when the pressure gradient is small enough such that inertial effects are negligible. To ensure such a small pressure gradient, we set the first pressure gradient to be applied as $'10 Pa/m'$. However, this value can be adapted using the keyword `Problem.MinimumPressureGradient` in params.input. it is recommended to use more than 10 pressure sample points which can be set in the input file. As mentioned before, considering a slight difference between Darcy (intrinsic) permeability and Forchheimer permeability, in many applications they can be used interchangeabely. Here, however, we distinguish between them, calculate and report them separately. + +The code documentation is structured as follows: ## Upscaling helper struct (`upscalinghelper.hh`) @@ -42,57 +52,114 @@ the pore network in flow direction in order to find the upscaled Darcy permeabil #include #include #include +#include +#include namespace Dumux { -struct UpscalingHelper +template +class UpscalingHelper { +public: ``` -### Calculate the intrinsic permeability +### Set data points to calculate intrinsic permeability and Forchheimer coefficient This function first evaluates the mass flux leaving the network in the direction of the applied pressure gradient. -Afterwards, the mass flux is converted into an area specify volume flux from which finally the intrinsic Darcy -permeability $`\mathbf{K}`$ [m$`^2`$] can be evaluated. +Afterwards, the mass flux is converted into an volume flux which is used to calculate the apparent velocity. +Then apparent permeability of the network is computed and stored for furthure calculations. ```cpp - template - static Scalar getDarcyPermeability(const Problem& problem, const Scalar totalMassFlux) + template + void setDataPoints(const Problem &problem, const Scalar totalMassFlux) { // get the domain side lengths from the problem auto sideLengths = problem.sideLengths(); - // create temporary stringstream with fixed scientific formatting without affecting std::cout - std::ostream tmp(std::cout.rdbuf()); - tmp << std::fixed << std::scientific; - static constexpr char dirNames[] = "xyz"; + + // get the applied pressure gradient + const auto pressureGradient = problem.pressureGradient(); + const auto pressureDrop = pressureGradient * sideLengths[problem.direction()]; + + // get the fluid properties + const auto liquidDensity = problem.liquidDensity(); + const auto liquidDynamicViscosity = problem.liquidDynamicViscosity(); // convert mass to volume flux - const auto volumeFlux = totalMassFlux / problem.liquidDensity();; + const auto volumeFlux = totalMassFlux / liquidDensity; + ; + // calculate apparent velocity sideLengths[problem.direction()] = 1.0; const auto outflowArea = std::accumulate(sideLengths.begin(), sideLengths.end(), 1.0, std::multiplies()); - const auto vDarcy = volumeFlux / outflowArea; - const auto K = vDarcy / problem.pressureGradient() * problem.liquidDynamicViscosity(); - tmp << "\n########################################\n" << std::endl; - tmp << dirNames[problem.direction()] << "-direction"; - tmp << ": Area = " << outflowArea << " m^2"; - tmp << "; Massflux = " << totalMassFlux << " kg/s"; - tmp << "; v_Darcy = " << vDarcy << " m/s"; - tmp << "; K = " << K << " m^2" << std::endl; - tmp << "\n########################################\n" << std::endl; - - return K; + const auto vApparent = volumeFlux / outflowArea; + + // compute apparent permeability + const auto KApparent = vApparent / pressureGradient * liquidDynamicViscosity; + // calculate rho v / mu, called inertia to viscous ratio in the rest of the code + const auto inertiaToViscousRatio = liquidDensity * vApparent / liquidDynamicViscosity; + + // store the required data for further calculations + totalPressureDrop_[problem.direction()].push_back(pressureDrop); + apparentVelocity_[problem.direction()].push_back(vApparent); + apparentPermeability_[problem.direction()].push_back(KApparent); + inertiaToViscousRatio_[problem.direction()].push_back(inertiaToViscousRatio); } ``` -### Determine the domain's side lengths automatically based on the bounding box of the network. +### Calculate intrinsic permeability and Forchheimer coefficient. +This function first calculate intrinsic permeability and Forchheimer coefficient using linear least squares regression method +and reports them. It also plot the apparent permeability of the porous medium versus Forchheimer number/pressure gradient in each +simulation. + +```cpp + template + void calculateUpscaledProperties(const Problem &problem, bool isCreepingFlow) + { + const auto sideLengths = problem.sideLengths(); + const auto liquidDynamicViscosity = problem.liquidDynamicViscosity(); + + for (const auto dirIdx : directions_) + { + // determine Darcy permeability as the maximum permeability of the domain + darcyPermeability_[dirIdx] = *max_element(apparentPermeability_[dirIdx].begin(), apparentPermeability_[dirIdx].end()); + if (!isCreepingFlow) + { + for (int i = 0; i < totalPressureDrop_[dirIdx].size(); i++) + { + // calculate the Darcy pressure drop. + const Scalar darcyPressureDrop = liquidDynamicViscosity * apparentVelocity_[dirIdx][i] * sideLengths[dirIdx] / darcyPermeability_[dirIdx]; + + // calculate the ratio of Dracy to total pressure drop + const Scalar pressureDropRatio = darcyPressureDrop / totalPressureDrop_[dirIdx][i]; + + // set sample points for upscaling of Forchheimer parameters. + // first, check the permability ratio to see if the flow regime is Forchheimer. + if (pressureDropRatio < 0.99) + { + samplePointsX_[dirIdx].push_back(inertiaToViscousRatio_[dirIdx][i]); + samplePointsY_[dirIdx].push_back(1 / apparentPermeability_[dirIdx][i]); + } + } + // determine regression line and accordingly the Forchheimer permeability and the Forchheimer coefficient + const auto [intercept, slope] = linearRegression(samplePointsX_[dirIdx], samplePointsY_[dirIdx]); + forchheimerPermeability_[dirIdx] = 1.0 / intercept; + forchheimerCoefficient_[dirIdx] = slope; + writePlotDataToFile(dirIdx); + } + } + } +``` + +### Determine the domain's side lengths + +We determine the domain side length by using the bounding box of the network ```cpp template - static auto getSideLengths(const GridGeometry& gridGeometry) + auto getSideLengths(const GridGeometry& gridGeometry) { using GlobalPosition = typename GridGeometry::GlobalCoordinate; - GlobalPosition result; + GlobalPosition result(0.0); std::cout << "Automatically determining side lengths of REV based on bounding box of pore network" << std::endl; for (int dimIdx = 0; dimIdx < GridGeometry::GridView::dimensionworld; ++dimIdx) @@ -102,13 +169,231 @@ permeability $`\mathbf{K}`$ [m$`^2`$] can be evaluated. } ``` +### Plot the data using Gnuplot + + +```cpp + void plot() + { + // plot permeability ratio vs. Forchheimer number + plotPermeabilityratioVsForchheimerNumber_(); + + // plot inverse of apparent permability vs. rho v / mu + plotInversePrmeabilityVsInertiaToViscousRatio_(); + } +``` + + +### Save the relevant data for plot + +```cpp + void writePlotDataToFile(std::size_t dirIdx) + { + // write permeability ratio vs. Forchheimer number + writePermeabilityratioVsForchheimerNumber_(dirIdx); + + // write inverse of apparent permability vs. rho v / mu + writeInversePrmeabilityVsInertiaToViscousRatio_(dirIdx); + } +``` + + +### Report the upscaled data + +```cpp + void report(bool isCreepingFlow) + { + // Report the results for each direction + for (const auto dirIdx : directions_) + { + std::cout << Fmt::format("\n{:#>{}}\n\n", "", 40) + << Fmt::format("{}-direction:\n", dirName_[dirIdx]) + << Fmt::format("-- Darcy (intrinsic) permeability = {:.3e} m^2\n", darcyPermeability_[dirIdx]); + + // Report non-creeping flow upscaled properties + if (!isCreepingFlow) + { + std::cout << Fmt::format("-- Forchheimer permeability = {:.3e} m^2\n", forchheimerPermeability_[dirIdx]); + std::cout << Fmt::format("-- Forchheimer coefficient = {:.3e} m^-1\n", forchheimerCoefficient_[dirIdx]); + } + + std::cout << Fmt::format("\n{:#>{}}\n", "", 40) << std::endl; + } + } +``` + + +### Compare with reference data provided in input file + +```cpp + void compareWithReference(std::vector referenceData) + { + for (const auto dirIdx : directions_) + { + const auto K = darcyPermeability_[dirIdx]; + static const Scalar eps = getParam("Problem.TestEpsilon", 1e-3); + if (Dune::FloatCmp::ne(K, referenceData[dirIdx], eps)) + { + std::cerr << "Calculated permeability of " << K << " in " + < directions) + { + directions_ = directions; + } + // [[codeblock]] + +private: + // ### Save the relevant data for plot of permeability ratio vs. Forchheimer number + // [[codeblock]] + void writePermeabilityratioVsForchheimerNumber_(std::size_t dirIdx) + { + // open a logfile + std::ofstream logfile(dirName_[dirIdx] + "-dir-PermeabilityratioVsForchheimerNumber.dat"); + + // save the data needed to be plotted in logfile + for (int i = 0; i < apparentPermeability_[dirIdx].size(); i++) + { + // compute the Forchheimer number + const Scalar forchheimerNumber = darcyPermeability_[dirIdx] * forchheimerCoefficient_[dirIdx] * inertiaToViscousRatio_[dirIdx][i]; + // ratio between apparrent permeability and darcy permeability + const Scalar permeabilityRatio = apparentPermeability_[dirIdx][i] / darcyPermeability_[dirIdx]; + + logfile << forchheimerNumber << " " << permeabilityRatio << std::endl; + } + } +``` + +### Save the relevant data for plot of inverse of apparent permability vs. rho v / mu + +```cpp + void writeInversePrmeabilityVsInertiaToViscousRatio_(std::size_t dirIdx) + { + // open a logfile and write inverese of apparent permeability given by the model vs. inertial to viscous ratio (rho v / mu) + std::ofstream logfile(dirName_[dirIdx] + "-dir-InversePrmeabilityVsInertiaToViscousRatio.dat"); + + // save the data needed to be plotted in logfile + for (int i = 0; i < apparentPermeability_[dirIdx].size(); i++) + { + const Scalar inertiaToViscousRatio = inertiaToViscousRatio_[dirIdx][i]; + const Scalar inverseAppPermeability = 1 / apparentPermeability_[dirIdx][i]; + + // compute inverse of apparent permeability using the Forchheimer permeability and coefficient + const Scalar inverseAppPermeabilityForchheimer = 1 / forchheimerPermeability_[dirIdx] + inertiaToViscousRatio * forchheimerCoefficient_[dirIdx]; + + logfile << inertiaToViscousRatio << " " << 1e-12 * inverseAppPermeability << " " << 1e-12 * inverseAppPermeabilityForchheimer << std::endl; + } + } +``` + + +### Plot permeability ratio vs. Forchheimer number using Gnuplot + + +```cpp + void plotPermeabilityratioVsForchheimerNumber_() + { + // using gnuplot interface + Dumux::GnuplotInterface gnuplot(true); + gnuplot.setOpenPlotWindow(true); + + for (const auto dirIdx : directions_) + { + gnuplot.resetAll(); + std::string title{}, option{}; + + // add the data in each direction for plot + gnuplot.addFileToPlot(dirName_[dirIdx] + "-dir-PermeabilityratioVsForchheimerNumber.dat", "notitle with lines"); + // set the properties of lines to be plotted + option += "set linetype 1 linecolor 1 linewidth 7\n"; + // report the darcy permeability in each direction as the title of the plot + title += Fmt::format("{}-direction, Darcy permeability= {:.3e} m^2 ", dirName_[dirIdx], darcyPermeability_[dirIdx]); + + option += "set title \"" + title + "\"\n"; + option += "set logscale x""\n"; + option += "set format x '10^{%L}'""\n"; + + gnuplot.setXlabel("Forchheimer Number [-]"); + gnuplot.setYlabel("Apparent permeability / Darcy permeability [-]"); + gnuplot.setOption(option); + gnuplot.plot("permeability_ratio_versus_forchheimer_number"); + } + } +``` + + +### Plot inverse of apparent permability vs. rho v / mu using Gnuplot + + +```cpp + void plotInversePrmeabilityVsInertiaToViscousRatio_() + { + // using gnuplot interface + Dumux::GnuplotInterface gnuplot(true); + gnuplot.setOpenPlotWindow(true); + + for (const auto dirIdx : directions_) + { + gnuplot.resetAll(); + std::string title{}, option{}; + std::string legend0 = "u 1:2 title \"Network model\" with lines"; + // add the data in each direction for plot, first set of data + gnuplot.addFileToPlot(dirName_[dirIdx] + "-dir-InversePrmeabilityVsInertiaToViscousRatio.dat", legend0); + + // set the properties of lines to be plotted + option += "set linetype 1 linecolor 1 linewidth 5\n"; + + std::string legend1 = "u 1:3 title \"Forchheimer equation\" with lines"; + // add the data in each direction for plot, second set of data + gnuplot.addFileToPlot(dirName_[dirIdx] + "-dir-InversePrmeabilityVsInertiaToViscousRatio.dat", legend1); + + // set the properties of lines to be plotted + option += "set linetype 2 linecolor 2 linewidth 5\n"; + + // report the darcy permeability in each direction as the title of the plot + title += Fmt::format("{}-direction, Darcy permeability= {:.3e} m^2 ", dirName_[dirIdx], darcyPermeability_[dirIdx]); + + option += "set title \"" + title + "\"\n"; + option += "set logscale x""\n"; + option += "set format x '10^{%L}'""\n"; + + gnuplot.setXlabel("{/Symbol r} v / {/Symbol m} [1/m]"); + gnuplot.setYlabel("1/ Apparent permeability [1/m^2] x 1e12"); + gnuplot.setOption(option); + gnuplot.plot("inverse_apppermeability_versus_rhovmu-1"); + } + } +``` + ```cpp + std::array, 3> samplePointsX_; + std::array, 3> samplePointsY_; + std::array, 3>totalPressureDrop_; + std::array, 3> apparentVelocity_; + std::array, 3> apparentPermeability_; + std::array, 3> inertiaToViscousRatio_; + std::array darcyPermeability_; + std::array forchheimerPermeability_; + std::array forchheimerCoefficient_; + const std::array dirName_ = {"X", "Y", "Z"}; + std::vector directions_; }; } // end namespace Dumux ``` +[[/codeblock]] diff --git a/examples/porenetwork_upscaling/doc/upscalinghelper_intro.md b/examples/porenetwork_upscaling/doc/upscalinghelper_intro.md index fbafad4fab072239e69b03207a3994e3bd298025..76c143622b1edb753a7bd4fd153ffafb0ace4b63 100644 --- a/examples/porenetwork_upscaling/doc/upscalinghelper_intro.md +++ b/examples/porenetwork_upscaling/doc/upscalinghelper_intro.md @@ -1,21 +1,31 @@ -# Part 3: Upscaling helper +The upscaling helper evaluates the pore-network simulation results for each direction $`i`$ and calculates the upscaled properties in this direction. Firstly, it evaluates the apparent velocity as: -The upscaling helper evaluates the pore-network simulation results for each direction $`i`$ and calculates the upscaled intrinsic permeability in this direction using: +```math + v_{\mathrm{Apparent},i} = \frac{q_{\mathrm{mass,tot},i} / \varrho}{A_{\mathrm{tot},i}} +``` + +where $`q_{\mathrm{mass,tot},i}`$ is the total mass flow leaving the network over the REV's boundary with area +$`A_{\mathrm{tot},i}`$ in $`i`$-direction. $`\varrho `$ is the fluid mass density. Then, we calculate upscaled permeability as: ```math - K_i = v_{\mathrm{Darcy},i} / \nabla p_i ~ \mu. + K_D = v_{\mathrm{Apparent},i} ~ \mu / \nabla p_i. ``` -$`\nabla p_i`$ is a given pressure gradient in $`i`$-direction and $`\mu`$ the fluid dynamic viscosity. - -We evaluate the Darcy velocity as +$`\nabla p_i`$ is a given pressure gradient in $`i`$-direction and $`\mu`$ the fluid dynamic viscosity. In the creeping flow simulation, the calculated permeability, $`K_D`$, is the Darcy (intrinsic) permeability of the system. +To simulate non-creeping flow, we use Forchheimer's equation to upscale the properties. ```math - v_{\mathrm{Darcy},i} = \frac{q_{\mathrm{mass,tot},i} / \varrho}{A_{\mathrm{tot},i}} + \nabla p_i = \frac{\mu}{K_f} v_{\mathrm{Apparent},i} + \varrho \beta v_{\mathrm{Apparent},i}^2, ``` +where $`K_f`$ is the Forchheimer permeability and $`\beta`$ is the Forchheimer coefficient. -where $`q_{\mathrm{mass,tot},i}`$ is the total mass flow leaving the network over the REV's boundary with area -$`A_{\mathrm{tot},i}`$ in $`i`$-direction. $`\varrho `$ is the fluid mass density. +Although it is sometimes assumed for the sake of simplicity that $`K_f = K_D`$, they are not exactly the same properties. As the velocity increases, the flow regime in a porous medium shifts from the Darcy to the Forchheimer regime. This change in the flow regime means that in addition to viscous dissipation (the first term in the Forchheimer equation), we also need to consider inertia effects (the second term in the Forchheimer equation). Moreover, the shift in the flow regime, means a different viscous dissipation than what the same porous medium experiences a Darcy flow. In other words, moving from the Darcy to the Forchheimer regime establishes a new velocity field in the porous medium which causes the difference in viscous dissipation in addition to inertia effects. The first and second terms in the Forchheimer equation are correlated. For more detail, we refer to the study conducted by [Dukhan and Minjeur (2010)](http://dx.doi.org/10.1007/s10934-010-9393-1). -The code documentation is structured as follows: +In this example the end of the creeping flow (Darcy regime) is defined as the moment when the pressure drop calculated by the Dracy (intrinsic) permeability $`\Delta p_i = v_{\mathrm{Apparent},i} \mu l/ K_D`$ becomes less than 99% of the total pressure drop [Muljadi et al.](https://doi.org/10.1016/j.advwatres.2015.05.019). +To calculate upscaled properties, we rearrange Forchehimer's equation as: -[[_TOC_]] +```math + \frac{\nabla p_i v_{\mathrm{Apparent},i}}{\mu} = \frac{1}{K_f} + \frac{\varrho v_{\mathrm{Apparent},i}}{\mu} \beta . +``` +Finding the linear regression line of $`\nabla p_i v_{\mathrm{Apparent},i}/\mu `$ versus $`\varrho v_{\mathrm{Apparent},i}/\mu `$ and using the intercept and the slope of the regression line, we can respectively calculate the Forchheimer permeability and coefficient. It should be noted that the calculation of the Forchheimer permeability can be affected by the pressure range applied to the porous medium as well as the number of sample points that are used in the regression process. We compute the Darcy (intrinsic) permeability as the maximum permeability of the sample of data of the system which happens when the pressure gradient is small enough such that inertial effects are negligible. To ensure such a small pressure gradient, we set the first pressure gradient to be applied as $'10 Pa/m'$. However, this value can be adapted using the keyword `Problem.MinimumPressureGradient` in params.input. it is recommended to use more than 10 pressure sample points which can be set in the input file. As mentioned before, considering a slight difference between Darcy (intrinsic) permeability and Forchheimer permeability, in many applications they can be used interchangeabely. Here, however, we distinguish between them, calculate and report them separately. + +The code documentation is structured as follows: diff --git a/examples/porenetwork_upscaling/img/inverse_apppermeability_versus_rhovmu-1.png b/examples/porenetwork_upscaling/img/inverse_apppermeability_versus_rhovmu-1.png new file mode 100644 index 0000000000000000000000000000000000000000..470e2935de64e0e2b060071ac0457abadbd322e5 Binary files /dev/null and b/examples/porenetwork_upscaling/img/inverse_apppermeability_versus_rhovmu-1.png differ diff --git a/examples/porenetwork_upscaling/img/permeability_ratio_versus_forchheimer_number.png b/examples/porenetwork_upscaling/img/permeability_ratio_versus_forchheimer_number.png new file mode 100644 index 0000000000000000000000000000000000000000..f597a47586dd9eb89ac82dc564975a6840f7d0cf Binary files /dev/null and b/examples/porenetwork_upscaling/img/permeability_ratio_versus_forchheimer_number.png differ diff --git a/examples/porenetwork_upscaling/main.cc b/examples/porenetwork_upscaling/main.cc index 2cedb4fdfbb97986b4475e8f816ea66d17fb2cd7..540109b475ddf409ebde8a375a9d49f5cb7d993e 100644 --- a/examples/porenetwork_upscaling/main.cc +++ b/examples/porenetwork_upscaling/main.cc @@ -27,7 +27,10 @@ #include +#include + #include // for floating point comparison +#include #include // for GetPropType #include // for getParam @@ -35,14 +38,14 @@ #include // for ILU0BiCGSTABBackend #include // for LinearPDESolver +#include #include -#include #include -#include +#include +#include #include // for pore-network grid -#include #include // for getting the total mass flux leaving the network #include "upscalinghelper.hh" @@ -50,22 +53,16 @@ // [[/codeblock]] // [[/details]] // -// ### Beginning of the main function -// [[codeblock]] -int main(int argc, char** argv) try -{ - using namespace Dumux; - - // maybe initialize MPI and/or multithreading backend - Dumux::initialize(argc, argv); - - // We parse the command line arguments. - Parameters::init(argc, argv); - - // Convenience alias for the type tag of the problem. - using TypeTag = Properties::TTag::PNMUpscaling; - // [[/codeblock]] +// ### The driver function +// +// The template argument `TypeTag` determines if we run the example assuming +// a creeping flow regime or not. Which regime is selected is set with the parameter +// `Problem.AssumeCreepingFlow` in the input file. +namespace Dumux { +template +void runExample() +{ // ### Create the grid and the grid geometry // [[codeblock]] // The grid manager can be used to create a grid from the input file @@ -84,14 +81,12 @@ int main(int argc, char** argv) try // ### Initialize the problem and grid variables // [[codeblock]] - using SpatialParams = GetPropType; - auto spatialParams = std::make_shared(gridGeometry); using Problem = GetPropType; - auto problem = std::make_shared(gridGeometry, spatialParams); + auto problem = std::make_shared(gridGeometry); // instantiate and initialize the discrete and exact solution vectors using SolutionVector = GetPropType; - SolutionVector x(gridGeometry->numDofs()); + SolutionVector x(gridGeometry->numDofs()); // zero-initializes // instantiate and initialize the grid variables using GridVariables = GetPropType; @@ -99,43 +94,45 @@ int main(int argc, char** argv) try gridVariables->init(x); // [[/codeblock]] - // ### Initialize VTK output - using VtkOutputFields = GetPropType; - using VtkWriter = PoreNetwork::VtkOutputModule, SolutionVector>; - VtkWriter vtkWriter(*gridVariables, x, problem->name()); - VtkOutputFields::initOutputModule(vtkWriter); - vtkWriter.addField(gridGeometry->poreVolume(), "poreVolume", Vtk::FieldType::vertex); - vtkWriter.addField(gridGeometry->throatShapeFactor(), "throatShapeFactor", Vtk::FieldType::element); - vtkWriter.addField(gridGeometry->throatCrossSectionalArea(), "throatCrossSectionalArea", Vtk::FieldType::element); - // ### Instantiate the solver - // We use the `LinearPDESolver` class, which is instantiated on the basis + // We use the `NewtonSolver` class, which is instantiated on the basis // of an assembler and a linear solver. When the `solve` function of the - // `LinearPDESolver` is called, it uses the assembler and linear - // solver classes to assemble and solve the linear system around the provided - // solution and stores the result therein. + // `NewtonSolver` is called, it uses the assembler and linear + // solver classes to assemble and solve the non-linear system. // [[codeblock]] - using Assembler = FVAssembler; + using Assembler = FVAssembler; auto assembler = std::make_shared(problem, gridGeometry, gridVariables); using LinearSolver = UMFPackBackend; auto linearSolver = std::make_shared(); - LinearPDESolver solver(assembler, linearSolver); - solver.setVerbose(false); // suppress output during solve() - // [[/codeblock]] - // ### Prepare the upscaling procedure. - // Specify the directions for which the permeability shall be determined (default: x, y, z for 3D). - // [[codeblock]] - using GridView = typename GetPropType::GridView; - const auto defaultDirections = GridView::dimensionworld == 3 ? std::vector{0, 1, 2} - : std::vector{0, 1}; - const auto directions = getParam>("Problem.Directions", defaultDirections); + using NewtonSolver = NewtonSolver; + NewtonSolver nonLinearSolver(assembler, linearSolver); // [[/codeblock]] + // ### Initialize VTK output + using VtkOutputFields = GetPropType; + using VtkWriter = PoreNetwork::VtkOutputModule, SolutionVector>; + VtkWriter vtkWriter(*gridVariables, x, problem->name()); + VtkOutputFields::initOutputModule(vtkWriter); + // specify the field type explicitly since it may not be possible + // to deduce this from the vector size in a pore network + vtkWriter.addField(gridGeometry->poreVolume(), "poreVolume", Vtk::FieldType::vertex); + vtkWriter.addField(gridGeometry->throatShapeFactor(), "throatShapeFactor", Vtk::FieldType::element); + vtkWriter.addField(gridGeometry->throatCrossSectionalArea(), "throatCrossSectionalArea", Vtk::FieldType::element); + + // ### Prepare the upscaling procedure. // Set up a helper class to determine the total mass flux leaving the network const auto boundaryFlux = PoreNetwork::BoundaryFlux(*gridVariables, assembler->localResidual(), x); + // ### The actual upscaling procedure + // #### Instantiate the upscaling helper + // [[codeblock]] + using Scalar = GetPropType; + using UpscalingHelper = UpscalingHelper; + UpscalingHelper upscalingHelper; + // [[/codeblock]] + // Set the side lengths used for applying the pressure gradient and calculating the REV outflow area. // One can either specify these values manually (usually more accurate) or let the UpscalingHelper struct // determine it automatically based on the network's bounding box. @@ -146,78 +143,101 @@ int main(int argc, char** argv) try if (hasParam("Problem.SideLength")) return getParam("Problem.SideLength"); else - return UpscalingHelper::getSideLengths(*gridGeometry); + return upscalingHelper.getSideLengths(*gridGeometry); }(); // pass the side lengths to the problem problem->setSideLengths(sideLengths); // [[/codeblock]] - // ### The actual upscaling procedure + // Get the maximum and minimum pressure gradient and the population of sample points specified in the input file + // [[codeblock]] + const Scalar minPressureGradient = getParam("Problem.MinimumPressureGradient", 1e1); + const Scalar maxPressureGradient = getParam("Problem.MaximumPressureGradient", 1e10); + + if (!(minPressureGradient < maxPressureGradient)) + DUNE_THROW(Dune::InvalidStateException, "Maximum pressure gradient must be greater than minimum pressure gradient"); + + const int numberOfSamples = getParam("Problem.NumberOfPressureGradients", 1); + // [[/codeblock]] - // Iterate over all directions specified before, apply the pressure gradient, calculated the mass flux - // and finally determine the permeability. + // Iterate over all directions specified before, apply several pressure gradient, calculated the mass flux + // and finally determine the the upscaled properties. // [[codeblock]] + const auto directions = getParam>("Problem.Directions", std::vector{0, 1, 2}); + upscalingHelper.setDirections(directions); for (int dimIdx : directions) { - // reset the solution - x = 0; - // set the direction in which the pressure gradient will be applied problem->setDirection(dimIdx); - // solve problem - solver.solve(x); + for (int i = 0; i < numberOfSamples; i++) + { + // reset the solution + x = 0; - // write the vtu file for the given direction - vtkWriter.write(dimIdx); + // set the pressure gradient to be applied + Scalar pressureGradient = maxPressureGradient * std::exp(i + 1 - numberOfSamples); + if (i == 0) + pressureGradient = std::min(minPressureGradient, pressureGradient); - // get the Scalar type - using Scalar = GetPropType; + problem->setPressureGradient(pressureGradient); - // calculate the permeability - const Scalar totalFluidMassFlux = boundaryFlux.getFlux(std::vector{problem->outletPoreLabel()})[0]; - const Scalar K = UpscalingHelper::getDarcyPermeability(*problem, totalFluidMassFlux); + // solve problem + nonLinearSolver.solve(x); - // optionally compare with reference data - static const auto referenceData = getParam>("Problem.ReferenceData", std::vector{}); - if (!referenceData.empty()) - { - static const Scalar eps = getParam("Problem.TestEpsilon"); - if (Dune::FloatCmp::ne(K, referenceData[dimIdx], eps)) - { - std::cerr << "Calculated permeability of " << K << " does not match with reference value of " << referenceData[dimIdx] << std::endl; - return 1; - } + // set the sample points + const Scalar totalFluidMassFlux = boundaryFlux.getFlux(std::vector{ problem->outletPoreLabel() })[0]; + upscalingHelper.setDataPoints(*problem, totalFluidMassFlux); } + + // write a vtu file for the given direction for the last sample + vtkWriter.write(dimIdx); } + // calculate and report the upscaled properties + constexpr bool isCreepingFlow = std::is_same_v; + upscalingHelper.calculateUpscaledProperties(*problem, isCreepingFlow); + upscalingHelper.report(isCreepingFlow); + + // compare the Darcy permeability with reference data if provided in input file and report in case of inconsistency + static const auto referenceData = getParam>("Problem.ReferencePermeability", std::vector{}); + if (!referenceData.empty()) + upscalingHelper.compareWithReference(referenceData); + + // plot the results just for non-creeping flow + // creeping flow would just result in a straight line (permeability is independent of the pressure gradient) + if (!isCreepingFlow) + upscalingHelper.plot(); +}; + +} // end namespace Dumux +// [[/codeblock]] + +// ### The main function +// [[details]] main +// [[codeblock]] +int main(int argc, char** argv) +{ + using namespace Dumux; + + // We parse the command line arguments. + Parameters::init(argc, argv); + + // Convenience alias for the type tag of the problem. + using CreepingFlowTypeTag = Properties::TTag::PNMUpscalingCreepingFlow; + using NonCreepingFlowTypeTag = Properties::TTag::PNMUpscalingNonCreepingFlow; + // // [[/codeblock]] + + // user decides whether creeping flow or non-creeping flow should be run + if (getParam("Problem.AssumeCreepingFlow", false)) + runExample(); + else + runExample(); + // program end, return with 0 exit code (success) return 0; } // [[/codeblock]] -// ### Exception handling -// In this part of the main file we catch and print possible exceptions that could -// occur during the simulation. -// [[details]] error handler -catch (const Dumux::ParameterException &e) -{ - std::cerr << std::endl << e << " ---> Abort!" << std::endl; - return 1; -} -catch (const Dune::DGFException & e) -{ - std::cerr << "DGF exception thrown (" << e << - "). Most likely, the DGF file name is wrong " - "or the DGF file is corrupted, " - "e.g. missing hash at end of file or wrong number (dimensions) of entries." - << " ---> Abort!" << std::endl; - return 2; -} -catch (const Dune::Exception &e) -{ - std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl; - return 3; -} // [[/details]] // [[/content]] diff --git a/examples/porenetwork_upscaling/params.input b/examples/porenetwork_upscaling/params.input index 186151998cf86fe93a341486a17d4db161a4be3d..12261739b5861cb6094b9bc2d780c6dfd5d6e38c 100644 --- a/examples/porenetwork_upscaling/params.input +++ b/examples/porenetwork_upscaling/params.input @@ -1,9 +1,12 @@ [Problem] EnableGravity = 0 # disable gravity Name = upscaling_pnm -PressureGradient = 10 -Directions = 0 1 2 -Epsilon = 1e-3 +MaximumPressureGradient = 1e10 #Pa/m +NumberOfPressureGradients = 10 # the number of sample points for regression process +Directions = 0 # the directions to be plotted 0 1 2 +Epsilon = 1e-13 +AssumeCreepingFlow = false +ReferencePermeability = 3.326e-12 1.394e-12 7.449e-13 [Grid] UpperRight = 4e-4 4e-4 4e-4 @@ -18,7 +21,7 @@ ParameterRandomNumberSeed = 1 BoundaryPoreLabels = xMin:1 xMax:2 yMin:3 yMax:4 zMin:5 zMax:6 MinThroatLength = 1e-10 -DeletionProbability = 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 +DeletionProbability = 0.0 0.6 0.98 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 RemoveThroatsOnBoundary = 0 1 2 3 4 5 CapPoresOnBoundaries = 0 1 2 3 4 5 Sanitize = true @@ -27,3 +30,6 @@ DeletionRandomNumberSeed = 33 [Component] LiquidKinematicViscosity = 1e-6 LiquidDensity = 1e3 + +[Assembly.NumericDifference] +PriVarMagnitude = 1e5 diff --git a/examples/porenetwork_upscaling/problem.hh b/examples/porenetwork_upscaling/problem.hh index 2a4203195b5163c02ada076b33172fdb3f8e273c..2b028111c5525f47bf1cc21f7170dd59699260a4 100644 --- a/examples/porenetwork_upscaling/problem.hh +++ b/examples/porenetwork_upscaling/problem.hh @@ -39,43 +39,55 @@ template class UpscalingProblem : public PorousMediumFlowProblem { // [[details]] convenience aliases + // [[codeblock]] using ParentType = PorousMediumFlowProblem; - using GridGeometry = GetPropType; - using FVElementGeometry = typename GetPropType::LocalView; - using SubControlVolume = typename FVElementGeometry::SubControlVolume; using Scalar = GetPropType; + using GridGeometry = GetPropType; + using FVElementGeometry = typename GridGeometry::LocalView; + using SubControlVolume = typename GridGeometry::SubControlVolume; + using Element = typename GridGeometry::GridView::template Codim<0>::Entity; + using GlobalPosition = typename Element::Geometry::GlobalCoordinate; using PrimaryVariables = GetPropType; using Indices = typename GetPropType::Indices; using BoundaryTypes = Dumux::BoundaryTypes; - using Element = typename GridGeometry::GridView::template Codim<0>::Entity; - using GlobalPosition = typename Element::Geometry::GlobalCoordinate; + // [[/codeblock]] // [[/details]] - - // ### The constructor of our problem. + // + // #### The constructor of our problem. // [[codeblock]] public: - template - UpscalingProblem(std::shared_ptr gridGeometry, std::shared_ptr spatialParams) - : ParentType(gridGeometry, spatialParams) + UpscalingProblem(std::shared_ptr gridGeometry) + : ParentType(gridGeometry) { // the applied pressure gradient - pressureGradient_ = getParam("Problem.PressureGradient"); + pressureGradient_ = getParam("Problem.PressureGradient", 1e5); // We can either use pore labels (given in the grid file) to identify inlet and outlet pores // or use the network's bounding box to find these pores automatically. Using labels is usually much // more accurate, so this is the default here. useLabels_ = getParam("Problem.UseLabels", true); - // an epsilon value for the bounding box approach + // an epsilon value for the floating point comparisons to determine inlet/outlet pores eps_ = getParam("Problem.Epsilon", 1e-7); } // [[/codeblock]] + // #### Temperature + // We need to specify a constant temperature for our isothermal problem. + // Fluid properties that depend on temperature will be calculated with this value. + // [[codeblock]] + Scalar temperature() const + { return 283.15; } + // [[/codeblock]] + + // Set the pressure gradient to be applied to the network + void setPressureGradient(Scalar pressureGradient) + { pressureGradient_ = pressureGradient; } + // #### Boundary conditions // This function is used to define the __type of boundary conditions__ used depending on the location. - // Here, we use Dirichlet boundary conditions (fixed pressures) at the inlet and outlet and Neumann - // boundary conditions at all remaining boundaries. Note that the PNM only supports Neumann no-flow boundaries. - // The specify a certain mass flux, we would have to use a source term on the boundary pores (which is not done in this example). + // Here, we use Dirichlet boundary conditions (fixed pressures) at the inlet and outlet. Note that the PNM does not support Neumann boundaries. + // To specify a certain mass flux on a boundary, we would have to use a source term on the boundary pores (which is not done in this example). // [[codeblock]] BoundaryTypes boundaryTypes(const Element &element, const SubControlVolume& scv) const { @@ -84,8 +96,6 @@ public: // fix the pressure at the inlet and outlet pores if (isInletPore_(scv)|| isOutletPore_(scv)) bcTypes.setAllDirichlet(); - else // Neumann (no-flow) for the remaining boundaries - bcTypes.setAllNeumann(); return bcTypes; } @@ -147,8 +157,9 @@ public: // Return the applied pressure gradient. Scalar pressureGradient() const { return pressureGradient_; } - - + // [[/codeblock]] + // [[/details]] + // // Return the label of inlet pores assuming a previously set direction. int inletPoreLabel() const { @@ -187,12 +198,8 @@ private: int direction_; GlobalPosition length_; bool useLabels_; - - // [[/codeblock]] - // [[/details]] }; } // end namespace Dumux -// [[/codeblock]] // [[/content]] #endif diff --git a/examples/porenetwork_upscaling/properties.hh b/examples/porenetwork_upscaling/properties.hh index 0359a2a1603c31ed35dfd60ce25abc66a545d873..0aa7ad8fe19808cedd855c1651c95477232b5199 100644 --- a/examples/porenetwork_upscaling/properties.hh +++ b/examples/porenetwork_upscaling/properties.hh @@ -34,6 +34,13 @@ // type tag, which we want to modify or for which no meaningful default can be set. #include // for `TTag::PNMOneP` +// The class that contains a collection of single-phase flow throat transmissibilities +// among them the transmisibility model to be used can be specified in AdvectionType class +#include + +// The class that provides specializations for both creeping and non-creeping advection types. +#include + // The local residual for incompressible flow is included. // The one-phase flow model (included above) uses a default implementation of the // local residual for single-phase flow. However, in this example we are using an @@ -46,18 +53,19 @@ #include #include - // The classes that define the problem and parameters used in this simulation #include "problem.hh" #include "spatialparams.hh" // [[/details]] // // ### `TypeTag` definition -// A `TypeTag` for our simulation is defined, which inherits properties from the -// single-phase flow model and the box scheme. +// Two `TypeTag` for our simulation are defined, one for creeping flow and another for non-creeping flow, +// which inherit properties from the single-phase pore network model. The non-creeping flow inherits +// all properties from the creeping flow simulation but sets an own property for the `AdvectionType`. namespace Dumux::Properties { namespace TTag { -struct PNMUpscaling { using InheritsFrom = std::tuple; }; +struct PNMUpscalingCreepingFlow { using InheritsFrom = std::tuple; }; +struct PNMUpscalingNonCreepingFlow { using InheritsFrom = std::tuple; }; } // ### Property specializations @@ -67,43 +75,49 @@ struct PNMUpscaling { using InheritsFrom = std::tuple; }; // [[codeblock]] // We use `dune-foamgrid`, which is especially tailored for 1D networks. template -struct Grid +struct Grid { using type = Dune::FoamGrid<1, 3>; }; // The problem class specifying initial and boundary conditions: template -struct Problem +struct Problem { using type = UpscalingProblem; }; -//! The spatial parameters to be employed. +//! The spatial parameters template -struct SpatialParams +struct SpatialParams { -private: using GridGeometry = GetPropType; using Scalar = GetPropType; public: using type = PoreNetwork::UpscalingSpatialParams; }; -//! The advection type. +//! The advection type for creeping flow template -struct AdvectionType +struct AdvectionType { -private: using Scalar = GetPropType; using TransmissibilityLaw = PoreNetwork::TransmissibilityPatzekSilin; public: using type = PoreNetwork::CreepingFlow; }; -// We use a single liquid phase consisting of a component with constant fluid properties. +//! The advection type for non-creeping flow (includes model for inertia effects) template -struct FluidSystem +struct AdvectionType { -private: using Scalar = GetPropType; + using TransmissibilityLaw = PoreNetwork::TransmissibilityPatzekSilin; public: + using type = PoreNetwork::NonCreepingFlow; +}; + +// We use a single liquid phase consisting of a component with constant fluid properties. +template +struct FluidSystem +{ + using Scalar = GetPropType; using type = FluidSystems::OnePLiquid >; }; // [[/codeblock]] @@ -111,7 +125,7 @@ public: // Moreover, here we use a local residual specialized for incompressible flow // that contains functionality related to analytic differentiation. template -struct LocalResidual +struct LocalResidual { using type = OnePIncompressibleLocalResidual; }; } // end namespace Dumux::Properties diff --git a/examples/porenetwork_upscaling/spatialparams.hh b/examples/porenetwork_upscaling/spatialparams.hh index a0267a274ad83d2032af2cc1eeeb1876cc320229..4faa56f3458808280bc6fdcf52e3c06eae9d748a 100644 --- a/examples/porenetwork_upscaling/spatialparams.hh +++ b/examples/porenetwork_upscaling/spatialparams.hh @@ -102,11 +102,24 @@ public: return M_PI * r * r; } + // dimensionless kinetic-energy coefficient which for non-creeping flow + template + Scalar kineticEnergyCoefficient(const Element& element, + const SubControlVolume& scv, + const ElementSolutionVector& elemSol) const + { return 1.0; } + + // dimensionless momentum coefficient which for non-creeping flow + template + Scalar momentumCoefficient(const Element& element, + const SubControlVolume& scv, + const ElementSolutionVector& elemSol) const + { return 1.0; } + private: std::vector poreShapeFactor_; }; - } // namespace Dumux::PoreNetwork #endif diff --git a/examples/porenetwork_upscaling/upscalinghelper.hh b/examples/porenetwork_upscaling/upscalinghelper.hh index 6f25be60322ca317846894b458008da52e0177d8..b5beebae82366c1c4cb52fa4da409b03377bac5e 100644 --- a/examples/porenetwork_upscaling/upscalinghelper.hh +++ b/examples/porenetwork_upscaling/upscalinghelper.hh @@ -30,53 +30,109 @@ #include #include #include +#include +#include namespace Dumux { -struct UpscalingHelper +template +class UpscalingHelper { - // ### Calculate the intrinsic permeability +public: + // ### Set data points to calculate intrinsic permeability and Forchheimer coefficient // This function first evaluates the mass flux leaving the network in the direction of the applied pressure gradient. - // Afterwards, the mass flux is converted into an area specify volume flux from which finally the intrinsic Darcy - // permeability $`\mathbf{K}`$ [m$`^2`$] can be evaluated. + // Afterwards, the mass flux is converted into an volume flux which is used to calculate the apparent velocity. + // Then apparent permeability of the network is computed and stored for furthure calculations. // [[codeblock]] - template - static Scalar getDarcyPermeability(const Problem& problem, const Scalar totalMassFlux) + template + void setDataPoints(const Problem &problem, const Scalar totalMassFlux) { // get the domain side lengths from the problem auto sideLengths = problem.sideLengths(); - // create temporary stringstream with fixed scientific formatting without affecting std::cout - std::ostream tmp(std::cout.rdbuf()); - tmp << std::fixed << std::scientific; - static constexpr char dirNames[] = "xyz"; + + // get the applied pressure gradient + const auto pressureGradient = problem.pressureGradient(); + const auto pressureDrop = pressureGradient * sideLengths[problem.direction()]; + + // get the fluid properties + const auto liquidDensity = problem.liquidDensity(); + const auto liquidDynamicViscosity = problem.liquidDynamicViscosity(); // convert mass to volume flux - const auto volumeFlux = totalMassFlux / problem.liquidDensity();; + const auto volumeFlux = totalMassFlux / liquidDensity; + ; + // calculate apparent velocity sideLengths[problem.direction()] = 1.0; const auto outflowArea = std::accumulate(sideLengths.begin(), sideLengths.end(), 1.0, std::multiplies()); - const auto vDarcy = volumeFlux / outflowArea; - const auto K = vDarcy / problem.pressureGradient() * problem.liquidDynamicViscosity(); - tmp << "\n########################################\n" << std::endl; - tmp << dirNames[problem.direction()] << "-direction"; - tmp << ": Area = " << outflowArea << " m^2"; - tmp << "; Massflux = " << totalMassFlux << " kg/s"; - tmp << "; v_Darcy = " << vDarcy << " m/s"; - tmp << "; K = " << K << " m^2" << std::endl; - tmp << "\n########################################\n" << std::endl; - - return K; + const auto vApparent = volumeFlux / outflowArea; + + // compute apparent permeability + const auto KApparent = vApparent / pressureGradient * liquidDynamicViscosity; + // calculate rho v / mu, called inertia to viscous ratio in the rest of the code + const auto inertiaToViscousRatio = liquidDensity * vApparent / liquidDynamicViscosity; + + // store the required data for further calculations + totalPressureDrop_[problem.direction()].push_back(pressureDrop); + apparentVelocity_[problem.direction()].push_back(vApparent); + apparentPermeability_[problem.direction()].push_back(KApparent); + inertiaToViscousRatio_[problem.direction()].push_back(inertiaToViscousRatio); + } + // [[/codeblock]] + + // ### Calculate intrinsic permeability and Forchheimer coefficient. + // This function first calculate intrinsic permeability and Forchheimer coefficient using linear least squares regression method + // and reports them. It also plot the apparent permeability of the porous medium versus Forchheimer number/pressure gradient in each + // simulation. + // [[codeblock]] + template + void calculateUpscaledProperties(const Problem &problem, bool isCreepingFlow) + { + const auto sideLengths = problem.sideLengths(); + const auto liquidDynamicViscosity = problem.liquidDynamicViscosity(); + + for (const auto dirIdx : directions_) + { + // determine Darcy permeability as the maximum permeability of the domain + darcyPermeability_[dirIdx] = *max_element(apparentPermeability_[dirIdx].begin(), apparentPermeability_[dirIdx].end()); + if (!isCreepingFlow) + { + for (int i = 0; i < totalPressureDrop_[dirIdx].size(); i++) + { + // calculate the Darcy pressure drop. + const Scalar darcyPressureDrop = liquidDynamicViscosity * apparentVelocity_[dirIdx][i] * sideLengths[dirIdx] / darcyPermeability_[dirIdx]; + + // calculate the ratio of Dracy to total pressure drop + const Scalar pressureDropRatio = darcyPressureDrop / totalPressureDrop_[dirIdx][i]; + + // set sample points for upscaling of Forchheimer parameters. + // first, check the permability ratio to see if the flow regime is Forchheimer. + if (pressureDropRatio < 0.99) + { + samplePointsX_[dirIdx].push_back(inertiaToViscousRatio_[dirIdx][i]); + samplePointsY_[dirIdx].push_back(1 / apparentPermeability_[dirIdx][i]); + } + } + // determine regression line and accordingly the Forchheimer permeability and the Forchheimer coefficient + const auto [intercept, slope] = linearRegression(samplePointsX_[dirIdx], samplePointsY_[dirIdx]); + forchheimerPermeability_[dirIdx] = 1.0 / intercept; + forchheimerCoefficient_[dirIdx] = slope; + writePlotDataToFile(dirIdx); + } + } } // [[/codeblock]] - // ### Determine the domain's side lengths automatically based on the bounding box of the network. + // ### Determine the domain's side lengths + // + // We determine the domain side length by using the bounding box of the network // [[codeblock]] template - static auto getSideLengths(const GridGeometry& gridGeometry) + auto getSideLengths(const GridGeometry& gridGeometry) { using GlobalPosition = typename GridGeometry::GlobalCoordinate; - GlobalPosition result; + GlobalPosition result(0.0); std::cout << "Automatically determining side lengths of REV based on bounding box of pore network" << std::endl; for (int dimIdx = 0; dimIdx < GridGeometry::GridView::dimensionworld; ++dimIdx) @@ -85,8 +141,212 @@ struct UpscalingHelper return result; } // [[/codeblock]] + + // ### Plot the data using Gnuplot + // + // [[codeblock]] + void plot() + { + // plot permeability ratio vs. Forchheimer number + plotPermeabilityratioVsForchheimerNumber_(); + + // plot inverse of apparent permability vs. rho v / mu + plotInversePrmeabilityVsInertiaToViscousRatio_(); + } + // [[/codeblock]] + // + // ### Save the relevant data for plot + // [[codeblock]] + void writePlotDataToFile(std::size_t dirIdx) + { + // write permeability ratio vs. Forchheimer number + writePermeabilityratioVsForchheimerNumber_(dirIdx); + + // write inverse of apparent permability vs. rho v / mu + writeInversePrmeabilityVsInertiaToViscousRatio_(dirIdx); + } + // [[/codeblock]] + // + // ### Report the upscaled data + // [[codeblock]] + void report(bool isCreepingFlow) + { + // Report the results for each direction + for (const auto dirIdx : directions_) + { + std::cout << Fmt::format("\n{:#>{}}\n\n", "", 40) + << Fmt::format("{}-direction:\n", dirName_[dirIdx]) + << Fmt::format("-- Darcy (intrinsic) permeability = {:.3e} m^2\n", darcyPermeability_[dirIdx]); + + // Report non-creeping flow upscaled properties + if (!isCreepingFlow) + { + std::cout << Fmt::format("-- Forchheimer permeability = {:.3e} m^2\n", forchheimerPermeability_[dirIdx]); + std::cout << Fmt::format("-- Forchheimer coefficient = {:.3e} m^-1\n", forchheimerCoefficient_[dirIdx]); + } + + std::cout << Fmt::format("\n{:#>{}}\n", "", 40) << std::endl; + } + } + // [[/codeblock]] + // + // ### Compare with reference data provided in input file + // [[codeblock]] + void compareWithReference(std::vector referenceData) + { + for (const auto dirIdx : directions_) + { + const auto K = darcyPermeability_[dirIdx]; + static const Scalar eps = getParam("Problem.TestEpsilon", 1e-3); + if (Dune::FloatCmp::ne(K, referenceData[dirIdx], eps)) + { + std::cerr << "Calculated permeability of " << K << " in " + < directions) + { + directions_ = directions; + } + // [[codeblock]] + +private: + // ### Save the relevant data for plot of permeability ratio vs. Forchheimer number + // [[codeblock]] + void writePermeabilityratioVsForchheimerNumber_(std::size_t dirIdx) + { + // open a logfile + std::ofstream logfile(dirName_[dirIdx] + "-dir-PermeabilityratioVsForchheimerNumber.dat"); + + // save the data needed to be plotted in logfile + for (int i = 0; i < apparentPermeability_[dirIdx].size(); i++) + { + // compute the Forchheimer number + const Scalar forchheimerNumber = darcyPermeability_[dirIdx] * forchheimerCoefficient_[dirIdx] * inertiaToViscousRatio_[dirIdx][i]; + // ratio between apparrent permeability and darcy permeability + const Scalar permeabilityRatio = apparentPermeability_[dirIdx][i] / darcyPermeability_[dirIdx]; + + logfile << forchheimerNumber << " " << permeabilityRatio << std::endl; + } + } + // [[/codeblock]] + + // ### Save the relevant data for plot of inverse of apparent permability vs. rho v / mu + // [[codeblock]] + void writeInversePrmeabilityVsInertiaToViscousRatio_(std::size_t dirIdx) + { + // open a logfile and write inverese of apparent permeability given by the model vs. inertial to viscous ratio (rho v / mu) + std::ofstream logfile(dirName_[dirIdx] + "-dir-InversePrmeabilityVsInertiaToViscousRatio.dat"); + + // save the data needed to be plotted in logfile + for (int i = 0; i < apparentPermeability_[dirIdx].size(); i++) + { + const Scalar inertiaToViscousRatio = inertiaToViscousRatio_[dirIdx][i]; + const Scalar inverseAppPermeability = 1 / apparentPermeability_[dirIdx][i]; + + // compute inverse of apparent permeability using the Forchheimer permeability and coefficient + const Scalar inverseAppPermeabilityForchheimer = 1 / forchheimerPermeability_[dirIdx] + inertiaToViscousRatio * forchheimerCoefficient_[dirIdx]; + + logfile << inertiaToViscousRatio << " " << 1e-12 * inverseAppPermeability << " " << 1e-12 * inverseAppPermeabilityForchheimer << std::endl; + } + } + // [[/codeblock]] + // + // ### Plot permeability ratio vs. Forchheimer number using Gnuplot + // + // [[codeblock]] + void plotPermeabilityratioVsForchheimerNumber_() + { + // using gnuplot interface + Dumux::GnuplotInterface gnuplot(true); + gnuplot.setOpenPlotWindow(true); + + for (const auto dirIdx : directions_) + { + gnuplot.resetAll(); + std::string title{}, option{}; + + // add the data in each direction for plot + gnuplot.addFileToPlot(dirName_[dirIdx] + "-dir-PermeabilityratioVsForchheimerNumber.dat", "notitle with lines"); + // set the properties of lines to be plotted + option += "set linetype 1 linecolor 1 linewidth 7\n"; + // report the darcy permeability in each direction as the title of the plot + title += Fmt::format("{}-direction, Darcy permeability= {:.3e} m^2 ", dirName_[dirIdx], darcyPermeability_[dirIdx]); + + option += "set title \"" + title + "\"\n"; + option += "set logscale x""\n"; + option += "set format x '10^{%L}'""\n"; + + gnuplot.setXlabel("Forchheimer Number [-]"); + gnuplot.setYlabel("Apparent permeability / Darcy permeability [-]"); + gnuplot.setOption(option); + gnuplot.plot("permeability_ratio_versus_forchheimer_number"); + } + } + // [[/codeblock]] + // + // ### Plot inverse of apparent permability vs. rho v / mu using Gnuplot + // + // [[codeblock]] + void plotInversePrmeabilityVsInertiaToViscousRatio_() + { + // using gnuplot interface + Dumux::GnuplotInterface gnuplot(true); + gnuplot.setOpenPlotWindow(true); + + for (const auto dirIdx : directions_) + { + gnuplot.resetAll(); + std::string title{}, option{}; + std::string legend0 = "u 1:2 title \"Network model\" with lines"; + // add the data in each direction for plot, first set of data + gnuplot.addFileToPlot(dirName_[dirIdx] + "-dir-InversePrmeabilityVsInertiaToViscousRatio.dat", legend0); + + // set the properties of lines to be plotted + option += "set linetype 1 linecolor 1 linewidth 5\n"; + + std::string legend1 = "u 1:3 title \"Forchheimer equation\" with lines"; + // add the data in each direction for plot, second set of data + gnuplot.addFileToPlot(dirName_[dirIdx] + "-dir-InversePrmeabilityVsInertiaToViscousRatio.dat", legend1); + + // set the properties of lines to be plotted + option += "set linetype 2 linecolor 2 linewidth 5\n"; + + // report the darcy permeability in each direction as the title of the plot + title += Fmt::format("{}-direction, Darcy permeability= {:.3e} m^2 ", dirName_[dirIdx], darcyPermeability_[dirIdx]); + + option += "set title \"" + title + "\"\n"; + option += "set logscale x""\n"; + option += "set format x '10^{%L}'""\n"; + + gnuplot.setXlabel("{/Symbol r} v / {/Symbol m} [1/m]"); + gnuplot.setYlabel("1/ Apparent permeability [1/m^2] x 1e12"); + gnuplot.setOption(option); + gnuplot.plot("inverse_apppermeability_versus_rhovmu-1"); + } + } + // [[/codeblock]] + std::array, 3> samplePointsX_; + std::array, 3> samplePointsY_; + std::array, 3>totalPressureDrop_; + std::array, 3> apparentVelocity_; + std::array, 3> apparentPermeability_; + std::array, 3> inertiaToViscousRatio_; + std::array darcyPermeability_; + std::array forchheimerPermeability_; + std::array forchheimerCoefficient_; + const std::array dirName_ = {"X", "Y", "Z"}; + std::vector directions_; }; } // end namespace Dumux +// [[/codeblock]] // [[/content]] #endif diff --git a/test/references/example_porenetwork_upscaling_noncreeping_flow_X-dir_InversePrmeabilityVsInertiaToViscousRatio.dat b/test/references/example_porenetwork_upscaling_noncreeping_flow_X-dir_InversePrmeabilityVsInertiaToViscousRatio.dat new file mode 100644 index 0000000000000000000000000000000000000000..b4d27680fa14119188aac2e3244b99a1410f9e38 --- /dev/null +++ b/test/references/example_porenetwork_upscaling_noncreeping_flow_X-dir_InversePrmeabilityVsInertiaToViscousRatio.dat @@ -0,0 +1,10 @@ +0.0332584 0.300676 0.31292 +11107 0.302029 0.31387 +29965.6 0.30431 0.315483 +79894.3 0.310254 0.319753 +207384 0.324902 0.330657 +512009 0.357721 0.356711 +1.17636e+06 0.423229 0.413531 +2.49648e+06 0.542104 0.526439 +4.93037e+06 0.74615 0.734604 +9.19349e+06 1.08773 1.09922 diff --git a/test/references/example_porenetwork_upscaling_noncreeping_flow_X-dir_PermeabilityratioVsForchheimerNumber.dat b/test/references/example_porenetwork_upscaling_noncreeping_flow_X-dir_PermeabilityratioVsForchheimerNumber.dat new file mode 100644 index 0000000000000000000000000000000000000000..e071f632ecafef7e6ad4621bdaaf907e47b97b29 --- /dev/null +++ b/test/references/example_porenetwork_upscaling_noncreeping_flow_X-dir_PermeabilityratioVsForchheimerNumber.dat @@ -0,0 +1,10 @@ +9.46045e-09 1 +0.00315941 0.995519 +0.00852378 0.988056 +0.0227261 0.969127 +0.0589909 0.925434 +0.145642 0.84053 +0.334619 0.710432 +0.710131 0.554645 +1.40246 0.402969 +2.61511 0.276426