On the quasistatic effective elastic moduli for elastic waves in three-dimensional phononic crystals

Journal of the Mechanics and Physics of Solids 61 (2013) 2260–2272

Contents lists available at ScienceDirect

Journal of the Mechanics and Physics of Solids journal homepage: www.elsevier.com/locate/jmps

On the quasistatic effective elastic moduli for elastic waves in three-dimensional phononic crystals A.A. Kutsenko a, A.L. Shuvalov a, A.N. Norris b,n a b

Institut de Mécanique et d'Ingénierie de Bordeaux, Université de Bordeaux, UMR CNRS 5469, Talence 33405, France Mechanical and Aerospace Engineering, Rutgers University, Piscataway, NJ 08854-8058, USA

a r t i c l e in f o

abstract

Article history: Received 16 December 2012 Received in revised form 19 June 2013 Accepted 24 June 2013 Available online 9 July 2013

Effective elastic moduli for 3D solid–solid phononic crystals of arbitrary anisotropy and oblique lattice structure are formulated analytically using the plane-wave expansion (PWE) method and the recently proposed monodromy-matrix (MM) method. The latter approach employs Fourier series in two dimensions with direct numerical integration along the third direction. As a result, the MM method converges much quicker to the exact moduli in comparison with the PWE as the number of Fourier coefficients increases. The MM method yields a more explicit formula than previous results, enabling a closed-form upper bound on the effective Christoffel tensor. The MM approach significantly improves the efficiency and accuracy of evaluating effective wave speeds for high-contrast composites and for configurations of closely spaced inclusions, as demonstrated by threedimensional examples. & 2013 Elsevier Ltd. All rights reserved.

Keywords: Effective moduli Plane-wave expansion Homogenization Monodromy matrix

1. Introduction Long-wave low-frequency dispersion of acoustic waves in periodic structures is of both fundamental and practical interest, particularly due to the current advances in manufacturing of metamaterials and phononic crystals. In this light, the leading order dispersion (quasistatic limit) has recently been under intensive study by various theoretical approaches such as plane-wave expansion (PWE) (Krokhin et al., 2003; Ni and Cheng, 2005; Kutsenko et al., 2011a), scaling technique (Parnell and Abrahams, 2008; Andrianov et al., 2011), asymptotics of multiple-scattering theory (Mei et al., 2007; Torrent and Sánchez-Dehesa, 2008; Wu and Zhan, 2009) and a newly proposed monodromy-matrix (MM) method (Kutsenko et al., 2011b). The cases treated were mostly confined to scalar waves in 2D (two-dimensional) structures. Regarding vector waves in 2D and especially in 3D phononic crystals, a variety of methods have been proposed for calculating the quasistatic effective elastic properties of 3D periodic composites containing spherical inclusions arranged in a simple cubic array. For the case of rigid inclusions an integral equation on the sphere surface was solved numerically to obtain the effective properties (Nunan and Keller, 1984). Spherical voids (Nemat-Nasser and Taya, 1981) and subsequently elastic inclusions were considered using a Fourier series approach (Nemat-Nasser et al., 1982). Alternative procedures for elastic spherical inclusions include the method of singular distributions (Sangani and Lu, 1987) and infinite series of periodic multipole solutions of the equilibrium equations (Kushch, 1987). The latter multipole expansion method has also been applied to cubic arrays of ellipsoidal inclusions (Kushch, 1997). A particular PWE-based method of calculation of quasistatic speeds in 3D

n

Corresponding author. E-mail addresses: [email protected] (A.A. Kutsenko), [email protected] (A.L. Shuvalov), [email protected] (A.N. Norris).

0022-5096/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jmps.2013.06.003

A.A. Kutsenko et al. / J. Mech. Phys. Solids 61 (2013) 2260–2272

2261

phononic crystals of cubic symmetry has been formulated and implemented in Ni and Cheng (2007). A review of numerical methods for calculating effective properties of composites can be found in Milton (2001, Sections 2.8, 14.11). Of all the methods available for calculating effective elastic moduli the PWE method is arguably the simplest and most straightforward to implement. It requires only Fourier coefficients of the inclusion in the unit cell, which makes it the method of choice for many problems. Unfortunately, PWE is not a very practical tool for the 3D case, where the vectors and matrices in the Fourier space are of very large algebraic dimension, especially if the phononic crystal is composed of highly contrasting materials (examples in Section 5 illustrate this critical drawback). The present paper provides the PWE and MM analytical formulations of the 21 components of the effective elastic stiffness for 3D solid–solid phononic crystals of arbitrary anisotropy and arbitrary oblique lattice. While the PWE method is widely used in some fields its formulation for general anisotropic static elasticity has, surprisingly, not been discussed before. The PWE is presented here in a compact form (see Eq. (15)) suitable for numerical implementation. The main thrust of the paper is concerned with the MM approach. The motivation for advocating this method as an alternative to the more conventional PWE technique is, first, that the MM method ‘spares’ Fourier expansion in one of the coordinates (this is particularly advantageous for the 3D numerics) and, second, that the MM method has much faster convergence than PWE. Comparison of the MM and PWE calculations provided in the paper confirms a markedly better efficiency of the MM method. The paper is organized as follows. In Section 2, the quasistatic perturbation theory is used to define the effective Christoffel equation in the form which serves as the common starting point for the PWE and MM methods. The PWE formulation of the effective elastic moduli follows readily and is also presented in Section 2. The MM formulation is described in Section 3: the derivation of the MM formula is in Section 3.1 (see also Appendix A), its numerical implementation is discussed in Section 3.2, generalization to the case of an oblique lattice is presented in Section 3.3 and the scheme for recovering the full set of effective elastic moduli is provided in Section 3.4. A closed-form estimate of the effective Christoffel matrix is presented in Section 4. Examples of the MM and PWE calculations are provided in Section 5. Concluding remarks are given in Section 6. 2. Background. PWE formula Consider a 3D anisotropic medium with density and elastic stiffness ρðxÞ ¼ ρðx þ ep Þ;

CðxÞ ¼ Cðx þ ep Þ;

ð1Þ

which are assumed to be 1-periodic, i.e. invariant to period or translation vectors ep ¼ ðδpq Þ (a cubic lattice, otherwise see Section 3.3). All roman indices run from 1 to 3. In the following, n and + mean complex and Hermitian conjugation respectively. Assume no dissipation so that the elements of C satisfy cijkl ¼ cnklij ð ¼ cklij for real case). Our goal is the quasistatic eff effective elastic stiffness Ceff with elements ceff ijkl that has the same symmetries as those of C, and matrices Cjl defined by analogy with Eq. (2). For compact writing, introduce the matrices Cjl ¼ ðcijkl Þ3i;k ¼ 1 ¼ Cþ lj

ð2Þ iωt

with components numbered by i; k. The elastodynamic equation for time-harmonic waves vðx; tÞ ¼ vðxÞe 2

∂j ðCjl ∂l vÞ ¼ ρω v;

is ð3Þ

where ∂j ≡∂=∂xj and repeated indices are summed. The differential operator in Eq. (3) is self-adjoint with respect to the Floquet condition vðxÞ ¼ uðxÞeik⋅x with 1-periodic uðxÞ ¼ uðx þ ep Þ and k ¼ kκ ðjκj ¼ 1Þ. Substituting this condition in (3) casts it into the form 2

ðC0 þ kC1 þ k C2 Þu ¼ ρω2 u

where

C0 u≡∂j ðCjl ∂l uÞ;

C1 u≡iðκ j Cjl ∂l þ κ l ∂j Cjl Þu;

C2 u≡κ j κ l Cjl u:

ð4Þ

All operators C are self-adjoint. We introduce for future use the linear, areal and volumetric averages over the unit-cell: 〈  〉j is the average over coordinate xj, 〈  〉j is the average over the section orthogonal to xj, and 〈  〉 is the complete average. These averages in turn define inner products of vector-valued functions. Thus, for a scalar function f and vector-functions f; h Z 1 Z þ þ 〈f 〉j ¼ f dxj ; 〈f 〉 ¼ f dx; ðf; hÞj ¼ 〈h f〉j ; ðf; hÞ ¼ 〈h f〉: ð5Þ ½0;13

0

Next we apply perturbation theory to (4) with a view to define the quasistatic effective Christoffel matrix whose eigenvalues yield the effective speeds cα ≡cα ðκÞ≡lim ωα ðkÞ=k; k-0

α ¼ 1; 2; 3:

ð6Þ

For k ¼ 0, the eigenvalue ω ¼ 0 has multiplicity 3 and corresponds to three constant linear independent eigenvectors u0α . Consider the asymptotics 2

3

ω2α ¼ 0 þ kλ1α þ k λ2α þ Oðk Þ; 2

ð7aÞ 3

uα ¼ u0α þ ku1α þ k u2α þ Oðk Þ:

ð7bÞ

2262

A.A. Kutsenko et al. / J. Mech. Phys. Solids 61 (2013) 2260–2272

Substituting (7) into (4) and collecting terms with the same power of k yield 1 : C0 u0α ¼ 0;

ð8aÞ

k : C1 u0α þ C0 u1α ¼ ρλ1α u0α ;

ð8bÞ

2

k : C2 u0α þ C1 u1α þ C0 u2α ¼ ρλ1α u1α þ ρλ2α u0α :

ð8cÞ

Scalar multiplying (8b) by u0α and (8c) by ek ; and using ðC1 u0α ; u0α Þ ¼ 0 together with self-adjointness of C0 lead to u1α ¼ C1 0 C 1 u0α ;

λ1α ¼ 0; where

λ2α ¼ c2α

α ¼ 1; 2; 3;

ð10Þ

3 Γ ¼ ððC2 ek ; ei ÞðC1 0 C 1 ek ; C 1 ei ÞÞi;k ¼ 1 . Γ ¼ ðΓ ik Þ3i;k ¼ 1 in the form

ΓðκÞ ¼ Cejl κj κ l The matrix

Cejl

ΓðκÞ ¼

ð9Þ

due to λ1α ¼ 0 and (6). Inserting u1α from (9)2 in (9)3 gives

Γu0α ¼ 〈ρ〉c2α u0α ; where matrix

ðC2 u0α ; ek Þ þ ðC1 u1α ; ek Þ ¼ 〈ρ〉λ2α ðu0α ; ek Þ; k ¼ 1; 2; 3;

Substituting C1 ; C2 from (4) defines the quasistatic effective 3  3 Christoffel þ

e 1 where Cejl ≡〈Cjl 〉〈ð∂p Cþ pj ÞC 0 ð∂q Cql Þ〉ð ¼ Clj Þ:

Ceff jl ,

is distinguished from   eff κj κl : ¼ 12 Ceff jl þ Clj

ð11Þ

in terms of which the Christoffel matrix is

Ceff jl κ j κ l

ð12Þ

Comparison of Eqs. (11)1 and (12) implies that Cejl and Ceff jl are related by eff e e Ceff jl þ Clj ¼ Cjl þ Clj

ð13Þ

and they are equal if j¼l. Eq. (11)2 does not yield Ceff jl explicitly, only in the combination (13); however, this connection is sufficient for the purpose of finding all elements of Ceff , as described in Section 3.4. For now we focus on methods to calculate Cejl . Eq. (11)2 is still an implicit formula for the matrix Cejl insofar as the operator C1 0 is not specified. One way to an explicit implementation of (11)2 is via its transformation to Fourier space which must be truncated to define C1 0 as a matrix inverse. Thus taking the 3D Fourier expansion b ðgÞe2πigx ; Cjl ðxÞ ¼ ∑ C jl

b ðgÞ ¼ 〈C ðxÞe2πigx 〉 C jl jl

ð14Þ

g∈Z3

and plugging it into (11)2 yields the PWE formula for Cejl as follows: b ð0Þqþ C1 q where q ¼ ðg C b Cejl ¼ C jl l j i ij ðgÞÞg∈Z3 \0 ; j 0

b ðgg′ÞÞ C0 ¼ ðg k g′p C kp g;g′∈Z3 \0 :

ð15Þ

This result corresponds to the quasistatic limit of the effective elastic coefficients for a dynamically homogenized periodic medium, see Norris et al. (2012, Eq. (2.11)). Unfortunately, although the PWE formula (15) is straightforward, its numerical use in 3D is complicated by the large algebraic dimensions of vectors and matrices of Fourier coefficients. This motivates using the monodromy matrix (MM) method which confines the PWE computation to 2D, see next section. The MM method will also be shown to have significantly faster convergence than PWE. 3. MM formula 3.1. Derivation The idea underlying the MM method is to reduce the three-dimensional problem of Eq. (11)2 to an equivalent onedimensional equation that can be integrated. This is achieved by focusing on a single coordinate and using Fourier transforms in the orthogonal coordinates. The MM formula may be deduced proceeding from Eq. (11)2; using integration by parts, let us recast it into the form Cejl ¼ 〈Cjl 〉〈ð∂p Cþ pj ÞU〉 ¼ 〈Cjp ∂p V〉;

ð16Þ

where we have denoted C1 0 ð∂q Cql Þ ¼ U and V ¼ U þ xl I3 with I3 standing for the 3  3 identity matrix. Thus we need to solve the equation C0 U ¼ ∂q Cql ⇔ C0 U ¼ C0 xl I3 ⇔ C0 V ¼ 0 ⇔ ∂q ðCqp ∂p VÞ ¼ 0

with V ¼ U þ xl I3 :

ð17Þ

In the following derivation the indices j; l are regarded as fixed, all repeated indices are summed and among them the indices a; b are specialized by the condition a; b≠l. The suffix 0 of the differential operators Q0 ; M0 below indicates their reference to ω; k ¼ 0 (similarly to C0 ).

A.A. Kutsenko et al. / J. Mech. Phys. Solids 61 (2013) 2260–2272

Eq. (17) can be rewritten as an ordinary differential equation in the designated coordinate xl ! ! V B1 B1 A1 Ξ′ ¼ Q0 Ξ with ′≡∂l ; Ξ ¼ ; Q0 ¼ 1 1 ; Clp ∂p V A2 Aþ A1 Aþ 1B 1B

2263

ð18Þ

where Ξ is a 6  3 matrix and the matrix operators B, A1 and A2 are defined by BV ¼ Cll V;

A1 V ¼ Clb ∂b V;

A2 V ¼ ∂a ðCab ∂b VÞ:

ð19Þ

Note that B and A2 are self-adjoint at any fixed xl with respect to the inner product ð; Þl (see (5)). Denote ΞðsÞ≡ΞðxÞjxl ¼ s . The solution to (18) with some initial matrix function Ξð0Þ has the form Zd xl ðI þ Q0 dxl Þ; ð20Þ Ξðxl Þ ¼ M0 ðxl ÞΞð0Þ with M0 ðxl Þ ¼ 0

R where M0 is a propagator matrix defined through the multiplicative integral b and I is the identity operator. Note the identities !   0 I3 ~ ~ ~ ~ for W0 ¼ Mþ ð21Þ : Q0 W0 ¼ 0; Qþ ; W0 ¼ 0 W 0 ¼ 0 ⇒ M0 W0 ¼ W 0 ; 0 W0 ¼ W0: I3 0 It is seen that for any value of xl the operator M0 I has no inverse but at the same time it is a one-to-one mapping from ~ 0 : By the definitions of V and Ξ in (17), (18) and due to the subspace orthogonal to W0 onto the subspace orthogonal to W periodicity of U, it follows that Ξð1Þ ¼ Ξð0Þ þ W0 ⇒ Ξð0Þ ¼ ðM0 ð1ÞI Þ1 W0 ≡S:

ð22Þ

Thus from (20), (22) and (18), (21) ~ þ M0 ðxl ÞS: Clp ∂p V ¼ W 0

Ξðxl Þ ¼ M0 ðxl ÞS ⇒ V ¼ Wþ 0 M0 ðxl ÞS;

Substituting V obtained from (23) into (16) yields the desired formula for Note that ~ 0 Þ ¼ ðS; Mþ ðxl ÞW ~ 0 Þ ¼ ðS; W ~ 0 Þ ¼ 〈W ~ þ S〉 : 〈Clp ∂p V〉l ¼ ðM0 ðxl ÞS; W 0 0 l l l l

ð23Þ Cejl ,

which is discussed further in Section 3.2.3. ð24Þ

The latter identity implies that the laterally averaged ‘traction’ component of Ξðxl Þ is independent of xl, and may be identified as a net ‘force’ acting on the faces xl ¼constant. e Further simplification can be achieved for the matrices Ceff ll ð ¼ Cll Þ. By (16) and (58) þ

þ

1 ~ ~ Ceff ll ¼ 〈Clp ∂p V〉 ¼ 〈W 0 S〉 ¼ 〈W 0 ðM0 ð1ÞI Þ W 0 〉l :

ð25Þ

Eq. (25) suffices to define the effective Christoffel tensor for k ¼ kκ parallel to the translation vector el in which case ΓðκÞ ¼ Ceff ll . Note that an alternative derivation of (25) is given in Appendix A. Finally it is noted that M0 ð1Þ which appears in the above expressions is formally a monodromy matrix (MM) relatively to the coordinate xl, for which reason this approach and its outcome formulas are referred to as the MM ones. The MM approach to finding effective speed of shear (scalar) waves in 2D structures was first presented in Kutsenko et al. (2011b), where Eq. (25) for vector waves was also mentioned but with neither derivation nor discussion. 3.2. Implementation 3.2.1. Propagator matrix in direction el As above, we fix l and keep a; b≠l. Introduce the 2D Fourier expansion b pq ðg; x Þ; Cpq ðxÞ ¼ ∑ e2πiga xa C l g∈Z2

b pq ðg; x Þ ¼ 〈Cpq ðxÞe2πiga xa 〉 : C l l

ð26Þ

Operators B; A1 ; A2 and matrix operators Q; M0 ðxl Þ defined in (18)–(20) are represented in the 2D Fourier space by the following infinite matrices truncated to a finite size in calculations: b ðgg′; x ÞÞ B↦B ¼ ðC ll l g;g′∈Z2 ;

ð27aÞ

b ðgg′; x Þg′ Þ A1 ↦A1 ¼ 2πiðC lb l b g;g′∈Z2 ;

ð27bÞ

b ðgg′; x Þg′ Þ A2 ↦A2 ¼ 4π 2 ðg a C ab l b g;g′∈Z2 ;

ð27cÞ

Q0 ↦Q 0 ¼

B1 A1

B1

1 A2 Aþ 1 B A1

1 Aþ 1B

! ;

ð27dÞ

2264

A.A. Kutsenko et al. / J. Mech. Phys. Solids 61 (2013) 2260–2272

M0 ðxl Þ↦M0 ðxl Þ ¼

Zd xl 0

ðI þ Q 0 dxl Þ:

ð27eÞ

3.2.2. Principal directions κ∥el eff Consider the effective matrices Ceff ll ð ¼ Γ if κ∥el Þ. In view of the above notations, formula (25) for Cll is expressed as þ

b ~ S Ceff ll ¼ W b 0

b ¼ ðM0 ð1ÞIÞ1 W ; with S b

ð28Þ

0

where Wb ¼ 0

! Eb 0 ; 0

~ ¼ W b 0

0 Eb

! with Eb ¼ ðδg0 I3 Þg∈Z2 : 0

0

ð29Þ

b in (28) can be performed by means of calculating M0 from its definition (27e), using any of the known Calculation of S methods for evaluating multiplicative integrals. However, this approach may suffer from numerical instabilities for M0 of large size due to exponential growth of some components of M0 . Another method rests on calculating the resolvent ðM0 αIÞ1 without calculating M0 . The advantage of doing so is due to the fact that growing dimension of M0 leads to a relatively moderate growth of the resolvent components. Let us consider this latter method in some detail. Denote the resolvent Rα ðxl Þ ¼ ðM0 ðxl ÞαIÞ1 where α is not an eigenvalue of M0 ðxl Þ. Since the matrix M0 ðxl Þ satisfies the differential equation M′0 ðxl Þ ¼ Q 0 ðxl ÞM0 ðxl Þ with the initial condition M0 ð0Þ ¼ I, it follows that Rα ðxl Þ with a randomly chosen α≠1 satisfies a Riccati equation with initial condition as follows: ( R′α ¼ Rα Q 0 ðI þ αRα Þ; ð30Þ Rα ð0Þ ¼ ð1αÞ1 I; where ′ ¼ ∂l . Since eigenvalues of M0 usually tend to lie close to the real axis and unit circle, it is recommended to take α∈C far from these sets. Integrating Eq. (30) numerically (we used the Runge–Kutta method of fourth order) provides R α ð1Þ ¼ ðM0 ð1ÞαIÞ1 where α≠1: To exploit it for finding Ceff ll given by (28) we use the identity 1 1 b Rα ð1ÞWb ¼ ðI þ ðα1ÞR α ð1ÞÞ1 S≡ðM 0 ð1ÞIÞ Wb ¼ ðI þ ðα1ÞR α ð1ÞÞ 0

0

Wb 0 : 1α

ð31Þ

b is found from the linear system Thus S b¼W : ð1αÞðI þ ðα1ÞR α ð1ÞÞS b

ð32Þ

0

~ þ T ¼ 0 with W ; W ~ from (29) and To solve (32), we first note that the matrix T≡ð1αÞðI þ ðα1ÞR α ð1ÞÞ satisfies TWb ¼ 0; W b b b 0 0 0 0 thus has three zero columns on the left of its vertical midline and three zero rows below the horizontal midline. Removing   these columns and rows yields the reduced matrix T~ ¼ T1 T2 where the square block T3 has three rows and three columns T3 T4

less than the block T2 . Since T3 is numerically large, while T1 ; T4 are medium and T2 is small, it is convenient to apply Schur's formula to arrive at the final relation in the form þ 1 1 : Ceff ll ¼ Eb ðT2 T1 T3 T4 Þ Eb 0 0

ð33Þ

3.2.3. Off-principal directions Consider k ¼ kκ of arbitrary orientation relative to the translations el . In order to find the effective Christoffel tensor ΓðκÞ of Eq. (12) for arbitrary κð∦el Þ, we need to calculate Cejl with j≠l. Applying the 2D Fourier expansion to Eqs. (16) and (23) yields ! b þ b b NÞ〉 b l where V ¼ M0 ðxl ÞS b l ¼ Ceff þ 〈Eþ ðC ll C jl ÞB1 ðA1 V Cejl ¼ 〈Eb ; ðA1 þ C jl ∂l ÞV〉 ð34Þ ll b b 0 0 N b jp ðgg′; x ÞÞ b and C jp ¼ ðC l g;g′∈Z2 with g ¼ ðg a g b Þ and a; b≠l. As detailed above, the matrix S can be calculated with a very good b l Þ and Nðx b l Þ is no longer that accurate, particularly for high-contrast structures which require precision. Evaluation of Vðx many terms in the Fourier series and hence need M0 of large algebraic dimension and therefore with large values of some components. Thus a negligible error in S may become noticeable after multiplying by M0 . In this regard, we present an alternative method of calculating Cejl with j≠l which circumvents (34). For the fixed functions ρðxÞ and cijkl ðxÞ with a given cubic lattice of periods ep ¼ ðδpq Þ, introduce the new periods ap ¼ Aep where A has integer components aij. The solutions vðxÞ ¼ uðxÞeik⋅x and ωðkÞ of the wave equation (3) remain unaltered, as they do not depend on the choice of periods. But now in order to define Ceff jl by the MM formula (which requires periods to coincide with base vectors) we should apply the change of variables x-Ax that leads, as explained in Section 3.3, to new

A.A. Kutsenko et al. / J. Mech. Phys. Solids 61 (2013) 2260–2272

2265

functions for the density and elasticity with the periods ep ¼ ðδpq Þ. According to Eqs. (49) and (56) of Section 3.3 eff þ eff ~ eff ~ eff Ceff pq ¼ apj C jl aql ¼ apj AC jl A aql ð⇔C jl ¼ bjp Cpq blq ; eff where C~ jl

eff

eff

þ C jl ¼ bjp BCeff pq B blq Þ

ð35Þ

are the effective matrices associated with the new profiles c~ ijkl ðxÞ ¼ bjp blq cipkq ðAxÞ and

and C jl

c ijkl ðxÞ ¼ bim bjp bkn blq cmpnq ðAxÞ, respectively, and bij stand of the transformations ap ¼ Aj ep , j ¼ 1; 2; 3 0 1 0 a1 ¼ e1 ; 1 0 0 1 B C B A2 ¼ @ 0 A1 ¼ @ 0 1 1 A⇔ a2 ¼ e2 þ e3 ; a3 ¼ e2 þ e3 ; 0 1 1 1

for components of A1 . We may apply (35) successively for each 1

0

1

1

C 0 A⇔…;

0

1

0

1

B A3 ¼ @ 1 0

1

0

1

1

C 0 A⇔…:

0

1

ð36Þ

Using the inverse forms in (35) leads to a variety of identities, for instance eff þ eff eff eff eff eff ~ eff Ceff 23 þ C32 ¼ 4ðC 22 ÞA1 C22 C33 ¼ 4ðAC 22 A ÞA1 C22 C33 ; eff þ eff eff eff eff eff ~ eff Ceff 31 þ C13 ¼ 4ðC 33 ÞA2 C33 C11 ¼ 4ðAC 33 A ÞA2 C33 C11 ; eff þ eff eff eff eff eff ~ eff Ceff 12 þ C21 ¼ 4ðC 11 ÞA3 C11 C22 ¼ 4ðAC 11 A ÞA3 C11 þ C22 :

ð37Þ

eff Inserting (37) in (12) eliminates Ceff jl þ Clj with j≠l and expresses the effective Christoffel tensor Γ fully in terms of the

~ ~ matrices Ceff ll and either of ðC ll ÞAn or ðC ll ÞAn n ¼1,2,3), which are defined by Eqs. (25), (28) with cijkl ðxÞ replaced by c ijkl ðxÞ or c ijkl ðxÞ. Thus all calculations have been reduced to the form (28) which ensures a numerically stable evaluation of Γ. Note that the transformed elasticity c ijkl ðxÞ retains the usual symmetries under the interchange of indices, while c~ ijkl ðxÞ does not (see Section 3.3). Also, the transformations defined by (36) for the cubic unit cell can be identified as rotations by virtue of the fact that p1ffiffi2Aj , j ¼ 1; 2; 3, are orthogonal matrices of unit determinant. Hence, apart from a factor of 14, ðc ijkl ðxÞÞAj is eff

eff

precisely the elasticity tensor represented in a coordinate system rotated about the axis ej by π=4 from the original. 3.3. The case of an oblique lattice 3.3.1. Equivalent problem on a cubic lattice Consider the problem of quasistatic asymptotics of the wave equation (3) for the general case of a 3D periodic elastic medium with ρðxÞ ¼ ρðx þ ap Þ;

cijkl ðxÞ ¼ cijkl ðx þ ap Þ;

ð38Þ

where the translation vectors ap form an oblique lattice. We will define the solution of this problem via the solution for a simpler case of a cubic lattice. The oblique lattice vectors are defined by a matrix A ð≠IÞ as ap ¼ Aep with A ¼ ðapq Þ3p;q ¼ 1 ¼ ðap  eq Þ3p;q ¼ 1 ;

B≡A1 ¼ ðbpq Þ3p;q ¼ 1

ð39Þ

where ep are the orthonormal vectors used previously. Define the new or transformed position variable x′ ¼ Bxð⇔x ¼ Ax′Þ, ~ ~ ðx′Þ ¼ cijkl ðxÞ, which are seen to be periodic the associated displacement vðx′Þ ¼ vðxÞ and material parameters ρðx′Þ ¼ ρðxÞ, cð1Þ ijkl ¼ ðcð1Þ Þ3 , the equation of motion (3) becomes in x′ with respect to the vectors ep . Setting Cð1Þ jl ijkl i;k ¼ 1 ~ ~ bjp blq ∂j′ ðCð1Þ ~ 2 v; pq ∂l′ vÞ ¼ ρω

ð40Þ

where ∂j′ ≡∂=∂xj ′. Using the fact that B is constant allows it to be removed explicitly from (40) by incorporation into a newly defined stiffness tensor. Thus, replacing x′-x we have ~ ¼ ρω ~ ~ 2 v; ∂j ðC~ jl ∂l vÞ

ð41Þ

~ where vðxÞ ¼ vðAxÞ and the material parameters are ~ ~ þ ep ÞÞ; ρðxÞ ¼ ρðAxÞ ð ¼ ρðx c~ ijkl ðxÞ ¼ bjp blq cipkq ðAxÞ ð ¼ c~ ijkl ðx þ ep ÞÞ; þ C~ jl ðxÞ ¼ ðc~ ijkl Þ3i;k ¼ 1 ¼ bjp blq Cpq ðAxÞ ¼ C~ lj ðxÞ;

ð42Þ

which are periodic in x with respect to the cubic lattice formed by vectors ep . Note that the tensor c~ ijkl for A≠I is of Cosserat type in that it is not invariant to permutations of indices i⇄j and k⇄l but retains the major symmetry c~ ijkl ¼ c~ nklij ð ¼ c~ klij for real case). The Floquet condition vðxÞ ¼ eikx uðxÞ with periodic uðxÞ ¼ uðx þ aj Þ satisfying ð∂l þ ikl ÞClq ð∂q þ ikq Þu ¼ ρω2 u

ð43Þ

2266

A.A. Kutsenko et al. / J. Mech. Phys. Solids 61 (2013) 2260–2272

~ ~ ~ ~ with periodic uðxÞ satisfying the equation that follows from (40) is equivalent to the condition vðxÞ ¼ eikx uðxÞ

~ ð∂l þ ik~ l ÞC~ lq ð∂q þ ik~ q Þu~ ¼ ρ~ ω~ 2 u;

ð44Þ

where ~ jκj ~ ¼ 1Þ; k~ ¼ Aþ kð ¼ k~ κ;

~ ¼ ωðkÞ; ωð ~ kÞ

~ ~ þ ep ÞÞ: uðxÞ ¼ uðAxÞð ¼ uðx

ð45Þ

Eq. (44), which is defined on a cubic lattice, has quasistatic asymptotics as described above. According to (10) ~ c~ 2α u~ 0α Γ~ u~ 0α ¼ 〈ρ〉

~ k; ~ α ¼ 1; 2; 3; ~ with c~ α ðκÞ≡lim ω~ α ðkÞ=

ð46Þ

~ k-0

eff ~ κÞ ~ ¼ κ~ j κ~ l C~ jl . Let us write a similar relation for the quasistatic asymptotics of (43) where Γð

Γu0α ¼ 〈ρ〉c2α u0α

with cα ðκÞ≡lim ωα ðkÞ=k;

ð47Þ

k-0

where k ¼ kκðjκj ¼ 1Þ and ΓðκÞ is to be determined. Comparing (46) and (47) with regard to (45) and making use of the R ~ equality 〈ρ〉≡ ½0;13 ρðAxÞ dx ¼ 〈ρ〉, we find that eff 1 2 1 ~2 ~ 〈ρ〉 k~ j k~ l ~ eff k Γ¼ k Γ ⇒ Γ¼ C ¼ κ p apj C~ jl aql κq ~ ~ k2 jl 〈ρ〉 〈ρ〉 〈ρ〉

ð48Þ

and hence eff eff ~ eff ~ eff ΓðκÞ ¼ Ceff pq κ p κ q : Cpq ¼ apj C jl aql ð⇔C jl ¼ bjp Cpq blq Þ:

ð49Þ

3.3.2. Alternative formulation using anisotropic mass density Premultiplication of Eq. (41) by B allows it to be reformulated as ∂j ðC jl ∂l vÞ ¼ ρω2 v;

ð50Þ

where v and the material parameters are ~ vðxÞ ¼ Aþ vðAxÞð ¼ Aþ vðxÞÞ; ~ ¼ ρ þ ðxÞð ¼ ρðx þ ep ÞÞ; ρðxÞ ¼ BBþ ρðAxÞ ¼ BBþ ρðxÞ c ijkl ðxÞ ¼ bim bjp bkn blq cmpnq ðAxÞð ¼ c ijkl ðx þ ep ÞÞ; þ C jl ðxÞ ¼ ðc ijkl Þ3i;k ¼ 1 ¼ bjp blq BCpq ðAxÞBþ ¼ BC~ pq ðxÞBþ ¼ C lj ðxÞ;

ð51Þ

which are periodic in x with respect to the cubic lattice of vectors ep . Note that the tensor c ijkl retains the major and minor symmetries of normal elasticity, c ijkl ¼ c nklij and c ijkl ¼ c jikl , respectively, while the mass density is no longer a scalar but becomes a symmetric tensor. The Floquet condition now becomes vðxÞ ¼ eikx uðxÞ with periodic uðxÞ satisfying ð∂l þ ik l ÞC lq ð∂q þ ik q Þu ¼ ρ ω 2 u;

ð52Þ

where k ¼ Aþ kð ¼ kκ; jκj ¼ 1Þ;

uðxÞ ¼ Aþ uðAxÞð ¼ uðx þ ep ÞÞ:

ωðkÞ ¼ ωðkÞ;

ð53Þ

Its quasistatic asymptotics are Γu 0α ¼ 〈ρ〉c 2α u 0α eff ΓðκÞ ¼ κ j κ l C jl .

where find that

with c α ðκÞ≡lim ω α ðkÞ=k; α ¼ 1; 2; 3; k-0

ð54Þ

~ ¼ BBþ 〈ρ〉, we Comparing (46) and (47) with regard to (53) and making use of the equality 〈ρ〉 ¼ BBþ 〈ρ〉 2

eff

ð55Þ

þ eff þ Ceff pq ¼ apj AC jl A aql ð⇔C jl ¼ bjp BCpq B blq Þ:

ð56Þ

〈ρ〉1 k Γ ¼ 〈ρ〉1 k AΓAþ ⇒ Γ ¼ κ p apj AC jl Aþ aql κ q 2

and hence eff

eff

3.3.3. Summary of the oblique case Given material constants ρðxÞ and cijkl ðxÞ on an oblique lattice we first identify the matrix A of (39). We may proceed in either of the two ways based on Cosserat elasticity with isotropic density or normal elasticity with anisotropic density. In each case we define material properties on a cubic lattice: ρðxÞ, ~ C~ jl ðxÞ from Eq. (42) or ρðxÞ, C jl ðxÞ from Eq. (51). Then use eff eff the formulas of Section 3 to obtain C~ ðxÞ-C~ or C ðxÞ-C , and finally insert the result into (49) or (56) to arrive at the jl

jl

jl

jl

A.A. Kutsenko et al. / J. Mech. Phys. Solids 61 (2013) 2260–2272

2267

sought effective Christoffel matrix Γ as a function of unit direction vector κ in the oblique lattice. Knowing ΓðκÞ yields the effective speeds cα ðκÞ according to (47). Note that although the formulation in Section 2 was restricted to isotropic density, the quasistatic effective elasticity is the same if one replaces the isotropic density tensor ρðxÞI by the anisotropic density ρðxÞJ with J ¼ Jþ constant positive definite. In the case of the anisotropic density formulation for the oblique lattice J ¼ BBþ . 3.4. Calculating the effective elastic moduli We return to the question of the full determination of ceff ijkl from the Christoffel tensor, or more specifically, from D defined by   eff dikjl ¼ 12 ceff ijkl þ cilkj :

ð57Þ

The elements of D satisfy the same symmetries as those of C ðdikjl ¼ djlik ¼ diklj Þ and they follow from Eq. (13) as       eff ¼ 12 Cejl þ Celj ≡Djl ¼ Dlj : ð58Þ ðdikjl Þ3i;k ¼ 1 ¼ 12 Ceff jl þ Clj   eff ;s eff eff ¼ 1 ceff Define the ‘totally symmetric’ part of C as ceff ijkl þ cikjl þ ciljk . This is seen to be equal to the totally symmetric part ijkl  3 s of D defined in (57), i.e. Ceff;s ¼ Ds where dijkl ¼ 13 dijkl þ dikjl þ diljk . Eq. (57) can then be rewritten as (Norris, 2006) D ¼ 32 Ceff;s 12Ceff ⇒ Ceff ¼ 3Ds 2D: The latter inverse 0 eff c11 ceff 12 B B ceff 22 B B B B B B B B S Y @

ð59Þ

relation may be simply represented in Voigt notation as 1 0 ceff ceff ceff ceff 2d66 d12 2d55 d13 2d56 d14 d 13 14 15 16 C B 11 C B ceff ceff ceff ceff d22 2d44 d23 d24 23 24 25 26 C C B B ceff ceff ceff ceff d33 d34 33 34 35 36 C B C¼B C B d23 ceff ceff ceff 44 45 46 C C B B eff C SYM M ceff c @ 55 56 A ceff 66

d15

1

d16

2d46 d25 d35 d36 d13

C C C 2d45 d36 C C C: d25 C C C d14 A d12 d26

ð60Þ

This one-to-one correspondence between the elements Ceff and D, combined with the identity (58), provides the means to find the effective moduli from Cejl . It remains to determine the full set of elements dijkl from Christoffel tensors for a given set of directions. It is known that data for at least six distinct directions are required (Van Buskirk et al., 1986; Norris, 1989). The necessary and sufficient condition that a given sextet fΓðαÞ ≡Γðκα Þ; α ¼ 1; ‥; 6g will yield the full elastic moduli tensor is that the six directions fκα g do not lie on a cone through the origin and cannot be contained in less than three distinct planes through the origin (Norris, 1989). The set fκα g ¼ fe1 ; e2 ; e3 ; p1ffiffi2 A1 e2 ; p1ffiffi2 A2 e3 ; p1ffiffi2A3 e1 g (see (36)) meets this requirement (and is in fact the set first proposed in Van Buskirk et al., 1986). Thus the complete set of dijkl follows from Norris (1989, Eq. (3.17)) (see Eq. (58))     3 3 1  α α αþ3 þ ∑ Γð9αβÞ κ αj κβl þ κ βj καl : ð61Þ Djl ¼ ∑ ΓðαÞ καj κ αl  pffiffiffi κ αþ3 κ þ κ κ j l l j 2 α¼1 α;β ¼ 1 β4α

The equivalent form of (61) in Voigt notation is 0 ð1Þ 0 1 Γ 11 d11 d12 d13 d14 d15 d16 B B C B Γ ð1Þ d22 d23 d24 d25 d26 C B 22 B B C B ð1Þ B B d33 d34 d35 d36 C B C B Γ 33 B C ¼ B ð1Þ d44 d45 d46 C B Γ 23 B B C B B ð1Þ B S Y M d55 d56 C A B Γ 31 @ @ d66 ð1Þ Γ 12

ð2Þ Γ 11

ð3Þ Γ 11

ð4Þ Γ 11

ð5Þ Γ 11

Γ ð6Þ 11

ð2Þ Γ 22 ð2Þ Γ 33 ð2Þ Γ 23 ð2Þ Γ 31 ð2Þ Γ 12

ð3Þ Γ 22 ð3Þ Γ 33 ð3Þ Γ 23 ð3Þ Γ 31 ð3Þ Γ 12

ð4Þ Γ 22 ð4Þ Γ 33 ð4Þ Γ 23 ð4Þ Γ 31 ð4Þ Γ 12

ð5Þ Γ 22 ð5Þ Γ 33 ð5Þ Γ 23 ð5Þ Γ 31 ð5Þ Γ 12

Γ ð6Þ 22 Γ ð6Þ 33 Γ ð6Þ 23 Γ ð6Þ 31 Γ ð6Þ 12

1

0 C 1 CB CB 0 CB CB CB 0 CB CB 0 CB CB C@ 0 A 0

0

0

0

 12

1

0

 12

0

0

1

 12

 12

0

0

1

0

0

0

0

1

0

0

0

0

 12

1

C  12 C C C 0 C C: 0 C C C 0 A 1

ð62Þ

The full set of ceff ijkl can therefore be obtained from Eqs. (60) and (62) with (see Eq. (37)) Γð1Þ ¼ Ceff 11 ; Γð4Þ ¼

Γð2Þ ¼ Ceff 22 ;

eff 2ðAC 22 Aþ ÞA1 ;

Γð3Þ ¼ Ceff 33 ; eff

Γð5Þ ¼ 2ðAC 33 Aþ ÞA2 ;

eff

Γð6Þ ¼ 2ðAC 11 Aþ ÞA3 :

ð63Þ

Other procedures for inverting a set of compatible Christoffel tensors to give the moduli can be found in Norris (1989). 4. Closed-form upper bound of the effective Christoffel tensor Let N≥0 be the truncation parameter of the 3D or 2D Fourier expansion of Cjl ðxÞ, meaning that the index g in (14) or (26) takes the values from the set ½N; N3 or ½N; N2 ;, respectively. Denote truncated approximations of the effective Christoffel

2268

A.A. Kutsenko et al. / J. Mech. Phys. Solids 61 (2013) 2260–2272

Fig. 1. (a) A cubic lattice of symmetric steel cubes in Epoxy at filling fraction f ¼ 1=8. (b) Effective wave speeds as a function of f. The PWE and MM calculated values are plotted by thin and thick lines (light blue and red online), respectively. The broad curves indicate the Hashin–Shtrikman lower bounds. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this article.)

tensor ΓðκÞ ¼ Ceff jl κ j κ l calculated from (11)1 via the PWE and MM formulas by Γ½NPWE and Γ½NMM , respectively. By analogy with Kutsenko et al. (2012), it can be proved that N1 ≤N2 ⇒ Γ½N1 MM ≥Γ½N 2 MM ;

Γ½N 1 PWE ≥Γ½N2 PWE ;

∀N ⇒ Γ ≤Γ½NMM ≤Γ½NPWE ;

ð64Þ

where the inequality sign between two matrices is understood in the sense that their difference is a sign definite matrix and hence the differences of their similarly ordered eigenvalues are sign definite. The inequalities (64)1 imply that the MM and PWE approximations of Γ obtained by truncating Eq. (11)1 are upper bounds which converge from above to the exact value with growing N. Furthermore, the MM approximation of Γ is more accurate than the PWE one at a given N, from (64)2. Taking N ¼0 in (64)2 yields Γ ≤Γ½0MM ≤Γ½0PWE ¼ 〈Γ〉;

ð65Þ

where Γ½0MM admits an explicit expression which is however rather cumbersome. Seeking a simpler result, consider the above inequalities for one of the principal directions κ∥el so that ΓðκÞ ¼ Ceff ll . Denote the PWE and MM approximations (15) eff eff and (28) of Ceff ll by Cll ½NPWE and Cll ½NMM : From (65) eff eff Ceff ll ≤Cll ½0MM ≤Cll ½0PWE ¼ 〈Cll 〉;

where

Ceff ll ½0MM

Q 0 ½0 ¼

ð66Þ

admits closed-form expression as follows. From (27d) and (27e) taken with N ¼0 (i.e. with g; g′ ¼ 0Þ ! ! I3 〈〈Cll 〉1 〉l 0 〈Cll 〉1 l l ⇒ M0 ð1Þ ¼ ; ð67Þ 0 I3 0 0

so that (28) with N ¼0 yields Ceff ll ½0MM ¼ ð0 I3 ÞS½0≡S2

with S½0 ¼

!1 

0

〈〈Cll 〉1 〉l l

I3

0

0

0

 ≡

S1 S2

! :

ð68Þ

Solving for S2 gives 0

〈〈Cll 〉1 〉l

0

0

l

!

S1

!

S2

 ¼

I3 0



〉1 ⇒ S2 ¼ 〈〈Cll 〉1 l : l

ð69Þ

Thus from (66), (68) and (69) 1 1 Ceff ll ≤〈〈Cll 〉l 〉l ≤〈Cll 〉

ð70Þ

where 〈Cll 〉 is identifiable as the Voigt average (Milton, 2001), known to provide an upper bound. The upper bound provided by the first inequality in (70) has not to our knowledge been presented before. Combining the new bound from (70) with the

A.A. Kutsenko et al. / J. Mech. Phys. Solids 61 (2013) 2260–2272

2269

Voigt inequality Ceff al ≤〈Cal 〉 ða≠lÞ yields 2 〉1 ΓðκÞ ≤ΓB ðκÞ≡ð〈Cal 〉κa κ l Þa≠l þ 〈〈Cll 〉1 l κl ; l

ð71Þ

where the subscript “B” implies bound. Denote the eigenvalues of the matrix ΓB by 〈ρ〉c2Bα and order them in the same way as the eigenvalues 〈ρ〉c2α of Γ, then it follows that cα ðκÞ ≤cBα ðκÞ;

α ¼ 1; 2; 3:

ð72Þ

It will be demonstrated in Section 5 that the upper bounds cBα of the effective speeds can also serve as their reasonable estimate. 5. Numerical examples We consider two examples of 3D phononic crystals composed of steel inclusions in epoxy matrix. The material parameters are c11 ðStÞ ¼ 170 GPa, c66 ðStÞ ¼ 80 GPa, ρðStÞ ¼ 7:7 g=cm3 for steel and c11 ðEpÞ ¼ 7:537 GPa, c66 ðEpÞ ¼ 1:482 GPa, ρðEpÞ ¼ 1:142 g=cm3 for epoxy. This implies cl ðStÞ ¼ 4:7 mm=μs, ct ðStÞ ¼ 3:22 mm=μs and cl ðEpÞ ¼ 2:57 mm=μs, ct ðEpÞ ¼ 1:14 mm=μs for the longitudinal and transverse speeds. The number of Fourier modes is ð2N þ 1Þ3 for the PWE method and ð2N þ 1Þ2 for the MM method; we performed the calculations for N ¼0,3,5. The first example assumes a cubic lattice of cubic steel inclusions (Fig. 1a). We present the effective longitudinal and transverse speeds respectively cl and ct in the principal direction as functions of the volume fraction of steel inclusions (Fig. 1b). The curves calculated by the PWE method are plotted by thin lines (light blue and red online), the curves calculated by the MM method are plotted by thick lines (dark blue and red online). For each method, we present three different data obtained with N ¼ 0; N ¼ 3 and N ¼5 (dotted, dashed and solid lines, respectively). The Hashin–Shtrikman lower bounds (Hashin and Shtrikman, 1963) are also plotted (the Hashin–Shtrikman upper bounds lie far above the other curves and are not displayed). It is observed from Fig. 1b that the results of both methods monotonically converge from above to the exact value with growing N in agreement with the general statement of Section 4. What is significant is that the convergence of the MM method is seen to be much faster than that of the PWE method. In fact, the explicit MM estimate for N ¼0 which follows from (70) in the form: c2l ¼

1 〈〈c11 〉1 〉1 ; 1 1 〈ρ〉

c2t ¼

1 〈〈c66 〉1 〉1 ; 1 1 〈ρ〉

ð73Þ

provides a much better estimate for cl and ct at f 4 0:5 than the PWE calculation with N ¼ 5, i.e. with matrices of about 4000  4000 size. Note that as f -1 in the example of Fig. 1, the bound (70) may be approximated, yielding ! 1=3 1 1 1f eff c11 ⪅ þ : ð74Þ c11 ðStÞ c11 ðEpÞ At the same time the geometry of the unit cell for f -1 indicates that the modulus ceff 11 can be estimated by an equivalent medium stratified in the 1-direction, for which the uniaxial strain assumption with constant stress s11 leads to the 1=3

1 1 approximation ceff , the same as the right hand side of (74) (as f -1Þ. This, combined with the fact that 1f tends 11 ≈〈c11 〉 to zero faster than 1f as f -1, explains the exceptional accuracy of the new bound as an estimate for the moduli. Note that, by comparison, the Hashin–Shtrikman bounds (upper or lower) do not provide a useful estimate in this case. The second example considers a cubic lattice of spheroidal steel inclusions in epoxy matrix. The shape of the inclusion evolves from formally a disk of unit diameter to a ball of unit diameter (that is, inscribed in a cubic unit cell) by means of elongating the radius along the x1 direction, see Fig. 2. We describe the dependence of the effective speeds cl and ct along x1

Fig. 2. A cubic lattice of Steel spheroids in Epoxy matrix. (a) The inclusions are oblate spheroids with minor axis a and unit major axes. (b) The periodic structure for a ¼0.5 and a ¼1 (spheres).

2270

A.A. Kutsenko et al. / J. Mech. Phys. Solids 61 (2013) 2260–2272

Fig. 3. (a) Effective wave speeds for the periodic structure of Fig. 2 as a function of the spheroid minor axis a calculated by the PWE and MM methods. (b) The corresponding elastic moduli.

and of the corresponding effective elastic moduli c11 and c66 on the shape of the spheroidal inclusion. Fourier coefficients for the PWE and MM methods are given in Appendix B. The results are obtained by the PWE method with N ¼3 and N ¼5 (open circles in Fig. 3) and by the MM method with N ¼ 0; N ¼ 3 and N ¼5 (dotted, dashed and solid lines in Fig. 3). The MM method is particularly efficient for the case in hand since it uses Fourier coefficients in the x2 x3 plane where the inclusions have circular cross-section and performs direct numerical integration of the Riccati equation (see Section 3.2) along the direction x1 where the shape is ‘distorted’. We observe an interesting feature of a drastic increase of the effective longitudinal speed cl and of the modulus c11 when the inclusions tend to touch each other, see Fig. 3a. This type of configuration where the inclusions are almost touching is known to be particularly difficult for numerical calculation of the effective properties (Nunan and Keller, 1984). MM appears to be particularly well suited to treating such problems with closely spaced inclusions since it explicitly accounts for the thin gap region via integration of the Riccati equation. The PWE method, on the other hand, clearly fails to capture the sharp increase in wave speed at N ¼5 (matrix size ≈4000  4000). In fact, the PWE for N ¼3 does not even satisfy the strict upper bounds (73) derived from MM at N ¼0.

6. Conclusion The PWE and MM methods of calculating quasistatic effective speeds in three-dimensional phononic crystals have been formulated and compared. The MM method can be viewed as a two-dimensional PWE combined with a one-dimensional propagator matrix approach. The propagator part of the MM scheme is calculated by numerical integration of a (nonlinear) Riccati differential equation to produce the monodromy matrix. It was shown both analytically and numerically that the MM method provides more accurate approximations than the PWE scheme. In particular, the closed form MM bounds (70) (see also (73)) using only one Fourier mode to estimate the effective speed gives better approximations than PWE bounds using more than a thousand (11 in each of xi, i ¼ 1; 2; 3) Fourier modes in the case of densely packed structures (see Fig. 1b for f 40:5). The speed-up of the MM method as compared with PWE via reduction in matrix size is particularly significant for the three-dimensional homogenization problem. Thus, numerical implementation of the PWE scheme needs a matrix of dimension 3ð2N þ 1Þ3  3ð2N þ 1Þ3 , requiring a considerable amount of computer memory even for small N. By contrast, the MM scheme uses matrices of dimension 6ð2N þ 1Þ2  6ð2N þ 1Þ2 . The reduced memory requirement for the MM method is at the cost of the computer time needed to solve the Riccati equation, a relatively small price to pay. In fact, the ability to set the step size in the Runge–Kutta scheme enables the MM method to efficiently and accurately solve configurations for which the PWE is particularly ill-suited, such as narrow gaps (see Fig. 1b for f -1) and closely spaced inclusions (Fig. 3 for a-1).

Acknowledgment A.A.K. acknowledges support from Mairie de Bordeaux. A.N.N. acknowledges support from Institut de Mécanique et d'Ingénierie, Université de Bordeaux.

A.A. Kutsenko et al. / J. Mech. Phys. Solids 61 (2013) 2260–2272

2271

eff

Appendix A. Alternative derivation of Eq. (25) for Cll

Let k be parallel to one of the translation vectors. Take the latter to be e1 ¼ ðδ1i Þ and so k ¼ ðk1 ; 0; 0Þ. Eq. (3) may be rewritten in the form ! ! v 0 0 2 η′ ¼ ðQ0 þ ω Q1 Þη where ′≡∂1 ; η ¼ ; ð75Þ ; Q1 ¼ C1p ∂p v ρI3 0 while Q0 is defined in (18) and (19) with j ¼ 1 and a; b ¼ 2; 3. Denote ηðx1 Þ≡ηðx1 ; x2 ; x3 Þ. The solution to (75) with some initial function ηð0Þ can be written via the matricant in the form ηðx1 Þ ¼ Mðx1 Þηð0Þ

with Mðx1 Þ ¼

Zd x1 0

ðI þ ðQ0 þ ω2 Q1 Þ dx1 Þ:

ð76Þ

Taking into account the assumed 1-periodicity in x1 and hence the Floquet condition v ¼ eik1 x1 u for the solution of (3) implies that the solution η of (75) must satisfy ηð1Þ ¼ eik1 ηð0Þ. Thus, with reference to (76), ωðk1 ; 0; 0Þ≡ωðk1 Þ is an eigenvalue of (3) iff there exists w≡wðx2 ; x3 Þ such that Mð1Þw ¼ eik1 w:

ð77Þ

Consider asymptotic expansion of (77) in small ω; k1 . By (76) Mð1Þ ¼ M0 þ ω2 M1 þ Oðω4 Þ where Z cb M0 ≡M0 ½1; 0; M0 ½b; a ¼ ðI þ Q0 dx1 Þ; a

Z M1 ¼

1

0

M0 ½1; x1 Q1 M0 ½x1 ; 0 dx1 :

ð78Þ

The identity M0 W0 ¼ W0 with the 6  3 matrix W0 ¼ ðI3 0Þþ (see (21)) implies triple multiplicity of the zero-order ω ¼ 0. Therefore we may write 2

ωα ðk1 Þ ¼ cα k1 þ Oðk1 Þ;

2

3

wα ¼ w0α þ k1 w1α þ k1 w2α þ Oðk1 Þ

with w0α ¼ W0 u0α ;

ð79Þ

where α ¼ 1; 2; 3 and  u 0α are some constant linear independent 3  1 vectors. Inserting (78)–(79) along with 2 3 eik1 ¼ 1 þ ik1 12k1 þ O k1 in (77) and equating the terms of the same order in k1 yield 1 : M0 w0α ¼ w0α ;

ð80aÞ

k1 : M0 w1α ¼ iw0α þ w1α ;

ð80bÞ

2

k1 : M0 w2α þ c2α M1 w0α ¼ 12w0α þ iw1α þ w2α :

ð80cÞ

~ 0 ¼ ð0 I3 Þ satisfying Express w1α from (80b) and substitute it in (80c), then scalar multiply the latter by the 6  3 matrix W ~ ~ the identity Mþ 0 ½b; aW 0 ¼ W 0 (see (21)). As a result, we obtain þ

þ

1 ~ Ceff 11 ¼ 〈W 0 ðM0 I Þ W 0 〉1

for κ ¼ e1 ¼ ðδ1i Þ:

ð81Þ

It is seen that (25) with M0 ð1Þ≡M0 and l¼1 is the same as (81). Appendix B. Fourier coefficients for spheroidal inclusions The coefficients for the spheroids of Fig. 2a are as follows: 1. MM method: Identity (27) yields 8 > > > Cpq ðEpÞ;   < b C pq g 2 ; g 3 ; x1 ¼ > > > : Cpq ðEpÞ þ ðCpq ðStÞCpq ðEpÞÞχ^ 1 ðg 2 ; g 3 ; x1 Þ;



1a 1 þ a ; 2 2

1a 1 þ a ; x1 ∈ 2 2 x1 ∉

with qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R J 1 ð2πR g 22 þ g 23 Þ   qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi χ^ 1 g 2 ; g 3 ; x1 ¼ ð1Þg2 þg3 ; g 22 þ g 23 where J1 is the first order Bessel function. 2. PWE method: Identity (13) yields b pq ðgÞ ¼ Cpq ðEpÞδg0 þ ðCpq ðStÞCpq ðEpÞÞχ^ ðgÞ; C 2

R2 ¼ 1

ð2x1 1Þ2 ; a2

2272

A.A. Kutsenko et al. / J. Mech. Phys. Solids 61 (2013) 2260–2272

where χ^ 2 ðgÞ ¼

    að1Þg1 þg2 þg3 ð sin πjga j πjga j cos πjga jÞ ; 2π 2 jga j3

jga j ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðag1 Þ2 þ g 22 þ g 23

and δ is a Kronecker symbol. References Andrianov, I.V., Awrejcewicz, J., Danishevs'kyy, V.V., Weichert, D., 2011. Higher order asymptotic homogenization and wave propagation in periodic composite materials. J. Comput. Nonlinear Dyn. 6, 011015. Hashin, Z., Shtrikman, S., 1963. A variational approach to the elastic behavior of multiphase minerals. J. Mech. Phys. Solids 11, 127–140. Krokhin, A.A., Arriaga, J., Gumen, L.N., 2003. Speed of sound in periodic elastic composites. Phys. Rev. Lett. 91 (26), 264302+. Kushch, V.I., 1987. Computation of the effective elastic moduli of a granular composite material of regular structure. Sov. Appl. Mech. 23, 362–365. Kushch, V.I., 1997. Microstresses and effective elastic moduli of a solid reinforced by periodically distributed spheroidal particles. Int. J. Solids Struct. 34, 1353–1366. Kutsenko, A.A., Shuvalov, A.L., Norris, A.N., Poncelet, O., 2011a. On the effective shear speed in 2D phononic crystals. Phys. Rev. B 84, 064305. Kutsenko, A.A., Shuvalov, A.L., Norris, A.N., 2011b. Evaluation of the effective speed of sound in phononic crystals by the monodromy matrix method. J. Acoust. Soc. Am. 130, 3553–3557. Kutsenko, A.A., Shuvalov, A.L., Norris, A.N., 2012. Converging bounds for the effective shear speed in 2D phononic crystals. J. Elasticity, doi: http://dx.doi.org/ 10.1007/s10659-012-9417-y2. Mei, J., Liu, Z., Wen, W., Sheng, P., 2007. Effective dynamic mass density of composites. Phys. Rev. B 76 (13), 134205+. Milton, G.W., 2001. The Theory of Composites, 1st ed. Cambridge University Press. Nemat-Nasser, S., Taya, M., 1981. On effective moduli of an elastic body containing periodically distributed voids. Q. Appl. Math. 39, 43–59. Nemat-Nasser, S., Iwakuma, T., Hejazi, M., 1982. On composites with periodic structure. Mech. Mater. 1, 239–267. Ni, Q., Cheng, J., 2005. Anisotropy of effective velocity for elastic wave propagation in two-dimensional phononic crystals at low frequencies. Phys. Rev. B 72, 014305. Ni, Q., Cheng, J., 2007. Long wavelength propagation of elastic waves in three-dimensional periodic solid–solid media. J. Appl. Phys. 101, 073515. Norris, A.N., 1989. On the acoustic determination of the elastic moduli of anisotropic solids and acoustic conditions for the existence of planes of symmetry. Q. J. Mech. Appl. Math. 42, 413–426. Norris, A.N., 2006. Elastic moduli approximation of higher symmetry for the acoustical properties of an anisotropic material. J. Acoust. Soc. Am. 119, 2114–2121. Norris, A.N., Shuvalov, A.L., Kutsenko, A.A., 2012. Analytical formulation of 3D dynamic homogenization for periodic elastic systems, Proc. R. Soc. A 468, 1629–1651. Nunan, K.C., Keller, J.B., 1984. Effective elasticity tensor of a periodic composite. J. Mech. Phys. Solids 32 (4), 259–280. Parnell, W.J., D Abrahams, I., 2008. Homogenization for wave propagation in periodic fibre-reinforced media with complex microstructure. I—theory. J. Mech. Phys. Solids 56, 2521–2540. Sangani, A.S., Lu, W., 1987. Elastic coefficients of composites containing spherical inclusions in a periodic array. J. Mech. Phys. Solids 35, 1–21. Torrent, D., Sánchez-Dehesa, J., 2008. Anisotropic mass density by two-dimensional acoustic metamaterials. New J. Phys. 10 (2), 023004+. Van Buskirk, W.C., Cowin, S.C., Carter Jr, R., 1986. A theory of acoustic measurement of the elastic constants of a general anisotropic solid. J. Mater. Sci. 21, 2759–2762. Wu, Y., Zhan, Z.-Q., 2009. Dispersion relations and their symmetry properties of electromagnetic and elastic metamaterials in two dimensions. Phys. Rev. B 79, 195111+.