A Localized Collocation Solver Based on T-Complete Functions for Anti-Plane Transverse Elastic Wave Propagation Analysis in 2D Phononic Crystals

: In this paper, we introduce a novel localized collocation solver for two-dimensional (2D) phononic crystal analysis. In the proposed collocation solver, the displacement at each node is expressed as a linear combination of T-complete functions in each stencil support and the sparse linear system is obtained by satisfying the considered governing equation at interior nodes and boundary conditions at boundary nodes. As compared with ﬁnite element method (FEM) results and the analytical solutions, the efﬁciency and accuracy of the proposed localized collocation solver are veriﬁed under a benchmark example. Then, the proposed method is applied to 2D phononic crystals with various lattice forms and scatterer shapes, where the related band structures, transmission spectra, and displacement amplitude distributions are calculated as compared with the FEM.


Introduction
In recent decades, more and more attention has been paid to a new kind of artificial periodic composite structures, which are well known as phononic crystals [1][2][3][4][5][6]. Due to the periodicity of phononic crystals, a wave cannot propagate through the structures within certain frequency band ranges, which is called the phononic bandgaps. Thanks to this bandgap characteristic, one can control the direction and path of wave propagation by designing a specific phononic crystal [7][8][9][10], and also can realize the shielding and reflection of a specific wave [11][12][13]. Since different shapes and placements of scatterers have a critical influence on wave propagation in phononic crystals, numerical simulation plays an important role in the computational design of phononic crystals.
In this study, we make the first attempt to apply the novel localized collocation Trefftz scheme to phononic crystal analysis, which includes a numerical solution of the Helmholtz eigenvalue equation for band structures and a numerical solution of Helmholtz equation for frequency responses and displacement amplitude distributions. This paper is organized as follows: The mathematical formulation of the anti-plane transverse elastic wave in the 2D phononic crystal is described in Section 2; the corresponding discretization formulation based on the localized collocation Trefftz method (LCTM) is clearly introduced in Section 3; in Section 4, as compared with the COMSOL simulation, the efficiency and accuracy of the proposed LCTM are verified under a benchmark example, and then the wave propagation behavior is investigated by calculation of band structures, transmission spectra, and displacement amplitude distributions in phononic crystals with various lattice forms and scatterer shapes; and in Section 5, our conclusions are summarized.

Mathematical Model of Shear Horizontal (SH) Wave Propagation by Phononic Crystals
Under the assumptions of the linear elastic wave theory [3], the propagation of timeharmonic SH waves in phononic crystals is under consideration. The corresponding governing equation [32] can be represented as follows: where ω stands for the angular frequency; µ and ρ denotes, respectively, the shear modulus and mass density of the solid; and u is the anti-plane displacement along the z-axis. The superscript j stands for the quantities related to the matrix (j = 0) or the scatterer (j = 1). In this study, two types of the phononic crystals, infinite periodic structure and semi-infinite periodic structure are investigated. For the perfectly phononic crystal with infinite periodic structure, as shown in Figure 1, we can only consider the unit cell instead of the whole structure. Γ 1 -Γ 4 (Γ 1 -Γ 6 ) denote the boundaries of the square (hexagonal) unit cell and Γ 0 represents the interface between the matrix and the scatterer. All the quantities of the anti-plane transverse elastic wave field on periodic boundaries Γ 1 -Γ 4 (Γ 1 -Γ 6 ) should satisfy the Bloch theorem, namely, where k = (k x 1 , k x 2 ) is the anti-plane Bloch wave vector; T = µ i ∂u i /∂n denotes the traction vector; n represents the unit normal vector; i = √ −1, a = m 1 a 1 + m 2 a 2 , in which m 1 , m 2 are arbitrary integers; and a 1 and a 2 stand for the fundamental translation vectors of the lattice. The following continuity conditions should be satisfied on the interface boundary: Γ 0  Therefore, the solution of the SH wave equation can be written as a superposition of anti-plane Bloch waves. The dispersion relation of all anti-plane Bloch waves forms the band structure which can be obtained by solving the following generalized eigenvalue problem, possibly subjected with surface or interface boundary conditions (2) and (3):  Therefore, the solution of the SH wave equation can be written as a superposition of anti-plane Bloch waves. The dispersion relation of all anti-plane Bloch waves forms the band structure which can be obtained by solving the following generalized eigenvalue problem, possibly subjected with surface or interface boundary conditions (2) and (3): where κ j = ω ρ j /µ j . The band structure can be obtained by evaluating the dispersion relations along the boundary of the irreducible Brillouin zone, in which λ = e ik x1 a , e ik x2 a = 1 from Γ − X; λ = e ik x2 a , e ik x1 a = −1 from X − M; λ = e ik x1 a = e ik x2 a from M − Γ in square lattice and λ = e ik x2 a , e ik x1 a = 1 from Γ − X; λ = e ik x2 a = e and λ = e ik x2 a = e i √ 3k x1 a from M − Γ in triangular lattice. Then, the corresponding eigen frequency, λ, can be obtained by solving the generalized eigenvalue equations with sweeping the angular frequency, ω, along the boundary of the irreducible Brillouin zone. For a perfectly phononic crystal with semi-infinite periodic structure in the y-direction, the periodic unit cell, as shown in Figure 2, can be considered, where Figure 2a,b represents two types of the SH wave propagation along ΓX− and ΓM− directions. Then, the displacements and tractions on the boundaries Γ * 2 , Γ * 4 satisfy the following relationships: and the left boundary Γ * 1 is excited by the following unit displacement: and the right boundary Γ * 3 is imposed on the following absorbing boundary condition: where wavenumber By solving Equation (9), the corresponding transmitted displacements on the right boundary

Localized Collocation Trefftz Method
In this section, the localized collocation Trefftz method (LCTM) is introduced to discretize the aforementioned generalized eigenvalue Equations (4) or (5) and the set of linear equations in Equation (9). By considering the interface continuity conditions (3) and boundary conditions (6)-(8), the following set of linear equations can be obtained: where wavenumber κ 0 = ω ρ 0 /µ 0 , κ 1 = ω ρ 1 /µ 1 . By solving Equation (9), the corresponding transmitted displacements on the right boundary u 0 (x Γ * 3 ) can be obtained for every given angular frequency ω. Then, the transmission spectra can be obtained by for different angular frequency ω.

Localized Collocation Trefftz Method
In this section, the localized collocation Trefftz method (LCTM) is introduced to discretize the aforementioned generalized eigenvalue Equations (4) or (5) and the set of linear equations in Equation (9).
To construct the LCTM numerical differentiation formulation, each given i-th node x i 0 with m nearest nodes x i 1 , x i 2 , . . . , x i m around x i 0 form a subdomain Ξ i , whose center can be chosen as Figure 3. Then, the approximated displacement u x i j inside the stencil support Ξ i can be represented by a linear combination of T-complete functions φ i k with unknown coefficients α i where the following T-complete functions are used to satisfy the governing Equation (1) in advance: where r i k , θ i k is constructed based on the polar coordinate system determined by the nodes in a subdomain Ξ i with the center x i as the origin and J i is Bessel function of the first kind of order i. Then, the vector of unknown coefficients α i can be expressed as follows:

Numerical Results
In this section, the efficiency and accuracy of the proposed LCTM are first verified under a benchmark example, and then the wave propagation behavior is numerically investigated by the calculation of band structures, transmission spectra, and displacement amplitude distributions in the phononic crystals with various lattice forms and scatterer shapes. Unless otherwise specified, aurum (Au) and epoxy are used as the materials of scatterers and matrix in phononic crystals, and the wave velocities and the densities of aurum and epoxy are, respectively, given by 0 1239 c ms = ,

Convergence and Numerical Efficiency Analysis
To verify the numerical efficiency and accuracy, first, the proposed LCTM is implemented to solve a benchmark example about SH wave scattering problem of a rigid circular cavity in an infinite elastic domain. Consider the incident SH wave as the form of 0 cos = ikr I e θ φ , the analytical solutions of the scattered field s u [59] for rigid boundary is as follows: where the k is wave number, the prime denotes differentiation with respect to ka, a stands for the radii of the circular scatterer, and the analytical solutions are numerically calculated by using the first 100 terms in the series solution Expression (16).
In the proposed LCTM implementation, the computational domain, as shown in Figure 4, needs to be truncated into finite domain with absorbing boundary ∞ Γ , and the weakened Sommerfeld boundary condition is used as the following absorbing boundary conditions (ABC) [59]: According to the conclusion from the literature [59] coupled with the extensive numerical investigation, the length of a reliable region RB is chosen as the half length of a whole truncated region TB. In this example, the region (a circle with radius TB = 10a) It should be mentioned, because the used T-complete function satisfies the governing equation in advance, the LCTM numerical differentiation formulation of u i 0 at x i 0 inside the domain of matrix and scatterer can be rewritten as follows: For a given angular frequency ω, substituting Equations (13)-(15) into the generalized eigenvalue Equations (4) or (5), the related system of linear algebra equations can be obtained to calculate the band structure, and substituting Equations (13)-(15) into the set of linear equations in Equation (9), the related system of linear algebra equations can be obtained to calculate the transmission spectra and displacement amplitude distributions.

Numerical Results
In this section, the efficiency and accuracy of the proposed LCTM are first verified under a benchmark example, and then the wave propagation behavior is numerically investigated by the calculation of band structures, transmission spectra, and displacement amplitude distributions in the phononic crystals with various lattice forms and scatterer shapes. Unless otherwise specified, aurum (Au) and epoxy are used as the materials of scatterers and matrix in phononic crystals, and the wave velocities and the densities of aurum and epoxy are, respectively, given by c 0 = 1239 m/s , ρ 0 = 19, 500 kg/m 3 , c 1 = 1161 m/s, ρ 1 = 1180 kg/m 3 ; in the simulation of transmission spectra and displacement amplitude distributions, the numbers of the unit cells are chosen as Na = 9, Nb = 4 and Nc = 7, the node number in each subdomain is m = 21.

Convergence and Numerical Efficiency Analysis
To verify the numerical efficiency and accuracy, first, the proposed LCTM is implemented to solve a benchmark example about SH wave scattering problem of a rigid circular cavity in an infinite elastic domain. Consider the incident SH wave as the form of φ I = e ikr cos θ 0 , the analytical solutions of the scattered field u s [59] for rigid boundary is as follows: where the k is wave number, the prime denotes differentiation with respect to ka, a stands for the radii of the circular scatterer, and the analytical solutions are numerically calculated by using the first 100 terms in the series solution Expression (16).
In the proposed LCTM implementation, the computational domain, as shown in Figure 4, needs to be truncated into finite domain with absorbing boundary Γ ∞ , and the weakened Sommerfeld boundary condition is used as the following absorbing boundary conditions (ABC) [59]:   In Table 5, the condition numbers of i Φ (Cf) and * A (Ca) in the proposed LCTM + ABC are given, and the CPU times of LCTM and FEM (COMSOL) without matrix generated are shown. Due to the highly ill-conditioning matrix i Φ , the Moore-Penrose inverse According to the conclusion from the literature [59] coupled with the extensive numerical investigation, the length of a reliable region RB is chosen as the half length of a whole truncated region TB. In this example, the region (a circle with radius TB = 10a) bounded by the absorbing boundary Γ ∞ can be divided into a reliable region (a circle with radius RB = 5a) and an error region. The numerical accuracy inside the reliable region is calculated by the L 2 relative error, which is defined as follow: where N is the number of total discretization nodes and u(x i ), u(x i ) denotes the numerical solution and analytical solution at x i . Tables 1 and 2 list the numerical errors (L 2 relative errors) in the reliable region obtained by the proposed LCTM and the FEM (COMSOL) coupled with the absorbing boundary conditions (17), where the symbol "/" denotes that the L 2 relative error is larger than 0.01 and ∆h stands for the average node spacing. It can be found from Tables 1 and 2 that both the proposed LCTM and the FEM provide acceptable results, and with an increasing wavenumber k, more refined nodes are required to achieve acceptable results. It should be mentioned, because of that the effect of the ABC (17), the numerical accuracy cannot be further improved with an increasing number of nodes/elements in the proposed LCTM/FEM. Therefore, more efficient additional techniques, i.e., the absorbing boundary conditions based on wave expansion functions (ABCWEF) and perfectly matched layer (PML), can be introduced. In this paper, to improve the numerical accuracy, the ABCWEF is employed in the proposed LCTM and the PML is employed in the FEM (COMSOL). One can find more details about the PML in the literature [59]. Here the ABCWEF is briefly introduced. By using the Hankel function, the numerical solution u x i j can be represented by a linear combination of Hankel functions as: in which γ 0 , γ 1 , · · · , γ 2n ABC −1 , γ 2n ABC are underdetermined coefficients and R i j , ϑ i j is constructed based on the polar coordinate system. To construct the ABCWEF formulation for each given i-th node x i 0 on the absorbing boundary Γ ∞ , its m nearest nodes x i 1 , x i 2 , . . . , x i m in the subdomain Ξ i are used, whose center is the center of the computational domain. By using the wave expansion formulation (19) at x i 0 on the absorbing boundary Γ ∞ and its m nearest nodes x i 1 , x i 2 , . . . , x i m in the subdomain Ξ i , the following system of linear equations can be obtained: and then the approximate displacement u x i 0 can be written as follow: substituting Equation (21) into Equation (13), the following formulation is obtained:   Table 3 lists the numerical errors (L 2 relative errors) in the reliable region obtained by the proposed LCTM coupled with the absorbing boundary conditions based on wave expansion functions (ABCWEF). It can be found that the proposed LCTM + ABCWEF can converge rapidly to the analytical solution with an increase in the node number, and it provides more accurate results than the results obtained by the LCTM + ABC. In addition, it is easy to introduce the PML to the FEM in the COMSOL software. As listed in Table 4, the FEM + PML provides more accurate results than the results obtained by the FEM + ABC and converge rapidly to the analytical solution with an increase in the element number. From Tables 3 and 4, as compared with the FEM + PML, the proposed LCTM + ABCWEF can provide more accurate results with less computational resources.  In Table 5, the condition numbers of Φ i (Cf) and A * (Ca) in the proposed LCTM + ABC are given, and the CPU times of LCTM and FEM (COMSOL) without matrix generated are shown. Due to the highly ill-conditioning matrix Φ i , the Moore-Penrose inverse technique is implemented to calculate Φ i −1 . We can observe that the condition numbers of the resultant matrix A * is relatively small. Additionally, it should be mentioned that the CPU times are in the same magnitude between the proposed LCTM and COMSOL, and the code of the proposed LCTM may be further improved to save the CPU time.

Phononic Crystals with Square Lattice
Next, we investigate the wave propagation behavior in the phononic crystals with square lattice and consider both the infinite periodic structure and semi-infinite periodic structure in the y-direction. For the infinite periodic structure, the related band structure is calculated. For the semi-infinite periodic structure in the y-direction, the related transmission spectra and displacement amplitude distributions are computed. All the numerical results obtained by the proposed LCTM + ABC are compared with the COMSOL simulation with the FEM + ABC and FEM + PML.

Circular Scatterer Case
In the first example, the unit cell structure of aurum/epoxy phononic crystal consists of a square lattice and a circular scatterer. The filling fraction of the circular scatterer is 0.4. Figure 5 shows the corresponding node distribution in the proposed LCTM. In the present numerical implementation, the discretization node number is N = 209 with node spacing ∆h = 1/15. The calculated band structures obtained by using the proposed LCTM and the FEM (COMSOL, 288 elements) are in good agreement, which is shown in Figure 6. In the first example, the unit cell structure of aurum/epoxy phononic crystal consists of a square lattice and a circular scatterer. The filling fraction of the circular scatterer is 0.4. Figure 5 shows the corresponding node distribution in the proposed LCTM. In the present numerical implementation, the discretization node number is The calculated band structures obtained by using the proposed LCTM and the FEM (COMSOL, 288 elements) are in good agreement, which is shown in Figure  6.    Then, the phononic crystal with semi-infinite periodic structure in the y-direction is investigated. The related node distribution is shown in Figure 7.  Then, the phononic crystal with semi-infinite periodic structure in the y-direction is investigated. The related node distribution is shown in Figure 7. Figure 8 presents the transmission spectra of anti-plane transverse elastic wave along with ΓXdirection. The band structures of the corresponding infinite periodic structure calculated by the proposed LCTM are also displayed in this figure. From Figure 8, it can be observed that the proposed LCTM with ∆h = 1/13 (N = 2729) can perform equally well with the FEM (COMSOL, 6507 elements), which reveals that the developed LCTM can also effectively calculate the SH wave transmission spectra of the phononic crystal. In addition, under the same materials and fill fraction, both the infinite periodic structure and semi-infinite periodic structure in the y-direction may hinder the SH wave propagation in the same angular frequency regions, namely, the bandgap of infinite periodic structure corresponds to the very low transmission coefficient of semi-infinite periodic structure. Therefore, in the design of phononic crystals, the use a nine-column finite period structure is recommended instead of an ideal infinite period structure to block the SH wave propagation at specific angular frequencies.
In order to validate and intuitively understand the band gap characteristics of phononic crystals, the displacement amplitude distributions computed by the proposed LCTM and FEM (ABC), FEM (PML) are presented at two specific normalized frequencies ( ωa 2πc 0 = 0.5, 1). All the numerical results obtained by the proposed LCTM and FEM (ABC), FEM (PML) are in good agreement, which reveals that the LCTM and FEM with the ABC are enough to obtain an acceptable solution for transmission spectra and displacement amplitude distributions. It is easy to observe from Figure 9 that the SH wave is blocked at the normalized frequency ( ωa 2πc 0 = 0.5) in the bandgap region, and it cannot hinder the SH wave propagation at the normalized frequency ( ωa 2πc 0 = 1) in the pass band region. calculate the SH wave transmission spectra of the phononic crystal. In addition, under the same materials and fill fraction, both the infinite periodic structure and semi-infinite periodic structure in the y-direction may hinder the SH wave propagation in the same angular frequency regions, namely, the bandgap of infinite periodic structure corresponds to the very low transmission coefficient of semi-infinite periodic structure. Therefore, in the design of phononic crystals, the use a nine-column finite period structure is recommended instead of an ideal infinite period structure to block the SH wave propagation at specific angular frequencies. FEM (PML) are in good agreement, which reveals that the LCTM and FEM with the ABC are enough to obtain an acceptable solution for transmission spectra and displacement amplitude distributions. It is easy to observe from Figure 9 that the SH wave is blocked at the normalized frequency (   To investigate the effect of the number of columns (scatterers) on the wave shielding, 7-column, 9-column and 11-column finite periodic structures are considered, as shown in Figure 10. Figure 11 shows the result of the transmission spectra of 7-column, 9-column, and 11-column finite periodic structures. As shown in Figure 11, with an increasing column number, the transmission coefficient tends to the one with infinite periodic structures. Since the transmission coefficients are almost the same between 9-column and 11column finite periodic structures, nine-column finite periodic structures are employed in the following numerical investigations.  To investigate the effect of the number of columns (scatterers) on the wave shielding, 7-column, 9-column and 11-column finite periodic structures are considered, as shown in Figure 10. Figure 11 shows the result of the transmission spectra of 7-column, 9-column, and 11-column finite periodic structures. As shown in Figure 11, with an increasing column number, the transmission coefficient tends to the one with infinite periodic structures. Since the transmission coefficients are almost the same between 9-column and 11-column finite periodic structures, nine-column finite periodic structures are employed in the following numerical investigations.
To investigate the effect of the number of columns (scatterers) on the wave shielding, 7-column, 9-column and 11-column finite periodic structures are considered, as shown in Figure 10. Figure 11 shows the result of the transmission spectra of 7-column, 9-column, and 11-column finite periodic structures. As shown in Figure 11, with an increasing column number, the transmission coefficient tends to the one with infinite periodic structures. Since the transmission coefficients are almost the same between 9-column and 11column finite periodic structures, nine-column finite periodic structures are employed in the following numerical investigations.

Cross-Shaped Scatterer Case
Here, the unit cell structure of aurum/epoxy phononic crystal consists of a square lattice and a cross-shaped scatterer and the measurements of the cross are 0.75 by 0.25. The filling fraction of the cross-shaped scatterer is 0.3125. Figures 12 and 13 show the corresponding node distributions in the proposed LCTM. In the present numerical implementation, the discretization node number is 238 N = with node spacing 1/13 h Δ = for band structure calculation, and the discretization node number is 3950 N = with the same node spacing 1/15 h Δ = for the calculation of transmission spectra and displacement amplitude distributions. The calculated band structures and transmission spectra obtained by using the proposed LCTM and the FEM (COMSOL, 264 elements for band structures, 6284 elements for transmission spectra) are in good agreement, which is shown in Figure 14.

Cross-Shaped Scatterer Case
Here, the unit cell structure of aurum/epoxy phononic crystal consists of a square lattice and a cross-shaped scatterer and the measurements of the cross are 0.75 by 0.25. The filling fraction of the cross-shaped scatterer is 0.3125. Figures 12 and 13 show the corresponding node distributions in the proposed LCTM. In the present numerical implementation, the discretization node number is N = 238 with node spacing ∆h = 1/13 for band structure calculation, and the discretization node number is N = 3950 with the same node spacing ∆h = 1/15 for the calculation of transmission spectra and displacement amplitude distributions. The calculated band structures and transmission spectra obtained by using the proposed LCTM and the FEM (COMSOL, 264 elements for band structures, 6284 elements for transmission spectra) are in good agreement, which is shown in Figure 14.

Cross-Shaped Scatterer Case
Here, the unit cell structure of aurum/epoxy phononic crystal consists of a square lattice and a cross-shaped scatterer and the measurements of the cross are 0.75 by 0.25. The filling fraction of the cross-shaped scatterer is 0.3125. Figures 12 and 13 show the corresponding node distributions in the proposed LCTM. In the present numerical implementation, the discretization node number is 238 N = with node spacing 1/13 h Δ = for band structure calculation, and the discretization node number is 3950 N = with the same node spacing 1/15 h Δ = for the calculation of transmission spectra and displacement amplitude distributions. The calculated band structures and transmission spectra obtained by using the proposed LCTM and the FEM (COMSOL, 264 elements for band structures, 6284 elements for transmission spectra) are in good agreement, which is shown in Figure 14.   As shown in Figure 14, it can be found that when the angular frequency is within the range of band gaps, the transmission coefficient is very low, namely, the low transmission coefficient region corresponds well to the calculated band gaps. Figure 15 shows the displacement amplitude distributions computed by the proposed LCTM and FEM + ABC/PML (COMSOL) presented at two specific normalized frequencies (   As shown in Figure 14, it can be found that when the angular frequency is within the range of band gaps, the transmission coefficient is very low, namely, the low transmission coefficient region corresponds well to the calculated band gaps. Figure 15 shows the displacement amplitude distributions computed by the proposed LCTM and FEM + ABC/PML (COMSOL) presented at two specific normalized frequencies (  As shown in Figure 14, it can be found that when the angular frequency is within the range of band gaps, the transmission coefficient is very low, namely, the low transmission coefficient region corresponds well to the calculated band gaps. Figure 15 shows the displacement amplitude distributions computed by the proposed LCTM and FEM + ABC/PML (COMSOL) presented at two specific normalized frequencies ( ωa 2πc 0 = 0.5, 0.8). It is easy to observe from Figure 15 that the SH wave is blocked at the normalized frequency ( ωa 2πc 0 = 0.5) in the bandgap region, and it cannot hinder the SH wave propagation at the normalized frequency ( ωa 2πc 0 = 0.8) in the pass band region.

Phononic Crystals with Triangular Lattice
In this subsection, another popular type of lattice, the triangular lattice, is investigated under the phononic crystals with both the infinite periodic structure and semi-infinite periodic structure in the y-direction. All the numerical results obtained by the proposed LCTM are compared with the solutions obtained by the FEM (ABC) and FEM (PML).

Circular Scatterer Case
In this case, the unit cell structure of aurum/epoxy phononic crystal consists of a hexagon matrix and a circular scatterer. The filling fraction of the circular scatterer is 0.4. Figure 16 shows the corresponding node distributions for the infinite periodic structure and semi-infinite periodic structure in the y-direction in the proposed LCTM. In the present numerical implementation, the discretization node number is N = 199 with node spacing ∆h = 1/13 for band structure calculation, and the discretization node number is N = 2808 with the same node spacing ∆h = 1/15 for the calculation of transmission spectra and displacement amplitude distributions. The calculated band structures and transmission spectra obtained by using the proposed LCTM and the FEM (COMSOL, 284 elements for band structures, 4071 elements for transmission spectra) are in good agreement, which can be found in Figure 17. The following similar conclusion can be drawn from Figure 17, that is, the low transmission coefficient region in the transmission spectra corresponds well with the calculated band gaps. Figure 18 shows the displacement amplitude distributions at two specific normalized frequencies ( ωa 2πc 0 = 0.5, 0.9) computed by the proposed LCTM and FEM + ABC/PML (COMSOL) to validate and intuitively understand the band gap characteristics of phononic crystals.

Phononic Crystals with Triangular Lattice
In this subsection, another popular type of lattice, the triangular lattice, is investigated under the phononic crystals with both the infinite periodic structure and semi-infinite periodic structure in the y-direction. All the numerical results obtained by the proposed LCTM are compared with the solutions obtained by the FEM (ABC) and FEM (PML).

Circular Scatterer Case
In this case, the unit cell structure of aurum/epoxy phononic crystal consists of a hexagon matrix and a circular scatterer. The filling fraction of the circular scatterer is 0.4. Figure 16 shows the corresponding node distributions for the infinite periodic structure and semi-infinite periodic structure in the y-direction in the proposed LCTM. In the present numerical implementation, the discretization node number is 199 N = with node spacing 1/13 h Δ = for band structure calculation, and the discretization node number is 2808 N = with the same node spacing 1/15 h Δ = for the calculation of transmission spectra and displacement amplitude distributions. The calculated band structures and transmission spectra obtained by using the proposed LCTM and the FEM (COMSOL, 284 elements for band structures, 4071 elements for transmission spectra) are in good agreement, which can be found in Figure 17. The following similar conclusion can be drawn from Figure 17, that is, the low transmission coefficient region in the transmission spectra corresponds well with the calculated band gaps. Figure 18 shows the displacement amplitude distributions at two specific normalized frequencies (

Square Scatterer Case
In this case, the unit cell structure of aurum/epoxy phononic crystal consists of a hexagon matrix and a square scatterer. The filling fraction of the square scatterer is 0.25. Figure 19 shows the corresponding node distributions for the infinite periodic structure and semi-infinite periodic structure in the y-direction in the proposed LCTM. In the present numerical implementation, the discretization node number is N = 212 with node spacing ∆h = 1/13 for band structure calculation, and the discretization node number is N = 2978 with the node spacing ∆h = 1/15 for the calculation of transmission spectra and displacement amplitude distributions. The calculated band structures and transmission spectra obtained by using the proposed LCTM and the FEM (COMSOL, 248 elements for band structures, 3357 elements for transmission spectra) are in good agreement, which is shown in Figure 20. The following similar conclusion can be drawn from Figure 20, that is, the low transmission coefficient region in the transmission spectra corresponds well with the calculated band gaps. Figure 21 shows the displacement amplitude distributions at two specific normalized frequencies ( ωa

Square Scatterer Case
In this case, the unit cell structure of aurum/epoxy phononic crystal consists of a hexagon matrix and a square scatterer. The filling fraction of the square scatterer is 0.25. Figure 19 shows the corresponding node distributions for the infinite periodic structure and semi-infinite periodic structure in the y-direction in the proposed LCTM. In the present numerical implementation, the discretization node number is 212 N = with node spacing 1/13 h Δ = for band structure calculation, and the discretization node number is 2978 N = with the node spacing 1/15 h Δ = for the calculation of transmission spectra and displacement amplitude distributions. The calculated band structures and transmission spectra obtained by using the proposed LCTM and the FEM (COMSOL, 248 elements for band structures, 3357 elements for transmission spectra) are in good agreement, which is shown in Figure 20. The following similar conclusion can be drawn from Figure 20, that is, the low transmission coefficient region in the transmission spectra corresponds well with the calculated band gaps. Figure 21 shows the displacement amplitude distributions at two specific normalized frequencies (   (b) Figure 19. Nodes distribution for (a) the unit cell of phononic crystal with a hexagon matrix including a square scatterer; (b) Node distribution with error region and element distribution with PML region for the phononic crystal with semiinfinite periodic structure in the y-direction.   Table 6. It can be observed from Table 6 that all the numerical calculations are less than 4 min, and the code of the proposed LCTM may be further improved to save CPU time.  Table 6. It can be observed from Table 6 that all the numerical calculations are less than 4 min, and the code of the proposed LCTM may be further improved to save CPU time.

Conclusions
In this paper, a novel localized collocation scheme based on T-complete functions is applied, for the first time, to calculate the band structures, transmission spectra, and displacement amplitude distribution for anti-plane transverse elastic waves in 2D solid phononic crystals. As compared with analytical solutions and COMSOL simulation, the proposed localized collocation Trefftz method (LCTM) can provide similar accurate results with less computational times for calculating the transmission spectra and displacement amplitude distribution of the simple/complicated-shaped scatterers in the square/triangular lattice. This is because the semi-analytical T-complete functions are employed in the proposed LCTM.
Under the present numerical investigations, it can be summarized that the computed wave transmission spectra for a phononic crystal with semi-infinite periodic structure are basically consistent with the corresponding band structures for a phononic crystal with infinite periodic structure. Additionally, in the design of phononic crystals, it is sufficient to use a nine-column finite period structure instead of an ideal infinite period structure to block the SH wave propagation at specific angular frequencies.
Overall, it is concluded that the proposed LCTM could be considered to be a competitive alternative for the calculation of band structures, transmission spectra, and displacement amplitude distribution in phononic crystals, after further extensive numerical and theoretical study. It is also possible, in principle, to extend the proposed LCTM to the calculation of band structures, transmission spectra, and displacement amplitude distribution in three-dimensional (3D) phononic crystals, which is under intense study and which we plan to report in a subsequent paper.

Conflicts of Interest:
The authors declare no conflict of interest.