Second order homogenization of quasi-Periodic structures

Ω;R3) which minimizes the energy functional J  over H10 (Ω;R3): J (u) = min v∈H10 (Ω;R3) J (v) with J (v) = ∫ Ω ( 1 2 Aεx(v) · εx(v)− f · v ) dx. (5) As usual for the asymptotic homogenization approaches, the micro-cell Y  is re-scaled by introducing the “unit” cell Y so that Y  = Y and we assume that |Y | = 1. In the quasi-periodic case, the properties of a point in the cell depends not only on its position but also the position of the cell in the medium. Specifically, the local elasticity tensor Aijkl in the unit cell can be represented by a function Aijkl(x,y) which is assumed to be smooth in the first variable x ∈ Ω and - eriodic in th second one y ∈ Y : Aijkl( ) ijkl(x, x  ). (6) Ω  Fig. 1. A typical quasi-periodic medium considered in the paper, where the inclusions have all the same size but a stiffness which varies smoothly with their location x in the domain Ω. Remark 1. This type of quasi-periodic medium is a particular case of a more general class of heterogeneous media where the microstructure varies smoothly. Indeed, here we assume that the cell Y is always the same, independent of x. A more general case would consist in considering that Y depends on x. For instance, one could consider that the dependence of the stiffness on the position reads as A(x) = A ( x, ϕ(x)  ) , where ϕ is a smooth diffeomorphism. In such a case, assuming that A is still periodic with respect to the second variable with period Y , the image of Y by ϕ−1 will depend on x. Fig. 1. A typical quasi-periodic mediu sidered in the pa er, where the inclusions have all th same size but a tiffness which varies smoothly wit their location x i the domain Ω Remark 1. This type of quasi-periodic medium is a particular case of a more general class of heterogeneous media where the microstructure varies smoothly. Indeed, here we assume that the cell Y is always the same, independent of x. A more general case would consist in considering that Y depends on x. For instance, one could consider that the dependence of the stiffness on the position reads as Ae(x) = A ( x, ϕ(x) e ) , whereϕ is a smooth diffeomorphism. In such a case, assuming that A is still periodic with respect to the second variable with period Y, the image of eY byϕ−1 will depend on x. Second order homogenization of quasi-periodic structures 329 By using the traditional multi-scale method like in [1, 15], the displacement field ue is represented through an asymptotic expansion in terms of the small parameter e ue(x) = u0 ( x, x e ) + eu1 ( x, x e ) + e2u2 ( x, x e ) + e3u3 ( x, x e ) + . . . (7) where each term ui(x,y) is a function of the two variables x and y, periodic with respect to the “fast” variable y = (y1, y2, y3) = x/e with period Y. Accordingly, the stress field σe admits the same type of expansion σe(x) = σ0 ( x, x e ) + eσ1 ( x, x e ) + e2σ2 ( x, x e ) + e3σ3 ( x, x e ) + . . . (8) Remark 2. From a regularity point of view, the solution ue must be of finite energy and hence must belong to H10 ( Ω;R3 ) . Therefore, one could expect that each term ui of the expansion has the same regularity and hence that ui ∈ H1 (Ω×Y;R3). Moreover, the ui’s must be periodic on Y and one requires that they vanish on ∂Ω× Y. That leads to introduce the following set U at which each ui a priori should belong U = { v ∈ H1 (Ω×Y;R3) : v = 0 on ∂Ω×Y,v is Y-periodic } . However, we will see during the asymptotic procedure that, except u0, the subsequent ui’s are in general less regular than expected and that they can satisfy the boundary condition on ∂Ω only under particular conditions on the microstructure. We will make some comments concerning this lack of regularity, but the study of the boundary layer effects due to the loss of the boundary conditions is outside the scope of the present paper, see Dumontet [4], Devries et al. [16] where boundary layer correctors are calculated in the periodic case. Inserting the expansion (7) into the energy functional Je (ue), and taking into account the following derivation et integration rules ∂xl → ∂xl + 1 e ∂yl , ∫ Ω (·)dx→ 1|Y| ∫ Ω×Y (·)dxdy. (9) Since |Y| = 1, the factor 1/|Y| in front of the integral may be omitted. The energy func- tional Je (ue) becomes a series in e Je (ue) = Je∗ ( u0,u1,u2, . . . ) := 1 e2 J(−2)∗ + 1 e J(−1)∗ + e (0) ∗ + eJ (1) ∗ + e2 J (2) ∗ + . . . (10) where J(p)∗ = J(p) ( u0, . . . ,up+1 ) , (11) the functionals J(p) being defined for p = −2 by J(−2) ( v0 ) = ∫ Ω×Y 1 2 Aεy ( v0 ) · εy (v0)dxdy, (12) 330 Duc Trung Le, Jean-Jacques Marigo and for p ≥ −1 by J(p) ( v0, . . . ,vp+1 ) := 1 2 ∑m+n=p ∫ Ω×Y A [ εx (vm) + εy ( vm+1 )] · [ εx (vn) + εy ( vn+1 )] dxdy − ∫ Ω×Y f · vpdxdy. (13) Accordingly, the minimization problem (5) leads to the following Sequence of ordered minimization problems. Find, for i ∈ N, the functions ui ∈ U , which minimize the energy functional Je ( v0,v1, . . . ) = ∑ p≥−2 ep J(p) ( v0, . . . ,vp+1 ) , over all vi ∈ U . In other words the formal asymptotic procedure leads to replace (5) by a sequence of minimization problems which are detailed in the next section. Let us note that the different terms of the stress expansion are related to those of the displacement expansion by the following chain rule σi = A ( εx ( ui ) + εy ( ui+1 )) , ∀i ≥ 0. (14) 3. SOLVING OF THE SEQUENCE OF ORDEREDMINIMIZATION PROBLEMS 3.1. J(−2) and J(−1) minimization problems The field u0 has to minimize J(−2)(v0) over all v0 ∈ U . Since the stiffness tensor A is definite positive, we have J(−2)(v0) ≥ 0 with the equality only if εy(v0) = 0 and hence only if v0 is independent of y. Therefore, the minimization of J(−2) requires that u0 is function of the macroscopic variable x only u0 = u0(x), (15) and hence u0 must belong to H10(Ω;R 3). In turn, considering v0 such that v0 = v0(x) leads to J(−1)(v0,v1) = 0, regardless of v1 Therefore, for such v0, the energy starts at the order 0 Je ( v0,v1, · · · ) = J(0) ( v0,v1 ) + eJ(1) ( v0,v1,v2 ) + . . . 3.2. J(0) minimization problem Considering v0 independent of y, the energy J(0)(v0,v1) reads as J(0) ( v0,v1 ) = ∫ Ω×Y 1 2 A [ εx ( v0 ) + εy ( v1 )] · [ εy ( v0 ) + εy ( v1 )] dxdy− ∫ Ω×Y f · v0dxdy. (16) We have to find u0 ∈ H10(Ω;R3) and u1 ∈ U which minimize J(0). We proceed in two steps: (1) First, for a given v0 ∈ H10(Ω;R3) , we minimize J(0)(v0, ·) over U . Second order homogenization of quasi-periodic structures 331 (2) The so-obtained minimum becomes a functional of v0 that we minimize to ob- tain u0. 3.2.1. Minimizing with respect to v1 at given v0 At given v0, since J(0)(v0,v1) is a convex functional of v1, its minimizer, denoted by v1 and which depends on v0, is such that the first derivative of J(0)(v0, ·) vanishes. That leads to the following variational equation for v1∫ Ω×Y A [ εx ( v0 ) + εy ( v1 )] · εy(v)dxdy = 0, (17) which must hold for any v ∈ U . Let us choose v under the form v(x,y) = w(x)φ(y) with w ∈ D(Ω) and φ ∈ H1#(Y;R3). Inserting this test function in (17) and using classic argument of calculus of variation, we get (at almost all x ∈ Ω)∫ {x}×Y A [ εx ( v0 ) + εy ( v1 )] · εy(φ)dy = 0, ∀φ ∈ H1# ( Y;R3 ) . (18) Since (18) brings into play the derivative of v1 with respect to y only, v1 is determined at this stage up to an arbitrary function of x. Moreover, by linearity, v1 can be written as follows (where the convention of summation is still used): v1(x,y) = εxpq ( v0 ) (x)χpq(x,y) +V1(x). (19) In (19), V is the arbitrary function of x and the six1 vector fields χpq are solution of the family of six elementary problems posed on Y and indexed by x, each one corresponding to the response of the unit cell to a prescribed macroscopic strain tensor. Specifically, at given x, the Y-periodic field y 7→ χpq(x,y) is solution of the following variational equation ∫ {x}×Y A [ ep  eq + εy (χpq) ] · εy(φ)dy = 0, ∀φ ∈ H1# (Y;R3) , (20) where  denotes the symmetrized tensorial product and ei stands for the ith basis vec- tor of the cartesian coordinates. Since (20) determines χpq up to a function of x only, the indetermination is removed by adding the condition that the average of χpq over Y vanishes: 〈χpq〉 = 0 with 〈·〉 = 1|Y| ∫ Y ·dy. (21) Accordingly, V1 corresponds to the average of v1 over Y but remains undetermined at this stage. Let us note that the fields χpq do not depend on v0 but depend on the quasi periodic repartition of the microstructure only. To compare with what happens in the case of a periodic medium, the solution of the cell problems now depends not only on the microscopic variable y but also on the macroscopic one x. 1The number of vector fields to be determined is reduced to six because of the obvious symmetry χpq = χqp. 332 Duc Trung Le, Jean-Jacques Marigo 3.2.2. Minimizing with respect to y0 Using (19), the energy J(0)(v0,V1) becomes the following functional J(0) defined for v0 ∈ H10(Ω;R3) by J(0) ( v0 ) = ∫ Ω 1 2 A(0)εx ( v0 ) · εX (v0)dx− ∫ Ω f · v0dx. (22) In (22) A(0)(x) is the (classical) homogenized stiffness tensor at x which is symmetric, pos- itive definite and depends on x because of the quasi-periodicity assumption. Specifically, owing to (20), the components of A(0)(x) can be written under the two following forms A (0) mnpq(x) = ∫ {x}×Y A [ em  en + εy (χmn) ] · [ep  eq + εy (χpq)]dy (23) = ∫ {x}×Y Aijpq [ δimδjn + εyij (χ mn) ] dy. (24) The symmetry and the positivity of A(0) appear more clearly from (23), but the form (24) will be useful later. Therefore J(0) is convex and coercive on H10 ( Ω;R3 ) . So its minimizer u0 is unique and such that∫ Ω ( A(0)εx ( u0 ) · εx(v)− f · v)dx = 0, ∀v ∈ H10 (Ω;R3) . (25) Accordingly, u0 is the unique solution of the following linear elastic problem defined on the “homogenized” body (that is to say, the body whose microstructure has been removed and replaced by its effective stiffness A(0)){ divx ( A(0)εx ( u0 )) + f = 0, in Ω u0 = 0, on ∂Ω (26) Once this problem is solved, u0 is known and we obtain from (19) that the second term u1 of the expansion of ue can read as u1(x,y) = εxpq ( u0 ) (x)χpq(x,y) + 〈 u1 〉 (x), (27) its average value 〈u1〉 only remaining unknown at this stage. Moreover, the first order stress field σ0 is also known and given by σ0(x,y) = A(x,y) ( εx ( u0 ) (x) + εxpq ( u0 ) (x)εy (χpq) (x,y) ) . (28) Let us note that σ0 depends linearly on the macroscopic strain εx(u0), σ0(x,y) = a0(x,y)εx ( u0 ) , (29) where the fourth-order stiffness tensor field a0 (sometimes called the stress localization field) depends only on the microstructure a0ijkl = Aijkl + Aijpqεypq ( χkl ) . (30) Second order homogenization of quasi-periodic structures 333 The macroscopic stress field of order 0 is the average value of σ0 on Y and, by virtue of (24), is associated to the strain field of order 0 by the effective relation of order 0〈 σ0 〉 = A(0)εx ( u0 ) . Therefore A(0) is nothing but that the average value of a0: A(0) = 〈a0〉. Let us make some comments concerning the regularity and the boundary conditions satisfied by u1 in the spirit of Remark 2. Remark 3. (i) In order that u1 is in H1 ( Ω;R3 ) , it is necessary that u0 is in H2(Ω;R3) which is the case if the effective stiffness tensor A(0) and the body forces f are smooth functions of x. (ii) The verification of the boundary condition u1 = 0 on ∂Ω × Y is more delicate. Indeed, in general εx(u0) does not vanish on ∂Ω. Therefore, the unique possibility for u1 to satisfy the boundary condition is that the six fields χpq vanish on ∂Ω× Y (otherwise boundary layer effects will appear). This is the case when the body is homogeneous in a neighborhood of the boundary, i.e. when A does not depend on y when x is close to ∂Ω. In such a case it suffices that 〈u1〉 = 0 on ∂Ω in order that the boundary condition is satisfied. In the sequel, we will assume that all the conditions are satisfied in order that u1 actually belongs to H10(Ω;R 3) . 3.3. Calculation of J(1)(u0,u1,u2) By the definition (13) one gets J(1) ( u0,u1,u2 ) = ∫ Ω×Y ( A [ εx ( u0 ) + εy ( u1 )] · [ εx ( u1 ) + εy ( u2 )]− f · u1)dxdy. (31) Using (18) with (u0,u1,u2) in place of (v0,v1,v) one obtains∫ Ω×Y A [ εx ( u0 ) + εy ( u1 )] · εy ( u2 ) dxdy = 0. (32) Using (25) with v = 〈u1〉 gives∫ Ω ( A(0)εx ( u0 ) · εx (〈u1〉)− f · 〈u1〉)dx = 0. Therefore, using the form (24) for A(0), the previous equality reads also as∫ Ω×Y ( A [ εx ( u0 ) + εy ( u1 )] · εx (〈 u1 〉) − f · 〈 u1 〉) dxdy = 0. (33) Inserting (32) and (33) into (31) shows that J(1)(u0,u1,u2) depends in fact only on u0, but not on 〈u1〉 and u2. Therefore, J(1)∗ := J(1)(u0,u1,u2) can be considered as known at this stage. Furthermore, let us remark that the energy J(1)∗ can be seen as a quadratic form of u0 which involves new homogenized stiffness tensors that we propose to identify. Indeed, 334 Duc Trung Le, Jean-Jacques Marigo since f does not depend on y, the last term in (31) can read as∫ Ω×Y f · u1dxdy = ∫ Ω f · 〈 u1 〉 dx. Then, using (28), (32) and (33), (31) becomes J(1)∗ = ∫ Ω×Y σ0 · εx ( εxij ( u0 ) χij ) dxdy. For a future use, it is more convenient to make an integration by parts of the integrand to obtain J(1)∗ = − ∫ Ω×Y divx σ0 · χijεxij ( u0 ) dxdy, the boundary term on ∂Ω×Y cancelling by virtue of Remark 3. Using (29) and expanding the derivative with respect to x lead to divx σ0 = ∂xq ( a0pqkl ) εxkl ( u0 ) ep + a0pqkl gxklq ( u0 ) ep, (34) where gx(u) denotes the gradient of the symmetrized gradient of u (all the derivatives are with respect to x) whose components are given by gxijk(u) = ∂xk ( εxij(u) ) = 1 2 ( ∂2ui ∂xj∂xk + ∂2uj ∂xi∂xk ) . (35) Inserting (34) into J(1)∗ leads to J(1)∗ = ∫ Ω ( A(1)εx ( u0 ) · εx (u0)dx+ B(0)gx (u0) · εx (u0))dx, (36) where appear the fourth-order tensor field A(1) and the fifth-order tensor field B(0) de- fined by A (1) ijkl(x) = − ∫ {x}×Y ∂xq ( a0pqkl ) χ ij pdy, B (0) ijklq(x) = − ∫ {x}×Y a0pqklχ ij pdy. (37) These two tensor fields depend only on the microstructure and can be considered as two other effective stiffness tensor fields. Let us note also that A(1) has the minor symmetries, not the major one (but its symmetric part only is involved in (36)) A (1) ijkl = A (1) jikl = A (1) ijlk, A (1) ijkl 6= A(1)klij, and that B(1) has the following symmetries B (0) ijklm = B (0) jiklm = B (0) ijkm. It remains to determine 〈u1〉. It will be given by minimizing J(2). Second order homogenization of quasi-periodic structures 335 3.4. J(2) minimization problem Let us consider J(2)(u0,v1,v2,v3) with v1 = εxpq ( u0 ) χpq + 〈 v1 〉 , 〈 v1 〉 ∈ H10 ( Ω;R3 ) and v2,v3 in U . Note that y1 differs from u1 only by its average value and hence that εy(v1) = εy(u1). By virtue of (13), one gets J(2) ( u0,v1,v2,v3 ) = ∫ Ω×Y 1 2 A [ εx ( v1 ) + εy ( v2 )] · [εx (v1)+ εy (v2)]dxdy+ + ∫ Ω×Y A [ εx ( u0 ) + εy ( u1 )] · [εx (v2)+ εy (v3)]dxdy− ∫ Ω×Y f · v2dxdy. (38) Using (18) with (u0,u1,v3) in place of (v0,v1,v) one obtains∫ Ω×Y A [ εx ( u0 ) + εy ( u1 )] · εy ( v3 ) dxdy = 0, (39) and hence J(2)(u0,v1,v2,v3) does not depend on the choice of v3 in U . By virtue of (25), see also (33), one has∫ Ω×Y ( A [ εx ( u0 ) + εy ( u1 )] · εx (〈 v2 〉)− f · 〈v2〉)dxdy = 0, and hence, since u0 is known and since εy(〈v2〉) = 0, J(2)(u0,v1,v2,v3) only depends on 〈v1〉 and v2 − 〈v2〉. Therefore, the problem consists in minimizing, for u0 given in H10(Ω;R 3) , the following quadratic functional J(2) defined on H10(Ω;R 3)×U0 J¯(2)(〈v〉,w) = ∫ Ω×Y A [ εx ( u0 ) + εxpq ( u0 ) εy (χ pq) ] · εx(w)dxdy+ + ∫ Ω×Y 1 2 A [ εx(〈v〉) + εy(w) + εx ( εxpq ( u0 ) χpq )] · [ εx(〈v〉) + εy(w) + εx ( εxpq ( u0 ) χpq )] dxdy, (40) where U0 denotes the linear subset of U made of fields defined on Ω× Y whose average value with respect to Y vanishes (almost) everywhere in Ω U0 = {v ∈ U : 〈v〉(x) = 0, ∀x ∈ Ω}. (41) The minimizer of J(2) is (〈u1〉,u2 − 〈u2〉) (which means that u2 will be determined up to a function of x at this stage). To minimize J(2), we proceed in three steps: (1) First, for a given 〈v〉 ∈ H10(Ω;R3) , we minimize J(2)(〈v〉, ·) over U0. (2) The so-obtained minimum becomes a functional of 〈v〉 defined on H10(Ω;R3) which involves six effective stiffness tensor fields. (3) Finally, minimizing that functional gives 〈u1〉 and hence u2 − 〈u2〉. 336 Duc Trung Le, Jean-Jacques Marigo 3.4.1. Minimizing with respect to w at given 〈v〉 At given 〈v〉, since J(2)(〈v〉,w) is a convex functional of w, its minimizer, denoted by W and which depends on u0 and 〈v〉, is such that the first derivative of J(2)(〈v〉, ·) vanishes. Using the relation (28) for σ0, that leads to the following variational equation for w 0 = ∫ Ω×Y ( A [ εy(w) + εx(〈v〉) + εx ( εxpq ( u0 ) χpq )] · εy(w) + σ0 · εx(w))dxdy, (42) which must hold for any w ∈ U0. Because of the presence of the gradient of w with respect to x in the last term of (42) we cannot directly use classical arguments to suppress the integral over Ω for obtaining local problems posed at almost all points x of Ω. We must before make an integration by parts of that term. So, integrating by parts the last term of (42) gives 0 = ∫ Ω×Y ( A [ εy(w) + εx(〈v〉) + εx ( εxpq ( u0 ) χpq )] · εy(w)− divxσ0 ·w)dxdy, (43) the boundary term on ∂Ω cancelling because w ∈ U . Let us choose w under the form w(x,y) = ϕ(x)φ(y) with ϕ ∈ D(Ω) and φ ∈ H1#0 ( Y;R3 ) , where H1#0(Y;R 3) denotes the Y-periodic fields with zero average value H1#0 ( Y;R3 ) = { φ ∈ H1# ( Y;R3 ) : 〈φ〉 = 0 } . (44) Inserting this test function into (43), we can now use classical argument of calculus of variation to get (at almost all x ∈ Ω):∫ {x}×Y A [ εy(w) + εx(〈v〉) + εx ( εxpq ( u0 ) χpq )] · εy(φ)dy = ∫ {x}×Y divxσ0 ·φdy, (45) which must hold for any φ ∈ H1#0(Y;R3). Let us note that both εx(εxpq(u0)χpq) and divxσ0 depend linearly on εx(u) and its gradient gx(u0) . The dependency of the latter is given by (34) whereas for the former one gets εxij ( εxkl ( u0 ) χkl ) = εxkl ( u0 ) εxij ( χkl ) + 1 2 ( gxkli ( u0 ) χklj + gxklj ( u0 ) χkli ) . (46) Using these properties of linearity allow us to decompose w as the combination of fields depending on the microstructure only. Specifically, w depends linearly on the six (inde- pendent) components of εx(〈v〉) , the six (independent) components of εx(u0) and the eighteen (independent) components of gx(u0), and hence can be read as w(x,y) = εxpq(〈v〉)(x)χpq(x,y) + εxpq ( u0 ) (x)ψpq(x,y) + gxpqr ( u0 ) (x)ξpqr(x,y). (47) In (47), the quasi periodic displacement fields associated with εx(v) are the same fields χpq as those defined by (20) and (21), whereas the other ones, ψpq and ξpqr, associated Second order homogenization of quasi-periodic structures 337 with the strain εx(u0) and its gradient gx(u0) respectively, are new quasi periodic fields which are solution in U0 of the following problems defined on the unit cell {x} ×Y: Find ψpq ∈ H1#0 ( Y;R3 ) such that ∀φ ∈ H1#0 ( Y;R3 ) ,∫ {x}×Y A [ εy (ψ pq) + εx (χ pq) ] · εy(φ)dy = ∫{x}×Y ∂xj ( a0ijpq ) φidy; (48) Find ξpqr ∈ H1#0 ( Y;R3 ) such that ∀φ ∈ H1#0 ( Y;R3 ) ,∫ {x}×Y A [ εy (ξ pqr) + χpq  er ] · εy(φ)dy = ∫{x}×Y a0irpqφidy. (49) Remark 4. Concerning the regularity and the boundary conditions satisfied by w: (i) The existence of ψpq and ξpqr is not guaranteed but depends on the regularity of the fields χpq on x and hence on the microstructure. This regularity must be checked on a case by case basis. On the other hand, the uniqueness is ensured once the existence is established. (ii) Even if a solution ψpq and ξpqr exists for each “elementary problem”, the field w belongs to H1(Ω×Y;R3) only if u0 and 〈v〉 are smooth enough. (iii) The verification of the boundary condition w = 0 on ∂Ω× Y is as delicate as for u1. Indeed, in general neither εx(〈v〉) , nor εx(u0) and nor gx(u0) vanish on ∂Ω. Therefore, the unique possibility for w to satisfy the boundary condition is that all the fields χpq,ψpq and ξpqr vanish on ∂Ω×Y. This is still the case when the body is homogeneous in a neighborhood of the boundary. In the sequel, we will assume that all the conditions are satisfied in order that w exists and belongs to H10(Ω;R 3). 3.4.2. Expression of J(2) Inserting (47) into (40) obtains a functional J (2) of 〈v〉 defined on H10 ( Ω,R3 ) . The functional J (2) contains several new effective stiffness tensor fields. Specifically, after some tedious calculations which are not detailed here, one eventually gets J (2) (〈v〉) = ∫ Ω ( 1 2 A(0)εx(〈v〉) · εx(〈v〉) + ( A(1)εx ( u0 ) + B(0)gx ( u0 )) · εx(〈v〉))dx + ∫ Ω ( 1 2 A(2)εx ( u0 ) · εx (u0)+ B(1)gx (u0) · εx (u0)+ 12C(0)gx (u0) · gx (u0) ) dx, (50) where A(0),A(1) and B(0) are the already defined stiffness tensor fields, see (23) and (37), whereas A(2),B(1) and C(0) are the new stiffness tensor fields given by A (2) ijkl = ∫ Y ( Aεx ( χij ) · εx ( χkl ) − Aεy ( ψij ) · εy ( ψkl )) dy, (51) B (1) ijklm = ∫ Y ( Aεx ( χij ) · ( χkl  em ) − Aεy ( ψij ) · εy ( ξklm )) dy, (52) C (0) ijklmn = ∫ Y ( A ( χij  ek ) · ( χlm  en ) − Aεy ( ξijk ) · εy ( ξlmn )) dy. (53) 338 Duc Trung Le, Jean-Jacques Marigo Note that one must use the elementary problems giving the χij’s, ψij’s and ξklm to obtain the expressions above. As the three previous ones, these stiffness tensor fields depend only on the microstructure and admit the same type of symmetries. Specifically, one has A (2) ijkl = A (2) ijlk = A (2) klij, B (1) ijklm = B (1) jiklm = B (1) ijlkm, C (0) ijklmn = C (0) jiklmn = C (0) lmnijk. The properties of positivity of A(2) and C(0) do not appear clearly from the expressions above. In fact, their positivity is not guaranteed as we will se in the example developed in Section 4. 3.4.3. Minimizing with respect to 〈v〉 Minimizing J (2) leads to the following variational equation for 〈u1〉∫ Ω ( A(0)εx (〈 u1 〉) + A(1)εx ( u0 ) + B(0)gx ( u0 )) · εx(〈v〉)dx = 0, ∀〈v〉 ∈ H10 (Ω,R2) (54) which admits a unique solution under the conditions stated in Remark 3. Accordingly, 〈u1〉 is solution of the following linear elastic problem defined on the “homogenized” body which involves the three effective tensor fields A(0), A(1) and B(0)divx ( A(0)εx (〈 u1 〉) + A(1)εx ( u0 ) + B(0)gx ( u0 )) = 0, in Ω〈 u1 〉 = 0, on ∂Ω (55) Once this problem is solved, u1 is entirely determined and we obtain from (47) that the third term u2 of the expansion of ue can read as u2 = εxpq (〈 u1 〉) χpq + εxpq ( u0 ) ψpq + gxpqr ( u0 ) ξpqr + 〈 u2 〉 , (56) its average value 〈u2〉 only remaining unknown at this stage. Moreover, the value of the energy J(2)∗ = J (2) (〈 u1 〉) is obtained from (50) J(2)∗ = ∫ Ω ( 1 2 A(0)εx (〈 u1 〉) · εx (〈 u1 〉) + ( A(1)εx ( u0 ) + B(0)gx ( u0 )) · εx (〈u1〉))dx + ∫ Ω ( 1 2 A(2)εx ( u0 ) · εx (u0)+ B(1)gx (u0) · εx (u0)+ 12C(0)gx (u0) · gx (u0) ) dx, (57) and the second order stress field σ1 is also known, σ1 = A ( εx ( u1 ) + εy ( u2 )) . (58) 3.5. To summarize At this stage, we have obtained the first two terms of the expansion both for the displacements and the stresses, that is to say u0,u1,σ0 and σ1. The next term u2 of the displacement expansion is known up to its average value 〈u2〉. Moreover the first three terms J(0)∗ , J (1) ∗ and J (2) ∗ of the expansion of the energy at equilibrium Je(ue) are also Second order homogenization of quasi-periodic structures 339 known and we could prove that the next term J(3)∗ is also determined. All these determi- nations require: (i) first, to solve three families of elastic problems posed on the current unit cell {x} × Y whose solutions χij,ψij and ξijk (with i, j and k running in {1, 2, 3}) de- pend only on the microstructure at x; (ii) then, to calculate the six effective tensor fields A(0),A(1),A(2),B(0),B(1),C(0) defined on the “homogenized” domain Ω; (iii) then, to solve the two macroscopic elastic problems whose solutions are u0 and 〈u1〉; (iv) finally, the fields u0,u1,σ0 and σ1 are obtained by linear combinations of the previous ones and the energetic terms J(0)∗ , J (1) ∗ and J (2) ∗ are given by integral calculus. We could continue the procedure to obtain the next terms. 4. AN ILLUSTRATIVE EXAMPLE To illustrate the procedure described in the previous sections, we consider the case of an heterogeneous cylinder made of a bi-layer laminate and submitted to the gravity. Specifically, the cylinder is Ω = S × (0, H) with height H and whose cross section S is a bounded, connected open subset of R2 with boundary ∂S. This cylinder is made of a finely bi-layered composite, the two materials constituting the layers having the same density ρ and being isotropically elastic with Lame´ coefficients (λ1, µ1) and (λ2, µ2) . The spatial distribution of the layers is such that the stiffness Ae depends on (x1, x2) only and not on x3, Ae(x)ijkl = λe (x1, x2) δijδkl + µe (x1, x2) ( δikδjl + δilδjk ) . (59) The lateral boundary of the cylinder ∂S× (0, H) is clamped whereas the sections S× {0} and S× {H} cannot move in the transversal directions but are free to move in the axial direction. Specifically the boundary conditions read as ue = 0 on ∂S× (0, H), { ue1 = u e 2 = 0 σe33 = 0 on S× {0, H}. This cylinder being submitted to the gravity −ge3, the body forces are f(x) = −ρge3 and hence are uniform. Under these conditions, it is easy to show that the (exact) displace- ment at equilibrium is anti-plane: ue(x) = ue3 (x1, x2) e3. Moreover, by linearity, its unique non-null component ue3 can read as ue3 (x1, x2) = ρgu e (x1, x2) , where ue is the unique solution of the following problem posed on the cross-section S{ div (µe∇ue) = 1, in S ue = 0, on ∂S (60) From an energetic point of view, ue is solution of the following minimization problem Variational formulation: Find ue ∈ H10(S) which minimizes the energy functional Je over 340 Duc Trung Le, Jean-Jacques Marigo H10(S) Je (ue) = min v∈H10 (S) Je(v) with Je(v) = ∫ S ( 1 2 µe∇v · ∇v− v ) dS. (61) Note that the solution depends only on the distribution of the shear modulus µe in the cross-section S. Concerning that distribution, one considers a quasi-periodic case charac- terized by the following assumption µe (x1, x2) = µ ( x1, x2, x1 e ) with µ (x1, x2, y) = { µ1 if 2|y| < θ (x1, x2) µ2 if θ (x1, x2) < 2|y| < 1 (62) where θ is a smooth function of (x1, x2) with values in [0, 1] and µ is periodic with re- spect to y with period 1. In other words, the two materials are layered in the direction 1, in proportion θ for the material 1 and 1− θ for the material 2, that proportion changing smoothly in the cross-section. So, we are in a situation corresponding to what is studied in the previous sections and even a little simpler. Indeed, here the body forces are con- stant, the unit cell is one-dimensional and corresponds to the interval Y = (−1/2,+1/2), the quasi-periodic repartition of the heterogeneities is characterized by the scalar func- tion θ, the unknown is a scalar field and the domain S is two-dimensional. In order to prevent boundary layer effects, it suffices to assume that θ is equal to 0 (or equivalently to 1) in a neighborhood of ∂S. We can follow the procedure described above to find the first terms of the expansion of the solution with respect to e. Throughout this section, the Greek indices α, β,γ, . . . run in {1, 2}. 4.1. Determination of χα3, a0α3β3,A (0) α3β3,A (1) α3β3 and B (0) α3β3γ Here all the displacement fields have only their third component which is not iden- tically zero and hence can be considered as scalar functions defined on S× (−1/2, 1/2) χij(x,y) = χij (x1, x2, y) e3. Furthermore, χα3 = χ3α only are to be determined for α = 1, 2. Starting from (20) and using the one-dimensional character of the unit cell, χ13 and χ23 are Y-periodic and such that  ∂y ( µ ( ∂yχ 13 + 1 )) = 0〈 χ13 〉 = 0 , { ∂y ( µ ( ∂yχ 23)) = 0〈 χ23 〉 = 0 (63) One easily deduces that χ23 = 0 and that χ13 is the following piecewise linear and odd function of y with slopes depending on x := (x1, x2) , see Fig. 2 χ13(x, y) = (µ2 − µ1) m(θ(x))  (1− θ(x))y when 2|y| ≤ θ(x) θ(x) ( sign(y) 2 − y ) when θ(x) ≤ 2|y| ≤ 1 (64) where m(θ) := (1− θ)µ1 + θµ2. (65) Note that χ13 = 0 in the homogeneous cases, i.e. when µ1 = µ2 or θ ∈ {0, 1}. Second order homogenization of quasi-periodic structures 341Hom genizati n of quasi-periodic struct res 17 12-12 Ȑ2-Ȑ2 -0.1 0.1 y Χ13 Fig. 2. Graph of χ13 when θ = 0.4 and µ2 = 3µ1 We will use also the derivatives ∂xαχ 13 of χ13 which read as ∂xαχ 13(x, y) = (µ2 − µ1) m ( θ(x) )2 ∂xαθ(x)  −µ2y when 2|y| < θ(x) µ1 ( sign(y) 2 − y ) when θ(x) < 2|y| < 1 . (66) Note that they are discontinuous at y = ±θ(x)/2 and that they are proportional to the gradient of θ. Note also that they are identically null when µ1 = µ2, but not when θ(x) = 0 or 1 and ∂xαθ(x) 6= 0. Using (30), one obtains the components a0α3β3 (and those which are deduced by symmetry) of the stress localization tensor field: a01313(x, y) = µ1µ2 m(θ(x)) , a01323(x, y) = a 0 2313(x, y) = 0, a 0 2323(x, y) = µ(x, y). (67) Then, their average value over Y gives the components of the effective stiffness tensor A(0): A (0) 1313(x) = µ1µ2 m(θ(x)) , A (0) 1323(x = A (0) 2313(x) = 0, A (0) 2323(x) = θ(x)µ1 + ( 1− θ(x))µ2 (68) which bring the harmonic mean value 1/〈1/µ〉 and the arithmetic mean value 〈µ〉 of the shear modulus into play. Since χ13 is an odd function of y and since χ23 vanishes, it is easy to deduce from (37) and (67) that all the relevant components of B(0) vanish: B (0) α3β3γ = B (0) 3αβ3γ = B (0) α33βγ = B (0) 3α3βγ = 0, ∀α, β, γ ∈ {1, 2}. (69) The relevant components of A(1) are A (1) α3β3 (and those deduced by symmetry). From (37) they read A (1) α3β3 = −〈∂xγ ( a0γ3β3 ) χα3〉. Since χ23 = 0, one immediately obtains A(1)2313 = A (1) 2323 = 0. Using (67) and noting that a 0 1313 does not depend on y, A (1) 1313 becomes A (1) 1313 = −〈∂x1 ( a01313 ) χ13〉 = −∂x1 ( a01313 ) 〈χ13〉 = 0. Fig. 2. Graph f 13 .4 and µ2 = 3µ1 e will use also the derivatives ∂xαχ 13 of χ13 which read as ∂xαχ 13(x, y) = (µ2 − µ1) m(θ(x))2 ∂xαθ(x)  µ2y hen 2|y| θ(x) µ1 ( sign(y) 2 − y ) when θ(x) < 2|y| < 1 (66) Note that they are discontinuous at y = ±θ(x)/2 and that they are proportional to the gradient of θ. Note also that they are identically null when µ1 = µ2, but not when θ(x) = 0 or 1 and ∂xαθ(x) 6= 0. Using (30), one obtains the components a0α3β3 (and those which are deduced by sym- metry) of the stress localization tensor field: a01313(x, y) = µ1µ2 m(θ(x)) , a01323(x, y) = a 0 2313(x, y) = 0, a 0 2323(x, y) = µ(x, y). (67) Then, their average value over Y gives the components of the effective stiffness tensor A(0) A(0)1313(x) = µ1µ2 m(θ(x)) , A(0)1323(x) = A (0) 2313(x) = 0, A (0) 2323(x) = θ(x)µ1 + (1− θ(x))µ2, (68) which bring the harmonic mean value 1/〈1/µ〉 and the arithmetic mean value 〈µ〉 of the shear modulus into play. Since χ13 is an odd function of y and since χ23 vanishes, it is easy to deduce from (37) and (67) that all the relevant components of B(0) vanish B (0) α3β3γ = B (0) 3αβ3γ = B (0) α33βγ = B (0) 3α3βγ = 0, ∀α, β,γ ∈ {1, 2}. (69) The relevant components of A(1) are A(1)α3β3 (and those deduced by symmetry). From (37) they read A(1)α3β3 = −〈∂xγ(a0γ3β3)χα3〉. Since χ23 = 0, one immediately obtains A(1)2313 = 342 Duc Trung Le, Jean-Jacques Marigo A (1) 2323 = 0. Using (67) and noting that a 0 1313 does not depend on y,A (1) 1313 becomes A (1) 1313 = − 〈 ∂x1 ( a01313 ) χ13 〉 = −∂x1 ( a01313 ) 〈 χ13 〉 = 0. The calculation of A(1)1323 requires some attention. Indeed, by (37) and (67), A (1) 1323 = −〈χ13∂x2µ〉 which brings the spatial derivative of µ into play. But since µ is piecewise constant, that derivative must be understood in a weak sense. Since 〈µχ13〉 = 0 by symmetry, one must read in fact A(1)1323 = 〈µ∂x2χ13〉 where the derivative of χ13 can be considered in the classical sense because χ13 is continuous and piecewise differentiable in S×Y. Accordingly, since ∂x2χ13 is an odd function of y and µ is an even function of y, one finally gets A(1)1323 = 0. Therefore, all the relevant components of A (1) in the present anti-plane context vanish. 4.2. Determination of u0, σ0α3, 〈u1〉, J(0)∗ , J(1)∗ The first term u0 of the expansion of ue depends on (x1, x2) only and is the unique solution of the following linear problem posed on the “homogenized” cross-section{ ∂x1 ( 〈1/µ〉−1∂x1 u0 ) + ∂x2 (〈µ〉∂x2 u0) = 1 in S u0 = 0 on ∂S (70) Therefore u0 depends only on the repartition of the proportion (x1, x2) 7→ θ(x1, x2) of the material 1 in the section S. Moreover u0 minimizes J(0) over H10(S) with J (0) given by J(0)(v) = 1 2 ∫ S ( 1 〈1/µ〉∂x1 v∂x1 v + 〈µ〉∂x2 v∂x2 v ) dS− ∫ S vdS. (71) By standard arguments, on deduces that the minimum is given by J(0)∗ := J (0) (u0) = −1 2 ∫ S u0dS. (72) Once u0 is determined, the first order stress field σ0 is obtained from (29) and (67) σ013(x) = 1 〈1/µ〉(x)∂x1 u 0(x), σ023(x, y) = µ(x, y)∂x2 u 0(x). (73) Since all the relevant components of A(1) and B(0) vanish, one immediately deduces from (36) and (55) that 〈 u1 〉 = 0, J(1)∗ = 0. (74) Therefore u1 is given by u1(x, y) = ∂x1 u 0(x)χ13(x, y). (75) 4.3. Determination of ψα3, ξα3β,A(2)α3β3,B (1) α3β3γ,C (0) α3βγ3ζ and J (2) ∗ The goal of this subsection is to obtain J(2)∗ , i.e. the term of the order of e2 in the energy expansion. For that, one must solve the other elementary problems and calculate the other effective stiffness tensors. Second order homogenization of quasi-periodic structures 343 4.3.1. Determination of ψα3 and ξα3β Let us first note that the effective stiffness tensors A(2),B(1) and C(0) bring into play the gradients with respect to y of the fields ψα3, ξα3β only, see (51)–(53). Therefore, it is suf- ficient to determine ∂yψα3 and ∂yξα3β, leaving undetermined the constant of integration (which could be obtained by the conditions 〈ψα3〉 = 〈ξα3β〉 = 0). From (48) and (67), one sees that ψ13 must satisfy〈 µ ( ∂yψ 13 + ∂x1χ 13 ) ∂yφ 〉 = 〈 ∂x1a 0 1313φ 〉 , ∀φ ∈ H1#(S), 〈φ〉 = 0. But since a01313 does not depend on y and 〈φ〉 = 0, the right hand side of the above variational equality vanishes and hence 〈µ(∂yψ13 + ∂x1χ13)∂yφ〉 = 0 for all φ ∈ H1#(S) such that 〈φ〉 = 0. Therefore µ(∂yψ13 + ∂x1χ13) does not depend on y and one gets ∂yψ 13(x, y) + ∂x1χ 13(x, y) = a(x) µ(x, y) , (76) with the constant a(x) to be determined. But since χ13 is Y-periodic and since χ13 is an odd function of y, one has 〈∂yψ13〉 = 〈∂x1χ13〉 = 0. Hence (76) gives a(x) = 0 and one finally obtains ∂yψ 13 = −∂x1χ13, 〈 ψ13 〉 = 0, (77) which determines ψ13 with the help of (64). The determination of ψ23 requires some attention. Indeed, since χ23 = 0, using (48) and (67), ψ23 should satisfy〈 µ∂yψ 23∂yφ 〉 = 〈φ∂x2µ〉 , ∀φ ∈ H1#(S), 〈φ〉 = 0. But since µ is piecewise constant, the derivative of µ must be understood in a weak sense. Specifically, let us consider that 〈φ∂x2µ〉 is a distribution, i.e. an element of the dual D′(S) of D(S). Then 〈φ∂x2µ〉 is defined by 〈φ∂x2µ〉 (ϕ) = − ∫ S×Y µ(x, y)φ(y)∂x2ϕ(x)dydS, ∀ϕ ∈ D(S). By Fubini theorem, one gets 〈φ∂x2µ〉 (ϕ) = − ∫ S ( µ1 ∫ 2|y|<θ(x) φ(y)dy + µ2 ∫ θ(x)<2|y|<1 φ(y)dy ) ∂x2ϕ(x)dS. Now, owing to the assumed regularity of x 7→ θ(x) , one can make the integration by parts with respect to x2 to obtain 〈φ∂x2µ〉 (ϕ) = − ∫ S 1 2 (µ2 − µ1) ∂x2θ(x)(φ(−θ(x)/2) + φ(+θ(x)/2))ϕ(x)dS. Therefore ψ23 must in fact satisfy the following variational equation (where the depen- dence on x is not indicated)∫ 1/2 −1/2 µ∂yψ 23∂yφdy = −12 (µ2 − µ1) ∂x2θ(φ(−θ/2) + φ(+θ/2)), (78) 344 Duc Trung Le, Jean-Jacques Marigo for all φ ∈ H1#(S) with 〈φ〉 = 0. It can be solved in a closed form and after some calcula- tions left to the reader one finally gets ∂yψ 23(x, y) = (µ1 − µ2) ∂x2θ(x)  y µ1 when 2|y| < θ(x) y µ2 − sign(y) 2µ2 when θ(x) < 2|y| < 1 (79) That determines ψ23 up to a constant which is fixed by the condition 〈ψ23〉 = 0. From (49) and (67), one gets that ξ131 is the unique solution of ξ 131 ∈ H1#(S), 〈 ξ131 〉 = 0〈 µ ( ∂yξ 131 + χ13 ) ∂yφ 〉 = 0, ∀φ ∈ H1#(S), 〈φ〉 = 0. Let us verify that ξ131 given by ∂yξ 131 = −χ13, 〈 ξ131 〉 = 0, (80) is solution. Indeed, ξ131 given by (80) verifies the variational equation above. From 〈χ13〉 = 0, one deduces that 〈∂yξ131〉 = 0 and hence that ξ131 is Y-periodic. The relation ∂yξ 131 = −χ13 determines ξ131 up to a constant which is fixed by the condition 〈ξ131〉 = 0. In a same manner, one gets from (49) and (67) that ξ132 = ξ231 = 0. It remains to find ξ232 which is the unique solution of{ ξ232 ∈ H1#(S), 〈 ξ232 〉 = 0,〈 µ∂yξ 232∂yφ 〉 = 〈µφ〉, ∀φ ∈ H1#(S), 〈φ〉 = 0. It can be solved in a closed form and after some calculations left to the reader one finally gets ∂yξ 232(x, y) = (µ2 − µ1)  1− θ µ1 y when 2|y| < θ(x) θ µ2 ( sign(y) 2 − y ) when θ(x) < 2|y| < 1 (81) That determines ξ232 up to a constant which is fixed by the condition 〈ξ232〉 = 0. Let us note that the displacement fields ψ13,ψ23, ξ131 and ξ232 enjoy the following properties: (1) They are identically null when the medium is homogeneous, i.e. when θ ∈ {0, 1} everywhere or when µ1 = µ2, but ψ13(x) or ψ23(x) are not identically null when the medium is only locally homogenous at x with θ(x) ∈ {0, 1} and (∂x1θ(x), ∂x2θ(x)) 6= (0, 0). (2) They are even, periodic and piecewise quadratic functions of y, see Fig. 3. 4.3.2. Calculation of A(2)α3β3,B (1) α3β3γ and C (0) α3βγ3ζ Using (51)–(53), the relevant components of these effective stiffness tensors read as A (2) α3β3 = 〈 µ∂xζχ α3∂xζχ β3 − µ∂yψα3∂yψβ3 〉 , (82) Second order homogenization of quasi-periodic structures 345 Homogenization of quasi-periodic structures 21 12-12 Ȑ2-Ȑ2 -0.02 0.02 y Ψ13 12-12 Ȑ2-Ȑ2 -0.02 0.02 y Ψ23 12-12 Ȑ2-Ȑ2 -0.02 0.02 y Ξ131 12-12 Ȑ2-Ȑ2 -0.02 0.02 y Ξ232 Fig. 3. Graph of ψ13, ψ23, ξ131 and ξ232 when θ = 0.4, µ2 = 3µ1 and ∂x1θ = ∂x2θ = 1 4.3.2. Calculation of A (2) α3β3, B (1) α3β3γ and C (0) α3βγ3ζ Using (51)–(53), the relevant components of these effective stiffness tensors read as A (2) α3β3 = 〈µ∂xζχα3∂xζχβ3 − µ∂yψα3∂yψβ3〉 (82) B (1) α3β3γ = 〈µ∂xγχα3χβ3 − µ∂yψα3∂yξβ3γ〉 (83) C (0) α3βγ3ζ = 〈µχα3χγ3δβζ − µ∂yξα3β∂yξγ3ζ〉 (84) where δ denotes the Kronecker symbol. Since all the elementary fields are known, it suffices to calculate the integrals over Y . Using (66), (77), (79) and (82), one gets for A(2): A (2) 1313 = 〈µ∂x2χ13∂x2χ13〉 = (µ2 − µ1)2µ1µ2 ( (1− θ)3µ1 + θ3µ2 ) 12 ( (1− θ)µ1 + θµ2 )4 (∂x2θ)2, (85) A (2) 1323 = 〈µ∂x1χ13∂yψ23〉 = (µ2 − µ1)2 ( (1− θ)3µ1 + θ3µ2 ) 12 ( (1− θ)µ1 + θµ2 )2 ∂x1θ ∂x2θ, (86) A (2) 2323 = −〈µ∂yψ23∂yψ23〉 = − (µ2 − µ1)2 12µ1µ2 ( (1− θ)3µ1 + θ3µ2 )( ∂x2θ )2 . (87) For B(1), using (83) and the expressions for the elementary fields, one gets B (1) 13131 = B (1) 13231 = B (1) 23132 = 0, (88) Fig. 3. Graph of ψ13,ψ23, ξ131 and ξ232 when θ = 0.4, µ2 = 3µ1 and ∂x1θ = ∂x2θ = 1 B (1) α3β3γ = 〈 µ∂xγχ α3χβ3 − µ∂yψα3∂yξβ3γ 〉 , (83) C (0) α3βγ ζ = 〈 µχα3χγ3δβζ ∂yξ α3β∂yξ γ3ζ 〉 , (84) where δ denotes the Kronecker symbol. Since all the elementary fields are known, it suffices to calculate the integrals over Y. Using (66), (77), (79) and (82), one gets for A(2) A (2) 1313 = 〈 µ∂x2χ 13∂x2χ 13 〉 = (µ2 − µ1 µ1µ2 ( (1− θ)3µ1 + θ3µ2 ) 12 ((1− θ)µ1 + θµ2)4 (∂x2θ) 2 , (85) A (2) 1323 = 〈 µ∂x1χ 13∂yψ 23 〉 = (µ2 − µ1)2 ( (1− )3µ1 + θ3µ2 ) 12 ((1− θ)µ1 + θµ2)2 ∂x1θ∂x2θ, (86) A (2) 2323 = − 〈 µ∂yψ 23∂yψ 23〉 = − (µ2 − µ1)2 12µ1µ2 ( (1− θ)3µ1 + θ3µ2 ) (∂x2θ) 2 . (87) For B(1), using (83) and the expressions for the elementary fields, one gets B (1) 13131 = B (1) 13231 = B (1) 23132 = 0, (88) B (1) 13132 = 〈 µχ13∂x2χ 13 〉 = (µ2 − µ1)2 µ1µ2θ(1− θ)(1− 2θ) 12 ((1− θ)µ1 + θµ2)3 ∂x2θ, (89) B (1) 23131 = 〈 µχ13∂yψ 23 〉 = (µ2 − µ1)2 θ(1− θ)(1− 2θ) 12 ((1− θ)µ1 + θµ2) ∂x2θ, (90) 346 Duc Trung Le, Jean-Jacques Marigo B (1) 13232 = 〈 µ∂x1χ 13∂yξ 232 〉 = (µ2 − µ1)2 θ(1− θ) ( (1− θ)2µ1 − θ2µ2 ) 12 ((1− θ)µ1 + θµ2)2 ∂x1θ, (91) B (1) 23232 = − 〈 µ∂yψ 23∂yξ 232〉 = − (µ2 − µ1)2 θ(1− θ) ((1− θ)2µ1 − θ2µ2) 12µ1µ2 ∂x2θ. (92) For C(0), using (84) and the expressions for the elementary fields, one gets that the unique non null components are C (0) 131232 = 〈 µχ13∂yξ 232 〉 = (µ2 − µ1)2 θ2(1− θ)2 12 ((1− θ)µ1 + θµ2) , (93) C (0) 132132 = 〈 µχ13χ13 〉 = (µ2 − µ1)2 θ2(1− θ)2 (θµ1 + (1− θ)µ2) 12 ((1− θ)µ1 + θµ2)2 , (94) C (0) 232232 = − 〈 µ∂yξ 232∂yξ 232〉 = − (µ2 − µ1)2 12µ1µ2 θ2(1− θ)2 ((1− θ)µ1 + θµ2) , (95) and those which are obtained by symmetry. 4.3.3. Calculation of J(2)∗ Inserting the expressions above into (57) gives the term of order e2 in the expansion of the energy, J(2)∗ = ∫ S 1 2 ( A (2) 1313∂1u 0∂1u0 + 2A (2) 1323∂1u 0∂2u0 + A (2) 2323∂2u 0∂2u0 ) dS + ∫ S ( B (1) 13132∂1u 0∂212u 0 + B (1) 13232∂1u 0∂222u 0 + B (1) 23131∂2u 0∂211u 0 + B (1) 23232∂2u 0∂222u 0 ) dS + ∫ S 1 2 ( C2132132∂ 2 12u 0∂212u 0 + 2C(0)131232∂ 2 11u 0∂222u 0 + C(0)232232∂ 2 22u 0∂222u 0 ) dS, where ∂αu0 and ∂2αβu 0 stand for ∂xαu 0 and ∂xα(∂xβu 0), respectively. Since the components of A(2) and B(1) contain the gradient of θ, it is more convenient to render this dependence explicit by expressing J(2)∗ as J(2)∗ = ∫ S 1 2 ( a11 ( ∂1θ∂1u0 )2 + 2a12∂1θ∂1u0∂2θ∂2u0 + a22 ( ∂2θ∂2u0 )2) dS + ∫ S ( b112∂2θ∂1u0∂212u 0 + b122∂1θ∂1u0∂222u 0 + b211∂2θ∂2u0∂211u 0 + b222∂2θ∂2u0∂222u 0)dS + ∫ S 1 2 ( c1212∂212u 0∂212u 0 + 2c1122∂211u 0∂222u 0 + c2222∂222u 0∂222u 0)dS. (96) In (96) the coefficients aαβ, bαβγ and cαβγζ depend only on θ and their expression can be easily obtained from (85)–(95). Remark 5. The second order energy term depends both on the second gradient of the displace- ment and on the gradient of the characteristic parameter of the microstructure. It turns out that the coefficients entering in the expression of this energy term have diffe rent signs: for instance, A (2) 1313,C (0) 131232 and C (0) 132132 are positive whereas A (2) 2323 and C (0) 232232 are negative. Second order homogenization of quasi-periodic structures 347 5. CONCLUSION In this paper, we developed a homogenization procedure based on the combination of asymptotic expansions and a variational approach in order to construct the higher or- der macroscopic energy of a quasi-periodic heterogeneous medium. The main assump- tion used in this approach to describe the quasi-periodicity is that the properties of a point in the cell depends not only on its position but also the position of the cell in the medium. By using the asymptotic expansion of the strain energy, the minimization problem of en- ergy microscopic becomes a series of successive minimization problems. The solutions of each problems give the corresponding macroscopic elastic energies, the components of these higher-order elastic energies are obtained by solving the cell problems where the solutions of previous orders become the entries for the next orders. We have shown that the effective energy density depends on the strain gradient and on the gradient of the microstructure as of the second order. Concerning the perspectives open by this work and what extensions could be inves- tigated, let us mention the following ones: - The method could be apply to obtain the effective thermo-mechanical properties of thermo-elastic quasi-periodic media. - It would be interesting to apply the same procedure for more general quasi-periodic media like those mentioned in Remark 1. - The method was applied here by assuming that the microstructure is given (and fixed). It would be interesting to remove this assumption by considering that the mi- crostructure depends on parameters that one has to optimize. For instance, we could extend the study made in Section 4 by considering that the proportion θ of the material 1 (which can be considered as the damaged material if µ1 < µ2) is governed by a principle of least energy like in [17]. - Another possible application would be to consider micro-cracked media where the length and the orientation of the micro-cracks vary smoothly in the domain. We could so justify some models of damage which are regularized by introducing gradient damage terms like in [18–20]. - The analysis was made here by considering situations where there is no boundary layer effects. Since these situations are more the exception than the rule, it would be important to extend our study to the cases where boundary layer effects exist. REFERENCES [1] A. Bensoussan, J. L. Lions, and G. Papanicolaou. Asymptotic analysis of periodic structures, Vol. 374. Elsevier Science, (1978). [2] F. Murat and L. Tartar. Calcul des variations et homoge´ne´isation. Les me´thodes del’homoge´ne´isation: the´orie et applications en physique, 57, (1985), pp. 319–369. [3] G. Allaire. Homogenization and two-scale convergence. SIAM Journal on Mathematical Anal- ysis, 23, (6), (1992), pp. 1482–1518. https://doi.org/10.1137/0523084. [4] H. Dumontet. Study of a boundary layer problem in elastic composite materials. ESAIM: Mathematical Modelling and Numerical Analysis, 20, (2), (1986), pp. 265–286. https://doi.org/10.1051/m2an/1986200202651. 348 Duc Trung Le, Jean-Jacques Marigo [5] G. A. Francfort and F. Murat. Homogenization and optimal bounds in linear elas- ticity. Archive for Rational Mechanics and Analysis, 94, (4), (1986), pp. 307–334. https://doi.org/10.1007/bf00280908. [6] R. Abdelmoula and J. J. Marigo. The effective behavior of a fiber bridged crack. Journal of the Mechanics and Physics of Solids, 48, (11), (2000), pp. 2419–2444. https://doi.org/10.1016/s0022- 5096(00)00003-x. [7] V. A. Marchenko and E. Y. Khruslov. Homogenization of partial differential equations, Vol. 46. Springer Science & Business Media, (2008). [8] U. Hornung. Homogenization and porous media, Vol. 6. Interdisciplinary Applied Mathematics, (1997). [9] B. Gambin and E. Kro¨ner. Higher-order terms in the homogenized stress-strain rela- tion of periodic elastic media. Physica Status Solidi (B), 151, (2), (1989), pp. 513–519. https://doi.org/10.1002/pssb.2221510211. [10] C. Boutin. Microstructural effects in elastic composites. International Journal of Solids and Structures, 33, (7), (1996), pp. 1023–105. https://doi.org/10.1016/0020-7683(95)00089-5. [11] I. V. Andrianov, J. Awrejcewicz, and A. A. Diskovsky. Homogenization of quasi- periodic structures. Journal of Vibration and Acoustics, 128, (4), (2006), pp. 532–534. https://doi.org/10.1115/1.2202158. [12] F. Su, Z. Xu, J. Z. Cui, and Q. L. Dong. Multi-scale method for the quasi-periodic structures of composite materials. Applied Mathematics and Computation, 217, (12), (2011), pp. 5847–5852. https://doi.org/10.1016/j.amc.2010.12.068. [13] L. Guillot, Y. Capdeville, and J. J. Marigo. 2-D non-periodic homogenization of the elas- tic wave equation: SH case. Geophysical Journal International, 182, (3), (2010), pp. 1438–1454. https://doi.org/10.1111/j.1365-246x.2010.04688.x. [14] Y. Capdeville, L. Guillot, and J. J. Marigo. 2-D non-periodic homogenization to upscale elastic media for P-SV waves. Geophysical Journal International, 182, (2), (2010), pp. 903–922. https://doi.org/10.1111/j.1365-246x.2010.04636.x. [15] E. S. Palencia. Non-homogeneous media and vibration theory, Vol. 127. Springer-Verlag Berlin, (1980). [16] F. Devries, H. Dumontet, G. Duvaut, and F. Le´ne´. Homogenization and damage for com- posite structures. International Journal for Numerical Methods in Engineering, 27, (2), (1989), pp. 285–298. https://doi.org/10.1002/nme.1620270206. [17] G. A. Francfort and J. J. Marigo. Stable damage evolution in a brittle continuous medium. European Journal of Mechanics Series a Solids, 12, (1993), pp. 149–149. [18] K. Pham and J. J. Marigo. Approche variationnelle de l’endommagement: I. Les concepts fondamentaux. Comptes Rendus Me´canique, 338, (4), (2010), pp. 191–198. https://doi.org/10.1016/j.crme.2010.03.009. [19] K. Pham and J. J. Marigo. Approche variationnelle de l’endommagement: II. Les mode`les a` gradient. Comptes Rendus Me´canique, 338, (4), (2010), pp. 199–206. [20] C. Dascalu, G. Bilbie, and E. Agiasofitou. Damage and size effects in elastic solids: a homog- enization approach. International Journal of Solids and Structures, 45, (2), (2008), pp. 409–430. https://doi.org/10.1016/j.ijsolstr.2007.08.025.