Phononic band structures and stability analysis using radial basis function method with consideration of different interface models

Physica B 489 (2016) 1–11

Contents lists available at ScienceDirect

Physica B journal homepage: www.elsevier.com/locate/physb

Phononic band structures and stability analysis using radial basis function method with consideration of different interface models Zhi-zhong Yan a,b,n, Chun-qiu Wei a,b, Hui Zheng c, Chuanzeng Zhang c a

School of Mathematics and Statistics, Beijing Institute of Technology, Beijing 100081, PR China Beijing Key Laboratory on MCAACI, Beijing Institute of Technology, Beijing 100081, PR China c Department of Civil Engineering, University of Siegen, D-57068 Siegen, Germany b

art ic l e i nf o

a b s t r a c t

Article history: Received 13 November 2015 Received in revised form 3 February 2016 Accepted 18 February 2016

In this paper, a meshless radial basis function (RBF) collocation method is developed to calculate the phononic band structures taking account of different interface models. The present method is validated by using the analytical results in the case of perfect interfaces. The stability is fully discussed based on the types of RBFs, the shape parameters and the node numbers. And the advantages of the proposed RBF method compared to the finite element method (FEM) are also illustrated. In addition, the influences of the spring-interface model and the three-phase model on the wave band gaps are investigated by comparing with the perfect interfaces. For different interface models, the effects of various interface conditions, length ratios and density ratios on the band gap width are analyzed. The comparison results of the two models show that the weakly bonded interface has a significant effect on the properties of phononic crystals. Besides, the band structures of the spring-interface model have certain similarities and differences with those of the three-phase model. & 2016 Elsevier B.V. All rights reserved.

Keywords: RBF collocation method Phononic crystals Stability Band structures Interface models

1. Introduction Phononic crystals (PCs) which are periodic structures and feature acoustic/elastic wave band gaps have attracted great attentions in the recent twenty years [1]. The wave band gaps can be used in acoustic filters, silent blocks, and so on. Early studies about PCs are mostly focused on the perfectly bonded interfaces. But due to all kinds of reasons, defects and damages are more or less presented at interfaces. Moreover, some specific interfaces may be required in practical applications. Generally speaking, the interface can be considered as a special medium layer which is different from the component materials, namely the “interface layer” [2]. Many existing researches about the composites have shown that different interface conditions may have significant effects on their mechanical behaviors. Thus the interface effect on the wave motion property of PCs becomes an important research topic [3]. Although a lot of studies on composite materials with imperfect interfaces have been carried out over the years [4], few researches about the effects of the interfaces on the phononic band structures have been performed. Recently, Zhu et al. [5] studied the band gap structures considering the interface effects n Corresponding author at: School of Mathematics and Statistics, Beijing Institute of Technology, Beijing 100081, PR China. E-mail address: [email protected] (Z.-z. Yan).

http://dx.doi.org/10.1016/j.physb.2016.02.026 0921-4526/& 2016 Elsevier B.V. All rights reserved.

based on the boundary element method (BEM). Liu et al. [6] developed the Mie scattering matrix and considered the boundary conditions. Zhen et al. [7] used the Dirichlet-to-Neumann map method to calculate the band structures in nano-scale phononic crystals. So far, the interaction between waves and imperfect interfaces in PCs is getting more and more attention in the research and engineering applications. In recent years, several numerical methods have been developed for the calculation of phononic band structures, such as the plane wave expansion (PWE) method [8,9], the multiple scattering theory (MST) method [10,11], the wavelet method [12,13], the finite difference time domain (FDTD) method [14–18] and the finite element method (FEM) [19,20], etc. However, most of these numerical methods regard a phononic crystal as a whole domain, for example, the PWE method, the FDTD method and the wavelet method. In this sense, they are not applicable in dealing with the imperfect interface conditions. The traditional FEM which is used to consider surface/interface effects is not convenient in the mesh generation and re-meshing, and its convergence speed is relatively slow in some ultra-thin structures. Therefore, it is necessary to develop some more efficient and convenient numerical methods to analyze the phononic band structures considering the imperfect interfaces. The RBF methods [21,22] have provided an additional dimension in solving partial differential equations (PDEs) in recent years. Compared to the traditional domain-type methods such as the

2

Z.-z. Yan et al. / Physica B 489 (2016) 1–11

FDTD and the FEM methods, the merits of the RBF methods are mathematically simple, independent of dimensionality and truly meshless. The RBF methods can also avoid tedious mesh generation for the high-dimensional problems involving the moving boundary-value problems, and the insertion of additional nodes is much simpler than the re-meshing in the domain-type methods with respect to the post-processing of the computational data. In this paper, we employ a strong-form RBF method to calculate the band structures of one-dimensional (1D) PCs with consideration of different interface models. Compared with the meshless methods constructed by the weak-form, the present approach appears to be attractive as it avoids domain and boundary integrations. Moreover, the advantages of the strong-form method are its simpler numerical implementation and lower computational cost. This article is organized as follows. In Section 2, the basic approximation idea of the RBF method is described, and the strongform collocation method using the RBFs is developed for phononic crystals. The feasibility, the stability and the efficiency of the RBF method are discussed in Section 3, two kinds of imperfect interface models (spring-interface model, three-phase model) are analyzed in detail. Moreover, the influences of the length ratios and the material density ratios on the band structures are also investigated and illustrated. Then some concluding remarks are given in Section 4.

2. Theoretical method

coefficients that we need to compute. Popular choices for RBFs include multi-quadric (MQ) RBFs

φ(r ) =

r 2 + cs2 ,

Inverse multi-quadric (inverse MQ) RBFs

φ(r ) =

1 r 2 + cs2

(6)

⎛ r2 ⎞ φ(r ) = exp⎜⎜ − 2 ⎟⎟, cs ⎠ ⎝

(7)

where r = ||x − x i|| is the Euclidean distance and cs is the shape parameter used to scale the basic functions. By using Eq. (4), we obtain the discrete boundary value problem n

∑ Ηφ(

xq − xi )αi = ℏ(xq),

xq ∈ Ω

q = 1, 2, ⋯n1,

i=1

(8)

n

∑ Κφ( xs − xi

)αi = ƛ(xs ),

xs ∈ I

s = 1, 2, ⋯n 2,

)αi = ℓ(xt ),

xt ∈ ℵ

t = 1, 2, ⋯n3.

i=1

(9)

n

2.1. Function approximation with RBFs for a boundary value problem In this section, we give a brief introduction to the RBF collocation method for initial-boundary value problems of PDEs. The basic idea (see, e.g., [21] or [22]) is to use a set of RBF basis functions to represent an unknown function (the approximate solution of the PDE) over a set of scattered nodes in a domain Ω . To illustrate this point, we consider the following one-dimensional boundary value problem

x∈Ω

(1)

where Η is the differential operator in Ω , u(x ) is the unknown field quantity, and ℏ(x ) is the source term. Besides, the boundary conditions are as follows

Κu(x) = ƛ(x), x ∈ I

(2)

Ρu(x) = ℓ(x), x ∈ ℵ

(3)

where I is the Dirichlet and ℵ is the Neumann boundary, Κ and Ρ are the differential operators on I and ℵ respectively, ƛ(x ) and ℓ(x ) are known source functions. In the frame of the RBF collocation method, the general solution u(x ) can be approximated as follows

i=1

∑ φ(‖x − xi‖)αi

(4)

i=1

where i is the index of the nodes in Ω , n is the total number of all nodes, φ is the RBF that we choose, and αi are the unknown

(10)

where n = n1 + n2 + n3 with n1, n2 and n3 being the numbers of the collocation nodes. From the Eqs. (8)–(10), the unknown field quantities and other related quantities can be evaluated when the unknown coefficients αi are calculated. 2.2. RBF collocation method for phononic crystals Here, the RBF collocation method as presented in Section 2.1 is applied for the band structure calculations of one-dimensional (1D) phononic crystals by taking into account the existence of imperfect interfaces. We consider a 1D phononic crystal composed of two different materials A and B as shown in Fig. 1. Let the length of the material A and B be d1 and d2, respectively, so the lattice constant is a = d1 + d2 as shown in Fig. 1. The longitudinal waves are assumed to propagate normally along the x-direction. The governing equation can be expressed as

d2uj (x) dx2

+

ρj ω2 λ j + 2μj

uj (x) = 0

(11)

where uj (x ) is the displacement in the jth (j¼1, 2) sublayer, ω is the circular frequency, ρj , λj and μj are the mass density, the Lamé constant and the shear modulus, respectively, clj =

n

u(x) =

,

Gaussian RBFs

∑ Ρφ ( x t − x i

Ηu(x) = ℏ(x),

(5)

(λj + 2μj ) /ρj is

the longitudinal sound velocity, where the quantities with the subscript “ j = 1” are referred to the material A, while those with the subscript “ j = 2” to the material B. In the following, different interface models are considered.

P wave Lattice constant a x Fig. 1. Schematic structure of a one-dimensional phononic crystal.

Z.-z. Yan et al. / Physica B 489 (2016) 1–11

Table 1 Material properties.

u1(0)eika−u2(a) = F⋅σ (0),

Material

ρ (kg/ m3 )

λ (kg/ ms2 )

μ (kg/ ms2 )

Al

2.73 × 103

6.82 × 1010

2.87 × 1010

Epoxy

1.18 × 103

4.43 × 109

1.59 × 109

2.2.1. Perfect interface model The perfect interface model is one of the most simple interface models, i.e., the displacement and stress are continuous at the interface. The boundary conditions are as follows

[u] = 0,

(12)

[σ ] = 0

where [•] represents the jump value across interface, [u] = u+ − u− denotes the displacement discontinuity value, [σ ] = σ + − σ −denotes the stress discontinuity value, u+ , u−, σ + and σ − represent the displacements and stresses at both sides of the interface, respectively. The component form of Eq. (12) can be written as

u1(0)eika = u2(a), u1(d1) = u2(d1),

σ1(0)eika = σ2(a) σ1(d1) = σ2(d1)

(13)

2.2.2. Imperfect interface models Here we mainly consider two kinds of weakly bonded interface models, i.e., the spring-interface model and the three-phase model. 2.2.2.1. Model 1: spring-interface model. The spring-interface model ignores the interface layer thickness and assumes the component materials on two sides of the interface could produce the relative displacement. Thus, it is replaced with the spring interface layer which satisfies the Hooke law and ignores its own quality. The interface boundary conditions can be expressed as

[u] = F • σ ,

[σ ] = 0,

(14)

The component form of Eq. (14) can be represented as follows

u1(d1)−u2(d1) = F⋅σ (d1)

3

σ1(0)eika = σ2(a) σ1(d1) = σ2(d1)

(15)

where F = h/(λ 0 + 2μ0 ) is the spring coefficient with λ 0 and μ0 are the material parameters of the boundary, respectively, h is the thickness of the interface layer [23]. 2.2.2.2. Model 2: three-phase model. The three-phase model regards the interface layer as a single medium, which is still perfectly bonded with the components. The interface conditions are denoted as follows

[u] = F • σ ,

[σ ] = 0,

(16)

The component form can be written as [23]

u+ − u− =

⎞ ⎞ h⎛ 1 1 h⎛ 1 1 ⎜⎜ ⎟⎟σ + − ⎜⎜ ⎟⎟σ − − − λ2 + 2μ2 ⎠ λ1 + 2μ1 ⎠ 2 ⎝ λ 0 + 2μ 0 2 ⎝ λ 0 + 2μ 0

=

⎞ ∂u+ ⎞ ∂u− h ⎛ λ2 + 2μ2 h ⎛ λ1 + 2μ1 ⎜⎜ − 1⎟⎟ + ⎜⎜ − 1⎟⎟ 2 ⎝ λ 0 + 2μ 0 2 ⎝ λ 0 + 2μ 0 ⎠ ∂x ⎠ ∂x

(17)

where λ 0 and μ0 are the average values of the elastic constants of materials Al and Epoxy, respectively, h is the thickness of the interface layer. Without loss of generality, the RBF collocation method for the PCs in the case of perfect interfaces is described. First, the 2n nodes x1, x2, ⋅⋅⋅, x n, x n + 1, ⋯ , x2n are distributed in the interval [0, a], where x1, x n , x n + 1 and x2n are the boundary points, namely, x1 = 0, x n = x n + 1 = d1 and x2n = a, respectively, the step size is 1/(2n − 2). Based on Eq. (4), the displacements u1(x ) and u2(x ) are approximated by using the radial basis functions, which can be expressed as n

u1(x) =

∑ αiφi(x)

x ∈ ⎡⎣ 0, d1⎤⎦

i=1 2n

u2(x) =

∑

βi φi(x)

x ∈ ⎡⎣ d1, a⎤⎦

i = n+ 1

(18)

where φi(x ) = 1/ (csri )2 + 1 (i = 1, 2, ⋯2n ) are the inverse MQ RBFs.

Fig. 2. The band structures using different inverse MQ numbers: (a) n¼19, (b) n¼ 11.

4

Z.-z. Yan et al. / Physica B 489 (2016) 1–11

Fig. 3. Influences of the shape parameter, the basis function number and the node number on the relative error.

Fig. 4. Relative error versus the shape parameter with the MQ function.

Fig. 5. Relative error versus the shape parameter with the Gauss function.

Z.-z. Yan et al. / Physica B 489 (2016) 1–11

5

Fig. 6. Band structures with their respectively best cs values for different kinds of RBFs.

0.1

According to the Bloch theory, the displacement can be written

Relative error

FEM RBF

as

0.01

uj (x + a) = eikau j (x)

1E-3

where k is the Bloch wave vector in the first Brillouin zone (BZ) [−π /a, π /a] of the reciprocal lattice, i = −1 . Because the displacement and the stress in the joints, i.e., x1, x2n are continuous, we have the following conditions

⎧ eikau1(x1) = u2(x2n) ⎪ ⎨ ∂u ∂u ⎪ eikaE1 1 (x1) = E2 2 (x2n) ⎩ ∂x ∂x

1E-4

10

100

Nodes/Degrees of freedom Fig. 7. Convergence rates of the RBF method and the FEM.

(19)

(20)

where E1 and E2 are the Young’s modulus of the materials A and B, respectively. Using the Eqs. (4), (18) and (20), we have

6

Z.-z. Yan et al. / Physica B 489 (2016) 1–11

18 16

FEM RBF

14

Time(s)

12 10 8 6 4 2

20

40

60

80

100

Nodes/Degrees of freedom Fig. 8. Comparisons of the computation time for the RBF method and the FEM. ⎧ n n d 2φi (x v ) ⎪ (λ1 + 2μ1) ∑ α i = −ρω 2 ∑ α iφi (x v ) v = 2, 3, ⋯ , n − 1 ⎪ 2 dx i =1 i =1 ⎪ ⎪ 2n 2 2n ⎪ (λ + 2μ ) ∑ α d φi (x v ) = −ρω 2 ∑ α φ (x ) v = n + 2, ⋯ , 2n − 1 2 i i i v 2 ⎪ 2 dx i =n +1 i =n +1 ⎪ 2n n ⎪ ⎪ ∑ α iφi (x n ) = ∑ α iφi (x n + 1) ⎪ ⎪ i =1 i =n +1 ⎨ 2n n ⎪ dφ (x n + 1) dφi (x n ) ⎪ E1 ∑ α i = E2 ∑ α i i dx dx ⎪ i =1 i =n +1 ⎪ 2n n ⎪ eika ∑ α iφi (x1) = ∑ α iφi (x2n ) ⎪ ⎪ i =n +1 i =1 ⎪ 2n n ⎪ dφ (x1) dφ ( x2 n ) eikaE1 ∑ α i i = E2 ∑ α i i ⎪ ⎪ dx dx ⎩ i =n +1 i =1

Fig. 9. Comparison of the band structures between the spring-interface model and the perfect interface.

(21)

The Eq. (21) can be recast into the following matrix form

⎡ α1⎤ ⎡ α1⎤ M⎢ ⎥ = ω2N⎢ ⎥ ⎣ α 2⎦ ⎣ α 2⎦

(22)

where M and N are n × n matrices, α1 and α2 are n × 1 vectors related to the undetermined coefficients. The specific expressions of M, N, α1 and α2 are given in Appendix. For the imperfect interface models, we follow the same procedure as for the perfect interface and consider the Eq. (19) as well as the imperfect interface conditions, i.e., Eqs. (15) and (17), the generalized eigenvalue equation can be obtained and solved numerically by using the eigensolver in MATLAB.

structures. The filling fraction of the Al/epoxy phononic crystal is f = d1: d2 = 1: 1, the normalized value of the lattice constant is a = 1 and cl2 is the longitudinal wave velocity of epoxy. Fig. 2 shows the band structures obtained by using the RBF method (scattered symbols) and the exact analytical solutions (solid lines) [24,25]. Fig. 2(a) is the results obtained by employing 19 inverse MQ basis functions. From Fig. 2(a), we can find that the RBF results agree with the exact analytical solutions very well. Hence it verifies the effectiveness of the present method developed above. However, the results with 11 basis functions shown in Fig. 2 (b) show that the poor consistency appears from the 9th band to the higher frequency bands, which indicates that the number of basis functions has a significant effect on the stability of the RBF method. Thus, a natural question to ask is, how would the stability of the RBF method be altered if considering different numbers of RBFs? Are there any other influence factors? In the following, we give the detailed discussions about the influences of the basis function numbers, the types of RBFs and the shape parameters. For convenience, we consider the relative error of the 12th band in Fig. 2(a), the relative error δ is defined as follows 21

Σ ω w − ω wexact

δ= 3. Numerical results and discussions In this section, the developed RBF collocation method is applied to the band structure calculations of the phononic crystals with imperfect interfaces. 3.1. Results with the perfect interface and stability analysis of the RBF method In this subsection, the phononic crystal composed of aluminum (Al) and epoxy with perfect interfaces is considered. The employed material parameters are listed in Table 1. First, to check the correctness of the present RBF method, we select the inverse MQ functions to calculate the phononic band

w= 1

21

Σ ω wexact

w= 1

, (23)

where the number of k is 21. Fig. 3 plots the influences of the basis function number, the shape parameter and the node number on the relative error. From Fig. 3(a) and (b) (locally enlarged graph of Fig. 3(a)), we can clearly see that the shape parameter could affect the relative error significantly. The best value of cs which causes the relative error to be minimum is 4, 5 and 5 when the basis function number n is equal to 11, 15 and 19, respectively, which shows that the appropriate cs must be chosen for different n in order to reduce the relative error, and the shape parameter is not as bigger as better. According to previous researches [26–29], an optimal shape parameter can be found based on the golden search method and so on. Here, by using Fig. 3, we can quickly find out the optimal value of cs for a certain n. On the other hand, from

Z.-z. Yan et al. / Physica B 489 (2016) 1–11

7

Fig. 10. Band structures under different interface spring stiffness.

Fig. 3(a) and (b), we find that the increase of the basis function number leads to the decrease of the relative error value. In fact, it is verified again by Fig. 3(c) which indicates that the relative error indeed decreases with the increase of the number of basis functions although the relative errors are different by choosing different shape parameters cs . Also, it should be pointed out that the number of nodes is not the more the better, which may cause the matrix singular phenomenon along with the increase of the numbers. For the MQ and Gauss RBFs, we will only show the best value of cs for the certain n through the above method because the other conclusions are the same as those of the inverse MQ. For the MQ basis function, according to Fig. 4, when n is equal to 11, 15 and 19, the best value of cs corresponds to 4.8, 5 and 6.8, respectively. For the Gauss basis function, Fig. 5 indicates that when n is equal to 11, 15 and 19, the best value of cs corresponds to 0.005, 0.012 and

0.005, respectively, and the best value of cs is not unique. According to the above tests, we can easily find that in the RBF simulation of the PCs, a large node number is necessary if we want to get more correct results in the high frequencies. However, when the number of nodes is increased, the matrix singular phenomenon will appear. In this paper, we have ascertained by trial that 19 nodes are sufficient to obtain the correct results. This issue needs to be investigated in depth. On the other hand, the shape parameter affections are obtained in the PCs, according to which an optimal shape parameter can be found. Of course, we can also find an optimal shape parameter based on other methods [28,29]. So, a pre-calculation will be done to find out the shape parameter cs that has the smallest relative errors, and the reasonable node number will be given in our numerical examples in the following. Fig. 6 shows the results obtained by using the above three kinds

8

Z.-z. Yan et al. / Physica B 489 (2016) 1–11

Fig. 11. Band gap widths under different interface spring coefficients and length ratios.

Fig. 13. Comparison of the band structures between the three-phase model (h ¼ 0.03) and the perfect interface.

freedom, the results of the RBF method are more accurate than those of the FEM. Besides, in the log10(⋅) scale the RBF collocation method converges faster than the FEM in terms of the slope of the linearly fitted error curve. The computation time of the RBF method and the FEM for the 12 bands in Fig. 2 is compared in Fig. 8. The calculation is implemented on the sample laptop with DELL-Inspiron-14-7437. The computation time of the RBF method is less than that of the FEM in the case of the same nodes or degrees of freedom. 3.2. Effects of the imperfect interfaces on the band structures

Fig. 12. Band gap widths under different interface spring coefficients and density ratios.

of RBFs with their respectively best cs values for 19 basis functions. It can be clearly seen that in the low frequency range, the inverse MQ can lead to the most ideal results among the three types of RBFs, although they are all somewhat inconsistent with the exact analytical solutions in the high frequencies. So in the following numerical calculations, we employ the inverse MQ and n ¼ 19 in the RBF method. It is necessary to point out that the number of basis functions must be increased if the higher frequency band structures are calculated, and the optimal shape parameters need to be re-determined at this time. Next, the computational efficiency of the RBF collocation method is analyzed and discussed. In Fig. 7, the convergence rate of the present RBF method for the average of all eigenvalues of the 8th band in Fig. 2 is given and it is compared with that of the FEM to illustrate the superiority of the RBF method for the investigation of the frequency band gap properties. For comparison's sake, the relative error is used and defined as (Sυ − Sς ) /Sυ , where Sυ is the result of the RBF method by using 101 nodes, and Sς is the result of the RBF method or the FEM with the actual degrees of freedom. From Fig. 7, we can find that the order of the accuracy can reach 10−3 with about 20 nodes, and with the same nodes or degrees of

We first consider the effect of the spring interface on the phononic band structures. Fig. 9 plots the results with the perfect interface and the spring-interface model. The spring stiffness is chosen as F = 5F0 ( F0 = 1.5014 × 10−11). From Fig. 9, we can observe that compared with the perfect interface, the frequency bands become more flat slightly in the case of the spring-interface model, which shows that the interface imperfection has a significant effect on the phononic band structures. Next, to account for the effect of the interface imperfection, we give the detailed analysis of the band structures under different interface spring stiffness. From Fig. 10(a) and (b), we can find that the band structures are nearly unchanged when the interface spring stiffness F is very small, which indicates that a slight interface imperfection has little effect on the wave propagation characteristics in the PNCs. However, the band structures become more and more downward and flat with the increase of F as shown in Fig. 10(c)–(e). For example, with the increase of F , more flat bands appear from F = F0 to F = 30F0 , namely, the second band becomes more flat and the upper edge of the first band becomes more downward with the increase of F which causes the low-frequency band gap width to be wider. These phenomena indicate that the imperfect interface has a significant impact on the wave behaviors. With the stiffness of the spring-interface increasing, the localized resonance flat bands appear, and the first band gap becomes wider. In the following, we consider the influences of different length ratios on the width of the first band gap. Fig. 11 shows the width of the first band gap under different length ratios and spring coefficients. We can see that the band gap width decreases with the increase of the length ratios no matter what the interface imperfection is. Moreover, the first band gap becomes wider with the

Z.-z. Yan et al. / Physica B 489 (2016) 1–11

Fig. 14. Band structures in the case of the three-phase model for (a) h = 0.001; (b) h = 0.005; (c) h = 0.02 ; (d) h = 0.04 ; (e) h = 0.089 .

Fig. 15. The width of the first band gap under different length ratios and h .

Fig. 16. The width of the first band gap under different density ratios and h .

9

10

Z.-z. Yan et al. / Physica B 489 (2016) 1–11

increasing of F . In addition, we also observe that the interface needs to be stronger in order to obtain the lower and wider frequency ban gap in the phononic crystals with larger length ratio. In the following, the width of the first band gap under different density ratios and spring parameters is shown in Fig. 12. It is indicated that the width of the first band gap increases with the increase of density ratios and spring coefficients. Interestingly, the band gap width significantly increases no matter whether the interface is strong or weak. However, when the density ratio reaches a certain high value, the increase of density ratios has little effect on the band gap width. In this part, we consider the three-phase model. Fig. 13 illustrates the comparison of the band structures between the threephase model and the perfect interface. We can clearly see that the overall bands present an upward trend. However, the width of the band gaps does not change obviously. These phenomena are not the same as those of the spring interface model. Corresponding to the spring-interface model, the effects of the three-phase interface parameter h on the band structures are also investigated. From Fig. 14, we can see that when h is small, the band structures have almost no changes. However, the top bands show an upward trend, with the increase of h, the phenomenon becomes more obvious, namely, all the bands are gradually upward as illustrated in Fig. 14(c)–(e). Unlike the spring-interface model, there are no multiple flat bands in the band structures of the three-phase model. Next, the width of the first band gap under different length ratios and h is shown in Fig. 15, similar to the spring-interface model, the width of the first band gap decreases with the increase of length ratios. However, at the initial stages of the length ratios, the band gap width increases with the increase of h, but when the length ratio reaches a certain high value, the phenomenon is just the opposite. It is different from the spring-interface model which shows that in order to obtain a wider low-frequency band gap in the phononic crystals with large length ratio, the smaller value of h is necessary. Fig. 16 shows the effect of the density ratios on the width of the first band gap, compared with Fig. 12, we can clearly see that h has little to do with the widening of the first band gap, the width of the first band gap increases with the increase of density ratios.

4. Conclusions In this paper a novel RBF collocation method is developed to calculate the band structures of the phononic crystals with different interface models. The spring-interface and the three-phase models are considered. The effectiveness and the stability of the present RBF method are discussed, and the effects of different interface models, the interface imperfection conditions, the length ratios and the density ratios on the band gaps are analyzed. The following important results have been obtained 1) The RBF collocation method is very effective in the calculation of the band structures. The stability of the method is influenced by the types of basis functions, the node numbers and the shape parameters. 2) The band structures are nearly unchanged when the interface imperfection is very small. However, when the interface is weak, the band structures changed, for the spring-interface model, with the spring stiffness increasing, the band structures become more and more downward and flat. And the localized resonant flat bands appear with the first band gap width increases. For the three-phase model, the top bands show an upward trend. With the increase of h, the phenomenon becomes more obvious, namely, all the bands are gradually

upward. Unlike the spring-interface model, there are no multiple flat bands in the band structures of the three-phase model. 3) For the spring-interface model, the band gap width decreases with the increase of length ratios no matter what the interface imperfection is. Moreover, the first band gap width increases with the increase of F . In addition, we also observe that the interface needs to be weaker in order to obtain the lower and wider frequency ban gaps in the phononic crystals with large length ratio. The width of the first band gap increases with the increase of density ratios and spring coefficients. Interestingly, the band gap width significantly increases no matter whether the interface is strong or weak. However, when the density ratio reaches a certain high value, the increase of the density ratios has little effect on the band gap width. 4) For the three-phase model, similar to the spring-interface model, the width of the first band gap decreases with the increase of the length ratios. However, at the initial stages of the length ratios, the band gap width increases with the increase of h, but when the length ratio reaches a certain high value, the phenomenon is just the opposite. It is different from the springinterface model which shows that in order to obtain a wider low-frequency band gap in the phononic crystals with large length ratio, the smaller value of h is necessary. We can clearly see that h has little to do with the widening of the band gaps, the width of the first band gap increases with the increase of the density ratios.

Acknowledgments Support by the National Natural Science Foundation of China, China (No. 11002026 and 11372039), Beijing Natural Science Foundation, China (No. 3133039), the Scientific Research Foundation for the Returned (No. 20121832001) are gratefully acknowledged.

Appendix. The matrices M and N in Eq. (22)

⎡ M1 0 ⎤ ⎢ ⎥ M = ⎢ 0 M2 ⎥ , where ⎢B B ⎥ ⎣ 1 2⎦

⎡ φ ′′(x2) φ ′′(x2) 2 ⎢ 1 φ2′′(x3) ⋅⋅⋅ ⎢ φ ′′(x ) φ ′′(x ) ⎣ 1 n n 2

M1 =

λ1 + 2μ1 ⎢ φ ′′(x ) 3 1 ρ1 ⎢

M2 =

λ2 + 2μ2 ⎢ φ ′′(x n + 3) n+1 ρ2 ⎢

⋅⋅⋅ φn + 1′′(x2) ⎤ ⎥ ⋅⋅⋅ φn + 1′′(x3) ⎥ , ⎥ ⋅⋅⋅ ⎥ ⋅⋅⋅ φn + 1′′(x n)⎦

⎡ φ ′′(x n + 2) φ ′′(x n + 2) n+2 ⎢ n+1 φn + 2′′(x n + 3) ⋅⋅⋅ ⎢ φ ′′(x ) φ ′′(x ) ⎣ n + 1 2n 2n n+2

⎡ φ (x n + 1) φ2(x n + 1) 1 ⎢ E E 1 ⎢ 1 φ′ (x ) φ′ (x ) E2 2 n + 1 ⎢ E2 1 n + 1 B1 = ⎢ ikΗ ikΗ e φ2(x1) ⎢ e φ1(x1) ⎢ ikΗ E1 ikΗ E1 ⎢⎣ e E2 φ′1(x1) e E2 φ′2 (x1)

⋅⋅⋅ φ2n + 1′′(x n + 2)⎤ ⎥ ⋅⋅⋅ φ2n + 1′′(x n + 3) ⎥ , ⎥ ⋅⋅⋅ ⎥ ⋅⋅⋅ φ2n + 1′′(x2n) ⎦

φn + 1(x n + 1) ⎤ ⎥ ⎥ ⋅⋅⋅ ⎥ , ⎥ ⋅⋅⋅ eikΗ φn + 1(x1) ⎥ ⎥ E ⋅⋅⋅ eikΗ E1 φ′n + 1(x1)⎥⎦ 2 4 ×(n + 1)

⋅⋅⋅

⎡ φ (x n + 1) φ (x n + 1) n+ 2 ⎢ n+ 1 ⎢ φ′n + 1(x n + 1) φ′n + 2 (x n + 1) B2 = − ⎢ ) φn + 2(x2n + 1) φ (x ⎢ n + 1 2n + 1 ⎢⎣ φ′n + 1(x2n + 1) φ′n + 2 (x2n + 1)

E1 φ′ (x ) E2 n + 1 n + 1

⋅⋅⋅ φ2n + 1(x n + 1) ⎤ ⎥ ⋅⋅⋅ φ′2n + 1(x n + 1) ⎥ , ⋅⋅⋅ φ2n + 1(x2n + 1) ⎥ ⎥ ⋅⋅⋅ φ′2n + 1(x2n + 1)⎥⎦ 4 ×(n + 1)

Z.-z. Yan et al. / Physica B 489 (2016) 1–11

⎡N1 0 ⎤ ⎢ ⎥ N = ⎢ 0 N2 ⎥, where ⎢⎣ 0 0 ⎥⎦

11

References

⎡ φ (x2) φ (x2) 1 ⎢ 1 λ1 + 2μ1 ⎢ φ (x ) φ (x ) 3 3 1 2 N1 = ⋅⋅⋅ ρ1 ⎢ ⎢ φ (x ) φ (x ) ⎣ 1 n 2 n

⋅⋅⋅ φn + 1(x2) ⎤ ⎥ ⋅⋅⋅ φn + 1(x3) ⎥ , ⎥ ⋅⋅⋅ ⎥ ⋅⋅⋅ φn + 1(x n)⎦

⎡ φ (x n + 2) φ (x n + 2) n+ 2 ⎢ n+ 1 λ2 + 2μ2 ⎢ φ (x ) φ (x ) n+ 1 n+ 3 n+ 2 n+ 3 N2 = ⎢ ⋅⋅⋅ ρ2 ⎢ φ (x ) φ (x ) ⎣ n + 1 2n n + 2 2n

⋅⋅⋅ φ2n + 1(x n + 2)⎤ ⎥ ⋅⋅⋅ φ2n + 1(x n + 3) ⎥ , ⎥ ⋅⋅⋅ ⎥ ⋅⋅⋅ φ2n + 1(x2n) ⎦

α1 = (α1, α2, ⋯, αn + 1)T , α2 = (αn + 1, αn + 2, ⋯, α2n + 1)T , 1

φi(x v ) = (cs2(x v − xi )2 + 1)− 2 ; 3

φi′(x v ) = − cs2(x v − xi )[cs2(x v − xi )2 + 1]− 2 , 3

φi′′(x v ) = − cs2[cs2(x v − xi )2 + 1]− 2 +

5 3 4 2 cs [cs (x v − xi )2 + 1]− 2 . 2

Appendix A. Supplementary material Supplementary data associated with this article can be found in the online version at http://dx.doi.org/10.1016/j.physb.2016.02.026.

[1] M.S. Kushwaha, P. Halevi, L. Dobrzynski, B. Djafari-Rouhani, Phys. Rev. Lett. 71 (1993) 2022. [2] J.P. Jones, J.S. Whittier, J. Appl. Mech. 34 (1967) 905. [3] Y.S. Wang, G.L. Yu, Z.M. Zhang, Adv. Mech. 30 (2000) 378. [4] H. Mahiou, A. Beakou, Compos. Part A: Appl. Sci. Manuf. 29 (1998) 1035. [5] X.Y. Zhu, S. Zhong, D.Q. Sun, A.K. Ye, F.W. Deng, Physica. B 450 (2014) 121. [6] W. Liu, J. Chen, Y. Liu, X. Su, Phys. Lett. A 376 (2012) 605. [7] N. Zhen, Y.S. Wang, C.Z. Zhang, Physica. E 54 (2013) 125. [8] M.M. Sigalas, J. Acoust. Soc. Am. 101 (1997) 1256. [9] X.Z. Zhou, Y.S. Wang, Ch. Zhang, J. Appl. Phys. 106 (2009) 014903. [10] M. Kafesaki, E.N. Economou, Phys. Rev. B 60 (1999) 11993. [11] Z. Liu, C.T. Chan, P. Sheng, A.L. Goertzen, J.H. Page, Phys. Rev. B 62 (2000) 2446. [12] Z.Z. Yan, Y.S. Wang, Phys. Rev. B 74 (2006) 224303. [13] Z.Z. Yan, Y.S. Wang, Ch Zhang, Phys. Rev. B 78 (2008) 094306. [14] Y. Tanaka, Y. Tomoyasu, S. Tamura, Phys. Rev. B 62 (2000) 7387. [15] A. Khelif, A. Choujaa, S. Benchabane, B. Djafari-Rouhani, V. Laude, Appl. Phys. Lett. 84 (2004) 4400. [16] X.H. Hu, C.T. Chan, J. Zi, Phys. Rev. E 71 (2005) 055601. [17] P.F. Hsieh, T.T. Wu, J.H. Sun, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 53 (2006) 148. [18] X.X. Su, Y.S. Wang, J. Syn. Crys. 38 (2009) 1375. [19] W. Axmann, P. Kuchment, J. Comput. Phys. 150 (1999) 468. [20] J.B. Li, Y.S. Wang, Ch. Zhang, IEEE. Int. Ultrason. Symp. (2008) 1468. [21] G.R. Liu, Y.T. Gu, Springer Publications, The Netherlands, 2011. [22] W. Chen, Z.J. Fu, C.S. Chen, Springer Verlag, 2014. [23] Y. Benveniste, J. Mech. Phys. Solids. 54 (2006) 708–734. [24] R.E. Camley, B. Djafari-Rouhani, L. Dobrzynski, A.A. Maradudin, Phys. Rev. B 27 (1983) 7318. [25] E.H. EI Boudouti, B. Djafari-Rouhani, A. Akjouj, L. Dobrzynski, Surf. Sci. Rep. 64 (2009) 471–594. [26] G.E. Fasshatier, J.G. Zhang, Numer. Algorithms 45 (2007) 345–368. [27] B. Fornberg, C. Piret, J. Comput. Phys. 227 (2008) 2758–2780. [28] C.H. Tsai, J. Kolibal, M. Li, Eng. Anal. Bound. Elem. 34 (2010) 738–746. [29] S. Hernandez, I. Flores, J.A. Vazquez, J. Appl. Res. Technol. 10 (2012) 388–397.