Skip to main content
  • Research Article
  • Open access
  • Published:

Weakly periodic boundary conditions for the homogenization of flow in porous media

An Erratum to this article was published on 30 April 2015

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 two-scale 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 two-scale 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 non-linear 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 so-called 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 saddle-point 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 in-flux 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 two-scale formulation of the saddle-point 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 non-periodic 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 phasea. 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 Ґ P F where the pressure p is prescribed and Ґ V 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 ∂Ф ( v ) v = σ v where σv is the deviatoric part of the Cauchy stress b. Thus, we seek vϵV that satisfies the constrained problem

minimize Ω F Ф ( v ) d V - Ґ P F t ^ v d S
(1)
subject to: v = 0 on Ω F
(2)

where t ^ =- p ^ n is the prescribed pressure on the boundary Ґ P F Ґ F and is defined below. This can equivalently be written as the inf-sup problem

inf v ϵ V sup p ϵ P Ω F Ф v d V - Ω F p ( v ) d V - Ґ P F t ^ v d S
(3)

where

V = v ϵ H 1 Ω F 3 : v = 0 on Ґ int , v = v ^ n n on Ґ V F
(4)
P = p ϵ L 2 Ω F
(5)

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.

Figure 1
figure 1

Topologically periodic substructure. (a) Fully resolved domain. (b) Division of fully resolved domain. (c) Homogenized macroscale domain. (d) Representative Volume Element.

We proceed by splitting the domain into a finite number of n domains Ω□,i such that Ω= i = 1 n Ω , 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 (v,p)ϵV×P is continuous on the whole Ω F. Rewriting Equation 2 as the sum of the energy contribution from all subdomains gives

inf v ϵ V sup p ϵ P Σ i = 1 n Ω , i F Ф ( v ) d V - Ω , i F p ( v ) d V - Ґ P F t ^ v d S
(6)

In order to separate the macro and subscales features, we split the pressure term p into a smooth part p M ϵ P M and a fluctuating part p S ϵ P S such that p=pM+pS and P= P M + P S , P S being the hierachial complement to P M . Integration by parts on pM in the continuity constraint in Equation 4 yields

Σ i = 1 n Ω , i F p M ( v ) d V = Σ i = 1 n - Ω , i F p M v d V + Ґ , i F n v p M d S
(7)

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 pM are continuous. Thus, after introducing the split in p, Equation 4 can be restated as

inf v ϵ V sup p M ϵ P M p S ϵ P S Σ i = 1 n Ω , i F Ф ( v ) d V - Ω , i F p S ( v ) d V + Ω , i F p M v d V - Ґ P F t ^ v d S - Ґ F n v p M d S
(8)

Furthermore, the last two terms can be rewritten as

- Ґ P F t ^ v d S - Ґ F n v p M d S = - Ґ P F t ^ v + n v p M d S - Ґ V F n v p M d S
(9)

where it is noted that the last term contains the prescribed velocity, v ^ n =vn. Under the assumption that t=- p ^ n=- p M n the integral over Ґ P F vanishes and Equation 6 is equivalent to

inf v ϵ V sup p M ϵ P M p S ϵ P S Σ i = 1 n Ω , i F Ф ( v ) d V - Ω , i F p S ( v ) d V + Ω , i F p M v d V - Ґ V F v ^ n p M d S
(10)

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 pS are periodic over, and continuous inside, each Ω F , thus replacing the condition on continuity over the boundaries Ґ F . As an intermediate step, we note that by removing continuity over Ґ 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; Ґ = Ґ + Ґ - where the +/- sign is the sign of the normal to that part of the boundaryc. Furthermore, we introduce the jump operator

where x is a point on Ґ + and x-(x) is the corresponding point on the opposite side of the RVE. The conditions for periodicity are given as

where tS+ and tS- are the subscale tractions along the edges Ґ + and Ґ - 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 [ pS]=0 and allow the constraints to be fulfilled in average, rather than confining the respective solution spaces, we get

where

V = v ϵ H 1 Ω F 3 : v = 0 on Ґ int , v = v ^ n n on Ґ V F
(11)
P S = p ϵ H 1 Ω F
(12)
P M = p ϵ H 1 Ω F : p = p ^ on Ґ P F
(13)
B = β ϵ L 2 Ґ F 3
(14)
G = γ ϵ L 2 Ґ F
(15)

The Lagrange multiplier β can be interpreted as the traction needed to maintain periodicity on v and γ as the flux needed to maintain periodicity on pS. The infimum on γ is further discussed in Remark 1. In order to allow for strong (essential) boundary conditions on Ґ P F in the subsequent macroscale problem, the function space P 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 pS (which already is a Lagrange multiplier) in Equation 8.

sup p S ϵ P S Σ i = 1 n - Ω , i F p S v d V
(16)

which, when adding the constraint [pS]=0 becomes, locally,

or in weak form

We now introduce the total energy potential П which is split into n RVE potentials П i int which, in turn, can be expressed using the RVE mean potential П , i = П i int Ω , i as

П p M , v , p S , β , γ = Σ i = 1 n П i int p M , v , p S , β , γ - П ext p M = Σ i = 1 n π , i p M , v , p S , β , γ Ω , i - П ext p M
(17)

Here, the RVE potential is given as

and the external load as

П ext p M = Ґ V F v ^ n p M d S
(18)

Furthermore, we note that by introducing separation of scales, i.e. for each coordinate x - ϵΩ there exist one RVE, thus, the RVE mean potential functions can be written as

π p M , v , p S , β , γ , x - = π , i p M , v , p S , β , γ
(19)

where i is the number of the RVE occupyed by coordinate x - . Here, we define the RVE such that 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 x - . We give the RVE mean potential π on explicit integral form as

We proceed by writing the total potential on compact form as

П p M , v , p S , β , γ = Ω π p M , v , p S , β , γ , x - d V - П ext p M
(20)

Nested saddle-point formulation

We proceed by introducing the macroscale potential function ψ(pM) as

inf v ϵ V sup p M ϵ P M p S ϵ P S β ϵ B inf γ ϵ G П p M , v , p S , β , γ = sup p M ϵ P M inf v ϵ V sup p S ϵ P S β ϵ B inf γ ϵ G П p M , v , p S , β , γ = sup p M ϵ P M ψ p M
(21)

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 p - and the macroscale pressure gradient g - as

p - = def p g - = def p -
(22)

and assume 1st order homogenization, i.e. the macroscale pressure pM varies linearily inside the RVE. Thus, we have

p M = p - ( x - ) + g - x - x - x - F
(23)

where x - F is the center of mass of Ω F . We can now express ψ{pM} in terms of the macroscale pressure p - as

ψ p - = Ω 𝜓 g - d Ω - П ext p -
(24)

where we have introduced the local macroscale potential

𝜓 g - : = inf v ϵ V sup p S ϵ P S β ϵ B inf γ ϵ G π g - , v , p S , β , γ
(25)

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

ѱ p - = Ω ѱ g - δ g - d V - Ґ V F v ^ n δ p - d S = 0
(26)

We now define the macroscale seepage as

w - = ѱ g - = 1 Ω Ω F v d V = ϕ v
(27)

where ϕ is the porosity defined as Ω F / Ω and 〈•〉 is the intrinsic averaging operator. We now recognize the weak form of the macroscale problem as that of finding all p - ϵ P M such that

Ω w - δ g - d V - Ґ V F v ^ n δ p - d S = 0 δ p - ϵ P M,0
(28)

where PM and PM,0 are the trial and test spaces respectively; now pertaining to the macroscale pressure p - .

The RVE problem

The local (subscale) problem for a given macroscale pressure gradient g - is produced by seeking the stationary point for variations of subscale quantities in Equation 20. The problem is stated as: Find v , p S , β , γ ϵ V × P S × β × G such that

a v ; δ v - b δ v , p S - c δ v , β = - e δ v , g -
(29)
- b v , δ p S - d δ p S , γ = 0
(30)
- c v , δ β = 0
(31)
- d p S , δγ = 0
(32)

for all δvϵ V , δ p S ϵ P S , δβϵ β and δγϵ G , where

Homogenization of velocity and the macroscale tangent

From the definition of seepage in Equation 28, we produce the possibly non-linear relation between the seepage and the macroscale pressure gradient by differentiation as

w - = w - g - d w - = d w - d g - d g - = - K - g - d g -
(33)

Remark 2.

Note that the minus sign on the positive definite permeability tensor 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

a v ; δ v , d v - b δ v , d p S - c δ v , d β = - e δ v , e i
(34)
- b d v , δ p S - d δ p S , d γ = 0
(35)
- c d v , δ β = 0
(36)
- d d p S , δγ = 0
(37)

for all δvϵ V , δ p S ϵ P S , δβϵ β and δγϵ G where a is the directional derivative of a. Following [3], we can express an arbitrary unit pressure gradient as

d g - = Σ i = 1 n dim e i e i d g ^
(38)

From here, we make an ansatz of the response dv as

d v = Σ i = 1 n dim v ( i ) e i d g ^ = Σ i = 1 n dim v ( i ) e i d g ^
(39)

Using the definition of seepage in Equation 28 on the above equation, we produce the relation between seepage and pressure gradient perturbations as

d w - = 𝜙 d v = 𝜙 Σ i = 1 n dim v ( i ) e i d g - = ϕ Σ i = 1 n dim v ( i ) e i d g - = - K - d g -
(40)

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 V and G . Furthermore, by choosing the function space V in such a way that periodicity is always fulfilled, the supremum on β is void, producing the following inequality

ψ g - = inf v ϵ V sup p S ϵ P S β ϵ B inf γ G π g - , v , p S , β , γ inf v ϵ V sup p S ϵ P S inf γ ϵ G π g - , v , p S , 0 , γ = ψ fU g -
(41)

Equivalently, a lower bound is produced by confining the spaces P S and β

ψ g - = inf v ϵ V sup p S ϵ P S β ϵ B inf γ ϵ G π g - , v , p S , β , γ inf v ϵ V sup p S ϵ P S β ϵ B π g - , v , p S , β , 0 = ψ L g -
(42)

By combining Equations 37 and 38, we get

ψ L g - ψ g - ψ U g -
(43)

We shall now consider the special case of linear flow, defined by Ф(v )= μ 2 [ v ] sym : [ v ] sym . Assuming that v, pS, β and γ satisfies Equation 30 for some g - , ψ ( g - ) is rendered stationary. Thus, the stationarity condition for Equation 30 is

a v ; δ v = - e δ v , g - δ v ϵ V
(44)

In the case of a linear flow, choosing δ v=v in the stationarity condition, Equation 40 is given as

Ω F μ v sym : [ v ] sym d V = - Ω F v d V p -
(45)

Inserting the stationary point into π and using 41, we see that the RVE mean potential is given as

1 Ω F Ω F μ 2 [ v ] sym : [ v ] sym d V + 1 Ω F Ω F g - v d V = 1 2 w - g - = - 1 2 g - K - g -
(46)

Thus, by bounding ѱ, we have also bounded K - . More specifically, we may represent Equation 39, in terms of the permeability tensor as

g - K - L g - g - K - g - g - K - U g -
(47)

where

ψ L g - = - 1 2 g - K - U g -
(48)
ψ g - = - 1 2 g - K - g -
(49)
𝜓 U g - = - 1 2 g - K - L g -
(50)

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.

B = β ϵ R 2 : β = Σ i = 0 n p b i s i l , G = γ ϵ R : γ = Σ i = 0 n p g i s i l
(51)

where n p are the polynomial order in the respective approximation, s is a parameterized coordinate along Ґ + , 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 V 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. N1 is a master and W1, W6, S1, S2, E1 and E6 its slaves. Finally, we connect opposite sides, i.e W2 is a slave to E2, S2 to N2 etc.

Figure 2
figure 2

Example of boundary elements approximating the velocity.

For a lower bound on the energy, we choose P S as

P S = p S ϵ P S : p S = 0 on Ґ F P S
(52)

Thus, pS 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 non-periodic 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 Taylor-Hood elements (linear pressure, quadratic velocity). The used fluid model is σv=μ lsym where the viscosity μ is chosen as unity and lsym 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 non-periodic mesh (see Figure 3(c)). As to the actual computation of the permeability 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 Non-periodic mesh

  1. (a)

    Upper bound (constant pressure, weakly periodic velocity)

  2. (b)

    Weak periodicity (weakly periodic pressure and velocity)

  3. (c)

    Lower bound (weakly periodic pressure, velocity is either constant or quadratic along Ґ F )

Figure 3
figure 3

Figures shows the unitcell (a) domain (b) periodic mesh and (c) mesh used during computations.

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 K - , we use the function spaces suggested in Equations 46 and 47 and choose G and β as in Equation 45. Figure 4 shows the first and second eigenvalues K I and K II (K II K I ) of 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.

Figure 4
figure 4

The Figure shows (a) the first eigenvalue of the permeability tensor K I versus n p (b) K II versus n p (c) K I K II versus n p .

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

n G = n p + n f + 1 2
(53)

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 Taylor-Hood 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 L2 norm of the error as

where [] represents the jump of • on Ґ F . As can be seen in Figure 5, both velocities converge quickly compared to the pressure pS. Indeed, since the pressure is discretized by piecewise linear polynomials, true periodicity can never be achieved on a non-periodic mesh, less in the cases of linear or constant pressure along Ґ F+ and Ґ 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 5
figure 5

The error e 2 versus n p .

Figure 6 shows the pressure pS 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.

Figure 6
figure 6

Pressure p S and flux γ along the vertical parts of Ґ F corresponding to Lagrange multipliers of order (a) order 0 (b) order 2 (c) order 4 (d) order 15.

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 7
figure 7

Velocity v and β 1 along the vertical parts of Ґ F corresponding to Lagrange multipliers of order (a) order 0 (b) order 2 (c) order 4 (d) order 15.

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.

Figure 8
figure 8

Velocity v and β 2 along the vertical parts of Ґ F corresponding to Lagrange multipliers of order (a) order 0 (b) order 2 (c) order 4 (d) order 15.

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

  1. (a)

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

II Ω contains pseudo-randomly placed, circular obstacles which can cross the boundary

  1. (a)

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

  2. (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

f i , j = L i - 1 2 + Ф μ , σ , L j - 1 2 + Ф ( μ , σ )
(54)

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 9-11 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 pseudo-randomly 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 9
figure 9

Example meshes for domains consisting of (a) periodic obstacles (b) pseudo-randomly placed obstacles.

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
figure 10

The graphs show (a) the first eigenvalue of the permeability tensor for an RVE with periodic unit cells and (b) the mean first eigenvalue of the permeability tensor for RVEs with random obstacles and constant Lagrange multipliers.

Figure 11
figure 11

Magnitude of velocity corresponding to (a) upper bound, (b) weak periodicity and (c) lower bound for a pressure gradient g - = 1 0 .

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 pseudo-random 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.

Figure 12
figure 12

Graphs showing (a) the mean first eigenvalue of the permeability tensor for RVEs with random obstacles and 4th order, piecewise applied Lagrange multipliers and (b) the standard deviation of K I .

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 Ω .

Figure 13
figure 13

Magnitude of velocity corresponding to a pressure gradient g - = 1 0 on RVEs with weak periodicity of order 0 for sizes (a) 4×4 (b) 8×8 (c) 20×20.

Figure 14
figure 14

Magnitude of velocity corresponding to a pressure gradient p M =[1 0] on RVEs with weak periodicity of order 0 for sizes (a) 4×4 (b) 8×8 (c) 20×20.

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 saddlepoint-nature 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, Ф(v )= μ 2 v sym : v 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

П ^ v , p M = sup p S ϵ P S β ϵ β inf γ ϵ G П v , p M , p S , β , γ = Ф ( v ) - b v , p M
(55)

From [23], we have the max-min inequality as

sup p M Θ p M = sup p M inf v П ^ v , p M inf v sup p M П ^ v , p M
(56)

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 Θ(pM) is concave in pM for all v. For future use, we note that a variation in v yields, by stationarity

Ф ( v ; δ v ) - b δ v , p M = 0
(57)

and a subsequent perturbation in v

Ф ′′ ( v ; δ v , d v ) - b δ v , d p M = 0
(58)

To show concavity of Θ, we make a perturbation in pM and choose δ v=dv in Equation 53, which yields

Θ p M ; d p M = Ф ( v ; d v ) - b d v , p M - b v , d p M = - b v , d p M
(59)

By yet another perturbation and the use of Equation 54

Θ ′′ p M ; d p M , d p M = - b d v , d p M = - Ф ′′ v ; d v , d v < 0
(60)

if Θ′′ is positive definite and b satisfies the required, classic inf-sup condition

inf q sup w b ( w , q ) w q γ b > 0
(61)

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.

    MATH  Google Scholar 

  2. Du X, Ostoja-Starzewski 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

    Article  MATH  MathSciNet  Google Scholar 

  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

    Article  Google Scholar 

  4. Sandström C, Larsson F, Runesson K, Johansson H: A two-scale 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

    Article  Google Scholar 

  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.

    MathSciNet  Google Scholar 

  6. Sánchez-Palencia E: Non-homogeneous Media and Vibration Theory. Springer, Berlin Heidelberg; 1980.

    MATH  Google Scholar 

  7. Allaire G: Homogenization of the Stokes flow in a connected porous medium. Asymptotic Anal 1989, 2: 203–222.

    MATH  MathSciNet  Google Scholar 

  8. Hornung U: Homogenization and porous media. Springer, Berlin Heidelberg; 1997.

    Book  MATH  Google Scholar 

  9. Nemat-Nasser S, Lori M, Datta SK (1996) Micromechanics: overall properties of heterogeneous materials. J Appl Mech 63(2): 561. Nemat-Nasser 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: Multi-scale computational homogenization: trends and challenges. J Comput Appl Math 2010,234(7):2175–2182. 10.1016/j.cam.2009.08.077

    Article  MATH  Google Scholar 

  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/s00466-006-0108-z

    Article  MathSciNet  Google Scholar 

  12. Löhnert S, Wriggers P: Homogenisation of microheterogeneous materials considering interfacial delamination at finite strains. Technische Mechanik 2003,23(2-4):167–177.

    Google Scholar 

  13. Ngo N, Tamma K: Microscale permeability predictions of porous fibrous media. Int J Heat Mass Transf 2001, 44: 3135–3145. 10.1016/S0017-9310(00)00335-5

    Article  MATH  Google Scholar 

  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

    Article  Google Scholar 

  15. Arbogast T, Lehr H: Homogenization of a Darcy-Stokes system modeling vuggy porous media. Comput Geosci 2006, 78712: 1–18.

    MathSciNet  Google Scholar 

  16. Ohman M, Runesson K, Larsson F: Computational Mesoscale modeling and homogenization of liquid-phase sintering of particle agglomerates. Technische Mechanik 2011, 32: 463–483.

    Google Scholar 

  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

    Article  Google Scholar 

  18. Ostoja-Starzewski M (2007) Microstructural randomness and scaling in mechanics of materials. Chapman & Hall/CRC, Boca Raton, FL. Ostoja-Starzewski 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

    Article  MATH  MathSciNet  Google Scholar 

  20. Larsson F, Runesson K, Saroukhani S, Vafadari R: Computational homogenization based on a weak format of micro-periodicity for RVE-problems. Comput Methods Appl Mech Eng 2011,200(1-4):11–26. 10.1016/j.cma.2010.06.023

    Article  MathSciNet  Google Scholar 

  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.

    Google Scholar 

  22. Patzák B, Bittnar Z: Design of object oriented finite element code. Adv Eng Softw 2001, 32: 759–767. 10.1016/S0965-9978(01)00027-8

    Article  MATH  Google Scholar 

  23. Boyd S, Vandenberghe L (2010) Convex Optimization. vol. 25. Cambridge University Press. pp 487-487. Chap. 1,10,11.

Download references

Acknowledgements

The project is financed by the Swedish Research Council (www.vr.se) under contract 2011-5388. The authors would also like to thank Håkan Johansson for fruitful discussions.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Carl Sandstöm.

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/s40323-015-0024-x.

Authors’ original submitted files for images

Below are the links to the authors’ original submitted files for images.

Authors’ original file for figure 1

Authors’ original file for figure 2

Authors’ original file for figure 3

Authors’ original file for figure 4

Authors’ original file for figure 5

Authors’ original file for figure 6

Authors’ original file for figure 7

Authors’ original file for figure 8

Authors’ original file for figure 9

Authors’ original file for figure 10

Authors’ original file for figure 11

Authors’ original file for figure 12

Authors’ original file for figure 13

Authors’ original file for figure 14

Authors’ original file for figure 15

Authors’ original file for figure 16

Authors’ original file for figure 17

Authors’ original file for figure 18

Authors’ original file for figure 19

Authors’ original file for figure 20

Authors’ original file for figure 21

Authors’ original file for figure 22

Authors’ original file for figure 23

Authors’ original file for figure 24

Authors’ original file for figure 25

Authors’ original file for figure 26

Authors’ original file for figure 27

Authors’ original file for figure 28

Authors’ original file for figure 29

Authors’ original file for figure 30

Authors’ original file for figure 31

Authors’ original file for figure 32

Authors’ original file for figure 33

Authors’ original file for figure 34

Authors’ original file for figure 35

Authors’ original file for figure 36

Authors’ original file for figure 37

Authors’ original file for figure 38

Authors’ original file for figure 39

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.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

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/s40323-014-0012-6

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s40323-014-0012-6

Keywords