On the quasistatic effective elastic moduli for elastic waves in three-dimensional phononic crystals
OOn the quasistatic effective elastic moduli for elastic waves inthree-dimensional phononic crystals
A.A. Kutsenko a , A.L. Shuvalov a , A.N. Norris b, ∗ a Institut de M´ecanique et d’Ing´enierie de Bordeaux,Universit´e de Bordeaux, UMR CNRS 5469, Talence 33405, France b Mechanical and Aerospace Engineering, Rutgers University,Piscataway, NJ 08854-8058, USA
Abstract
Effective elastic moduli for 3D solid-solid phononic crystals of arbitrary anisotropy and oblique latticestructure are formulated analytically using the plane-wave expansion (PWE) method and the recently pro-posed monodromy-matrix (MM) method. The latter approach employs Fourier series in two dimensionswith direct numerical integration along the third direction. As a result, the MM method converges muchquicker 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 boundon the effective Christoffel tensor. The MM approach significantly improves the efficiency and accuracyof evaluating effective wave speeds for high-contrast composites and for configurations of closely spacedinclusions, as demonstrated by three-dimensional examples.
1. Introduction
Long-wave low-frequency dispersion of acoustic waves in periodic structures is of both fundamental andpractical interest, particularly due to the current advances in manufacturing of metamaterials and phononiccrystals. In this light, the leading order dispersion (quasistatic limit) has recently been under intensivestudy by various theoretical approaches such as plane-wave expansion (PWE) [1, 2, 3], scaling technique[4, 5], asymptotics of multiple-scattering theory [6, 7, 8] and a newly proposed monodromy-matrix (MM)method [9]. 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 beenproposed for calculating the quasistatic effective elastic properties of 3D periodic composites containingspherical inclusions arranged in a simple cubic array. For the case of rigid inclusions an integral equationon the sphere surface was solved numerically to obtain the effective properties [10]. Spherical voids [11] andsubsequently elastic inclusions were considered using a Fourier series approach [12]. Alternative proceduresfor elastic spherical inclusions include the method of singular distributions [13] and infinite series of periodicmultipole solutions of the equilibrium equations [14]. The latter multipole expansion method has also beenapplied to cubic arrays of ellipsoidal inclusions [15]. A particular PWE-based method of calculation ofquasistatic speeds in 3D phononic crystals of cubic symmetry has been formulated and implemented in[16]. A review of numerical methods for calculating effective properties of composites can be found in [17, § § ∗ Corresponding author
Email addresses: [email protected] (A.A. Kutsenko), [email protected] (A.L. Shuvalov), [email protected] (A.N. Norris)
Preprint submitted to Elsevier November 13, 2018 a r X i v : . [ m a t h - ph ] J un ractical tool for the 3D case, where the vectors and matrices in the Fourier space are of very large algebraicdimension, especially if the phononic crystal is composed of highly contrasting materials (examples in § § § § §
2. Background. PWE formula
Consider a 3D anisotropic medium with density and elastic stiffness ρ ( x ) = ρ ( x + e p ) , C ( x ) = C ( x + e p ) , (1)which are assumed to be -periodic, i.e. invariant to period or translation vectors e p = ( δ pq ) (a cubiclattice, otherwise see § ∗ and + mean complexand Hermitian conjugation. Assume no dissipation so that the elements of C satisfy c ijkl = c ∗ klij (= c klij for real case). Our goal is the quasistatic effective elastic stiffness C eff with elements c eff ijkl that have thesame symmetries as those of C , and matrices C eff jl defined by analogy with Eq. (2). For compact writing,introduce the matrices C jl = ( c ijkl ) i,k =1 = C + lj (2)with components numbered by i, k . The elastodynamic equation for time-harmonic waves v ( x , t ) = v ( x ) e − iωt is ∂ j ( C jl ∂ l v ) = − ρω v , (3)where ∂ j ≡ ∂/∂x j and repeated indices are summed. The differential operator in Eq. (3) is self-adjoint withrespect to the Floquet condition v ( x ) = u ( x ) e i k · x with -periodic u ( x ) = u ( x + e p ) and k = kκκκ ( | κκκ | = 1).Substituting this condition in (3) casts it into the form( C + k C + k C ) u = ρω u , where C u ≡ − ∂ j ( C jl ∂ l u ) , C u ≡ − i ( κ j C jl ∂ l + κ l ∂ j C jl ) u , C u ≡ κ j κ l C jl u . (4)All operators C are self-adjoint. We introduce for future use linear, areal and volumetric averages over theunit-cell: (cid:104)·(cid:105) j is the average over coordinate x j ; (cid:104)·(cid:105) j is the average over the section orthogonal to x j , and (cid:104)·(cid:105) is the complete average. These averages in turn define inner products of vector-valued functions. Thus, fora scalar function f and vector-functions f , h , (cid:104) f (cid:105) j = (cid:90) f dx j , (cid:104) f (cid:105) = (cid:90) [0 , f d x , ( f , h ) j = (cid:104) h + f (cid:105) j , ( f , h ) = (cid:104) h + f (cid:105) . (5)2ext we apply perturbation theory to (4) with a view to defining the quasistatic effective Christoffelmatrix whose eigenvalues yield the effective speeds c α ≡ c α ( κκκ ) ≡ lim k → ω α ( k ) /k, α = 1 , , . (6)For k = 0 , the eigenvalue ω = 0 has multiplicity 3 and corresponds to three constant linear independenteigenvectors u α . Consider the asymptotics ω α = 0 + kλ α + k λ α + O ( k ) , (7a) u α = u α + k u α + k u α + O ( k ) . (7b)Substituting (7) into (4) and collecting terms with the same power of k yields1 : C u α = , (8a) k : C u α + C u α = ρλ α u α , (8b) k : C u α + C u α + C u α = ρλ α u α + ρλ α u α . (8c)Scalar multiplying (8b) by u α and (8c) by e k , and using ( C u α , u α ) = 0 together with self-adjointness of C leads to λ α = 0; u α = −C − C u α ; ( C u α , e k ) + ( C u α , e k ) = (cid:104) ρ (cid:105) λ α ( u α , e k ) , k = 1 , , , (9)where λ α = c α due to λ α = 0 and (6). Inserting u α from (9) in (9) gives Γu α = (cid:104) ρ (cid:105) c α u α , α = 1 , , , (10)where Γ = (( C e k , e i ) − ( C − C e k , C e i )) i,k =1 . Substituting C , C from (4) defines the quasistatic effective3 × Γ = (Γ ik ) i,k =1 in the form Γ ( κ ) = C e jl κ j κ l where C e jl ≡ (cid:104) C jl (cid:105) − (cid:104) ( ∂ p C + pj ) C − ( ∂ q C ql ) (cid:105) (cid:0) = C e + lj (cid:1) . (11)The matrix C e jl is distinguished from C eff jl , in terms of which the Christoffel matrix is Γ ( κ ) = C eff jl κ j κ l = 12 (cid:0) C eff jl + C eff lj (cid:1) κ j κ l . (12)Comparison of Eqs. (11) and (12) implies that C e jl and C eff jl are related by C eff jl + C eff lj = C e jl + C e lj (13)and they are equal if j = l . Equation (11) does not yield C eff jl explicitly, only in the combination (13);however, this connection is sufficient for the purpose of finding all elements of C eff , as described in § C e jl .Equation (11) is still an implicit formula for the matrix C e jl in so far as the operator C − is not specified.One way to an explicit implementation of (11) is via its transformation to Fourier space which must betruncated to define C − as a matrix inverse. Thus taking the 3D Fourier expansion C jl ( x ) = (cid:88) g ∈ Z (cid:98) C jl ( g ) e πi g · x , (cid:98) C jl ( g ) = (cid:104) C jl ( x ) e − πi g · x (cid:105) (14)and plugging it into (11) yields the PWE formula for C e jl as follows C e jl = (cid:98) C jl ( ) − q + j C − q l , where q j = ( g i (cid:98) C ij ( g )) g ∈ Z \ , C = ( g k g (cid:48) p (cid:98) C kp ( g − g (cid:48) )) g , g (cid:48) ∈ Z \ . (15)This result corresponds to the quasistatic limit of the effective elastic coefficients for a dynamically homog-enized periodic medium, see [18, eq. (2.11)].Unfortunately, although the PWE formula (15) is straightforward, its numerical use in 3D is complicatedby the large algebraic dimensions of vectors and matrices of Fourier coefficients. This motivates using themonodromy matrix (MM) method which confines the PWE computation to 2D, see next Section. The MMmethod will also be shown to have significantly faster convergence than PWE.3 . MM formula The idea underlying the MM method is to reduce the 3-dimensional problem of Eq. (11) to an equivalent1-dimensional equation that can be integrated. This is achieved by focusing on a single coordinate and usingFourier transforms in the orthogonal coordinates. The MM formula may be deduced proceeding from Eq.(11) ; using integration by parts, let us recast it into the form C e jl = (cid:104) C jl (cid:105) − (cid:104) ( ∂ p C + pj ) U (cid:105) = (cid:104) C jp ∂ p V (cid:105) , (16)where we have denoted C − ( ∂ q C ql ) = U and V = U + x l I with I standing for the 3 × C U = ∂ q C ql ⇔ C U = −C x l I ⇔ C V = 0 ⇔ ∂ q ( C qp ∂ p V ) = 0 with V = U + x l I . (17)In the following derivation the indices j, l are regarded as fixed, all repeated indices are summed and amongthem the indices a, b are specialized by the condition a, b (cid:54) = l. The suffix 0 of the differential operators Q , M below indicates their reference to ω, k = 0 (similarly to C ).Equation (17) can be rewritten as an ordinary differential equation in the designated coordinate x l Ξ (cid:48) = Q Ξ with (cid:48) ≡ ∂ l , Ξ = (cid:18) VC lp ∂ p V (cid:19) , Q = (cid:18) −B − A B − A − A +1 B − A A +1 B − (cid:19) , (18)where Ξ is a 6 × B , A and A are defined by B V = C ll V , A V = C lb ∂ b V , A V = − ∂ a ( C ab ∂ b V ) . (19)Note that B and A are self-adjoint at any fixed x l with respect to the inner product ( · , · ) l (see (5)). Denote Ξ ( s ) ≡ Ξ ( x ) | x l = s . The solution to (18) with some initial matrix function Ξ (0) has the form Ξ ( x l ) = M ( x l ) Ξ (0) with M ( x l ) = (cid:100)(cid:90) x l ( I + Q dx l ) , (20)where M is a propagator matrix defined through the multiplicative integral (cid:98)(cid:82) and I is the identity operator.Note the identitiesfor W = (cid:18) I (cid:19) , (cid:102) W = (cid:18) (cid:19) : Q W = , Q +0 (cid:102) W = ⇒ M W = W , M +0 (cid:102) W = (cid:102) W . (21)It is seen that for any value of x l the operator M − I has no inverse but at the same time it is aone-to-one mapping from the subspace orthogonal to W onto the subspace orthogonal to (cid:102) W . By thedefinitions of V and Ξ in (17), (18) and due to periodicity of U , it follows that Ξ (1) = Ξ (0) + W ⇒ Ξ (0) = ( M (1) − I ) − W ≡ S . (22)Thus from (20), (22) and (18), (21), Ξ ( x l ) = M ( x l ) S ⇒ V = W +0 M ( x l ) S , C lp ∂ p V = (cid:102) W +0 M ( x l ) S . (23)Substituting V obtained from (23) into (16) yields the desired formula for C e jl , which is discussed furtherin § M is self-adjoint with respect to ( · , · ) l , (cid:104) C lp ∂ p V (cid:105) l = ( M ( x l ) S , (cid:102) W ) l = ( S , M +0 ( x l ) (cid:102) W ) l = ( S , (cid:102) W ) l = (cid:104) (cid:102) W +0 S (cid:105) l . (24)The latter identity implies that the laterally averaged ’traction’ component of Ξ ( x l ) is independent of x l ,and may be identified as a net ’force’ acting on the faces x l = constant.4urther simplification can be achieved for the matrices C eff ll (= C e ll ). By (16) and (58) C eff ll = (cid:104) C lp ∂ p V (cid:105) = (cid:104) (cid:102) W +0 S (cid:105) ≡ (cid:104) (cid:102) W +0 ( M (1) − I ) − W (cid:105) l . (25)Equation (25) suffices to define the effective Christoffel tensor for k = kκ parallel to the translation vector e l in which case Γ ( κ ) = C eff ll . Note that an alternative derivation of (25) is given in Appendix 1.Finally it is noted that M (1) which appears in the above expressions is formally a monodromy matrix(MM) relatively to the coordinate x l , for which reason this approach and its outcome formulas are referredto as the MM ones. The MM approach to finding effective speed of shear (scalar) waves in 2D structureswas first presented in [9], where Eq. (25) for vector waves was also mentioned but with neither derivationnor discussion .3.2. Implementation3.2.1. Propagator matrix in direction e l As above, we fix l and keep a, b (cid:54) = l . Introduce the 2D Fourier expansion C pq ( x ) = (cid:88) g ∈ Z e πig a x a (cid:98) C pq ( g , x l ) , (cid:98) C pq ( g , x l ) = (cid:104) C pq ( x ) e − πig a x a (cid:105) l . (26)Operators B , A , A and matrix operators Q , M ( x l ) defined in (18)-(20) are represented in the 2D Fourierspace by the following infinite matrices truncated to a finite size in calculations: B (cid:55)→ B = ( (cid:98) C ll ( g − g (cid:48) , x l )) g , g (cid:48) ∈ Z , (27a) A (cid:55)→ A = 2 πi ( (cid:98) C lb ( g − g (cid:48) , x l ) g (cid:48) b ) g , g (cid:48) ∈ Z , (27b) A (cid:55)→ A = 4 π ( g a (cid:98) C ab ( g − g (cid:48) , x l ) g (cid:48) b ) g , g (cid:48) ∈ Z , (27c) Q (cid:55)→ Q = (cid:18) − B − A B − A − A +1 B − A A +1 B − (cid:19) , (27d) M ( x l ) (cid:55)→ M ( x l ) = (cid:100)(cid:90) x l ( I + Q dx l ) . (27e) κ (cid:107) e l Consider the effective matrices C eff ll (= Γ if κ (cid:107) e l ). In view of the above notations, formula (25) for C eff ll is expressed as C eff ll = (cid:102) W + (cid:98) (cid:98) S with (cid:98) S = ( M (1) − I ) − W (cid:98) , (28)where W (cid:98) = (cid:18) E (cid:98) (cid:19) , (cid:102) W (cid:98) = (cid:18) (cid:98) (cid:19) with E (cid:98) = ( δ g0 I ) g ∈ Z . (29)Calculation of (cid:98) S in (28) can be performed by means of calculating M from its definition (27e), usingany of the known methods for evaluating multiplicative integrals. However, this approach may suffer fromnumerical instabilities for M of large size due to exponential growth of some components of M . Anothermethod rests on calculating the resolvent ( M − α I ) − without calculating M . The advantage of doingso is due to the fact that growing dimension of M leads to a relatively moderate growth of the resolventcomponents. Let us consider this latter method in some detail.Denote the resolvent R α ( x l ) = ( M ( x l ) − α I ) − where α is not an eigenvalue of M ( x l ). Since the matrix M ( x l ) satisfies the differential equation M (cid:48) ( x l ) = Q ( x l ) M ( x l ) with the initial condition M (0) = I , itfollows that R α ( x l ) with a randomly chosen α (cid:54) = 1 satisfies a Riccati equation with initial condition asfollows (cid:40) R (cid:48) α = − R α Q ( I + α R α ) , R α (0) = (1 − α ) − I , (30)5here (cid:48) = ∂ l . Since eigenvalues of M usually tend to lie close to the real axis and unit circle, it is rec-ommended to take α ∈ C far from these sets. Integrating Eq. (30) numerically (we used the Runge-Kuttamethod of fourth order) provides R α (1) = ( M (1) − α I ) − where α (cid:54) = 1 . To exploit it for finding C eff ll givenby (28) we use the identity (cid:98) S ≡ ( M (1) − I ) − W (cid:98) = ( I + ( α − R α (1)) − R α (1) W (cid:98) = ( I + ( α − R α (1)) − W (cid:98) − α . (31)Thus (cid:98) S is found from the linear system(1 − α ) (cid:0) I + ( α − R α (1) (cid:1)(cid:98) S = W (cid:98) . (32)To solve (32), we first note that the matrix T ≡ (1 − α ) (cid:0) I + ( α − R α (1) (cid:1) satisfies TW (cid:98) = , (cid:102) W + (cid:98) T = with W (cid:98) , (cid:102) W (cid:98) from (29) and thus has 3 zero columns on the left of its vertical midline and 3 zero rows belowthe horizontal midline. Removing these columns and rows yields the reduced matrix (cid:101) T = (cid:16) T T T T (cid:17) wherethe square block T has 3 rows and 3 columns less than the block T . Since T is numerically large, while T , T are medium and T small, it is convenient to apply Schur’s formula to arrive at the final relation inthe form C eff ll = E + (cid:98) (cid:0) T − T T − T (cid:1) − E (cid:98) . (33) Consider k = kκ of arbitrary orientation relatively to the translations e l . In order to find the effectiveChristoffel tensor Γ ( κ ) of Eq. (12) for arbitrary κ ( ∦ e l ) , we need to calculate C e jl with j (cid:54) = l. Applying the2D Fourier expansion to Eqs. (16) and (23) yields C e jl = (cid:104) E + (cid:98) ( A + C jl ∂ l ) (cid:98) V (cid:105) l = C eff ll + (cid:104) E + (cid:98) ( C ll − C jl ) B − ( A (cid:98) V − (cid:98) N ) (cid:105) l , where (cid:32) (cid:98) V (cid:98) N (cid:33) = M ( x l ) (cid:98) S (34)and C jp = ( (cid:98) C jp ( g − g (cid:48) , x l )) g , g (cid:48) ∈ Z with g = ( g a g b ) and a, b (cid:54) = l . As detailed above, the matrix (cid:98) S can be calculated with a very good precision. Evaluation of (cid:98) V ( x l ) and (cid:98) N ( x l ) is no longer that accurate,particularly for high-contrast structures which require many terms in the Fourier series and hence need M of large algebraic dimension and therefore with large values of some components. Thus a negligible errorin S may become noticeable after multiplying by M . In this regard, we present an alternative method ofcalculating C e jl with j (cid:54) = l which circumvents (34).For the fixed functions ρ ( x ) and c ijkl ( x ) with a given cubic lattice of periods e p = ( δ pq ), introduce thenew periods a p = Ae p where A has integer components a ij . The solutions v ( x ) = u ( x ) e i k · x and ω ( k ) ofthe wave equation (3) remain unaltered, as they do not depend on the choice of periods. But now in orderto define C eff jl by the MM formula (which requires periods to coincide with base vectors) we should apply thechange of variables x → Ax that leads, as explained in § e p = ( δ pq ). According to Eqs. (49) and (56) of § C eff pq = a pj (cid:101) C eff jl a ql = a pj AC eff jl A + a ql (cid:16) ⇔ (cid:101) C eff jl = b jp C eff pq b lq , C eff jl = b jp BC eff pq B + b lq (cid:17) (35)where (cid:101) C eff jl and C eff jl are the effective matrices associated with the new profiles (cid:101) c ijkl ( x ) = b jp b lq c ipkq ( Ax )and c ijkl ( x ) = b im b jp b kn b lq c mpnq ( Ax ), respectively, and b ij stand for components of A − . We may apply(35) successively for each of the transformations a p = A j e p , j = 1 , , A = −
10 1 1 ⇔ a = e , a = e + e , a = − e + e , A = − ⇔ ..., A = − ⇔ ... . (36) C eff23 + C eff32 = 4( (cid:101) C eff22 ) A − C eff22 − C eff33 = 4( AC eff22 A + ) A − C eff22 − C eff33 , C eff31 + C eff13 = 4( (cid:101) C eff33 ) A − C eff33 − C eff11 = 4( AC eff33 A + ) A − C eff33 − C eff11 , C eff12 + C eff21 = 4( (cid:101) C eff11 ) A − C eff11 − C eff22 = 4( AC eff11 A + ) A − C eff11 + C eff22 . (37)Inserting (37) in (12) eliminates C eff jl + C eff lj with j (cid:54) = l and expresses the effective Christoffel tensor Γ fullyin terms of the matrices C eff ll and either of ( (cid:101) C eff ll ) A n or ( C eff ll ) A n ( n = 1 , , c ijkl ( x ) replaced by (cid:101) c ijkl ( x ) or c ijkl ( x ), respectively. Thus all calculations have been reduced tothe form (28) which ensures a numerically stable evaluation of Γ .Note that the transformed elasticity c ijkl ( x ) retains the usual symmetries under the interchange ofindices, while (cid:101) c ijkl ( x ) does not (see § √ A j , j = 1 , ,
3, are orthogonal matrices of unitdeterminant. Hence, apart from a factor of , (cid:0) c ijkl ( x ) (cid:1) A j is precisely the elasticity tensor represented in acoordinate system rotated about the axis e j by π from the original. Consider the problem of quasistatic asymptotics of the wave equation (3) for the general case of a 3Dperiodic elastic medium with ρ ( x ) = ρ ( x + a p ) , c ijkl ( x ) = c ijkl ( x + a p ) , (38)where the translation vectors a p form an oblique lattice. We will define the solution of this problem via thesolution for a simpler case of a cubic lattice.The oblique lattice vectors are defined by a matrix A ( (cid:54) = I ) as a p = Ae p with A = ( a pq ) p,q =1 = ( a p · e q ) p,q =1 ; B ≡ A − = ( b pq ) p,q =1 (39)where e p are the orthonormal vectors used previously. Define the new or transformed position variable x (cid:48) = Bx ( ⇔ x = Ax (cid:48) ), the associated displacement (cid:101) v ( x (cid:48) ) = v ( x ) and material parameters (cid:101) ρ ( x (cid:48) ) = ρ ( x ), c (1) ijkl ( x (cid:48) ) = c ijkl ( x ), which are seen to be periodic in x (cid:48) with respect to the vectors e p . Setting C (1) jl = (cid:0) c (1) ijkl (cid:1) i,k =1 , the equation of motion (3) becomes b jp b lq ∂ j (cid:48) (cid:0) C (1) pq ∂ l (cid:48) (cid:101) v (cid:1) = − (cid:101) ρω (cid:101) v , (40)where ∂ j (cid:48) ≡ ∂/∂x (cid:48) j . Using the fact that B is constant allows it to be removed explicitly from (40) byincorporation into a newly defined stiffness tensor. Thus, replacing x (cid:48) → x we have ∂ j ( (cid:101) C jl ∂ l (cid:101) v ) = − (cid:101) ρω (cid:101) v , (41)where (cid:101) v ( x ) = v ( Ax ) and the material parameters are (cid:101) ρ ( x ) = ρ ( Ax ) (= (cid:101) ρ ( x + e p )) , (cid:101) c ijkl ( x ) = b jp b lq c ipkq ( Ax ) (= (cid:101) c ijkl ( x + e p )) , (cid:101) C jl ( x ) = ( (cid:101) c ijkl ) i,k =1 = b jp b lq C pq ( Ax ) = (cid:101) C + lj ( x ) , (42)which are periodic in x with respect to the cubic lattice formed by vectors e p . Note that the tensor (cid:101) c ijkl for A (cid:54) = I is of Cosserat type in that it is not invariant to permutations of indices i (cid:29) j and k (cid:29) l but retainsthe major symmetry (cid:101) c ijkl = (cid:101) c ∗ klij (= (cid:101) c klij for real case) . v ( x ) = e i k · x u ( x ) with periodic u ( x ) = u ( x + a j ) satisfying − ( ∂ l + ik l ) C lq ( ∂ q + ik q ) u = ρω u (43)is equivalent to the condition (cid:101) v ( x ) = e i (cid:101) k · x (cid:101) u ( x )with periodic (cid:101) u ( x ) satisfying the equation that follows from(40), − ( ∂ l + i (cid:101) k l ) (cid:101) C lq ( ∂ q + i (cid:101) k q ) (cid:101) u = (cid:101) ρ (cid:101) ω (cid:101) u , (44)where (cid:101) k = A + k ( = (cid:101) k (cid:101) κ, | (cid:101) κ | = 1) , (cid:101) ω ( (cid:101) k ) = ω ( k ) , (cid:101) u ( x ) = u ( Ax ) (= (cid:101) u ( x + e p )) . (45)Equation (44), which is defined on a cubic lattice, has quasistatic asymptotics as described above. Accordingto (10), (cid:101) Γ (cid:101) u α = (cid:104) (cid:101) ρ (cid:105) (cid:101) c α (cid:101) u α with (cid:101) c α ( (cid:101) κ ) ≡ lim (cid:101) k → (cid:101) ω α ( (cid:101) k ) / (cid:101) k, α = 1 , , , (46)where (cid:101) Γ ( (cid:101) κ ) = (cid:101) κ j (cid:101) κ l (cid:101) C eff jl . Let us write a similar relation for the quasistatic asymptotics of (43), Γu α = (cid:104) ρ (cid:105) c α u α with c α ( κ ) ≡ lim k → ω α ( k ) /k, (47)where k = k κ ( | κ | = 1) and Γ ( κ ) is to be determined. Comparing (46) and (47) with regard for (45) andmaking use of the equality (cid:104) (cid:101) ρ (cid:105) ≡ (cid:82) [0 , ρ ( Ax ) d x = (cid:104) ρ (cid:105) , we find that1 (cid:104) ρ (cid:105) k Γ = 1 (cid:104) (cid:101) ρ (cid:105) (cid:101) k (cid:101) Γ ⇒ Γ = (cid:104) ρ (cid:105)(cid:104) (cid:101) ρ (cid:105) (cid:101) k j (cid:101) k l k (cid:101) C eff jl = κ p a pj (cid:101) C eff jl a ql κ q (48)and hence Γ ( κ ) = C eff pq κ p κ q : C eff pq = a pj (cid:101) C eff jl a ql (cid:16) ⇔ (cid:101) C eff jl = b jp C eff pq b lq (cid:17) . (49) Premultiplication of Eq. (41) by B allows it to be reformulated as ∂ j ( C jl ∂ l v ) = − ρ ω v , (50)where v and the material parameters are v ( x ) = A + v ( Ax ) (cid:0) = A + (cid:101) v ( x ) (cid:1) , ρ ( x ) = BB + ρ ( Ax ) = BB + (cid:101) ρ ( x ) = ρ + ( x ) (= ρ ( x + e p )) ,c ijkl ( x ) = b im b jp b kn b lq c mpnq ( Ax ) (= c ijkl ( x + e p )) , C jl ( x ) = ( c ijkl ) i,k =1 = b jp b lq BC pq ( Ax ) B + = B (cid:101) C pq ( x ) B + = C + lj ( x ) , (51)which are periodic in x with respect to the cubic lattice of vectors e p . Note that the tensor c ijkl retains themajor and minor symmetries of normal elasticity, c ijkl = c ∗ klij and c ijkl = c jikl , while the mass density is nolonger a scalar but becomes a symmetric tensor.The Floquet condition now becomes v ( x ) = e i k · x u ( x ) with periodic u ( x ) satisfying − ( ∂ l + ik l ) C lq ( ∂ q + ik q ) u = ρ ω u , (52)where k = A + k ( = kκ, | κ | = 1) , ω ( k ) = ω ( k ) , u ( x ) = A + u ( Ax ) (= u ( x + e p )) . (53)Its quasistatic asymptotics are Γu α = (cid:104) ρ (cid:105) c α u α with c α ( κ ) ≡ lim k → ω α ( k ) /k, α = 1 , , , (54)8here Γ ( κ ) = κ j κ l C eff jl . Comparing (47) and (54) with regard for (53) and making use of the equality (cid:104) ρ (cid:105) = BB + (cid:104) (cid:101) ρ (cid:105) = BB + (cid:104) ρ (cid:105) , we find that (cid:104) ρ (cid:105) − k Γ = (cid:104) ρ (cid:105) − k AΓA + ⇒ Γ = κ p a pj AC eff jl A + a ql κ q (55)and hence C eff pq = a pj AC eff jl A + a ql (cid:16) ⇔ C eff jl = b jp BC eff pq B + b lq (cid:17) . (56) Given material constants ρ ( x ) and c ijkl ( x ) on an oblique lattice we first identify the matrix A of (39).We may proceed in either of two ways based on Cosserat elasticity with isotropic density, or normal elasticitywith anisotropic density. In each case we define material properties on a cubic lattice: (cid:101) ρ ( x ), (cid:101) C jl ( x ) fromEq. (42) or ρ ( x ), C jl ( x ) from Eq. (51), respectively. Then use the formulas of § (cid:101) C jl ( x ) → (cid:101) C eff jl or C jl ( x ) → C eff jl , and finally insert the result into (49) or (56) to arrive at the sought effective Christoffelmatrix Γ as a function of unit direction vector κ in the oblique lattice. Knowing Γ ( κ ) yields the effectivespeeds c α ( κ ) according to (47). Note that although the formulation in § ρ ( x ) I by theanisotropic density ρ ( x ) J with J = J + constant positive definite. In the case of the anisotropic densityformulation for the oblique lattice J = BB + . We return to the question of the full determination of c eff ijkl from the Christoffel tensor, or more specifically,from D defined by d ikjl = 12 (cid:0) c eff ijkl + c eff ilkj (cid:1) . (57)The elements of D satisfy the same symmetries as those of C ( d ikjl = d jlik = d iklj ) and they follow fromEq. (13) as ( d ikjl ) i,k =1 = 12 (cid:0) C eff jl + C eff lj (cid:1) = 12 (cid:0) C e jl + C e lj (cid:1) ≡ D jl (cid:0) = D lj (cid:1) . (58)Define the ’totally symmetric’ part of C eff as c eff , s ijkl = (cid:0) c eff ijkl + c eff ikjl + c eff iljk (cid:1) . This is seen to be equal to thetotally symmetric part of D defined in (57), i.e. C eff , s = D s where d s ijkl = (cid:0) d ijkl + d ikjl + d iljk (cid:1) . Equation(57) can then be rewritten as [19] D = 32 C eff , s − C eff ⇒ C eff = 3 D s − D . (59)The latter inverse relation may be simply represented in Voigt notation as c eff11 c eff12 c eff13 c eff14 c eff15 c eff16 c eff22 c eff23 c eff24 c eff25 c eff26 c eff33 c eff34 c eff35 c eff36 c eff44 c eff45 c eff46 S Y M c eff55 c eff56 c eff66 = d d − d d − d d − d d d d d − d d d − d d d d d d − d d d d SY M d d d . (60)This one-to-one correspondence between the elements C eff and D , combined with the identity (58), providesthe means to find the effective moduli from C e jl .It remains to determine the full set of elements d ijkl from Christoffel tensors for a given set of directions.It is known that data for at least six distinct directions are required [20, 21]. The necessary and sufficient9ondition that a given sextet { Γ ( α ) ≡ Γ ( κ α ) , α = 1 , .., } will yield the full elastic moduli tensor is thatthe six directions { κ α } do not lie on a cone through the origin and cannot be contained in less than threedistinct planes through the origin [21]. The set { κ α } = { e , e , e , √ A e , √ A e , √ A e } (see (36))meets this requirement (and is in fact the set first proposed in [20]). Thus the complete set of d ijkl followsfrom [21, Eq. (3.17)] (see Eq. (58)) D jl = (cid:88) α =1 Γ ( α ) (cid:0) κ αj κ αl − √ ( κ α +3 j κ αl + κ αj κ α +3 l ) (cid:1) + (cid:88) α,β =1 β>α Γ (9 − α − β ) ( κ αj κ βl + κ βj κ αl ) . (61)The equivalent form of (61) in Voigt notation is d d d d d d d d d d d d d d d d d d S Y M d d d = Γ (1)11 Γ (2)11 Γ (3)11 Γ (4)11 Γ (5)11 Γ (6)11 Γ (1)22 Γ (2)22 Γ (3)22 Γ (4)22 Γ (5)22 Γ (6)22 Γ (1)33 Γ (2)33 Γ (3)33 Γ (4)33 Γ (5)33 Γ (6)33 Γ (1)23 Γ (2)23 Γ (3)23 Γ (4)23 Γ (5)23 Γ (6)23 Γ (1)31 Γ (2)31 Γ (3)31 Γ (4)31 Γ (5)31 Γ (6)31 Γ (1)12 Γ (2)12 Γ (3)12 Γ (4)12 Γ (5)12 Γ (6)12 − − − − − −
00 0 0 1 0 00 0 0 0 1 00 0 0 0 0 1 . (62)The full set of c eff ijkl can therefore be obtained from Eqs. (60) and (62) with (see Eq. (37)) Γ (1) = C eff11 , Γ (2) = C eff22 , Γ (3) = C eff33 , Γ (4) = 2( AC eff22 A + ) A , Γ (5) = 2( AC eff33 A + ) A , Γ (6) = 2( AC eff11 A + ) A . (63)Other procedures for inverting a set of compatible Christoffel tensors to give the moduli can be found in[21].
4. Closed-form upper bound of the effective Christoffel tensor
Let N ≥ C jl ( x ), meaning that theindex g in (14) or (26) takes the values from the set [ − N, N ] or [ − N, N ] , respectively. Denote truncatedapproximations of the effective Christoffel tensor Γ ( κ ) = C eff jl κ j κ l calculated from (11) via the PWE andMM formulas by Γ [ N ] PWE and Γ [ N ] MM , respectively. By analogy with [22], it can be proved that N ≤ N ⇒ Γ [ N ] MM ≥ Γ [ N ] MM , Γ [ N ] PWE ≥ Γ [ N ] PWE ; ∀ N ⇒ Γ ≤ Γ [ N ] MM ≤ Γ [ N ] PWE , (64)where the inequality sign between two matrices is understood in the sense that their difference is a sign defi-nite matrix and hence the differences of their similarly ordered eigenvalues are sign definite. The inequalities(64) imply that the MM and PWE approximations of Γ obtained by truncating Eq. (11) are upper boundswhich 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) .Taking N = 0 in (64) yields Γ ≤ Γ [0] MM ≤ Γ [0] PWE = (cid:104) Γ (cid:105) , (65)where Γ [0] MM admits an explicit expression which is however rather cumbersome. Seeking a simpler result,consider the above inequalities for one of the principal directions κ (cid:107) e l so that Γ ( κ ) = C eff ll . Denote thePWE and MM approximations (15) and (28) of C eff ll by C eff ll [ N ] PWE and C eff ll [ N ] MM . From (65), C eff ll ≤ C eff ll [0] MM ≤ C eff ll [0] PWE = (cid:104) C ll (cid:105) , (66)10here C eff ll [0] MM admits closed-form expression as follows. From (27d) and (27e) taken with N = 0 (i.e.with g , g (cid:48) = ), Q [0] = (cid:18) (cid:104) C ll (cid:105) − l (cid:19) ⇒ M (1) = (cid:18) I (cid:104)(cid:104) C ll (cid:105) − l (cid:105) l (cid:19) , (67)so that (28) with N = 0 yields C eff ll [0] MM = (cid:0) (cid:1) S [0] ≡ S with S [0] = (cid:18) (cid:104)(cid:104) C ll (cid:105) − l (cid:105) l (cid:19) − (cid:18) I (cid:19) ≡ (cid:18) S S (cid:19) . (68)Solving for S gives (cid:18) (cid:104)(cid:104) C ll (cid:105) − l (cid:105) l (cid:19) (cid:18) S S (cid:19) = (cid:18) I (cid:19) ⇒ S = (cid:104)(cid:104) C ll (cid:105) − l (cid:105) − l . (69)Thus from (66), (68) and (69), C eff ll ≤ (cid:104)(cid:104) C ll (cid:105) − l (cid:105) − l ≤ (cid:104) C ll (cid:105) (70)where (cid:104) C ll (cid:105) is identifiable as the Voigt average [17], known to provide an upper bound. The upper boundprovided by the first inequality in (70) has not to our knowledge been presented before.Combining the new bound from (70) with the Voigt inequality C eff al ≤ (cid:104) C al (cid:105) ( a (cid:54) = l ) yields Γ ( κ ) ≤ Γ B ( κ ) ≡ ( (cid:104) C al (cid:105) κ a κ l ) a (cid:54) = l + (cid:104)(cid:104) C ll (cid:105) − l (cid:105) − l κ l , (71)where the subscript ”B” implies bound. Denote the eigenvalues of the matrix Γ B by (cid:104) ρ (cid:105) c α and order themin the same way as the eigenvalues (cid:104) ρ (cid:105) c α of Γ , then it follows that c α ( κ ) ≤ c B α ( κ ) , α = 1 , , . (72)It will be demonstrated in § c B α of the effective speeds can also serve as theirreasonable estimate. 11) fc c l (Ep) c t (Ep) c t (St) c l (St) PWE MM N =5 N =5 N =3 N =3 N =0 N =0HS bounds b) Figure 1: (a) A cubic lattice of symmetric steel cubes in Epoxy at filling fraction f = 1 /
8. (b) Effective wave speeds as afunction 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.
5. Numerical examples
We consider two examples of 3D phononic crystals composed of steel inclusions in epoxy matrix. Thematerial parameters are c (St) = 170 GPa, c (St) = 80 GPa, ρ (St) = 7 . for steel and c (Ep) =7 .
537 GPa, c (Ep) = 1 .
482 GPa, ρ (Ep) = 1 .
142 g/cm for epoxy. This implies c l (St) = 4 . µ s, c t (St) = 3 .
22 mm/ µ s and c l (Ep) = 2 .
57 mm/ µ s, c t (Ep) = 1 .
14 mm/ µ s for the longitudinal and transversespeeds. The number of Fourier modes is (2 N + 1) for the PWE method and (2 N + 1) for the MM method;we performed the calculations for N = 0 , , c l and c t in the principal direction as functions of the volume fractionof steel inclusions (Fig. 1b). The curves calculated by the PWE method are plotted by thin lines (lightblue and red online), the curves calculated by the MM method are plotted by thick lines (dark blue andred 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 [23] are also plotted(the Hashin-Shtrikman upper bounds lie far above the other curves and are not displayed). It is observedfrom Fig. 1b that the results of both methods monotonically converge from above to the exact value withgrowing N in agreement with the general statement of §
4. What is significant is that the convergence of theMM 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 c l = 1 (cid:104) ρ (cid:105) (cid:68) (cid:104) c (cid:105) − (cid:69) − , c t = 1 (cid:104) ρ (cid:105) (cid:68) (cid:104) c (cid:105) − (cid:69) − , (73)provides a much better estimate for c l and c t at f > . N = 5 , i.e. withmatrices of about 4000 × f → c eff11 (cid:47) (cid:16) c (St) + 1 − f / c (Ep) (cid:17) − . (74)At the same time the geometry of the unit cell for f → c eff11 can be estimated byan equivalent medium stratified in the 1 − direction, for which the uniaxial strain assumption with constantstress σ leads to the approximation c eff11 ≈ (cid:104) c − (cid:105) − , the same as the right hand side of (74) (as f → − f / tends to zero faster than 1 − f as f →
1, explains the exceptionalaccuracy of the new bound as an estimate for the moduli. Note that, by comparison, the Hashin-Shtrikmanbounds (upper or lower) do not provide a useful estimate in this case.a) b)
Figure 2: A cubic lattice of Steel spheroids in Epoxy matrix. (a) The inclusions are oblate spheroids with minor axis a andunit major axes. (b) The periodic structure for a = 0 . a = 1 (spheres). The second example considers a cubic lattice of spheroidal steel inclusions in epoxy matrix. The shapeof the inclusion evolves from formally a disk of unit diameter to a ball of unit diameter (that is, inscribedin a cubic unit cell) by means of elongating the radius along the x direction, see Fig. 2. We describe thedependence of the effective speeds c l and c t along x and of the corresponding effective elastic moduli c and c on the shape of the spheroidal inclusion. Fourier coefficients for the PWE and MM methods aregiven in Appendix 2. The results are obtained by the PWE method with N = 3 and N = 5 (open circlesin 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 x x plane where the inclusions have circular cross-section and performs direct numerical integration of theRiccati equation (see § x where the shape is ’distorted’. We observe an interestingfeature of a drastic increase of the effective longitudinal speed c l and of the modulus c when the inclusionstend to touch each other, see Fig. 3a. This type of configuration where the inclusions are almost touchingis known to be particularly difficult for numerical calculation of the effective properties [10]. MM appearsto be particularly well suited to treating such problems with closely spaced inclusions since it explicitlyaccounts for the thin gap region via integration of the Riccati equation. The PWE method, on the otherhand, clearly fails to capture the sharp increase in wave speed at N = 5 (matrix size ≈ × N = 3 does not even satisfy the strict upper bounds (73) derived from MM at N = 0.13 c c l (Ep) c t (Ep) c t (St) c l (St) cl ct N =5 N =3 N =0MMPWE N =3 , a) ac jj c (Ep) c (Ep) c c N =5 N =3 N =0MMPWE N =3 , b) Figure 3: (a) Effective wave speeds for the periodic structure of Fig. 2 as a function of the spheroid minor axis a calculated bythe PWE and MM methods. (b) The corresponding elastic moduli.
6. Conclusion
The PWE and MM methods of calculating quasistatic effective speeds in three-dimensional phononiccrystals have been formulated and compared. The MM method can be viewed as a two-dimensional PWEcombined with a one-dimensional propagator matrix approach. The propagator part of the MM scheme iscalculated by numerical integration of a (nonlinear) Riccati differential equation to produce the monodromymatrix.It was shown both analytically and numerically that the MM method provides more accurate approxi-mations than the PWE scheme. In particular, the closed form MM bounds (70) (see also (73)) using onlyone Fourier mode to estimate the effective speed gives better approximations than PWE bounds using morethan a thousand (eleven in each of x i , i = 1 , ,
3) Fourier modes in the case of densely packed structures(see Fig. 1b for f > . N + 1) × N + 1) , requiring a considerable amount of computermemory even for small N . By contrast, the MM scheme uses matrices of dimension 6(2 N + 1) × N + 1) .The reduced memory requirement for the MM method is at the cost of the computer time needed to solve theRiccati equation, a relatively small price to pay. In fact, the ability to set the step size in the Runge-Kuttascheme enables the MM method to efficiently and accurately solve configurations for which the PWE isparticularly ill-suited, such as narrow gaps (see Fig. 1b for f →
1) and closely spaced inclusions (Fig. 3 for a → ppendix Appendix 1. Alternative derivation of Eq. (25) for C eff ll . Let k be parallel to one of the translation vectors. Take the latter to be e = ( δ i ) and so k = ( k , , . Equation (3) may be rewritten in the form η (cid:48) = ( Q + ω Q ) η where (cid:48) ≡ ∂ , η = (cid:18) vC p ∂ p v (cid:19) , Q = (cid:18) − ρ I (cid:19) , (75)while Q is defined in (18) and (19) with j = 1 and a, b = 2 , . Denote η ( x ) ≡ η ( x , x , x ). The solutionto (75) with some initial function η (0) can be written via the matricant in the form η ( x ) = M ( x ) η (0) with M ( x ) = (cid:100)(cid:90) x ( I + ( Q + ω Q ) dx ) . (76)Taking into account assumed 1-periodicity in x and hence the Floquet condition v = e ik x u for thesolution of (3) implies that the solution η of (75) must satisfy η (1) = e ik η (0). Thus, with reference to (76), ω ( k , , ≡ ω ( k ) is an eigenvalue of (3) iff there exists w ≡ w ( x , x ) such that M (1) w = e ik w . (77)Consider asymptotic expansion of (77) in small ω, k . By (76), M (1) = M + ω M + O ( ω ) where M ≡ M [1 , , M [ b, a ] = (cid:98)(cid:82) ba ( I + Q d x ) , M = (cid:90) M [1 , x ] Q M [ x ,
0] d x . (78)The identity M W = W with the 6 × W = ( I ) + (see (21)) implies triple multiplicity of thezero-order ω = 0. Therefore we may write ω α ( k ) = c α k + O ( k ) , w α = w α + k w α + k w α + O ( k ) with w α = W u α , (79)where α = 1 , , u α are some constant linear independent 3 × e ik = 1 + ik − k + O ( k ) in (77) and equating the terms of the same order in k yields1 : M w α = w α , (80a) k : M w α = i w α + w α , (80b) k : M w α + c α M w α = − w α + i w α + w α . (80c)Express w α from (80b) and substitute it in (80c), then scalar multiply the latter by the 6 × (cid:102) W = ( ) + satisfying the identity M +0 [ b, a ] (cid:102) W = (cid:102) W (see (21)). As a result, we obtain C eff11 = (cid:104) (cid:102) W +0 ( M − I ) − W (cid:105) for κ = e = ( δ i ) . (81)It is seen that (25) with M (1) ≡ M and l = 1 is the same as (81), QED. Appendix 2. Fourier coefficients for spheroidal inclusions
The coefficients for the spheroids of Fig. 2a are as follows:
1. MM method.
Identity (27) yields (cid:98) C pq ( g , g , x ) = (cid:40) C pq (Ep) , x (cid:54)∈ (cid:2) − a , a (cid:3) , C pq (Ep) + ( C pq (St) − C pq (Ep)) ˆ χ ( g , g , x ) , x ∈ (cid:2) − a , a (cid:3) χ ( g , g , x ) = ( − g + g RJ (2 πR (cid:112) g + g ) (cid:112) g + g , R = 1 − (2 x − a , where J is the first order Bessel function.
2. PWE method.
Identity (13) yields (cid:98) C pq ( g ) = C pq (Ep) δ g0 + ( C pq (St) − C pq (Ep)) ˆ χ ( g ) , where ˆ χ ( g ) = a ( − g + g + g π | g a | (sin( π | g a | ) − π | g a | cos( π | g a | )) , | g a | = (cid:113) ( ag ) + g + g and δ is a Kronecker symbol. Acknowledgment
A.A.K. acknowledges support from Mairie de Bordeaux. A.N.N. acknowledges support from Institut deM´ecanique et d’Ing´enierie, Universit´e de Bordeaux.
References [1] A. A. Krokhin, J. Arriaga, and L. N. Gumen. Speed of sound in periodic elastic composites.
Phys. Rev. Lett. ,91(26):264302+, 2003.[2] Q. Ni and J. Cheng. Anisotropy of effective velocity for elastic wave propagation in two-dimensional phononic crystals atlow frequencies.
Phys. Rev. B , 72:014305, 2005.[3] A. A. Kutsenko, A. L. Shuvalov, A. N. Norris, and O. Poncelet. On the effective shear speed in 2D phononic crystals.
Phys. Rev. B , 84:064305, 2011.[4] W. J. Parnell and I. D Abrahams. Homogenization for wave propagation in periodic fibre-reinforced media with complexmicrostructure. I - Theory.
J. Mech. Phys. Solids , 56:2521–2540, 2008.[5] I.V. Andrianov, J. Awrejcewicz, V.V. Danishevs’kyy, and D. Weichert. Higher order asymptotic homogenization and wavepropagation in periodic composite materials.
J. Comput. Nonlinear Dynam. , 6:011015, 2011.[6] J. Mei, Z. Liu, W. Wen, and P. Sheng. Effective dynamic mass density of composites.
Phys. Rev. B , 76(13):134205+,2007.[7] D. Torrent and J. S´anchez-Dehesa. Anisotropic mass density by two-dimensional acoustic metamaterials.
New J. Phys. ,10(2):023004+, 2008.[8] Y. Wu and Z.-Q. Zhan. Dispersion relations and their symmetry properties of electromagnetic and elastic metamaterialsin two dimensions.
Phys. Rev. B , 79:195111+, 2009.[9] A. A. Kutsenko, A. L. Shuvalov, and A. N. Norris. Evaluation of the effective speed of sound in phononic crystals by themonodromy matrix method.
J. Acoust. Soc. Am. , 130:3553–3557, 2011.[10] K. C. Nunan and J. B. Keller. Effective elasticity tensor of a periodic composite.
J. Mech. Phys. Solids , 32(4):259 – 280,1984.[11] S. Nemat-Nasser and M. Taya. On effective moduli of an elastic body containing periodically distributed voids.
Q. Appl.Math. , 39:43–59, 1981.[12] S. Nemat-Nasser, T. Iwakuma, and M. Hejazi. On composites with periodic structure.
Mech. Mater. , 1:239–267, 1982.[13] A. S. Sangani and W. Lu. Elastic coefficients of composites containing spherical inclusions in a periodic array.
J. Mech.Phys. Solids , 35:1–21, 1987.[14] V. I. Kushch. Computation of the effective elastic moduli of a granular composite material of regular structure.
Sov. Appl.Mech. , 23:362–365, 1987.[15] V. I. Kushch. Microstresses and effective elastic moduli of a solid reinforced by periodically distributed spheroidal particles.
Int. J. Solids Struct. , 34:1353–1366, 1997.[16] Q. Ni and J. Cheng. Long wavelength propagation of elastic waves in three-dimensional periodic solid-solid media.
J.Appl. Phys. , 101:073515, 2007.[17] G. W. Milton.
The Theory of Composites . Cambridge University Press, 1st edition, 2001.[18] A. N. Norris, A. L. Shuvalov, and A. A. Kutsenko. Analytical formulation of 3D dynamic homogenization for periodicelastic systems.
Proc. R. Soc. A , doi:10.1098/rspa.2011.0698, 2012.[19] A. N. Norris. Elastic moduli approximation of higher symmetry for the acoustical properties of an anisotropic material.
J. Acoust. Soc. Am. , 119:2114–2121, 2006.[20] W.C. Van Buskirk, S. C. Cowin, and R. Carter, Jr. A theory of acoustic measurement of the elastic constants of a generalanisotropic solid.
J. Mater. Sci. , 21:2759–2762, 1986.[21] A. N. Norris. On the acoustic determination of the elastic moduli of anisotropic solids and acoustic conditions for theexistence of planes of symmetry.
Q. J. Mech. Appl. Math. , 42:413–426, 1989.
22] A. A. Kutsenko, A. L. Shuvalov, and A. N. Norris. Converging bounds for the effective shear speed in 2D phononiccrystals.
J. Elasticity , doi:10.1007/s10659-012-9417-y, 2012.[23] Z. Hashin and S. Shtrikman. A variational approach to the elastic behavior of multiphase minerals.
J. Mech. Phys. Solids ,11:127–140, 1963.,11:127–140, 1963.