Effective speed of sound in phononic crystals
EEffective speed of sound in phononic crystals
A.A. Kutsenko a , A.L. Shuvalov a , A.N. Norris a,b a Universit´e de Bordeaux, Institut de M´ecanique et d’Ing´enierie de Bordeaux,UMR 5295, Talence 33405, France, b Mechanical and Aerospace Engineering,Rutgers University, Piscataway, NJ 08854, USA
A new formula for the effective quasistatic speed of sound c in 2D and 3D periodic materials isreported. The approach uses a monodromy-matrix operator to enable direct integration in one ofthe coordinates and exponentially fast convergence in others. As a result, the solution for c has amore closed form than previous formulas. It significantly improves the efficiency and accuracy ofevaluating c for high-contrast composites as demonstrated by a 2D example with extreme behavior. PACS numbers: 62.65.+k, 43.20.+g, 02.70.Hm, 43.90.+v
I. INTRODUCTION
Long-standing interest in modelling effective elasticproperties of composites with microstructure has sub-stantially intensified with the emerging possibility of de-signing periodic structures in air and in solids to formphononic crystals and other exotic metamaterials, whichopen up exciting application prospects ranging from neg-ative index lenses to small scale multiband phononicdevices . This new prospective brings about the need forfast and accurate computational schemes to test ideas insilico . The most common numerical tool is the Fourieror plane-wave expansion method (PWE). It is widelyused for calculating various spectral parameters includ-ing the effective quasistatic speed of sound in acoustic and elastic phononic crystals. At the same time, thePWE calculation is known to face problems when ap-plied to high-contrast composites , which are of especialinterest for applications. Particularly riveting is the casewhere a soft ingredient is embedded in a way breakingthe connectivity of densely packed regions of stiff ingre-dient. Physically speaking, the speed of sound, which islarge in a homogeneously stiff medium, should fall dra-matically when even a small amount of soft componentforms a ’quasi-insulating network’. Note that this case,which implies a strong effect of multiple interactions, isalso ungainly for the multiple-scattering approach .The purpose of present Letter is to highlight a newmethod for evaluating the quasistatic effective soundspeed c in 2D and 3D phononic crystals. The idea isto recast the wave equation as a 1st-order ’ordinary’ dif-ferential system (ODS) with respect to one coordinate(say x ) and to use a monodromy-matrix operator de-fined as a multiplicative (or path) integral in x . Bythis means, we derive a formula for c whose essentialadvantages are an explicit integration in x and an expo-nentially small error of truncation in other coordinate(s).Both these features of the analytical result are shown tosignificantly improve the efficiency and accuracy of itsnumerical implementation in comparison with the con-ventional PWE calculation, which is demonstrated fora 2D steel/epoxy square lattice. The power of the newapproach is especially apparent at high concentration f of steel inclusions, where the effective speed c displays asteep, near vertical, dependence for f ≈
1, a feature notcaptured by conventional techniques like PWE.
II. EFFECTIVE SPEED: 2D ACOUSTIC WAVES
A. SETUP.
Consider the scalar wave equation ∇∇∇ · ( µ ∇∇∇ v ) = − ρω v, (1)for time-harmonic shear displacement v ( x , t ) = v ( x ) e − iωt in a 2D solid continuum with T -periodic density ρ ( x )and shear coefficient µ ( x ). Assume a square unit cell T = (cid:80) i t i a i = [0 , with unit translation vectors a ⊥ a taken as the basis for x = (cid:80) i x i a i . Imposing the Floquetcondition v ( x ) = u ( x ) e i k · x where u ( x ) is periodic and k = kκκκ ( | κκκ | = 1), Eq. (1) becomes( C + C + C ) u = ρω u with C u = −∇∇∇ ( µ ∇∇∇ u ) , C u = − i k · ( µ ∇∇∇ u + ∇∇∇ ( µu )) , C u = k µu. (2)Regular perturbation theory applied to (2) yields the ef-fective speed c ( κκκ ) = lim ω,k → ω ( k ) /k in the form c ( κκκ ) = µ eff ( κκκ ) / (cid:104) ρ (cid:105) , µ eff ( κκκ ) = (cid:104) µ (cid:105) − M ( κκκ ) with (3) M ( κκκ ) = (cid:80) i,j =1 M ij κ i κ j , M ij = (cid:0) C − ∂ i µ, ∂ j µ (cid:1) = M ji , where ∂ i ≡ ∂/∂x i , spatial averages are defined by (cid:104) f (cid:105) ≡ (cid:82) T f ( x )d x (cid:0) = (cid:104)(cid:104) f (cid:105) (cid:105) , (cid:104) f (cid:105) i ≡ (cid:82) f ( x )d x i (cid:1) , (4)and ( · , · ) denotes the scalar product in L ( T ) so that( f, h ) = (cid:104) f h ∗ (cid:105) ( ∗ means complex conjugation). The dif-ficulty with (3) is that it involves the inverse of a par-tial differential operator C . One solution is to apply adouble Fourier expansion to C − and ∂ i µ in (3). Thisleads to the PWE formula for the effective speed whichis expressed via infinite vectors and the inverse of theinfinite matrix of Fourier coefficients of µ ( x ). Numeri-cal implementation of the PWE formula requires dealingwith large dense matrices, especially in the case of high-contrast composites for which the PWE convergence isslow (see § IV). An alternative ”brute force” procedure of a r X i v : . [ m a t h - ph ] J un the scaling approach is to numerically solve the partialdifferential equation C h = ∂ i µ for the -periodic func-tion h ( x ) (e.g. via the boundary integral method ).The new approach proposed here leads to a more effi-cient formula for c based on direct analytical integrationin one coordinate direction. There are two ways of doingso. The first proceeds from the ODS form of the waveequation (1) itself, which means ’skipping’ (3). This isconvenient for deriving c ( κκκ ) in the principal directions κκκ (cid:107) a , , see § IIB. The second method is more closelyrelated to the conventional PWE and scaling approachesin that it also starts from (3) but treats it differently,namely, the equation C h = ∂ i µ is cast in ODS form andanalytically integrated in one coordinate. This is basi-cally equivalent to the former method, but enables aneasier derivation of the off-diagonal component M forthe anisotropic case, see § IIC.
B. Wave speed in the principal directions.
Thewave equation (1) may be recast as ηηη (cid:48) = Q ηηη with A = − ∂ ( µ∂ ) , Q = (cid:18) µ − A − ρω (cid:19) , ηηη ( x ) = (cid:18) vµv (cid:48) (cid:19) , (5)where (cid:48) stands for ∂ . The solution to Eq. (5) for initialdata ηηη (0 , x ) ≡ ηηη (0 , · ) at x = 0 is ηηη ( x , · ) = M [ x , ηηη (0 , · ) with M [ a, b ] = (cid:98)(cid:82) ab ( I + Q d x ) . (6)The operator M [ x ,
0] is formally the matricant, or prop-agator, of (5) defined through the multiplicative integral (cid:98)(cid:82) (with I denoting the identity operator). It is assumedfor the moment that ρ ( x ) and µ ( x ) are smooth to en-sure the existence of M . The matricant over a period, M [1 , k = ( k T so that v ( x ) = u ( x ) e ik x and ηηη (1 , · ) = ηηη (0 , · ) e ik . By (6) , this implies the eigenproblem M [1 , w ( k ) = e ik w ( k ) . (7)where M [1 ,
0] depends on ω . Eq. (7) defines k = k ( ω )and hence ω = ω ( k ) , where ω is the eigenvalue of (1)with v ( x ) = u ( x ) e ik x . The effective speed c ( κ ) =lim ω,k → ω/k can therefore be determined by applyingperturbation theory to (7) as ω, k →
0. The asymptoticform of M [1 ,
0] follows from definitions (5) and (6) as M [1 ,
0] = M + ω M + O ( ω ) where : M ≡ M [1 , , M [ a, b ] = (cid:98)(cid:82) ab ( I + Q d x ) with Q ≡ Q ω =0 = (cid:18) µ − A (cid:19) , (8) M = (cid:90) M [1 , x ] (cid:18) − ρ (cid:19) M [ x ,
0] d x . Note the identities Q w = , Q +0 (cid:101) w = and hence M [ a, b ] w = w , M +0 [ a, b ] (cid:101) w = (cid:101) w ( ∀ a, b )for w = (1 0) T , (cid:101) w = (0 1) T . (9)By (9) w is an eigenvector of M with the eigenvalue 1,and it can be shown to be a single eigenvector. Therefore w ( k ) = w + k w + k w + O ( k ) and ω = ck + O ( k ).Insert these expansions along with (8) in (7) and collectthe first-order terms in k to obtain M w = w + i w ⇒ w = i ( M − I ) − w . (10)According to (9), M − I has no inverse but is a one-to-one mapping from the subspace orthogonal to w ontothe subspace orthogonal (cid:101) w ; hence, w exists and (cid:101) w · w is uniquely defined. The terms of second-order in k in(7) then imply M w + c M w = i w + w . (11)Scalar multiplication on both sides by (cid:101) w leads, with ac-count for (9) and (8) , to c (cid:104) ρ (cid:105) = − i (cid:104) (cid:101) w · w (cid:105) , whenceby (10 ) c ( κ ) = (cid:104) ρ (cid:105) − (cid:10) (cid:101) w · ( M − I ) − w (cid:11) , (12)where the notation (cid:104)·(cid:105) is explained in (4). Interchangingvariables x (cid:29) x in the above derivation yields a similarresult for c ( κ ) as follows c ( κ ) = (cid:104) ρ (cid:105) − (cid:68) (cid:101) w · ( (cid:102) M − I ) − w (cid:69) where (cid:102) M = (cid:98)(cid:82) ( I + (cid:101) Q d x ) , (13) (cid:101) Q = (cid:18) µ − (cid:101) A (cid:19) , (cid:101) A = − ∂ µ∂ . The result for a rectangular lattice with T = [0 , T ] × [0 , T ] is obtained by replacing x i with x i /T i . C. The full matrix M ij . The anisotropy of theeffective speed c ( κκκ ) , i.e. its dependence on the wavenormal κκκ ≡ k /k , is determined by the quadratic form M ( κκκ ) = (cid:80) i,j =1 M ij κ i κ j (see Eq. (3)), and representedby the ellipse of (squared) slowness c − ( κ ). Eqs. (12)and (13) , which define c ( κ i ) and so M ii , suffice forthe case where T is rectangular and µ ( x ) is even in (atleast) one of x i so that the effective-slowness ellipse is c − ( κκκ ) = (cid:80) i =1 , c − ( κ i ) κ i with the principal axes par-allel to a ⊥ a . Otherwise c ( κκκ ) for arbitrary κκκ requiresfinding the off-diagonal component M . For this pur-pose, with reference to (3), consider the equation C h = ∂ µ (14)for -periodic h ( x ). With the above notations this canbe written as − ( µh (cid:48) ) (cid:48) + A h = µ (cid:48) or, more conveniently,( µ (cid:101) h (cid:48) ) (cid:48) = A (cid:101) h with (cid:101) h = h + x . The latter is equivalent to ξξξ (cid:48) = Q ξξξ where ξξξ = (cid:18) h + x µ ( h (cid:48) + 1) (cid:19) (15)and Q is given in (8) . The general solution to (15) is ξξξ ( x , · ) = M [ x , ξξξ (0 , · ) , (16)where M [ x ,
0] is defined in (8) , and ξξξ (0 , · ) is the initialdata at x = 0. The periodicity of h implies ξξξ (1 , · ) = ξξξ (0 , · ) + w , while ξξξ (1 , · ) = M ξξξ (0 , · ) by (16). Hence ξξξ (0 , · ) = ( M − I ) − w and so (14) is solved by ξξξ ( x , · ) = M [ x ,
0] ( M − I ) − w . (17)Substituting (17) into the definition of M in (3) yields M = ( C − ∂ µ, ∂ µ ) = (cid:104) h∂ µ (cid:105) = (cid:104) ∂ µ w · ξξξ (cid:105) = (cid:10) ∂ µ w · M [ x , M − I ) − w (cid:11) . (18)Note that the formula (18) for M requires more com-putation than the formulas (12) and (13) for M ii . Inter-estingly, if the unit cell T is square, then, for an arbitrary(periodic) µ ( x ), Eq. (18) can be circumvented by usingthe identity M = ( (cid:102) M − (cid:102) M ) /
2, where (cid:102) M ii follow fromEqs. (12) and (13) applied to the square lattice obtainedfrom the given one by turning it 45 ◦ . D. Discussion.
The two lines of attack outlinedmen-tioned in § II.A are equivalent in that the formula (12)for the effective speed c ( κ ) in the principal directioncan also be inferred from Eq. (3). Inserting the solution(17) of (14) defines the component M as M = (cid:0) C − ∂ µ, ∂ µ (cid:1) = (cid:104) hµ (cid:48) (cid:105) = (cid:104) µ (cid:48) w · ξξξ (cid:105)−(cid:104) x µ (cid:48) (cid:105) . (19)Integrating by parts each term in the last identity andusing the periodicity of µ ( x ) along with Eqs. (8) , (9),(15)-(17) (see also the notation (4)) yields (cid:104) µ (cid:48) w · ξξξ (cid:105) = −(cid:104) (cid:101) w · ( M − I ) − w (cid:105) + (cid:104) µ (0 , x ) (cid:105) , − (cid:104) x µ (cid:48) (cid:105) = (cid:104) µ (cid:105) − (cid:104) µ (1 , x ) (cid:105) = (cid:104) µ (cid:105) − (cid:104) µ (0 , x ) (cid:105) . (20)Thus, M = (cid:104) µ (cid:105)− (cid:10) (cid:101) w · ( M − I ) − w (cid:11) which leads to(12), QED. Note that Eq. (18) is also obtainable via themonodromy matrix of the wave equation (1) (the ap-proach of § IIB) with v ( x ) = u ( x ) e i k · x and k ∦ a i , but thismethod of derivation of M is lengthier than in § IIC.As another remark, it is instructive to recover a knownresult for the case where µ ( x ) is periodic in one coordi-nate and does not depend on the other, say µ ( x , x ) = µ ( x ). Using (8) , (8) and (13) gives( M − I ) (cid:18) (cid:104) µ − (cid:105) − (cid:19) = w , ( (cid:102) M − I ) (cid:18) µ ( x ) (cid:19) = w . (21)Therefore, by (12) and (13) , c ( κ ) = (cid:104) µ − (cid:105) − / (cid:104) ρ (cid:105) and c ( κ ) = (cid:104) µ (cid:105) / (cid:104) ρ (cid:105) while M = 0 by (18) with ∂ µ = 0.Finally, we note that, while the above evaluation ofquasistatic speed c is exact, using the same monodromy-matrix approach also provides a closed-form approxima-tion of c. For the isotropic case, it is as follows (see formore details): c ≈ (cid:104) ρ (cid:105) (cid:18)(cid:68)(cid:10) µ − (cid:11) − (cid:69) + (cid:68) (cid:104) µ (cid:105) − (cid:69) − (cid:19) . (22) III. EFFECTIVE SPEEDS IN PRINCIPALDIRECTIONS FOR 3D ELASTIC WAVES
The equation for time-harmonic elastic wave motion v ( x , t )= v ( x ) e − iωt is, with repeated suffices summed, − ∂ j ( c ijkl ∂ l v k ) = ρω v i ( i, j, k, l = 1 , , , (23)where density ρ ( x ) and compliances c ijkl ( x ) are T -periodic in a 3D periodic medium. Assume a cubic unitcell T = (cid:80) i t i a i = [0 , and refer the components x i , v i and c ijkl to the orthogonal basis formed by the transla-tion vectors a i . Impose the condition v ( x ) = u ( x ) e i k · x with periodic u ( x ) = ( u i ) and take k parallel to one of a i , e.g. to a . Eq. (23) may be rewritten in the form ηηη (cid:48) = Q ηηη with ηηη ( x ) = (cid:18) ( u i )( c i kl ∂ l u k ) (cid:19) , Q = (cid:18) −C − A C − ω ρδ ij + A − A C − A A C − (cid:19) (24)where the self-adjoint matrix operators C and A , are C = ( c i k ) , A ( u i ) = ( c i ka ∂ a u k ) , A ( u i ) = ( ∂ a ( c iakb ∂ b u k )) with a, b = 2 , . (25)Like in the 2D case, denote the monodromy matrix for(24) at ω = 0 by M = (cid:98)(cid:82) ( I + Q dx ) where Q = Q ω =0 ,and also introduce the 6 × W = ( δ ij ) T and (cid:102) W = ( δ ij ) T . Reasoning similar to that in § II.C leadsus to the conclusion that the effective speeds c α ( κ ) =lim ω,k → ω/k ( α = 1 , ,
3) of the three waves with k ≡ kκκκ parallel to a are the eigenvalues of the 3 × (cid:68)(cid:68) (cid:102) W · ( M − I ) − W (cid:69) (cid:69) (cid:0) with (cid:104)·(cid:105) i ≡ (4) (cid:1) . (26) IV. NUMERICAL IMPLEMENTATION
There are several ways to use the above analytical re-sults for calculating the effective speed. One approachis to transform to Fourier space with respect to coordi-nate(s) other than the coordinate of integration in themonodromy matrix. Consider the 2D case and apply theFourier expansion f ( x , x ) = (cid:80) n ∈ Z (cid:98) f n ( x ) e πinx in x for the functions f = µ and µ − . Then the operator ofmultiplying by the function µ − ( x , · ) and the differentialoperator A ( x )= − ∂ ( µ ( x , · ) ∂ ) become matrices µ − (cid:55)−→ µµµ − ( x ) = ( (cid:100) µ − n − m ) = ( (cid:98) µ n − m ) − , A (cid:55)−→ A ( x ) = 4 π ( nm (cid:98) µ n − m ) , n, m ∈ Z , (27)and Eq. (12) reduces to following form c ( κ ) = (cid:104) ρ (cid:105) − (cid:101) w (cid:98) · ( M − I ) − w (cid:98) with M = (cid:98)(cid:82) ( I + Q dx ) , Q ( x ) = (cid:18) µµµ − A 0 (cid:19) , (28) (cid:101) w (cid:98) = ( δ n ) T , w (cid:98) = ( δ n ) T , where c ( κ ) = c =const. for any κκκ in the isotropic case.The above vectors and matrices are, strictly speaking, ofinfinite dimension, which needs to be truncated for nu-merical purposes. In this sense there is no loss of general-ity in assuming a smooth µ ( x ) in the course of derivationsin § II. Implementation of Eq. (28) consists of two steps. fc c Ep c St EpSt c MM num.: N = 1 N = 10 N = 20 c PWE num.: N = 1 N = 10 N = 20 c MM ≈ (22): FIG. 1: Effective speed c versus concentration f of square rodsfor 2D St/Ep and Ep/St lattices (see details in the text). Step 1.
Calculate the multiplicative integral (28) defin-ing M . For an arbitrary µ ( x ), one way is to use adiscretization scheme. Divide the segment x ∈ [0 , N intervals [ x ( i )1 , x ( i +1)1 ) ≡ ∆ i , i = 1 ..N , of smallenough length. Calculate 2 N + 1 Fourier coefficients (cid:98) µ n ( x ( i )1 ), n = − N..N and the (2 N +1) × (2 N +1) matrices Q ( x ( i )1 ) for each i = 1 ..N , and then use the approximateformula M = (cid:81) i = N exp (cid:104) ∆ i Q ( x ( i )1 ) (cid:105) . Recall that (cid:98)(cid:82) satisfies the chain rule and is exactly equal to exp(∆ i Q )for x ∈ ∆ i if µ ( x ) does not depend on x within ∆ i .Therefore the calculation is much simpler in the com-mon case of a piecewise homogeneous unit cell with onlya few inclusions of simple shape (see the example below). Step 2.
Solve the system ( M − I ) w (cid:98) = i w (cid:98) for unknown w (cid:98) . First remove one zero row and one zero column inthe matrix M − I (see the remark below (10)). Then thevector w (cid:98) is uniquely defined and may be found by anystandard method. Note that only a single component of h is needed to evaluate (cid:101) w (cid:98) · w (cid:98) . Finally dividing by (cid:104) ρ (cid:105) yields the desired result (28) .As an example, we calculate the effective shear-wavespeed c versus the volume fraction f of square rods pe-riodically embedded in a matrix material forming a 2Dsquare lattice with translations parallel to the inclusionedges. A high-contrast pair of materials is chosen suchas steel ( ≡ St, with ρ = 7 . kg/m , µ = 80 GPa) andepoxy ( ≡ Ep, with ρ = 1 .
14 10 kg/m , µ = 1 .
48 GPa).We consider two conjugated St/Ep and Ep/St configu-rations, where the matrix and rod materials are eitherSt and Ep or Ep and St, respectively. The results aredisplayed in Fig. 1. The curves c MM ( f ) are computed by the present monodromy-matrix (MM) method, Eq.(28) , they are complemented by the approximation (22).Also shown for comparison are the curves c PWE ( f ) com-puted from the truncated formula of the conventionalPWE method based on a 2D Fourier transform of (3).Calculations are performed for a different fixed number2 N + 1 ≡ d of the 1D Fourier coefficients of µ ( x ), whichimplies 2 d × d monodromy matrix in (28) and, by con-trast, d × d matrix in the PWE formula . Apart fromthis advantage of the MM calculation, it is also seen tobe remarkably more stable - with a reasonable fit pro-vided already at N = 1. The difference between the MMand PWE numerical curves is especially notable for thecase of densely packed steel rods. Interestingly, the MMcomputation and estimate both predict a steep fall for c ( f ) when a small concentration 1 − f of epoxy forms a’quasi-insulating network’. The PWE fails to capture thisimportant physical feature for reasons described next.The far superior stability and accuracy of the MMmethod observed in Fig. 1 can be explained as fol-lows. The PWE formula implies calculating M ≈ (cid:80) | g | This work has been supportedby the grant ANR-08-BLAN-0101-01 and the projectSAMM. A.N.N. acknowledges support from the CNRS.———— D. Torrent, A. H˚akansson, F. Cervera, and J. S´anchez-Dehesa, Phys. Rev. Lett. , 204302 (2006); D. Torrent, J.S´anchez-Dehesa, Phys. Rev. B 74, 224305 (2006); J. Phys.9, 323 (2007); ibid. 10, 023004 (2008). J. Mei, Z. Liu, W. Wen, P. Sheng, Phys. Rev. Lett. ,024301 (2006); Phys. Rev. B , 134205 (2007). J.O. Vasseur, P.A. Deymier, B. Djafari-Rouhani, Y. Pennec,A.-C. Hladky-Hennion, Phys. Rev. B , 085415 (2008); Y.Pennec, J.O. Vasseur,B. Djafari-Rouhani, L. Dobrzynski,P.A. Deymier, Surf. Sci. Rep. , 229-291 (2010). A. A. Krokhin, J. Arriaga, L. N. Gumen, Phys. Rev. Lett. , 264302 (2003). Q. Ni and J. Cheng, J. Appl. Phys. , 073515 (2007). A. A. Kutsenko, A. L. Shuvalov, A. N. Norris, O. Poncelet,submitted. I.V. Andrianov, J. Awrejcewicz, V.V. Danishevs’kyy, D.Weichert, J. Comput. Nonlinear Dynam. , 011015 (2011) The subsequent results are equally valid for acoustic wavesin fluid-like phononic crystals under the standard inter-change of ρ and µ for solids by K − and ρ −1