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

Abstract

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 finite element method (FEM) results and the analytical solutions, the efficiency and accuracy of the proposed localized collocation solver are verified 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.

1. 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.

Up to now, various numerical methods have been proposed in the simulation of wave propagation in phononic crystals, such as finite difference time domain (FDTD) [14,15,16], plane wave expansion (PWE) [17,18,19,20], wavelet method [21], multiple scattering theory (MST) [22,23], finite element method (FEM) [24,25,26], boundary element method (BEM) [27,28], singular boundary method (SBM) [29], meshless methods [30], and generalized finite difference method (GFDM) [31,32], just to mention a few. Although the aforementioned numerical methods have been successfully applied to phononic crystal analysis, each of them has advantages and disadvantages. Hence, it is necessary to develop a robust and effective numerical method to analyze wave propagation in phononic crystals.

It should be mentioned that localized collocation schemes, such as the localized radial basis function collocation method (LRBFCM) [33,34,35,36,37,38], generalized finite difference method (GFDM) [39,40,41,42], RBF-based finite difference method (RBF-FD) [43,44,45,46], and localized collocation methods based on semi-analytical basis functions [47,48,49,50,51], etc., have been proposed for a numerical solution of PDEs. The common feature in these localized collocation schemes is to inherit the simple discretization formulation from the FDM [52] and also the meshless property [53,54,55,56] without generating the ill-conditioning resultant matrices. Among these localized collocation schemes, localized semi-analytical collocation methods introduced the semi-analytical basis functions to construct numerical differentiation formulations to approximate the partial derivative of physical quantities at each discretization node, which have better accuracy and efficiency as compared with the other localized collocation methods. Here, we focus on a novel localized collocation scheme based on T-complete functions, called the localized collocation Trefftz method (LCTM) and, so far, the LCTM has been successfully applied to Laplace and biharmonic equations [57,58].

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.

2. Mathematical Model of Shear Horizontal (SH) Wave Propagation by Phononic Crystals

Under the assumptions of the linear elastic wave theory [3], the propagation of time-harmonic SH waves in phononic crystals is under consideration. The corresponding governing equation [32] can be represented as follows:

$${\mu }^j\Delta u^j+{\rho }^j{\omega }^2{u}^j=0,$$

(1)

where $\omega$ stands for the angular frequency; $\mu$ and $\rho$ 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. ${\Gamma }_1$–${\Gamma }_4$(${\Gamma }_1$–${\Gamma }_6$) denote the boundaries of the square (hexagonal) unit cell and ${\Gamma }_0$ represents the interface between the matrix and the scatterer. All the quantities of the anti-plane transverse elastic wave field on periodic boundaries ${\Gamma }_1$–${\Gamma }_4$(${\Gamma }_1$–${\Gamma }_6$) should satisfy the Bloch theorem, namely,

$$\begin{array}{l}{u}^0(x+a)=e^{ika}{u}^0(x)\\ T^0(x+a)=e^{ika}{T}^0(x)\end{array},$$

(2)

where $k=(k_{{\mathbf{x}}_1},k_{{\mathbf{x}}_2})$ is the anti-plane Bloch wave vector; $T={\mu }_i\partial u_i/\partial \mathbf{n}$ denotes the traction vector; $n$ represents the unit normal vector; $i=\sqrt{-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: ${\Gamma }_0$

$$\begin{array}{l}{u}^1(x)=u^0(x),\hspace{0.17em}x\in {\Gamma }_0\\ T^1(x)=T^0(x),\hspace{0.17em}x\in {\Gamma }_0\end{array}$$

(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):

$$\mathbf{AU}=\left[\begin{array}{cc}{u}^0\left({\mathbf{x}}_{\Gamma 3}\right)& 0\\ u^0\left({\mathbf{x}}_{\Gamma 4}\right)& 0\\ \partial u^0\left({\mathbf{x}}_{\Gamma 3}\right)/\partial \mathbf{n}& 0\\ \partial u^0\left({\mathbf{x}}_{\Gamma 4}\right)/\partial \mathbf{n}& 0\\ u^0\left({\mathbf{x}}_{\Gamma 0}\right)& -u^1\left({\mathbf{x}}_{\Gamma 0}\right)\\ \partial u^0\left({\mathbf{x}}_{\Gamma 0}\right)/\partial \mathbf{n}& -\partial u^1\left({\mathbf{x}}_{\Gamma 0}\right)/\partial \mathbf{n}\\ \left(\Delta +{\kappa }_0^2\right)u^0& 0\\ 0& \left(\Delta +{\kappa }_1^2\right)u^1\end{array}\right]=\left[\begin{array}{cc}{e}^{i{k}_{\mathbf{x}1}\mathbf{a}}{u}^0\left({\mathbf{x}}_{\Gamma 1}\right)& 0\\ e^{i{k}_{\mathbf{x}2}\mathbf{a}}{u}^0\left({\mathbf{x}}_{\Gamma 2}\right)& 0\\ e^{i{k}_{\mathbf{x}1}\mathbf{a}}\partial u^0\left({\mathbf{x}}_{\Gamma 1}\right)/\partial \mathbf{n}& 0\\ e^{i{k}_{\mathbf{x}2}\mathbf{a}}\partial u^0\left({\mathbf{x}}_{\Gamma 2}\right)/\partial \mathbf{n}& 0\\ 0& 0\\ 0& 0\\ 0& 0\\ 0& 0\end{array}\right]=\lambda \mathbf{BU},$$

(4)

$$\mathbf{AU}=\left[\begin{array}{cc}\begin{array}{c}{u}^0\left({\mathbf{x}}_{\Gamma 4}\right)\\ u^0\left({\mathbf{x}}_{\Gamma 5}\right)\\ u^0\left({\mathbf{x}}_{\Gamma 6}\right)\end{array}& \begin{array}{c}0\\ 0\\ 0\end{array}\\ \partial u^0\left({\mathbf{x}}_{\Gamma 4}\right)/\partial \mathbf{n}& 0\\ \partial u^0\left({\mathbf{x}}_{\Gamma 5}\right)/\partial \mathbf{n}& 0\\ \partial u^0\left({\mathbf{x}}_{\Gamma 6}\right)/\partial \mathbf{n}& 0\\ u^0\left({\mathbf{x}}_{\Gamma 0}\right)& -u^1\left({\mathbf{x}}_{\Gamma 0}\right)\\ \partial u^0\left({\mathbf{x}}_{\Gamma 0}\right)/\partial \mathbf{n}& -\partial u^1\left({\mathbf{x}}_{\Gamma 0}\right)/\partial \mathbf{n}\\ \left(\Delta +{\kappa }_0^2\right)u^0& 0\\ 0& \left(\Delta +{\kappa }_1^2\right)u^1\end{array}\right]=\left[\begin{array}{cc}{e}^{i\left(\frac{\sqrt{3}}{2}{k}_{\mathbf{x}1}\mathbf{a}+\frac{1}{2}{\mathbf{k}}_{\mathbf{x}2}\mathbf{a}\right)}{u}^0\left({\mathbf{x}}_{\Gamma 1}\right)& 0\\ e^{i{k}_{\mathbf{x}2}\mathbf{a}}{u}^0\left({\mathbf{x}}_{\Gamma 2}\right)& 0\\ e^{i\left(-\frac{\sqrt{3}}{2}{k}_{\mathbf{x}1}\mathbf{a}+\frac{1}{2}{k}_{\mathbf{x}2}\mathbf{a}\right)}{u}^0\left({\mathbf{x}}_{\Gamma 3}\right)& 0\\ e^{i\left(\frac{\sqrt{3}}{2}{k}_{\mathbf{x}1}\mathbf{a}+\frac{1}{2}{k}_{\mathbf{x}2}\mathbf{a}\right)}\partial u^0\left({\mathbf{x}}_{\Gamma 1}\right)/\partial \mathbf{n}& 0\\ \begin{array}{c}{e}^{i{k}_{\mathbf{x}2}\mathbf{a}}\partial u^0\left({\mathbf{x}}_{\Gamma 2}\right)/\partial \mathbf{n}\\ e^{i\left(-\frac{\sqrt{3}}{2}{k}_{\mathbf{x}1}\mathbf{a}+\frac{1}{2}{k}_{\mathbf{x}2}\mathbf{a}\right)}\partial u^0\left({\mathbf{x}}_{\Gamma 1}\right)/\partial \mathbf{n}\\ 0\end{array}& \begin{array}{c}0\\ 0\\ 0\end{array}\\ 0& 0\\ 0& 0\\ 0& 0\end{array}\right]=\lambda \mathbf{BU},$$

(5)

where ${\kappa }_j=\omega \sqrt{{\rho }^j/{\mu }^j}$. The band structure can be obtained by evaluating the dispersion relations along the boundary of the irreducible Brillouin zone, in which $\lambda =e^{i{k}_{\mathbf{x}1}\mathbf{a}},\hspace{0.17em}{e}^{i{k}_{\mathbf{x}2}\mathbf{a}}=1$ from $\Gamma -\mathrm{X}$; $\lambda =e^{i{k}_{\mathbf{x}2}\mathbf{a}},\hspace{0.17em}{e}^{i{k}_{\mathbf{x}1}\mathbf{a}}=-1$ from $\mathrm{X}-\mathrm{M}$; $\lambda =e^{i{k}_{\mathbf{x}1}\mathbf{a}}=e^{i{k}_{\mathbf{x}2}\mathbf{a}}$ from $\mathrm{M}-\Gamma$ in square lattice and $\lambda =e^{i{k}_{\mathbf{x}2}\mathbf{a}},\hspace{0.17em}{e}^{i{k}_{\mathbf{x}1}\mathbf{a}}=1$ from $\Gamma -\mathrm{X}$; $\lambda =e^{i{k}_{\mathbf{x}2}\mathbf{a}}=e^{i(\frac{4\pi }{3\mathbf{a}}-\frac{1}{\sqrt{3}}{k}_{\mathbf{x}1})\mathbf{a}}$ from $\mathrm{X}-\mathrm{M}$; and $\lambda =e^{i{k}_{\mathbf{x}2}\mathbf{a}}=e^{i\sqrt{3}{k}_{\mathbf{x}1}\mathbf{a}}$ from $\mathrm{M}-\Gamma$ in triangular lattice. Then, the corresponding eigen frequency, $\lambda$, can be obtained by solving the generalized eigenvalue equations with sweeping the angular frequency, $\omega$, 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 $\Gamma \mathrm{X}-$ and $\Gamma \mathrm{M}-$ directions. Then, the displacements and tractions on the boundaries ${\Gamma }_2^{*},\hspace{0.17em}\hspace{0.17em}{\Gamma }_4^{*}$ satisfy the following relationships:

$$u_{{\Gamma }_4^{*}}^0=u_{{\Gamma }_2^{*}}^0,\hspace{0.17em}{T}_{{\Gamma }_4^{*}}^0=T_{{\Gamma }_2^{*}}^0,$$

(6)

and the left boundary ${\Gamma }_1^{*}$ is excited by the following unit displacement:

$$u_{{\Gamma }_1^{*}}^0=1,$$

(7)

and the right boundary ${\Gamma }_3^{*}$ is imposed on the following absorbing boundary condition:

$$\mathbf{n}\cdot \nabla u^0=-i{\kappa }_0{u}^0.$$

(8)

By considering the interface continuity conditions (3) and boundary conditions (6)–(8), the following set of linear equations can be obtained:

$${\mathbf{A}}^{*}\mathbf{U}=\left[\begin{array}{cc}{u}^0\left({\mathbf{x}}_{{\Gamma }_1^{*}}\right)& 0\\ u^0\left({\mathbf{x}}_{{\Gamma }_4^{*}}\right)-u^0\left({\mathbf{x}}_{{\Gamma }_2^{*}}\right)& 0\\ \partial u^0\left({\mathbf{x}}_{{\Gamma }_3^{*}}\right)/\partial \mathbf{n}+i{\kappa }_0{u}^0\left({\mathbf{x}}_{{\Gamma }_3^{*}}\right)& 0\\ \partial u^0\left({\mathbf{x}}_{{\Gamma }_4^{*}}\right)/\partial \mathbf{n}-\partial u^0\left({\mathbf{x}}_{{\Gamma }_2^{*}}\right)/\partial \mathbf{n}& 0\\ u^0\left({\mathbf{x}}_{{\Gamma }_0^{*}}\right)& -u^1\left({\mathbf{x}}_{{\Gamma }_0^{*}}\right)\\ \partial u^0\left({\mathbf{x}}_{{\Gamma }_0^{*}}\right)/\partial \mathbf{n}& -\partial u^1\left({\mathbf{x}}_{{\Gamma }_0^{*}}\right)/\partial \mathbf{n}\\ \left(\Delta +{\kappa }_0^2\right)u^0& 0\\ 0& \left(\Delta +{\kappa }_1^2\right)u^1\end{array}\right]=\left[\begin{array}{c}1\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\end{array}\right]=\mathbf{b},$$

(9)

where wavenumber ${\kappa }_0=\omega \sqrt{{\rho }^0/{\mu }^0}$, ${\kappa }_1=\omega \sqrt{{\rho }^1/{\mu }^1}$. By solving Equation (9), the corresponding transmitted displacements on the right boundary $u^0({\mathbf{x}}_{{\Gamma }_3^{*}})$ can be obtained for every given angular frequency $\omega$. Then, the transmission spectra can be obtained by computing the transmission coefficient $\mathrm{T}=20\times \log \left(u^0\left({\mathbf{x}}_{{\Gamma }_3^{*}}\right)/\left(2\times {10}^{-5}\times u^0\left({\mathbf{x}}_{{\Gamma }_1^{*}}\right)\right)\right)$ for different angular frequency $\omega$.

3. 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 ${\mathbf{x}}_0^i$ with m nearest nodes $\left({\mathbf{x}}_1^i,\hspace{0.17em}{\mathbf{x}}_2^i,\hspace{0.17em}\dots \hspace{0.17em},{\mathbf{x}}_m^i\right)$ around $x_0^i$ form a subdomain ${\Xi }_i$, whose center can be chosen as ${\tilde{\mathbf{x}}}^i=\frac{1}{m+1}{\displaystyle \sum _{j=0}^m{\mathbf{x}}_j^i}$ shown in Figure 3. Then, the approximated displacement $\overline{u}\left({\mathbf{x}}_j^i\right)$ inside the stencil support ${\Xi }_i$ can be represented by a linear combination of T-complete functions ${\varphi }_k^i$ with unknown coefficients ${\alpha }_k^i$

$$\overline{u}\left({\mathbf{x}}_j^i\right)={\displaystyle \sum _{k=0}^m{\varphi }_k^i{\alpha }_k^i}\mathrm{with}\mathrm{its}\mathrm{matrix}\mathrm{form}\overline{u}={\mathbf{\Phi }}^i{\alpha }^i,$$

(10)

where the following T-complete functions are used to satisfy the governing Equation (1) in advance:

$$\begin{array}{l}{\varphi }_k^i=\{J_0\left({\kappa }_1{r}_k^i\right),\hspace{0.17em}{J}_1\left({\kappa }_1{r}_k^i\right)\cos \left({\theta }_k^i\right),J_1\left({\kappa }_1{r}_k^i\right)\sin \left({\theta }_k^i\right),\hspace{0.17em}\dots \hspace{0.17em},\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{J}_{m/2}\left({\kappa }_1{r}_k^i\right)\cos \left(m{\theta }_k^i/2\right),J_{m/2}\left({\kappa }_1{r}_k^i\right)\sin \left(m{\theta }_k^i/2\right)\}\end{array}$$

(11)

where $\left(r_k^i,{\theta }_k^i\right)$ is constructed based on the polar coordinate system determined by the nodes in a subdomain ${\Xi }_i$ with the center ${\tilde{\mathbf{x}}}^i$ as the origin and $J_i$ is Bessel function of the first kind of order i. Then, the vector of unknown coefficients ${\alpha }^i$ can be expressed as follows:

$${\alpha }^i={\left({\mathbf{\Phi }}^i\right)}^{-1}{\overline{u}}_{}^i.$$

(12)

Next, the approximate solution ${\overline{u}}_0^i$ and traction ${\overline{T}}_0^i$ at ${\mathbf{x}}_0^i$ can be represented as follows:

$${\overline{u}}_0^i={\mathbf{\Phi }}_0^i{\alpha }^i={\mathbf{\Phi }}_0^i{\left({\mathbf{\Phi }}_{}^i\right)}^{-1}{\overline{u}}_{}^i={\mathbf{W}}^i{\overline{u}}_{}^i={\displaystyle \sum _{j=0}^m{w}_j^i{\overline{u}}_j^i},$$

(13)

$${\overline{T}}_0^i={\mu }^i\frac{\partial {\overline{u}}_0^i}{\partial \mathbf{n}}={\mu }^i\frac{\partial {\mathbf{\Phi }}_0^i}{\partial \mathbf{n}}{\alpha }^i={\mu }^i\frac{\partial {\mathbf{\Phi }}_0^i}{\partial \mathbf{n}}{\left({\mathbf{\Phi }}_{}^i\right)}^{-1}{\overline{u}}_{}^i={\overline{\mathbf{W}}}^j{\overline{u}}_{}^i={\displaystyle \sum _{j=0}^m{\overline{w}}_j^i{\overline{u}}_j^i}.$$

(14)

It should be mentioned, because the used T-complete function satisfies the governing equation in advance, the LCTM numerical differentiation formulation of ${\overline{\mathbf{u}}}_0^i$ at ${\mathbf{x}}_0^i$ inside the domain of matrix and scatterer can be rewritten as follows:

$${\overline{\mathbf{u}}}_0^i-{\displaystyle \sum _{j=0}^m{w}_j^i{\overline{u}}_j^i=0}.$$

(15)

For a given angular frequency $\omega$, 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.

4. 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\hspace{0.17em}m/s$, ${\rho }_0=19500\hspace{0.17em}\hspace{0.17em}kg/m^3$, $c_1=1161\hspace{0.17em}m/s$, ${\rho }_1=1180\hspace{0.17em}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.

4.1. 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 ${\varphi }_I=e^{ikr\cos {\theta }_0}$, the analytical solutions of the scattered field $u^s$ [59] for rigid boundary is as follows:

$$u^s(r,\theta )=-\frac{{J^{\prime }}_0(ka)}{{H^{\prime }}_0^{}(ka)}{H}_0^{}(kr)-2\sum _{n=1}^{\infty }{i}^n\frac{{J^{\prime }}_n(ka)}{{H^{\prime }}_n^{}(ka)}{H}_n^{}(kr)\cos n\theta ,$$

(16)

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 ${\Gamma }_{\infty }$, and the weakened Sommerfeld boundary condition is used as the following absorbing boundary conditions (ABC) [59]:

$$\frac{\partial u}{\partial r}-iku=0,\left(r,\theta \right)\in {\Gamma }_{\infty }.$$

(17)

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 ${\Gamma }_{\infty }$ 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₂ relative error, which is defined as follow:

$${\mathrm{L}}_2=\sqrt{{\displaystyle \sum _{i=1}^N{\epsilon }_i^2}/{\displaystyle \sum _{i=1}^N{u}_i^2}},\mathrm{in}\mathrm{which}{\epsilon }_i=\left|\left(\overline{u}\left(x_i\right)-u\left(x_i\right)\right)\right|,$$

(18)

where $N$ is the number of total discretization nodes and $\overline{u}\left(x_i\right),u\left(x_i\right)$ denotes the numerical solution and analytical solution at $x_i$.

Table 1. Numerical errors (L₂) in the reliable region obtained by the proposed localized collocation Trefftz method (LCTM) + absorbing boundary conditions (ABC) with respect to the increasing node number and wavenumber.

Node Number (Node Density)	k = 1	k = 5	k = 10
998 ($\Delta h$ = 0.57)	7.54 × 10−3	/	/
2025 ($\Delta h$ = 0.395)	6.73 × 10−3	/	/
5032 ($\Delta h$ = 0.25)	5.99 × 10−3	/	/
10,840 ($\Delta h$ = 0.17)	5.99 × 10−3	1.71 × 10−3	/
25,069 ($\Delta h$ = 0.1117)	6.00 × 10−3	1.62 × 10−3	/
50,010 ($\Delta h$ = 0.079)	5.97 × 10−3	1.62 × 10−3	1.93 × 10−3
80,108 ($\Delta h$ = 0.0624)	5.95 × 10−3	1.62 × 10−3	1.30 × 10−3
119,863 ($\Delta h$ = 0.051)	5.94 × 10−3	1.62 × 10−3	1.30 × 10−3

and

Table 2. Numerical errors (L₂) in reliable region obtained by the FEM + ABC (COMSOL) with respect to the increasing element number and wavenumber.

Element Number	k = 1	k = 5	k = 10
1016	7.18 × 10−3	/	/
2010	7.03 × 10−3	/	/
5079	6.73 × 10−3	/	/
9878	6.57 × 10−3	7.35 × 10−3	/
25,805	6.42 × 10−3	2.05 × 10−3	/
51,742	6.37 × 10−3	1.99 × 10−3	9.57 × 10−3
81,503	6.37 × 10−3	1.99 × 10−3	3.90 × 10−3
120,266	6.37 × 10−3	1.99 × 10−3	1.95 × 10−3

list the numerical errors (L₂ 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₂ relative error is larger than 0.01 and $\Delta h$ stands for the average node spacing. It can be found from Table 1 and Table 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 $\overline{u}\left({\mathbf{x}}_j^i\right)$ can be represented by a linear combination of Hankel functions as:

$$\overline{u}(R_j^i,{\vartheta }_j^i)={\gamma }_0{H}_0\left(k{R}_j^i\right)+{\displaystyle \sum _{n=1}^{n_{ABC}}\left({\gamma }_{2n}\cos n{\vartheta }_j^i+{\gamma }_{2n-1}\sin n{\vartheta }_j^i\right)H_n\left(k{R}_j^i\right)},$$

(19)

in which ${\gamma }_0,{\gamma }_1,\cdots ,{\gamma }_{2n_{ABC}-1},{\gamma }_{2n_{ABC}}$ are underdetermined coefficients and $\left(R_j^i,{\vartheta }_j^i\right)$ is constructed based on the polar coordinate system. To construct the ABCWEF formulation for each given $i$-th node ${\mathbf{x}}_0^i$ on the absorbing boundary ${\Gamma }_{\infty }$, its m nearest nodes $\left({\mathbf{x}}_1^i,\hspace{0.17em}{\mathbf{x}}_2^i,\hspace{0.17em}\dots \hspace{0.17em},{\mathbf{x}}_m^i\right)$ in the subdomain ${\Xi }_i$ are used, whose center is the center of the computational domain. By using the wave expansion formulation (19) at ${\mathbf{x}}_0^i$ on the absorbing boundary ${\Gamma }_{\infty }$ and its m nearest nodes $\left({\mathbf{x}}_1^i,\hspace{0.17em}{\mathbf{x}}_2^i,\hspace{0.17em}\dots \hspace{0.17em},{\mathbf{x}}_m^i\right)$ in the subdomain ${\Xi }_i$, the following system of linear equations can be obtained:

$${\overline{u}}^i={\mathbf{H}}_{}^i{\gamma }_{}^i,$$

(20)

and then the approximate displacement $\overline{u}\left({\mathbf{x}}_0^i\right)$ can be written as follow:

$$\overline{u}\left(x_0^i\right)=\left[H_0\left(k{R}_0^i\right),\cdots ,H_{2n_{ABC}}\left(k{R}_0^i\right)\cos n_{ABC}{\vartheta }_j^i,H_{2n_{ABC}}\left(k{R}_0^i\right)\sin n_{ABC}{\vartheta }_j^i\right]{\left({\mathbf{H}}_{}^i\right)}^{-1}{\overline{u}}^i,$$

(21)

substituting Equation (21) into Equation (13), the following formulation is obtained:

$$\left(\left[H_0\left(k{R}_0^i\right),\cdots ,H_{2n_{ABC}}\left(k{R}_0^i\right)\cos n_{ABC}{\vartheta }_j^i,H_{2n_{ABC}}\left(k{R}_0^i\right)\sin n_{ABC}{\vartheta }_j^i\right]{\left({\mathbf{H}}_{}^i\right)}^{-1}-{\mathbf{W}}^i\right){\overline{u}}^i=0.$$

(22)

Table 3. Numerical errors (L₂) in the reliable region obtained by the proposed LCTM + absorbing boundary conditions based on wave expansion functions (ABCWEF) with respect to the increasing node number and wavenumber.

Node Number (Node Density)	k = 1	k = 5	k = 10
998 ($\Delta h$ = 0.57)	8.59 × 10−3	/	/
2025 ($\Delta h$ = 0.395)	1.26 × 10−3	/	/
5032 ($\Delta h$ = 0.25)	1.75 × 10−4	/	/
10,840 ($\Delta h$ = 0.17)	1.02 × 10−4	3.18 × 10−4	/
25,069 ($\Delta h$ = 0.1117)	9.74 × 10−5	4.00 × 10−5	/
50,010 ($\Delta h$ = 0.079)	5.26 × 10−5	1.85 × 10−5	4.88 × 10−5
80,108 ($\Delta h$ = 0.0624)	2.89 × 10−5	1.34 × 10−5	1.23 × 10−5
119,863 ($\Delta h$ = 0.051)	2.81 × 10−5	6.42 × 10−6	5.78 × 10−6

lists the numerical errors (L₂ 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. Numerical errors (L₂) obtained by the FEM + PML (COMSOL) with respect to the increasing element number and wavenumber.

Element Number	k = 1	k = 5	k = 10
1016	1.99 × 10−3	/	/
2010	7.03 × 10−4	/	/
5079	2.50 × 10−4	/	/
9878	1.25 × 10−4	7.78 × 10−3	/
25,805	3.36 × 10−5	1.20 × 10−3	/
51,742	2.50 × 10−5	3.31 × 10−3	8.39 × 10−3
81,503	2.35 × 10−5	1.34 × 10−4	3.63 × 10−3
120,266	2.32 × 10−5	6.38 × 10−5	1.68 × 10−3

, 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 Table 3 and Table 4, as compared with the FEM + PML, the proposed LCTM + ABCWEF can provide more accurate results with less computational resources.

In

Table 5. Condition number Cf and Ca and time cost of LCTM and FEM for each node number.

Element Number	Ca	Cf	CPU Time of LCTM	CPU Time of FEM
998 ($\Delta h$ = 0.57)	6.26 × 105	5.64 × 1050	1 s	1 s
2025 ($\Delta h$ = 0.395)	7.81 × 105	9.14 × 1052	1 s	1 s
5032 ($\Delta h$ = 0.25)	6.84 × 106	1.06 × 1055	2 s	2 s
10,840 ($\Delta h$ = 0.17)	1.19 × 107	1.05 × 1060	5 s	2 s
25,069 ($\Delta h$ = 0.1117)	8.99 × 106	3.81 × 1063	13 s	4 s
50,010 ($\Delta h$ = 0.079)	6.04 × 107	4.08 × 1064	22 s	9 s
80,108 ($\Delta h$ = 0.0624)	3.02 × 108	6.38 × 1066	42 s	11 s
119,863 ($\Delta h$ = 0.051)	9.28 × 108	2.79 × 1068	68 s	16 s

, the condition numbers of ${\mathbf{\Phi }}^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 ${\mathbf{\Phi }}^i$, the Moore–Penrose inverse technique is implemented to calculate ${\left({\mathbf{\Phi }}^i\right)}^{-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.

4.2. 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.

4.2.1. 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 $\Delta 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.

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 $\mathsf{\Gamma }\mathrm{X}$- direction. 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 $\Delta 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 ($\frac{\omega a}{2\pi 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 ($\frac{\omega a}{2\pi c_0}=0.5$) in the bandgap region, and it cannot hinder the SH wave propagation at the normalized frequency ($\frac{\omega a}{2\pi c_0}=1$) in the pass band region.

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.

4.2.2. 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. Figure 12 and Figure 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 $\Delta h=1/13$ for band structure calculation, and the discretization node number is $N=3950$ with the same node spacing $\Delta 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.

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 ($\frac{\omega a}{2\pi c_0}=0.5,0.8$). It is easy to observe from Figure 15 that the SH wave is blocked at the normalized frequency ($\frac{\omega a}{2\pi c_0}=0.5$) in the bandgap region, and it cannot hinder the SH wave propagation at the normalized frequency ($\frac{\omega a}{2\pi c_0}=0.8$) in the pass band region.

4.3. 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).

4.3.1. 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 $\Delta h=1/13$ for band structure calculation, and the discretization node number is $N=2808$ with the same node spacing $\Delta 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 ($\frac{\omega a}{2\pi 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.

4.3.2. 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 $\Delta h=1/13$ for band structure calculation, and the discretization node number is $N=2978$ with the node spacing $\Delta 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 ($\frac{\omega a}{2\pi c_0}=0.6,1.48$) computed by the proposed LCTM and FEM + ABC/PML (COMSOL) to validate and intuitively understand the band gap characteristics of phononic crystals.

Finally, computational times of the proposed LCTM and the FEM in the solution of the cases in Section 4.2 and Section 4.3 are listed in

Table 6. CPU times by using the proposed LCTM + ABC and FEM + ABC in Section 4.2 and Section 4.3.

Cases	Transmission by LCTM	Transmission by COMSOL	Band Structure by LCTM	Band Structure by COMSOL
4.2.1	133 s	32 s	68 s	58 s
4.2.2	197 s	29 s	120 s	47 s
4.3.1	146 s	18 s	58 s	56 s
4.3.2	134 s	16 s	95 s	62 s

. 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.

5. 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.