1 Introduction

There is a steady growing interest in both experimental and theoretical investigations on the propagation of elastic/acoustic waves in phononic crystals in the last decades because of their remarkable features^{[1, 2, 3, 4, 5, 6, 7, 8, 9]} . The most interesting and important aspect of these features is that such artificial composites can exhibit elastic/acoustic wave bandgaps, in which elastic and acoustic waves are all forbidden regardless of the polarization and propagation direction of the elastic/acoustic waves. This behavior leads to many interesting physical properties and potential applications. The existence of elastic/acoustic bandgaps is significant for better understanding the propagation properties of elastic and acoustic waves in composite media^{[10, 11]}, aswellas various innovative applications such as acoustic filters, ultrasonic silent blocks, and acoustic mirrors. Therefore, one key problem in this research field is to compute the band structures (i.e., dispersion curves) and search for the so-called complete phononic bandgaps. Up to now, several numerical methods have been proposed in literature. Roughly speaking, these methods may be divided into two main categories: eigenfunction expansion methods (meshless methods) and domain/boundary discretization methods (domain or boundary methods). The former category includes the plane-wave expansion (PWE) method^{[1, 2, 12, 13]}, the wavelet method^{[8, 9]}, the multiple scattering theory (MST) method^{[14, 15]}, the Dirichlte-to-Neumann map (DtN-map) method^{[16, 17, 18]}, etc. The latter class includes the finite element method (FEM)^{[19, 20]}, the lumpedmass method (LMM)^[21], the boundary element method (BEM)^{[22, 23]}, the finite difference time domain (FDTD) method^[24], the spectral element method^{[25, 26]}, etc. Each of these methods has its own advantages and disadvantages. For example, the widely used PWE method, the wavelet method, and the FDTD method cannot handle the fluid-solid interface conditions properly. The MST and DtN-map methods are generally convenient, and have much small system matrices for almost all cases including the large acoustic mismatch case. However, these two methods are only suitable for circular and spherical inclusions. The BEM involves singular integrals and may miss some roots (eigenfrequencies)^[22]. The FEM is convenient for complex phononic systems with arbitrarily shaped inclusions, but the boundary conditions on the fluid-solid interaction interface are weakened^[20], and generally need a large number of elements.

Recently, the present authors have proposed a new numerical strategy to calculate the phononic band structures based on the generalized multipole technique (GMT)^[27]. However, only scalar wave propagation in two-dimensional (2D) phononic crystals was investigated^[27]. In this paper, the GMT is extended to calculate the band structures of in-plane waves in 2D phononic crystals which are composed of arbitrarily shaped solid or fluid cylinders (including vacuum holes) embedded in a solid matrix medium. According to the basic idea of the GMT^{[28, 29, 30]}, the wave fields in the matrix and the inclusion in one representative unit-cell are expanded by using the fundamental solutions with multiple origins. This makes it possible to analyze the phononic crystals with non-circular inclusions. An extra monopole source is introduced which acts as the external excitation. By varying the frequency of the excitation, the eigenvalues can be localized as the extreme points of an appropriately chosen function. By sweeping the frequency range of interest and the boundary of the irreducible first Brillouin zone (FBZ), the band structures can be obtained. In addition, in order to make the extreme points of the response curve more evident, two regularizing procedures are developed^[30]. Depending on the use of different basis functions, the GMT can be divided into two categories: the multiple multipole method (MMP)^{[27, 28, 29]} and the multiple monopole (MMoP) method^{[27, 30]}. The present authors have demonstrated in Ref. ^[27] that the MMoP method can yield satisfactory results. Therefore, in this paper the MMoP method is also used.

The remainder of this paper is organized as follows. In Section 2, the basic vector wave equations and the boundary conditions in a unit-cell are described. A brief overview of the GMT and the derivation of the proposed method for the computation of bandgaps of the phononic crystals are given in Section 3. Then, some numerical examples are presented and discussed in Section 4 to demonstrate the accuracy and efficiency of the method. Finally, Section 5 concludes the paper. 2 Basic equations and boundary conditions

Consider a 2D phononic crystal which is composed of straight and infinite cylinders (inclusions) forming a square or triangular array with the lattice constant a in a matrix (see Fig. 1). The inclusions may be solid, fluid, or vacuum. The cross-section of the inclusions may be of an arbitrary shape. Γ₁-Γ₄ (Γ₁-Γ₆) represent the edges of the square (hexagonal) unit-cell, and Γ₀ denotes the interface between the matrix and the inclusion. The subscript indices 0 and 1 are used to distinguish the matrix (Ω₀) and the inclusion (Ω₁) materials. Time-harmonic mixed plane wave modes propagating in the xy-plane are considered.

The governing equations for time-harmonic waves in isotropic and linear elastic solids can be described by the following vector equations^[31]:

where x = (x, y) ∈ R², j = 0, 1, u_j(x) = (u_{x, j}(x), u_{y, j}(x)) is the displacement vector at the position x = (x, y), ω is the angular frequency, λ_j , μ_j , ρ_j are the Lamé constant, the shear modulus, and the mass density in Ω₀ (j=0) and Ω₁ (j=1), respectively, and is the Laplace operator. For convenience, we introduce the Lamé potentials φj and ψj such that Then, the stress-potential relations are given by In terms of these potentials, the wave equation (1) is reduced to the following scalar Helmholtz equation: where k_{l, j} = ω/v_{l, j} and t_{t, j} = ω/v_{t, j} are the wave numbers of the longitudinal and transverse waves with v_{l, j} = √ (λ_j + 2μ_j)/ρ_j and v_{t, j} = √ μ_j/ρ_j being the velocities of the longitudinal and transverse waves, respectively.

Due to the periodicity of the system under consideration, the computation can be restricted to a single unit-cell. On the interface Γ₀ between the matrix and the inclusion, the displacement and traction vectors are continuous for the solid/solid system. Thus, we have

where n(x) is the unit normal vector of Γ₀, and σ_j (x) is the stress tensor on Γ₀. In addition, the periodicity of the structure implies that the displacements and tractions should satisfy the Bloch theorem (see Appendix A) on the boundaries of the unit-cell as follows: where i = √ -1 , x_m = m₁e₁+m₂e₂ with m = (m₁, m₂) ∈ Z² and e₁ and e₂ being the primitive translation vectors, and k = (k_x, k_y) is a real Bloch wave vector belonging to the FBZ.

If the inclusion Ω₁ is an ideal fluid, the displacement field of a propagating wave can be described by a displacement potential η₁(x), i.e.,

The pressure p₁(x) is then determined by where η1(x) satisfies the Helmholtz equation which is similar to (6) in the frequency domain, i.e., where k₁ = ω/v₁ is the wave number with v₁ = √ λ₁/ρ₁ being the longitudinal wave velocity. In this case, the boundary conditions on the interface Γ₀ require the continuity of the normal displacement and the normal traction as well as the zero tangential traction. Therefore, we have In the case of a vacuum cavity, the traction vector vanishes on the cavity surface, i.e., 3 GMT and calculation of band structures

The GMT^{[32, 33]} is a versatile and reliable numerical method for wave scattering problems in composites with arbitrarily shaped inclusions. It relies on a simple idea, namely, the scattered fields in various domains of the considered composite structure are expressed as a superposition of the fields radiated by adequate sources located outside the corresponding domains. These sources have no physical existence (i.e., fictitious sources) and can be multipole sources or monopole sources. The weights of these sources are unknowns to be determined. To obtain these unknowns, over-determined inhomogeneous linear equations are generated and solved by the least-squares method, and the inhomogeneous term of the equations stems from the incident wave field. In this paper, we use the generalized point-matching technique (GPMT), in which more collocation points than necessary are applied, i.e., the number of the collocation points is larger than the number of the sources or unknowns. Then, we can get the optimal weights of the sources. Once the weights of all the sources are obtained, we can express the scattered wave field in any domain.

For obtaining the band structures of a phononic crystal, an eigenvalue problem without the incident wave field has to be generally solved. However, in the present method, it is difficult to obtain a solution from the over-determined homogeneous linear equations. As suggested in Ref. ^[27] for anti-plane waves, we can introduce an external fictitious excitation that is applied in such a way that all wave modes are thereby excited. This excitation is equivalent to the incident wave field. This concept mimics the measurement of the resonance frequencies, where one always needs an excitation of the resonator and some test points (or ports) at which the signals are measured. By sweeping the frequency of the excitation, the peaks of an appropriately chosen function in the test points can be readily assigned to the resonance frequencies of the different modes, yielding the eigenvalues we are looking for. Next, we give a brief description of the method. For more details, the readers can refer to Refs. ^[30]-^[33].

According to the above mentioned basic idea of the GMT, instead of the homogeneous wave equation (6), we consider the following inhomogeneous wave equation in the matrix domain Ω₀ of the unit-cell:

where ƒ_{l, 0}(x) and ƒ_{t, 0}(x) are the excitation terms induced by the source at the position x_ex = (x_ex, yex) with the wave numbers k_{l, 0} and k_{t, 0}, respectively. The equation in the inclusion domain Ω₁ remains unchanged as (6). If we note that φ_ex(x) and ψ_ex(x) are indeed particular solutions corresponding to ƒ_{l, 0}(x) and ƒ_{t, 0}(x), respectively, the displacement potentials in Ω₀ can be expressed according to the GMT as follows: The displacement potentials in Ω₁ can be expressed as follows: For the ideal fluid inclusions, the displacement potential η1(x) in Ω₁ can be expressed as In (19)-(23), the superscript “approx” denotes the approximation to the actual fields, a_{j, n}, b_{j, n}, and c_{j, n} are understood as the unknown weights of the sources S_{j, n} which are located outside the corresponding domains, F_{φj, n}(x) and F_{ψj, n}(x) are the wave fields radiated from the sources S_{j, n} and satisfy the governing equation (6), and N₀ and N₁ are the numbers of the sources for the wave fields in the matrix and inclusion, respectively. According to the Vekua theory^[29], F_{φj, n}(x), F_{ψj, n}(x), and F_{η1, n}(x) form the complete basis of the series (19)-(23). In the MMP method^{[29, 30, 31, 32, 33, 34, 35]}, F_{φj, n}(x), F_{ψj, n}(x), and F_{η1, n}(x) take the following identical forms: where F_{j, n}^s is the position vector of the nth source S_{j, n}, and H_m ⁽¹⁾ (·) is the Hankel function of the first kind and mth order. In practical calculations, the truncation m ∈ [−M, M] is taken in (24)-(26). In the MMoP method (also known as the method of fundamental solutions^[30] and the method of fictitious sources^[33]), F_{φj, n}(x), F_{φj, n}(x), and F_{η1, n}(x) take the following simple forms: In this paper, we use the MMoP method, i.e., (27)-(29). The particular solutions φ_ex(x) and ψ_ex(x) take the same forms as (27) and (28)^[29], i.e., To determine the unknown weights, the interface conditions ((7)-(8) for the solid/solid system, (14)-(16) for the fluid/solid system, and (17) for the vacuum/solid system) are imposed on a set of collocation points x_k(k = 1, 2, · · · , N_c1) on the interface Γ₀. In addition, the Bloch boundary conditions (9) and (10) are imposed on another set of collocation points x_k(k = 1, 2, ··· , N_c2) on the boundary of the unit-cell. Then, we can get the following matrix equation: where A_r is an over-determined rectangular matrix with the number of the rows larger than that of the columns, the vector a contains the unknowns ((a_{0, n}, b_{0, n}, a_{1, n}, b_{1, n}) for the solid/solid system, (a_{0, n}, b_{0, n}, c_{1, n}) for the fluid/solid system, and (a_{0, n}, b_{0, n}) for the vacuum/solid system), and the vector b stems from the external excitation source. The specific forms of (32) for the square and triangular lattices are listed in Appendix A. Then, we can get an optimal solution (a_opt⁰, b_opt⁰, ab_opt¹, b_opt¹)T for the solid/solid system, (a_opt⁰, b_opt⁰, cb_opt¹)T for the fluid/solid system, and (a_opt⁰, b_opt⁰)T for the vacuum/solid system from (32) based on the least-squares method for a given frequency ω and a Bloch wave vector (k_x, k_y).

Substituting the optimal solutions a_opt⁰and b_opt⁰into (19) and (20) and subsequently into (2), we can define a norm as

where x_t¹ = (x_t¹, y_t¹) denotes the test point randomly distributed in the matrix domain Ω₀, and Nt is the number of the test points. For a Bloch vector (k_x, k_y) given along the boundary of the irreducible FBZ, the maximum values of the searching function ƒ₁(ω, k_x, k_y) will be found when we sweep the frequency range of interest which does not necessarily start from zero. These maximum values correspond to the eigenvalues we are looking for. Sweeping the boundary of the irreducible FBZ, we obtain the band structures.

It is noted that residual errors are inevitable in the least-squares method. According to (3), we can get the residual error vector in the following form:

Then, we can define another searching function as In order to improve the calculation accuracy, we also introduce an alternative searching function as follows: where p ≥ 2 is the order of this searching function. The eigenvalues to be determined are characterized by the minima of the searching function ƒ₂ or ƒ₃. In this paper, we use the searching function ƒ₃ for the solid/solid system and the vacuum/solid system. For the fluid/solid system, the vibration is localized in the fluid inclusion generally, especially for the fluid with a low density. Thus, we substitute the optimal solution c_opt¹ first into (23) and then into (12), and define another norm as where x_t¹ = (x_t¹, y_t¹) denotes the test point randomly distributed in the matrix domain Ω₁. However, if the above method is used directly, the response curve of ƒ₃ (or ƒ₁, ƒ₂, ƒ₄)

as a function of the frequency ω may have no extreme point^{[27, 30]}. To overcome this difficulty, Reutskiy^[30] proposed two regularizing procedures which can mark out the extreme points clearly. The first one is to introduce a damping term in the governing equations (6) and (13). Then instead of (6) and (13), we have

where ε is a small parameter. In this case, the basis functions (27)-(29) should be modified as In the second regularizing procedure, a shift δ_k between the wave numbers of the external exciting source and the studied mode is introduced, and then we take the particular solutions as follows: These two regularizing techniques are named as the ε-procedure and δ-procedure, respectively^[30]. In general, accurate results can be obtained by using any one of the two procedures. However, the latter one is generally used because the calculations of the Hankel functions with a complex argument in the first procedure are computationally demanding. 4 Numerical examples and discussion

In this section, the above described method based on the GMT is used to calculate the band structures for the mixed modes of the in-plane waves in some typical 2D phononic crystals. In order to validate the present method, the results are compared with those of the FEM. 4.1 Solid/solid system

We first check the accuracy and efficiency of the present method by computing the band structures of the periodic systems composed of circular aurum (Au) cylinders and the epoxy matrix in a square or triangular lattice with the filling fraction ƒ = 0.40. The densities and the wave velocities of the component materials are ρ₀ = 1 180 kg/m³, v_{l, 0} = 2 535 m/s, and v_{t, 0} = 1 157 m/s for the epoxy, and ρ₁ = 19 500 kg/m³, v_{l, 1} = 3 360 m/s, and v_{t, 0} = 1 239 m/s for Au. The wave impedance ratios of the systems are ρ₁v_{l, 1}/(ρ₀v_{l, 0}) = 21.9 and ρ₁v_{t, 1}/(ρ₀v_{t, 0}) = 17.7. It is well known that some other methods, such as the PWE method^[34], may yield poor results for the periodic systems with large acoustic mismatch. In the present method, the distribution of the monopole sources for both the square and triangular lattices is shown in Fig. 2.

For the square lattice as depicted in Fig. 2(a), the radius of the inclusion is r = a√ f/π = 0.36a, the sources inside the inclusion, outside the inclusion, and outside the unit-cell are distributed along circles with the radii r₁ = 0.19a, r₂ = 0.52a, and r₃ = 0.75a, respectively. We choose the numbers of the sources within the inclusion, outside the inclusion, and outside the unit-cell as N₁ = 28 , N₂ = 32, and N₃ = 40, and the numbers of the collocation points on the interface and the boundary of the unit-cell are N_c1 = 36 and N_c2 = 60, respectively. The number of the test points is N_t = 2 unless otherwise stated. In fact this method is not very sensitive to the number of the test points^{[29, 30]}. The position of the external excitation source is set at x_ex = (0.2a, 0.1a), and the regularizing parameter in the δ-procedure is δ_ω = δka/(2π) = 1.0. For the triangular lattice, we choose r₁ = 0.18a, r₂ = 0.50a, r₃ = 0.70a, N₁ = 24, N₂ = 30, N₃ = 42, N_c1 = 36, N_c2 = 60, x_ex = (0.2a, 0.1a), and δ_ω = 1.0. The eigenvalues are obtained from the searching function ƒ₃.

Figure 3 illustrates the band structures calculated by the present method and the FEM. The vertical axis is the normalized frequency ωa/(2πv_{t, 0}), and the horizontal axis represents the edge of the irreducible FBZ. The results of the FEM are obtained by the acoustics module of COMSOL Multiphysics 3.5a. The numbers of the degrees of freedom are 2 960 for the square lattice and 2 082 for the triangular lattice, respectively. Here, we use the “Direct (SPOOLES)” solver to solve the linear system of algebraic equations. The constraint handling method we used is the “Lagrange multipliers” method. The corresponding computational parameters in the MMoP and FEM models are summarized in Table 1. We can see that the results from the two methods are in good agreement, and they are also in good agreement with those from the DtN-map method^[16].

From this example, it is demonstrated that the developed MMoP method can compute the band structures of in-plane waves in 2D solid/solid phononic crystals accurately and efficiently. It is noted that the convergence and the accuracy of the method are affected by many factors including the number and the position of the fictitious sources, the number of the collocation points, the regularizing parameter, and the position of the external excitation soured. Next, we discuss the effects of these factors on the convergence and the accuracy of the present method by taking the square lattice as an example. Of course, similar investigations can be done for the triangular lattice. 4.1.1 Effects of monopole sources and collocation points

We first discuss the influence of the monopole sources by computing the eigenfrequency ω_e of the fifth band at k = (π/(2a), 0) for the square lattice (see Fig. 3(a)). Figure 4(a) shows the convergence of the computation with the number of the monopole sources when the number of the collocation points is sufficiently large. Fast and good convergence is observed although some response curves exhibit slight oscillations with the amplitude less than 5 before convergence, as shown in Fig. 4(a1). From these figures, we can find that the minimum numbers of the monopole sources are taben as N₁ = 22, N₂ = 22, and N₃ = 30. Here, we fix the position of the sources, but there is still a certain degree of freedom for setting the position of the sources. Generally, the sources should be sufficiently far from the boundaries^[30]. Figure 4(b) shows the convergence of the results with the increase of the distance between the sources and the interface when the numbers of the sources are taken as N₁ = 32, N₂ = 32, and N₃ = 40, and the numbers of the collocation points are chosen as N_c1 = 80 and N_c2 = 120. From this figure, we can conclude that the results are convergent when the numbers of the sources and the collocation points are large enough. However, no convergent result can be obtained when the sources are too close to the expanded domain. For example, no convergence is observed when r₁ = 0.30a for the sources inside the inclusion (see Fig. 4(b1)). This is because that some numerical dependencies are induced when the distance between the adjacent sources is too small^[29]. Similar behaviors are shown in Fig. 4(b2) and 4(b3) for the sources outside the inclusion and outside the unit-cell.

In the previous examples, we assume that the number of the collocation points is large enough in order to guarantee the convergence. Moreover, the number of the collocation points should be kept to a minimum to reduce the computing time. In the following, we discuss the influence of the number of the collocation points when the number and the position of the fictitious sources are fixed. Figures. 4 4(c1) and 4(c2) show the convergence of the computation with the number of the collocation points for fixed sources. The collocation points are uniformly distributed on the interface and the boundary of the unit-cell. We can observe here that the results are convergent with the increasing number of the collocation points. The minimum numbers of the collocation points for the convergence of the computation are N_c1 = 32 and N_c2 = 44. In practical calculations, we can take N_c1 = 36 and N_c2 = 60 in order to guarantee the convergence.

4.1.2 Effects of regularizing parameter and external excitation

In the proposed method, the regularizing parameter δ_ω plays an important role in making the eigenvalues more prominent. Figure 5(a) presents the influence of this parameter on the response curves for the square lattices. From this figure, it can be seen that the minimum points (values) of the response curves become evident with the increasing regularizing parameter. In general, we take the regularizing parameter in the range of 0.1 ≤ δ_ω ≤ 2.5. Of course, there are other possibilities for regularizing the response curves, for example, double-precision calculations are used, and the number of the collocation points are increasing. However, they will increase the computational effort.

In order to search the eigenvalues, we introduce a fictitious external excitation which is applied in such a way that all wave modes are thereby excited. This concept mimics the measurement of the resonance frequency. This procedure was first introduced in the method of auxiliary sources^[36], and later extended to MMP by Moreno et al.^[34] and Smajic et al.^[35], and to the method of fundamental solutions by Reutskiy^[30]. Interestingly, the type and the location of the external excitation have several options. For convenience, we choose the excitation having the same form as the fictitious sources (see (30) and (31)). The effects of the location of the external excitation on the response curves and the convergence are shown in Fig. 5(b). From this figure, we can observe that the convergence is not sensitive to the location of the excitation. However, the excitation should not be located on the high-symmetry points of the unit-cell. Otherwise, some wave modes cannot be excited, leading to the loss of the roots or eigenfrequencies (see Figs. 5(b3) and 5(b6)).

As mentioned before, the present MMoP method is flexible to deal with arbitrarily shaped inclusions. To demonstrate this advantage, we conduct the computation for a square lattice of Au cylinders with a square (see Fig. 6(a)) or elliptical (see Fig. 6(b)) cross-section embedded in the epoxy matrix. For all cases under consideration, the filling ratio is chosen as ƒ = 0.40. In the case of square inclusions, the side-length of the square inclusion is r = a√ƒ = 0.632 5a. The sources inside the inclusion are distributed along a square with the side-length r₁ = 0.50a, the sources outside the inclusion are distributed along a circle with the radius r₂ = 0.50a, and the sources outside the unit-cell are distributed along a circle with the radius r₃ = 0.75a. The numbers of the sources are taken to be N₁ = 48, N₂ = 52, and N₃ = 60, and the numbers of the collocation points to be N_c1 = 80 and N_c2 = 100. The regularizing parameter is chosen as δ_ω = 1.0, and the position of the external excitation source is selected as x_ex = (0.2a, 0.04a), as shown in Fig. 6(a). In the system with elliptical inclusions, the ratio of the semi-major axis and the semi-minor axis of the inclusion is assumed as b : c = 0.8 : 1.0. The sources inside the elliptical inclusion are distributed along an ellipse with the semi-major and semi-minor axes being b₁ = 0.28a and N₁ = 0.36a, the sources outside the elliptical inclusion are distributed along an ellipse with the semi-major and semi-minor axes being b2 = 0.40a and c2 = 0.50a, and the sources outside the unit-cell are distributed along a circle with the radius r₃ = 0.75a.

We choose the numbers of the sources as N₁ = 44, N₂ = 52, and N₃ = 60, and the numbers of the collocation points as N_c1 = 68 and N_c2 = 96. The regularizing parameter is taken as δ_ω = 1.0, and the position of the external excitation source is x_ex = (0.2a, 0.1a) (see Fig. 6(b)). It is important to point out that the numbers of the sources and the collocation points are larger than those in the circular inclusion case. The numbers of the degrees of freedom are 2 810 for the square inclusion case and 2 834 for the elliptical inclusion case, resdentudy. The corresponding computational parameters in the MMoP and FEM models are summarized in Table 1. Figure 7 shows the band structures calculated by the MMoP method (scattered dots) and the FEM (solid lines). It can be seen here that the results obtained by the two methods are in good agreement. The investigated examples verify the applicability of the MMoP method for non-circular inclusions, even for the inclusions with sharp corners.

4.2 Vacuum/solid system

In this subsection, we analyze both the square and the triangular lattices of circular cylindrical vacuum cavities in an aluminum (Al) matrix with the filling fraction ƒ = 0.55. The densities and the wave velocities of the Al matrix are ρ₀ = 2 697 kg/m³, v_{l, 0} = 6 410m/s, and v_{t, 0} = 3 110 m/s, respectively. The MMoP models for this system are shown in Fig. 2 (but without the sources “T” in the inclusion). The interface condition (17) is used. Here, we choose r₁ = 0.20a, r₃ = 0.80a, N₁ = 40, N₃ = 56, N_c1 = 40, N_c2 = 96, δ_ω = 1.0, and x_ex = (0.2a, 0.1a) for the square lattice, and r₁ = 0.20a, r₃ = 0.75a, N₁ = 36, N₃ = 48, N_c1 = 42, N_c2 = 90, δ_ω = 1.0, and x_ex = (0.2a, 0.1a) for the triangular one. The numbers of the degrees of freedom are 1 400 for the square lattice and 1 032 for the triangular lattice, respectively. The solver of the linear system of algebraic equations and the constraint handling method we used in COMSOL Multiphysics for vacuum/solid system are the same as those for the previously investigated solid/solid systems. The corresponding computational parameters in the MMoP method and FEM models are summarized in Table 1. Figure 8 illustrates the band structures calculated by the MMoP method (scattered dots) and the FEM (solid lines). Here again, we can see that the results from the two methods are almost identical, and that these results are also in good agreement with those of the DtN-map method^[16] and the BEM^[23]. This validates the applicability of the present GMT for vacuum/solid systems.

4.3 Fluid/solid system

As we know, some numerical methods, such as the PWE method and the wavelet method, do not work well for the fluid/solid system^{[8, 37]}. However, the present MMoP method can calculate the band structures of this mixed system by considering the different wave modes in the fluid and solid as well as the proper fluid-solid interface conditions (14)-(16). In this subsection, we first consider the examples in Subsection 4.2 where the circular cavities are filled now with air. The density and the wave velocity of air are ρ₁ = 1.21 kg/m³ and v₁ = 340 m/s. This system has a very small acoustic impedance ratio ρ₁v₁/(ρ₀v_{l, 0}) = 0.000 023 8. The MMoP models for this system are shown in Fig. 2. Here, we choose r₁ = 0.20a, r₂ = 0.56a, r₃ = 0.80a, N₁ = 36, N₂ = 26, N₃ = 52, N_c1 = 40, N_c2 = 72, δ_ω = 1.0, and x_ex = (0.2a, 0.1a) for the square lattice, and r₁ = 0.20a, r₂ = 0.52a, r₃ = 0.75a, N₁ = 36, N₂ = 24, N₃ = 48, N_c1 = 40, N_c2 = 72, δ_ω = 1.0, and x_ex = (0.2a, 0.1a) for the triangular lattice. The numbers of the degrees of freedom are 2 173 for the square lattice and 1 569 for the triangular lattice, respectively. It is noted here that the “Direct (SPOOLES)” solver is used to solve the linear system of algebraic equations. And the used constraint handling method is again the “Lagrange multipliers” method for fluid/solid systems. The corresponding computational parameters in the MMoP and FEM models are summarized in Table 1, and the eigenvalue searching function ƒ₄ is used. Figure 9 shows the band structures calculated by the MMoP method (scattered dots) and the FEM (solid lines). From these figures, it can be recognized that the results obtained from the two methods are almost identical. Moreover, these results are also in good agreement with the results of Ref. ^[20] (see Fig. 4 therein). A distinguishing feature shown in these figures is the appearance of many flat bands. These flat bands, which are completely different from those appearing when the PWE method is used^[37] and did not appear in the predicted results by the DtN-map method^[16], are physical bands and associated with the localized vibration modes of the air inside the circular cavity^[20]. To better understand these localized vibration modes, we enlarge the curve of the searching function ƒ₄ in the normalized frequency range [0.05, 0.25] with the corresponding wave vector k = (π/a, π/a) for the square lattice in Fig. 10. From this figure, we can identify that the flat bands correspond to very small peaks (A, C, D, and E). This implies that these localized vibration modes do not affect the wave propagation behavior in the air/Al system.

As the second example, we consider a water/iron system. The material properties of the water are ρ₁ = 1 025 kg/m³ and v₁ = 1 531 m/s, and those of the iron are ρ₀ = 7 670 kg/m³, v_{l, 0} = 6 010 m/s, and v_{t, 0} = 3 230 m/s. The acoustic impedance ratio of the system is ρ₁v₁/(ρ₀v_{l, 0}) = 0.034 which is much larger than that of the air/Al system. The filling fraction of this system is ƒ = 0.50. We choose the parameters of the computational model in the MMoP method for this system as r₁ = 0.20a, r₂ = 0.56a, r₃ = 0.80a, N₁ = 36, N₂ = 26, N₃ = 52, N_c1 = 40, N_c2 = 72, δ_ω = 1.0, x_ex = (0.2a, 0.1a). Here, we choose all the same computational parameters as those in the air/Al system, even though the filling ratio is not exactly the same (see Table 1). Here and in the following examples, we use the searching function ƒ₃. The computed band structures are shown in Fig. 11, which are in good agreement with those calculated by the MST method^[15], the DtN-map method^[16], the FEM^[20], and the BEM^[23]. It can be seen from Fig. 11 that there is no complete band-gap, and that the number of the flat bands is drastically reduced. This is because that the water has a higher density than the air, and its interaction with the solid is much stronger than that of the air. This example verifies the applicability of the MMoP method for fluid/solid systems with a large acoustic impedance ratio. It also shows that all computational parameters in the proposed method are independent of the material properties, but only depend on the geometry of the periodic structures under consideration. This makes it more convenient to deal with periodic systems having various acoustic impedance ratios. In contrast, the number of the degree of freedom in the FEM is increased drastically (see Table 1), which means that the computational cost is higher.

It is noted that the numbers of the sources and collocation points needed in the vacuum/solid or the fluid/solid system are larger than those used in the solid/solid system.

Next, we consider the system of circular mercury (Hg) cylinders in the Al matrix in a square lattice. The density and the wave velocity of Hg are ρ₁ = 13 600 kg/m³ and v₁ = 1 450 m/s, respectively. This system has a much larger acoustic impedance ratio ρ₁v₁/(ρ₀v_{l, 0}) = 1.14. The filling fraction is ƒ = 0.40. Similar to the water/iron system, we choose

The corresponding computational parameters in the MMoP and FEM models are summarized in Table 1. The results of the MMoP method and FEM are shown in Fig. 12. They are also in good agreement with those calculated by the DtN-map method^[16], the FEM^[20], and the BEM^[23]. No apparent flat bands can be observed here, because Hg is a special fluid with a very high density at the room temperature, and its interaction with the solid is very strong. From Table 1, we can find that the number of the degree of freedom in the FEM for the Hg/Al system is larger than that for the air/Al and water/iron systems.

Finally, we consider the Hg/rubber system which has a very large acoustic impedance ratio ρ₁v₁/(ρ₀v_{l, 0}) = 10.776. The material properties of the rubber are ρ₀ = 1 000 kg/m³, v_{l, 0} = 1 830 m/s, and v_{t, 0} = 500 m/s, and the filling fraction is ƒ = 0.40. The computational parameters are chosen as

The band structures computed by the MMoP method and FEM are presented in Fig. 13. From this figure, it can be observed that the results from the two methods are in good agreement, and that no physical flat band appears in this system. The corresponding computational parameters in the MMoP and FEM models are summarized in Table 1. It is noted that the computational parameters of the MMoP model for this system are the same as those for the Hg/Al system with the square lattice. However, the FEM needs a larger number of the degree of freedom for the Hg/rubber system, even though this system has the same geometrical structure as the Hg/Al system. This fact shows again that the MMoP method is more convenient to deal with periodic systems with various acoustic impedance ratios. 5 Conclusions

In this paper, we present an accurate and efficient numerical method for the band structure calculations of in-plane waves in 2D phononic crystals based on the GMT. The method expands the wave fields using the fundamental solutions with multiple sources, which makes it possible to deal with periodic composite systems containing non-circular inclusions. An additional excitation source is introduced, and thus the original eigenvalue problem is transformed into a response problem. The extreme values of the response curves correspond to the eigenvalues which we search for. Numerical calculations are performed for typical phononic crystal systems with various acoustic impedance ratios. The following conclusions can be drawn from the present numerical results.

(I) The MMoP method based on the GMT can yield accurate results for the band structures of solid/solid, vacuum/solid, and fluid/solid phononic crystals with various acoustic impedance ratios, including the localized vibration modes. They exhibit some advantages such as good stability, fast convergence, and direct computation in an arbitrary frequency range. In this method, the fluid-solid interface conditions are taken into account strictly for periodic systems with fluid inclusions.

(II) All parameters of the computational model in the MMoP method are independent of the material properties for the same type of periodic systems (solid/solid, fluid/solid, or vacuum/solid), but depend on the geometry of the periodic structure under consideration. This fact makes the method more convenient to deal with phononic crystal systems with various acoustic impedance ratios.

(III) An important feature of the MMoP method is that it can handle periodic systems with arbitrarily shaped inclusions, which is difficult for other methods based on the expansions of cylindrical wave functions, such as the MST and DtN-map methods. Appendix A

Application of the Bloch theorem (see (9) and (10)) to the boundary of a unit-cell of the square lattice leads to the following equations:

Combining (7) and (8) with (A1)-(A6) and taking collocation points N_c1 and N_c2 on the interface Γ₀ and the external boundary of the unit-cell, respectively, we have

where U_{Γlk, j, l, x}, U_{Γlk, j, l, y}, U_{Γlk, j, t, x}, U_{Γlk, j, t, y}, σ_{Γlk, j, l, x}, σ_{Γlk, j, l, y}, σ_{Γlk, j, t, x}, and σ_{Γlk, j, t, y} stem from the displacements and tractions on the corresponding boundaries Γij in the x- and y-directions, respectively. The subscript "ex” denotes the external excitation. (A7) is the specific form of (32) for the solid/solid system in a square lattice.

For the triangular lattice, the application of the Bloch theorem to the boundary of a unit-cell leads to similar results which are given by

Click Browse formula

By combining (7) and (8) with (A8)-(A16) and taking collocation points N_c1 and N_c2 on the interface Γ₀ and the boundary of the unit-cell, respectively, the following matrix equation can be obtained:

Click Browse formula

which is the specific form of (32) for the solid/solid system in a triangular lattice.

Similarly, combining (14)-(16) with (A1)-(A6) and taking collocation points N_c1 and N_c2 on the interface Γ₀ and the external boundary of the unit-cell, respectively, we can obtain the specific form of (32) for the fluid/solid system in a square lattice as follows:

Click Browse formula

where the superscripts “n” and “t” denote the normal and tangential directions, respectively.

Similarly, for the vacuum/solid system in a square lattice, we can obtain the following equation by combining (17) with (A1)-(A6):

Click Browse formula

The corresponding equations for the triangular lattices of the fluid/solid and vacuum/solid systems can be obtained similarly but are not given here.