On bounded approximations of periodicity for computational homogenization of Stokes flow in porous media

By separation of scales and the homogenization of a flow through porous media, a two‐scale problem arises where a Darcy‐type flow is present on the macroscale and a Stokes flow on the subscale. In this paper, the problem is given as the minimization of a potential. Additional constraints imposing periodicity in a weak sense are added using Lagrange multipliers. In particular, the upper and lower energy bounds for the corresponding strongly periodic problem are produced, quantifying the accuracy of the weakly periodic boundary conditions. A numerical example demonstrates the evaluation of energy bounds and the performance of weakly periodic boundary conditions on a representative volume element. © 2016 The Authors. International Journal for Numerical Methods in Engineering Published by John Wiley & Sons Ltd


INTRODUCTION
In this paper, we consider the two-scale problem where a Stokes flow is present on the subscale and a Darcy flow on the macroscale. Thus, the macroscale is considered homogeneous, while the subscale is considered porous, that is, it consists of a solid and a fluid phase. The two-scale approach used to simulate the behavior of materials is common for a large number of engineering applications including oil geology [1], sintering [2], resin transfer modeling [3,4], and transportation of matter [5] to name a few.
The objective for models of this type, known as multiscale models, is to capture the subscale effects without the computational effort involved in resolving the complete microstructure. To that end, computational homogenization [6] is employed. Here, a representative volume element (RVE), containing a suitably large volume of the microstructure, is used as a replacement for the classic phenomenological model pertinent to the macroscale problem, that is, a Darcy flow. A Stokes flow is present on the RVE, which is driven by the macroscale pressure gradient. The resulting macroscale seepage is computed using volume averaging of the velocity on the RVE. In this work, the solid part of the porous material is considered rigid.
The derivation of a Darcy flow from a subscale Stokes flow can be found in, for example, [7][8][9]. Previous work concerning the upper and lower energy bounds of a Darcy flow can be found in [10] and [11] where various kinds of correlation functions, describing the characteristics of the microstructure, are used to estimate the bounds. In this paper, the bounds produced using weak periodicity [12], which allows for handling non-periodic meshes. 308 C. SANDSTRÖM AND F. LARSSON We consider the variationally consistent homogenization (VCH) method as a framework for the multiscale model. VCH is an extension of variational multiscale method (VMS) proposed by Hughes and coworkers [13,14]. The more classic format of VMS has also been used on a Stokes flow as shown in [15]. The extension of VMS to VCH was first proposed by Larsson et al. [16], where the pertinent generalized macrohomogeneity condition was introduced for guaranteeing energy equivalence for the homogenized response.
Following along the lines of Sandström et al. [12], we consider a viscous potential on a fully resolved domain consisting of a fluid phase contained within a rigid and geometrically periodic open pore solid. By splitting the domain into a finite number of subdomains and assuming separation of scales, we define the geometry and size of the pertinent RVE. From here, we use, as an approximation, periodic boundary conditions that, as shown in [17], satisfy the variationally consistent macrohomogeneity condition. In other words, we ensure energy equivalence on macroscale and subscale [16]. In particular, we impose periodic boundary conditions in weak form, resulting in additional unknown Lagrange multipliers on the RVE. The unknowns constitute external loads necessary to uphold periodicity on the pertinent fields.
By minimizing a viscous potential on the RVE, a saddle point problem is produced, allowing for the analysis of upper and lower energy bounds for the strongly periodic problem. The bounds are produced by confining pertinent solution spaces to contain only strongly periodic velocities or pressures. In a 2D setting, this can be achieved using quadratic interpolation over each opening on the boundary because the boundary is 1D [12]. In this way, we allow for the fluid to enter and exit the RVE while satisfying the no-slip condition on obstacles crossing the boundary. However, as the geometry of the surface of the RVE in 3D is of arbitrary 2D shape, the choice of global base functions is not as straightforward.
As a suggestion for the choice of basis functions on the 2D surface, we introduce the concept of solution-based shape functions. To compute the shape functions, we use a possibly non-periodic solution of a Stokes flow to produce an analytically periodic function. The function is then projected onto the boundary, producing a function that is periodic up to the error introduced by the discretization of the domain while satisfying no-slip conditions. We note that this procedure is used on the velocity to produce an upper bound and on the pressure to produce a lower bound. In the case of velocity, a method for compensating for possible compressibility artifacts due to the prescribed inflow on all open boundaries is also discussed. We note that an important part of this paper concerns the weak periodic boundary condition [12]. In [18], an alternative approach to weak boundary conditions is discussed.
The outline of the remainder of this paper is as follows. In Section 2, we start out from the energy formulation and derive the weak form of the problem. In Section 3, the upper and lower energy bounds are discussed and for the special case of linear flow, the permeability. In Section 4, a general method for producing function spaces pertinent to the bounds is shown along with special handling of incompressibility. The function spaces for the bounds are studied in more detail in Section 4.2 for the upper energy bounds and in Section 4.3, for the lower bound. Section 5 contains a numerical example of a simple RVE where the bounds are computed and the performance of the weak periodic boundary condition is evaluated. In Section 6, we finish the paper with some concluding remarks and suggestions for future work.

Multiscale Stokes flow -the saddle point problem
In this section, we provide the reader with the basic background on the two scale problem and the pertinent boundary conditions.
For the subsequent two-scale problem, we introduce the macroscale domain and the RVE . The boundary D @ of the macroscale domain is split into two parts, D v [ p , where p is the part where the pressure is prescribed and v is the part where the out-flux is prescribed. The RVE consists of two subdomains D F [ S where F is the volume containing the fluid and S the volume containing the solid phase, as shown in Figure 1. The solid phase is considered rigid, and the interface between the two phases is denoted int . F is the part of @ where fluid can enter and exit the domain. For reasons given momentarily, we also denote the part of F having a normal in a positive direction as FC and for part having a normal in the negative direction as F . Finally, we also define the boundary between FC and F as F0 D FC \ F as shown in Figure 1(c).
In order to establish a coupling between the macroscale and the subscale, we first define the pressure p on both scales. The pressure p on the RVE is then split into two parts, p D p M C p S , where p M is the smooth macroscale part and p S is the fluctuating subscale part. Furthermore, we introduce the macroscale pressure N p and the macroscale pressure gradient N g. As first-order homogenization is assumed, p M is allowed to vary linearly within the RVE as where x is the coordinate and N x F is the centroid of the fluid part of the RVE. From Equation (1), we conclude the macroscale-subscale coupling as N g D r p M . For a discussion concerning the uniqueness of the split, we refer to [17].
The macroscale problem is given on strong form as where N w is the seepage, n is the outward pointing normal to the boundary, O p is a prescribed pressure, and O v n is the prescribed flux in the outward pointing normal direction of v . Furthermore, ¹ º denotes implicit dependence on .
We define the seepage N w in terms of the stationary mean RVE potential ¹ N gº In the case of a laminar, incompressible, and strongly periodic flow, the potential is given as a stationary point to the optimization problem 310 C. SANDSTRÖM AND F. LARSSON where v is the velocity of a fluid particle, the subscale pressure p S is a Lagrange multiplier due to the incompressibility condition, andˆ.l/ is a viscous potential (expressed in terms of the velocity gradient l D v˝r ), such that @@ l D V being the deviatoric part of the Cauchy stress. The solution spaces in Equation (4) are defined as As a final remark on the potential, we note that t S D . V p S I/ n is the traction pertinent to the reaction forces due to split of into a finite number of RVEs. Henceforth, we will neglect this term as it is zero when the macrohomogeneity condition is fulfilled. For further reading on this topic, we refer to [17].
For the special case of linear flow, we have the viscous stress V defined as V D 2 OEv˝r sym (7)

Weak periodicity on the RVE using Lagrange multipliers
We will now weaken the conditions on v and p S by using the larger function spaces V and P S as replacements for V P and P S;P and instead impose periodicity on the pressure and velocity fields on the RVE in a weak sense. For this purpose, we introduce the jump operator as where x 2 FC and x .x/ are the corresponding points on F . When enforcing periodicity, we do so by manipulating the stationary mean RVE potential in Equation (4). Thus, we add the conditions where t SC and t S are tractions along FC and F , respectively. As the conditions enter into the optimization problem, they give rise to Lagrange multiplierˇfor the velocity constraint and for the pressure constraint. The stationary mean RVE potential now takes the form where is the mean RVE potential, defined as v; p S ;ˇ; ; N g For the situation in Equation (10), we note that the seepage satisfies due to the stationarity of (cf. [12]). Here, 0 denotes the directional derivative of w.r.t. N g in the direction of ı N g.

Remark 1
The solution to Equation (10) can be stated as if we assume a unique solution. Consequently, the explicit description of the stationary RVE problem becomes We stress that the existence of a unique solution v ¹ N gº ; p S ¹ N gº ;ˇ¹ N gº ; ¹ N gº is a stronger requirement than the formulation in Equation (10).

Remark 2
We note that by using the definition of the split and first-order homogenization, the periodicity condition on the subscale pressure p S is equivalent to the condition Thus, we can rewrite Equation (11) as Here, the split of the domain into RVEs gives the traction t D t S N g OEx x F n on the boundary. Assuming that the integral containing the subscale traction t S vanishes because of boundary conditions satisfying the variationally consistent macrohomogeneity condition, we obtain the last integral in Equation (17). This formulation allows for the computation of the full pressure up to a constant rather than the subscale pressure p S . 312 C. SANDSTRÖM AND F. LARSSON

Weak form
For the weak form of the macroscale equation, we refer to [12] and proceed by taking variations in all subscale quantities in Equation (11), producing the following weak form of the subscale problem: In Equation (18a), the known macroscale pressure gradient N g acts as a driving load for the flow.

Remark 3
In a discrete setting, the tangent stiffness of Equation (18) is a block matrix where submatrices C and D constitute the weak boundary conditions. In the case where relevant base functions are global, the submatrices can be considered dense. As a result, the computational cost increases rapidly as the Lagrange multiplier approximations are refined.

Confining the solution spaces
From Equation (10), we note that by confining the function spaces in certain ways, we can compute the upper and lower energy bounds. In the special case of a linear fluid material, we are also able to compute upper and lower bounds for the permeability. By replacing V with V 0 V and G with G 0 G , we produce an upper energy bound. More specifically, by choosing V 0 to contain only strongly periodic functions, we produce an upper bound for strongly periodic velocities. This also voids the supremum onˇ. The subspace G 0 contains the discretized , giving the inequality

313
By the same reasoning, we replace the function spaces P S and B with P S 0 and B 0 , where the particular choice of P S 0 only contains strongly periodic pressures. This yields the inequality pertinent to the lower energy bound for strong periodicity as (21) and (22) gives the bounds on strong periodicity as

Permeability and energy relation in the linear case
In the special case where we have a linear fluid, we have from [12] the following relations between the energy potential on the RVE and the permeability N K Here, the permeability N K is computed as where v .i / is the velocity field pertinent to a unit pressure gradient N g in the ith spatial direction. h i is the intrinsic averaging operator, defined as

Computing discrete shape functions
For the analysis of upper and lower bounds for strong periodicity according to Section 3, we require the velocity or pressure function spaces to contain only strongly periodic functions. In the case of a 2D flow, this can be achieved using global quadratic shape functions along the fluid part of the boundary (cf. [12]). However, in 3D, there is no apparent way of producing global shape functions for that purpose, as the arbitrary geometry is in 2D. To that end, we introduce the concept of solutionbased shape functions being global base functions that are produced by manipulating a possibly non-periodic solution of a Stokes flow. If these base functions are strongly periodic, we have ensured that the resulting flow is indeed also strongly periodic.
As the discussed method is used to create solution spaces for both velocity and pressures and because the procedure is the same for both, we introduce the loosely defined function spaces Q and Q 0 , which we can switch for its velocity or pressure counterparts later. Q D ® q W q is sufficiently regular and satisfies Dirichlet conditions on F ¯( 27a) We also introduce the function ' h 2 Q being either the discrete pressure or velocity field. Here, we restrict ourselves to use the same number of solution-based shape functions as the number of spatial dimensions of the problem, n dim .
For clarity, the following procedure is shown in the flowchart in Figure 2. As a starting point, we solve a Stokes flow on F where a unit pressure gradient N g D e i is applied. The index i D 1; : : : ; n dim indicates the spatial direction of the pressure gradient. The discrete solution to this problem is denoted ' .i/ h , being either the velocity solution v .i / or the pressure solution p S.i/ . Naturally, ' .i/ h must satisfy all Dirichlet conditions in order to qualify as a proper global base function.
From here, we use the discrete ' .i/ h (cf. Figure 3(a)) to produce an analytically periodic function Q ' .i/ 2 Q 0 (cf. Figure 3 where n is the number of mirror points for a point on F and x m j is an operator such that x m j .x/ is the j th mirror point of x. We note that the number of mirror points varies over the cubic domain. For instance, on a surface, there are two mirror points, one on FC and one on the opposite side F . On an edge of a 3D cube, there are four mirror points, and on a corner, there are eight mirror points, as shown in Figure 4.
By projecting the analytically periodic function Q ' .i/ onto the discretized boundary F , we create a quasi-periodic function Q ' .i/ h with an error in periodicity up to that introduced by the discretization of the boundary. More specifically, the value of Q ' .i/ h in coordinates of nodes is exactly Q ' .i/ and is otherwise interpolated using the base functions pertinent to the discretization for that particular function. That is, in the case of Taylor-Hood elements, if Q ' .i/ represents the velocity, it is interpolated using quadratic base functions, and if Q ' .i/ represents the pressure, it is interpolated using linear base functions.  By producing the solution-based shape functions Q v h;i and Q p h;i , we can now give the velocity and the pressure as and where a i 2 R and b i 2 R are arbitrary weights. Furthermore, we note that because both Q v h;i and Q p h;i tend to strong periodicity as the fineness of the mesh increases, so does both Q v h and Q p h . As a concluding remark, we note that as v in Equation (21) and p S in Equation (22) tend to strong periodicity, the respective bounds approach the strongly periodic energy content. Thus, we can use weak periodic boundary conditions on the RVE problem pertinent to the production of ' .i/ h in order to enhance the bound.

Special case of incompressibility
In the case where we choose to use solution-based shape functions on the velocity field, projecting the analytically periodic function Q v onto F may introduce compressibility into the model because of non-periodicity of the mesh. This implies a violation of the incompressibility condition. In order to compensate for such an error, we compute correction factors that are used to scale the prescribed velocity on FC and F . As a first step, we split the solution-based shape function Q v h into three parts as where Q v C h and Q v h are the parts of Q v h on FC and F , respectively, while Q v 0 h is the part on F0 . We can now introduce variables a C and a such that we weight the terms in Equation (31) as We want to find a C and a such that N Q v h is as close to Q v as possible while maintaining global incompressibility on F , that is, the same amount of fluid enters and exits the domain. Thus, we state the optimization problem as subject to: In weak form, we state the problem as follows: we can solve a C , a , and from the following system of equations: 2 For completeness, we give the explicit expressions for the unknowns in (35a)

Solution-based shape functions on the velocity
For the computation of the upper energy bound discussed in Section 3, we introduce the function space where b i is an unknown coefficient and N Q v h;i is the solution-based shape function arising from a pressure gradient in the direction e i . In order to avoid numerical difficulties relevant to possible compressibility, we define the velocity N Q v h;i on @ F as where the correction factors a C and a are defined in Equation (37). We note that, because N Q v h;i is strongly periodic, so is v 2 V 0 and c .v;ˇ/ D 0 for anyˇ2 B . The formulation of the problem is now as follows: Find .v; p S ; / 2 V 0 P S G such that a .vI ıv/ b .ıv; p S / D e .ıv; In the finite element setting, V 0 is simply constructed by defining the global basis functions pertaining to N Q v h;i on the nodes residing on F .

Solution-based shape functions for the pressure
For the lower energy bound, we follow the procedure described in Section 4.3 and introduce the function space where a i is an unknown coefficient and Q p h;i is the solution-based shape function pertinent to the ith spatial direction.
As the pressure is strongly periodic, we have d .p S ; / D 0 for any p S 2 P S 0 and 2 G . Thus, we state the problem on weak form, as follows: Find .v; p S ;ˇ/ 2 V P S 0 B such that a .vI ıv/ b .ıv; p S / c .ıv;ˇ/ D e .ıv; Similarly as for the velocities in Section 4.3, we note that the finite element counterpart of P S 0 can be constructed by constraining the nodal values of the pressure on F pertinent to Q p h;i .

Preliminaries
In the numerical example, we will investigate the permeability and its upper and lower bounds on a geometrically periodic structure. More specifically, we will investigate how the weak periodic boundary conditions, and the proposed bounds, perform compared with strong periodicity. First, we consider the approximations on fixed discretization (Section 5.2) and later also the behavior of the bounds during refinement of the finite element discretization (Section 5.3).
To this end, we will use a 'quasi-periodic' mesh, which simply means that the surface elements on opposite sides are of similar size. The permeability of the RVE with the quasi-periodic mesh and weak periodic boundary conditions will be compared with two alternative discretizations of the same RVE: (i) a mesh where surface elements on opposing sides are of very different sizes and (ii) a strongly periodic mesh where the sizes of surface elements are approximately the same as those of the quasi-periodic mesh.
As to the geometry of the subscale, we choose to use a body centered cubic (BCC) structure of slightly intersecting spheres as a solid, around which a Stokes flow is present. The material model used for the fluid is the standard linear model presented in Equation (7). The RVE is chosen as a unit cell as seen in Figure 5(a). The space occupied by the RVE is chosen in such a way that the (x; y) and the (x;´) planes are planes of symmetry and the RVE window is slightly shifted along the x-axis in order to create a non-symmetric geometry. The distance from the center sphere to the boundary of the RVE is 0.5 in the y and´directions and 0.46 in the x direction at its closest point (cf. Figure 5(b)). The distance between the centers of spheres in the same plane is 1 in the (x; y), (y;´), and (x;´) planes, respectively, and the radius of each spheres is 0.45. The solid phase in the RVE is shown in Figure 5(c).
The shape of the elements in all simulations is tetrahedral, and the discretization is P2-P1, that is, a Taylor-Hood element (quadratic velocity and linear pressure).
For the following numerical example, we define the function spaces for the discretization of the Lagrange multipliers as where n p is the order of the polynomial, b ij 2 R 2 and g ij 2 R are coefficients, and and Á are the local coordinates on the surface of the RVE. As the polynomial depends on two coordinates, note that the number of terms in the expression is ..n p C 1/ 2 C n p C 1/=2 if we consider all cross-terms. Furthermore, we choose solution spaces for the velocity and pressure fields according to Section 4.

Remark 4
In order to maintain compatibility between the velocity and pressure spaces, we avoid tetrahedral elements located such that all velocity degrees of freedom are either prescribed or hanging on a solution-based shape function. This is solved by flipping the edge between the tetrahedral element and its neighbor such that one node is located inside the volume.

Weak approximations and bounds on fixed meshes
We shall now investigate how the weak periodic boundary conditions and the pertinent bounds converge on fixed finite element discretizations (meshes). The behavior will be studied for fully non-  periodic as well as quasi-periodic meshes. Here, quasi-periodic means having elements of similar sizes on opposing sides of the RVE, whereas the fully non-periodic mesh has differently sized elements on opposing sides. For reference, we also consider a strongly periodic mesh with conventional periodic boundary conditions enforced.
As to the discretization of the domain, three meshes are used: (1) A 'quasi-periodic' mesh with approximately equally sized element on all sides (cf. Figure 7(a)). (2) A non-periodic mesh with a coarse mesh on the F and a fine mesh on FC (cf. Figure 7(b)).
(3) A periodic mesh with the same element sizes as in 1.
Detailed information on the number of nodes and elements for the various meshes is given in Table I.
In order to estimate how the weak periodic boundary conditions converge towards strong periodicity as we increase the order of the polynomial, we evaluate the upper and lower bounds for the permeability. Thus, if the gap between the upper and lower bounds is small, the weak periodicity will perform well. Figure 8 shows how the first eigenvalue K I of the permeability tensor changes for the quasi-periodic and non-periodic meshes. The same quantity for the strong periodic boundary conditions using the periodic mesh is shown in both plots. Note that the order of all polynomials is the same for each choice of n p , both when producing the shape functions and for imposing weak periodicity.
We note that the shape of the curves in both plots is similar but the values in Figure 8(a) are higher than that for the strongly periodic solution, while the values in Figure 8(b) are lower. These differences arise from the different approximations induced by the individual meshes (cf. the subsequent convergence study in Section 5.3). We see that K I for the strongly periodic solution is closer to K I for the quasi-periodic mesh with high n p as compared with the non-periodic one. This is expected because the periodic and quasi-periodic meshes are of similar size. Figure 9 shows the velocity field in a cross-section of the RVEs with weak periodic boundary conditions. The pressure fields for the same solutions are shown in Figure 10. Finally, for the non- where the upper bound does not decrease monotonically. The same quantity for the strongly periodic boundary conditions using the periodic mesh is shown in both plots. Again, we note the fact that the strongly periodic solution is not within all bounds due to different meshes used.

Convergence of bounds with mesh refinement
Finally, to demonstrate how the bounds behave as the mesh is refined, we compute K I for meshes with various choices of largest allowed length h of an element edge and varying orders of weak periodicity when computing the solution-based shape functions. We compare the resulting bounds with the strongly periodic solutions of a mesh of similar size. The results are shown in Figure 13. Note that the meshes are of similar size only; hence, there is no guarantee that the strong periodic solution should be contained within the bounds (as indeed is the case for some choices of h). However, as h decreases, the solution tends to a value within the bounds.

CONCLUSIONS AND OUTLOOK
In this paper, we have shown a novel approach on how to compute upper and lower energy bounds for flow through porous materials from a saddle point problem. In the case of a linear flow, upper and lower bounds on the permeability are computed, and with the use of these bounds, we have also evaluated the performance on weak periodicity and shown that they perform well, at least for higher-order approximations of the Lagrange multipliers.
The bounds were computed using global base functions for the velocity and pressure using solutions-based shape functions, pertaining from the solution of a Stokes flow that is then processed to produce a strongly periodic base function. Correction factors due to possible compressibility in the otherwise incompressible flow are also discussed.
In the numerical section, we have shown that the use of solution-based shape functions is suitable to use on both quasi-periodic as well as non-periodic meshes and that the bounds converge as the flux and traction required for weak periodicity increase are refined. This implies that we can replace the strongly periodic boundary conditions with weak ones of sufficiently high-order of approximation.
As a suggestion for future work using solutions-based shape functions, an evaluation of the use of this technique for the Dirichlet boundary condition for the Stokes flow in this type of application should be done. This would enable us to compute upper and lower bounds for a non-periodic microstructure.
An important issue in computations involving porous materials is that of the coupling between permeability and deformation; a framework for a variationally consistent model that takes this coupling into consideration is of great interest.