Fluid flow through porous media is typically modeled with the Darcy equations Darcy1856; Whitaker1986
. When there are large cracks and voids in the porous medium, then a homogenization of the material into a single permeability tensor is no longer appropriate. The creeping flow in those domains can more accurately be modeled with the Stokes equations. The coupling relations between the Darcy and Stokes flow regimes were first studied by Beavers and JosephBeavers:67.1, and later supplemented by Saffman Saffman:71.1. The resulting Stokes/Darcy equations coupled via Beavers-Joseph-Saffman interface conditions play an important role in many disciplines, for instance modeling groundwater flow Discacciati:02.1, petroleum and karst reservoirs Popov2009; Wu2006, in-situ leach mining Forbes1994), perfusion of blood through tissue Stoter2017; Mosharaf2019, filtration devices Nassehi1998; Hanspal2006, and chemical reactors Schneider2005. The wide range of applications has lead to significant interest in methods that computationally approximate the coupled Stokes/Darcy problem, including their numerical analysis and algorithmic treatment Layton:02.1; Discacciati:07.1; Chen:11.1; Baber:12.1; Discacciati:09.1; Burman:07.1; Badia:09.1; Chidyagwai:09.1; Cao2010. Many of these applications are characterized by incomplete or uncertain data in terms of the geometry and topology of the two flow domains (e.g., uncertain subsurface soil characteristics, limited resolution of CT or MRI scans). At the same time, they generally involve a wide range of length scales (e.g., from small rock cavities to a complete mining site, from small capillaries to a full organ).
In this article, our objective is to enable the efficient simulation of such multiscale systems with parametrically defined free-flow Stokes domains at the lowest scale. To achieve this objective, we make use of model order reduction by means of reduced basis methods Chinesta2017; quarteroni2015; Hesthaven2016
on small repetitive subdomains. A reduced basis method replaces a computationally expensive high-dimensional finite element discretization of a parametrized partial differential equation (PDE) by a small low-dimensional set of basis functions that have high approximation power with respect to the solution manifold of the parametrized problemWillcox2002; Binev2011. We develop a reduced basis method for the coupled Stokes/Darcy model on a variable internal geometry, and model the large-scale system consisting of many repetitive subdomains by coupling together many such reduced basis functions. Similar approaches to coupling reduced basis functions on repetitive subdomains have been used in Antil2010; Antil2011; Baiges2013; Rademacher2014; Ferrero2018; Xiao2019.
In a reduced basis context, the reduced basis functions themselves are defined on a high-dimensional finite element approximation space. Consequently, the high-resolution finite element mesh must remain fixed. This requirement conflicts with our aim of flexibly handling the geometry of the internal Stokes and Darcy domains. Usually, parametric geometry is handled by mapping back to a reference domain quarteroni2015; Hesthaven2016, but such an approach does not permit changes of topology as this would degenerate the Jacobian of the mapping. We therefore require a method that is able to model a topologically flexible Stokes domain that may merge and disperse within the Darcy domain. Diffuse interface methods Maury:01.1; Raetz:06.1; Ramiere:07.1; Li:09.1; Lervaag:14.1, also known as diffuse domain or phase-field methods, offer such a flexible framework for solving coupled boundary value problems on non-boundary fitted meshes. The internal geometry is implicitly represented by a phase-field function, which smoothly transitions from zero to one. The phase-field indicator function and its gradient can be leveraged to replace integrals on subdomains or interfaces by weighted volumetric integrals on the complete domain. The resulting phase-field formulation is equivalent to the sharp-boundary interface problem when the width of the diffuse interface (controlled by a characteristic length-scale parameter) limits to zero. Phase-field geometry representations have been widely applied, for instance in growth modeling of tissues and crystals Oden2010; Golovin1998, for tracking the evolution of crack patterns Miehe:10.2; Borden:12.1; Schillinger:15.2; Mikelic:15.1, for enforcing boundary conditions in imaging data based analysis Nguyen:17.1; Nguyen:18.1, for modeling multi-phase flow Boyer:99.1; Aland:10.1; Teigen:11.1
, and for modeling phase transition and segregation processesAnders:12.1; Teigen:09.1; Liu:15.1; Zhao:15.1.
A crucial aspect of a reduced basis method is the affine decomposition of the linear and bilinear forms. Affine decomposition enables the precomputation of reduced-basis stiffness matrices such that the final reduced order model is completely independent of the size of the original high-fidelity model. However, the (bi)linear forms resulting from the diffuse representation of the internal geometry do not satisfy this criterion. This is a common issue for many relevant parametrized PDEs, such as those that feature nonlinearities Chaturantabut2010; Chaturantabut2012; Drohmann2012; Amsallem2012, complex material laws Ryckelynck2009; Kerfriden2011; Fauque2018, or in general, spatially varying model coefficients Barrault2004; Grepl2007. For those cases, the discrete empirical interpolation method (DEIM) may be used Barrault2004; Maday2015. With DEIM, all non-affine parameter dependent fields are replaced by a low-dimensional interpolation on DEIM modes. For our diffuse representation, this method solves the problem of non-affine parameter dependence, but the interpolation of a domain indicator function is likely to produce Gibb’s-type oscillations at regions with high gradients (i.e., the interfaces). Oscillations necessarily imply negative values and values larger than one, which in turn could produce nonphysical domain representations and unstable system matrices. In this paper, we mitigate this issue by introducing a non-negativity preserving version of DEIM.
Our article is structured as follows: in Section 2, we derive the diffusely coupled Stokes/Darcy equations. We show that all three Beavers-Joseph-Saffman conditions can be treated naturally within a diffuse interface framework. In Section 3, we propose a non-negativity preserving variation of the discrete empirical interpolation method for the dimensional reduction of the phase-field geometry representation and show its effectiveness for three benchmark problems. In Section 4, we use the same three benchmark problems to study the relation between the required number of phase-field DEIM modes and the number of reduced basis functions in the reduced order model. In Section 5, we apply our methodology to efficiently characterize the subsurface flow characteristics of an in-situ leach (ISL) mining site that consists of a large number of repetitive hexagonal units. In Section 6, we summarize our work and draw some conclusions.
2 Diffusely coupled Stokes/Darcy equations on parametrically defined domains
We consider an incompressible fluid moving through a partially porous medium at velocities that are sufficiently small to neglect the convective components in the material derivative. Steady state equilibrium of the fluid is then described by:
with the Cauchy stress tensor, the fluid velocity, the body force, and the -dimensional spatial domain.
2.1 Weak formulation of the coupled Stokes/Darcy equations
Our porous medium in contains voids and cracks where the flow is unobstructed. We denote the union of these (potentially disconnected) subdomains , and in these subdomains we use the constitutive relation of an incompressible Newtonian fluid:
with the viscosity and the pressure field. Hence, in , the governing equations are the Stokes equations.
In the remaining domain, , the interaction between the porous medium and the fluid produces a flow resistance that is assumed to be sufficiently high such that the viscous effects in the Newtonian fluid may be neglected. The closing relations become:
where is the second order permeability tensor. These assumptions lead to the Darcy equations in .
The interface that couples the Stokes and Darcy domains is denoted . On
, the unit normal vectorsand point out of the Stokes and Darcy domains, respectively. Across the interface, the solution fields are coupled due to the required balance of mass and balance of linear momentum. Additionally, the porous material introduces a shearing resistance to the fluid on the Stokes domain. These considerations lead to the so-called Beavers-Joseph-Saffman coupling conditions:
where subscripts and refer to the solution fields in the Stokes and Darcy domains, respectively. The tensor is the tangential projector and is the shear resistance parameter.
On the boundary of , denoted , we permit the following boundary conditions:
The first two conditions represent no-inflow with free-slip and the last condition is a pressure condition. These specific boundary conditions are chosen because they are valid essential and natural conditions for both the Stokes and the Darcy equations. The inflow condition is chosen homogeneous for simplicity, but it may just as well be set to a non-zero value.
We obtain a weak formulation by multiplying Eq. 1 by test functions and , integrating over the domain , substituting the constitutive relations Eq. 2 and Eq. 3 on their respective domains, performing integration by parts, and by substituting the coupling and boundary conditions of Eqs. 5 and 4. This leads to the statement:
where the function space is defined as:
This space ensures sufficient regularity on the vector valued functions in both the Stokes and Darcy domains.
2.2 Diffuse interface representation
Next, we consider a parametrically defined domain (and thus also ), where is a point in the parameter space . The parameter space is finite dimensional and encodes in some sense the set of potential geometries of . For example, parameters in could denote the number, position, size and/or orientation of voids and cracks through the domain . To handle this extensive geometric flexibility, we introduce a diffuse representation of the interface geometry in the weak statement. The strength of such an implicit interface representation is that the computational mesh does not have to fit the internal interface.
Let be a “phase-field” indicator function with values between 1 and 0. The function tends to 1 in the Stokes domain and to 0 in the Darcy domain, and it monotonically decreases from 1 to 0 along straight lines that cross a thin region around the Stokes/Darcy interface. The thickness of the diffuse interface is characterized by the parameter , as illustrated in Fig. 1. Based on this indicator function, volume integrals on or and surface integrals on can be approximated by volume integrals on according to:
In the following, we indicate dependencies on only when relevant.
The normal vector may be approximated as , such that:
Finally, we assume that the solutions are such that , whereby on can be approximated as twice the solution in the diffuse interface.
By making use of these approximations in Eq. 6, the new weak statement becomes:
where we introduce the new second order tensor field , defined as:
and the function space , defined as:
2.3 Relation with the Brinkman model
Diffuse geometry representations are often plagued by non-physical behavior in the diffuse interface region. The only way to address this issue is to reduce the length-scale parameter and thereby reduce the interface width. As a result, significant adaptive mesh refinement is required to ensure that the induced error is below a tolerance threshold, which can easily lead to a prohibitive increase in computational cost. We would like to point out, however, that the weak formulation (10) corresponds to model equations with physical relevance even for a non-vanishing interface width.
Assuming that the obtained solution pair is sufficiently smooth, we can perform integration by parts to arrive at the following statements:
which hold for all and all . The corresponding strong form equations are:
with and . These are the Brinkman equations Brinkman1949, which are used frequently for describing a sufficiently fast moving fluid in porous media Popov2009; Juntunen2010; Mosharaf2019. The field , which was originally a domain indicator, is now a material parameter, and the third Beavers-Joseph-Saffman coupling condition (4) has emerged as an additional orthotropic contribution to the material porosity.
2.4 Finite element approximation
A finite element formulation based on the weak formulation (10) no longer requires a mesh that fits the interface. Instead, we are able to use the same mesh for computing the solution for all domain configurations in the parameter space . This is a crucial requirement for obtaining a reduced basis representation of the parametric problem later on.
The mixed nature of the equations, however, warrants careful selection of the finite element spaces. This is particularly challenging for the coupled Stokes/Darcy equations: the well-posedness requirements, i.e., the Ladyzhenskaya-Babuška-Brezzi (LBB) conditions, in the Stokes limit are different from those in the Darcy limit Boffi2013. Since our domain is parametric, but our mesh and approximation spaces are fixed, we seek a pressure/velocity pair of elements that is stable for both a pure Stokes problem and a pure Darcy problem. So, we require approximation spaces that are subspaces of the relevant function spaces for all domain configurations:
and satisfy the LBB conditions in the limiting cases of both the Stokes and the Darcy equations. One such pair is the combination of linear nodal elements for the construction of and so-called MINI elements for the construction of Mardal2002; Juntunen2010. The MINI element is a simplicial element and consist of linear interpolation functions enriched with a bubble function Arnold1984. We use this combination of elements throughout the remainder of this article. For all simulations, we make use of the FEniCS finite element library AlnaesBlechta2015a for computing the large stiffness matrices and the PETSc library for manipulating those matrices petsc-web-page.
3 Non-negativity preserving discrete empirical interpolation of the phase-field
To ensure efficient formation of the stiffness matrix in the reduced order model, it is imperative that the bilinear and linear forms under consideration depend on the parameter in an affine sense. That is, we require:
where and are the bilinear and linear forms corresponding to Eq. 10. The (bi)linear forms and are parameter independent, and the parameter dependency of and is captured through multiplication with the scalar-valued functions and .
3.1 The discrete empirical interpolation method
In our case, the complex parameter dependency of the phase-field does not naturally permit a decomposition of the bilinear form as shown in Eq. 16a. To remedy this deficiency, the affine decomposition can be approximated with the discrete empirical interpolation method (DEIM). DEIM produces the following low-dimensional representation of a field:
where are the parameter independent functions in space (the DEIM interpolation modes) and are the corresponding parameter dependent weighting functions. In summary, the DEIM algorithm for obtaining such a decomposition is a two-step procedure:
The manifold is explored and the modes that best represent the complete manifold in some norm are identified. In practice, a discrete representation of the manifold is obtained as a snapshot matrix of solution vectors of the projection of the functions onto a finite element space for a sampling of the parameter space. The modes
follow from a principal orthogonal decomposition or singular value decomposition of the snapshot matrix.
The functions are designed such that Eq. 17 induces an interpolation of by the modes at a selection of points. These points are iteratively determined: the -th point is the location where the interpolation of the -th mode by the first modes produces the largest error. When the interpolation locations are defined for , the weighting values are computed by sampling and by solving the interpolation system of equations:
The reader is referred to Barrault2004; Grepl2007 for more details.
3.2 Positivity preserving DEIM
The following properties of the phase-field and tensor field guarantee that the stiffness matrix is positive semi-definite:
These conditions are naturally satisfied by the true phase-field. To guarantee stability of the reduced basis method, they must also be satisfied by the DEIM approximation of these fields after the reconstruction strategy of Eq. 17. This is challenging due to the nature of the possible fields : they are mostly constants of 1 or 0 except at thin regions where they have steep gradients. The non-local support of the DEIM interpolation functions would produce Gibbs-like oscillations in the neighborhood of the interface.
To remedy this issue, we propose to rewrite the fields (19) in the following form:
We will now approximate the new fields , and with Eq. 17 and the DEIM algorithm. The affine representation of the original fields then follows by expanding the summation multiplication:
where, in Eq. 21c, is taken for simplicity of notation.
Since the final approximations are constructed from squares, these series expansions will satisfy the requirements (19). The downside of this approach is the increase of the cost for the summations involved.
3.3 Benchmark problems
In the following, we will investigate the relation between the quality of the approximation and the required number of modes for three benchmark problems. For each of these problems, we illustrate some of the DEIM modes and , and we show reconstructions of the true fields and for representative parameter points. For the first benchmark, we also compare our positivity preserving DEIM method to a standard DEIM method.
3.3.1 Circular domain with an inwards spiraling channel
The first benchmark problem consists of a circular porous domain with a circular hole in the center. The phase-field, representing Stokes flow, is spiraling inwards. The parameter space is one-dimensional and denotes the number of full rotations of the spiral. The problem set-up is illustrated in Fig. (a)a, and Figure (b)b shows the high-fidelity mesh that is used for all subsequent computations.
First, we obtain the DEIM approximation modes of the fields , and . We use a homogeneously spaced discrete approximation of the parameter space for the parameter sampling: . For each we construct discrete approximations , and
as interpolants of the true fields on the high-fidelity mesh. We collect the solution vectors corresponding to the interpolated fields in a matrix, and performing a singular value decomposition on this snapshot matrix. The first four left-singular-vectors of the fieldsare shown in Fig. 3.
As an example, we use these interpolation modes to reconstruct the field at the parameter point . The original field is shown in Figure (a)a. The other subfigures show the reconstructions with the non-negativity preserving DEIM method for different numbers of modes. We observe that the reconstruction of is non-negative everywhere, guaranteed via our reconstruction strategy of Eq. 21. To demonstrate the importance of non-negativity preserving reconstruction, we employ the standard DEIM method for the reconstruction of the same field with the same interpolation bases. The result for the minimal and maximal number of modes is shown in Fig. 5. As the color bar indicates, there are banded regions around the interface in the approximated phase-field that are negative. Such undershoots would result in a negative diffusive term, rendering the weak formulation (10) unstable.
The first four interpolation modes for the field are shown in Fig. 6, and an example reconstruction for the case with various numbers of modes is shown in Fig. 7. We observe that a higher number of interpolation modes is required to retrieve a reasonable approximation of the original field. This is to be expected; the localized support and the vector valued nature of make interpolation challenging. Nevertheless, with 15 modes a qualitatively reasonable approximation of the initial field can be obtained.
The singular values of the singular value decomposition of the snapshot matrix can be leveraged to quantify the approximation power of a certain number of modes of either one of the fields , and with respect to the entire (discrete) parameter set. The Eckart-Young-Mirskey theorem Eckart1936 suggests the following measure:
where are the singular values. The value for a given number of modes represents the optimal relative error (in the Frobenius norm) that can be achieved when representing the entire snapshot matrix with the first interpolation modes. As such, is a useful indicator for the required number of modes to reach a certain error tolerance.
Figure 8 shows for the different fields , and . The graphs confirm that out of the three fields, is the most challenging to represent accurately with a finite-dimensional linear approximation. Still, with a total of 50 modes, the measure for the relative error drops below 1%, which is more than sufficient for applications where interface geometries are described by a phase-field. The figure also shows a much quicker drop-off for the fields and . This indicates that a computationally optimal reduced basis method will require different numbers of interpolation modes to represent the different phase-field quantities, a concept that we will explore in more depth in Section 4.
3.3.2 Filtering device
We now consider a channel of Stokes flow is obstructed by three blocks of porous material. The location of these blocks is variable. Each block is permitted to displace by a distance of a single width compared to the base configuration as illustrated in Fig. 9. The three blocks can thus partially overlap to form two, or even just one single block. This benchmark illustrates the capability of describing internal geometry variations with changing topology. Using our phase-field representation of the geometry, we can naturally change the topological properties of the Stokes and Darcy domains. As emphasized in the introduction, this is not the case for the standard approach to incorporating geometric flexibility in a reduced basis method, which is based on a reinterpretation of the geometric variability as a mapping from a reference domain quarteroni2015; Hesthaven2016.
The parameter space of this benchmark problem is three-dimensional. We construct the sampling set as a tensor product of 21 uniformly spaced samples along each axis, yielding a set of 9261 samples: . When we use this sampling set for the DEIM algorithm, we obtain the DEIM modes for the fields and depicted in Figs. 12 and 10. Figures 13 and 11 then show reconstructions of the fields and of the parameter point . Note that this parameter point produces an overlap of the left two blocks of porous material.
The reconstructions in Figs. 12 and 10 make use of more interpolation modes compared to benchmark problem in order to obtain acceptable approximations. This is due to the higher dimensionality of the parameter space. This should be reflected in the distribution of the singular values. Figure 14 plots the error measure from Eq. 22. When compared to Fig. 8, we indeed observe a larger overall . Similar to before, it is the field of that exhibits the slowest decay of the error measure.
3.3.3 Rotating channel
An angled straight channel of Stokes flow is embedded in a rectangular domain of porous material. The parameter is the angle of rotation of the channel, such that the parameter space is one-dimensional. The problem is illustrated in Fig. (a)a, and the corresponding high-fidelity mesh is illustrated in Fig. (b)b. While this may seem like the simplest of the three benchmark problems, it is actually the one where the DEIM approximation is the least effective. This is due to the large Kolmogorov -width of the -manifold: for different parameter points , i.e., different angles of the rotating channel, the phase-field is largely uncorrelated. This makes the DEIM compression relatively ineffective. The purpose of this benchmark is to investigate the required number of DEIM modes for problems that are not particularly suitable for the suggested model order reduction technique.
We use as a discrete approximation of the full parameter space for the parameter sampling. Figures 18 and 16 shows the DEIM modes of the fields and , and examples of reconstructed fields and are shown in Figs. 19 and 17. The error measure is plotted in Fig. 20 for each of the three fields. These results confirm that this benchmark is particularly difficult to handle with the DEIM representation of the geometry. The reconstructed fields still exhibit quite significant errors away from the center of the channel, even for 20 modes. The error measure also indicates higher errors for the same number of modes compared to the previous benchmarks. Figure 20 indicates that modes are required to reach an error measure of for the -field, whereas and modes were required to reach the same threshold for benchmark problems 1 and 2 respectively (as shown in Figs. (d)d and (d)d).
4 Reduced basis method of the coupled equations on the parametric domain
From Section 2, we have the following parameter dependent finite element formulation:
where bilinear and linear forms are those of Eq. 10. The mixed finite element space is that described in Section 2.4 and could be kept independent of . This approximation space is assumed to be sufficiently refined such that the finite element solutions can for all intents and purposes be considered the true solution.
Now that we have affine representations of all parameter dependent fields, an affine decomposition of the bilinear form would follow as:
We now wish to find a low dimensional subspace (with ) that has high approximation power with respect to the full solution manifold . Such a reduced basis can again be obtained from a principal orthogonal decomposition or a singular value decomposition of a snapshot matrix generated from solution vectors of high-fidelity finite element solutions. We perform this computation for all three of the benchmark problems, and use the same parameter samplings as before for the generation of the snapshot matrices. The first four reduced basis functions are shown in Figs. 23, 22 and 21 for each of the three benchmarks. The figures show pressure fields overlapped by velocity fields. This is indicative of the nature of these functions: since we perform a singular value decomposition on the snapshot matrix of velocity-pressure pairs, each of the left singular vectors represents both a velocity and a pressure field.
The singular values corresponding to the reduced basis functions again relate to the approximation power of the cumulative basis compared to the full solution manifold. This is quantified by the measure from Eq. 22, which is plotted for all three benchmark problems in Fig. 24. We observe the same trends as we did for the DEIM geometry representation in Section 3: benchmark 1 performs best, followed by benchmark 2, and benchmark 3 is worst. Again, this can be motivated from the complexity of the problems, and the related Kolmogorov -width of the solution manifold. The higher dimensionality of benchmark 2 produces a higher Kolmogorov -width, but the irregularity of benchmark 3 has a higher impact.
where represents solution tuple that combines the velocity and the pressure solution. The matrix representation of this problem becomes:
where each and .
The reduced order model is operated by taking the following steps:
In the offline phase (i.e., before operation), the matrices corresponding to the bilinear forms , and are precomputed.
Also in the offline phase, the interpolation matrix from Eq. 18 corresponding to each of the weight vectors , and is inverted and stored.
In the online phase (i.e., during operation), the weight vectors , and are computed for the given parameter .
From those weights, we compute the weights , and for to the positivity preserving approximation of Eq. 21.
The summation from Section 4 is carried out and the resulting system of equation is be solved.
Figures 27, 26 and 25 show example results of the complete reduced order model for various numbers of reduced basis functions (i.e., various ). In all cases, 50 modes were used for the DEIM approximation of the different phase-field quantities to ensure that this did not visually affect the solution graphs.
The remaining open issue is the relation between the number of DEIM modes for each of the three fields, the number of solution modes in the reduced basis, and the effectiveness of the reduced order model throughout the parameter space. Recall that plotted in Figs. 24, 20, 14 and 8 represents the lowest possible relative error that may be achieved when approximating all solutions in the snapshot matrix with the first modes. The question remains open whether the DEIM approximation and the reduced order model yield solutions close to that optimum. At the same time, the measure concerns an average over the parameter space, which says little about the maximum relative error, which is arguably the more important measure. Nevertheless, the computed values may guide the design of the reduced order model in terms of numbers of modes for the various reduced order approximations.
As a first indicator of the quality of the reduced order model, and its dependence on the DEIM approximation of the phase-field geometry representation, we plot the maximum -error of the reduced order model over the complete (discrete) parameter space as a function of for various numbers of DEIM interpolation modes in Figs. (c)c, (b)b and (a)a for each of the three benchmark problems. The general trend in each of the figures is the same: for small values, the maximum -error is independent of the number of modes for the DEIM approximations. As increases above a certain threshold, the -error starts to drop. This drop plateaus at different values for the different numbers of DEIM modes. For all three benchmark problems, an acceptable maximal relative -error of is achieved for solution modes, where the relative error is computed with respect to the largest -norm in the discrete parameter space.
In Section 3, we show that the DEIM approximation of the three phase-field related quantities performs differently well. At the same time, the error in the reconstruction of the three fields may have different impact on the final solution error. These observations imply that different numbers of modes may be required for the DEIM approximation of the three fields. A more targeted and optimized choice for the different numbers of DEIM interpolation modes could reduce the computational expense of operating the reduced order model without adversely affecting the quality of the final solution. We propose to equate the -value for the particular choice of number of solution modes (i.e., those in Fig. 24) to a weighted sum of the -values for the DEIM interpolation modes (i.e., those in Figs. 20, 14 and 8
). As a weighting factor, we choose the maximal eigenvalueof the matrix corresponding to each of the fields (i.e., those of Eqs. 25d, 25c and 25b), for the central parameter point in the parameter space. That is, we use the minimal values , and for the DEIM reconstruction of the respective fields such that: