S Brisard

L Dormieux
email: luc.dormieux@enpc.fr

K Sab

A variational form of the equivalent inclusion method for numerical homogenization

Keywords: Homogenization, Elasticity, Equivalent Inclusion Method, Galerkin Discretization, Green Operator

Due to its relatively low computational cost, the equivalent inclusion method is an attractive alternative to traditional full-field computations of heterogeneous materials formed of simple inhomogeneities (spherical, ellipsoidal) embedded in a homogeneous matrix. The method can be seen as the discretization of the Lippmann-Schwinger equation with piecewise polynomials. Contrary to the original approach of Moschovidis and Mura, who discretized the strong form of the Lippmann-Schwinger equation through Taylor expansions, we propose in the present paper a Galerkin discretization of the weak form of this equation. Combined with the new, mixed boundary conditions recently introduced by the authors, the resulting method is particularly well-suited to homogenization. It is shown that this new, variational approach has a number of benefits: (i) the resulting linear system is well-posed, (ii) the numerical solution converges to the exact solution as the maximum degree of the polynomials tends to infinity and (iii) the method can provide rigorous bounds on the apparent properties of the statistical volume element, provided that the matrix is stiffer (or softer) than all inhomogeneities. This paper presents the formulation and implementation of the new, variational form of the equivalent inclusion method. Its efficiency is investigated through numerical applications in 2D and 3D elasticity.

Introduction

A wide range of tools are available for the determination of the macroscopic properties of heterogeneous materials. Standard micromechanical techniques, based on the solution of [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problems[END_REF], are found at one end of this range. Examples of such techniques are the scheme [START_REF] Mori | Average stress in matrix and average elastic energy of materials with misfitting inclusions[END_REF] [see also [START_REF] Benveniste | A new approach to the application of Mori-Tanaka's theory in composite materials[END_REF]] and the self-consistent scheme [START_REF] Walpole | On the overall elastic moduli of composite materials[END_REF][START_REF] Kröner | Bounds for effective elastic moduli of disordered materials[END_REF]. Although not based on the solution of [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problems[END_REF], the generalized self-consistent scheme [START_REF] Christensen | Solutions for effective shear properties in three phase sphere and cylinder models[END_REF][START_REF] Hervé | n-layered inclusion-based micromechanical modelling[END_REF]) also falls into this category. The main asset of these schemes is their flexibility, as they apply to linear as well as non-linear behaviours [START_REF] Suquet | Continuum Micromechanics[END_REF], while leading to analytical or semi-analytical models; various types of couplings can also be included. The weakness of such models lies in the limited amount of microstructural information they can account for in a quantitative way: volume fractions and possibly distributions of orientations, but no higher-order correlations.

Full-field calculations can be found at the other end of this range. These computation return the 'exact' (up to some numeri-cal error) local strains and stresses of a specific realization of the material. Standard (finite elements, boundary elements) as well as non-standard [based for example on the Fast Fourier Transform [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF]] homogenization techniques can be used for this computation. Contrary to the micromechanical technique, the amount of microstructural information that full-field calculations can account for is limited by the resolution (or fineness of the mesh) only, thus leading to models with greater accuracy. However, such type of computations are also much more demanding both in terms of memory and processing time, to the effect that advanced techniques such as parallelization or multi-threading are usually invoked. This can be problematic for stastical analyses requiring computations on numerous microstructures [START_REF] Kanit | Determination of the size of the representative volume element for random composites: statistical and numerical approach[END_REF][START_REF] Ostoja-Starzewski | Material spatial randomness: From statistical to representative volume element[END_REF].

The equivalent inclusion method is an intermediate technique which is both more accurate than micromechanical approaches, and less costly than full-field computations. Like the latter, it requires a realization of the material under investigation, and produces estimates of the local strain and stress fields. Because the total number of degrees of freedom remains limited, the required memory and processing time are lowered. It should also be noted that the preparation of the computation is straightforward, as no mesh of the microstructure is required. The price to pay is, of course, a loss of accuracy with respect to full-field computations. Still, the accuracy gain over micromechanical approaches can be significant.

The equivalent inclusion method finds its roots in the pioneering work of [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problems[END_REF], who derived the exact solution to the problem of a single ellipsoidal inhomogeneity, embedded in an infinite matrix. To this end, he introduced a strictly equivalent ellipsoidal inclusion, with appropriate eigenstrain. [START_REF] Moschovidis | Two-ellipsoidal inhomogeneities by the equivalent inclusion method[END_REF] then extended this work to find an approximate solution to the problem of multiple inhomogeneities, embedded in an infinite matrix. Similarly, each inhomogeneity is susbtituted with an inclusion, which is equivalent only in an approximate sense. The eigenstrains to which the inclusions are submitted are found by solving a linear system, resulting from a discretization of the so-called consistency equation by means of Taylor expansions. Since this consistency equation is equivalent to the Lippmann-Schwinger equation [START_REF] Korringa | Theory of elastic constants of heterogeneous media[END_REF][START_REF] Zeller | Elastic constants of polycrystals[END_REF][START_REF] Kröner | On the physics and mathematics of self-stresses[END_REF], the equivalent inclusion method can be seen as a particular discretization of this integral equation.

Although many applications of the equivalent inclusion method found in the literature deal with pairs of inhomogeneities [START_REF] Moschovidis | Two-ellipsoidal inhomogeneities by the equivalent inclusion method[END_REF][START_REF] Rodin | On the problem of linear elasticity for an infinite region containing a finite number of non-intersecting spherical inhomogeneities[END_REF][START_REF] Shodja | Interacting cracks and ellipsoidal inhomogeneities by the equivalent inclusion method[END_REF], some authors have successfully applied this technique to problems involving many inclusions [START_REF] Fond | Effects of mechanical interactions on the hydrostatic stress in randomly distributed rubber particles in an amorphous polymer matrix[END_REF][START_REF] Benedikt | On elastic interactions between spherical inclusions by the equivalent inclusion method[END_REF]. It has also recently been used to enrich the shape functions of finite element models with non-conforming meshes [START_REF] Novák | A micromechanics-enhanced finite element formulation for modelling heterogeneous materials[END_REF].

The equivalent inclusion method is limited to composites with ellipsoidal inhomogeneities and homogeneous matrix. While this might be too stringent a restriction for many real materials, this method is a valuable tool for the quantification -on model microstructures-of the influence of some specific microstructural parameters, such as particle-size distribution, or local orientational order.

In this paper, a variational form of the equivalent inclusion method is introduced. It is based on a Galerkin discretization (through piecewise polynomials) of a modified Lippmann-Schwinger equation, which is better-suited to numerical homogenization than the standard Lippmann-Schwinger equation (Brisard et al., 2013b). The new, variational form of the equivalent inclusion method improves upon its original, Taylor-based form in many respects. Indeed, the resulting linear system is well-posed, and convergence with respect to the degree of the polynomial approximations can be proved. By contrast, none of these results has been established for the Taylor-based equivalent inclusion method, and examples can be found where increasing the degree of the polynomials actually lowers the accuracy of the approximations [START_REF] Fond | Mechanical interaction between spherical inhomogeneities: an assessment of a method based on the equivalent inclusion[END_REF]. Furthermore, the principle of Hashin and Shtrikman (1962a), extended to the modified Lippmann-Schwinger equation (Brisard et al., 2013b), can be used to show that the Galerkin-based equivalent inclusion method provides rigorous bounds on the macroscopic properties.

The present paper is organized as follows. Sec. 2 provides some background information regarding the original (Taylorbased) equivalent inclusion method and its relation to the Lippmann-Schwinger equation, as well as its application to numerical ho-mogenization. Sec. 3 introduces the variational framework for the proposed alternative version of the equivalent inclusion method. The discretized equations are derived from a Galerkin approach. Using the principle of Hashin and Shtrikman (1962a) [see also Hill (1963b) and [START_REF] Willis | Bounds and self-consistent estimates for the overall properties of anisotropic composites[END_REF]], it is then shown that the newly introduced method can provide bounds on the macroscopic properties. Sec. 4 provides technical details on the implementation of the method. In particular, the calculation of the self-influence and influence pseudotensors is presented. Finally, some applications are proposed in Sec. 5 to illustrate the efficiency of the method.

In the remainder of this paper, we refer to the original (Taylor-based) form of the equivalent inclusion method as EIM-T. Likewise, EIM-G refers to the new (Galerkin-based) form of the equivalent inclusion method.

Background

The present section provides background information on the equivalent inclusion method and its application to numerical homogenization. The presentation below differs from the original paper by [START_REF] Moschovidis | Two-ellipsoidal inhomogeneities by the equivalent inclusion method[END_REF], in that the Green operator for strains Γ ∞ 0 is used in place of the potentials Φ i j and Ψ i j . Furthermore, derivation of the EIM-T is traditionally based on eigenstrains [START_REF] Moschovidis | Two-ellipsoidal inhomogeneities by the equivalent inclusion method[END_REF]. By contrast, the EIM-G presented in this paper is most conveniently derived using prestresses (more precisely, polarizations). As a consequence, slightly altering the terminology introduced by [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problems[END_REF], an inclusion is defined in the present paper as a region in a homogeneous medium (no material mismatch), subjected to eigenstrains or prestresses. By contrast, an inhomogeneity is a region in a heterogeneous medium which is occupied by a homogeneous material whose mechanical properties differ from those of the matrix.

Sec. 2 is organized as follows. In Sec. 2.1, the tuple notation is first introduced, which will allow for more compact formulas in the remainder of this paper. Then, in Sec. 2.2, the fundamentals of the equivalent inclusion (EIM-T) method are briefly recalled. Following [START_REF] Moschovidis | Two-ellipsoidal inhomogeneities by the equivalent inclusion method[END_REF], the derivation of the method is based on eigenstrains; unlike these authors, our formulation of the consistency equation makes explicit use of the Green operator for strains. This reveals important connections between the EIM-T and the well-known Lippmann-Schwinger equation [START_REF] Korringa | Theory of elastic constants of heterogeneous media[END_REF][START_REF] Zeller | Elastic constants of polycrystals[END_REF][START_REF] Kröner | On the physics and mathematics of self-stresses[END_REF]. Finally, applicability of this method to numerical homogenization is discussed in Sec. 2.3.

On the tuple notation

In the remainder of this paper, quantities indexed with a dot (e.g. k

• , l • , x • , y • , . . .) denote d-tuples; in particular, 0 • denotes the null tuple 0 • = (0, . . . , 0) d times .
Depending on the context, the components of the tuple can be naturals or reals. For example, k

• = (k i ) 1≤i≤d ∈ N d denotes a multi-index, while x • = (x i ) 1≤i≤d ∈ R d denotes
the coordinates of the vector x with respect to the global basis (e i ) 1≤i≤d . Similarly, x ′ • = (x ′ j ) 1≤ j≤d ∈ R d denotes the coordinates of x with respect to an auxiliary basis (e ′ j ) 1≤ j≤d . For any tuple a • , a + (resp. a × ) denotes the sum (resp. product) of its components

a + = a 1 + • • • + a d , and a × = a 1 • • • a d .
In particular, dx × is the following volume element

dx × = dx 1 • • • dx d .
It is understood that all binary operators (+, ×, ≤,. . . ) should apply component-wise to tuples. For example,

a • + b • = c • , where c i = a i + b i for all i = 1, . . . , d, similarly a • ≤ b • if, and only if, a i ≤ b i for all i = 1, . . . , d.
Exponentiation of tuples is then defined as follows

x n • • = x n 1 1 • • • x n d d ,
and it should be noted that the above quantity is not intrinsic, as it depends on the components x 1 , . . . , x d of the vector x in a specified basis. In particular, x n • would not be a meaningful notation, since

x ′n • • x n • • . From the binomial theorem (a • + b • ) n • = 0 • ≤l • ≤n • n • l • a l • • b n • -l • • , (1)
where the following multi-binomial coefficient has been defined

n • l • = n 1 k 1 • • • n d k d = n 1 ! k 1 ! (n 1 -k 1 )! • • • n d ! k d ! (n d -k d )!
.

The equivalent inclusion method

The N-inhomogeneity problem. The equivalent inclusion method allows the computation of the approximate strain field in an assembly of inhomogeneities, embedded in an infinite matrix, and submitted to a uniform strain at infinity. The microstructure under consideration is formed of N inhomogeneities Ω 1 , . . . , Ω N ⊂ R d (d: dimension of the physical space), embedded in an infinite matrix

Ω 0 = R d \ (Ω 1 ∪ • • • ∪ Ω N ).
Inhomogeneity α is centered at x α (α = 1, . . . , N), and χ α denotes the indicator function of inhomogeneity α, translated back to the origin. In other words, x ∈ Ω α if, and only if, χ α (xx α ) = 1; besides, χ α = χ β if inhomogeneities α and β are identical (up to a translation).

Assuming a linear elastic behaviour for all constituants, C α (resp. C 0 ) denotes the elastic stiffness of inhomogeneity α (resp. the matrix). The local stiffness, defined over the whole space R d then reads

C(x) = C 0 + N α=1 χ α (x -x α ) (C α -C 0 ) . (2)
With this notation at hand, the problem to be solved by the equivalent inclusion method reads (Fig. 1, left)

∇ x • σ = 0 (x ∈ R d ), (3a) σ(x) = C(x) : ε(x) (x ∈ R d ), (3b) ε(x) = ∇ s x u (x ∈ R d ), (3c) u(x) ∼ E ∞ • x ( x → +∞), (3d)
under the constraint that ε -E ∞ be square-integrable to ensure well-posedness (Brisard et al., 2013b). In Eqs.

(3), u (resp. ε, σ) denotes the local displacement (resp. strain, stress); furthermore, ∇ s u denotes the symmetric gradient of the displacement. Boundary conditions (3d) will be called KUBC ∞ [Kinematic Uniform Boundary Conditions at infinity, see Brisard et al. (2013b)].

The N-inclusion problem. Prior to solving problem (3), [START_REF] Moschovidis | Two-ellipsoidal inhomogeneities by the equivalent inclusion method[END_REF] consider the following auxiliary N-inclusion problem, in which Ω 1 , . . . , Ω N are inclusions rather than inhomogeneities, submitted to the uniform strain at infinity E ∞ and the eigenstrain η (Fig. 1, right)

∇ x • σ = 0 (x ∈ R d ), (4a) σ(x) = C 0 : ε(x) -η(x) (x ∈ R d ), (4b) ε
(x) = ∇ s x u (x ∈ R d ), (4c) u(x) ∼ E ∞ • x ( x → +∞). (4d)
It should be noted that η(x) = 0 for all x ∈ Ω 0 ; besides, contrary to the single-inclusion problem of [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problems[END_REF], the eigenstrain η is allowed to vary spatially within each inclusion. The solution to the N-inclusion problem (4) depends linearly on E ∞ and η. It is conveniently written in terms of the Green operator for strains Γ ∞ 0 of the homogeneous, infinite medium C 0 . It is recalled that this operator returns the opposite of the strain induced in the reference, unbounded medium by any prestress ̟ [START_REF] Korringa | Theory of elastic constants of heterogeneous media[END_REF][START_REF] Zeller | Elastic constants of polycrystals[END_REF][START_REF] Kröner | On the physics and mathematics of self-stresses[END_REF]. In other words

Γ ∞ 0 * ̟ = -∇ s u, (5)
where ' * ' stands for the standard convolution product, and u is the solution to the following partial differential equation in

R d ∇ • C 0 : ∇ s u + ̟ = 0, (6)
with ∇ s u square integrable and u(x) → 0 as x → +∞. From Eqs. ( 5) and ( 6), it is readily found that the local strain ε solution to the N-inclusion problem (4) is given by the convolution product

ε = E ∞ + Γ ∞ 0 * C 0 : η . ( 7
)
From inhomogeneities to inclusions. To solve the N-inhomogeneity problem (3), [START_REF] Moschovidis | Two-ellipsoidal inhomogeneities by the equivalent inclusion method[END_REF] stress induced in the heterogeneous medium by E ∞ and the material mismatch, they derive the so-called consistency equation, with the eigenstrain η as unknown. As a second step, they then propose a procedure to discretize the consistency equation, thus allowing the numerical computation of approximate solutions to the N-inhomogeneity problem (3). These two steps are briefly summarized below. An equivalent formulation of the N-inhomogeneity problem (3) is first introduced

σ = C 0 : ε u(x) ∼ E ∞ • x, x → +∞ Ω α Ω 1 Ω 2 Ω 3 Ω 4 Ω 5 σ = C α : ε σ = C 0 : ε u(x) ∼ E ∞ • x, x → +∞ Ω α Ω 1 Ω 2 Ω 3 Ω 4 Ω 5 σ = C 0 : ε -η
∇ x • σ = 0 (x ∈ R d ), (8a) σ
(x) = C 0 : ε(x) -η(x) (x ∈ R d ), (8b) η
(x) = -C -1 0 : (C(x) -C 0 ) : ε(x) (x ∈ R d ), (8c) ε
(x) = ∇ s x u (x ∈ R d ), (8d) u(x) ∼ E ∞ • x ( x → +∞), (8e)
where the quantity η [defined by Eq. (8c)] can be viewed as an eigenstrain applied to Ω 1 , . . . , Ω N , which are no longer inhomogeneities, but inclusions. It should be noted that the eigenstrain η does not affect the matrix Ω 0 since by construction, η(x) = 0 for x ∈ Ω 0 [see Eq. ( 8c)]. The solution to Eqs. (8a), (8b), (8d), (8e) is given by Eq. ( 7), which, upon substitution in Eq. ( 8c), leads to the so-called consistency equation [START_REF] Moschovidis | Two-ellipsoidal inhomogeneities by the equivalent inclusion method[END_REF])

(C -C 0 ) -1 : C 0 : η = -E ∞ -Γ ∞ 0 * C 0 : η . (9)
The consistency equation is an integral equation with the local eigenstrain η as unknown. Its discretization is addressed in the next paragraph. It should be noted that if the local strain ε is chosen as main unknown, Eq. ( 9) reduces to

ε + Γ ∞ 0 * ((C -C 0 ) : ε) = E ∞ ,
which is known as the Lippmann-Schwinger equation [START_REF] Korringa | Theory of elastic constants of heterogeneous media[END_REF][START_REF] Zeller | Elastic constants of polycrystals[END_REF][START_REF] Kröner | On the physics and mathematics of self-stresses[END_REF].

Finally, it will be seen in Sec. 3 that the most natural choice of main unknown for the derivation of the EIM-G is neither the local eigenstrain η, nor the local strain ε, but the local polarization τ, defined as follows [see [START_REF] Willis | Bounds and self-consistent estimates for the overall properties of anisotropic composites[END_REF] among others]

τ = (C -C 0 ) : ε = -C 0 : η. ( 10
)
It is readily observed that τ is homogeneous to a prestress, and that τ(x) = 0 for all x ∈ Ω 0 . With τ as main unknown, Eq. ( 9) reads

(C -C 0 ) -1 : τ + Γ ∞ 0 * τ = E ∞ . ( 11
)
Discretization of the consistency equation. Following [START_REF] Moschovidis | Two-ellipsoidal inhomogeneities by the equivalent inclusion method[END_REF], an approximate solution to Eq. ( 9) is found by use of piecewise polynomial expansions of degree p for the eigenstrain η

η(x) = N α=1 k • ∈I p χ α (x -x α )(x • -x α,• ) k • η k • α , (12)
for all x ∈ R d , where I p is the following set of multi-indices k

• I p = k • ∈ N d , k + ≤ p . For α = 1, . . . , N and k • ∈ I p , η k • α is a set of d 2 constants η k • α,i j , such that η k • α,i j = η k • α, ji (i, j = 1, . . . , d); however, η k • α
should not be considered as a true second-rank, symmetric tensor. Indeed, Eq. ( 12) is not intrinsic (because exponentiation is not an intrinsic operation). In particular, η k • α does not follow the required transformation laws under a change of basis. Such objects will be called pseudotensors in the remainder of this paper.

Substituting the assumed form ( 12) for η into Eq. ( 7), it is readily seen that the local strain ε resulting from the imposed uniform strain at infinity E ∞ and the piecewise polynomial eigenstrain η defined by Eq. ( 12) reads

ε(x) = E ∞ + N β=1 l • ∈I p D l • β (x -x β ) : η l • β , (13)
where the fourth-rank pseudotensor D l • β is defined as follows

D l • β (x) = y∈R d χ β (y)y l • • Γ ∞ 0 (x -y) : C 0 dV y . (14)
In the above expression, dV y denotes the d-dimensional volume element at y ∈ R d . D l • β (x) depends on the shape of inclusion β, but not on its position in space; for ellipsoidal inclusions, these pseudotensors can be calculated explicitely [see [START_REF] Moschovidis | Two-ellipsoidal inhomogeneities by the equivalent inclusion method[END_REF] and references therein].

Substitution of Eq. ( 14) into Eq. ( 9) leads to the following set of equations which must hold for all α = 1, . . . , N and x ∈ Ω α

(C α -C 0 ) -1 : C 0 : η(x) = -E ∞ - N β=1 l • ∈I p D l • β (x -x β ) : η l • β . (15)
Clearly, these equations can only be solved for the η l • β in an approximate sense. Indeed, the left-hand side of Eq. ( 15) is polynomial over Ω α , while the right-hand side is not. Following [START_REF] Moschovidis | Two-ellipsoidal inhomogeneities by the equivalent inclusion method[END_REF], it is natural (since the left-hand side is polynomial) to approximate the right-hand side of Eq. ( 15) with a Taylor expansion in the neighborhood of the center

x α of Ω α D l • β (x -x β ) ≃ k • ∈I p x • -x α,• k • D k • l • αβ : η l • β , (16)
where

D k • l • αβ = 1 k 1 ! • • • k d ! ∂ k 1 • • • ∂ k d D l • β ∂x k 1 1 • • • ∂x k d d x α -x β ,
and, upon substitution into Eq. ( 15)

k • ∈I p (x • -x α,• ) k • (C α -C 0 ) -1 : C 0 : η k • α = -E ∞ - N β=1 k • ,l • ∈I p x • -x α,• k • D k • l • αβ : η l • β .
The above identity must be satisfied for all x ∈ Ω α . Equating all monomials results in the following set of algebraic equations with unknowns

η k • α (α = 1, . . . , N, k • ∈ I p ) (C α -C 0 ) -1 : C 0 : η k • α + N β=1 l • ∈I p D k • l • αβ : η l • β =        -E ∞ if k • = 0 • , 0 otherwise. ( 17
)
The system of equations ( 17) defines the original form of the equivalent inclusion method (EIM-T), first derived by [START_REF] Moschovidis | Two-ellipsoidal inhomogeneities by the equivalent inclusion method[END_REF]. From the solution to these equations, an approximate expression of the local strain ε(x), solution to the N-inhomogeneity problem (3) can be computed by means of Eq. ( 13).

Since the seminal paper of [START_REF] Moschovidis | Two-ellipsoidal inhomogeneities by the equivalent inclusion method[END_REF], alternative linear systems of equations have been proposed. Instead of the inclusions centers x α , [START_REF] Benedikt | On elastic interactions between spherical inclusions by the equivalent inclusion method[END_REF] compute Taylor expansions of D l • β [see Eq. ( 16)] at selected points of interest (where the stresses are seeked); they show that the resulting estimates are more accurate. However, this approach requires a new inversion of the linear system (17) if a new set of points of interest is considered. [START_REF] Shodja | Interacting cracks and ellipsoidal inhomogeneities by the equivalent inclusion method[END_REF] use point collocation in place of Taylor expansions, which makes the assembly of the linear system much simpler.

Closed-form expressions of the pseudotensors D k • l • αβ are available [START_REF] Mura | Micromechanics of defects in solids[END_REF]; alternatively, they can be computed numerically [START_REF] Nakasone | Numerical equivalent inclusion method: a new computational method for analyzing stress fields in and around inclusions of various shapes[END_REF]. Then the linear system of equations (17) can be assembled and solved numerically for the coefficients of the polynomial expansion of the eigenstrain, η k • α [see Eq. ( 12)].

The EIM-T presented above has two shortcomings. First, there is no guarantee that the linear system given by Eq. ( 17) is indeed invertible; the EIM-T therefore lacks robustness, as it might fail in some circumstances. Second, increasing the degree p of the expansions does not necessarily improve the quality of the approximate solution given by Eqs. ( 12) and ( 13). This has already been reported by many authors [START_REF] Rodin | On the problem of linear elasticity for an infinite region containing a finite number of non-intersecting spherical inhomogeneities[END_REF][START_REF] Fond | Mechanical interaction between spherical inhomogeneities: an assessment of a method based on the equivalent inclusion[END_REF][START_REF] Benedikt | On elastic interactions between spherical inclusions by the equivalent inclusion method[END_REF]. Also, convergence (as p → +∞) of this numerical method is not established.

In Sec. 3, an alternative to Taylor expansions is proposed to carry out the discretization of the consistency equation, and overcome these problems. Following a Galerkin approach based on the weak form of Eq. ( 11), a new system of equations is derived for the polynomial expansion of the unknown local polarization τ. This new system replaces the system obtained by Moschovidis and Mura [see Eq. (17)]; for this new system, wellposedness and convergence (as p → +∞) to the true polarization can then be investigated with the help of standard mathematical tools [see e.g. [START_REF] Ern | Theory and Practice of Finite Elements[END_REF]].

Before this variational approach is introduced in Sec. 3, application of the EIM to homogenization problems is first dicussed in Sec. 2.3. It is shown that a slight modification [introduced by Brisard et al. (2013b)] of the standard Lippmann-Schwinger equation ( 11) makes it much better-suited to this kind of problems. Galerkin discretization will therefore be carried out on the modified Lippmann-Schwinger equation [see Eq. ( 21)].

The EIM for numerical homogenization

Provided that the degree p of the polynomial expansions is not too high, each inhomogeneity has relatively few degrees of freedom; in other words, the linear system (17) remains small, even for large assemblies of inhomogeneities. The EIM is therefore an attractive tool for numerical homogenization, where hundreds to thousands of inhomogeneities must be considered simultaneously. Attention must however be paid to the boundary conditions, as was already noted by [START_REF] Fond | Effects of mechanical interactions on the hydrostatic stress in randomly distributed rubber particles in an amorphous polymer matrix[END_REF].

We seek to determine (numerical estimates of) the apparent stiffness of a statistical volume element Ω ⊂ R d [SVE, following the terminology of [START_REF] Ostoja-Starzewski | Material spatial randomness: From statistical to representative volume element[END_REF]]. As previously, the microstructure is formed of inhomogeneities Ω 1 , . . . , Ω N embedded in a homogeneous matrix. However, contrary to Sec. 2.2, the SVE Ω is now a bounded domain.

It is recalled that the apparent stiffness C app relates the macroscopic stress σ to the macroscopic strain ε σ = C app : ε,

where overlined quantities denote volume averages over the bounded domain

Ω ε = 1 |Ω| Ω ε and σ = 1 |Ω| Ω σ.
In Eq. ( 18), the macroscopic strain and stress are computed from the solution to the following auxiliary problem, which states that the SVE is in (elastic) equilibrium (Fig. 2, left)

∇ x • σ = 0 (x ∈ Ω), σ(x) = C(x) : ε(x) (x ∈ Ω), ε(x) = ∇ s x u (x ∈ Ω).
In addition, appropriate boundary conditions must be specified. The boundary conditions most frequently met are kinematic uniform (KUBC), static uniform (SUBC) and periodic (PBC)

σ = C 0 : ε u(x) = E • x Ω α Ω 1 Ω 2 Ω 3 Ω 4 Ω 5 σ = C α : ε Ω σ = C 0 : ε u(x) ∼ E ∞ • x, x → +∞ Ω α Ω 1 Ω 2 Ω 3 Ω 4 Ω 5 σ = C α : ε Ω Figure 2
: Left: evaluation of the apparent stiffness traditionally requires the determination of the elastic equilibrium of the SVE Ω, with appropriate boundary conditions imposed at the boundary ∂Ω (boundary conditions of type KUBC are represented here). Right: in order to evaluate the apparent stiffness by means of the EIM, the SVE Ω must be embedded in an infinite homogeneous medium with same stiffness as the matrix, with imposed strain E ∞ at infinity (boundary conditions of type KUBC ∞ ). In this case, the loading parameter E ∞ is not the macroscopic strain ε.

boundary conditions [see among others Hill (1963a[START_REF] Hill | The essential structure of constitutive laws for metal composites and polycrystals[END_REF], [START_REF] Mandel | Plasticité classique et viscoplasticité[END_REF] and [START_REF] Gusev | Representative volume element size for elastic composites: A numerical study[END_REF]]. With kinematic uniform and periodic boundary conditions, the loading parameter is the macroscopic strain (E = ε is specified), while with static uniform boundary conditions, the loading parameter is the macroscopic stress (Σ = σ is specified). For a finite-size SVE, these boundary conditions produce three different estimates of the apparent stiffness:

C app KUBC , C app SUBC , C app PBC .
However, under statistical homogeneity and ergodicity assumptions, all estimates converge to the effective stiffness C eff as the size of the SVE tends to infinity [START_REF] Sab | On the homogenization and the simulation of random materials[END_REF].

As previously mentioned, the EIM is an attractive technique to compute numerical estimates of the apparent stiffness of the SVE Ω. This would of course require to embed the bounded domain Ω in an infinite, homogeneous medium with same stiffness as the matrix (Fig. 2, right). However, the EIM-T is a discretization of the N-inhomogeneity problem (3) with boundary conditions of type KUBC ∞ , which is not strictly speaking a standard auxiliary problem, for two reasons. First, elastic equilibrium is specified in the whole space R d (as opposed to the bounded domain Ω ⊂ R d ); second, neither macroscopic strain nor macroscopic stress can be specified. Indeed, there is no simple relationship between the sole loading parameter E ∞ [see Eq. (3d)] and the macroscopic strain or stress (in particular, ε E ∞ ).

Notwithstanding, the N-inhomogeneity problem (3) can still be used to define an apparent stiffness C app KUBC ∞ of the SVE Ω. Indeed, from the linearity of this problem with respect to E ∞ , both σ and ε depend linearly on the strain at infinity

E ∞ ε = A KUBC ∞ : E ∞ , σ = B KUBC ∞ : E ∞ ,
where A KUBC ∞ and B KUBC ∞ are fourth-rank localization tensors. The apparent stiffness results from the elimination of the loading parameter E ∞ [see e.g. [START_REF] Fond | Mechanical interaction between spherical inhomogeneities: an assessment of a method based on the equivalent inclusion[END_REF][START_REF] Fond | Effects of mechanical interactions on the hydrostatic stress in randomly distributed rubber particles in an amorphous polymer matrix[END_REF]]

C app KUBC ∞ = B KUBC ∞ : A -1 KUBC ∞ , (19)
and this definition is consistent in the sense that C app KUBC ∞ → C eff as |Ω| → +∞ (under statistical homogeneity and ergodicity assumptions). To sum up, the macroscopic strain and stress must be computed for six independent values of the loading parameter E ∞ (strain at infinity); from these six computations, the localization tensors A KUBC ∞ and B KUBC ∞ can be computed, and the apparent stiffness C app KUBC ∞ can be retrieved. However, evaluating ε and σ from the EIM estimate of the prestress involves complex surface integrals at the boundary ∂Ω of the SVE [START_REF] Fond | Mechanical interaction between spherical inhomogeneities: an assessment of a method based on the equivalent inclusion[END_REF][START_REF] Fond | Effects of mechanical interactions on the hydrostatic stress in randomly distributed rubber particles in an amorphous polymer matrix[END_REF]. These operations are costly and potentially inaccurate. Furthermore, there is no guarantee that the apparent stiffness thus defined is symmetric, positive definite for finite-size SVEs. Using an energy approach to define the apparent stiffness partially resolves these issues [START_REF] Rodin | On the problem of linear elasticity for an infinite region containing a finite number of non-intersecting spherical inhomogeneities[END_REF].

Mixed boundary conditions (MBC) were introduced and analyzed in detail by Brisard et al. (2013b) as a way to circumvent the above-mentioned shortcomings of the boundary conditions of type KUBC ∞ . The resulting system (20) of partial differential equations is very similar to the initial problem (3) with boundary conditions of type KUBC ∞ . Yet, it is better-suited to homogenization problems, since the loading parameter E is now the macroscopic strain [in other words, the solution to problem (20) satisfies ε = E]. Besides, the macroscopic stress can readily be derived from the average polarization, the computation of which does not involve complex surface integrals. Furthermore, the corresponding apparent stiffness is symmetric, positive definite, and we have

C app SUBC ≤ C app MBC ≤ C app KUBC ,
where inequalities should be understood in the sense of quadratic forms.

The N-inhomogeneity problem with mixed boundary conditions is depicted in Fig. 3. The SVE Ω is embedded in an infinite, homogeneous medium with same stiffness as the matrix. It is then submitted to a uniform strain E at infinity, as well as a surface load t • n applied to the boundary ∂Ω of the SVE, where t is a constant, second-rank, symmetric tensor and n is the outward normal to ∂Ω

∇ x • σ = 0 (x ∈ R d ), (20a)
σ(x) = C(x) : ε(x) (x ∈ R d ), (20b)
ε(x) = ∇ s x u (x ∈ R d ), (20c) [[σ]](x) • n(x) = -t • n(x) (x ∈ ∂Ω), (20d) u(x) ∼ E • x ( x → +∞), (20e) ε = E, (20f)
where [[σ]] denotes the stress-jump across the boundary ∂Ω of the SVE Ω. From the above set of equations, it is apparent that t is not a loading parameter, but should be chosen so as to ensure that Eq. (20f) holds. In other words, the unique loading parameter in the N-inhomogeneity problem with boundary conditions of type MBC is the strain at infinity E, which is forced to coincide with the average strain over the SVE. At first sight, problem (20) might seem more complex than problem (3) because of the additional unknown t. However, it can be shown (Brisard et al., 2013b) that for ellipsoidal SVEs Ω,

σ = C 0 : ε u(x) ∼ E • x, x → +∞ Ω α Ω 1 Ω 2 Ω 3 Ω 4 Ω 5 σ = C α : ε Ω n [[σ]] • n = -t • n Figure 3:
The N-inhomogeneity problem with mixed boundary conditions (MBC). Similarly to the situation depicted in Fig. 2 (right), the SVE Ω is embedded in an infinite medium with same stiffness as the matrix. It is submitted to a uniform strain at infinity E, as well as a surface load t • n applied at the boundary ∂Ω of the SVE Ω. The constant, second-rank, symmetric tensor t is a free parameter which is chosen so as to ensure that the unique loading parameter E coincides with the macroscopic strain ε.

Eqs. ( 20) are equivalent to the following modified Lippmann-Schwinger equation, from which the supplementary unknown, constant tensor t is absent

(C -C 0 ) -1 : τ + Γ ∞ 0 * (τ -χτ) = E, (21)
where χ denotes the indicator function of Ω and τ again denotes the polarization [see Eq. ( 10)]. Furthermore, the local strain ε and stress σ are readily retrieved from the solution τ to Eq. ( 21) (Brisard et al., 2013b)

ε = E -Γ ∞ 0 * τ -χτ and σ = C 0 : ε + τ. (22)
The above formulation as an integral equation shows that the N-inhomogeneity problems with boundary conditions of type MBC [see Eqs. ( 21) and ( 22)] and KUBC ∞ [see Eqs. ( 11) and ( 7)] are very similar, and can be discretized through the same techniques [including Taylor expansions used by [START_REF] Moschovidis | Two-ellipsoidal inhomogeneities by the equivalent inclusion method[END_REF] for Eq. ( 9)]. From Eq. ( 22) 2 and the constraint ε = E, the following identity is readily found

τ = C app MBC -C 0 : E, (23)
where C app MBC denotes the apparent stiffness associated with mixed boundary conditions. As a consequence, unlike the original Lippmann-Schwinger equation ( 11), estimation of the apparent stiffness from a numerical solution to the modified Lippmann-Schwinger equation ( 21) does not involve complex surface integrals. Indeed, only the volume average of the polarization is required, which is trivial to compute provided a simple form is assumed for the numerical estimate of the polarization [see Eq. ( 26)]. Slightly modifying the equation to be solved numerically therefore allows to overcome one of the shortcomings of the EIM-T.

Adopting a Galerkin approach (rather than Taylor expansions) for the discretization of this equation allows to overcome all other shortcomings of the EIM-T listed in Sec. 2. This will be shown in Sec. 3, where a variational form of the EIM is proposed, based on the modified Lippmann-Schwinger equation (21), and piecewise polynomial approximations of the local polarization [see Eq. ( 26)].

To close this section, it should be noted that other types of mixed boundary conditions have been proposed in the past by various authors [see among others [START_REF] Hazanov | Order relationships for boundary conditions effect in heterogeneous bodies smaller than the representative volume[END_REF]; [START_REF] Coenen | Novel boundary conditions for strain localization analyses in microstructural volume elements[END_REF]; [START_REF] Salmi | Apparent and effective mechanical properties of linear matrix-inclusion random composites: Improved bounds for the effective behavior[END_REF]]. However, such boundary conditions usually apply at the boundary ∂Ω of the (bounded) SVE. This precludes the use of the Green operator for strains Γ ∞ 0 of the whole space to solve the underlying system of partial differential equations. By contrast, the mixed boundary conditions introduced in Brisard et al. (2013b) and used in the present paper apply partly at infinity. Eq. ( 21) shows that this approach allows to associate a very simple integral equation of the Lippmann-Schwinger type to the underlying (seemingly complex) system of partial differential equations ( 20). It should be emphasized that with these boundary conditions, the Green operator for strains Γ ∞ 0 of the whole space arises naturally, without any approximation. This is essential, since rigorous bounds on the apparent stiffness can then be produced through this approach, as shown in Sec. 3.2. Finally, the fact that the Green operator of the whole space (rather than the bounded SVE Ω) appears in the integral equation ( 21) is what makes the EIM practical. Indeed, Γ ∞ 0 has very simple expressions; besides, it is translation invariant, which makes assembly of the underlying linear system (37) [or indeed ( 17)] less costly.

Derivation of the proposed method

It is first noted that the equivalence between the N-inhomogeneity problem (20) with boundary conditions of type MBC and the modified Lippmann-Schwinger equation ( 21) holds for ellipsoidal SVEs only (Brisard et al., 2013b). Therefore, it will be assumed in the remainder of this paper that Ω is an ellipsoid.

Galerkin discretization of the consistency equation

The point of departure is the following weak form of the modified Lippmann-Schwinger equation ( 21), which is found by premultiplying with a test function ̟ ∈ V, and averaging over the SVE Ω

Find τ ∈ V such that a(τ, ̟) = f (̟) for all ̟ ∈ V, ( 24
)
where V is the set of square integrable, second-rank, symmetric tensors, supported in R d \ Ω 0 ; a and f are bilinear and linear forms, respectively, given by

a(τ, ̟) = ̟ : (C -C 0 ) -1 : τ + ̟ : Γ ∞ 0 * τ -χτ , f (̟) = E : ̟.
It should be noted that since the SVE Ω is an ellipsoid, the theorem of [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problems[END_REF] applies, and for all x ∈ Ω Γ ∞ 0 * χτ (x) = P Ω : τ, where P Ω denotes the Hill tensor of the domain Ω with respect to the reference medium C 0 . Therefore, the bilinear form a reduces to

a(τ, ̟) = ̟ : (C -C 0 ) -1 : τ + ̟ : Γ ∞ 0 * τ -̟ : P Ω : τ.
The above continuous variational problem is then discretized following a standard Galerkin procedure. More precisely, Eq. ( 24) is replaced with

Find τ p ∈ V p such that a(τ p , ̟ p ) = f (̟ p ) for all ̟ p ∈ V p , ( 25
)
where V p is the finite dimension set of piecewise polynomial polarizations τ p defined as follows

τ p (x) = N α=1 k • ∈I p χ α (x -x α )(x • -x α,• ) k • τ k • α . (26)
Comparison with Eqs. ( 10) and ( 12) shows that the present variational form of the EIM uses the same piecewise polynomial approximation of the polarization as the classical EIM of [START_REF] Moschovidis | Two-ellipsoidal inhomogeneities by the equivalent inclusion method[END_REF].

The discrete variational problem (25) reduces to a linear system, with the coefficients τ k • α of the polynomial expansion (26) as unknowns. Introducing the following expansion of the test function

̟ p ∈ V p ̟ p (x) = N α=1 k • ∈I p χ α (x -x α )(x • -x α,• ) k • ̟ k • α , (27)
the linear system to be solved is derived from the substitution of Eqs. ( 26) and ( 27) into Eq. ( 25). Recalling that all inclusions are completely included in Ω, so that

Ω χ α = R d χ α (α = 1, . . . N),
it is readily found that

̟ p = 1 |Ω| N α=1 k • ∈I p M k • α ̟ k • α , (28)
where

M k • α denotes the moment of order k • of inclusion α M k • α = x∈R d χ α (x)x k • • dV x . (29)
Consequently, the value of the linear form f at

̟ p ∈ V p reads f (̟ p ) = 1 |Ω| N α=1 k • ∈I p M k • α ̟ k • α : E. ( 30
)
Similarly

̟ p : (C -C 0 ) -1 : τ p = 1 |Ω| N α=1 k • ,l • ∈I p M k • +l • α ̟ k • α : (C α -C 0 ) -1 : τ l • α , (31)
and

̟ p : Γ ∞ 0 * τ p = 1 |Ω| N α,β=1 k • ,l • ∈I p ̟ k • α : T k • l • αβ (r αβ ) : τ l • β , (32)
where r αβ = x βx α and T k • l • αβ (r) denotes the following influence pseudotensor

T k • l • αβ (r) = x,y∈R d χ α (x)χ β (y)x k • • y l • • Γ ∞ 0 (r + y -x)dV x dV y . ( 33
)
For k • = l • = 0 • , these pseudotensors coincide with T II and T I J introduced by [START_REF] Berveiller | The problem of two plastic and heterogeneous inclusions in an anisotropic medium[END_REF] and [START_REF] Molinari | The problem of elastic inclusions at finite concentration[END_REF]. It is assumed that the inclusions do not overlap. Therefore, if α β, then r + yx is always non-zero in the above integral, and the integrand is never singular. By contrast, if α = β, then r = 0, and the integrand is singular. This suggests to single out so-called self-influence pseudotensors, defined as follows

S k • l • α = T k • l • αα (0), ( 34
)
and it is noted that

|Ω α | -1 S 0 • 0 • α
is the Hill tensor of inclusion α with respect to the reference medium C 0 . Eq. ( 32) then reads

̟ p : Γ ∞ 0 * τ p = 1 |Ω| N α,β=1 α β k • ,l • ∈I p ̟ k • α : T k • l • αβ (r αβ ) : τ l • β + 1 |Ω| N α=1 k • ,l • ∈I p ̟ k • α : S k • l • α : τ l • α , (35)
Evaluation of the bilinear form a for τ p , ̟ p ∈ V p is then carried out from Eqs. ( 28), ( 31) and ( 35)

a(τ p , ̟ p ) = 1 |Ω| N α=1 k • ,l • ∈I p M k • +l • α ̟ k • α : (C α -C 0 ) -1 : τ l • α + 1 |Ω| N α=1 k • ,l • ∈I p ̟ k • α : S k • l • α : τ l • α + 1 |Ω| N α,β=1 α β k • ,l • ∈I p ̟ k • α : T k • l • αβ (r αβ ) : τ l • β - 1 |Ω| 2 N α,β=1 k • ,l • ∈I p M k • α M l • β ̟ k • α : P Ω : τ l • β . ( 36
)
Gathering Eqs. ( 30) and ( 36), it is found that the discretized variational problem ( 25) is equivalent to the following linear system with unknowns τ k • α (one equation for each value of α = 1, . . . , N and k

• ∈ I p ) l • ∈I p       M k • +l • α (C α -C 0 ) -1 + S k • l • α - M k • α M l • α |Ω| P Ω       : τ l • α + N β=1 β α l • ∈I p         T k • l • αβ (r αβ ) - M k • α M l • β |Ω| P Ω         : τ l • β = M k • α E. ( 37
)
For inhomogeneities with simple geometries, closed-form expressions of the moments

M k • α , self-influence pseudotensors S k • l •
α and influence pseudotensors T k • l • αβ (r) can be derived (see Sec. 4); assembly of the above system is therefore straightforward. Following an approach similar to that of [START_REF] Brisard | Combining Galerkin approximation techniques with the principle of Hashin and Shtrikman to derive a new FFT-based numerical method for the homogenization of composites[END_REF] (see in particular Sec. 4.2 in this reference), it can be shown that this system is well-posed.

To sum up, the EIM-G is defined through the discretization (26) and the linear system (37). The solution to this system provides the coefficients of the expansion of the approximate polarization, τ p . In turn, estimates of the strain and stress in the inclusions can be computed as follows

ε p (x) = (C α -C 0 ) -1 : τ p (x) (x ∈ Ω α ), (38)
and

σ p (x) = C 0 : ε p (x) + τ p (x) (x ∈ Ω α ).
For points belonging to the matrix, Eq. ( 38) is meaningless (since C -C 0 is singular), and the approximate strains must be computed from the convolution product ε p = E -Γ ∞ 0 * τ p ; this involves complex numerical integration, which are not necessary for the determination of the approximate apparent stiffness, C app,p MBC . Indeed, from Eq. ( 23)

C app,p MBC : E = C 0 : E + τ p = C 0 : E + 1 |Ω| N α=1 k • ∈I p M k • α τ k • α . (39)
To close this section, it should be noted that for p = 0 (constant polarization in each inclusion), the equations of the EIM-G are very close to those found by previous authors [START_REF] Rodin | On the problem of linear elasticity for an infinite region containing a finite number of non-intersecting spherical inhomogeneities[END_REF][START_REF] Rodin | The overall elastic response of materials containing spherical inhomogeneities[END_REF][START_REF] Molinari | The problem of elastic inclusions at finite concentration[END_REF][START_REF] El Mouden | Thermoelastic properties of composites containing ellipsoidal inhomogeneities[END_REF]. In the present work however, periodization of the unit-cell is not necessary owing to the introduction of the mixed boundary conditions and the corresponding correction term involving the Hill tensor P Ω of the domain Ω [see Eq. ( 36)].

Bounds on the apparent stiffness

Provided that the matrix is stiffer (resp. softer) than all inhomogeneities, the above approach leads to an upper bound (resp. lower bound) on the apparent stiffness of the SVE Ω. Indeed, an extremum principle of the Hashin and Shtrikman (1962a) type has been proved for the modified Lippmann-Schwinger equation (21) (Brisard et al., 2013b). More precisely, introducing for ̟ ∈ V the following functional

H(̟) = f (̟) -1 2 a(̟, ̟),
it can readily be shown that H is stationary at the solution τ to Eq. ( 21), and

H(τ) = 1 2 E : C app MBC -C 0 : E. ( 40
)
It was further shown in (Brisard et al., 2013b) that H is minimum (resp. maximum) at τ, provided C α ≤ C 0 (resp. C α ≥ C 0 ) for all α = 1, . . . , N.

In particular, if one of these conditions is fulfilled, H(τ) can be compared to H(τ p ), where τ p is the EIM-G estimate of τ. Since τ p is the solution to the discretized variational problem (25), it verifies a(τ p , τ p ) = f (τ p ), and from Eq. ( 39), it is readily found that

H(τ p ) = 1 2 E : C app,p MBC -C 0 : E. (41)
If the matrix is stiffer than all inhomogeneities, then H is minimum at τ, and H(τ) ≤ H(τ p ). Gathering Eqs. ( 40) and ( 41)

C app MBC ≤ C app,p MBC ,
the EIM-G estimate of the apparent stiffness is in fact an upper bound. Conversely, if the matrix is softer than all inhomogeneities, then the EIM-G estimate of the apparent stiffness is a lower bound.

To sum up, if the matrix is stiffer (resp. softer) than the inhomogeneities, then the variational form of the EIM returns upper (resp. lower) bounds on the apparent stiffness. This result further guarantees that the quality of the solution will not deteriorate as the order p of the method is increased. Indeed, the discretized variational problem (25) can be seen as the optimization of H over the finite-dimension subspace V p of V. As p grows, so does the subspace V p , and the quality of the approximate optimum is therefore improved. This desirable property is not necessarily observed with the EIM-T [START_REF] Fond | Mechanical interaction between spherical inhomogeneities: an assessment of a method based on the equivalent inclusion[END_REF], which is not based on a variational setting.

Implementation of the method

The key point for the implementation of the method is the computation of the self-influence and influence pseudotensors. The complex integrals involved [see Eqs. ( 33) and ( 34)] must be evaluated with the help of a computer algebra system1 . Examples of such computations can be found in Brisard et al. (2013a).

Computation of the self-influence pseudotensors

The self-influence pseudotensors S k • l • α are defined as follows [see Eqs. ( 33) and (34)]

S k • l • α = x,y∈R d χ α (x)χ α (y)x k • • y l • • Γ ∞ 0 (y -x)dV x dV y .
In the above expression, principal values are taken according to Eq. (A.1) in order to remove the singularity at x = y

S k • l • α = x∈R d χ α (x)x k • • lim δ→0 y∈R d y-x >δ χ α (y)y l • • Q 0 (y -x)dV y dV x + x∈R d χ α (x)x k • +l • • P 0 dV x . (42)
The last term reduces to M k • +l • α P 0 . For convex inclusions, the following change of variables is performed to compute the first term (see Fig. 4 andAppendix B.1 for more details) where R α (x, n) is defined as the largest value of r such that χ α (x + rn) = 1 for χ α (x) = 1. It can then be shown that

y = x + rn, with n = 1 and 0 ≤ r ≤ R α (x, n), (43)
S k • l • α = M k • +l • α P 0 + x,y∈R d χ α (x)χ α (y)x k • • y l • • -x l • • Q 0 (y -x)dV y dV x + x∈R d n =1 χ α (x)x k • +l • • log R α (x, n)Q 0 (n)dS n dV x , (44)
where dS n denotes the surface element at n on the unit sphere.

For spherical inclusions, the last term vanishes, while the second term can be computed analytically (see Appendix B.2).

Computation of the influence pseudotensors

The influence pseudotensors T k • l • αβ (r) are given by Eq. ( 33) (see also Fig. 5). As inclusions α and β do not overlap, only the regular part Q 0 of Γ ∞ 0 is involved in this computation [see Eq. (A.1)].

Direct evaluation of integral (33) is complex (if possible at all). In particular, the method proposed by [START_REF] Berveiller | The problem of two plastic and heterogeneous inclusions in an anisotropic medium[END_REF] -based on Fourier transforms-for the evaluation of the lowest-order influence pseudotensor T k • l • αβ (r) with k • = l • = 0 • does not generalize to higher orders. In the present paragraph, an alternative approach is therefore proposed, based on multipole expansions.

Let r = rn ( n = 1) be a fixed vector, and r ′ = r (n + ξ). From the homogeneity of degree -d of the regular part Q 0 of Γ ∞ 0 , it is found that

Q 0 (r ′ ) = r -d Q 0 (n + ξ),
wich is expanded in a Taylor series with respect to the powers of the components of ξ. The resulting expansion is then inserted in Eq. ( 33), with ξ = (yx) /r. Each term of the series being a monomial in the variables x i /r and y j /r (i, j = 1, . . . , d), its integral for x ∈ Ω α and y ∈ Ω β can readily be computed (provided Ω α and Ω β have simple shapes). This procedure returns the influence pseudotensors T k • l • αβ (r) as series of the negative powers of r. In practice, this series must be truncated, and the above method ought to be viewed as approximate. However, experience shows that for spherical inclusions, all computed expansions have a finite number of nonzero terms. This suggests that the present computation would in fact be exact in that case, although this conjecture remains unproved for the time being.

The above procedure involves rather tedious, but very systematic, algebra. As such, it lends itself to straightforward implementation in a computer algebra system; the resulting expressions are then imported into a numerical code, in charge of assembling and inverting the global matrix.

While it is possible to compute literal expressions of the influence tensors in any basis, the resulting expressions are intractable, and lead to cluttered code. αβ is not a tensor (see Appendix C).

Applications

In the present section, two applications of the variational form of the equivalent inclusion method (EIM-G) derived in Secs. 3 and 4 are proposed. The first application (see Sec. 5.1) is a plane strain elasticity application, while the second application (see Sec. 5.2) is a 3D elasticity application.

It is recalled that the EIM-G requires the SVE to be of ellipsoidal shape (Brisard et al., 2013b). Therefore, circular (resp. spherical) SVEs are considered in Sec. 5.1 (resp. 5.2).

Monodisperse assemblies of circular pores in plane strain elasticity

The present example deals with porous media in plane strain elasticity. The circular SVE Ω contains N circular pores of radius a. The porosity is φ = 0.4; with N = 160 pores, the radius R of the SVE Ω is R = 20a (see Fig. 6). The shear modulus and Poisson ratio of the matrix are µ 0 (arbitrary value) and ν 0 = 0.3.

To account for statistical fluctuations of the apparent mechanical properties of each individual SVE, 1 000 configurations were considered. The (mean) apparent shear modulus of these Table 1: Upper bounds on the apparent shear modulus of a monodisperse assembly of circular pores in plane strain elasticity, for increasing orders p of the EIM-G. For each value of p, the corresponding number of degrees of freedom (DOFs) is also reported.

microstructures is then estimated through the EIM-G, and the results are reported in Table 1 for various values of p (it is recalled that p is the maximum degree of the polynomial expansions used to approximate the polarization τ). Due to the large number of independent configurations considered here, the amplitude of the 99 % confidence interval is smaller than one unit in the last place of each value reported in Table 1.

Since the reference medium (the matrix) is stiffer than the inclusions, results presented in Sec. 3.2 apply, and the estimates of µ app are in fact upper bounds on this quantity. This is consistent with the fact that this bound decreases as p increases, as expected (the functional H is minimized on sub-spaces of increasing dimension).

Observation of the results presented in Table 1 shows that increasing p significantly improves the upper bound on µ app . Indeed, from p = 0 to p = 3, the upper-bound is reduced by approximately 20 %, while the total number of degrees of freedom is multiplied by a factor 10.

In order to quantify the error on the apparent shear modulus, the above results were compared to finite element (FEM) estimates computed on the same 1 000 configurations. Strictly speaking, EIM-G and FEM computations are not equivalent. Indeed, kinematic uniform boundary conditions were adopted for convenience for the FEM models, while EIM-G models require mixed boundary conditions (see Sec. 2.3 for a definition of these two types of boundary conditions). Since the SVEs under consideration are very large (R/a = 20), finite-size effects should be negligible, and the apparent shear moduli resulting from these two sets of boundary conditions are expected to coincide (Hill, 1963a).

The FEM estimate of the apparent shear modulus was found to be 0.244 µ 0 ; again, the 99 % confidence interval is narrow enough to ensure that this value is correct up to one unit in the last place. Using this value as a reference, Fig. 7 shows in semi-log scale the relative error on µ app as a function of the number of degrees of freedom. Observing Fig. 7, it seems that the EIM-G converges exponentially with respect to the number of degrees of freedom. This empirical result was expected, as similar asymptotic behaviors are also observed with the pversion of the FEM, which is very similar in spirit to EIM-G [START_REF] Szabó | The p-version of the finite element method[END_REF].

Comparison of the respective sizes of the EIM-G and FEM models emphasizes the efficiency of the EIM-G. Indeed, each of the 1 000 FEM model contained about 320 000 degrees of freedom (this figure was variable from one configuration to another), while with only 4 800 degrees of freedom, the EIM-G achieves a relative error of approximately 1.2 %. In other words, the EIM-G can provide at a relatively low cost estimates of some quantities of interest with a small (but finite) error. Of course, if high accuracy is required, then the FEM should be preferred to the EIM-G.

To close this section, it should also be noted that the Hashin and Shtrikman (1962b) upper-bound on µ app reads in this case: µ HS+ = 0.349 µ 0 . Clearly, this bound is a poor estimate in plane strain elasticity of the effective shear modulus of a composite with circular inclusions. Even the 0-th order EIM-G bound leads to an improvement of 11 %.

Polydisperse assembly of spherical pores in 3D elasticity

The present example deals with porous media in 3D elasticity. The spherical SVE Ω contains N 1 = 20 (resp. N 2 = 40, N 3 = 140) spherical pores of radius a 1 (resp. a 2 = 0.7 a 1 , a 3 = 0.4 a 1 ). The total porosity is φ = 0.45, so that the radius R of the SVE Ω is R = 4.56 a 1 (see Fig. 6, right). The shear modulus and Poisson ratio of the matrix are µ 0 (arbitrary value) and ν 0 = 0.3.

To account for statistical fluctuations, 100 such SVEs were generated. It should be noted that fluctuations were smaller in the present, 3D case than in the previous, plane strain one. In both cases, the amplitude of the statistical error (99 % confidence interval) was identical, even if the number of generated configurations was ten times smaller in the 3D application. The (mean) apparent shear modulus µ app of these microstructures is then estimated through the EIM-G, and the results are reported in Table 2 for various values of p. Again, the results are accurate (with probability 99 %) within one unit in the last place. For the same reasons as previously, they can be considered as upper bounds on µ app (which is again consistent with the fact that this bound decreases as p increases).

However, observing Table 2, it seems that the relative improvement of the bounds is much lower in the present case than in the previous, plane strain case, for the same increase of the number of degrees of freedom. It is necessary to thoroughly investigate this apparent loss of efficiency; this investigation will be carried out in future work.

Conclusion

In the present paper, we proposed a new form of the equivalent inclusion method (EIM). Recognizing that the EIM can be seen as an approximation of the Lippmann-Schwinger equation, we proposed a Galerkin discretization of the weak form of this integral equation as an alternative to the discretization of its strong form, based on Taylor expansions, initially proposed by [START_REF] Moschovidis | Two-ellipsoidal inhomogeneities by the equivalent inclusion method[END_REF]. All details are provided for the implementation of this new method.

Combined with the new, mixed boundary conditions recently introduced by the authors (Brisard et al., 2013b), the resulting numerical method is an attractive tool for homogenization. In particular, provided that the matrix is stiffer or softer than all inhomogeneities, it can provide rigorous bounds on the macroscopic properties, at a relatively low cost. Numerical examples show that the method is extremely efficient in plane strain elasticity, less so in three dimensional elasticity. We intend to investigate in future work the origins of this loss of efficiency, and how to overcome it. We will also extend the method to ellipsoidal inhomogeneities, which should allow us to capture local orientational order effects.

Gathering Eqs. (B.7), (B.9) and (B.10), the following expression of S k • l • α is finally found for spherical inclusions

S k • l • α = M k • +l • P 0 + 0 • ≤h • ≤l • h • 0 • 0≤2 j≤h + (-1) h + h + l • h • h + 2 j         x∈R d n =1 χ α (x) (x • n) h + -2 j a 2 α -x 2 + (x • n) 2 j x k • +l • -h • • n h • • Q 0 (n)dS n dV x         .
The above sum reduces to integrals of trigonometric polynomials. Its computation therefore involves only simple (though tedious) algebra, and can easily be implemented in a computer algebra system.

Appendix C. Change of basis for the influence pseudotensors

As already argued in Sec. 4, literal expressions of the influence pseudotensors are computed in a convenient, local basis (e ′ j ) 1≤ j≤d . In order to retrieve the pseudotensors in the global basis (e i ) 1≤i≤d , a change of basis must be carried out, keeping in mind that the influence pseudotensors are not tensors. This is detailed below.

In the remainder of this section, it will be convenient to introduce the following notation: For x ∈ R d , let x • (resp. x ′ • ) be the coordinates of x in the global (resp. local) frame; x • is related to x ′

• through the classical formula

x i = R i j x ′ j ,
where R is the rotation tensor which maps the global basis onto the local basis (e ′ i = R • e i ). Then, from the multinomial theorem

x p i = R i j x ′ j p = k + =p p! k 1 ! • • • k d ! R i1 x ′ 1 k 1 • • • R id x ′ d k d ,
and

x p • • is expressed as a function of x ′ • x p • • = k •+ =p • p 1 ! • • • p d ! k 11 ! • • • k dd ! R k 11 11 • • • R k dd dd x ′k +• • .
Using the above identity, the variables x ′ • and y ′ • are substituted to x • and y • in Eq. ( 33)

T p • q • αβ (r) = h •+ =p • k •+ =q • p 1 ! • • • p d ! h 11 ! • • • h dd ! q 1 ! • • • q d ! k 11 ! • • • k dd ! R h 11 +k 11 11 • • • R h dd +k dd dd x,y∈R d χ α (x)χ β (y)x ′h +• • y ′k +• • Γ ∞ 0 (r + y -x)dx ′ × dy ′ × ,
where it is recalled that det R = 1. The integrals in the above sums are the influence pseudotensors T ′h +• k +• αβ (r), expressed in the local basis

T ′p • q • αβ (r) = x,y∈R d χ α (x)χ β (y)x ′p • • y ′q • • Γ ∞ 0 (r + y -x)dx ′ × dy ′ × ,
which shows that the influence pseudotensors in the global basis can be expressed as a linear combination of the influence pseudotensors in the local basis

T p • q • αβ (r) = r + =p + s + =q + Λ p • q • r • s • T ′r • s • αβ ,
where

Λ p • q • r • s • = h •+ =p • h +• =r • k •+ =q • k +• =s • p 1 ! • • • p d ! h 11 ! • • • h dd ! q 1 ! • • • q d ! k 11 ! • • • k dd ! R h 11 +k 11 11 • • • R h dd +k dd dd .
Finally, recalling that Γ ∞ 0 is a true tensor, and introducing the components T p • q • αβ,i jkl (r) [resp. T ′r • s • αβ,I JKL (r)] of T p • q • αβ (r) [resp. T ′r • s • αβ (r)] in the global (resp. local) basis, it is found that

T r • s • αβ,i jkl (r) = R iI R jJ R kK R lL r + =p + s + =q + Λ p • q • r • s • T ′r • s • αβ,I JKL .
The above formulas can be used to compute the influence tensor in any basis, provided all influence tensors have been computed in the local basis.

Figure 1 :

 1 Figure 1: Left: the N-inhomogeneity problem. Right: the N-inclusion problem.

Figure 4 :

 4 Figure 4: Variables for the calculation of the self-influence pseudotensor of inclusion Ω α .

Figure 5 :

 5 Figure 5: Variables for the calculation of the influence pseudotensor of inclusions Ω α and Ω β .

  For the sake of clarity and simplicity, closed-form expressions of the influence tensors were computed in a local basis (e ′ 1 , . . . , e ′ d ), where the d-th direction e ′ d coincides with the line of centers of the two inclusions (r = re ′ d ); the resulting pseudotensors are denoted T ′k • l • αβ (r). A change of basis must be performed to retrieve the components of T k • l • αβ in the global basis (e 1 , . . . , e d ); this is performed numerically, during the assembly of the global matrix. It should be noted that this change of basis does not reduce to a standard linear algebra operation, since T k • l •

Figure 6 :

 6 Figure 6: Left: graphical representation of one out of the 1 000 SVEs considered in Sec. 5.1. Circular, monodisperse pores are embedded in a homogeneous matrix; the porosity is φ = 0.4. Right: graphical representation of one out of the 100 SVEs considered in Sec. 5.2. Spherical, polydisperse pores are embedded in a homogeneous matrix (only the pores are shown in this image); the porosity is φ = 0.45. Order p Bound on µ app DOFs 0 0.310 µ 0 480 1 0.278 µ 0 1 440 2 0.257 µ 0 2 880 3 0.247 µ 0 4 800

Figure 7 :

 7 Figure 7: Relative error on µ app as a function of the number of degrees of freedom of the EIM-G, in semi-log scale. Convergence seems to be exponential.

  k •• denotes doubly indexed tuples of size d × d, k •• = (k i j ) 1≤i, j≤d . Given k •• , two d-tuples are further introduced k •+ = (k i1 + • • • + k id ) 1≤i≤d and k +• = (k 1 j + • • • + k d j ) 1≤ j≤d ,which are the horizontal and vertical projections of k •• , respectively.

Table 2 :

 2 Upper bounds on the apparent shear modulus of a polydisperse assembly of spherical pores in 3D elasticity, for increasing orders p of the EIM-G. For each value of p, the corresponding number of degrees of freedom (DOFs) is also reported.

	Order p	Bound on µ app	DOFs
	0	0.381	1 200
	1	0.371	4 800
	2	0.363	12 200

Preprint submitted to International Journal of Solids and StructuresDecember 30, 2013

The present work was carried out with the Maxima computer algebra system, version 5.30.0, http://maxima.sourceforge.net/ (last visited

2013-06- 05).

Appendix A. On the infinite body Green's operator for strains It is recalled that the Green's operator for strains Γ ∞ 0 decomposes into a regular part Q 0 and a singular part P 0 (Γ ∞ 0 * τ)(x) = lim δ→0 y∈Ω y-x ≥δ Q 0 (x-y) : τ(y)dV y +P 0 : τ(x), (A.1)

where P 0 is the Hill tensor of d-dimensional spheres. Litteral expressions of P 0 and Q 0 (r) can be found in reference textbooks [see e.g. [START_REF] Torquato | Random heterogeneous materials: microstructure and macroscopic properties[END_REF]; [START_REF] Buryachenko | Micromechanics of Heterogeneous Materials[END_REF]; [START_REF] Kanaun | Self-Consistent Methods for Composites[END_REF]]. For d = 3,

and

where µ 0 (resp. ν 0 ) is the shear modulus (resp. Poisson ratio) of the reference medium, r = r , n = r/r and J (resp. K) is the spherical (resp. deviatoric) fourth-rank projection tensor. More precisely, J = 1 d i ⊗ i and K = I -J, where i (resp. I) is the second-(resp. fourth-) rank identity tensor. For d = 2 (plane strain elasticity),

(A.5)

It should be observed that Q 0 is homogeneous of degree -d. Besides, the following identity will prove useful in the present paper

where dS n denotes the surface element at n on the unit sphere.

Appendix B. Computation of the self-influence pseudotensors

In this section, a procedure is presented for the computation of the self-influence pseudotensors of spherical inclusions. The general expression [see Eq. ( 42)] of the self-influence pseudotensors is first transformed so as to handle the principal values (see Appendix B.1). It is then shown that the resulting integral can be computed analytically in the case of spherical inclusions (see Appendix B.2).

Appendix B.1. Proof of Eq. (44)

These pseudotensors are defined by Eq. ( 42), where the principal value must be computed. Let

Assuming inclusion α is convex, it is possible to perform in Eq. (B.1) the change of variables specified by Eq. ( 43)

where dS n denotes the surface element at n on the unit sphere.

From the homogeneity of Q 0 (see Appendix A)

where the singularity is removed with the help of the following substitution

Indeed, Eq. ( 1) readily shows that the first term is regular as r → 0

Therefore,

and, from Eq. (A.6), the integral in the last term reduces to

while the change of variables specified by Eq. ( 43) can be reversed in the first term

This completes the proof, since Eq. ( 44) is readily retrieved from Eqs. (B.2), (B.3), (B.4) and (B.5).

Appendix B.2. Application to spherical inclusions

Let a α denote the radius of the spherical inclusion Ω α . It can readily be verified that

In expanded form, the self-influence pseudotensor of convex inclusions reads [see Eq. ( 44)]

The third term in Eq. (B.7) vanishes. Indeed, from the identities Q

From Eq. (B.6), R α (x, -n)R α (x, n) = a 2 αx 2 , and the above integral reduces to

where Eq. (A.6) has been used. Turning now to the second term of Eq. (B.7), Eq. (B.8) can again be used to transform the integral over the unit sphere

Using Eq. (B.6), it is readily seen that all odd powers of a 2 αx 2 + (x • n) 2 vanish in the integrand .10)


Theodoros Katsaounis
email: thodoros@tem.uoc.gr

Charalambos Makridakis

Chiara Simeoni
email: simeoni@math.unice.fr

Stability and convergence of relaxation finite element schemes for the incompressible Navier-Stokes equations

We are interested in the numerical modelling of the incompressible Navier-Stokes equations, in two space dimensions. The approach considered consists in developing semi-discrete finite element schemes for appropriate relaxation models, which formally arise as hyperbolic approximations of the Navier-Stokes equations. Stability properties of the relaxation finite element schemes are derived from estimating suitable modifications of the standard energy functional, that is suggested by the presence of relaxation terms. These techniques are also applied to prove error estimates and deduce the convergence of relaxation finite element schemes to the (smooth) solutions of the incompressible Navier-Stokes equations.

Introduction

We consider the incompressible Navier-Stokes equations, in two space dimensions,

∂ t u + u • ∇u -ν ∆u + ∇p = 0 , ∇ • u = 0 , t ∈ R + , x ∈ T 2 , (1)
which describe the flow of an incompressible fluid with kinematic viscosity ν > 0, where u(t, x) ∈ R 2 and p(t, x) ∈ R are the velocity and the pressure field, respectively. For the sake of simplicity, we assume the fluid moves in a two-dimensional torus T 2 ⊂ R 2 , to reduce to homogeneous boundary conditions.

We introduce the relaxation approximation of (1) considered in [START_REF] Brenier | On a relaxation approximation of the incompressible Navier-Stokes equations[END_REF], for ε > 0,

∂ t u ε + ∇ • U ε + ∇p ε = 0 , ∂ t U ε + ν ε ∇u ε = -1 ε (U ε -u ε ⊗ u ε ) , ∇ • u ε = 0 , t ∈ R + , x ∈ T 2 , (2)
which is obtained by singular perturbations of the nonlinear term in (1) through a supplementary (matrix valued) variable U ε (t, x) ∈ R 2 × R 2 , or rather, by means of diffusive scaling from standard relaxation approximations of the Euler equations, most in the spirit of [START_REF] Jin | The relaxation schemes for systems of conservation laws in arbitrary space dimensions[END_REF]. We mention that similar relaxation models are also proposed in [START_REF] Banda | Compressible and incompressible limits for hyperbolic systems with relaxation[END_REF] for numerical purposes, as reduced discrete velocity kinetic models, and these approximations are systematically used to study the relaxation limits of hyperbolic conservation laws (see [START_REF] Donatelli | Convergence of singular limits for multi-D semilinear hyperbolic systems to parabolic systems[END_REF] and [START_REF] Tzavaras | Materials with internal variables and relaxation to conservation laws[END_REF], for instance). We recall the one-field equation associated with the relaxation system (2), namely

∂ t u ε + u ε • ∇u ε -ν ∆u ε + ε ∂ tt u ε + ∇p ε + ε ∂ t ∇p ε = 0 , ∇ • u ε = 0 , t ∈ R + , x ∈ T 2 , (3)
which turns out to be a hyperbolic approximation of the Navier-Stokes equations (1) with a small dissipative correction (that is viewed as part of a wave operator). We refer to [START_REF] Brenier | On a relaxation approximation of the incompressible Navier-Stokes equations[END_REF] for the rigorous justification of the convergence, as ε → 0, toward the incompressible Navier-Stokes equations, by the "modulated energy" method [START_REF] Brenier | Convergence of the Vlasov-Poisson system to the incompressible Euler equations[END_REF].

We deal with a numerical approach to the Navier-Stokes equations (1) based on standard finite element schemes for the relaxation approximation (3). Let T h be a (uniformly regular) partition of R 2 , with characteristic parameter h > 0. The discrete velocity u h and pressure p h are determined in finite element spaces denoted by V h ⊂ H 1 (T 2 ) and Π h ⊂ L 2 (T 2 ), respectively. These spaces are assumed to satisfy typical approximability properties, namely for V h consisting of piecewise polynomials (globally continuous) on T h such that [START_REF] Brenier | On a relaxation approximation of the incompressible Navier-Stokes equations[END_REF] inf

v∈V h [ ||u -v|| L 2 + h ||u -v|| H 1 ] ≤ C h k ||u|| H k , u ∈ H k (T 2 ) , k ≥ 2 .
We focus on spatial discretizations of (3), the so-called "semi-discrete schemes" (for the computational issues, a semi-implicit discretization in time is eventually considered in [START_REF] Ch | Stability and convergence of a class of finite element schemes for hyperbolic systems of conservation laws[END_REF], for which some analogue of the estimates below can also be proved).

In a first attempt, we only address the question of establishing error estimates for the velocity field, therefore we introduce the subspace Z h ⊂ V h of approximate divergence-free functions, given by ( 5)

Z h = { v ∈ V h ; ∇ • v , q = 0 , ∀ q ∈ Π h } ,
where the inner-product . , . is defined for the L 2 -norm.

According to [START_REF] Brenner | The mathematical theory of finite element methods Second Edition[END_REF], we take the variational formulation of (3) into account and the corresponding finite element scheme, for the numerical solution u h ∈ Z h , reads

∂ t u h , v + N (u h , u h ), v + ν ∇u h , ∇v + ε ∂ tt u h , v = 0 , ∀ v ∈ Z h , (6)
N (u h , u h ) = u h • ∇u h + 1 2 (∇ • u h ) u h , (7)
and the modified numerical nonlinear term ( 7) is introduced to keep the conservation property N (u h , u h ), u h = 0 in the finite element framework (we remark that the divergence-free constraint holds only weakly in the finite element space Z h ).

"Sub-characteristic" condition and stability estimates

The qualitative analysis of finite element schemes is performed by applying similar techniques as for the corresponding continuous problems (we refer to [START_REF] Girault | Finite element methods for Navier-Stokes equations[END_REF], [START_REF] Temam | Navier-Stokes equations. Theory and numerical analysis[END_REF] and [START_REF] Heywood | Finite element approximation of the nonstationary Navier-Stokes problem. I. Regularity of solutions and second-order error estimates for spatial discretization[END_REF] for standard finite element approximations of the Navier-Stokes equations).

To deduce stability properties of ( 6)-( 7), we study an appropriate modification of the standard L 2 -energy functional, whose form is suggested by the presence of the relaxation terms (see [START_REF] Brenier | On a relaxation approximation of the incompressible Navier-Stokes equations[END_REF], [START_REF] Ch | Stability and convergence of a class of finite element schemes for hyperbolic systems of conservation laws[END_REF] and [START_REF] Tzavaras | Materials with internal variables and relaxation to conservation laws[END_REF], for instance).

Lemma 2.1. Under the following assumption on the numerical solution,

(8) u h L ∞ ≤ 2 ν 5 ε , ∀ t ∈ R + ,
the relaxation finite element scheme (6)-( 7) satisfies the dissipation estimate

(9) 1 2 u h + ε ∂ t u h 2 L 2 + 1 2 ε 2 ∂ t u h 2 L 2 + ε ν ∇u h 2 L 2 ≤ C(u 0 h ) ,
for some constant C(u 0 h ) only depending on the approximation of initial data of (1).

Proof. We take v = u h ∈ Z h in (6) to obtain (10) 1 2

d dt u h 2 L 2 + ν ∇u h 2 L 2 + ε d dt (∂ t u h , u h ) -ε ∂ t u h 2 L 2 = 0 .
To estimate the time derivative, we take v = ∂ t u h ∈ Z h in ( 6) and we obtain (11)

∂ t u h 2 L 2 + N (u h , u h ), ∂ t u h + 1 2 ν d dt ∇u h 2 L 2 + 1 2 ε d dt ∂ t u h 2 L 2 = 0 .
Then, we multiply (11) by 2 ε (as a correction by lower order effects) and we combine with [START_REF] Th | Relaxation methods and finite element schemes for the incompressible Navier-Stokes equations[END_REF], rearranging the resulting equality as follows, 1 2

d dt u h + ε ∂ t u h 2 L 2 + 1 2 ε 2 d dt ∂ t u h 2 L 2 + ε ν d dt ∇u h 2 L 2 = -ν ∇u h 2 L 2 -ε ∂ t u h 2 L 2 -2 ε N (u h , u h ), ∂ t u h = -ε ∂ t u h + N (u h , u h ) 2 L 2 -ν ∇u h 2 L 2 + ε N (u h , u h ) 2 L 2 . (12)
Recalling [START_REF] Girault | Finite element methods for Navier-Stokes equations[END_REF], the last two terms of ( 12) are bounded as follows,

(13) -ν ∇u h 2 L 2 + ε N (u h , u h ) 2 L 2 ≤ -ν - 5 2 ε u h 2 L ∞ ∇u h 2 L 2 .
The result ( 9) is thus achieved, after integrating [START_REF] Tzavaras | Materials with internal variables and relaxation to conservation laws[END_REF] with respect to time, by using the uniform bound (8) to conclude the estimate of (13).

We remark that the a priori L ∞ -bound ( 8) is equivalent to the "sub-characteristic condition" required in [START_REF] Ch | Stability and convergence of a class of finite element schemes for hyperbolic systems of conservation laws[END_REF] to derive stability estimates for relaxation finite element approximations of hyperbolic conservation laws, that reveals the regularization role of the wave operator in (3).

The "sub-characteristic condition" in ( 8) is actually satisfied for the relaxation finite element scheme ( 6)- [START_REF] Girault | Finite element methods for Navier-Stokes equations[END_REF], because of some "inverse-inequalities" which hold in the finite element spaces, namely, in two space dimensions, ( 14)

u h L ∞ ≤ C h -1 u h L 2 , ∀ t ∈ R + .
Such inequalities, obviously degenerate with respect to the characteristic parameter h, are proved for several classes of finite element spaces (refer to [START_REF] Girault | Finite element methods for Navier-Stokes equations[END_REF] and [START_REF] Brenner | The mathematical theory of finite element methods Second Edition[END_REF], for instance), also including the framework [START_REF] Brenner | The mathematical theory of finite element methods Second Edition[END_REF]. Therefore, it follows from ( 14) that a sufficient condition for (8) to be satisfied is

(15) u h L 2 ≤ C h 2 ν 5 ε , ∀ t ∈ R + ,
that is guaranteed by the stability estimate [START_REF] Jin | The relaxation schemes for systems of conservation laws in arbitrary space dimensions[END_REF], provided that (15) holds for an approximation of the initial data of (1), and with (at least) ε = O(h 2 ) as the ratio between the diffusion parameters.

Error estimate for the velocity

We consider standard arguments to derive error estimates for finite element schemes.

According to [START_REF] Girault | Finite element methods for Navier-Stokes equations[END_REF], the technique is based on splitting the error between the solution to the Navier-Stokes equations (1) and the numerical solution obtained by applying the relaxation finite element scheme ( 6)-( 7) into two parts, as follows,

-u h = (u -P h u) + (P h u -u h ) , (16) u
where P h is some appropriate projection operator into the finite element space [START_REF] Brenner | The mathematical theory of finite element methods Second Edition[END_REF].

The consistency error η(u) = (u-P h u) is interpreted in terms of approximability properties of the finite element spaces, according to the definition of P h u ∈ Z h . In a first attempt, we introduce the "elliptic projection" associated with the linear stationary part of the Navier-Stokes equations (1), namely (17)

∇ P h u , ∇v = ∇u , ∇v , ∀ v ∈ Z h ,
for which some analogue of (4) in the framework (5) also holds (see [START_REF] Brenner | The mathematical theory of finite element methods Second Edition[END_REF], for instance). We thus deduce that

(18) η(u) L 2 ≤ C h k ||u|| H k , u ∈ H k (T 2 ) , k ≥ 2 .
Taking the variational formulation of ( 1) and ( 6) into account, together with (17), we establish a variational equation for the stability error ξ(u, u h ) = (P h u-u h ) ∈ Z h in (16), namely

∂ t ξ , v + u • ∇u -N (u h , u h ) , v + ν ∇ξ , ∇v + ε ∂ tt ξ, v = -∂ t η , φ + ε ∂ tt P h u , φ -∇p, v , ∀ v ∈ Z h . (19)
We make similar computations as for the proof of the L 2 -stability estimate (9) in Lemma 2.1 and, in the right-hand side of (19), we use the consistency properties (18) of the elliptic projection, to deduce the following convergence result. Theorem 3.1. Under the assumption (8) on the numerical solution, the relaxation finite element scheme (6)-( 7) satisfies the following error estimate, for k ≥ 2,

(20) u -u h L 2 ≤ C exp( 1 ν ) h k , ∀ t ∈ R + ,
with some constant C depending on smoothness properties of the solution to (1). We remark that is the treatment of the incompressibility condition which leads to the singular dependence of the error constants in (20) with respect to the diffusion parameters. That issue is particular important in the numerical analysis of the incompressible Navier-Stokes equations, as the question of dealing with the divergence-free condition rigorously is still an open problem. The main idea of relaxation finite element schemes for the incompressible Navier-Stokes equations is to achieve stability also in the limit of small viscosity, when the transport terms become dominant and we recover the framework of hyperbolic conservation laws (we actually consider relaxation approximations of the incompressible Euler equations), for which these techniques have been successfully applied in [START_REF] Ch | Stability and convergence of a class of finite element schemes for hyperbolic systems of conservation laws[END_REF]. In that sense, further works are in progress for some different relaxation models (see [START_REF] Donatelli | Convergence of singular limits for multi-D semilinear hyperbolic systems to parabolic systems[END_REF], for instance).

We refer to [START_REF] Th | Relaxation methods and finite element schemes for the incompressible Navier-Stokes equations[END_REF] for the details of the proof of Theorem (3.1) and the more general statement of the subjects covered in these notes.

Acknowledgments

The authors would like to thank Prof. Athanasios E. Tzavaras and Prof. P.E. Jabin for several helpful discussions. This work was partially supported by the HYKE European Network HPRN-CT-2002-00282 and by the MCWAVE Marie Curie Host Fellowship at the I.A.C.M. -F.O.R.T.H. Heraklion, Crete (Greece).