 Research Article
 Open Access
 Published:
Weakly periodic boundary conditions for the homogenization of flow in porous media
Advanced Modeling and Simulation in Engineering Sciences volume 1, Article number: 12 (2014)
Abstract
Background
Seepage in porous media is modeled as a Stokes flow in an open pore system contained in a rigid, impermeable and spatially periodic matrix. By homogenization, the problem is turned into a twoscale problem consisting of a Darcy type problem on the macroscale and a Stokes flow on the subscale.
Methods
The pertinent equations are derived by minimization of a potential and in order to satisfy the Variationally Consistent Macrohomogeneity Condition, Lagrange multipliers are used to impose periodicity on the subscale RVE. Special attention is given to the bounds produced by confining the solutions spaces of the subscale problem.
Results
In the numerical section, we choose to discretize the Lagrange multipliers as global polynomials along the boundary of the computational domain and investigate how the order of the polynomial influence the permeability of the RVE. Furthermore, we investigate how the size of the RVE affect its permeability for two types of domains.
Conclusions
The permeability of the RVE depends highly on the discretization of the Lagrange multipliers. However, the flow quickly converges towards strong periodicity as the multipliers are refined.
Background
We consider the classical problem of flow in porous media. On the macroscale, this phenomenon is often modeled as seepage governed by Darcy's law. Such seepage occur in a vast amount of natural as well as engineered materials, and applications include geomechanics, biomechanics and foam materials designed for energy absorption.
On the subscale, where the details of the pore system are resolved, Stokes' flow is an accurate description of the problem. In order to capture the effective properties of the subscale, the Stokes flow is solved on a Representative Volume Element (RVE) which should be large enough to represent the true subscale yet small enough to be as computationally efficient as possible [1]. In order to allow for the use of RVEs, the microstructure of the pertinent material should be ergodic and statistically spatially homogeneous. For further reading on the size of the RVE for homogenization of Stokes flow, we refer to [2].
Following the work by Sandström et al. [3],[4], we consider a twoscale problem where the subscale is represented as a Stokes flow on a strongly heterogeneous domain, consisting of a fluid within an open pore system. By adopting the concept of Variationally Consistent Homogenization [5], a macroscale representation in the form of a Darcy flow is produced. In the special case of linear flow, the homogenized tangent represents the permeability tensor whereas in the nonlinear case it serves as the consistent tangent in the macroscale Newton iterations.
In mechanics, homogenization is used to capture microstructural effects in a material subject to some load by deriving smooth effective properties on the structural scale. The model defined by the RVE can be used as a constitutive relation in itself, in a concurrent manner, or as a tool for calibrating an existing macroscale model. Several types of homogenization exists, such as asymptotic expansion which can be used to determine the macroscale properties in an analytical manner, see e.g [6][8]. In recent years, computational homogenization, where a local boundary value problem is solved on an RVE, has been subject to intense research, see e.g [9][12]. In the context of computational homogenization of porous materials, important areas of application are Resin Transfer Molding (RTM) [13],[14], oil geology [15], sintering [16] and transportation of matter [17].
Assuming separation of scales, we may adopt homogenization to derive the problem on two scales: the macroscale, representing the global structure, and the microscale, where the microstructure of the material is resolved. Classical homogenization concerns average theorems for the macroscale (effective) fluxes and primal variables (including possible gradients). Enforcing energy or work equivalence for the formulations on the different scales defines the socalled macrohomogeneity condition, cf. [9].
For single field problems, such as e.g. elasticity or heat conduction, there are three classical boundary conditions that satisfy macrohomogeneity: Dirichlet, Neumann and Periodic. However, in the case of a Stokes flow, it is not obvious how to choose the Dirichlet and Neumann conditions since there exists two primary fields of unknowns, namely velocity and pressure. Suggestions on how to choose Dirichlet and Neumann boundary conditions are given in [18]. It should also be mentioned that, in most cases, the periodic boundary condition performs better than Dirichlet and Neumann boundary conditions in terms of convergence with the size of the RVE [19]. In Sandström and Larsson [3], it was shown that periodic boundary conditions on the subscale (fluctuating) pressure and the (total) velocity defines a prolongation condition that satisfies the generalized macrohomogeneity presented in [5] thus ensuring no energy production on the subscale.
In this paper, we consider homogenization of the saddlepoint problem pertinent to the fully resolved Stokes’ problem within an open pore system. In contrast to the derivation by Sandström and Larsson [3], we thus carry out computational homogenization on pertinent potentials, rather than balance equations. We shall consider the particular choice of periodic boundary conditions, whereby the end result will be identical to that in [3]. However, we present this alternative derivation with the motivation that the arising subscale potentials will be utilized for computing upper and lower bounds on the effective properties, cf. below.
The classical approach in Finite Element Analysis of RVEs is to enforce periodic boundary conditions by treating two degrees of freedom on opposing sides of a domain as one single degree of freedom. Although computationally effective, this approach calls for a mesh which has identical discretization on either opposite side, which is a severe difficulty in 3D in the case of unstructured meshes.
The purpose of this work is to void the dependence on mesh periodicity for the periodic boundary conditions and instead impose periodicity in a weak sense, cf. [20],[21] where the momentum equation has been solved on the subscale for an elasticity problem. We note that in the former paper, the Lagrange multipliers are discretized with piecewise polynomials and in the latter, the displacement is interpolated by polynomials. With a minimization problem as point of departure, constraints pertinent to the boundary condition are added and as a result, Stokes flow with additional terms containing Lagrange multipliers is produced. The Lagrange multipliers can be identified as the required influx and traction necessary to maintain periodicity in pressure and velocity. From the minimization problem, bounds for the effective permeability are produced.
The remainder of this paper is organized as follows: In the Section “Methods”, the twoscale formulation of the saddlepoint problem pertinent to Stokes flow is derived in detail. Two numerical examples are presented in Section “Results and discussion”. The first example concerns an RVE in the form of a unit cell with a nonperiodic mesh. In the second example, we investigate how the size of the RVE affects the macroscale permeability. Finally, the conclusions and an outlook to future work are presented in Section “Conclusions”.
Methods
The single scale problem
Consider a fully resolved porous domain Ω=Ω^{F}⋃Ω^{S}, such as the one depicted in Figure 1(a). The domain consists of a topologically periodic substructure where Ω^{F} is the part of Ω occupied by the fluid phase and Ω^{S} the part occupied by the solid phase^{a}. The interface between the solid and fluid phases is denoted Ґ ^{int} and the part of Ґ where fluid can enter and exit the domain is denoted Ґ^{F}=∂ Ω^{F}\ Ґ^{int}=Ґ ⋂∂ Ω^{F}. The fluid part of the boundary Ґ^{F} is further divided into ${\mathrm{\u0490}}_{\mathrm{P}}^{\mathrm{F}}$ where the pressure p is prescribed and ${\mathrm{\u0490}}_{\mathrm{V}}^{\mathrm{F}}$ where the velocity v is prescribed. We hereby restrict ourselves to flows with low Reynolds numbers and purely viscous, incompressible fluids, whereby the fluid velocity field v can be found by minimizing the energy potential pertaining to a local viscous potential Ф(v⨂∇), defined such that $\frac{\mathrm{\partial \u0424}(\mathit{v}\u2a02\nabla )}{\partial \mathit{v}\u2a02\nabla}={\mathit{\sigma}}^{\mathrm{v}}$ where σ^{v} is the deviatoric part of the Cauchy stress ^{b}. Thus, we seek $\mathit{v}\u03f5\mathcal{V}$ that satisfies the constrained problem
where $\widehat{\mathit{t}}=\widehat{p}\mathit{n}$ is the prescribed pressure on the boundary ${\mathrm{\u0490}}_{\mathrm{P}}^{\mathrm{F}}\subset {\mathrm{\u0490}}^{\mathrm{F}}$ and is defined below. This can equivalently be written as the infsup problem
where
and p is a Lagrange multiplier resulting from the continuity condition. Note that due to the fact that σ^{v} is the deviatoric part of the Cauchy stress σ, p is interpreted as the pressure.
We proceed by splitting the domain into a finite number of n domains Ω_{□,i} such that $\Omega ={\bigcup}_{i=1}^{n}{\Omega}_{\square ,i}$ and such that each subdomain retains geometric periodicity (cf. Figure 1(a) and the periodic cutout in Figure 1(d)). By the choice of function spaces and , all functions $(\mathit{v},p)\u03f5\mathcal{V}\times \mathcal{P}$ is continuous on the whole Ω ^{F}. Rewriting Equation 2 as the sum of the energy contribution from all subdomains gives
In order to separate the macro and subscales features, we split the pressure term p into a smooth part ${p}^{\mathrm{M}}\u03f5{\mathcal{P}}^{\mathrm{M}}$ and a fluctuating part ${p}^{\mathrm{S}}\u03f5{\mathcal{P}}^{\mathrm{S}}$ such that p=p^{M}+p^{S} and $\mathcal{P}={\mathcal{P}}^{\mathrm{M}}+{\mathcal{P}}^{\mathrm{S}}$, ${\mathcal{P}}^{\mathrm{S}}$ being the hierachial complement to ${\mathcal{P}}^{\mathrm{M}}$. Integration by parts on p^{M} in the continuity constraint in Equation 4 yields
where n is the outward pointing normal. The boundary integral on the right hand side in Equation 5 vanish on all internal boundaries as v and p^{M} are continuous. Thus, after introducing the split in p, Equation 4 can be restated as
Furthermore, the last two terms can be rewritten as
where it is noted that the last term contains the prescribed velocity, ${\widehat{v}}_{n}=\mathit{v}\bullet \mathit{n}$. Under the assumption that $\mathit{t}=\widehat{p}\mathit{n}={p}^{\mathrm{M}}\mathit{n}$ the integral over ${\mathrm{\u0490}}_{\mathrm{P}}^{\mathrm{F}}$ vanishes and Equation 6 is equivalent to
Computational homogenization
Up to this point, nothing has changed since the original problem except the formulation. To proceed, we now assume separation of scales, i.e. that the subscale feature has a length scale much smaller than that of the macroscale. Furthermore, we also make the assumption that v and p^{S} are periodic over, and continuous inside, each ${\Omega}_{\square}^{\mathrm{F}}$, thus replacing the condition on continuity over the boundaries ${\mathrm{\u0490}}_{\square}^{\mathrm{F}}$. As an intermediate step, we note that by removing continuity over ${\mathrm{\u0490}}_{\square}^{\mathrm{F}}$, reaction forces arise, which eventually will contribute to the subsequent macrohomogeneity condition. In [3] it is shown that periodic boundary conditions satisfy the aforementioned condition. In order to impose periodicity (either in weak or strong form), we start out by following along the lines of [3] and split the subscale boundary Ґ _{□} into two parts; ${\mathrm{\u0490}}_{\square}={\mathrm{\u0490}}_{\square}^{+}\bigcup {\mathrm{\u0490}}_{}^{\square}$ where the +/ sign is the sign of the normal to that part of the boundary^{c}. Furthermore, we introduce the jump operator
where x is a point on ${\mathrm{\u0490}}_{\square}^{+}$ and x^{}(x) is the corresponding point on the opposite side of the RVE. The conditions for periodicity are given as
where t^{S+} and t^{S} are the subscale tractions along the edges ${\mathrm{\u0490}}_{\square}^{+}$ and ${\mathrm{\u0490}}_{\square}^{}$ respectively. By imposing the periodicity constraints in a weak sense, i.e. introducing the Lagrange multiplier β for the constraint [ v ]=0 and γ for the constraint [ p^{S}]=0 and allow the constraints to be fulfilled in average, rather than confining the respective solution spaces, we get
where
The Lagrange multiplier β can be interpreted as the traction needed to maintain periodicity on v and γ as the flux needed to maintain periodicity on p^{S}. The infimum on γ is further discussed in Remark 1. In order to allow for strong (essential) boundary conditions on ${\mathrm{\u0490}}_{\mathrm{P}}^{\mathrm{F}}$ in the subsequent macroscale problem, the function space ${\mathcal{P}}^{\mathrm{M}}$ is confined, replacing the former integral formulation of the condition.
Remark1.
In order to motivate the infimum on γ, consider the supremum of the term containing p^{S} (which already is a Lagrange multiplier) in Equation 8.
which, when adding the constraint [p^{S}]=0 becomes, locally,
or in weak form
We now introduce the total energy potential П which is split into n RVE potentials ${\u041f}_{i}^{\text{int}}$ which, in turn, can be expressed using the RVE mean potential ${\u041f}_{\square ,i}=\frac{{\mathrm{\u041f}}_{i}^{\text{int}}}{\u2502{{\Omega}_{\square}}_{,i}\u2502}$ as
Here, the RVE potential is given as
and the external load as
Furthermore, we note that by introducing separation of scales, i.e. for each coordinate $\stackrel{}{\mathit{x}}\u03f5\Omega $ there exist one RVE, thus, the RVE mean potential functions can be written as
where i is the number of the RVE occupyed by coordinate $\stackrel{}{\mathit{x}}$. Here, we define the RVE such that $\stackrel{}{\mathit{x}}$ is the centroid of Ω_{□}. By the assumption that the RVE is small compared to the macroscale, we identify │Ω_{□}_{,i}│ as a volume element on the macroscale and rewrite the sum in Equation 16 as an integral. It should be noted that the term │Ω_{□}_{,i}│ in the definition of the RVE mean potential is left unchanged during the transition from sum to integral as we are interested in the mean potential in the vicinity of $\stackrel{}{\mathit{x}}$. We give the RVE mean potential π_{□} on explicit integral form as
We proceed by writing the total potential on compact form as
Nested saddlepoint formulation
We proceed by introducing the macroscale potential function ψ(p^{M}) as
where we refer to Appendix "Commutativity of inf and sup" for a proof on the commutativity of the inf and sup operators. We now introduce the macroscale pressure $\stackrel{}{p}$ and the macroscale pressure gradient $\stackrel{}{\mathit{g}}$ as
and assume 1st order homogenization, i.e. the macroscale pressure p^{M} varies linearily inside the RVE. Thus, we have
where ${\stackrel{}{\mathit{x}}}^{\mathrm{F}}$ is the center of mass of ${\Omega}_{\square}^{\mathrm{F}}$. We can now express ψ{p^{M}} in terms of the macroscale pressure $\stackrel{}{p}$ as
where we have introduced the local macroscale potential
Weak form of the macroscale problem
Although this paper mainly focuses on the subscale problem, we choose to present the macroscale equation in order to achieve completeness. By taking the directional derivative of the the global macroscale potential, we produce
We now define the macroscale seepage as
where ϕ is the porosity defined as $\u2502{\Omega}_{\square}^{\mathrm{F}}\u2502/\u2502{\Omega}_{\square}\u2502$ and 〈•〉_{□} is the intrinsic averaging operator. We now recognize the weak form of the macroscale problem as that of finding all $\stackrel{}{p}\u03f5{\mathcal{P}}^{\mathrm{M}}$ such that
where P^{M} and P^{M,0} are the trial and test spaces respectively; now pertaining to the macroscale pressure $\stackrel{}{p}$.
The RVE problem
The local (subscale) problem for a given macroscale pressure gradient $\stackrel{}{g}$ is produced by seeking the stationary point for variations of subscale quantities in Equation 20. The problem is stated as: Find $\left(\mathit{v},{p}^{\mathrm{S}},\mathit{\beta},\gamma \right)\u03f5{\mathcal{V}}_{\square}\times {\mathcal{P}}_{\square}^{\mathrm{S}}\times {\mathcal{\beta}}_{\square}\times {\mathcal{G}}_{\square}$ such that
for all $\delta \mathit{v}\u03f5{\mathcal{V}}_{\square}$, $\delta {p}^{\text{S}}\u03f5{\mathcal{P}}_{\square}^{\mathrm{S}}$, $\delta \mathit{\beta}\u03f5{\mathcal{\beta}}_{\square}$ and $\mathrm{\delta \gamma}\u03f5{\mathcal{G}}_{\square}$, where
Homogenization of velocity and the macroscale tangent
From the definition of seepage in Equation 28, we produce the possibly nonlinear relation between the seepage and the macroscale pressure gradient by differentiation as
Remark 2.
Note that the minus sign on the positive definite permeability tensor $\stackrel{}{\mathit{K}}$ in Equation 32 is to ensure positive dissipation due to drag interaction between the solid and fluid phases.
From Equation 30, we see that the unit sensitivity field is given as
for all $\delta \mathit{v}\u03f5{\mathcal{V}}_{\square}$, $\delta {p}^{\text{S}}\u03f5{\mathcal{P}}_{\square}^{\mathrm{S}}$, $\delta \mathit{\beta}\u03f5{\mathcal{\beta}}_{\square}$ and $\mathrm{\delta \gamma}\u03f5{\mathcal{G}}_{\square}$ where a^{′} is the directional derivative of a. Following [3], we can express an arbitrary unit pressure gradient as
From here, we make an ansatz of the response dv as
Using the definition of seepage in Equation 28 on the above equation, we produce the relation between seepage and pressure gradient perturbations as
whereby the macroscale tangent is identified.
Bounds on effective properties for strong periodicity
According to Equation 26, an upper bound is produced by confining the function spaces ${\mathcal{V}}_{\square}$ and ${\mathcal{G}}_{\square}$. Furthermore, by choosing the function space ${\mathcal{V}}_{\square}$ in such a way that periodicity is always fulfilled, the supremum on β is void, producing the following inequality
Equivalently, a lower bound is produced by confining the spaces ${\mathcal{P}}_{\square}^{\mathrm{S}}$ and ${\mathcal{\beta}}_{\square}$
By combining Equations 37 and 38, we get
We shall now consider the special case of linear flow, defined by $\u0424(\mathit{v}\u2a02\nabla )=\frac{\mu}{2}{[\mathit{v}\u2a02\nabla ]}^{\text{sym}}:{[\mathit{v}\u2a02\nabla ]}^{\text{sym}}$. Assuming that v, p^{S}, β and γ satisfies Equation 30 for some $\stackrel{}{\mathit{g}}$, ${\psi}_{\square}\left(\phantom{\rule{0.3em}{0ex}}\stackrel{}{\mathit{g}}\right)$ is rendered stationary. Thus, the stationarity condition for Equation 30 is
In the case of a linear flow, choosing δ v=v in the stationarity condition, Equation 40 is given as
Inserting the stationary point into π_{□} and using 41, we see that the RVE mean potential is given as
Thus, by bounding ѱ_{□}, we have also bounded $\stackrel{}{\mathit{K}}$. More specifically, we may represent Equation 39, in terms of the permeability tensor as
where
Discretization of solutions spaces on the RVE boundary
As to the specific choice of solution spaces for the Lagrange multipliers we note that which is the most efficient depends on both the discretization and the geometry of the subscale domain. One example of a feasible discretization of the Lagrange multipliers is the one presented in [20] where the pertinent unknown functions are discretized on a mesh consisting of the union of all nodes on opposite sides of the domain. Here, however, we choose to discretize the Lagrange multipliers β and γ as global polynomials, i.e.
where n_{ p } are the polynomial order in the respective approximation, s is a parameterized coordinate along ${\mathrm{\u0490}}_{\square}^{+}$, b_{ i } and g_{ i } are the respective coefficients and l_{□} is the side length of the RVE.
For an upper bound of the energy, we choose ${\mathcal{V}}_{\square}$ such that the velocity is always periodic, removing the supremum on β.
where 𝜙_{ i } are basis functions for the N_{ v } velocity degrees of freedom a_{ i }. It should be noted that, if 𝜙_{ i } is represented in polynomial base, the constraint [ 𝜙_{ i }]=0 requires approximations of order higher than 1 in the case where obstacles cross the boundary of the RVE. The reason for this is simply that the no slip condition on the obstacle surface implies zero velocity on the RVE boundary if the velocity approximation is constant or linear. For the same reason, the velocity approximation is applied patchwise between obstacles along the boundary. In practice, we use global quadratic 1D element along the boundary as shown in the example in Figure 2 and make all nodes along the boundary hang on the global element. Furthermore, we connect all nodes located on a corner, i.e. N_{1} is a master and W_{1}, W_{6}, S_{1}, S_{2}, E_{1} and E_{6} its slaves. Finally, we connect opposite sides, i.e W_{2} is a slave to E_{2}, S_{2} to N_{2} etc.
For a lower bound on the energy, we choose ${\mathcal{P}}_{\square}^{\mathrm{S}}$ as
Thus, p^{S} is trivially periodic.
Results and discussion
In this section we present two numerical examples. The ambition of the first example is to investigates how the the order of a polynomial approximation of the Lagrange multiplier affect the solution and what order is required to reach convergence in terms of seepage, i.e. when the velocity field has converged to periodicity. This example is performed on a unit cell containing one single, circular, obstacle. The result on a nonperiodic mesh is compared the results from the corresponding problem with strong periodicity. Upper and lower bounds for the permeability are also presented. In the second numerical example, a quantitative convergence study is performed on the size effect on seepage of an RVE containing a set of random obstacles or periodic unit cells for a give order.
In all examples, the Stokes flow is solved using the Finite Element Method on triangular TaylorHood elements (linear pressure, quadratic velocity). The used fluid model is σ^{v}=μ l^{sym} where the viscosity μ is chosen as unity and l^{sym} is the symmetric velocity gradient.
All numerical simulations are performed using the open source software OOFEM [22].
Influence of polynomial order on permeability
The analysis in this section aim at evaluating how the order of the Lagrange multiplier approximation affect the periodicity of the solution and how the weak periodicity differ from strong periodicity. The simulations are performed on a unit cell containing a circular obstacle with radius 0.25 which is located at (0.26,0.5) in order to produce a nonperiodic mesh (see Figure 3(c)). As to the actual computation of the permeability $\stackrel{}{\mathit{K}}$, we use the method presented in [4]. In the following example, we compute the permeability on a domain using two different discretizations where each discretization has a set of subproblem as described below.
I Periodic mesh and strong periodicity
II Nonperiodic mesh

(a)
Upper bound (constant pressure, weakly periodic velocity)

(b)
Weak periodicity (weakly periodic pressure and velocity)

(c)
Lower bound (weakly periodic pressure, velocity is either constant or quadratic along ${\mathrm{\u0490}}_{\square}^{\mathrm{F}}$)
Since we aim at replicating the behavior of strong periodicity, I is used as a reference solution.
To compute the respective bounds of the permeability $\stackrel{}{\mathit{K}}$, we use the function spaces suggested in Equations 46 and 47 and choose ${\mathcal{G}}_{\square}$ and ${\mathcal{\beta}}_{\square}$ as in Equation 45. Figure 4 shows the first and second eigenvalues K_{ I } and K_{ II } (K_{ II }≤ K_{ I }) of $\stackrel{}{\mathit{K}}$ and their respective upper and lower bounds for orders ranging from 0 to 15. As the lower bound pertinent to the quadratic velocity profile gives significanly tighter bounds than the constant velocity profile, the later is omitted in the remaining parts of this paper. Since a unit cell in a periodic pattern is isotropic, K_{ I } and K_{ II } should tend to the same value as the solution approaches periodicity [3], which is indeed the case as shown in Figure 4(c). Note that none of the bounds reach isotropy, although the upper bound is significantly closer than the lower bound. According to the results, an order 4 is sufficient.
We choose the discretization of the unknown functions as proposed in Section "Discretization of solutions spaces on the RVE boundary" and note that in practice, the load is applied piecewise and in this case we have two sets of polynomials for each unknown function, one for the horizontal and one for the vertical part of the RVE boundary.
In the case where polynomials are used for the discretization of the Lagrange multipliers, the number of Gauss points, n_{ G }, needed to perform an exact integration of the integrals d_{□}(•,•) and c_{□}(•,•) are computed as
where n_{ f } is the order of the approximation of the pertinent field and n_{ p } is the order of the Lagrange multiplier approximation. For a TaylorHood element, n_{ f }=1 for the pressure and n_{ f }=2 for the velocity.
As an indicator of how close to strong periodicity the fields are, we compute the L^{2} norm of the error as
where [•] represents the jump of • on ${\mathrm{\u0490}}_{\square}^{\mathrm{F}}$. As can be seen in Figure 5, both velocities converge quickly compared to the pressure p^{S}. Indeed, since the pressure is discretized by piecewise linear polynomials, true periodicity can never be achieved on a nonperiodic mesh, less in the cases of linear or constant pressure along ${\mathrm{\u0490}}_{\square}^{\text{F+}}$ and ${\mathrm{\u0490}}_{\square}^{\text{F}}$. The same hold for a quadratic discretization but due to the larger number of degrees of freedom, a solution closer to periodicity is achieved. We would, however, like to point out that as the main interest lies in computing the effective permeability and/or seepage, we consider convergence in terms of the mean velocity.
Figure 6 shows the pressure p^{S} along the vertical parts of the RVE along with the pertinent Lagrange multiplier for orders 0, 2, 4 and an overkill solution with polynomials of order 15. The effect of linear elements can be seen in these graphs, as even in the overkill solution, relatively large differences between the pressure functions on either side of the RVE are present. However, the large values of the Lagrange multiplier in the corners of the overkill solution suggests that the pressure field is close to periodicity in that area.
Comparing the pressure curves in Figure 6 with the corresponding curves for the velocity u in Figure 7 we again notice that the velocity is closer to periodicity at order 4.
Figure 8 shows the velocity u along the vertical part of the boundary. Notice the even functions in the Lagrange multipliers γ and β_{1} and the odd functions in β_{2} due to the symmetric shape of the RVE.
Impact of RVE size on permeability
When imitating a material using homogenization on RVEs, two main sources of errors are introduced; boundary conditions and the statistical representation of the microstructure. If the microstructure is truly periodic, the periodic type boundary condition introduces no error, but this is often an approximation of the materials subscale geometry. However, as this paper aims at producing periodicity in a weak sense, we assume that the true solution is indeed periodic. In this case, the error introduced by boundary conditions are the order of the approximation of the Lagrange multipliers as a perfectly periodic solution is not guaranteed. As to the error introduced in terms of statistical representation of the subscale, this can be overcome by either increasing the number of RVEs or increase the size such that all geometrical effects are captured in one RVE. In fact, the relative error introduced by the order of approximation can also be decreased by increasing the size of the RVE.
In order to study the influence of RVE size on the effective permeability, we choose to perform homogenization on two types of domain:
I Ω_{□} contains a perfectly aligned, circular obstacles where no obstacles cross the boundary

(a)
Weak periodicity on pressure and velocity. Order of approximation is 0.
II Ω_{□} contains pseudorandomly placed, circular obstacles which can cross the boundary

(a)
Weak periodicity on pressure and velocity. Order of approximation is 0.

(b)
Weak periodicity on pressure and velocity. Order of approximation is 4 and the load is applied piecewise.
For type II domains, the RVE is generated according to Algorithm 1.
where L_{□} is the integer length of one side of a rectangular RVE and ∆ set to 1% of the radius of the obstacles. Furthermore, we choose the function f(i,j) as
where Ф is a normally distributed random variable, μ is the mean and σ is the standard deviation. Here, we choose μ=0 and σ=0.2. Note that lines 911 in Algorithm 1 implies geometric periodicity and constant porosity.
We use the function spaces suggested in Equations 46 and 47 when computing the bounds of the permeability. In cases where obstacles cross the boundary, the load pertinent to the periodicity condition is applied piecewise, thus, the solution space of the Lagrange multiplier approximation becomes larger.
It should be noted that the highest possible order of the polynomial approximation is limited by the number of boundary elements subject to that constraint. In the case of randomly placed obstacles, it is possible that one element only, separates two obstacles. If that is the case, depending on the discretization of the pertinent function and Lagrange polynomial, the subscale tangent becomes singular. In such cases, a simple rule is used; the number of unknowns in the Lagrange multiplier cannot exceed the number of unknowns belonging to the pertinent function, eg. if only one linear element is present, a linear approximation is used. This is taken into account when producing the RVEs. However, this situation rarely occurs.
In order to produce a reliable estimate for the permeability of an RVE containing pseudorandomly placed obstacles, a sufficiently large number of RVEs is generated and the permeability computed according to [4]. In the case of an RVE with one periodically repeating unit cell, one realization of each size of the RVE suffices. Examples of meshed RVEs with periodic unit cells and random obstacles are shown in Figure 9.
Figure 10 shows how the permeability changes as the size of the RVE increase for RVE types I (a) and II (a) along with their respective bounds. From Figure 10(a) we can see that the lower bound differs from the weak periodic solution while the upper bound coincides. Comparing to the results in Section "Influence of polynomial order on permeability", this is true for 0th order approximation, but as the order increase, the permeability of the weakly periodic solution decrease. We also emphasize that the error introduced by the 0th order approximation increase, the error in the homogenized result decrease as the RVE grows, since the solution approach the strongly periodic solution. The differences between the solutions are further illustrated in Figure 11 where a pressure gradient in the x direction is imposed on an RVE containing 3×3 periodic unit cells. The image shows the magnitude of the velocity field v. The solutions are similar inside the RVE but differs on the boundaries.
Figure 10(b) shows how μ{K_{ I }} behaves as the size of the RVE increases. As in the previous case, the upper bound and the weakly periodic solution produce similar solutions whereas the lower bound yields a significantly lower permeability. By increasing the size of the RVE while keeping the order of the approximation fixed, the error introduced by the approximation increase while the total error decrease. The bump at RVE sizes 2 and 3 are due to two mechanisms; The permeability increase in the direction of the pressure gradient, if the distance between the obstacle orthogonal to the pressure gradient, increase (the volume fluid passing through per second increases quadratically with the distance) and the permeability decrease as the distance parallel to the pressure gradient increase (as this implies a longer distance). The first mechanism is dominant for small RVEs but as the size increase, the two evens out.
In the final example, shown in Figure 12(a), the order of the approximating polynomial has been increased to 4 and the load pertinent to the weak periodic boundary condition is applied patchwise. In order to allow for proper comparison of the two last examples, the pseudorandom domains used in the last example are the same as in the previous. It is apparent that the permeability is dependent on the order of the approximation, even for large RVEs.
A comparison between Figures 13(a) and 14(a) illustrates the impact of additional terms and piecewise application of the loads on the periodicity, especially on the north and south edges of Ω _{□}.
To conclude the numerical section, we note that the lower bound is closer to the strongly periodic solution in all examples. Furthermore, we also note that by increasing the resolution of the Lagrange multiplier approximations, the solution approaches strong periodicity fast. The cost of the enriched approximations are the additional degrees of freedom and a stiffness matrix with dense sub matrices pertinent to the boundary integrals in Equation 33.
Conclusions
In this paper, we produce a multiscale problem for a Stokes flow by minimization of an energy potential. During the minimization process, a split in the pressure term is introduced, after which the weak form of the problem is introduced by variations of the pertinent quantities. The result is a Stokes flow subscale problem and a Darcy flow macroscale problem.
In order to satisfy the macrohomogeneity condition on a non periodic mesh, weak periodic boundary conditions are imposed using Lagrange multipliers. These are of two types: unknown tractions maintain periodicity on the velocity and unknown fluxes maintain periodicity on the pressure. Due to the saddlepointnature of the problem, bounds on the macroscale permeability are produced by confining the pertinent function spaces.
The numerical examples have shown the rapid convergence of periodicity using polynomial approximation of relevant Lagrange multipliers. Furthermore, the expected asymptotic convergence of macroscale permeability due to RVE size has been verified.
Concerning future developments, the primary goal is to be able to couple permeability with deformation of the porous material. Since a 2D representation of an open pore system is unable to carry static load when deformed, it is necessary to extend the study to a more relistic 3D representation.
Endnotes
^{a} As this work only concerns the fluid phase, the solid phase is considered rigid and is modeled as impermeable obstacles in the Ω domain.
^{b} In the case of linear flow, $\u0424(\mathit{v}\u2a02\nabla )=\frac{\mu}{2}{\left[\mathit{v}\u2a02\nabla \right]}^{\text{sym}}:{\left[\mathit{v}\u2a02\nabla \right]}^{\text{sym}}$.
^{c} It is assumed that all RVEs have parallel edges/surfaces.
Appendix
Commutativity of inf and sup
For the subsequent proof, we define the potential
From [23], we have the maxmin inequality as
which for a saddle point implies equality. What remains to be shown is the validity of the right hand side of the equation, i.e that the term Θ(p^{M}) is concave in p^{M} for all v. For future use, we note that a variation in v yields, by stationarity
and a subsequent perturbation in v
To show concavity of Θ, we make a perturbation in p^{M} and choose δ v=dv in Equation 53, which yields
By yet another perturbation and the use of Equation 54
if Θ^{′′} is positive definite and b satisfies the required, classic infsup condition
Authors’ contributions
All authors have carried out the theoretical parts described in this paper. CS has carried out the necessary software development and has had the major responsibility for preparing the manuscript. All authors have read and approved the final manuscript.
References
 1.
Zohdi TI, Wriggers P: An introduction to computational micromechanics. Springer, Berlin Heidelberg; 2008.
 2.
Du X, OstojaStarzewski M: On the size of representative volume element for Darcy law in random media. Proc R Soc A: Math Phys Eng Sci 2006,462(2074):2949–2963. 10.1098/rspa.2006.1704
 3.
Sandström C, Larsson F: Variationally consistent homogenization of stokes flow in porous media. Int J Multiscale Comput Eng 2013,11(2):117–138. 10.1615/IntJMultCompEng.2012004069
 4.
Sandström C, Larsson F, Runesson K, Johansson H: A twoscale finite element formulation of Stokes flow in porous media. Comput Methods Appl Mech Eng 2013, 261–262: 96–104. 10.1016/j.cma.2013.03.025
 5.
Larsson F, Runesson K, Su F: Variationally consistent computational homogenization of transient heat flow. Int J Numerical Methods Eng 2009,81(13):1659–1686.
 6.
SánchezPalencia E: Nonhomogeneous Media and Vibration Theory. Springer, Berlin Heidelberg; 1980.
 7.
Allaire G: Homogenization of the Stokes flow in a connected porous medium. Asymptotic Anal 1989, 2: 203–222.
 8.
Hornung U: Homogenization and porous media. Springer, Berlin Heidelberg; 1997.
 9.
NematNasser S, Lori M, Datta SK (1996) Micromechanics: overall properties of heterogeneous materials. J Appl Mech 63(2): 561. NematNasser S, Lori M, Datta SK (1996) Micromechanics: overall properties of heterogeneous materials. J Appl Mech 63(2): 561.
 10.
Geers MGD, Kouznetsova VG, Brekelmans WAM: Multiscale computational homogenization: trends and challenges. J Comput Appl Math 2010,234(7):2175–2182. 10.1016/j.cam.2009.08.077
 11.
Larsson F, Runesson K: RVE computations with error control and adaptivity: the power of duality. Comput Mech 2006,39(5):647–661. 10.1007/s004660060108z
 12.
Löhnert S, Wriggers P: Homogenisation of microheterogeneous materials considering interfacial delamination at finite strains. Technische Mechanik 2003,23(24):167–177.
 13.
Ngo N, Tamma K: Microscale permeability predictions of porous fibrous media. Int J Heat Mass Transf 2001, 44: 3135–3145. 10.1016/S00179310(00)003355
 14.
Pillai KM, Advani SG: Numerical and analytical study to estimate the effect of two length scales upon the permeability of a fibrous porous medium. Transp Porous Media 1995,21(1):1–17. 10.1007/BF00615332
 15.
Arbogast T, Lehr H: Homogenization of a DarcyStokes system modeling vuggy porous media. Comput Geosci 2006, 78712: 1–18.
 16.
Ohman M, Runesson K, Larsson F: Computational Mesoscale modeling and homogenization of liquidphase sintering of particle agglomerates. Technische Mechanik 2011, 32: 463–483.
 17.
Nilenius F, Larsson F: Macroscopic diffusivity in concrete determined by computational homogenization. Int J Numerical Anal Methods Geomech (May 2012) 2012,37(11):1535–1551. 10.1002/nag.2097
 18.
OstojaStarzewski M (2007) Microstructural randomness and scaling in mechanics of materials. Chapman & Hall/CRC, Boca Raton, FL. OstojaStarzewski M (2007) Microstructural randomness and scaling in mechanics of materials. Chapman & Hall/CRC, Boca Raton, FL.
 19.
Yue X: The local microscale problem in the multiscale modeling of strongly heterogeneous media: Effects of boundary conditions and cell size. J Comput Phys 2007,222(2):556–572. 10.1016/j.jcp.2006.07.034
 20.
Larsson F, Runesson K, Saroukhani S, Vafadari R: Computational homogenization based on a weak format of microperiodicity for RVEproblems. Comput Methods Appl Mech Eng 2011,200(14):11–26. 10.1016/j.cma.2010.06.023
 21.
Nguyen VD, Béchet E, Geuzaine C, Noels L: Imposing periodic boundary condition on arbitrary meshes by polynomial interpolation. Electrical Eng 2011,55(November):390–406.
 22.
Patzák B, Bittnar Z: Design of object oriented finite element code. Adv Eng Softw 2001, 32: 759–767. 10.1016/S09659978(01)000278
 23.
Boyd S, Vandenberghe L (2010) Convex Optimization. vol. 25. Cambridge University Press. pp 487487. Chap. 1,10,11.
Acknowledgements
The project is financed by the Swedish Research Council (www.vr.se) under contract 20115388. The authors would also like to thank Håkan Johansson for fruitful discussions.
Author information
Affiliations
Corresponding author
Additional information
Competing interests
The authors declare that they have no competing interests.
An erratum to this article is available at http://dx.doi.org/10.1186/s403230150024x.
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 2.0 International License ( https://creativecommons.org/licenses/by/2.0 ), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
About this article
Cite this article
Sandstöm, C., Larsson, F. & Runesson, K. Weakly periodic boundary conditions for the homogenization of flow in porous media. Adv. Model. and Simul. in Eng. Sci. 1, 12 (2014). https://doi.org/10.1186/s4032301400126
Received:
Accepted:
Published:
Keywords
 Multiscale modeling
 Computational homogenization
 Stokes flow
 Weak periodicity
 Porous media