Homogenization of Stokes System using Bloch Waves

In this work, we study the Bloch wave homogenization for the Stokes system with periodic viscosity coefficient. In particular, we obtain the spectral interpretation of the homogenized tensor. The presence of the incompressibility constraint in the model raises new issues linking the homogenized tensor and the Bloch spectral data. The main difficulty is a lack of smoothness for the bottom of the Bloch spectrum, a phenomenon which is not present in the case of the elasticity system. This issue is solved in the present work, completing the homogenization process of the Stokes system via the Bloch wave method.


Introduction and Main Result
We consider the Stokes system in which the viscosity is a periodically varying function of the space variable with small period ǫ > 0. Many physical phenomena (boiling flows, porous media, oil reservoirs, etc.) lead to mixture of fluids with different viscosities. For incompressible slow or creeping flows, such a situation is modeled by the system (1.1) for a Stokesian fluid with variable viscosity which is further assumed to be a periodic function. From the point of view of application, it is difficult to realize such a periodic distribution of droplets of one fluid in another without deforming the periodic structure, and (1.1) may seem as too much of an idealized system. Therefore, we also treat another model, which is a variant of the Stokes system and is physically more relevant. Namely, we consider the so-called incompressible elasticity system (1.11) which corresponds to a mixture of incompressible elastic phases in a composite material (this situation is quite common for rubber or elastomers).
We introduce now our first model. Assuming that the viscosity is a periodic function, the goal is to capture the effective viscosity of the mixture. To write down the model we start with a 1-periodic function µ = µ(y) ∈ L ∞ (T d ) or equivalently, a Y -periodic function where Y =]0, 1[ d . Here T d is the unit torus in R d . We assume µ(y) ≥ µ 0 > 0 a.e in T d .
(1. 6) We are interested here in the homogenization limit of (1.1), that is the asymptotic limit of the solution (u ǫ , p ǫ ) as ǫ → 0. This problem is very classical and its solution by means of a combination of two-scale asymptotic expansions and the method of oscillating test functions was provided in various references, including [2], [9], [13]. We recall their main results and follow the notations of [2] (cf. chapter I, section 10). The homogenized tensor (A * ) kl αβ , which represents "effective viscosity", is defined by µ(y)∇(χ k α + y α e k ) : ∇(χ l β + y β e l ) dy, (1.7) in which figure the cell test functions {χ k α ; α, k = 1 . . . d} solutions of the following problem in the torus T d : We impose T d χ k α dy = T d π k α dy = 0 to obtain uniqueness of the solutions. It is easy to see that the above homogenized tensor possesses the following "simple" symmetry, for any indices 1 ≤ α, β, k, l ≤ d, (A * ) kl αβ = (A * ) lk βα , (1.9) which corresponds to the fact that the fourth-order tensor A * is a symmetric linear map from the set of all matrices (or second-order tensors) into itself. Since we follow the notations of [2], the simple symmetry (1.9) seems a bit awkward since it mixes Latin and Greek indices but it is just the usual symmetry for a pair of indices (k, α) and (l, β) in a fourth-order tensor. In other words, (1.9) holds for a simultaneous permutation of k, l and α, β.
More precisely, we have the convergence of solutions: Note that the simple symmetry (1.9) does not imply that A * is symmetric in k, l or in α, β. However, in the homogenized equation (1.10), since A * is constant, only its symmetric version, obtained by symmetrizing in both k, l and α, β, plays a role. Let us next consider the second model of incompressible elasticity : Here the strain rate tensor is given by As before, there exists a unique solution (u ǫ s , p ǫ s ) of the above problem (1.11) in (H 1 0 (Ω)) d × L 2 0 (Ω) and using Korn's inequality and the inf-sup inequality, the following uniform bound can be proved : where the constant C does not depend on ǫ. Here the homogenized tensor (A * s ) kl αβ is given by where the cell test functions χ k α are now solutions in the torus T d of We impose T d χ k α = T d π k α = 0. It is known [3] that the above homogenized tensor possesses the following "full" symmetry, for any indices 1 ≤ α, β, k, l ≤ d, which corresponds to the fact that the fourth-order tensor A * s is a symmetric linear map from the set of all symmetric matrices into itself (the conditions (1.15) are the usual symmetry conditions for Hooke's laws in linearized elasticity). The homogenization limit of the problem (1.11) is again of the form (1.10) with A * s replacing A * . The first goal of this paper is to give an alternate proof of Theorem 1.1 using the Bloch Wave Method instead of two-scale asymptotic expansions and the method of oscillating test functions. The notion of Bloch waves is well-known in physics and mathematics [2], [4], [12], [16]. Bloch waves are eigenfunctions of a family of "shifted "spectral problems in the unit cell Y for the corresponding differential operator. Its link with homogenization theory was first explored in [2], [6], [10], [14]. The key point is that the homogenized operator can be defined in terms of differential properties of the bottom of the Bloch spectrum. The second goal of this paper is to explore this issue which is especially delicate in the case of Stokes equations. Indeed, it was discovered in [1] that the Bloch spectrum for the Stokes equations is not regular enough at the origin because of the incompressibility constraint. Therefore, its differential properties are all the more intricate to establish. Here we complete the task started in [1] and in particular we prove a conjecture of [1] on the homogenization of the Stokes system (1.1). Since the treatment of the incompressible elasticity system (1.11) is almost analogous to that of (1.1), we focus on (1.1) and we content ourselves in highlighting the main differences for (1.11) throughout the sequel.
The Bloch wave method for scalar equations and systems without differential constraints (like the incompressibility condition) was studied in [5,6,7,15]. In such cases, this approach gives a spectral representation of the homogenized tensor A * = (A * ) kl αβ in terms of the lowest energy Bloch waves and their behaviour for small momenta (what we call the bottom of the spectrum). For instance, the homogenized matrix in the scalar case was found to be equal to one -half of the Hessian of the ground energy (or first eigenvalue) at zero momentum. For a system, several bottom eigenvalues play a role and they are merely directionally differentiable by lack of simplicity. In the present case of the Stokes system, the situation is more complicated. The main characteristic of the Stokes system is the presence of the differential constraint expressing incompressibility of the fluid. One of its effects is that the Bloch energy levels are degenerate and the corresponding eigenfunctions are discontinuous at zero momentum. Even though energy levels are continuous at zero momentum, the second order derivatives are not (cf. Theorem 3.1). Thus, we cannot really make sense of the eigenvalue Hessian at zero momentum. Further, it is not clear if the homogenized tensor can be fully recovered from the Bloch spectral data. In fact, this issue is left open in [1]. In the non-self adjoint case treated in [15], only the symmetric part of the homogenized matrix is determined by Bloch spectral data and this is enough to determine the homogenized operator uniquely. Combining all these difficulties, the homogenization of Stokes system using Bloch waves is an interesting issue which is not a direct extension of previous results. Our work, roughly speaking, shows that Bloch spectral data does not determine the homogenized tensor uniquely, but determines the homogenized operator uniquely. This is in sharp contrast with the linear elasticity system treated in [7] in which the homogenized tensor was uniquely determined from Bloch spectral data. We see thus the effect of differential constraints (the incompressibility condition in the case of Stokes equations) on the homogenization process via Bloch wave method. For further discussion on this point, see Section 4. Bloch wave method of homogenization presented in Section 5 consists of localizing (1.1), taking its Bloch transform and passing to the limit to get the localized version of homogenized system in the Fourier space. Passage to the limit in the Bloch method is straight forward, though arguments are long. We do not run into the classical difficulty of having a product of two weakly convergent sequences. In fact, we use the Taylor approximation of Bloch spectral elements which gives strongly convergent sequences. This is one of the known features of the method. The required homogenized system is obtained by making a passage to the physical space from the Fourier space. Extracting macro constitutive relation and macro balance equation from the localized homogenized equation in the Fourier space turns out to be not very straight forward because of differential constraints.
The plan of this paper is as follows. In section 2, we recall from [1] the properties of Bloch waves associated with the Stokes operator. It turns out that the Bloch waves and their energies can be chosen to be directionally regular, upon modifying the spectral cell problem at zero momentum. Bloch transform using eigenfunctions lying at the bottom of the spectrum is also introduced in this section. Its asymptotic behaviour for low momenta is also described. Next, Section 3 is devoted to the computation of directional derivatives of Bloch spectral data. Even though these results are essentially borrowed from [1], some new ones are also included because of their need in the sequel. In particular we derive the so-called propagation relation linking the homogenized tensor A * with Bloch spectral data, and the extent to which it determines homogenized tensor is studied in Section 4. Using this information, we prove Theorem 1.1 in Section 5 following the Bloch wave homogenization method.

Bloch waves
In this section, we introduce Bloch waves associated to the Stokes operator following the lead of [1]. The Bloch waves are defined by considering the shifted (or translated) eigenvalue problem in the torus T d parametrized by elements in the dual torus which we take as T d again. We denote by y the points of the original torus and by η the points of the dual torus. The spectral Bloch problem amounts to find The solutions of (2.1) are a priori complex valued, so all functional spaces are complex valued too. Here, we denote D(η) = ∇ + iη the shifted gradient operator, with i the imaginary root √ −1. Its action on a vector function φ yields a matrix: (D(η)φ) kl = ∂φ l ∂y k + iη k φ l for all k, l = 1, . . . , d. The corresponding divergence operation yields a scalar: D(η) · φ = ∂φ k ∂y k + iη k φ k . Analogously, if φ is a matrix function then its shifted divergence D(η) · φ is a vector function obtained by acting D(η) on the column vectors of φ.
The main feature of (2.1) is that the state space keeps varying with η due to the differential constraints defined by the incompressibility of the fluid. That is why, the standard spectral theory for elliptic operators does not apply as such; it has to be modified. This is accomplished in [11]. Secondly, it is easily seen that when η = 0, the corresponding eigenvalue λ(0) is equal to zero and its multiplicity is d. In fact, we can take e k , k = 1 . . . d as eigenvectors (with corresponding eigen-pressure being zero). Because of this degeneracy, spectral elements of (2.1) are not guaranteed to be smooth at η = 0. Lack of regularity of the Bloch spectrum at η = 0 is an issue because the representation of the homogenized tensor in terms of Bloch spectral elements is then not clear. To overcome this difficulty, the idea is to consider directional regularity as we approach η = 0 [7]. Accommodating the directional limit at η = 0 requires a modification of the above shifted problem with the addition of a new constraint and corresponding Lagrange multiplier in the equation [1]. Fixing a direction e ∈ R d , |e| = 1 and taking η = δe, with δ > 0, we consider the modified problem: find Note that if δ = 0 then the relation e· It is natural to consider the system (2.2) with δ small as a perturbation from the following one which corresponds to δ = 0. We fix a unit vectorη ∈ S d−1 and we consider the eigenvalue problem: Existence of eigenvalues and eigenvectors for either (2.2) or (2.3) is proved in [1]. Let us recall their result, by specializing to the eigenvalue ν(η) = 0 of (2.3). Note that ν(η) = 0 is clearly an eigenvalue of multiplicity (d − 1) of (2.3) with corresponding eigenfunctions being constants, namely q 0 Doing perturbation analysis of the above situation, the following result was proved in [1].
The above theorem says that there are (d − 1) smooth curves emanating out of the zero eigenvalue as δ varies in an interval (−δ 0 , δ 0 ). We call them Rellich branches. Using them, for m = 1, . . . , (d−1), we can define the corresponding m th Bloch transform of For later purposes we need the Bloch transform for (H −1 (R d )) d elements also. Let us con- (2.5) Definition (2.5) is independent of the representation used for F ∈ (H −1 (R d )) d in terms of {g j , j = 0, ..., d} and is consistent with the previous definition (2.4) whenever F ∈ (L 2 (R d )) d .
In fact, by considering ∇ψ ) d then we can take g 0 = 0 and g j = ψe j j = 1, . . . , d, in (2.5) to obtain (2.6). That is, Bloch transform of gradient field is zero. Therefore the kernel of the Bloch transform B ǫ m,η : Roughly speaking since Bloch waves satisfy incompressibility condition the Bloch transform on gradient field vanish. Thus we may anticipate that the pressure effects may not be captured in the Bloch method. This impression is not correct. Indeed, as shown Section 5, by means of localization via a cut-off function, we manage to keep the pressure term.
Our next result is concerned with the asymptotic behavior of these Bloch transforms as ǫ → 0. Since φ m,η (y; 0) is a fixed unit vector (= φ 0 m,η ) orthogonal toη and independent of y (see Theorem 2.1), we have whereĝ denotes the Fourier transform of g and we recall thatη = ξ |ξ| .
By using Cauchy-Schwarz, the second term on the above right hand side can be estimated by the quantity Recall that δ is a function of (ǫ, ξ), namely δ = ǫ|ξ|. This quantity is easily seen to converge to zero as ǫ → 0 for each fixed ξ because of the directional continuity of φ m,η (., δ) → φ 0 m,η in (L 2 (T d )) d as δ → 0. We merely use the continuity of the m th Rellich branch at δ = 0 with values in (L 2 (T d )) d . On the other hand, thanks to our normalization, the integral on K is bounded by a constant independent of (ǫ, ξ). The proof is completed by a simple application of the Dominated Convergence Theorem which guarantees that the second term on the above right hand side converges strongly to 0 in L 2 loc (R d ξ ) as ǫ → 0. Since compactly supported elements are dense in (L 2 (R d )) d , we have the following : COROLLARY 2.1. In the setting of Theorem 2.2, if g ǫ be a sequence in (L 2 (R d )) d such that its support is contained in a fixed compact set K ⊂ R d , independent of ǫ and g ǫ → g in L 2 (R d ) d then we have the following strong convergence We recall the classical orthogonal decomposition : Let us denote By our choice, {φ 0 1,η , . . . , φ 0 d−1,η ,η} forms an orthonormal basis in R d , and so we can deduce the following : Proof. The proof is immediate, as {φ 0 1,η , . . . , φ 0 d−1,η } forms an orthogonal basis in R d−1 and φ 0 m,η · g = 0 for all m = 1, . . . , d − 1, so g(ξ) = c(ξ)ξ for some scalar c ∈ L 2 (R d ). Now if c = 0, it contradicts the hypothesis g ∈ X ⊥ . Thus c = 0. Consequently, g = 0.
COROLLARY 2.2. In the setting of Theorem 2.2 and Proposition 2.1 if g ǫ be a sequence in Proof. The proof simply follows as X ⊥ is a closed subspace of L 2 (R d ) d , so the limit g ∈ X ⊥ and the result follows by applying Proposition 2.1.
Bloch waves being incompressible are transversal. Longitudinal direction is missing and it has to be added to get the full basis. Naturally, asymptotics of the Bloch transform contains information of the Fourier transform only in transversal directions. It contains no information in the longitudinal direction. Because of this features, in the homogenization limit also, there is no information in the longitudinal direction. This is however proved to be enough to complete the homogenization process because the limiting velocity field is incompressible. See Section 5.
Using the above information in (3.2), we find that (φ ′ m,η (y; 0), q ′ m,η (y; 0)) is a solution of the following cell problem : Comparing this with (1.8), it can be seen that that φ ′ m,η (y; 0) is given by (see [1]) : where ζ m,η ∈ C d is a constant vector (independent of y), orthogonal toη. In other words, the y-dependence of φ ′ m,η (y; 0) is completely determined by the cell test function χ r α (y), solution of problem (1.8).
In a similar manner, the derivative of the eigenpressure q m,η (y; 0) is given by (see [1]): That is, the y-dependence of q ′ m,η (y; 0) is completely determined by the cell test function π r α (y) , solution of problem (1.8).
Second order derivatives : Next we differentiate (3.2) with respect to δ to obtain : (3.8) We consider (3.7) at δ = 0 and by integrating over T d , we get where M(η, A * ) is the symmetric matrix whose entries are given by This is nothing but a contraction of the homogenized tensor A * . As a simple consequence of (3.9), we get It is also follows that M(η, A * )φ 0 m,η ⊥ φ 0 m ′ ,η for all m = m ′ .
By summarizing the above computations, we have where ζ m,η ∈ C d is a constant vector (independent of y), orthogonal toη.
(iv) The derivative of the eigenfunction q m,η (y; δ) at δ = 0 satisfies: (v) The second derivative of the eigenvalue λ m,η (δ) and q 0,m,η (δ) at δ = 0 satisfy the relation where M(η, A * ) is the symmetric matrix whose entries are given by M(η, A * ) kl = (A * ) kl αβη αηβ . REMARK 3.1. The above matrix M(η, A * ) is precisely that which must be positive definite in the Legendre-Hadamard definition of ellipticity. A relation analogous to (3.10) is called "propagation relation" in [7] in the study of linearized elasticity system and it shows how the homogenized tensor A * enters into the Bloch wave analysis. The above relation (3.10) generalizes the relation (22) in [1].
In the linearized elasticity system, the propagation relation is an eigenvalue relation. Here, relation (3.10) can again be seen as an eigenvalue problem, posed in the (d − 1)-dimensional subspace orthogonal toη. More precisely, 1/2λ ′′ m,η (0) is an eigenvalue and φ 0 m,η (which is orthogonal toη) is an eigenvector of the restriction of the matrix M(η, A * ) to the subspaceη ⊥ . In (3.10) 1/2q ′′ 0,m,η (0) is the Lagrange multiplier corresponding to the constraint that the eigenvalue problem is posed in the (d−1)-dimensional subspace orthogonal toη.

Case of Symmetrized gradient :
We recall the incompressible elasticity system (1.11) with the symmetrized gradient introduced in Section 1.
We introduce Bloch waves associated to the Stokes operator defined in (3.11).

Recovery of homogenized tensor from Bloch waves
In the scalar self-adjoint case, it is known that the homogenized matrix is equal to one-half the Hessian of the first Bloch eigenvalue at zero momentum [6]. In the general (non-symmetric) scalar case, treated in [15], it was shown that only the symmetric part of the homogenized matrix is determined by the Bloch spectrum and it is given again by the same one-half of the Hessian of the first Bloch eigenvalue (which exists by virtue of the Krein-Rutman theorem). The fact that only the symmetric part of the homogenized matrix plays a role is not a big surprise since, the homogenized tensor A * being constant, the differential operator depends only on the symmetric part of A * .
In the case of systems, another phenomenon takes place. For example, the linearized elasticity system (in which there are no differential constraints) was treated in [7] where it was recognized that not only Bloch eigenvalues but also Bloch eigenfunctions at zero momentum are needed to determine the homogenized tensor. More precisely, this connection between Bloch eigenvalues and eigenfunctions, on the one hand, and the homogenized tensor, on the other hand, was expressed via a relation called propagation relation in [7] which uniquely determines the homogenized tensor.
In the case of Stokes system, a new phenomenon arises because of the presence of a differential constraint (the incompressibility condition). Even though there is an analogue of the propagation relation (see (3.10) above), it does not determine uniquely the homogenized tensor. In fact the propagation relation (3.10) is unaltered if we add a multiple of I ⊗I (where I is the d × d identity matrix) to the homogenized tensor. The homogenized Stokes operator clearly remains the same under such an addition since it corresponds to adding a gradient of the velocity divergence which vanishes because of the incompressibility constraint. The authors in [1] conjectured that the homogenized Stokes tensor is uniquely characterized by the propagation relation up to the addition of a term c(I ⊗ I) (where c is a constant). We prove this assertion in the case of the Stokes system (1.11) with a symmetrized gradient. For the other Stokes system (1.1), the homogenized tensor is not uniquely determined by the propagation relation (3.10). In this section, we investigate this non-uniqueness. Neverheless, we shall prove that for both Stokes systems the homogenized operators (1.10), and its equivalent for the symmetric gradient case of (1.11), are uniquely determined.
Our concern now is the following question: to what extent do the Bloch spectral elements determine the homogenized tensor A * via the propagation relation (3.10) ? Since λ ′′ m,η (0), q ′′ 0,m,η (0), φ 0 m,η are known from Bloch spectral data, it follows that M(η, A * )φ 0 m,η is uniquely determined via the relation (3.10). But it may happen that different tensors A * give rise to the same matrix M(η, A * ). Three main results are proved in this section and they are stated in the following three propositions.  Conversely, let us assume that there are two fourth-order tensors A * and B * , possessing the simple symmetry (1.9) and such that M(η, A * )φ 0 m,η = M(η, B * )φ 0 m,η , m = 1, ..., d − 1, for allη ∈ S d−1 . We must then deduce (4.1). For convenience, the proof is divided into five steps.
Step 1. First of all, we check that the matrix M(η, A * ) is symmetric. By interchanging the dummy indices α and β and using the simple symmetry (1.9) of the homogenized coefficients, (A * ) jl αβ = (A * ) lj βα , we get which shows the required symmetry.
(iii) (i = j, k = j). In this case, Now together with i = l we have N jl jj + N jl ij + N jl ji = cδ lj . (4.11) Then both j = l or j = l cases lead to verify (4.5) and (4.6) respectively.
(iv) Similarly, for (k = l and i = l) N ji ii + N ji ik + N ji ki = c(δ ji + δ jk ). (4.12) Together with k = j we have Then both i = j or i = j cases lead to verify (4.5) and (4.6) respectively.
Step 4. Now we consider the two remaining cases not covered in (4.6).
(i) (i, k) = (j, l). Then from (4.9) we have For i = k it gives using (4.8) (4.14) (ii) Similarly, for (k, i) = (j, l), together with i = k we have Step 5. Let us set N = N − c(I ⊗ I). Thanks to the properties (4.5) and (4.6), we can easily check that N is an anti-symmetric tensor in the sense that it satisfies N jl ik = −N jl ki = −N lj ik . whenever, (i, k) = (j, l) and (k, i) = (j, l) (4.16) From its very definition N also possesses the symmetry N jl ik = N lj ki . Thus N has all the properties listed in (4.2).
Next we extend Proposition 4.1 to the Stokes system (1.11), featuring a symmetric gradient tensor. In this case the propagation relation (3.10) is replaced by (3.14) and the homogenized tensor is denoted by A * s .  This symmetry combined with the anti-symmetry established in the previous step implies that N = 0. Note that antisymmetry property holds precisely for the interchange of those pairs of indices for which symmetry property does not hold. This can be seen as follows: whenever (i, k) = (j, l) and (k, i) = (j, l) Similarly, whenever (i, k) = (j, l) or (k, i) = (j, l) together with i = k; from (4.14), (4.15) we have N ik ik + N ik ki = c = N ki ik + N ki ki . Then using (4.18) and (4.8) we clearly have REMARK 4.1. The conclusion of the above proposition was conjectured in [1] and it is proved here to be true whenever we are working with the system (1.11) with symmetrized gradient. However, it is not true with the full gradient Stokes system (1.1) as shown by Proposition 4.1. However, in both of these cases the propagation relation fixes the homogenized operator (1.10) uniquely, as is stated in the following proposition. Proof. We have to check that A * and B * define the same Stokes differential operator for divergence-free vector fields. Indeed the Fourier symbol of the operator is (A * − B * ) kl αβ ξ α ξ β which, by virtue of (4.4), is equal to cξ k ξ l which is precisely the symbol of the operator u → −c∇(∇ · u) which vanishes on the space of divergence free functions.

Homogenization result
This section is devoted to a proof of Theorem 1.1, our main homogenization result stated in the first section. It is based on the tools that we have introduced so far. A similar proof is given for the linear elasticity problem in [15]. However, the presence of a pressure and a differential constraint in the Stokes system seriously complexifies the analysis and has a non-trivial effect in the homogenization process. Besides, we also bring some simplifications to the proof given in [15].
There are several steps in the proof. First, we localize the Stokes system (1.1) by applying a cut-off function technique to the velocity u in order to get the equation (5.2) in the whole R d . Next, by taking the Bloch transformation B ǫ m,η (1 ≤ m ≤ d − 1) of the equation (5.2) and passing to the limit, we arrive at the homogenized equation in the Fourier space. Finally, we take the inverse Fourier transform to go back to the physical space which gives our desired result.
Notation: in the sequel L.H.S. stands for left hand side, and R.H.S. for right hand side.
Step 1. Localization of the velocity u : Let v ∈ D(Ω) be arbitrary. Then vu ǫ and p ǫ satisfy (for l = 1, . . . , d) where, Note that, g ǫ l and h ǫ l correspond to terms containing zero and first order derivatives of µ ǫ respectively. In the sequel, we extend u ǫ and p ǫ by zero outside Ω and such extensions are denoted by the same letters.
They satisfy the following system because of (3.1) : Below, we treat each term of (5.5) one by one.