Localization of elastic deformation in strongly anisotropic, porous, linear materials with periodic microstructures: Exact solutions and dilute expansions

ARTICLE IN PRESS

Journal of the Mechanics and Physics of Solids 56 (2008) 1245–1268 www.elsevier.com/locate/jmps

Localization of elastic deformation in strongly anisotropic, porous, linear materials with periodic microstructures: Exact solutions and dilute expansions Franc- ois Willota,b,c, Yves-Patrick Pellegrinia,, Pedro Ponte Castan˜edab,c De´partement de Physique The´orique et Applique´e, Commissariat a` l’E´nergie Atomique, BP12, 91680 Bruye`res-le-Chaˆtel, France b De´partement de Me´canique, E´cole Polytechnique, 91128 Palaiseau Cedex, France c Department of Mechanical Engineering and Applied Mechanics, University of Pennsylvania, Philadelphia, PA 19104-6315, USA a

Received 2 May 2007; received in revised form 2 October 2007; accepted 7 October 2007

Abstract Exact solutions are derived for the problem of a two-dimensional, infinitely anisotropic, linear-elastic medium containing a periodic lattice of voids. The matrix material possesses either one infinitely soft, or one infinitely hard loading direction, which induces localized (singular) field configurations. The effective elastic moduli are computed as functions of the porosity in each case. Their dilute expansions feature half-integer powers of the porosity, which can be correlated to the localized field patterns. Statistical characterizations of the fields, such as their first moments and their histograms are provided, with particular emphasis on the singularities of the latter. The behavior of the system near the void close-packing fraction is also investigated. The results of this work shed light on corresponding results for strongly non-linear porous media, which have been obtained recently by means of the ‘‘second-order’’ homogenization method, and where the dilute estimates also exhibit fractional powers of the porosity. r 2007 Elsevier Ltd. All rights reserved. Keywords: Voids and inclusions; Anisotropic material; Constitutive behaviour; Energy methods; Localization

1. Introduction 1.1. Context Most of the currently available homogenization methods for strongly non-linear composites make use of an underlying linear homogenization estimate, whether aimed at computing dielectric (Willis, 1986; Zeng et al., 1988; Ponte Castan˜eda et al., 1992), or elastic–plastic transport properties (Ponte Castan˜eda, 1991, 1996, 2002; Masson et al., 2000; Lahellec and Suquet, 2004). The best results are obtained using an anisotropic linear estimate, the anisotropy of which is consistently determined, by the means and variances of the fields in each Corresponding author. Tel.: +33 1 69 26 73 49; fax: +33 1 69 26 70 77.

E-mail addresses: [email protected] (F. Willot), [email protected] (Y.-P. Pellegrini), [email protected] (P.P. Castan˜eda). 0022-5096/$ - see front matter r 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.jmps.2007.10.002

ARTICLE IN PRESS 1246

F. Willot et al. / J. Mech. Phys. Solids 56 (2008) 1245–1268

phase of the composite (Pellegrini, 2001; Ponte Castan˜eda, 2001). The anisotropy of the underlying linear medium accounts for the privileged direction imposed by the driving field in the non-linear medium. Applying the ‘‘second-order’’ method (Ponte Castan˜eda, 2002) to a random power-law material weakened by aligned cylindrical voids, it has been found that in the so-called ‘‘dilute’’ limit (i.e., that of a vanishingly small volume fraction of porosity f), and in the limit of infinite exponent where the matrix becomes ideally plastic with a yield threshold, the leading correction to the yield stress in pure shear behaves as f 2=3 . This result holds for both Hashin–Shtrikman and self-consistent estimates. The second-order method of Pellegrini (2001), suitably modified by replacing a Gaussian ansatz for the field distributions by a Heaviside ansatz (so as to cutoff high fields in a way compatible with threshold-type materials), provides in some cases similar predictions (Pellegrini and Ponte Castan˜eda, 2001, unpublished). This f 2=3 dilute behavior is unusual in the context of effective-medium theories for random media, where the low-order terms in dilute expansions from the literature usually consist in integer powers of f. With the f 2=3 correction, the derivative of the yield stress with respect to the porosity is infinite at f ¼ 0, due to the exponent being lower than 1. Thus, a vanishingly small volume fraction of voids induces a dramatic weakening of the porous medium. Ponte Castan˜eda (1996, 2002) interprets this phenomenon as an indication of localizing behavior in a regime where shear bands pass through the pores, remarking that the limit analysis of Drucker (1966) produces an upper bound for the yield threshold with a dilute correction behaving as f 1=2 in 2D. Numerical calculations and tests on perforated plates with periodically distributed holes (Francescato and Pastor, 1998) are consistent with these predictions in a plane stress situation. However, for a square network of circular holes, kinematic and static limit analyses result in the following analytical bounds on the effective yield stress Ye (Francescato et al., 2004): pffiffiffi 1  2ðf =pÞ1=2 pYe =Y pð2= 3Þ½1  2ðf =pÞ1=2 , (1) where Y is the yield threshold of the matrix. Numerical limit analysis on a hollow disk model in plane strain leads to somewhat different conclusions, with a correction behaving roughly as f 2=3 for uniform strain boundary conditions, and as f 1=2 for uniform stress boundary conditions (Pastor and Ponte Castan˜eda, 2002). However, the calculations were not sufficiently accurate to be completely definitive. On the other hand, the exponent is quite certainly associated with localizing behavior, and its presence seems to be independent on whether the system is random, periodic, or comprises a unique void (although the actual values of the exponent may be dependent on the specific microstructure considered). Because the second-order theory with exponent 23 inherits its properties from an underlying anisotropic linear effective-medium theory, whose anisotropy is consistently determined by the field variances, a natural question concerns the ‘‘localizing’’ properties of the anisotropic linear theory itself, and its relation with possible non-analyticities of the effective shear moduli in the dilute limit. Bearing these considerations in mind, the present paper focuses on a new exact solution for the special, but representative, case of a two-dimensional (2D) array of voids embedded in an anisotropic linear elastic matrix, in the singular limit of infinite anisotropy. We emphasize that no exact solution of a similar type is available for general non-linear periodic media, which provides additional motivation for undertaking this study. 1.2. Organization of the paper The constitutive laws, the prescribed loading conditions and the notations are defined in Section 2. The word ‘‘loading’’ refers to the prescription of either overall conditions of homogeneous stress, or of homogeneous strain in the linear medium. Both are equivalent since they are related by the effective moduli. The matrix in the composite is compressible and possesses anisotropic properties in shear. Two special limits of infinite anisotropy which lead to exact results are then considered (Sections 3 and 4). For each of these limits, simple shear (SS), pure shear (PS), and equibiaxial loading situations are examined, and solutions are provided in each case (for some loading and anisotropy conditions, however, only the limit of an incompressible matrix is considered; but this limit is the one relevant to non-linear homogenization of plastic porous media). The effective moduli and their dilute expansions, together with the mean and variance of the strain and stress fields are computed in each case. All these quantities are relevant to non-linear homogenization theories. A useful way of condensing the information contained in the solutions is by using field distributions which can in principle directly be used to compute the effective energy (Pellegrini, 2001) or, more trivially, the moments of the fields. Local extrema or

ARTICLE IN PRESS F. Willot et al. / J. Mech. Phys. Solids 56 (2008) 1245–1268

1247

saddle points in the field maps induce singularities in the distributions such as power or logarithmic divergences, or discontinuities, of the types observed in the density of states of a crystal (Van Hove, 1953; Abrikosov et al., 1993). These singularities are examined in the context of our exact results. Cule and Torquato (1998) computed analytically the distributions of the electric field in a particular dielectric/conductor composite. They found the singularities of the Van Hove type, but no extended singularities of the type put forward by Abrikosov et al. We extract below some singularities of this type, while singularities of yet a different type are encountered in Section 4.2.3. Obviously, in the periodic case, the field distributions are redundant with the exact solutions derived hereafter. However, for randomly disordered situations where exact solutions are not available, they remain the only way to collect useful informations on the fields. To ease their interpretation, a knowledge of their features in the periodic case is desirable, which provides a justification for the studies of Sections 3.2.4 and 4.2.3. Moreover, we find that the infinite field variances obtained for some particular loadings are correlated to field singularities directly linked to localization patterns, which have counterparts as characteristic features in the distributions. A summary of our findings and a conclusion close the paper (Section 5). 2. Problem formulation 2.1. Material constitutive law and microstructure We consider a periodic porous medium, with a square unit cell of size L made of a linear-elastic matrix (denoted by a phase index b ¼ 1), containing a single disk-shaped void of radius a (phase index b ¼ 2). Cartesian reference axes Ox, Oy are defined as in Fig. 1(a), such that the void center lies at ðx; yÞ ¼ ð0; 0Þ. The unit cell is the square region ½L=2; L=2  ½L=2; L=2. The porosity is the surface concentration of voids f ¼ pða=LÞ2 . The close-packing value where the voids touch is f c ¼ p=4 ’ 0:78, for a ¼ L=2. Hereafter, all lengths are rescaled by L, so that L  1. The following notations are employed throughout the text: the upperright quadrant (URQ) of the unit cell denotes the region ðx; yÞ 2 ½0; 122 , the lower-right quadrant (LRQ) stands for ðx; yÞ 2 ½0; 12  ½12; 0. Similar abbreviations ULQ and LLQ stand for the corresponding left quadrants. Small deformations and plane strain conditions are assumed, so that the strain e derives from the 2D (in-plane) displacement field u, and ij ¼ 0 if i or j ¼ z. In-plane stress equilibrium requires qi sij ¼ 0, i; j ¼ 1; 2 and we take the constitutive relation such that sxz ¼ syz ¼ 0, so that the problem is 2D. The linear constitutive relation in the medium is rðxÞ ¼ LðxÞ : eðxÞ. In the voids, LðxÞ  0. In the anisotropic matrix phase, LðxÞ has components (Latin indices henceforth vary over 1 and 2): SS PS Lð1Þ ij;kl ¼ 2kJ ij;kl þ 2lE ij;kl þ 2mE ij;kl ,

(2)

where k is the bulk compressibility modulus, where l, m are in-plane anisotropic shear moduli, and where the operators J, ESS , EPS are mutually orthogonal projectors of components (no summation over repeated indices here): J ij;kl ¼ 12dij dkl ;

1 E SS ij;kl ¼ 2ð1  dij Þð1  dkl Þ;

1 E PS ij;kl ¼ dij dkl ðdik  2Þ.

(3)

Fig. 1. Left, periodic porous medium with unit cell and reference axes (a). Right, unit cell with (b) 0-type fibers; and (c) 45-type fibers. The black arrows depict eigenmodes of strain: simple shear (SS) in (b), and pure shear (PS) in (c).

ARTICLE IN PRESS F. Willot et al. / J. Mech. Phys. Solids 56 (2008) 1245–1268

1248

A symmetric 2D tensor a is expanded as a ¼ am I þ aSS eSS þ aPS ePS where I is the 2  2 identity matrix and where     0 1 1 0 ; ePS ¼ , (4) eSS ¼ 1 0 0 1 the eigenvector sets of which are related by a 45 rotation. Corresponding strain modes are represented by arrows in Fig. 1(b) and (c). The equibiaxial component of a is am  ðaxx þ ayy Þ=2. We call the second and third components the SS and PS components, with aSS  axy and aPS  ðaxx  ayy Þ=2. Then J : a ¼ am I, ESS : a ¼ aSS eSS and EPS : a ¼ aPS ePS . Accordingly, the constitutive relations decompose into: sSS ¼ 2lSS ;

sPS ¼ 2mPS ;

sm ¼ 2km .

(5)

Elastic tensor (2) is of a special type of orthotropy in the plane with L1111 ¼ L2222 aL1122 aL1212 . When l4m, this type of anisotropic response could correspond to a fiber-reinforced material with fibers aligned with the cell axes (see Fig. 1(b)). On the other hand, when m4l, it would correspond to a material with reinforcing fibers aligned with the 45 direction (see Fig. 1(c)). We introduce a shear anisotropy ratio, a, and normalized shear elastic moduli, m and ‘: a ¼ l=m;

m ¼ m=k;

‘ ¼ l=k.

(6)

The principal directions of the matrix being ‘‘aligned’’ with dense lines of voids, a reinforcement of anisotropyinduced effects is expected. Besides, the loading modes considered below are aligned with the eigenstrains. These choices, motivated by future applications of this work to non-linear homogenization, focus on situations most relevant to plastic localization in porous media. Indeed, in the non-linear theory, the loading determines the anisotropy of the background linear medium and always coincides with one of its eigendirections. Moreover, shear bands in random porous media quite generally link neighboring voids together. 2.2. Overall behavior In this work, we are concerned with the problem of determining the overall behavior of the periodic, porous material described in the previous subsection. This overall behavior is defined as the relation between the average stress hri and the average strain hei. Under the assumption of separation of length scales, the overall behavior of the porous material may be determined from the effective strain potential (see Torquato, 2002) W ¼ ð1  f Þ minhwð1Þ ðeÞið1Þ , e2K

(7)

where hið1Þ denotes an average over the matrix phase, wð1Þ ðeÞ ¼ 12e : Lð1Þ : e is the strain potential in the matrix and K ¼ feje ¼ 12½ru þ ðruÞT , u ¼ heix þ u , u periodicg is the set of kinematically admissible strain fields. The overall constitutive relation is then given by hri ¼

qW . qhei

(8)

Because of linearity, we define the overall elasticity tensor e L of the porous material via the relation eij;kl defined by hri ¼ e L : hei. This tensor can be shown to take on the same form as (2) with components L e effective moduli l, e m and e k. In the analyses below, it will sometimes be more convenient to work with the effective stress potential U, such that the overall constitutive relation may be equivalently written hei ¼

qU ; qhri

U ¼ ð1  f Þ minhuð1Þ ðrÞið1Þ , r2S

(9)

where uð1Þ ðrÞ ¼ ð1=2Þr : ðLð1Þ Þ1 : r is the stress potential in the matrix and S denotes the set of periodic stresses that are divergence free in the unit cell, with prescribed average hri, and traction free on the boundaries of the pores. It should be noted that the case of an isotropic matrix (a ¼ 1) has been addressed by McPhedran and Movchan (1994) in the more general framework of arbitrary contrast between the matrix and the square array of isotropic inclusions.

ARTICLE IN PRESS F. Willot et al. / J. Mech. Phys. Solids 56 (2008) 1245–1268

1249

2.3. Limits of infinite anisotropy and loading modes Hereafter, only the strong anisotropy limits a ¼ 0 and a ¼ þ1, amenable to an exact solution, are considered. Therefore: (i) when a ¼ 0 (i.e. l ¼ 0 or m ¼ 1) the medium is soft for SS loading, and resists PS loading; (ii) when a ¼ 1 (i.e. l ¼ 1 or m ¼ 0) the medium is soft for PS loading, and resists SS, loading. Equibiaxial, SS, and PS loading modes will be considered separately hereafter, so that only one of the three averaged components of the strain or stress is non-zero at a time. It is denoted by  (resp. s). Depending on the loading mode considered, strain and stress components and displacements enjoy various symmetry properties summarized in Appendix A. For clarity, we detail the correspondence between the loading modes and the anisotropy properties of the material. Consider, e.g., the limiting case a ¼ l=m ! 0, attained by two different types of materials: (i) l ! 0 and m40 (Material 1); (ii) m ! 1 and lo1 (Material 2). We seek solutions with strain and stress finite almost everywhere (a.e.), i.e., except on points or on lines in the plane. In particular, sPS ¼ 2mPS is finite a.e., so that PS  0 a.e. in Material 2. Likewise, SS ¼ sSS =ð2lÞ is finite a.e., so that sSS  0 a.e. in Material 1. Hence, Material 1 is infinitely soft in the SS direction and can be finitely loaded in the PS direction, whereas Material 2 is infinitely rigid in the PS direction and can be finitely loaded only in the SS direction. Similar considerations apply for a ¼ 1, leading to the following correspondence: a¼0: a¼1:

PS loading2l ¼ 0; SS loading2m ¼ 0;

SS loading2m ¼ 1, PS loading2l ¼ 1.

In equibiaxial loading, both types of materials should be considered for each of the limits a ¼ 0 or 1. As is well known, the elasticity tensor (2) is positive definite when l, m, and k are all strictly positive, and we have existence and uniqueness of solutions. As mentioned above, however, the interest in this work is for the limiting cases where one of the eigenvalues of the elasticity tensor (2) (or, of its inverse, the compliance tensor) tend to zero. For these cases, care must be exercised when interpreting the solutions. As will be seen below, the relevant potential energy, or complementary energy functionals can still be shown to be strictly convex in the subspace of allowable strains, or stresses, and the standard theorems (see, for example, Proposition 1.2 of Ekeland and Temam, 1974) would still ensure existence and uniqueness of the solutions. However, appropriate components of the strain (or stress) field can develop discontinuities in these limiting cases. In this connection, it is also relevant to note that the governing equations lose ellipticity, leading to hyperbolic behavior (see below). This phenomenon is well known in 2D problems for ideally plastic materials (Kachanov, 1974), and has been exploited in a recent study of shape-memory polycrystals (Chenchiah and Bhattacharya, 2005). As was shown in this later reference, the hyperbolicity of the equations can be very helpful in interpreting the solutions obtained, and such connections will be made for the specific cases to be considered below. In any case, it should be kept in mind that the solutions of interest here are limiting cases of standard elasticity problems (with positive-definite elasticity tensors) for which solutions have been obtained by numerical methods (Willot et al., 2007). 2.4. Hyperbolicity conditions and characteristics Making use of the constitutive relations (5), the stress equilibrium conditions may be written in terms of the displacement field as a closed system of two second-order partial differential equations (PDEs): !     l 0 mþk 0 2 0 1 qx u þ q2y u þ ðl þ k  mÞ qx qy u ¼ 0. (10) 0 mþk 0 l 1 0 Following Otto et al. (2003), this system is decoupled as Dui ¼ 0;

Dij ¼ 0;

Dsij ¼ 0;

i ¼ x; y; j ¼ x; y,

(11)

ARTICLE IN PRESS F. Willot et al. / J. Mech. Phys. Solids 56 (2008) 1245–1268

1250

where the PDEs for the strain and stress fields follow immediately, and where D is the fourth-order differential operator: D ¼ q4x þ q4y þ 2rq2x q2y ;

r¼

lðm  kÞ þ 2km . lðm þ kÞ

(12)

This PDE has the symbolic form fðx; yÞ ¼ x4 þ y4 þ 2rx2 y2 , with characteristic lines (e.g., Zachmanoglou and Thoe, 1986) defined by the equation fðdx; dyÞ ¼ 0, such that: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi dx ¼  r  ðr2  1Þ1=2 . (13) dy Taking l, m, kX0, r lies on the segment line r 2 ½1; 1Þ. At any interior point 1oro1, the quantity dx=dy is strictly complex; there are no characteristic lines in the real plane and the problem is elliptic. This corresponds to either (i) 0ol; mo1; (ii) 0omo1, l ¼ 1, ko1; (iii) 0olo1, m ¼ 1, ko1. At the end points r ¼ 1 and r ¼ 1, the form fðdx; dyÞ has two distinct real roots for dx=dy, each of them of multiplicity two. Thus, the problem becomes hyperbolic with straight characteristics x=y ¼ cst: Case 1 : r ¼ 1 3 l ¼ 0 or m ¼ k ¼ 1 3 x; y ¼ cst, Case 2 : r ¼ 1 3 m ¼ 0 or l ¼ k ¼ 1 3 x ¼ y þ cst.

ð14Þ ð15Þ

In this study, only hyperbolic problems are considered, i.e. the medium will be taken incompressible (k ¼ 1) when l or m ¼ 1. All characteristics encountered are straight lines aligned with either one of the Cartesian axes (a ¼ 0) or with one of the two diagonals of the unit cell (a ¼ 1). As will be seen, the solutions found for the stress, strain, or displacement fields actually verify simpler, first-, or second-order PDEs which can be used instead of (11). 2.5. Averages, standard deviations, and field distributions For any stress or strain component aðxÞ produced by a loading a, normalized phase averages and standard deviations in phase b are defined as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi M ðbÞ ðaÞ  haiðbÞ =a; S ðbÞ ðaÞ  ha2 iðbÞ  hai2ðbÞ =a. (16) with the The distribution (i.e. the histogram) of a in the matrix M, Pa ðtÞ, is obtained by counting occurrences R Dirac distribution d. All individual averages in the matrix result from the definition hgðaÞið1Þ ¼ dtgðtÞPa ðtÞ, where Z Pa ðtÞ  ð1  f Þ1 d2 xdðaðxÞ  tÞ. (17) M

3. Material with anisotropy ratio a ¼ 0 3.1. Loading in pure shear In this section, a ¼ 0 with l ¼ 0. A PS stress loading s ¼ hsPS i is applied. 3.1.1. Stress fields Using (5), l ¼ 0 implies sSS ¼ sxy  0. Then, from stress equilibrium, sxx ðx; yÞ ¼ gðyÞ;

syy ðx; yÞ ¼ gðxÞ,

(18)

where g is an unknown 1-periodic function (using symmetry A1). Hence sPS;m ¼ ½gðyÞ  gðxÞ=2,

(19)

where the plus (resp. minus) sign applies to sPS (resp. sm ). This sign convention for PS loading, repeatedly employed hereafter, only holds for Section 3.1. The opposite convention will apply in Sections 3.2 and 3.3.

ARTICLE IN PRESS F. Willot et al. / J. Mech. Phys. Solids 56 (2008) 1245–1268

1251

We note that each stress component sii (no summation) obeys a simple first-order PDE, where x ¼ cst and y ¼ cst are characteristic lines, in agreement with (14). The stress vanishes in the void so that gðxÞ  0 for x 2 ½a; a. The structure of the solution is most easily grasped referring to Fig. 2(a). It is organized in three types of zones: (A), (B), and (D), separated by ‘‘frontiers’’ marked by solid lines. The PS stress reduces to sPS ¼ gðxÞ=2 or gðyÞ=2 in zone (B), to sPS ¼ 0 in zone (D), and to sPS ¼ ½gðxÞ þ gðyÞ=2 in zone (A). It thus vanishes in the square of length 2a consisting in the union of (D) and of the void (V ). This is more easily understood in terms of characteristics. Since the transverse stress component sxy is zero, the void boundary conditions for the stress reduces to sxx ¼ syy ¼ 0 at any point of the void–matrix interface. Hence, each of these components is zero along one of the two families of characteristics (vertical or horizontal) passing through a void. At the intersection of characteristic lines passing through a void (i.e. region D), all stress components must be zero. Thus, as far as stress is concerned, the voids behave as square voids. In particular, the effective modulus (34a) must be a natural function of a rather than of f (i.e., must not contain p when expressed in terms of a). The unknown function g is obtained by minimizing the complementary elastic energy functional (9), which is a strictly convex problem in the subspace of periodic stress fields with vanishing SS component (i.e., sSS ¼ 0):   Z 1 2 1 2 ð1Þ sPS ðx; yÞ þ sm ðx; yÞ dx dy, ð1  f Þhu ðrÞið1Þ ¼ 2k 1=2px;yp1=2 2m #2   Z 1=2  "Z 1=2 1 1 1 1 1 2 þ  ¼ g ðxÞ dx þ gðxÞ dx . 2 m k a m k a Functionally extremizing this integral with respect to g under the constraint Z 1=2 hsPS i ¼ s ¼ 2 gðxÞ dx

(20)

a

provides gðxÞ  s=ð1  2aÞ for x 2 ½12; a [ ½a; 12. Denoting by y½x1 ;x2  ðzÞ the characteristic function of the interval ½x1 ; x2 , we introduce wðzÞ  1  y½a;a ðzÞ, the characteristic function of the domain ½12; a [ ½a; 12. Then, sPS and sm are completely determined by sPS;m ðx; yÞ ¼ s

wðyÞ  wðxÞ . 2ð1  2aÞ

(21)

Referring to Fig. 2(a), sPS meets its highest (constant) value in zones (A), is of half this value in zones (B), and zero elsewhere. The discontinuous stress pattern obtained here as a solution for infinite anisotropy is of the type used by Drucker (1966) in his limit analysis of the ideally plastic porous medium. 3.1.2. Strain field From (6) and (21), the strain components PS and m in the matrix read PS ðx; yÞ ¼ s

wðyÞ þ wðxÞ ; 4mð1  2aÞ

m ðx; yÞ ¼ sm

wðyÞ  wðxÞ . 4mð1  2aÞ

(22)

Fig. 2. Structure of field patterns in the unit cell, in situations of infinite anisotropy: (a) pattern for a ¼ 0 and (b) a ¼ 1 (cf. Section 4).

ARTICLE IN PRESS F. Willot et al. / J. Mech. Phys. Solids 56 (2008) 1245–1268

1252

The relationship between s and  ¼ hPS i is now obtained. To compute , we use u . Expressions (22) equivalently read [cf. Equation (A1e)]: xx ðx; yÞ ¼ yy ðy; xÞ ¼

1 s ½ð1  mÞwðxÞ þ ð1 þ mÞwðyÞ. 4m 1  2a

(23)

In PS loading, the admissibility (i.e., compatibility) conditions in (7) imply that qx u x ¼ xx   and qy u y ¼ yy þ . The displacement component u x (resp. u y ) is odd and 1-periodic wrt. x (resp. y, see Eqs. (A1)). This requires u to be tangent to the boundary of the unit cell: u x ð12; yÞ  u y ðx; 12Þ  0. Hence, the previous equations are integrated as Z x Z y 0 0 ux ðx; yÞ ¼ dx ½xx ðx ; yÞ  ; uy ðx; yÞ ¼ dy0 ½yy ðx; y0 Þ þ . (24) 1=2

Furthermore, the periodicity condition Z 1=2 dx½xx ðx; yÞ   ¼ 0.

1=2

u x ð12; Þ

 0 yields: (25)

1=2

Choosing, e.g., y ¼ 12 (in the matrix) and inserting (23) in (25) provides the desired relation between s and : s¼

2mð1  2aÞ . 1 þ ðm  1Þa

(26)

Since m40 and 0pap12, the denominator of (26) is always 40. The still unknown xy  SS is computed from the admissibility (i.e., compatibility) conditions in (7), and from the expressions of u x and u y in (24). We obtain for ðx; yÞ in the matrix: "Z # Z y x 1 0 q 0 0 q 0 xx ðx ; yÞ þ yy ðx; y Þ . xy ðx; yÞ ¼ dx dy (27) 2 1=2 qy qx 1=2 Inserting solutions (23), the integration is first carried out in the LLQ for integration paths in the matrix. The result is extended to the whole unit cell appealing to identities (A1f), (A1g). For instance, in the URQ:       1þm 1 1 SS ðx; yÞ ¼ x  dðy  aÞ  y  dðx  aÞ , (28) 4 1 þ ðm  1Þa 2 2 where the Dirac distributions stem from the discontinuities in displacements (24) due to the w functions in (23). Hence, in PS loading, SS is localized on the four ‘‘frontier’’ lines in the unit cell. Points ðx; yÞ ¼ ða; 0Þ and ð0; aÞ are singular hot spots. At each tangency point of a ‘‘frontier’’ line with the void boundary, there occurs a jump of SS along the line. E.g., in the vicinity of ða; 0Þ, (28) and the symmetry xy ðx; yÞ ¼ xy ðx; yÞ provide: SS ðx; yÞ ’

 1þm dðx  aÞ sign y. 8 1 þ ðm  1Þa

(29)

3.1.3. Displacement field and ‘‘hot spots’’ From (24) and identity u x ðx; yÞ ¼ u x ðx; yÞ, the displacement u corresponding to the above strains reads in the URQ (y is the Heaviside function):     1 1  x ðm  1Þa þ ðm þ 1Þyða  yÞ ux ðx; yÞ ¼ 1 þ ðm  1Þa 2 2  1  ðm  1Þða  xÞyða  xÞ . ð30Þ 2 Moreover, u y ðx; yÞ ¼ u x ðy; xÞ. Again referring to Fig. 2(a), u is a linear function of x, y in (A), (B), and (D). Along the ‘‘frontiers’’ between (A) and (B) or between (B) and (D), its component normal to the ‘‘frontier’’ is

ARTICLE IN PRESS F. Willot et al. / J. Mech. Phys. Solids 56 (2008) 1245–1268

1253

continuous, whereas its tangential component is discontinuous. For instance, on the ‘‘frontier’’ y ¼ a for aoxo12, u x undergoes a jump ½½u x y  u x ðx; aþ Þ  u x ðx; a Þ given by   1 ð1 þ mÞ  x ; u1  . (31) ½½u x y ðx; aÞ ¼ u1 2 2½1 þ ðm  1Þa As expected, the displacement field can develop discontinuities along characteristic lines. This feature is reminiscent of rigid, ideally plastic bodies in isotropic 2D materials where discontinuities may only occur tangentially to slip lines (Kachanov, 1974, Proposition 39.4). The structure of the full deformed configuration is schematized in Fig. 3: overall deformation occurs by block sliding. Fig. 3 is drawn in the particular case of equal bulk and shear moduli, k ¼ m (m ¼ 1). Then, u vanishes strictly in zone (A)—but is finite there if kam—and is aligned with the axes in (B). The expression of u in (D) clarifies the nature of the hot spots: uðx; yÞ ¼ u1 ðsign x; sign yÞ;

aox; yoa.

(32)

We stress that this incomplete expression, valid for aox; yoa only, does not allow one to recover the singularity (29) of xy by taking a derivative. Since u is constant, save for orientation changes, in the four quadrants of zone (D), the hot spots are either points of extreme matter separation (at ðx; yÞ ¼ ð0; aÞ for 40; white elliptic markings in Fig. 3) or of matter crushing (at ðx; yÞ ¼ ð0; aÞ; dark markings). Fig. 3 also makes conspicuous four voided squares (in black) generated by the block-sliding pattern, at the intersections ðx; yÞ ¼ ða; aÞ of the ‘‘frontier’’ lines (the phenomenon is most remarkable for m ¼ 1). The immediate vicinity of each of these points is divided into four regions where the displacement vector locally takes on distinct values. For instance, around ðx; yÞ ¼ ða; aÞ, we have u ðaþ ; aþ Þ ¼ ðu0 ; u0 Þ, u ðaþ ; a Þ ¼ ðu0 þ Du; u0 Þ, u ða ; aþ Þ ¼ ðu0 ; u0  DuÞ, u ða ; a Þ ¼ ðu0 þ Du; u0  DuÞ, where u0 ¼ 2u1 ½ð1  mÞ=ð1 þ mÞaða  12Þ;

Du ¼ u1 ð12  aÞ.

(33)

The quantity Du represents the size of the voided square, to first order in . Remark that in terms of s, it reads Du ¼ sð1 þ mÞ=ð8mÞ, a void-independent expression (the dimensioning factor is the cell size L ¼ 1). eðf Þ 3.1.4. Distributions, moments and effective shear modulus m In the PS case, the distributions Pz ðtÞ or Psz ðtÞ consist uniquely of a sum of Dirac-type components, plus unusual singular components at z ! 1, due to the Dirac singularities in the fields. The latter cannot be accounted for by probability densities unless inconvenient limiting processes are employed. This may constitute a limitation of the use of distributions in a localizing regime. However, thepfirst ffiffiffiffiffiffiffiffi and second moments are readily computed. The effective shear modulus e m is read from (26). With a  f =po12, we have e 1  2a m ¼ 1  ð1 þ mÞðf =pÞ1=2 þ Oðf Þ, ¼ m 1 þ ðm  1Þa

(34a)

Fig. 3. Structure of the deformed matrix and of u (arrows) for an anisotropy ratio a ¼ 0, in the particular case of equal bulk and shear moduli k ¼ m. Loading is PS. Voids are in black. White (resp. gray) elliptic rings tag zones of extreme matter separation (resp. crushing).

ARTICLE IN PRESS F. Willot et al. / J. Mech. Phys. Solids 56 (2008) 1245–1268

1254

M ð1Þ ðPS Þ ¼

1  2a ; ð1  f Þ½1 þ ðm  1Þa

M ð2Þ ðPS Þ ¼

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð1  2aÞ½ð1 þ f Þa  f  S ðPS Þ ¼ ; ð1  f Þ½1 þ ðm  1Þa ð1Þ

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi að1  2aÞ pffiffiffiffiffiffiffiffiffiffiffi ; S ðm Þ ¼ ½1 þ ðm  1Þa 1  f ð1Þ

Sð1Þ ðSS Þ ¼ 1;

m

ðm þ 1Þa , f ½1 þ ðm  1Þa

(34b)

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð1 þ f Þa  f pffiffiffiffiffiffiffiffiffiffiffiffiffiffi , S ðsPS Þ ¼ ð1  f Þ 1  2a ð1Þ

Sð1Þ ðsm Þ ¼

(34c)

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a , ð1  f Þð1  2aÞ

(34d)

Sð1Þ ðsSS Þ ¼ 0.

(34e)

The mean strain in the pore stems from M ð2Þ ðPS Þ ¼ ½1  ð1  f ÞM ð1Þ ðPS Þ=f . Moreover, Sð1Þ ðsPS Þ ¼ ðm=e mÞSð1Þ ðPS Þ. The Dirac singularities are responsible for the transverse variance S ð1Þ ðSS Þ being infinite. Also, remark that M ð2Þ ðPS Þ blows up as f 1=2 when f ! 0. The curve e mðf Þ=m is displayed in Fig. 4 in the incompressible case m ¼ 0. The power singularity at f ¼ 0 with infinite negative slope, due to the f 1=2 term in the dilute expansion of (34a), indicates that an infinitesimal void dramatically weakens the medium. This is reminiscent of the situation encountered in plasticity (see Introduction). In the linear material considered here, it is a direct consequence of hyperbolicity. Each of the two components sxx and syy have constant values along one of the two families of characteristics. Due to the boundary conditions at the void–matrix interface, the parallel component sPS is smaller by half in region (B) than that in (A). Hence, each void lowers the stress field over large regions in the material, which are projections of the voids along the characteristic directions and involve infinitely long distances. The normalized standard deviations Sð1Þ ðPS Þ, Sð1Þ ðm Þ, S ð1Þ ðsPS Þ, and S ð1Þ ðsm Þ all behave as ðf =pÞ1=4 as f ! 0, which indicates that they grow more rapidly with f than in a linear isotropic medium. The effective modulus e mðf Þ vanishes linearly with ðf c  f Þ for f near f c ¼ p=4, the void close-packing fraction (see the comment in Section 3.2.2), along with the moments of the strain components in the matrix in the loading direction (M ð1Þ ðPS Þ ¼ S ð1Þ ðPS Þ ¼ 0). 3.2. Loading in SS (incompressible case only) According to Section 2.3, SS loading goes along with a ¼ 0 and m ¼ 1. Matrix incompressibility renders the problem hyperbolic and is assumed for simplicity (k ¼ 1, m ¼ 0). 1 ∼ λ/λ

α=0

0.8 ∼ μ/μ

0.6

0.4

0.2

0 0

0.2

0.4 f

0.6

0.8

e=m and e Fig. 4. Normalized shear effective moduli m l=l vs void concentration f, for an incompressible matrix with anisotropy ratio a ¼ l=m ¼ 0. Both moduli vanish at the close-packing threshold f ¼ f c ¼ p=4 ’ 0:78.

ARTICLE IN PRESS F. Willot et al. / J. Mech. Phys. Solids 56 (2008) 1245–1268

1255

3.2.1. Displacement, strain, and stress fields Taking m ¼ 1 under finite stress implies xx  yy ¼ 0. Incompressibility then requires xx ¼ yy ¼ 0. The only non-zero component is xy , and a variational calculation is simpler to carry out wrt. the strain rather than to the stress. Condition xx ¼ yy ¼ 0 implies u x ðx; yÞ ¼ u x ðyÞ and u y ðx; yÞ ¼ u y ðxÞ. By symmetry, see (A2c), u x ðzÞ ¼ u y ðzÞ  GðzÞ  z, where G is unknown. Introducing the derivative g  G 0 , one obtains xy ðx; yÞ ¼ ½gðxÞ þ gðyÞ=2.

(35)

The energy of the unit cell (7) is expressed as an integral over the matrix phase M (represented by the unit square minus the void): Z ð1  f Þhwð1Þ ið1Þ ¼ 2l d2 x2xy (ZM ) Z l ¼ dx dy  dx dy ½gðxÞ þ gðyÞ2 . ð36Þ 2 ½1=2;1=2½1=2;1=2 V After using identity (A2f), we expand it into separate integrals over the intervals ½0; a and ½a; 1=2, and we split g into independent functions gA and gB supported by these intervals. Subscripts A, B refer to the zones in Fig. 2(a): gðzÞ ¼ gB ðzÞy½0;a ðzÞ þ gA ðzÞy½a;1=2 ðzÞ.

(37)

The strain energy (7) is functionally minimized wrt. gA , gB under the constraint: "Z # Z a 1=2 hxy i   ¼ 2 gA ðzÞ dz þ gB ðzÞ dz . a

(38)

0

The system obtained determines gA as a constant. Introducing rðzÞ  equation for gB ðzÞ: gA ðzÞ  gB ðaÞ; z 2 ½a; 12, Z a 2 gB ðyÞ dy ¼  þ 2agB ðaÞ þ ½2rðzÞ  1gB ðzÞ;

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a2  z2 , it provides an integral ð39aÞ

z 2 ½0; a.

ð39bÞ

rðzÞ

In particular, for z ¼ a, gB ðaÞ is expressed in terms of gB ðzÞ, z 2 ½0; a½: Z a dzgB ðzÞ ¼  þ ð2a  1ÞgB ðaÞ. 2

(40)

0

Taking z ¼ 0 in (39b) yields a relationship between gB ð0Þ and gB ðaÞ: gB ð0Þ ¼ ½ þ 2agB ðaÞ=ð1  2aÞ.

(41)

Eq. (39b) can be transformed into a second-order differential equation with no simple analytical solution.1 Instead we express the strain, the displacement and the stress in terms of gB , and we directly deduce from (39b) some asymptotic results for the fields and for their distributions. We also compute numerically the whole solutions, the distributions, and the effective moduli by discretizing (39b) into a linear system. The solution for gB ðzÞ is represented vs. 0oz=ao1 for a 2 ½0; 12 in Fig. 5. The inset displays gB ð0Þ as a function of a. At fixed a, the highest values lie at z ¼ 0, the strain being higher on the cartesian axes of the unit cell. Two regimes are encountered as a increases: the strain first develops in zone (B) and concentrates around the axes; then, for a\0:4 (i.e. f \0:5), gB ð0Þ blows up as ðac  aÞ1=2 near ac ¼ 12 (see Section 3.2.3), while the strain localizes on the cartesian axes. In the limit, gB ðzÞ / dðzÞ. The field enhancement between voids near close-packing is boosted by anisotropy. The close-packing void fraction f c ¼ p=4 is, for periodic media, tantamount to the Rx With QðxÞ  0 gB ðyÞ dy, (39b) relates linearly Qða  rðzÞÞ to Q0 ðzÞ. Differentiating it wrt. z yields another equation connecting 0 0 Q ða  rðzÞÞ, Q ðzÞ and Q0 ðzÞ. The substitution y ¼ rðzÞ admits z ¼ rðyÞ as an inverse, and turns the first equation into one linking QðzÞ to 0 Q ða  rðzÞÞ and the second one into a differential equation for QðzÞ. 1

0

ARTICLE IN PRESS F. Willot et al. / J. Mech. Phys. Solids 56 (2008) 1245–1268

1256

6 gB (0) / ε

6

4

4

gB (z) /ε

2 1 0

0.2

2

0.4 a

1 0 -1 0

0.2

0.4

0.6

0.8

1

z/a Fig. 5. Plots of the solution gB ðzÞ of (39b) as a function of z=a, 0oz=ao1, for values of the void radius (from bottom to top at z ¼ 0) a ¼ 0:025 (black), 0:124 (red), 0:249 (green), 0:275 (blue), 0:373 (gray), 0:435 (violet), 0:485 (cyan), 0:498 (rose), 0:4995 (orange), 0:499985 (dark green). Color online. Line and symbol thicknesses increase with a.

void mechanical percolation threshold in random porous media (e.g., Torquato, 2002). When a ! ac , gB ðzÞ becomes discontinuous, with gB ðaÞ ¼ 1 and gB ða Þ ¼ 0, see Fig. 5. The solution in terms of gB is as follows. With (39a), gðzÞ in (37) becomes completely determined in terms of gB ðzÞ, z 2 ½0; a. In the URQ: xy ðx; yÞ ¼

gB ðaÞ 1 ½y½a;1=2 ðxÞ þ y½a;1=2 ðyÞ þ ½gB ðxÞy½0;a ðxÞ þ gB ðyÞy½0;a ðyÞ. 2 2

(42)

The integral equation (39b) admits a continuous solution. Then, gB and in turn xy , u and u are continuous, too. The continuity and oddness wrt. x (resp. y) of u x (resp. u y ) implies u x ð0; yÞ  u y ðx; 0Þ  0. Hence, in the URQ: Z y ux ðx; yÞ ¼ dzgðzÞ  y; u y ðx; yÞ ¼ u x ðy; xÞ ð0ox; yo12Þ. (43) 0

These expressions are extended to the unit cell using (A2a)–(A2c). From (A2d) the transverse components sxx and syy are odd in x and y, so that periodicity requires sxx ð12; yÞ  syy ðx; 12Þ  0. With sxy ðx; yÞ ¼ 2lxy ðx; yÞ given by (42), we deduce from stress equilibrium in the URQ of the matrix: Z sxx ðx; yÞ ¼ syy ðy; xÞ ¼

1=2

dx0 qy sxy ðx0 ; yÞ ¼ ðl=2Þð1  2xÞg0 ðyÞ.

(44)

x

Then in this region: sm;PS ðx; yÞ ¼ ðl=4Þ½ð1  2xÞg0 ðyÞ  ð1  2yÞg0 ðxÞ, where the plus (resp. minus) sign applies to sm (resp. sPS ). Hence, because gB ðzÞ ¼ const. for z 2 ½a; 12, the stresses sm and sPS vanish in zones (A) of Fig. 2(a). In zones ðBÞ, sm and sPS are closely related, with ( aoxo12 0 sm ðx; yÞ ¼ sPS ðx; yÞ ¼ ðl=4Þð1  2xÞgB ðyÞ for , ð45aÞ 0oyoa ( 0oxoa 0 . ð45bÞ sm ðx; yÞ ¼ sPS ðx; yÞ ¼ ðl=4Þð1  2yÞgB ðxÞ for aoyo12

ARTICLE IN PRESS F. Willot et al. / J. Mech. Phys. Solids 56 (2008) 1245–1268

1257

Taking the derivative of (39b) with respect to z, we deduce g0B ðzÞ ¼ 2z

gB ðzÞ þ gB ðrðzÞÞ , rðzÞ½1  2rðzÞ

(46)

so that g0B ð0Þ ¼ 0. Hence besides zones (A), regions of weak sm and sPS lie along the axes x ¼ 0 or y ¼ 0, and along the axes x ¼ 12, y ¼ 12. Near the ‘‘frontiers’’ in Fig. 2(a), sm and sPS blow up. Indeed, let us focus on the ‘‘frontier’’ y ¼ a: From (46), we deduce for yta pffiffiffiffiffi g0B ðyÞ ’  2a½gB ðaÞ þ gB ð0Þða  yÞ1=2 , (47) and the latter stresses behave as d 1=2 , where d is the distance to the ‘‘frontier’’ line. The variances involving an integral with integrand / s2 , this implies: Sð1Þ ðsPS Þ ¼ S ð1Þ ðsm Þ ¼ 1.

(48)

These divergences are the only ones encountered in the means and variances. The mean stress readily derives from Eqs. (39)–(40) and (42): "Z Z 1=2 Z a Z rðxÞ # 1=2 s ¼ hsxy i ¼ 8l dx dy  dx dy xy ¼ l½1 þ gB ðaÞ=, (49) 0

0

0

0

The effective shear modulus is therefore e l ¼ ðl=2Þ½1 þ gB ðaÞ=. 3.2.2. Moments and effective shear modulus e lðf Þ in the dilute limit Series expansions are now carried out for a51. Using the Taylor expansion: gB ðazÞ ¼ 

1 X

qn ðzÞan

(50)

(51)

ð0ozo1Þ,

n¼0

we solve Eq. (39b). Note that a being a parameter, we should have denoted gB ðzÞ by gB ðz; aÞ, so that in the above equation, gB ðazÞ would stand for gB ðaz; aÞ, expanded in powers of a. The following recursion is obtained: Z 1 pffiffiffiffiffiffiffiffiffiffiffiffiffi qnþ1 ðzÞ ¼ qn ð1Þ þ 1  z2 qn ðzÞ  pffiffiffiffiffiffiffiffi dyqn ðyÞ; q0 ðzÞ  1. (52) 2 1z2 Numerically, we find that series (51) converges for aoRa ’ 0:4. In the convergence region, this solution matches the numerical one obtained in the previous section. We carry out analytically the recursion up to fourth order to estimate gB ðazÞ=. The integrations are readily performed with a symbolic calculator. Using the result with z ¼ 1 provides: 3 2 4 5 gB ðaÞ= ¼ 1  2pa2  64 3 a þ 2ðp  6p  8Þa þ Oða Þ.

With expansion (53), we find in the dilute limit:   e l 32f 3=2 6 8 2 ¼1f   þ 1  f þ Oðf 5=2 Þ, l p p2 3p3=2 M ð1Þ ðSS Þ ¼ 1 

32f 3=2 þ Oðf 2 Þ; 3p3=2

M ð2Þ ðSS Þ ¼ 1 þ

(53)

(54a) 32f 1=2 þ Oðf Þ, 3p3=2

(54b)

pffiffiffi 4 2f 3=4 S ðSS Þ ¼ pffiffiffi þ Oðf 5=4 Þ; 3p3=4

pffiffiffi 4 2f 3=4 S ðsSS Þ ¼ pffiffiffi þ Oðf 5=4 Þ, 3p3=4

(54c)

Sð1Þ ðPS Þ ¼ 0;

Sð1Þ ðm Þ ¼ 0;

(54d)

ð1Þ

Sð1Þ ðsPS Þ ¼ 1;

ð1Þ

S ð1Þ ðsm Þ ¼ 1.

ARTICLE IN PRESS F. Willot et al. / J. Mech. Phys. Solids 56 (2008) 1245–1268

1258

A plot of e lðf Þ=l obtained from the full numerical solution of gB ðzÞ is displayed in Fig. 4. Compared to e mðf Þ=m, the less singular character of the solution goes along with a harder material at small porosities. But e l falls down to zero much faster than e m for higher porosities, except near the close-packing threshold. 3.2.3. Scaling in the close-packing limit Day et al. (1992) examined the linear elastic behavior of an isotropic material containing a honeycomb lattice of voids: the effective compressibility and shear moduli (the overall medium is isotropic) vanish as ðf c  f Þ1=2 near the void-packing threshold. A similar behavior is found here in shear for the anisotropic material. With the numerical solution of (39b) near f c , we observe that the function gB ðzÞ obeys a scaling of the form gB ðazÞ ðf c  f Þ1=2 geðzðf c  f Þ1=2 Þ;

f ! f c,

(55) t

where the master curve geðzÞ const. for z51 and geðzÞ z for zb1, where t is an exponent close to 2. Fig. 6 provides an illustration for three different values of a near ac ¼ 12. Two consequences are drawn from (55): first, the scaling shows that the strain is concentrated in a band of width x aðf c  f Þ1=2 , vanishing as the close-packing threshold is approached. Next, relation (55) used at z ¼ 0 together with (41) and (50) provides near f c : e l=l ðf c  f Þ1=2 ,

(56a)

as for an isotropic material. The means are R readily computed. Moreover, the variances are obtained with the help of the surface integral in the matrix M dx dy2SS ðx; yÞ ¼ 2 ½1 þ gB ðaÞ==2. The following behaviors near the close-packing threshold ensue: M ð1Þ ðSS Þ ðf c  f Þ1=2 ; Sð1Þ ðPS Þ ¼ S ð1Þ ðm Þ  0;

S ð1Þ ðSS Þ ðf c  f Þ1=4 ;

Sð1Þ ðsSS Þ ðf c  f Þ1=4 ,

S ð1Þ ðsPS Þ ¼ S ð1Þ ðsm Þ  1.

(56b)

3.2.4. Distributions and Van Hove singularities The distributions of the stress fields in the matrix for SS loading, at void concentration f ¼ 0:1, are displayed in Figs. 7 and 8 as thick solid lines. In these figures, straight vertical lines represent Dirac components. The distribution PsSS ðtÞ of the stress ‘‘parallel’’ to the applied loading (Fig. 7) comprises one Dirac component, proportional to ½ð1  2aÞ2 =ð1  f Þdðt  gB ðaÞÞ, which represents the contribution of zones (A) of constant sSS , and a divergence followed by a discontinuity of infinite amplitude. The latter comes from the (B) zones. The right ‘‘foot’’ (magnified in the inset) is produced by zones (D) of maximal a=1/2-2.5 10-4 a=1/2-7.6 10-6 a=1/2-10-6 (1/2-a)1/2gB (az)/ ε

(1/2-a)1/2 gB (az) / ε

0.3

0.2

0.1

10-1 10-2 10-3 10-4 10-5 10-6 10-7 10-2 10-1 100 101 102 103 z / (1/2-a)1/2

0.0 0

10

20

30

z / (1/2-a)1/2

Fig. 6. Data collapse illustrating the scaling properties of the function gB ðzÞ, for the three different values of the radius a indicated in the legend, near the close-packing value ac ¼ 12. Inset: same plot in log–log scale; the curves break down when gB ðazÞ become negative.

ARTICLE IN PRESS F. Willot et al. / J. Mech. Phys. Solids 56 (2008) 1245–1268

1259

12.0 f = 0.1

10.0 Pσ (t)

P (t)

Pσ (t)

PS

8.0

SS

PS load.

SS load.

α=∞

α=0

6.0

0.10

4.0

Pσ

PS

(t)

0.05

2.0

0.00 1.6 1.7 1.8

0.0 0.8

1.0

1.2

1.4 1.6 t = σPS, SS, m

1.8

2.0

Fig. 7. Distributions of the stress components parallel to the applied loading, at void concentration f ¼ 0:1, in simple shear loading (hsSS i ¼ 1) with a ¼ 0 (PsSS , solid), and in pure shear loading (hsPS i ¼ 1) with a ¼ 1 (PsPS , dashed). Straight vertical lines near t ¼ 0:8 indicate Dirac components. Inset: magnification of the ‘‘foot’’ in PsSS ðtÞ.

0.20 f = 0.1 Pσ (t)

0.15

Pσ (t)

SS

P (t)

PS

PS load. α=∞

0.10

SS load. α=0

Pσ (t) m

0.05

Pσ (t) m

PS load. α=∞

SS load. α=0

0.00 -3.0

-2.0

-1.0

0.0 1.0 t = σPS, SS, m

2.0

3.0

Fig. 8. Distributions of the stress components transverse to the applied loading, and of sm , at void concentration f ¼ 0:1, in simple shear loading (hsSS i ¼ 1) with a ¼ 0 (PsPS ðtÞ ’ Psm ðtÞ for jtjt1:5), and in pure shear loading (hsPS i ¼ 1) with a ¼ 1 (PsSS ðtÞ ’ Psm ðtÞ). Dirac components at t ¼ 0 should also be present in each of the four sets, due to the void and to zone A, but are omitted for clarity.

stress. The distribution has finite support since SS is bounded. The distributions of the transverse stresses sPS and of sm are displayed in Fig. 8. They are even (these fields average to 0). A Dirac contribution at s ¼ 0 (that of zone (A)) should be present in all plots but is omitted for clarity. The main components blow up at s ¼ 0 and slowly decrease at s ! 1 (the contribution of zones (B), essentially). The distribution Psm differs only slightly from PsPS by the contribution of zones (D), the difference being small at f ¼ 0:1, for the smallest values of the stress. The above singular features are understood as follows. Tails of the distributions at high values of the fields are determined by the square-root singularities responsible for the infinite variances (48). Definition (17) of the distribution generates Van Hove singularities (VHS), in the vicinity of values t ¼ t0 for which there exists extremal or saddle integration points xi such that aðxi Þ ¼ t0 and qx aðxi Þ ¼ 0. In 2D, depending on the eigenvalues of the matrix q2xx aðxi Þ being of same or of opposite signs, the singularity in P is either a discontinuity or a logarithmic divergence (Van Hove, 1953). ‘‘Extended’’ Van Hove singularities (EVHS) (Abrikosov et al., 1993) arise when in addition one or more eigenvalues of q2xx aðxi Þ vanish. Consider first the strain distribution. Using (39b), we obtain gB ðzÞ ¼ gB ð0Þ 

 þ gB ðaÞ 2 z þ Oðz4 Þ, að1  2aÞ2

(57)

ARTICLE IN PRESS F. Willot et al. / J. Mech. Phys. Solids 56 (2008) 1245–1268

1260

so that in zone (B), for aoxo12, 0oyoa, xy ðx; yÞ ’

 þ gB ðaÞ  þ gB ðaÞ 2  y þ , 2ð1  2aÞ 2að1  2aÞ2

(58)

of the generic type hðx; yÞ ¼ h0  h1 y2 , h1 40. This generates an EVHS near t ¼ h0 ¼ ð þ gB ðaÞÞ=½2ð1  2aÞ, of the form Z dx dy Lx pffiffiffiffiffi ðh0  tÞ1=2 yðh0  tÞ, dðh0  h1 y2  tÞ ¼ (59) Ph ðtÞ ’ V 2V h1 where Lx is the size of the integration domain on x values. Only values near y ¼ 0 contribute and the range of integration over y need not be completely specified. The four quadrants account for Lx ¼ 2ð1  2aÞ, so that pffiffiffiffiffi 2að1  2aÞ2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi jt  h0 j1=2 yðh0  tÞ. Pxy ðtÞ ’ (60) ð1  f Þ  þ gB ðaÞ Hence the divergence of PsSS , of square-root nature, and the drop of infinite amplitude in Fig. 7 are captured by this approximation, since PsSS ðtÞ ¼ Pxy ðt=ð2lÞÞ=ð2lÞ.

(61)

Consider next the field distributions PsPS ðtÞ and Psm ðtÞ of Fig. 8, and examine the contributions near t ¼ 0. According to (45) and to the remark following, weak field values lie on the axes x ¼ 0, y ¼ 0, x ¼ 12, y ¼ 12, where the four points ðx; yÞ ¼ ð12; 0Þ and ðx; yÞ ¼ ð0; 12Þ contribute most. Focus, for instance, on ðx; yÞ ¼ ð12; 0Þ. From (45), (46), and (49), and for y\0: sm ðx; yÞ ¼ sPS ðx; yÞ ’ s

ð1=2  xÞy . að1  2aÞ2

(62)

pffiffiffiffiffiffiffiffiffi This is a VHS of type hðx; yÞ ¼ h1 x2  h2 y2 (rotated 45 ) with h1 , h2 40 so that Ph ðtÞ  log jtj=ð2V h1 h2 Þ for t 0, where V is the integration surface. Using (62), the contributions of the above four points gather into:   4að1  2aÞ2 2 jtj ðBÞ ðBÞ PsPS ðtÞ Psm ðtÞ  log 2ð1  2aÞ ; t ! 0. (63) sð1  f Þ s The slowly decaying tails at high t, of PsPS ðtÞ and Psm ðtÞ in Fig. 8 are due to singularity (47). Again focusing on the domain aoxo12, 0oyoa in zone (B), and expanding near y ¼ a, we obtain from (45), (46), and (49): pffiffiffiffiffi l s 2a 1  2x 0 sPS ðx; yÞ ¼ ð1  2xÞgB ðyÞ  pffiffiffiffiffiffiffiffiffiffiffi ; yta. 4 4ð1  2aÞ a  y pffiffiffiffiffiffiffiffiffiffiffi Of type hðx; yÞ ¼ h0 ð1  2xÞ= a  y, h0 40, this expression contributes for   Z Z a 1 1=2 h20 ð1  2aÞ3 ð1  2aÞh0 pffiffiffi dx dydðhðx; yÞ  tÞ ¼ yðtÞy jtj  , (64) V a 3V a jtj3 0 where the lower integration bound on y provides the restriction on jtj which delimits the domain of validity of the approximation. Adding contributions from negative values, and multiplying by 4 due to square symmetry leads to: 2 PðBÞ sPS ðtÞ ’ s

að1  2aÞ 3 jtj ; 6ð1  f Þ

s jtjb pffiffiffi . 2 2

(65)

This is only the contribution of ðBÞ. The tail should also include a contribution from ðDÞ at high stresses, also decaying as jtj3 , but harder to handle due to sign changes at the points ða; pffiffiffi aÞ. Other singularities exist which have not been fully investigated: values of sm around x; y ¼ a= 2 produce a singularity in its histogram (extra wings starting at sm ’ 1:7 with vertical slope, see Fig. 8); moreover, for porosities f higher than 0.6, two inflection points show up in gB ðzÞ. They give rise to two local extrema of sm and sPS along the cartesian axes zone (B), also producing extra wings of the above type (not shown).

ARTICLE IN PRESS F. Willot et al. / J. Mech. Phys. Solids 56 (2008) 1245–1268

1261

3.3. Equibiaxial loading and limit l ! 0 We close the study of the case a ¼ 0 with equibiaxial loading  ¼ hm ia0, hPS i ¼ hSS i ¼ 0. The anisotropic limit a ¼ 0 follows either from letting l ! 0, or m ! 1, so that two different solutions may a priori arise. Only the case l ! 0, with finite k, is detailed in this section. Case m ! 1 is briefly addressed in Section 3.4. A method similar to that of Section 3.1 for PS loading is used. The differences relative to the PS case can be related to the symmetries at play. 3.3.1. Stress and strain fields The limit l ¼ 0 under finite macroscopic mean stress requires that sxy  0. Equilibrium and the equivalence of the Ox and Oy axes imply that sxx ðx; yÞ ¼ gðyÞ;

syy ðx; yÞ ¼ gðxÞ,

(66)

where gðzÞ is defined on the interval ½12; 12. Similar steps apply as in Section 3.1. Energy minimization under the same macroscopic stress constraint (20)—with s ¼ hsm i replacing hsPS i, provides sSS ðx; yÞ ¼ 0, and sm;PS ðx; yÞ ¼ s

wðyÞ  wðxÞ , 2ð1  2aÞ

(67)

where the plus (resp. minus) sign applies to sm (resp. sPS ). Associated strains follow from (5). Periodicity conditions on u entail: s¼

2mð1  2aÞ . m þ ð1  mÞa

(68)

Eq. (67) reproduces the PS result (21), with the indices PS and m, and with x and y, exchanged. Furthermore, the PS and equibiaxial loading modes exchange the parts played by k and m, which corresponds to interchanging m21=m. This circumstance almost allows for a one-to-one mapping of the results of this section onto those of Section 3.1. For instance, in the URQ:      1 1  1þm . (69) SS ¼ u2 x  dðy  aÞ þ y  dðx  aÞ ; u2 ¼ 2 2 4 m þ ð1  mÞa Eq. (69) differs from (28) by a substitution m21=m, and by a minus sign due to the difference between (A1g) and (A3g), or between (A1c) and (A3c). 3.3.2. Displacement field and singularities The displacement u% in the URQ is such that ux ðx; yÞ is there of the form (30), with m replaced by 1=m, but with now u y ðx; yÞ ¼ u x ðy; xÞ. Hence, u% is again piecewise linear, and tangentially discontinuous at the ‘‘frontiers’’ between ðAÞ, ðBÞ, and ðDÞ. In particular, the jump ½½u x y ðx; aÞ along the y ¼ a ‘‘frontier’’ in the URQ is given by (31) with m replaced by 1=m. With m ¼ 1, the overall pattern, however, differs from that of Fig. 3, by the fact that for o0 (resp. 40), the displacements u% in all parts of ðBÞ are directed towards the void (resp. outwards), the medium being subjected to compression (resp. extension). Moreover, owing to symmetry u y ðx; yÞ ¼ u x ðy; xÞ, the PS mechanism responsible for the formation of the square voids at ðx; yÞ ¼ ða; aÞ changes into one which induces a compressive singularity near these locations, independently of whether the equibiaxial loading mode is compressive or extensive, over a region of size Du ¼ ðm þ 1Þða=2  14Þ=½m þ ð1  mÞa. This is readily seen from the following expression valid in the vicinity of ðx; yÞ ¼ ða; aÞ: uðx; yÞ ¼ u2 ½ða þ 12Þð1; 1Þ þ ða  12Þðsignðy  aÞ; signðx  aÞÞ.

(70)

The displacement in ðDÞ now reads uðx; yÞ ¼ u2 ðsign x; sign yÞ, so that the four ‘‘hot spots’’ ðx; yÞ ¼ ð0; aÞ, ða; 0Þ now all undergo a similar local singularity, either of a compressive, or of an extensive nature.

ARTICLE IN PRESS F. Willot et al. / J. Mech. Phys. Solids 56 (2008) 1245–1268

1262

e 3.3.3. Moments and effective compressibility modulus k The effective compressibility modulus and the moments are provided by expressions (34a)–(34e), with m and e m replaced by k and e k, with m replaced by 1=m, and with the indices PS and m interchanged. For definiteness:   e ð1  2aÞ 1 1þm k ¼ 1 ðf =pÞ1=2 þ Oðf Þ . ¼ (71) m m m þ ð1  mÞa m Thus, e k possesses a dilute-limit correction f 1=2 at finite k, but blows up as f 1=2 in the incompressible limit m ¼ 0. It decays as ðf c  f Þ near f c . 3.4. Equibiaxial loading and limit m ! 1 Consider now the limit m ! 1, with finite l. A finite stress loading is consistent with this limit only if PS  0. This situation is that of a ¼ 0 under SS loading, the solution being of like complexity. In particular, for finite compressibility, it is governed by elliptic equations (see Section 2.4) and is much more complicated than the one of Section 3.3. Accordingly, we restrict ourselves to the incompressible case k ¼ 1. The solution then simplifies and coincides with that of Section 3.3 with m ¼ 0, see Eq. (6). Note that e k is infinite, due to m ¼ 1. 4. Material with anisotropy ratio a ¼ 1 This section is devoted to the case a ¼ 1, where the material is soft along the diagonals. As Fig. 2(b) illustrates, the field patterns are now rotated by 45 . Accordingly, the solutions obtained below could alternatively be derived in a frame where the void lattice is rotated. The corresponding rotation, R45o , would exchange, in the matrix, ESS and EPS , and l and m. In the rotated frame, the following correspondence should then be used: R45 ðvoid latticeÞ 3 fa21=a; and PS loading2SS loadingg.

(72)

Characteristics for the strain, stress or displacement fields are aligned with the two diagonals of the unit cell (see Section 2.4). Another difference with the case a ¼ 0 resides in the appearance of new regions in the matrix, denoted by ðCÞ in the figure, where the ðBÞ bands cross, with no holes in it. 4.1. Loading in SS We impose m ¼ 0, a ¼ 1 and consider SS loading conditions. Due to Eq. (5), m ¼ 0 implies sPS ¼ ðsxx  syy Þ=2  0, so that sxx ðx; yÞ ¼ syy ðx; yÞ  sðx; yÞ, an unknown function. Stress equilibrium implies q2 s =qx2  q2 s =qy2 ¼ 0. Using the identities in Appendix A, we deduce that sðx; yÞ ¼ ½gðx  yÞ gðx þ yÞ=2, where g is an unknown, even, 1-periodic function. The stress sxy ¼ sSS stems from integrating the stress equilibrium equations. The integration constant is zero, since s ¼ 0 in the void. Eventually, with gð0Þ ¼ 0, sSS ðx; yÞ ¼ ½gðx  yÞ þ gðx þ yÞ=2;

sm ðx; yÞ ¼ ½gðx  yÞ  gðx þ yÞ=2.

(73) pffiffiffi Due pffiffito ffi the pffiffiffi void one has, along the main diagonal, sm ðx; xÞ ¼ gð2xÞ=2  0 for xoa= 2. Hence, gðzÞ  0 on ½ 2a; 2a. In turn, stresses (73) vanish in the whole square made of zones ðV Þ and ðDÞ in Fig. 2(b). Due to 1-periodicity and to g being even, the stress pattern possesses a mirror symmetry with respect to the dashed lines of Fig. 2(b). Its fundamental unit cell is therefore the gray square delimited by them. Symmetry (72) makes this tantamount to the PS case of Section 3.1 with a ¼ 0 and l ¼ 0, with a cell size diminished by a pffifficase ffi factor 1= 2 at constant pore radius a and, consequently, with a porosity twice bigger. Effective moduli and e moments are then obtained from expressions (34) of Section 3.1, with m and e m replaced pffiffiffi by l and l, with m replaced by ‘, with the subscripts SS and PS exchanged, and with a replaced by 2a. Deformation in the non-voided crossing zones ðCÞ of zero stress is understood as follows. The obtained solutions for the displacement and strain fields in the ðDÞ zones are trivially continued to the voided zone ðV Þ. These admissible continuations are non-physical except on the void boundary. However, the above symmetry considerations

ARTICLE IN PRESS F. Willot et al. / J. Mech. Phys. Solids 56 (2008) 1245–1268

1263

make clear that, modulo an appropriate translation, these continuations also provide the physical deformation fields in ðCÞ. Mutatis mutandis, the obtained expressions for the fields and the moduli are valid only up to f ¼ f ð1Þ c ¼ p=8, due to the rescaling of a. This concentration acts as a ‘‘mechanically driven percolation’’ threshold, at half the close-packing value. It corresponds to the configuration where the summits ðDÞ in Fig. 2b come into contact with the cell boundaries. Then, ðAÞ and ðBÞ vanish, leaving only ðCÞ and ðDÞ for f up to f c ¼ p=4. Thus, for f 4f ð1Þ c , sSS vanishes in the matrix. The stress SS also vanishes, so that strain only takes place along the void boundaries. Consequently, e l vanishes for f Xf ð1Þ and the fields are such that: M ð1Þ ðSS Þ ¼ 0, c ð1Þ ð1Þ e M ðPS Þ ¼ 1=f , S ðSS Þ ¼ 0. Fig. 9 displays lðf Þ=l vs. f. Except for the threshold, the curve is alike that of e mðf Þ=m in Fig. 4. In words, we just showed that special field configurations in linear anisotropic periodic media can lower the geometric ‘‘percolation’’ threshold of the porous material (i.e. the close-packing threshold, for a periodic pore lattice) into one determined by effective ‘‘porous’’ zones (as far as they undergo vanishing stresses) at intersections of void-generated stress bands. This is unusual: for instance, in homogenization methods, percolation-like thresholds are usually thought (or found) to be independent of the constitutive law. 4.2. Loading in pure shear (incompressible case only) In the case a ¼ 1 with l ¼ 1 and m finite, under PS loading, the solution has the pattern of Fig. 2(b). It is ð1Þ of the type studied in Section 3.2, at least for f up to f ð1Þ c ¼ p=8. For f 4f c the problem becomes more involved due to a complex geometry, and has not been investigated in full. Only the range f of ð1Þ c is discussed hereafter. 4.2.1. Displacements, strain, and stress fields The solution obeys SS ¼ xy  0, and the incompressibility constraint xx þ yy  0. Expressed in terms of u x , u y , such that u x ð0; yÞ ¼ u y ðx; 0Þ  0 as implied by (A1a) and (A1b), these equations entail with the help of (A1c): u x ðx; yÞ ¼ ½Gðx þ yÞ þ Gðx  yÞ=2;

u y ðx; yÞ ¼ ½Gðx þ yÞ þ Gðx  yÞ=2,

(74)

and eventually with gðzÞ ¼ G0 ðzÞ an even and 1-periodic function: PS ðx; yÞ ¼ ½gðy  xÞ þ gðx þ yÞ=2.

(75)

1 ∼ μ/μ

α=∞

0.8

0.6 ∼ λ/λ

0.4

0.2

0 0

0.1

0.2

0.3

0.4

f Fig. 9. Normalized shear effective moduli e m=m and e l=l vs. void concentration f, for an incompressible matrix with anisotropy ratio a ¼ 1. The curve for e m=m is computed up to f ¼ p=8 only. Modulus e l vanishes at f cð1Þ ¼ p=8 ’ 0:39, half the close-packing threshold value.

ARTICLE IN PRESS F. Willot et al. / J. Mech. Phys. Solids 56 (2008) 1245–1268

1264

pffiffiffi pffiffiffi Introducing the 45 counterclockwise-rotated coordinates x0 ¼ ðy þ xÞ= 2 and y0 ¼ ðy  xÞ= 2, such that pffiffiffi p ffiffi ffi 1=ð2 2Þpx0 ; y0 p1=ð2 2Þ, functions gB and gA are introduced, such that pffiffiffi gð 2x0 Þ ¼ gB ðx0 Þy½0;a ðx0 Þ þ gA ðx0 Þy½a;1=ð2pffiffi2Þ ðx0 Þ. (76) The total elastic energy in the unit cell (7) is provided by the integral: Z d2 x2PS ð1  f Þhwð1Þ ið1Þ ¼ 2m (M Z ) Z 0 0 0 0 ¼ 2m 2 dx dy 2PS ðx0 ; y0 Þ. pffiffi pffiffi dx dy  ½1=ð2 2Þ;1=ð2 2Þ2

ð77Þ

V

One extra contribution of the void V (of no elastic energy) has been subtracted from the contribution of the gray square in Fig. 2(b), counted twice. Anticipating on the fact that on its interval of definition gA ðzÞ  gB ðaÞ, as in Section 3.2.1, energy (7) is functionally minimized under the constraint: pffiffiffi pffiffiffi Z a  ¼ hPS i ¼ 2 2 gB ðyÞ dy þ ð1  2 2aÞgB ðaÞ. (78) 0

This average is computed on one gray square, with the strain field (75) continued inside the void. Setting pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi rðzÞ  a2  z2 , an integral equation results: pffiffiffi pffiffiffi Z rðzÞ gB ðyÞ dy ¼ ½1  2rðzÞgB ðzÞ  gB ðaÞ; z 2 ½0; a. (79) 2 0

Its solution is obtained following Section 3.2.1. In particular, with z ¼ 0 and (78): pffiffiffi  þ ð1 þ 2 2aÞgB ðaÞ pffiffiffi gB ð0Þ ¼ . 2ð1  2aÞ

(80)

The displacement is everywhere continuous. The stress sPS is deduced from (75) and from sPS ¼ 2mPS in the matrix, with sPS ¼ 0 in the voids. The transverse components sm , sSS are computed analogously to (44), with integrations and differentiations carried out over x0 and y0 . Stress equilibrium provides: qx0 sPS ¼ qy0 ðsm  sSS Þ;

qy0 sPS ¼ qx0 ðsm þ sSS Þ.

(81)

Due to symmetry, sm ðx; xÞ ¼ sSS ðx; xÞ  0 and sm ðx; 1  xÞ ¼ sSS ðx; 1  xÞ  0, so that with 0ox; yo1 in the matrix (for convenience, a unit cell translated by a vector t ¼ ð12; 12Þ is considered):   Z x0   Z y0 1 1 0 0 sm;SS ðx0 ; y0 Þ ¼ dzq s ðz; y Þ  dzqx0 sPS ðx0 ; zÞ. (82) y PS 2 1=pffiffi2 2 0 The plus (resp. minus) sign applies to sm (resp. sSS ). Then, in the matrix: m (83) sm;SS ðx; yÞ ¼ ½ðx þ y  1Þg0 ðy  xÞ  ðy  xÞg0 ðx þ y  1Þ 2 pffiffiffi pffiffiffi pffiffiffi for 0ox; yo1. For z in ½0; a 2, the derivative g0 reduces to g0 ðzÞ ¼ g0B ðz= 2Þ= 2 where g0B is the derivative of gB . The stress and strain singularities are qualitatively the same as in the a ¼ 0, SS loading case. In particular, pffiffiffi from (79), gB ðzÞ ’ gB ðaÞ þ 2 a½gB ð0Þ þ gB ðaÞja  zj1=2 , for zta, so that sPS blows up along the ‘‘frontiers’’ of ðAÞ, ðBÞ, ðCÞ, and ðDÞ as d 1=2 where d is the distance to the ‘‘frontier’’. 4.2.2. Moments and effective modulus e mðf Þ in the dilute limit Expanding gB ðzÞ in powers of a, we obtain pffiffiffi    32 2 3 5 2 2 a þ gB ðaÞ= ¼ 1  2pa  p  6p  8 a4 þ Oða5 Þ. 3 2

(84)

ARTICLE IN PRESS F. Willot et al. / J. Mech. Phys. Solids 56 (2008) 1245–1268

Dilute expansions in f for the effective modulus and the moments ensue:   e 32 3 4 m ¼ 1  f  pffiffiffi f 3=2 þ 1   2 f 2 þ Oðf 5=2 Þ, p p m 3 2p3=2 32f 3=2 M ð1Þ ðPS Þ ¼ 1  pffiffiffi þ Oðf 2 Þ; 3=2 3 2p 4 21=4 f 3=4 þ Oðf 5=4 Þ; Sð1Þ ðPS Þ ¼ pffiffiffi 3p3=4 Sð1Þ ðSS Þ ¼ 0;

S ð1Þ ðsSS Þ ¼ 1;

32f 1=2 M ð2Þ ðPS Þ ¼ 1 þ pffiffiffi þ Oðf Þ, 3 2p3=2 4 21=4 f 3=4 Sð1Þ ðsPS Þ ¼ pffiffiffi þ Oðf 5=4 Þ, 3p3=4

S ð1Þ ðm Þ ¼ 0;

S ð1Þ ðsm Þ ¼ 1.

1265

(85a)

(85b)

(85c) (85d)

In the range 0of op=8, e mðf Þ is computed from a numerical solution of (79). Contrarily to the SS case where e lðf Þ vanishes at p=4, no ‘‘mechanically driven percolation’’ is found here, since the zones ðCÞ where bands cross are not zones of vanishing stress. Regular close-packing behavior must occur at f c ¼ p=4. 4.2.3. Distributions Field distributions are studied as in Section 3.2.4. Quite similar features are observed, summarized hereafter, along with some differences. The density PPS is non-symmetric, and supported by the interval ½gB ðaÞ; gB ð0Þ. A rescaling similar to (61), using m instead of l, provides PsPS , plotted in Fig. 7 (dashed curve). Due to the gradient of PS ðx; yÞ vanishing along the lines y ¼ x, where the field is maximum in (B), PPS ðtÞ blows up as jt  t0 j1=2 at t0 ¼ ½gB ðaÞ þ gB ð0Þ=2. A notable difference with the case a ¼ 0 lies in the existence of the extended ‘‘foot’’ on the right side of PsPS . It ends with a discontinuous VHS of finite amplitude, due to the maximum strain field value, max PS ¼ gB ð0Þ, being reached at ðx; yÞ ¼ ð12; 12Þ, in (C). In the case a ¼ 0, s ¼ 0 in this zone. The discontinuity is computed by remarking that according to (75), (76) and to a Taylor expansion of gB ðzÞ at z ¼ 0, the strain field PS at points ðx; yÞ ¼ ð12; 12Þ is a quadratic function with negative eigenvalues. Furthermore, a field hðx;pyÞ ’ ðh1 x2 þ h2 y2 Þ with h1 , h2 40, on a domain V , produces in Ph ðtÞ a ffiffiffiffiffiffiffiffiffi jump at the origin ½½Ph ð0Þ ¼ p=ð4 h1 h2 V Þ. Densities Psm , PsSS stem from Eq. (83). They blow up at the origin due to points of vanishing derivatives qx sm;SS , qy sm;SS . According to (83) this happens for ðx; yÞ ¼ ð12; 12Þ. In the vicinity of this point, the expansion of sm is quadratic of a regular type and leads to logarithmic divergent VHS of Psm ðtÞ at the origin. On the other hand, sSS ðx; yÞ is of an unusual cubic form: sSS ðx; yÞ / ðy  xÞðx þ y  1Þðx  12Þ.

(86)

This generates a singular contribution to PsSS ðtÞ at the origin, of a type not previously encountered: indeed, the distribution Ph ðtÞ of a function hðx; yÞ ¼ h0 ðy2  x2 Þx behaves near t ¼ 0 as jtj1=3 . Finally, the tails of Psm ðtÞ, PsSS ðtÞ when jtj ! 1 decay unsurprisingly as jtj3 , due to the blowing up of the concerned stress components along the ‘‘frontiers’’ between zones ðAÞ, ðBÞ, ðCÞ, ðDÞ as an inverse square root of the distance. 4.3. Equibiaxial loading We consider here the limit a ¼ l=m ¼ 1, achieved by taking m ! 0, under equibiaxial loading. The limit l ! 1 is not dealt with for the same reason as in Section 3.4. Moderate porosities f of c are assumed. The  method used pffiffifor ffi the case a ¼ 0 is applied to the void 0lattice rotated pffiffiffi by 45 ,0 and undergoing pffiffiffi symmetry (72). With as ¼ 2a and with the rotated coordinates x ¼ ðy þ xÞ= 2 and y ¼ ðy  xÞ= 2, the stresses are sPS ðx; yÞ  0, and sm;SS ðx; yÞ ¼ s½wðy0 Þ  wðx0 Þ=½2ð1  2as Þ,

(87)

in the rotated gray domain of Fig. 2(b). Hence, with $  ð=2Þ=½‘ þ ð1  ‘Þas , m ðx; yÞ ¼ ‘$½wðy0 Þ þ wðx0 Þ;

SS ðx; yÞ ¼ $½wðy0 Þ  wðx0 Þ.

(88)

ARTICLE IN PRESS F. Willot et al. / J. Mech. Phys. Solids 56 (2008) 1245–1268

1266

Likewise, the PS component of the strain field takes the Dirac-localized form for x0 ; y0 in the URQ of the rotated gray square:      1 1 PS ðx; yÞ ¼ ð1 þ ‘Þ$ x0  pffiffiffi dðy0 þ aÞ þ y0  pffiffiffi dðx0  aÞ . (89) 2 2 Then, the following effective compressibility modulus is obtained: e 1  2as k ¼ . l ‘ þ ð1  ‘Þas

(90)

‘‘Mechanically driven percolation’’ is again observed here. The moments of the fields are read from expressions (34b)–(34e), by carrying out the following substitutions: (i) replace a by as ; (ii) replace m by 1=‘; (iii) replace index PS by index m; (iv) replace index m by index SS; (v) replace index SS by index PS.

5. Discussion and conclusion Two essentially different types of solutions, in an infinitely anisotropic linear-elastic medium containing a periodic distribution of voids, have been exhibited. The anisotropy directions of the matrix—a ‘‘hard’’ one and a ‘‘soft’’ one—being aligned with the lattice directions, the fields are arranged in patterns of bands. On the one hand, for shear loading along a ‘‘hard’’ direction, the following observations were made: (i) the ‘‘parallel’’ components of the stress and strain are discontinuous, piecewise constant, with a bandwidth of one pore diameter; (ii) the deformation pattern has a ‘‘block-sliding’’ structure; (iii) the shear strain orthogonal to the loading has infinite variance, being localized as Dirac distributions along the band ‘‘frontiers’’, which are tangent to the voids (with jumps at the tangency points); (iv) the ‘‘perpendicular’’ component of the stress is zero. The presence of discontinuities in the tangential component of the displacement field is analogous to similar observation in ideal plasticity, including the Hencky plasticity model (Suquet, 1981). Moreover, in situations where the bands cross, the medium develops fictitious porous zones of zero parallel stress which lower the close-packing threshold. Also, a first-order correction f 1=2 is induced in the effective shear modulus. The latter decays linearly with the void concentration near the close-packing threshold. On the other hand, for loading along a ‘‘soft’’ direction, another type of solution is produced when the problem is hyperbolic, i.e. assumed to be incompressible. In this case, the following observations were made: (i) the strains are continuous everywhere; (ii) the leading-order correction to the effective shear modulus in the loading direction is now / f , i.e. less sensitive to increasing void concentration for small concentrations; (iii) the variances of the ‘‘parallel’’ stress and strain in the matrix are / f 3=4 and are again less sensitive to increasing void concentration; (iv) the parallel strain enjoys scaling properties near the close-packing threshold f c , operating within a band of width x aðf c  f Þ1=2 , and becomes fully localized only for f near f c . This situation is at best one of ‘‘very soft’’ void-induced strain localization. However, stresses accumulate in the ‘‘hard’’ directions (transverse and equibiaxial) and blow up along the ‘‘frontiers’’ of the field pattern as the inverse square root of the distance (the corresponding strains being zero). This particular stress localization configuration may be relevant to strain-locking materials, such as the shape-memory polycrystals that have been considered recently by Bhattacharya and Suquet (2005) and Chenchiah and Bhattacharya (2005). In conclusion, situations responsible for fractional exponents showing up in anisotropic linear theories, of relevance to non-linear homogenization methods, have been clarified. Comparisons with fully numerical results at moderate anisotropies, and for isotropic viscoplasticity, will be presented elsewhere (Willot et al., 2007). As a final remark, the fact that strong singularities are found (isolated points of matter overlapping or of matter separation) in one case is obviously of relevance for breakdown studies (e.g., in brittle materials), since these points may act as initiators. This may suggest that mechanically sounder solutions could be looked for in a large deformation framework. However, in view of our more modest aim to help better understand possible hallmarks of localization in non-linear homogenization theories, the model considered here is adequate.

ARTICLE IN PRESS F. Willot et al. / J. Mech. Phys. Solids 56 (2008) 1245–1268

1267

Acknowledgements The authors wish to thank the anonymous referee for remarks pertaining to the interpretation of the solutions in terms of characteristics. The work of P.P.C. was supported in part by the N.S.F. through Grant OISE-02-31867 in the context of collaborative project with the C.N.R.S. Appendix A. Symmetry properties for the strain, stress, and displacement fields Since the periodic anisotropic medium obeys the symmetries of the square, the symmetries of the displacement u are determined by the applied loading. They are easily deduced from Figs. 1(b) in SS and 1(c) in PS, and from an obvious flow pattern in equibiaxial loading. Differentiations wrt. x and y provide those enjoyed by ij . The symmetry group of the constitutive law carries them unchanged over sij . For v ¼ u or u and for a ¼ e or r:





In PS loading: vx ðx; yÞ ¼ vx ðx; yÞ ¼ vx ðx; yÞ, vy ðx; yÞ ¼ vy ðx; yÞ ¼ vy ðx; yÞ,

ðA1aÞ ðA1bÞ

vx ðx; yÞ ¼ vy ðy; xÞ, aii ðx; yÞ ¼ aii ðx; yÞ ¼ aii ðx; yÞ;

ðA1cÞ ðA1dÞ

i ¼ x; y,

axx ðx; yÞ ¼ ayy ðy; xÞ,

ðA1eÞ

axy ðx; yÞ ¼ axy ðx; yÞ ¼ axy ðx; yÞ, axy ðx; yÞ ¼ axy ðy; xÞ;

ðA1fÞ ðA1gÞ

In SS loading: vx ðx; yÞ ¼ vx ðx; yÞ ¼ vx ðx; yÞ,

ðA2aÞ

vy ðx; yÞ ¼ vy ðx; yÞ ¼ vy ðx; yÞ, vx ðx; yÞ ¼ vy ðy; xÞ,

ðA2bÞ ðA2cÞ

aii ðx; yÞ ¼ aii ðx; yÞ ¼ aii ðx; yÞ;



i ¼ x; y,

ðA2dÞ

axx ðx; yÞ ¼ ayy ðy; xÞ, axy ðx; yÞ ¼ axy ðx; yÞ ¼ axy ðx; yÞ,

ðA2eÞ ðA2fÞ

axy ðx; yÞ ¼ axy ðy; xÞ;

ðA2gÞ

In equibiaxial loading: vx ðx; yÞ ¼ vx ðx; yÞ ¼ vx ðx; yÞ, vy ðx; yÞ ¼ vy ðx; yÞ ¼ vy ðx; yÞ,

ðA3aÞ ðA3bÞ

vx ðx; yÞ ¼ vy ðy; xÞ, aii ðx; yÞ ¼ aii ðx; yÞ ¼ aii ðx; yÞ;

ðA3cÞ ðA3dÞ

i ¼ x; y,

axx ðx; yÞ ¼ ayy ðy; xÞ,

ðA3eÞ

axy ðx; yÞ ¼ axy ðx; yÞ ¼ axy ðx; yÞ, axy ðx; yÞ ¼ axy ðy; xÞ.

ðA3fÞ ðA3gÞ

References Abrikosov, A.A., Campuzano, J.C., Gofron, K., 1993. Experimentally observed extended saddle point singularity in the energy spectrum of YBa2 Cu4 O6 :9 and YBa2 Cu4 O8 and some of the consequences. Physica C: Superconductivity 214 (1–2), 73–79.

ARTICLE IN PRESS 1268

F. Willot et al. / J. Mech. Phys. Solids 56 (2008) 1245–1268

Bhattacharya, K., Suquet, P.M., 2005. A model problem concerning recoverable strains of shape memory polycrystals. Proc. R. Soc. London A 461, 2797–2816. Chenchiah, I.V., Bhattacharya, K., 2005. Examples of non-linear homogenization involving degenerate energies. I. Plane strain. Proc. R. Soc. London A 461, 3681–3703. Cule, D., Torquato, S., 1998. Electric field distribution in composite media. Phys. Rev. B 58, 11829–11832. Day, A.R., Snyder, K.A., Garboczi, E.J., Thorpe, M.F., 1992. The elastic moduli of a sheet containing circular holes. J. Mech. Phys. Solids 40, 1031–1051. Drucker, D.C., 1966. The continuum theory of plasticity on the macroscale and the microscale. J. Mater. 1, 873–910. Ekeland, I., Temam, R., 1974. Analyse convexe et proble`mes variationnels. Dunod, Paris. Francescato, P., Pastor, J., 1998. Re´sistance de plaques multiperfore´es: comparaison calcul-experience. Rev. Eur. E´le´ments Finis 7, 421–437. Francescato, P., Pastor, J., Riveill-Reydet, B., 2004. Ductile failure of cylindrically porous materials. Part I: plane stress problem and experimental results. Eur. J. Mech. A/Solids 23, 181–190. Kachanov, L.M., 1974. Fundamentals of the Theory of Plasticity. Mir Publishers, Moscow. Lahellec, N., Suquet, P., 2004. Nonlinear composites: a linearization procedure exact to second order in contrast and for which the strainenergy and affine formulations coincide. C. R. Me´canique 332, 693–700. Masson, R., Bornert, M., Suquet, P., Zaoui, A., 2000. An affine formulation for the prediction of the effective properties of non-linear composites and polycrystals. J. Mech. Phys. Solids 48, 1203–1227. McPhedran, R.C., Movchan, A.B., 1994. The Rayleigh multipole method for linear elasticity. J. Mech. Phys. Solids 42, 711–727. Otto, M., Bouchaud, J.-P., Claudin, P., Socolar, J.E.S., 2003. Anisotropy in granular media: classical elasticity and directed-force chain network. Phys. Rev. E 67, 031302. Pastor, J., Ponte Castan˜eda, P., 2002. Yield criteria for porous media in plane strain: second-order estimates versus numerical results. C. R. Me´canique 330, 741–747. Pellegrini, Y.P., 2001. Self-consistent effective-medium approximation for strongly nonlinear media. Phys. Rev. B 64, 134211. Ponte Castan˜eda, P., 1991. The effective mechanical properties of nonlinear isotropic composites. J. Mech. Phys. Solids 39 (1), 45–71. Ponte Castan˜eda, P., 1996. Exact second-order estimates for the effective mechanical properties of nonlinear composite materials. J. Mech. Phys. Solids 44, 827–862. Ponte Castan˜eda, P., 2001. Second-order theory for nonlinear dielectric composites incorporating field fluctuations. Phys. Rev. B 64, 214205. Ponte Castan˜eda, P., 2002. Second-order homogenization estimates for non-linear composites incorporating field fluctuations. I Theory. J. Mech. Phys. Solids 50, 737–757, and Part II—Applications, J. Mech. Phys. Solids 50, 759–782. Ponte Castan˜eda, P., DeBotton, G., Li, G., 1992. Effective properties of non-linear inhomogeneous dielectrics. Phys. Rev. B 46 (8), 4387–4394. Suquet, P.M., 1981. Sur les e´quations de la plasticite´: existence et re´gularite´ des solutions. J. Me´canique 20, 3–39. Torquato, S., 2002. Random Heterogeneous Materials: Microstructure and Macroscopic Properties. Springer, New York. Van Hove, L., 1953. The occurrence of singularities in the elastic frequency distribution of a crystal. Phys. Rev. 89 (6), 1189–1193. Willis, J.R., 1986. Variational estimates for the overall response of an inhomogeneous non-linear dielectric. In: Ericksen, J.L., et al. (Eds.), Homogenization and Effective Moduli Materials and Media. Springer, New York, pp. 247–263. Willot, F., Pellegrini, Y.-P., Idiart, M., Ponte Castan˜eda, P., 2007. In preparation. Zachmanoglou, E.C., Thoe, D.W., 1986. Introduction to Partial Differential Equations. Dover, New York. Zeng, X.C., Bergman, D.J., Hui, P.M., Stroud, D., 1988. Effective-medium theory for weakly non-linear composites. Phys. Rev. B 38, 10970–10973.