Coupling staggered-grid and vertex-centered finite-volume methods for coupled porous-medium free-flow problems

by   Martin Schneider, et al.
University of Stuttgart

In this work, a new discretization approach for coupled free and porous-medium flows is introduced, which uses a finite volume staggered-grid method for the discretization of the Navier-Stokes equations in the free-flow subdomain, while a vertex-centered finite volume method is used in the porous-medium flow domain. The latter allows for the use of unstructured grids in the porous-medium subdomain, and the presented method is capable of handling non-matching grids at the interface. In addition, the accurate evaluation of coupling terms and of additional nonlinear velocity-dependent terms in the porous medium is ensured by the use of basis functions and by having degrees of freedom naturally located at the interface. The available advantages of this coupling method are investigated in a series of tests: a convergence test for various grid types, an evaluation of the implementation of coupling conditions, and an example using the velocity dependent Forchheimer term in the porous-medium subdomain.



There are no comments yet.


page 11

page 14

page 15

page 16

page 18

page 19


Analysis of the Stokes-Darcy problem with generalised interface conditions

Fluid flows in coupled systems consisting of a free-flow region and the ...

Effective coupling conditions for arbitrary flows in Stokes-Darcy systems

Boundary conditions at the interface between the free-flow region and th...

A global-in-time domain decomposition method for the coupled nonlinear Stokes and Darcy flows

We study a decoupling iterative algorithm based on domain decomposition ...

An asynchronous incomplete block LU preconditioner for computational fluid dynamics on unstructured grids

We present a study of the effectiveness of asynchronous incomplete LU fa...

Uncertainty-aware Validation Benchmarks for Coupling Free Flow and Porous-Medium Flow

A correct choice of interface conditions and useful model parameters for...

Numerical implementation of efficient grid-free integral wall models in unstructured-grid LES solvers

Two zonal wall-models based on integral form of the boundary layer diffe...

Comparison of cell- and vertex-centered finite-volume schemes for flow in fractured porous media

Flow in fractured porous media is of high relevance in a variety of geot...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Descriptions of flow and transport processes in our natural, biological, and engineered environment exist for a diverse spectrum of flow conditions. Isolated from each other, these descriptive models can be used to effectively simulate flow within their domain of analysis. On the contrary, most relevant applications contain multiple coupled flow regimes. One common example of multi-domain systems that are challenging to resolve with an isolated flow model are coupled free-flow and multi-phase porous-medium flows. Mass, momentum, and heat exchange will occur across the interface between these domains dependent on the conditions in each subdomain. The flow fields on either side of the interface will change with this exchanged content, further driving the exchange processes between them.

Applications that rely on a detailed description of this exchange are abundant, ranging from large scale environmental applications such as soil-water evaporation and salt precipitation [37, 23], to the small scale exchange processes within a PEM fuel cell [16]. Evaluating these applications via direct numerical simulation of the flow fields, like those shown in [8], is usually infeasible due to the immense computational cost and the effort involved in procuring pore geometries. For the analysis of most applications, evaluations using the concept of a representative elementary volume (REV) scale [2] can be used, where averaging techniques create upscaled descriptions of flow in a porous domain [39].

Evaluating coupled free-flow and porous-medium flow systems at this REV scale can then be done using a single- or two-domain approach. A single-domain approach evaluates both domains with a single set of equations, and is possible using the Brinkman equations [6]; this method has been expanded upon significantly in [7]. With a two-domain approach, momentum conservation in each domain can be described using a model designed for its context, in this case, the Navier-Stokes equations for describing free flow, and the multi-phase Darcy’s law for describing multi-phase flow in a porous medium [28]. In order to couple these two flow domains, thermodynamically consistent interface conditions have to be used to enforce mass, momentum, and energy exchanges between these two flow domains [19].

Numerous discretization methods for these coupled two-domain models have been developed. Some models, e.g. [14, 33], couple cell-centered finite volumes, used to discretize the porous medium, with a Marker-and-Cell (staggered grid) finite volume scheme [18], used for the discretization of the free-flow domain. Other approaches utilize the same discretization method in each domain, using either staggered grids [22], finite elements [11], or a vertex-centered finite volume method, otherwise known as the box scheme [28, 21] for both domains. Other models have been developed to use an extra lower-dimensional grid, or interface mortar elements, at the interface where the exchange fluxes can be calculated [1, 5].

For many porous media based applications, variations of the core mathematical models have been developed to describe specific physical phenomena by considering additional terms in the balance equations. These terms include the Forchheimer term, describing viscous stresses within a porous medium [29], and dispersion terms, describing transport due to subscale fluctuations in velocity [31]. Both of these terms require a reconstruction of velocity for collocated discretization methods, typically requiring an expansion of local stencils [34]. Vertex-centered finite volume methods, such as the box scheme, show their strength in this reconstruction, where the basis functions can be used to determine the velocities used in these terms [21]. In this work, the free flow is modeled using the Navier-Stokes equations discretized with the staggered-grid method. The porous-medium flow is modeled using the multi-phase Darcy’s law including the Forchheimer term, and is discretized using a vertex-centered finite volume method.

The structure of this work is as follows: Section 2 reviews the governing model equations, Section 3 gives an overview of the discretization methods used in both domains, and Section 4 examines the performance of the proposed coupling method by means of numerical test cases. A review of this work, and an outlook towards further goals is provided in Section 5.

2 Governing Equations

The following section provides an outline of the assumptions made regarding geometry while developing this work, and reviews the notation. Furthermore, the governing equations for the free-flow and the porous-medium flow domains are presented, followed by the conditions used to couple them.

The domain in question, , , is an open, connected Lipschitz domain with boundary and -dimensional measure . This domain is then partitioned disjointly to and representing the free-flow and porous-medium flow subdomains, respectively. The boundaries of each subdomain are defined as the interface as well as the remainders and . In the following, the superscripts are omitted if no ambiguity arises.

Using a subscript to denote a prescribed boundary condition for either pressure, , or velocity, , the external boundary is further decomposed such that and disjointly. Further, we assume that for a subset of the boundary , a pressure boundary condition is imposed with positive measure, i.e.

, ensuring that the resulting system can be uniquely solved. The vector

denotes the outwardly oriented unit normal vector on with respect to .

Molecular diffusion and heat conductions are modeled after the following laws;

Fick’s law: (1)
Fourier’s law: (2)

where is the mass density of phase , is the diffusion coefficient of component in phase , and is the mass fraction of component in phase . is the effective thermal conductivity [37], which in is equal to the gas-phase thermal conductivity, and is the temperature, assuming local thermal equilibrium in .

These laws form an essential part of the governing equations in the subdomains and at the interface, which are presented in the subsequent sections.

2.1 Free Flow

Within the coupled model, the free-flow subdomain consists of a gas phase with multiple components . Here, we drop the subscript which relates all quantities to the gas phase. The mass balance (3a) for each component , the energy balance (3b), and the momentum balance equations (3c) are given as: equationparentequation


together with appropriate boundary conditions. The unknown variables are the gas-phase velocity , pressure and temperature , and one of the mass or mole fractions . Here, and denote the potentially solution-dependent gas-phase density and viscosity, respectively, while are source (or sink) terms, and describes the influence of gravity. is an additional momentum sink or source term.

is the identity tensor in

. and denote the enthalpy and internal energy, respectively.

2.2 Porous-Medium Flow

The equations governing compositional non-isothermal two-phase (water and air) flow in the porous-medium subdomain, , are given by: equationparentequation


together with appropriate boundary conditions. is the saturation of phase while is the porosity of the permeable medium. and are the solid-phase density and heat capacity.

Eq. 4c states that the momentum balance in the porous medium is given by Darcy’s law, with being the permeability tensor. Similar to the free-flow equations (3), denotes source or sink terms. The relative permeability, , as well as the capillary pressure, , used as closure relations, are modeled using the van Genuchten approach [36, 27], where and denote the gas and water phase, respectively.

2.3 Coupling Conditions

In order to derive a thermodynamically consistent formulation of the coupled problem, the conservation of mass, energy, and momentum must be guaranteed at the interface between the porous medium and the free-flow domain. We therefore impose the following interface conditions: equationparentequation


The momentum transfer normal to the interface is given by (5e) [26]. Condition (5f) is the commonly used Beavers-Joseph-Saffman slip condition [3, 30]. Here, denotes any unit vector from the tangent bundle of and is a parameter to be determined experimentally. We remark that this condition is technically a boundary condition for the free flow, not a coupling condition between the two flow regimes. Furthermore, it has been developed for free flow strictly parallel to the interface and might lose its validity for other flow configurations [40]. We are nevertheless going to use Eq. 5f in the upcoming numerical examples which also include non-parallel flow at the interface, thus accepting a potential physical inconsistency for the sake of numerical verification. Equation 5f could be replaced by a different condition in future work, such as those outlined in [13, 35], as this condition is not essential to the numerical model.

3 Discretization

This section provides a brief overview of the notation regarding the partition of , before the numerical schemes used in the subdomains and the incorporation of the coupling conditions (5) are outlined.

Definition 1 (Grid discretization).

The tuple denotes the grid discretization, in which

  1. is the set of primary grid elements such that . For each element , and denote its volume and barycenter.

  2. is the set of faces such that each face is a

    -dimensional hyperplane with measure

    . For each primary grid element , is the subset of such that . Furthermore, denotes the barycenter and the unit vector that is normal to and outward to .

  3. is the set of control volumes (dual grid elements) such that . For each control volume , denotes its volume.

  4. is the set of faces such that each face is a -dimensional hyperplane with measure . For each , is the subset of such that . Furthermore, denotes the face evaluation points and the unit vector that is normal to and outward to .

  5. is the set of cell centers such that and is star-shaped with respect to .

  6. is the set of vertices of the grid, corresponding to the corners of the primary grid elements .

For ease of exposition, we assume , however, this is not a limitation and the model can be readily extended to three dimensions. In fact, our implementation is capable of handling three-dimensional systems, which is demonstrated in the section containing the numerical examples. To distinguish geometrical quantities between the individual subdomains, the superscripts and are used, which are omitted if no ambiguity arises.

3.1 Staggered-grid scheme

A staggered-grid finite volume scheme, also known as MAC scheme [18, 33], is used in the free-flow subdomain. Scalar quantities, such as pressure and density, are defined at the cell centers while degrees of freedom corresponding to the components of the velocity vector are located on the primary control volumes’ faces . As opposed to collocated schemes, where all variables are defined at the same location, this resulting scheme is stable, guaranteeing oscillation-free solutions without any additional stabilization techniques [38].

The current implementation is restricted to grid cells that are aligned with the coordinate axes (hyperrectangles) such that it is assumed that is a structured conforming quadrilateral mesh. This is a restriction of this scheme and general interface topologies between the free-flow region and the porous medium cannot be resolved by using unstructured grids. Although extensions to this scheme, allowing general unstructured staggered schemes, are being investigated in ongoing research, these extensions are not within the scope of this paper and will not be discussed. The coupling approach presented in this paper can be seen as the basis for more general coupling concepts.

Equations to fully outline the staggered-grid discretization scheme are listed in [33] and unchanged in this work. For incorporation of boundary conditions, we refer the reader to [38].

3.2 Vertex-centered finite volume scheme

For typical cell-centered finite volume schemes, see [12, 32] for an overview, elements and control volumes usually coincide, . On the contrary, the box scheme [17, 21] introduces control volumes that do not coincide with the primary grid elements, but are defined around the grid vertices. In order to outline their construction, let us introduce with the set of primary grid faces adjacent to the vertex , i.e. for each it is . To continue, let us introduce the set of adjacent primary grid elements , i.e. we have for each . The control volume associated with the vertex of the grid is then constructed by connecting the centers of all adjacent elements and the centers of all adjacent faces . Please note that in three-dimensional settings, the construction furthermore considers the centers of the adjacent grid edges. For an illustration of the control volumes constructed at the interface to the free-flow domain, see Fig. 1.

Each control volume is partitioned into sub-control volumes , and with we denote the set of sub-control volumes embedded in the control volume such that . Each sub-control volume is defined as the overlap of the control volume with an adjacent element , that is, .

In the following, we want to derive the discrete formulation of 4. For the sake of readability, let us consider a general balance equation, representative for Eqs. 4b and 4a, of the form:


Here, is the conserved quantity and is the flux term. Integration of Eq. 6 over a control volume , making use of the divergence theorem, yields


Let us now split the integral over the control volume into integrals over the sub-control volumes, and let us introduce the discrete flux over a face :


Furthermore, let us denote with and the conserved quantity and the source term evaluated for the sub-control volume . Assuming that the time derivative is sufficiently smooth, this yields:


Finally, we approximate the time derivative by a finite difference:


where is the current time level and denotes the time step size. In this work, we use an implicit time discretization scheme, i.e. we evaluate the discrete fluxes and source terms in (10) on the time level . Recall that each control volume can be uniquely associated to a primary grid vertex. As a consequence, the same also holds for the sub-control volumes . A discrete value of is associated to each grid vertex, which we denote with , , and we understand , and to refer to the same discrete value. In this work, we use the mass lumping technique proposed by [21] to improve the stability of the scheme. Thus, we use the nodal values for storage term evaluations:


The above equation is the discrete form of (6) used in the box scheme. In the remainder of this section, the expressions for the discrete fluxes are presented. In the box scheme, fluxes are computed using continuous piecewise linear basis functions on the primary discretization. A basis function is associated to each vertex of the discretization, which means the basis functions are also uniquely associated to the control volumes . Moreover, the basis function associated with the vertex is equal to one at the vertex itself and zero on all other vertices of the grid:


where is the position of the vertex and denotes the Kronecker-Delta. As seen in the above condition, the basis functions have a small support region, i.e. they are only non-zero within the elements adjacent to :


The basis functions further describe a partition of unity, and we can express the discrete approximation of at a position with:


Extending from (14), the expression for the discrete gradient of is:


Let us further introduce , or the set of vertices contained in the primary grid element , which enables us to express the discrete approximation of as well as the gradient at a position by:


Making use of the above definitions, we define the discrete fluxes of the component across a face , which describes a part of the boundary of the control volume , and which is contained in the primary grid element , as follows: equationparentequation


Here, we introduced with the permeability defined for the primary grid element , and with and the phase density and diffusion coefficient associated with the face

, respectively. While the density at the face is defined by the interpolation of the nodal values

, the diffusion coefficient on the face is taken as the harmonic mean of the diffusion coefficients in the two adjacent sub-control volumes.

Similarly, we define the discrete heat fluxes as: equationparentequation


Here, the effective thermal conductivity on the face is again defined as the harmonic mean of the conductivities in the two sub-control volumes adjacent to it. In the discrete fluxes stated in (17b) and (18a) we have used the notation to denote quantities that are evaluated in the upstream or downstream sub-control volume of face , i.e. (used as default value) corresponds to a full upwinding and to an averaging. With the discrete conservation Eq. 11 and the discrete fluxes (17) and (18), we can state the discrete form of the mass balance Eq. 4a for each control volume :


and correspondingly for the energy balance Eq. 4b:


where .

3.3 Coupling

This section is devoted to the incorporation of the coupling conditions 5 in the discrete setting. The discretizations of the sub-domains are allowed to be non-matching at the interface, which results in arbitrary overlap between grid faces of the free-flow and the porous-medium domain, respectively. Let us partition the interface into a set of facets (line segments for the two-dimensional case), for which we note that . Each facet is the intersection between a face of the discretization of the porous-medium domain and a face of the discretization of the free-flow domain, depicted as coupling segments in Fig. 1. Each face of the discretizations of the subdomains overlaps with one or more such facets, which we collect in the set and where we note that .

Figure 1: Control volume of the box method at the interface to the free-flow subdomain. The unknowns of the box method are located at grid vertices, whereas pressures, mass fractions, and temperatures of the staggered scheme are located at cell centers, and velocities at faces. The basis functions of the box method can be used to interpolate solution values at the interface.

The coupling between the two subdomains is realized as follows: Eqs. 5d and 5b represent Neumann (flux) boundary conditions, where only the free-flow terms are calculated and set accordingly as boundary fluxes. This guarantees the conservation of mass and energy. To calculate these fluxes, we use the advantage of the vertex-centered finite volume scheme, where degrees of freedom are naturally located at the interface. In the following, we present the construction of the free-flow fluxes at the interface by using a two-point flux approximation (Tpfa) [32] to discretize diffusive and conductive fluxes given by Fick’s and Fourier’s law (1)-(2). This choice is based on the fact that, in the free-flow domain, we are restricted to grid cells that are aligned with the coordinate axes, as mentioned in Section 3.1. In this case, the Tpfa yields an accurate flux approximation.

Let be a face of some control volume of the free-flow grid (see Fig. 1) such that . Integrating the free-flow related terms of Eqs. 5d and 5b over yields the following discrete approximations:


with flux approximations


which are related to the control volume and the corresponding sub-control volume of the porous-medium domain such that . The interface density is evaluated with the following average, but the diffusion coefficient is evaluated as , using only information from the free-flow control volume. The upwind quantities are analogous to the previous definition and the default is such that only upstream information, i.e. either from the free-flow or the the porous-medium domain, is used.

For the construction of diffusive/conductive fluxes 24 we use the continuity of the primary variables and , which is provided via conditions 5c and 5a. This allows an evaluation of these quantities using information from the discrete porous-medium solution. This is done with two types of projection operators :


can be seen as a simplification of 25 by using a simple quadrature rule, which has the advantage of yielding smaller coupling stencils. By using the same fluxes also in the porous domain, i.e. , the conservation of advective/convective and diffusive/conductive fluxes can directly be observed.

The remaining coupling conditions (5e), for normal momentum transfer, and the slip condition (5f), are needed only for the free-flow subdomain. Equation 5e yields a Neumann boundary condition for the momentum balance by evaluating at each free-flow face. This is again realized by using the projections introduced above. When using the Saffman simplification of the Beavers-Joseph condition, as shown in Equation 5f, the porous-medium velocity is neglected meaning only free-flow unknowns are required and the condition is independent of the chosen porous-medium discretization scheme.

4 Numerical results

The coupling approach presented in the previous section was implemented within the the open-source simulation environment DuMu

[25], which exists as an additional module of the DUNE project [4]. We employ a monolithic approach, where both sub-problems are assembled into one system of equations and use an implicit Euler method for the time discretization. To accommodate for the nonlinearities in the systems of equations, Newton’s method is employed and the direct solver UMFPack [10] is used to solve the resulting linear systems of equations produced in each Newton iteration. The Dune-Subgrid [15] and Dune-UGGrid modules are used to represent the structured grids in the free-flow and the unstructured grids used in the porous-medium subdomain, respectively.

Within this section, a series of three test cases will be used to examine the performance of the proposed coupling approach. Section 4.1 provides a convergence study, while Section 4.2 presents a comparison between the proposed method and the coupling approaches introduced in [14]. Finally, a showcase simulation is given in Section 4.3 to demonstrate the flexibility of the method at application scale.

4.1 Convergence Study

In this section we investigate the convergence behavior of the proposed coupling scheme on a test case that was originally presented in [33]. It assumes incompressible single-phase flow, neglecting non-isothermal, compositional, and gravitational effects. Consider a domain , decomposed into and with interface . The material parameters are given by


where, depending on the parameter , the permeability tensor depends on the spatial coordinate and varies with the angular frequency . Furthermore, let the solution be given by


Neglecting time derivatives, we have in the free-flow domain


and in the porous medium


The validity of the coupling conditions can also be easily verified. In order to guarantee that is positive definite and a full tensor with non-negligible off-diagonal entries (depending on ), we set and . For calculating the fluxes in the porous domain, the permeability tensors are evaluated at the cell centers. The functions , and are numerically integrated using a 5th-order quadrature rule.

We use the following discrete pressure error norm,


where and where the parameter indicates the -th refinement level. Furthermore, and denote the control volume-wise constant discrete and exact pressure values, obtained by evaluating both in the center of the control volumes .

For the staggered-grid scheme, we additionally quantify the errors for the velocity unknowns. These errors are calculated as


where the velocity unknowns are associated with the dual control volumes , which are constructed around the faces (see Section 3). The exact velocity values at the face centers are denoted as .

Three different grid configurations in the porous-medium domain are considered: a grid that is conforming to the grid used in the free-flow domain (fig:convtest_grids_conforming), an unstructured simplex grid with non-matching interface (fig:convtest_grids_nonconforming), and a structured grid for which the dual grid used in the box scheme is conforming to the grid used in the free-flow domain (fig:convtest_grids_boxconforming). These cases are refered to here as the conforming, non-conforming, and box-conforming cases, respectively. Figure 2 shows all configurations after the first grid refinement. Moreover, we consider both projection operators (see Eq. 25) and (see Eq. 26).

Figure 2: Three grid configurations used in the convergence test: (fig:convtest_grids_conforming) conforming, (fig:convtest_grids_nonconforming) non-conforming, and (fig:convtest_grids_boxconforming) box-conforming. All grids are shown after the first refinement.

Figure 3 shows the error norms plotted over grid refinement for all combinations of grids and projection operators. Second order convergence of the considered discrete error norms is observed in all combinations, with only minor differences in the error values. The largest differences appear for , where the smallest errors are observed with on the box-conforming grid configuration, and the largest errors occur on the non-conforming grid configuration. Note that in this case, the errors observed with the projection operator are slightly smaller than with , despite being computationally cheaper. Note also that on the box-conforming grid, the projections and are equivalent, which is why we only show the result for here.

Figure 3: Convergence test. Error norms vs grid refinement.

4.2 Porous obstacle

The purpose of this test is to show that the presented coupling approach is able to accurately capture the interface exchange processes described by the continuous coupling conditions Eq. 5. By comparing different coupling methods, [14] have demonstrated that for an accurate description of the processes at the interface, additional interface unknowns have to be introduced. In their publication, a cell-centered Tpfa scheme on a structured grid has been used in the porous-medium domain coupled to a staggered-grid finite volume scheme used in the free-flow domain. When solving physically complex flow processes, described by nonlinear PDEs, the additional interface unknowns have to be calculated locally within a sub-routine that requires a nonlinear solver. As mentioned before, the advantage of the presented approach is that no interface solver is needed because the numerical schemes allow for an evaluation of all required quantities directly at the interface.

In this section, it is shown that the presented approach yields results that are comparable to those presented in [14] with the coupling method “CM4”, which involves local nonlinear solves for the interface quantities. In order to draw a fair comparison, we replicate a test case that was used in the original study and which is illustrated in Fig. 4. It considers a two-dimensional free-flow channel in which a flow of air is enforced by setting a parabolic velocity profile at the top boundary, while allowing outflow across the bottom boundary. A porous medium is placed in the center of the channel, leaving narrower channels on the left and right for the air to potentially bypass the porous medium. Due to the symmetry of the setup, the test case only considers the left half of the domain while applying symmetry boundary conditions on the right boundary. The left boundary is set to no-flow. For indications on the domain dimensions, see Fig. 4.

Figure 4: Porous obstacle. Setup and boundary conditions for the test case taken from [14].

A two-phase, two-component model is used in the porous medium (see Eq. 4), while a single-phase, two-component model is considered in the free-flow domain (see Eq. 3). Initially, the porous medium is partially saturated with water () and at a temperature of , while the temperature of the air in the free-flow domain is . We then simulate of air flowing through the channel, which leads to a partial desaturation of the porous medium, especially in its upper region where the highest gas velocities occur. Figure 5 illustrates the state at the final simulation time. We observe that the air in the free-flow channel is cooled down by the contact and infiltration of the cooler porous medium. Additional cooling occurs due to the evaporation of water inside the porous medium, which leads to a decrease in temperature in the upper left region of the porous medium, despite the fact that warmer air is flowing through it. The evaporated water is transported away with the gas phase as depicted in the central image of Fig. 5. The right image shows the water saturation in the porous medium, which also shows a slight decrease towards the upper left boundary, caused by some of the water evaporating into the gas phase.

Figure 5: Porous obstacle. The temperature and the mass fraction of water in the gas phase are depicted in (a) and (b), respectively. In (c), the velocity of air is shown in the free-flow domain while the water saturation is depicted in the porous medium. All snapshots were taken at the final simulation time.

In [14], four different coupling methods were compared on the basis of plots of along the interface. Figure 6 shows the data of the original study, plotted against the results obtained with the proposed method on a grid with the same discretization length. The coupling method “CM4” corresponds to the approach that involves solving local nonlinear problems, and is considered the most accurate approach in [14]. The plots show that the method presented in this work leads to very similar results without the need for an interface solver, while large deviations to the other coupling methods can be observed.

Figure 6: Comparison to [14]. The figures show and the temperature as obtained with the proposed method, plotted along the line - against the solutions presented in [14] with four different coupling methods. In (b) and (d), a close-up on the upper left part of the interface is shown.

Influence of interface grid conformity

So far, the grid in the porous medium was chosen such that its dual grid is conforming to the grid used in the free-flow domain. Unlike all other test cases that we have investigated so far, we have observed an oscillatory behavior for this test case when deviating from a box-conforming grid configuration. These oscillations also occur when neglecting non-isothermal and compositional effects. Therefore, to avoid unnecessary complexity, we will consider a single-phase system throughout the remainder of this section. The corresponding results are shown in Fig. 8, where is plotted over the line located on the upper part of the free-flow porous-medium interface. When using a box-conforming grid (see Fig. 7) such oscillations do not occur, as discussed before and shown in Fig. 6.

Figure 7: Conforming, box-conforming, and simplicial grid used for the test case where air flows around a porous obstacle. The interface factor is shown on the right. For the simplicial grids, faces that are located at the corners of the porous-medium domain are discretised with an interface factor of 0.5, for all other faces another factor can be set, e.g. 0.95 for the grid shown in this figure.

By neglecting non-isothermal and compositional effects and by considering only single-phase flow, the mass balance equation Eq. 4a for a stationary state simplifies to and the flux coupling condition Eq. 5a reduces to . Thus, the local discrete balance equation (see Section 3.2) associated with each box degree of freedom located at the interface is composed of interface fluxes, given by , and of interior fluxes given by Eq. 17. This means that derivatives of these local discrete balance equations (local residuals) at the interface with respect to free-flow velocity degrees of freedom are in the range of in contrast to derivatives with respect to box pressure Dofs which additionally scale with the permeability . Therefore, entries in the Jacobian matrix corresponding to local box residuals at the interface significantly vary. It is therefore not surprising that numerical artifacts may occur if this is additionally amplified by grid effects. These points are further discussed in the following and the oscillatory behavior is investigated in more detail.

Similar as for the convergence test case, three grid configurations, i.e. conforming, non-conforming, and box-conforming (see Fig. 7), are considered. For the non-conforming configuration, a triangular grid (simplices) is used in the porous-medium domain, which allows to specify the discretization lengths at , see Fig. 7.

Figure 8: The figures in the left column show plotted along the line - obtained on a conforming grid (upper row) and on a simplicial (non-conforming) grid with an interface factor of . The right column shows on the upper boundary of the porous medium. The discretization lengths for the free-flow grid is set to , (two times refined compared to Fig. 7).

To quantify the oscillations that occur at the interface (as exemplarily shown in Fig. 8), we introduce the total variaion


which accumulates the differences of vertical face velocities between neighboring free-flow faces along . In the following, we will use , that is, we only consider the upper part of the free-flow porous-medium interface. For simplicity, we omit the subscript and denote this measure as .

Figure 9 shows the total variation for different permeability values and Reynolds numbers for the grid configurations shown in Fig. 7. For comparison, the results when using the Tpfa method in the porous medium are also shown. Here, as before, CM1 denotes the simplified interface solver, where the cell pressures are used at the interface, whereas CM4 calculates an interface pressure from the coupling conditions by solving local problems, see [14] for more details. Since the Tpfa implementation currently available in DuMu is restricted to conforming meshes, only these results are shown. The integer values above the data points show the total number of detected sign changes, which can be calculated as , where denotes the sign function.

Figure 9: Total variation at upper part of the free-flow porous-medium interface, i.e. for different permeability values and Reynolds numbers. is chosen as the maximum absolute velocity value that occurs on this interface. The integer values above the data points are calculated as the total number of sign changes. For the box method, three grid configurations are considered, i.e. conforming, box-conforming, and non-conforming (simplices). The simplicial grid corresponds to an interface factor of , i.e. for all faces , see Fig. 7.
Figure 10: Total variation at upper part of the free-flow porous-medium interface, i.e. , plotted over various interface factors for different refinement levels when using the simplicial grids. is chosen as the maximum absolute velocity value at this interface. The integer values above the data points show the total number of detected sign changes. An interface factor of, e.g., 0.9 means that the length of a porous medium face at the interface if 0.9 times the length of a free-flow face.

Figure 9 shows that the largest total variation occurs for the conforming case with a maximum of 63 sign changes. Only for the box scheme on the box-conforming grid and for the Tpfa scheme using CM4 there are no sign changes. For the Tpfa CM1 scheme, two sign changes occur close to the upper left corner of the porous medium (). The results of the box scheme for simplices with interface factor is slightly worse. It can also be seen that the box scheme behaves better for larger permeability values. As mentioned before, derivatives with respect to box pressure Dofs directly scale with the permeability. Therefore, for larger permeability values, the magnitude of derivatives with respect to are closer to derivatives with respect to .

When changing the Reynolds numbers, as shown in the right plot of Fig. 9, the behavior of the different schemes does not significantly change. This is most likely again related to the fact that derivatives of the local box residuals do not scale with the Reynolds number, because boundary fluxes linearly depend on the free-flow velocities. Note that the area-weighted projection yielded very similar results, which are therefore not shown here.

To investigate the influence of the interface grid refinement ratio, the interface factor is varied. An interface factor of corresponds to the conforming case for all faces except for those that are located at the corners of the porous medium where is set (see Fig. 7). Figure 10 shows the results of the box scheme on the simplicial grids for different interface factors and for different grid refinement levels. Refinement 0 corresponds to discretization lengths of in the free-flow domain, which is halved for each refinement. The worst results are observed for (conforming case). For , which has also been used for the results of Fig. 9, the total variation and the number of sign changes already strongly decreases. For a factor of no sign changes are observed independent of the grid refinement level.

To summarize, by using a box-conforming grid or slightly refined grid () the oscillatory behavior can be prevented. With the box method used in the porous-medium subdomain, this can be easily achieved by using unstructured grids with local refinement towards the interface. Furthermore, it should be highlighted again that for all other test cases that we have considered so far, we did not observe any oscillatory behavior. In future work, we will investigate this in more detail and also compare our results with mortar approaches, where it is well-known that oscillations may occur if the discretization length of the mortar space is too fine compared to the coupled sub-domains.

4.3 Example: Cooling Filter

In order to demonstrate the flexibility and function of the coupling system discussed above, an additional case is developed and shown here. This case highlights the model’s flexibility for three dimensional geometries with non-conforming interfaces, the extension to non-isothermal models, and the simple incorporation of velocity dependent balance equations in the porous medium, here demonstrated with the Forchheimer term.

The setup, as seen in Fig. 11, incorporates two L-shaped free-flow channels, joined by a centrally located curved porous medium. Within the free-flow channels a structured grid of rectangular elements is used, where within the porous medium an unstructured tetrahedral grid is used.

Figure 11: To the left the full simulated domain is displayed. Highlighted to the right is the difference in grid elements between the free-flow channel and the porous-medium flow section, as well as the thermal boundary conditions set at the inlet and at the left surface of the porous-medium subdomain.

A pressure difference, , of is applied across the system using Dirichlet boundaries for pressure at both the inlet and outlet boundaries; flow will travel through the upper L-shaped channel, through the curved porous medium, and exit through the lower L-shaped channel. All other non-interface walls within the free flow are constrained with zero mass flux and no-slip conditions. Within the porous medium, all non-interface walls are also set to allow zero mass flux. The fluid parameters simulated here are based on DuMu ’s gaseous Air class defined within the H2OAir fluid system. The parameters porosity and permeability are set to and , respectively.

In addition, thermal boundary conditions are designated in the model. An increased temperature, , of is set at the free-flow domain inflow, shown in red in Fig. 11, as opposed to the ambient temperature, , . One external porous-medium boundary surface, shown in blue, will act as a heat sink. This flux is defined as a conductive flux, through a casing with a prescribed thickness, , and conductivity , and with a fixed external cooling temperature . In this case, the external cooling temperature is set to be , the insulation thickness to , and the conductivity . All other external boundaries are set to no-temperature flux.

In addition, the Forchheimer nonlinear extension to Darcy’s law [29], a velocity dependent term shown in Eq. 37, has been added to the momentum balance within the porous medium. The permeability, , is assumed to be a scalar value within this example. In this case, the coefficient is chosen as a fixed [29]. Using the box method, incorporation of this term is relatively simple and does not require any stencil extensions or velocity reconstructions, as would be the case for standard cell-centered finite volume methods. In this case, velocity values at the integration points are directly available from the shape functions.


As is seen in Figure 11, the interfaces between the free-flow and the porous-medium flow subdomains are non-conforming. This enables refinement to the interface and the use of unstructured grids in the porous-medium flow domain, which is required e.g. to model curved domains. Seen in red is an overlay of the free-flow domain’s cell boundaries at the interface. In black are the porous-medium domain’s element outlines. Across the two L-shaped free-flow domains the temperature decreases, and heat is further removed from the system across the cooling interface, as seen in Figure 12. Although this does not represent a concrete application case, these physics are representative of applications such as CPU cooling, where heat is removed from a porous system embedded within a larger flow system [41].

Figure 12: Flow fields are depicted in both domains using arrows, and the temperature field within the porous-medium domain is shown in the background. The cooling effect of the left cooling surface is seen to reach into the porous medium domain.

5 Conclusions

This work explores coupled two-domain free-flow and multi-phase porous medium flow models using a vertex-centered finite volume method to discretize the porous-medium flow domain. Use of this box discretization method in the porous-medium flow subdomain allows for the use of unstructured non-conforming grids, more flexible inclusion of interfacial coupling conditions, and a simpler incorporation of velocity dependent terms in the balance equations such as the Forchheimer extension of Darcy’s law, as shown in this work.

The convergence behavior of the proposed method is investigated in Section 4.1. Therein, conforming, non-conforming and box-conforming grid configurations are considered and convergence orders are determined across levels of refinement in comparison to an analytical solution for two projection operators and . For both projection operators and all considered grids, second order convergence is observed.

The performance of the proposed method in a non-isothermal, multi-component and multi-phase environment is examined in 4.2. The accuracy of the coupling system was then further investigated to provide insight into limitations due to oscillations when dealing with non box-conforming unstructured grids. Recommendations for grid constraints are made to ensure non-oscillatory behavior. In practice, it suffices to locally refine the grid used in the porous medium towards the interface.

Finally, an exemplary 3D application is presented, which illustrates the capabilities of the proposed method: The porous-medium subdomain is curved and is discretized with an unstructured tetrahedral grid, while the interfaces to the neighboring free-flow channels are non-conforming. Moreover, within the porous-medium subdomain, the nonlinear velocity dependent Forchheimer extension of Darcy’s law is incorporated, and heat transport is evaluated across the full coupled system.

As the use and development of coupled porous medium free-flow models are a field of active research, future work in this area will include both extensions to the discretization methods as well as developments to the incorporated mathematical models. At the domain interface, further coupling conditions should be investigated to better resolve the near interface conditions [13]

. Additionally, at the moment the two domain system is solved monolithically. Comparing the existing monolithically solved method with a partitioned model where domains are solved individually

[24], or to other domain decomposition models using interface mortars [5], could also provide insight to other solution methods. For expansions to the free-flow domain, turbulence models and radiation based expansions to the energy balance could be introduced to better simulate atmospheric free-flow conditions [20, 9]. To expand on the discretization in the free-flow domain, unstructured and adaptive staggered finite volume methods are being developed in order to efficiently resolve more complex geometries and flow cases.


We thank the German Research Foundation (DFG) for supporting this work by funding SFB 1313, Project Number 327154368, Research Project A02, and by funding SimTech via Germany’s Excellence Strategy (EXC 2075 – 390740016).

Appendix A

The following tables list the errors obtained in the convergence test of Section 4.1 on the different grids and different projections used.

free flow darcy flow
0 1.44e-01 - 2.54e-03 - 8.97e-03 - 2.52e-02 -
1 3.76e-02 1.94e+00 5.35e-04 2.25e+00 2.10e-03 2.10e+00 6.11e-03 2.05e+00
2 9.52e-03 1.98e+00 1.27e-04 2.07e+00 5.15e-04 2.03e+00 1.50e-03 2.03e+00
3 2.39e-03 1.99e+00 3.15e-05 2.02e+00 1.28e-04 2.01e+00 3.72e-04 2.01e+00
4 5.98e-04 2.00e+00 7.84e-06 2.00e+00 3.20e-05 2.00e+00 9.26e-05 2.00e+00
5 1.49e-04 2.00e+00 1.96e-06 2.00e+00 7.99e-06 2.00e+00 2.31e-05 2.00e+00
Table 1: Convergence study. Errors and rates for the primary variables (conforming case, l2-projection).
free flow darcy flow
0 1.44e-01 - 2.54e-03 - 8.97e-03 - 2.52e-02 -
1 3.76e-02 1.94e+00 5.35e-04 2.25e+00 2.10e-03 2.10e+00 6.11e-03 2.05e+00
2 9.52e-03 1.98e+00 1.27e-04 2.07e+00 5.15e-04 2.03e+00 1.50e-03 2.03e+00
3 2.39e-03 1.99e+00 3.15e-05 2.02e+00 1.28e-04 2.01e+00 3.72e-04 2.01e+00
4 5.98e-04 2.00e+00 7.84e-06 2.00e+00 3.20e-05 2.00e+00 9.26e-05 2.00e+00
5 1.49e-04 2.00e+00 1.96e-06 2.00e+00 7.99e-06 2.00e+00 2.31e-05 2.00e+00
Table 2: Convergence study. Errors and rates for the primary variables (conforming case, weighted dof-projection).
free flow darcy flow
0 1.64e-01 - 2.29e-03 - 8.69e-03 - 2.61e-02 -
1 4.64e-02 1.82e+00 4.67e-04 2.29e+00 2.02e-03 2.10e+00 6.58e-03 1.99e+00
2 1.21e-02 1.94e+00 1.08e-04 2.11e+00 4.89e-04 2.05e+00 1.66e-03 1.99e+00
3 3.11e-03 1.96e+00 2.69e-05 2.01e+00 1.21e-04 2.01e+00 4.17e-04 1.99e+00
4 7.81e-04 2.00e+00 6.67e-06 2.01e+00 3.02e-05 2.01e+00 1.04e-04 2.00e+00
5 1.96e-04 1.99e+00 1.68e-06 1.99e+00 7.56e-06 2.00e+00 2.61e-05 2.00e+00
Table 3: Convergence study. Errors and rates for the primary variables (simplices, l2-projection).
free flow darcy flow
0 1.53e-01 - 2.19e-03 - 8.26e-03 - 2.61e-02 -
1 4.35e-02 1.81e+00 4.60e-04 2.25e+00 2.01e-03 2.04e+00 6.58e-03 1.99e+00
2 1.12e-02 1.96e+00 1.06e-04 2.12e+00 4.85e-04 2.05e+00 1.66e-03 1.99e+00
3 2.92e-03 1.94e+00 2.62e-05 2.01e+00 1.20e-04 2.01e+00 4.17e-04 1.99e+00
4 7.29e-04 2.00e+00 6.35e-06 2.04e+00 2.97e-05 2.02e+00 1.04e-04 2.00e+00
5 1.84e-04 1.98e+00 1.62e-06 1.97e+00 7.49e-06 1.99e+00 2.62e-05 2.00e+00
Table 4: Convergence study. Errors and rates for the primary variables (simplices, weighted dof-projection).
free flow darcy flow