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 multidomain systems that are challenging to resolve with an isolated flow model are coupled freeflow and multiphase porousmedium 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 soilwater 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 freeflow and porousmedium flow systems at this REV scale can then be done using a single or twodomain approach. A singledomain 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 twodomain approach, momentum conservation in each domain can be described using a model designed for its context, in this case, the NavierStokes equations for describing free flow, and the multiphase Darcy’s law for describing multiphase 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 twodomain models have been developed. Some models, e.g. [14, 33], couple cellcentered finite volumes, used to discretize the porous medium, with a MarkerandCell (staggered grid) finite volume scheme [18], used for the discretization of the freeflow domain. Other approaches utilize the same discretization method in each domain, using either staggered grids [22], finite elements [11], or a vertexcentered finite volume method, otherwise known as the box scheme [28, 21] for both domains. Other models have been developed to use an extra lowerdimensional 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]. Vertexcentered 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 NavierStokes equations discretized with the staggeredgrid method. The porousmedium flow is modeled using the multiphase Darcy’s law including the Forchheimer term, and is discretized using a vertexcentered 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 freeflow and the porousmedium 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 freeflow and porousmedium 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 gasphase 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 freeflow 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
(3a)  
(3b)  
(3c) 
together with appropriate boundary conditions. The unknown variables are the gasphase velocity , pressure and temperature , and one of the mass or mole fractions . Here, and denote the potentially solutiondependent gasphase 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 PorousMedium Flow
The equations governing compositional nonisothermal twophase (water and air) flow in the porousmedium subdomain, , are given by: equationparentequation
(4a)  
(4b)  
(4c) 
together with appropriate boundary conditions. is the saturation of phase while is the porosity of the permeable medium. and are the solidphase 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 freeflow 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 freeflow domain. We therefore impose the following interface conditions: equationparentequation
(5a)  
(5b)  
(5c)  
(5d)  
(5e)  
(5f) 
The momentum transfer normal to the interface is given by (5e) [26]. Condition (5f) is the commonly used BeaversJosephSaffman 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 nonparallel 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

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

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 . 
is the set of control volumes (dual grid elements) such that . For each control volume , denotes its volume.

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 .

is the set of cell centers such that and is starshaped with respect to .

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 threedimensional 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 Staggeredgrid scheme
A staggeredgrid finite volume scheme, also known as MAC scheme [18, 33], is used in the freeflow 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 oscillationfree 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 freeflow 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.
3.2 Vertexcentered finite volume scheme
For typical cellcentered 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 threedimensional 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 freeflow domain, see Fig. 1.
Each control volume is partitioned into subcontrol volumes , and with we denote the set of subcontrol volumes embedded in the control volume such that . Each subcontrol 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:
(6) 
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
(7) 
Let us now split the integral over the control volume into integrals over the subcontrol volumes, and let us introduce the discrete flux over a face :
(8) 
Furthermore, let us denote with and the conserved quantity and the source term evaluated for the subcontrol volume . Assuming that the time derivative is sufficiently smooth, this yields:
(9) 
Finally, we approximate the time derivative by a finite difference:
(10) 
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 subcontrol 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:
(11) 
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:
(12) 
where is the position of the vertex and denotes the KroneckerDelta. As seen in the above condition, the basis functions have a small support region, i.e. they are only nonzero within the elements adjacent to :
(13) 
The basis functions further describe a partition of unity, and we can express the discrete approximation of at a position with:
(14) 
Extending from (14), the expression for the discrete gradient of is:
(15) 
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:
(16) 
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
(17a)  
(17b)  
(17c) 
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 subcontrol volumes.
Similarly, we define the discrete heat fluxes as: equationparentequation
(18a)  
(18b) 
Here, the effective thermal conductivity on the face is again defined as the harmonic mean of the conductivities in the two subcontrol 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 subcontrol 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 :
(19) 
and correspondingly for the energy balance Eq. 4b:
(20) 
where .
3.3 Coupling
This section is devoted to the incorporation of the coupling conditions 5 in the discrete setting. The discretizations of the subdomains are allowed to be nonmatching at the interface, which results in arbitrary overlap between grid faces of the freeflow and the porousmedium domain, respectively. Let us partition the interface into a set of facets (line segments for the twodimensional case), for which we note that . Each facet is the intersection between a face of the discretization of the porousmedium domain and a face of the discretization of the freeflow 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 .
The coupling between the two subdomains is realized as follows: Eqs. 5d and 5b represent Neumann (flux) boundary conditions, where only the freeflow 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 vertexcentered finite volume scheme, where degrees of freedom are naturally located at the interface. In the following, we present the construction of the freeflow fluxes at the interface by using a twopoint 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 freeflow 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 freeflow grid (see Fig. 1) such that . Integrating the freeflow related terms of Eqs. 5d and 5b over yields the following discrete approximations:
(21)  
(22) 
with flux approximations
(23)  
(24) 
which are related to the control volume and the corresponding subcontrol volume of the porousmedium domain such that . The interface density is evaluated with the following average, but the diffusion coefficient is evaluated as , using only information from the freeflow 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 freeflow or the the porousmedium 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 porousmedium solution. This is done with two types of projection operators :
(25)  
(26) 
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 freeflow subdomain. Equation 5e yields a Neumann boundary condition for the momentum balance by evaluating at each freeflow face. This is again realized by using the projections introduced above. When using the Saffman simplification of the BeaversJoseph condition, as shown in Equation 5f, the porousmedium velocity is neglected meaning only freeflow unknowns are required and the condition is independent of the chosen porousmedium discretization scheme.
4 Numerical results
The coupling approach presented in the previous section was implemented within the the opensource simulation environment DuMu
[25], which exists as an additional module of the DUNE project [4]. We employ a monolithic approach, where both subproblems 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 DuneSubgrid [15] and DuneUGGrid modules are used to represent the structured grids in the freeflow and the unstructured grids used in the porousmedium 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 singlephase flow, neglecting nonisothermal, compositional, and gravitational effects. Consider a domain , decomposed into and with interface . The material parameters are given by
(27) 
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
(28)  
(29) 
Neglecting time derivatives, we have in the freeflow domain
(30)  
(31) 
and in the porous medium
(32)  
(33) 
The validity of the coupling conditions can also be easily verified. In order to guarantee that is positive definite and a full tensor with nonnegligible offdiagonal 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 5thorder quadrature rule.
We use the following discrete pressure error norm,
(34) 
where and where the parameter indicates the th refinement level. Furthermore, and denote the control volumewise constant discrete and exact pressure values, obtained by evaluating both in the center of the control volumes .
For the staggeredgrid scheme, we additionally quantify the errors for the velocity unknowns. These errors are calculated as
(35) 
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 porousmedium domain are considered: a grid that is conforming to the grid used in the freeflow domain (fig:convtest_grids_conforming), an unstructured simplex grid with nonmatching 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 freeflow domain (fig:convtest_grids_boxconforming). These cases are refered to here as the conforming, nonconforming, and boxconforming 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 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 boxconforming grid configuration, and the largest errors occur on the nonconforming 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 boxconforming grid, the projections and are equivalent, which is why we only show the result for here.
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 cellcentered Tpfa scheme on a structured grid has been used in the porousmedium domain coupled to a staggeredgrid finite volume scheme used in the freeflow domain. When solving physically complex flow processes, described by nonlinear PDEs, the additional interface unknowns have to be calculated locally within a subroutine 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 twodimensional freeflow 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 noflow. For indications on the domain dimensions, see Fig. 4.
A twophase, twocomponent model is used in the porous medium (see Eq. 4), while a singlephase, twocomponent model is considered in the freeflow 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 freeflow 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 freeflow 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.
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.
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 freeflow 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 boxconforming grid configuration. These oscillations also occur when neglecting nonisothermal and compositional effects. Therefore, to avoid unnecessary complexity, we will consider a singlephase 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 freeflow porousmedium interface. When using a boxconforming grid (see Fig. 7) such oscillations do not occur, as discussed before and shown in Fig. 6.
By neglecting nonisothermal and compositional effects and by considering only singlephase 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 freeflow 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, nonconforming, and boxconforming (see Fig. 7), are considered. For the nonconforming configuration, a triangular grid (simplices) is used in the porousmedium domain, which allows to specify the discretization lengths at , see Fig. 7.
To quantify the oscillations that occur at the interface (as exemplarily shown in Fig. 8), we introduce the total variaion
(36) 
which accumulates the differences of vertical face velocities between neighboring freeflow faces along . In the following, we will use , that is, we only consider the upper part of the freeflow porousmedium 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 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 boxconforming 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 freeflow velocities. Note that the areaweighted 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 freeflow 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 boxconforming grid or slightly refined grid () the oscillatory behavior can be prevented. With the box method used in the porousmedium 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 wellknown that oscillations may occur if the discretization length of the mortar space is too fine compared to the coupled subdomains.
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 nonconforming interfaces, the extension to nonisothermal 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 Lshaped freeflow channels, joined by a centrally located curved porous medium. Within the freeflow channels a structured grid of rectangular elements is used, where within the porous medium an unstructured tetrahedral grid is used.
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 Lshaped channel, through the curved porous medium, and exit through the lower Lshaped channel. All other noninterface walls within the free flow are constrained with zero mass flux and noslip conditions. Within the porous medium, all noninterface 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 freeflow domain inflow, shown in red in Fig. 11, as opposed to the ambient temperature, , . One external porousmedium 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 notemperature 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 cellcentered finite volume methods. In this case, velocity values at the integration points are directly available from the shape functions.
(37) 
As is seen in Figure 11, the interfaces between the freeflow and the porousmedium flow subdomains are nonconforming. This enables refinement to the interface and the use of unstructured grids in the porousmedium flow domain, which is required e.g. to model curved domains. Seen in red is an overlay of the freeflow domain’s cell boundaries at the interface. In black are the porousmedium domain’s element outlines. Across the two Lshaped freeflow 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].
5 Conclusions
This work explores coupled twodomain freeflow and multiphase porous medium flow models using a vertexcentered finite volume method to discretize the porousmedium flow domain. Use of this box discretization method in the porousmedium flow subdomain allows for the use of unstructured nonconforming 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, nonconforming and boxconforming 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 nonisothermal, multicomponent and multiphase 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 boxconforming unstructured grids. Recommendations for grid constraints are made to ensure nonoscillatory 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 porousmedium subdomain is curved and is discretized with an unstructured tetrahedral grid, while the interfaces to the neighboring freeflow channels are nonconforming. Moreover, within the porousmedium 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 freeflow 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 freeflow domain, turbulence models and radiation based expansions to the energy balance could be introduced to better simulate atmospheric freeflow conditions [20, 9]. To expand on the discretization in the freeflow domain, unstructured and adaptive staggered finite volume methods are being developed in order to efficiently resolve more complex geometries and flow cases.Acknowledgements
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.44e01    2.54e03    8.97e03    2.52e02   
1  3.76e02  1.94e+00  5.35e04  2.25e+00  2.10e03  2.10e+00  6.11e03  2.05e+00 
2  9.52e03  1.98e+00  1.27e04  2.07e+00  5.15e04  2.03e+00  1.50e03  2.03e+00 
3  2.39e03  1.99e+00  3.15e05  2.02e+00  1.28e04  2.01e+00  3.72e04  2.01e+00 
4  5.98e04  2.00e+00  7.84e06  2.00e+00  3.20e05  2.00e+00  9.26e05  2.00e+00 
5  1.49e04  2.00e+00  1.96e06  2.00e+00  7.99e06  2.00e+00  2.31e05  2.00e+00 
free flow  darcy flow  

0  1.44e01    2.54e03    8.97e03    2.52e02   
1  3.76e02  1.94e+00  5.35e04  2.25e+00  2.10e03  2.10e+00  6.11e03  2.05e+00 
2  9.52e03  1.98e+00  1.27e04  2.07e+00  5.15e04  2.03e+00  1.50e03  2.03e+00 
3  2.39e03  1.99e+00  3.15e05  2.02e+00  1.28e04  2.01e+00  3.72e04  2.01e+00 
4  5.98e04  2.00e+00  7.84e06  2.00e+00  3.20e05  2.00e+00  9.26e05  2.00e+00 
5  1.49e04  2.00e+00  1.96e06  2.00e+00  7.99e06  2.00e+00  2.31e05  2.00e+00 
free flow  darcy flow  

0  1.64e01    2.29e03    8.69e03    2.61e02   
1  4.64e02  1.82e+00  4.67e04  2.29e+00  2.02e03  2.10e+00  6.58e03  1.99e+00 
2  1.21e02  1.94e+00  1.08e04  2.11e+00  4.89e04  2.05e+00  1.66e03  1.99e+00 
3  3.11e03  1.96e+00  2.69e05  2.01e+00  1.21e04  2.01e+00  4.17e04  1.99e+00 
4  7.81e04  2.00e+00  6.67e06  2.01e+00  3.02e05  2.01e+00  1.04e04  2.00e+00 
5  1.96e04  1.99e+00  1.68e06  1.99e+00  7.56e06  2.00e+00  2.61e05  2.00e+00 
free flow  darcy flow  

0  1.53e01    2.19e03    8.26e03    2.61e02   
1  4.35e02  1.81e+00  4.60e04  2.25e+00  2.01e03  2.04e+00  6.58e03  1.99e+00 
2  1.12e02  1.96e+00  1.06e04  2.12e+00  4.85e04  2.05e+00  1.66e03  1.99e+00 
3  2.92e03  1.94e+00  2.62e05  2.01e+00  1.20e04  2.01e+00  4.17e04  1.99e+00 
4  7.29e04  2.00e+00  6.35e06  2.04e+00  2.97e05  2.02e+00  1.04e04  2.00e+00 
5  1.84e04  1.98e+00  1.62e06  1.97e+00  7.49e06  1.99e+00  2.62e05  2.00e+00 
free flow  darcy flow  

Comments
There are no comments yet.