*Article* **A Numerical Method for Flexural Vibration Band Gaps in A Phononic Crystal Beam with Locally Resonant Oscillators**

#### **Xu Liang 1, Titao Wang 1,\*, Xue Jiang 2,\*, Zhen Liu 1, Yongdu Ruan <sup>1</sup> and Yu Deng <sup>1</sup>**


Received: 27 April 2019; Accepted: 3 June 2019; Published: 5 June 2019

**Abstract:** The differential quadrature method has been developed to calculate the elastic band gaps from the Bragg reflection mechanism in periodic structures efficiently and accurately. However, there have been no reports that this method has been successfully used to calculate the band gaps of locally resonant structures. This is because, in the process of using this method to calculate the band gaps of locally resonant structures, the non-linear term of frequency exists in the matrix equation, which makes it impossible to solve the dispersion relationship by using the conventional matrix-partitioning method. Hence, an accurate and efficient numerical method is proposed to calculate the flexural band gap of a locally resonant beam, with the aim of improving the calculation accuracy and computational efficiency. The proposed method is based on the differential quadrature method, an unconventional matrix-partitioning method, and a variable substitution method. A convergence study and validation indicate that the method has a fast convergence rate and good accuracy. In addition, compared with the plane wave expansion method and the finite element method, the present method demonstrates high accuracy and computational efficiency. Moreover, the parametric analysis shows that the width of the 1st band gap can be widened by increasing the mass ratio or the stiffness ratio or decreasing the lattice constant. One can decrease the lower edge of the 1st band gap by increasing the mass ratio or decreasing the stiffness ratio. The band gap frequency range calculated by the Timoshenko beam theory is lower than that calculated by the Euler-Bernoulli beam theory. The research results in this paper may provide a reference for the vibration reduction of beams in mechanical or civil engineering fields.

**Keywords:** phononic crystal; locally resonant; band gap; differential quadrature method

#### **1. Introduction**

Phononic crystals (PCs) are periodic composites or structures which modify the band structure in some way. The band gap is a frequency range in which the propagation of elastic waves in phononic crystal structures is suppressed. By adjusting the parameters of the artificial periodic structure, the position and width of the band gap and its ability to suppress wave propagation can be artificially regulated. In engineering practice, structures can be designed as phononic crystals, and then the band gap characteristics can be used in vibration and noise reduction.

There are two formation mechanisms of the elastic band gap in PCs: one is the Bragg scattering mechanism [1–18], and the other is the locally resonant mechanism. The wave length corresponding to the elastic band gap formed by Bragg scattering is generally equal to the lattice size or lattice constant, which restricts its application in engineering practice. In 1999, Liu et al. [19] proposed the locally resonant mechanism of the elastic band gap. It was found that the corresponding wave length of this locally resonant phononic crystal is far larger than the lattice size, which breaks through the limitation of the Bragg scattering mechanism and opens up a wide potential for application in low-frequency wave band. Since then, investigations on the locally resonant mechanism of phononic crystals have been carried out continuously [20–23].

The study of the band gap mechanism of PCs depends on the solving methods of elastic wave band gaps. At present, computational methods of vibration band gaps mainly include the transfer-matrix method (TM), multiple-scattering theory (MST), finite difference time-domain method (FDTD), lumped-mass method (LM), plane wave expansion method (PWE), differential quadrature method (DQM), finite element method (FEM), extended plane wave expansion method (EPWE) [24], wave finite element method [25], and boundary element method (BE) [26]. The transfer-matrix method [27,28] is widely applied to calculate the band gap characteristics of 1-D PCs. Although the analytical solution can be obtained quickly, it is not suitable for the study of the vibration dispersion relations of 2-D and 3-D periodic structures, since the transfer matrix can usually only be transmitted in one direction. The PWE method [29,30] is the most basic and common method to calculate the band gap characteristics of phononic crystals. It can be used to solve the elastic band gap of all-dimensional PCs. However, when the material parameters vary greatly, or the filling rate is too high or too low, it is difficult to achieve convergence [31]. It is worth mentioning that in order to overcome the limitation of the convergence of PWE, the improved plane wave expansion (IPWE) method [26] is proposed. The FDTD method [32,33] can calculate the transmission, reflection, and energy band characteristics of infinite structures. However, the computational amount is large, and the large elastic constant difference may lead to numerical instability and divergence. The multiple scattering theory [34,35] can not only calculate the dispersion curves of periodic structures, but also disordered structures. It has fast convergence and high accuracy and is easy to deal with the problem of elastic mismatch, but only some kinds of scatterers with regular shape can be dealt with. At present, only phononic crystals composed of cylindrical (two-dimensional) and spherical (three-dimensional) scatterers can be calculated. The lumped mass method [36,37] converges fast, but there are some difficulties in dealing with multi-field coupling problems. The FEM method [38–40] is a commonly used method in engineering with good applicability, wide application range, and good convergence. There are various kinds of mature commercial software, such as Comsol, ANSYS, etc., which facilitate the modeling and analysis of complex periodic structures. The finite element method can be used to accurately calculate the band gap characteristics of PCs of various dimensions and various shapes of scatterers.

The differential quadrature method can approximate the value of each derivative of the function at the node with the weighted sum of the function values at all nodes in the domain and transforms the problem of the continuous system into a discrete problem. It has strong convergence and high precision. To the authors' knowledge, the DQM was only applied to solve the elastic wave band gap of a beam or plate structure with the Bragg scattering mechanism [8,9]. However, the Bloch boundary conditions for locally resonant structures lead to nonlinear fundamental equations when using DQM, which causes difficulties in solving. Hence, numerical solutions of band gaps in LR structures based on DQM have not been reported up to now.

In this work, through applying an unconventional matrix-partitioning method and a variable substitution method, we transform the problem of solving the dispersion relation into a quadratic eigenvalue problem and propose a numerical method based on the differential quadrature method to calculate the bending vibration band gap of locally resonant beams. The purpose is to improve calculation accuracy and computational efficiency. Moreover, a parametric study is undertaken to investigate the effects of shear deformation and rotary inertia, the lumped mass, the spring stiffness coefficient, and the lattice size on the 1st band gap. Some major novelties of the contribution are pointed out as follows:


#### **2. Method**

#### *2.1. Di*ff*erential Quadrature Method*

The basic idea of the differential quadrature method is to use the sum of the weighted values at all discrete points in the computational domain to approximate the unknown function values and their derivatives at any discrete points, so that the continuous system can be transformed into a discrete system for solving as follows.

$$f(\boldsymbol{\varepsilon}) = \sum\_{j=1}^{N\_{\boldsymbol{\varepsilon}}} p\_j(\boldsymbol{\varepsilon}) f\_{j\boldsymbol{\varepsilon}} \tag{1}$$

$$\frac{\partial^r f(\varepsilon\_i)}{\partial \mathbf{x}^r} = \sum\_{i=1}^{N\_\mathbf{x}} A\_{ij}^{(r)} \cdot f(\varepsilon\_i), \tag{2}$$

where *Nx* is the total number of discrete points in the calculation domain; *i*, *j* = 1, 2, ... , *Nx*; *f*(ε*j*) = *fj*; *pj*(ε) are the Lagrange interpolation polynomials, and *A<sup>r</sup> ij* is the weight coefficient of the *<sup>r</sup>*th derivative defined by [41]:

$$\mathcal{A}\_{ij}^{(1)} = \frac{\prod\_{r=1, r \neq i}^{N\_x} (\varepsilon\_i - \varepsilon\_r)}{(\varepsilon\_i - \varepsilon\_j) \prod\_{r=1, r \neq j}^{N\_x} (\varepsilon\_j - \varepsilon\_r)},\tag{3}$$

$$A\_{ij}^{(r)} = r \left\{ A\_{ij}^{(r-1)} A\_{ij}^1 - \frac{A\_{ij}^{(r-1)}}{x\_i - x\_j} \right\} \tag{4}$$

where *r* = 1, 2, ... , *Nx*−1; *i*,*j* = 1, 2, ... , *Nx*, (*i <sup>j</sup>*); and *<sup>A</sup>*(*r*) *ii* are defined as:

$$A\_{ii}^{(r)} = -\sum\_{j=1, j \neq i}^{N\_x} A\_{ij}^{(r)}{}\_{,} \tag{5}$$

In order to obtain higher accuracy and faster convergence, the non-uniform mesh partition (Gauss-Lobatto-Chebyshev pattern) [42] is adopted in this paper. The coordinates of discrete points are as follows:

$$\varepsilon\_{i} = -\cos\left(\frac{\pi(i-1)}{N\_{\text{x}}-1}\right), i = 1, 2, \dots, N\_{\text{x}} \tag{6}$$

#### *2.2. Unconventional Matrix-Partitioning Method & Variable Substitution Method*

In the process of using the DQM method to calculate the band gaps of a locally resonant structure, a non-linear term of frequency exists in the matrix equation, which makes it impossible to solve the

dispersion relationship by using the conventional matrix-partitioning method. Hence, we propose an unconventional matrix-partitioning method and a variable substitution method to calculate the flexural band gap of a locally resonant beam. The idea is to separate the nonlinear term from the matrix by unconventional matrix partitioning and variable substitution, and then transform the matrix equation into a quadratic eigenvalue problem. The proposed method is applicable to both the Euler–Bernoulli beam model and the Timoshenko beam model. The detailed application process is shown in Section 3.2. The following is a concise and general formulation of the proposed method.

$$
\begin{bmatrix}
\mathbf{K}\_{qq}(\omega) & \mathbf{K}\_{qb} & \mathbf{K}\_{qd} \\
\mathbf{K}\_{bq} & \mathbf{K}\_{bb} & \mathbf{K}\_{bd} \\
\mathbf{K}\_{dq} & \mathbf{K}\_{db} & \mathbf{K}\_{dd}
\end{bmatrix} \cdot \begin{pmatrix}
\iota \mathbf{U}\_{q} \\
\mathbf{U}\_{b} \\
\mathbf{U}\_{d}
\end{pmatrix} - \omega^{2} \begin{bmatrix}
\mathbf{0} & \mathbf{0} & \mathbf{0} \\
\mathbf{0} & \mathbf{0} & \mathbf{0} \\
\mathbf{0} & \mathbf{0} & \mathbf{M}\_{dd}
\end{bmatrix} \cdot \begin{pmatrix}
\iota \mathbf{U}\_{q} \\
\mathbf{U}\_{b} \\
\mathbf{U}\_{d}
\end{pmatrix} = \mathbf{0},\tag{7}
$$

where **U** is a vector of independent variables. Subscripts "q", "b", and "d" refer to the shear boundary condition, other boundary conditions, and the domain, respectively. K*qq*(ω) is a scalar with the non-linear term. By variable substitutions and matrix operations, Equation (7) can be transformed into a quadratic eigenvalue problem.

$$r\_z^2 \mathbf{H}\_2 \cdot \mathbf{U}\_d + r\_z \mathbf{H}\_1 \cdot \mathbf{U}\_d + \mathbf{H}\_0 \cdot \mathbf{U}\_d = 0,\tag{8}$$

$$\mathbf{H}\_f = -\begin{bmatrix} \mathbf{0} & \mathbf{H}\_2 \\ -\mathbf{I} & \mathbf{0} \end{bmatrix}^{-1} \cdot \begin{bmatrix} \mathbf{H}\_0 & \mathbf{H}\_1 \\ \mathbf{0} & \mathbf{I} \end{bmatrix} \tag{9}$$

$$
\omega^2 = (-1/m\_z c\_4) r\_z + (k\_z/m\_z + c\_1/m\_z c\_4),
\tag{10}
$$

For any given wave vector *k* in the first Brillouin zone, *rz* can be obtained by calculating the eigenvalue of matrix *Hf*. One can get the circular frequency ω due to Equation (10), thus the dispersion relation and bending vibration band gaps can be plotted.

#### **3. Locally Resonant (LR) Beam Models and Solutions**

Figure 1 shows the configuration of a straight elastic metamaterial beam with periodical locally resonant (LR) oscillator structures. Harmonic locally resonant oscillators are periodically connected along the x-axis to the infinite Euler beam. Each oscillator is formed by a lumped mass *mz* and a spring with stiffness *kz*, taking the distance between two adjacent LR oscillators as lattice size *a*. By extending the Euler-Bernoulli beam theory and the Timoshenko beam theory, one can obtain the theoretical model of the LR beam. Since the structure has infinite periodicity, we can apply the Floquet-Bloch theorem to simplify the whole model into a unit cell.

**Figure 1.** (**a**) Configuration of a straight elastic metamaterial beam with locally resonant (LR) oscillators. (**b**) Diagram of the force equilibrium of the (*n* + 1)th LR oscillator.

#### *3.1. Euler–Bernoulli Model & Solution Procedures*

When the length of each beam unit cell is much larger than its height and width, the Euler-Bernoulli approximation is satisfied, thus the influences of shear force and rotary inertia can be ignored. The governing equation for the bending vibration of the *n*th cell is shown below [43]:

$$EI\frac{\partial^4 w\_n(\mathbf{x}, t)}{\partial \mathbf{x}^4} + \rho A \frac{\partial^2 w\_n(\mathbf{x}, t)}{\partial t^2} = 0,\tag{11}$$

where ρ is the density; *E* is Young's modulus; *A* is the cross-section area; *I* is the area moment of inertia with respect to the axis perpendicular to the beam axis, and *w*(*x*,*t*) is the lateral displacement at *x*. By assuming *w*(*x*, *t*) = *W*(*x*) exp(−*i*ω*t*), where *W*(*x*) is the vibration amplitude of the beam at *x*, and ω is the circular frequency, Equation (11) can be rewritten as follows:

$$EI\frac{\partial^4 \mathcal{W}\_\mathcal{U}(\mathbf{x})}{\partial \mathbf{x}^4} - \omega^2 \rho A \mathcal{W}\_\mathcal{U}(\mathbf{x}) = 0,\tag{12}$$

For the (*n* + 1)th locally resonant oscillator, consider the balance of forces in the *y*-axis direction, we can get:

$$\dot{Z}\_{n+1}(t) - m\_z \ddot{Z}\_{n+1}(t) = 0,\tag{13}$$

where *fn*+1(*t*) is the interaction force between the beam and the oscillator at node *xn*+1, *Zn*+1(*t*) = *Vn*+<sup>1</sup> exp(−*i*ω*t*) is the displacement of the lumped mass of the (*n* + 1)th oscillator, and the vibration amplitude of the (*n* + 1)th oscillator is denoted by the absolute value of *Vn*+1.

According to Hooke's law, *fn*+1(*t*) can be expressed as follows:

$$f\_{n+1}(t) = k\_z[w(\mathbf{x}\_{n+1}, t) - Z\_{n+1}(t)] = k\_z[W\_{n+1}(0) - V\_{n+1}] \exp(-i\omega t) \triangleq F\_{n+1} \exp(-i\omega t), \tag{14}$$

Substituting Equation (14) into Equation (13), one can get:

$$W\_{n+1} = \frac{k\_z}{k\_z - m\_z \omega^2} W\_{n+1}(0),\tag{15}$$

Ignoring the stress concentration between two adjacent units, the following boundary conditions can be listed according to the continuity of displacement, angle of rotation, bending moment, and shear force at the node *xn*+1.

$$\mathcal{W}\_{n+1}(0) = \mathcal{W}\_n(a),\tag{16}$$

$$\mathcal{W}'\_{n+1}(0) = \mathcal{W}'\_n(a),\tag{17}$$

$$EI\mathcal{W}\_{n+1}^{\prime\prime}(0) = EI\mathcal{W}\_{n}^{\prime\prime}(a),\tag{18}$$

$$EIN\_{n+1}^{\prime\prime\prime}(0) - F\_{n+1} = EIN\_n^{\prime\prime\prime}(a),\tag{19}$$

According to the Floquet-Bloch theorem [44], Equations (16)–(19) can be rewritten as follows:

$$
\epsilon^{\rm ika} \mathcal{W}\_n(0) = \mathcal{W}\_n(a),
\tag{20}
$$

$$
\epsilon^{\rm iku} \mathcal{W}'\_n(0) = \mathcal{W}'\_n(a),
\tag{21}
$$

$$
\epsilon^{\text{jka}} E I \mathcal{W}\_{\text{n}}^{\prime\prime}(0) = E I \mathcal{W}\_{\text{n}}^{\prime\prime}(a), \tag{22}
$$

$$e^{i\text{kar}}EI\mathbb{W}\_{n}^{\prime\prime\prime}(0) - e^{i\text{ka}}F\_{n} = EI\mathbb{W}\_{n}^{\prime\prime\prime}(a),\tag{23}$$

where *k* is the Bloch wave vector, also known as the wave number.

To facilitate the application of DQM conveniently, the computational domain of each cell needs to be converted to a standardized computational domain [[−1,1] by the reversible transformation below:

$$\varepsilon = \frac{\mathbf{x} - \mathbf{x}\_n^I}{L\_x} - \mathbf{1},\tag{24}$$

where *n* is the number of the cell; *x<sup>l</sup> <sup>n</sup>* represents the coordinates of the left side of the *n*th beam unit; *Lx* is equal to 0.5*a*, and ε is the local coordinate in the normalized computational domain.

*Crystals* **2019**, *9*, 293

Substituting Equation (24) into Equation (12) and Equations (20)–(23), the governing equation and boundary conditions in the normalized computational domain can be obtained.

$$\frac{EI}{L\_x^4} \frac{\partial^4 W\_\text{H}(\varepsilon)}{\partial \varepsilon^4} - \omega^2 \rho A W\_\text{n}(\varepsilon) = 0,\tag{25}$$

$$
\varepsilon^{\text{ik}\text{ra}}\mathcal{W}\_{\text{n}}(-1) - \mathcal{W}\_{\text{n}}(1) = 0,\tag{26}
$$

$$
\varepsilon^{\text{jka}} \mathcal{W}\_n'(-1) - \mathcal{W}\_n'(1) = 0,\tag{27}
$$

$$e^{i\mathbf{k}\cdot\mathbf{u}}\mathcal{W}\_{\textit{n}}^{\prime\prime}(-1) - \mathcal{W}\_{\textit{n}}^{\prime\prime}(1) = 0,\tag{28}$$

$$e^{i\mathbf{k}\mathbf{z}}\mathcal{W}\_{n}^{\prime\prime\prime}(-1) - \mathcal{W}\_{n}^{\prime\prime\prime}(1) - e^{i\mathbf{k}\mathbf{z}}\frac{L\_{\mathbf{x}}^{3}}{EI} \cdot \frac{m\_{\mathbf{z}}k\_{z}\omega^{2}}{k\_{z} - m\_{z}\omega^{2}}\mathcal{W}\_{n}(-1) = 0,\tag{29}$$

Substituting Equations (1)–(6) into Equations (25)–(29), the boundary conditions and governing equations discretized by DQM can be obtained as below:

$$\frac{EI}{L\_x^4} \sum\_{j=1}^{N\_x} A\_{ij}^{(4)} \mathcal{W}\_j - \omega^2 \rho A \mathcal{W}\_i = 0, \; i = 3, 4, \dots, N\_x - 2,\tag{30}$$

$$e^{i\mathbf{k}\mathbf{x}}\mathcal{W}\_1 - \mathcal{W}\_{\mathbf{N}\mathbf{x}} = \mathbf{0},\tag{31}$$

$$\epsilon^{ika}\sum\_{j=1}^{N\_x} A\_{1j}^{(1)}\mathcal{W}\_j - \sum\_{j=1}^{N\_x} A\_{Nxj}^{(1)}\mathcal{W}\_j = 0,\tag{32}$$

$$e^{ika}\sum\_{j=1}^{N\_x} A\_{1j}^{(2)}\mathcal{W}\_j - \sum\_{j=1}^{N\_x} A\_{Nxj}^{(2)}\mathcal{W}\_j = 0,\tag{33}$$

$$\epsilon^{\mathrm{ik}}\sum\_{j=1}^{N\_{\mathrm{x}}}A\_{1j}^{(3)}\mathcal{W}\_{j} - \sum\_{j=1}^{N\_{\mathrm{x}}}A\_{N\mathrm{x}j}^{(3)}\mathcal{W}\_{j} - \epsilon^{\mathrm{ik}}\frac{L\_{\mathrm{x}}^{3}}{EI} \cdot \frac{m\_{\mathrm{z}}k\_{2}\alpha^{2}}{k\_{z} - m\_{z}\alpha^{2}}\mathcal{W}\_{1} = 0,\tag{34}$$

Equations (30)–(34) can be expressed as a matrix equation form as shown below. Because the third term on the left side of Equation (34) is the non-linear term of ω2, it is impossible to solve the relationship between wave vector *k* and ω by using the conventional matrix-partitioning method. Next, we use the proposed unconventional matrix-partitioning method and the variable substitution method to solve this problem.

$$
\begin{bmatrix}
\mathbf{K}\_{qq}(\omega) & \mathbf{K}\_{qb} & \mathbf{K}\_{qd} \\
\mathbf{K}\_{bq} & \mathbf{K}\_{bb} & \mathbf{K}\_{bd} \\
\mathbf{K}\_{dq} & \mathbf{K}\_{db} & \mathbf{K}\_{dd}
\end{bmatrix} \cdot \begin{pmatrix}
\boldsymbol{\ell}\boldsymbol{I}\_{q} \\
\mathbf{U}\_{b} \\
\mathbf{U}\_{d}
\end{pmatrix} - \omega^{2} \begin{bmatrix}
\mathbf{0} & \mathbf{0} & \mathbf{0} \\
\mathbf{0} & \mathbf{0} & \mathbf{0} \\
\mathbf{0} & \mathbf{0} & \mathbf{M}\_{dd}
\end{bmatrix} \cdot \begin{pmatrix}
\boldsymbol{\ell}\boldsymbol{I}\_{q} \\
\mathbf{U}\_{b} \\
\mathbf{U}\_{d}
\end{pmatrix} = \mathbf{0},\tag{35}
$$

$$\mathbf{M}\_{dd} = \rho A \mathbf{I},\tag{36}$$

where *Uq* = *W*1; **U***<sup>b</sup>* = (*W*2, *WN*−1, *WN*) *<sup>T</sup>*; *Ud* <sup>=</sup> (*W*3, *<sup>W</sup>*4, ... , *WN*−2) *<sup>T</sup>*; **<sup>K</sup>***qb* is a 1 <sup>×</sup> 3 vector; **<sup>K</sup>***qd* is a 1 × (*N* − 4) vector; **K***bq* is a 3 × 1 vector; **K***bb* is a 3 × 3 matrix, which is always invertible; **K***bd* is a 3 × (*N* − 4) matrix; **K***dq* is a (*N* − 4) × 1 vector; **K***db* is a (*N* − 4) × 3 matrix; **K***dd* is a (*N* − 4) × (*N* − 4) matrix; **M***dd* is a (*N* − 4) × (*N* − 4) matrix; K*qq* is a scalar as shown below:

$$K\_{qq} = \frac{c\_1}{p\_z} + c\_{2\prime} \tag{37}$$

where *pz* = *kz* <sup>−</sup> *mz*ω2; *<sup>c</sup>*<sup>1</sup> = <sup>−</sup>*eikak*<sup>2</sup> *zL*3 *<sup>x</sup>*/*EI*; *<sup>c</sup>*<sup>2</sup> <sup>=</sup> *<sup>e</sup>ikaA*(3) 1,1 <sup>−</sup> *<sup>A</sup>*(3) *<sup>N</sup>*,1 <sup>+</sup> *eikakzL*<sup>3</sup> *<sup>x</sup>*/*EI*; Performing a partitioned matrix operation on Equation (35), the following three matrix equations can be obtained:

$$(\frac{c\_1}{p\_z} + c\_2)lI\_q + \mathbf{K}\_{qb} \cdot \mathbf{U}\_b + \mathbf{K}\_{qd} \cdot \mathbf{U}\_d = 0,\tag{38}$$

$$\mathbf{K}\_{h\eta} \mathbf{U}\_{\eta} + \mathbf{K}\_{hb} \cdot \mathbf{U}\_{b} + \mathbf{K}\_{hd} \cdot \mathbf{U}\_{d} = \mathbf{0},\tag{39}$$

$$
\mathbf{K}\_{\rm dq} \mathbf{U}\_q + \mathbf{K}\_{\rm db} \cdot \mathbf{U}\_b + \mathbf{K}\_{\rm dd} \cdot \mathbf{U}\_d - \alpha^2 \rho A \mathbf{I} \cdot \mathbf{U}\_d = 0,\tag{40}
$$

Performing matrix operation and simplification on Equation (39), and using **U***<sup>d</sup>* and *Uq* to represent **U***b*.

$$\mathbf{U}\_b = \mathbf{S}\_1 \cdot \mathbf{U}\_q + \mathbf{S}\_2 \cdot \mathbf{U}\_{d\prime} \tag{41}$$

where **<sup>S</sup>**<sup>1</sup> = <sup>−</sup>**K**−<sup>1</sup> *bb* · **<sup>K</sup>***bq*; **<sup>S</sup>**<sup>2</sup> <sup>=</sup> <sup>−</sup>**K**−<sup>1</sup> *bb* · **K***bd*, substituting Equation (41) into Equation (38) for calculation and simplification, the following equation can be obtained:

$$
\Delta U\_{\eta} = \left(\frac{1}{r\_z} \mathbf{S}\_5 + \mathbf{S}\_6\right) \cdot \mathbf{U}\_{d\prime} \tag{42}
$$

where **S**<sup>5</sup> = (*c*1/*c*4)**S**4; **S**<sup>6</sup> = (−1/*c*4)**S**4; *rz* = *c*4*pz* + *c*1; *c*<sup>4</sup> = *c*<sup>2</sup> + **K***qb* · **S**1; and **S**<sup>4</sup> = **K***qd* + **K***qb* · **S**2. Substituting Equation (42) into Equation (41), one can get:

$$\mathbf{U}\_b = \left(\frac{1}{r\_z}\mathbf{S}\_7 + \mathbf{S}\_8\right) \cdot \mathbf{U}\_{d'} \tag{43}$$

where **S**<sup>7</sup> = **S**<sup>1</sup> · **S**5; and **S**<sup>8</sup> = **S**<sup>1</sup> · **S**<sup>6</sup> + **S**2. Substituting Equation (42) and (43) into Equation (40), one can obtain a standard quadratic eigenvalue equation after simplification.

$$r\_z^2 \mathbf{H}\_2 \cdot \mathbf{U}\_d + r\_z \mathbf{H}\_1 \cdot \mathbf{U}\_d + \mathbf{H}\_0 \cdot \mathbf{U}\_d = 0,\tag{44}$$

$$\mathbf{H}\_f = -\begin{bmatrix} \mathbf{0} & \mathbf{H}\_2 \\ -\mathbf{I} & \mathbf{0} \end{bmatrix}^{-1} \cdot \begin{bmatrix} \mathbf{H}\_0 & \mathbf{H}\_1 \\ \mathbf{0} & \mathbf{I} \end{bmatrix} \tag{45}$$

$$
\omega^2 = (-1/m\_{\overline{z}}c\_4)r\_{\overline{z}} + (k\_{\overline{z}}/m\_{\overline{z}} + c\_1/m\_{\overline{z}}c\_4),
\tag{46}
$$

where **H**<sup>2</sup> = (ρ*A*/*mzc*4)**I**; **H**<sup>1</sup> = **K***dq* · **S**<sup>6</sup> + **K***db* · **S**<sup>8</sup> + **K***dd* − ρ*A*(*kz*/*mz* + *c*1/*mzc*4)**I**; and **H**<sup>0</sup> = **K***dq* · **S**<sup>5</sup> + **K***db* · **S**7. For any given wave vector *k* in the first Brillouin zone, *rz* can be got by calculating the eigenvalue of matrix **H***f*. One can get the circular frequency ω due to Equation (46), thus the dispersion relation and bending vibration band gaps can be plotted.

#### *3.2. Timoshenko Model & Solution Procedures*

For deep beams, the effects of transverse shear deformation and rotary inertia must be considered. Based on the Timoshenko beam theory, the governing equation for the bending vibration of the nth cell is shown below [8]:

$$k\_s G A \left(\frac{\partial \rho\_n(\mathbf{x})}{\partial \mathbf{x}} - \frac{\partial^2 \mathcal{W}\_n(\mathbf{x})}{\partial \mathbf{x}^2}\right) - \alpha^2 \rho A \mathcal{W}\_n(\mathbf{x}) = 0,\tag{47}$$

$$k\_{\rm s}GA\Big(\boldsymbol{\varrho}\_{\rm n}(\mathbf{x}) - \frac{\partial \mathcal{W}\_{\rm n}(\mathbf{x})}{\partial \mathbf{x}}\Big) - EI\frac{\partial^2 \boldsymbol{\varrho}\_{\rm n}(\mathbf{x})}{\partial \mathbf{x}^2} - \omega^2 \rho \boldsymbol{I} \boldsymbol{\varrho}\_{\rm n}(\mathbf{x}) = \mathbf{0},\tag{48}$$

where *ks* is the shear coefficient; *G* is the shear modulus, and ϕ is the rotation of the cross-section. Ignoring the stress concentration between two adjacent units, the following boundary conditions can be listed according to the continuity of displacement, angle of rotation, bending moment, and shear force at the node *xn*+1.

$$\mathcal{W}\_{n+1}(0) = \mathcal{W}\_n(a),\tag{49}$$

$$
\varphi\_{n+1}(0) = \varphi\_n(a),
\tag{50}
$$

$$\begin{cases} EI\varphi'\_{n+1}(0) = EI\varphi'\_n(a), \end{cases} \tag{51}$$

$$k\_s GA \Big[q\_{n+1}(0) - \mathcal{W}'\_{n+1}(0)\Big] - F\_{n+1} = k\_s GA \big[q\_n(a) - \mathcal{W}'\_n(a)\Big],\tag{52}$$

According to the Floquet–Bloch theorem [44], Equations (49)–(52) can be rewritten as follows:

$$
\epsilon^{\rm iku} \mathcal{W}\_n(0) = \mathcal{W}\_n(a),
\tag{53}
$$

$$
\epsilon^{\rm ika} \varphi\_n(0) = \varphi\_n(a),
\tag{54}
$$

$$\epsilon^{\text{ikr}} EI \varphi\_n'(0) = EI \varphi\_n'(a),\tag{55}$$

$$e^{i\mathbf{k}\mathbf{r}}k\_{\mathbf{s}}GA[\varphi\_n(0) - \mathcal{W}\_n'(0)] - e^{i\mathbf{k}\mathbf{r}}F\_n = k\_{\mathbf{s}}GA[\varphi\_n(a) - \mathcal{W}\_n'(a)],\tag{56}$$

Substituting Equation (24) into Equation (47), (48) and Equations (53)–(56), the governing equation and boundary conditions in the normalized computational domain can be obtained.

$$k\_s GA \left(\frac{1}{L\_\chi} \frac{\partial \rho\_n(\varepsilon)}{\partial \varepsilon} - \frac{1}{L\_\chi^2} \frac{\partial^2 \mathcal{W}\_n(\varepsilon)}{\partial \varepsilon^2} \right) - \alpha^2 \rho A \mathcal{W}\_n(\varepsilon) = 0,\tag{57}$$

$$k\_{\varepsilon}GA\Big(\eta\_{\text{fl}}\left(\varepsilon\right)-\frac{1}{L\_{\text{x}}}\frac{\partial \mathcal{W}\_{\text{fl}}\left(\varepsilon\right)}{\partial \varepsilon}\Big)-\frac{EI}{L\_{\text{x}}^{2}}\frac{\partial^{2}\rho\_{\text{fl}}\left(\varepsilon\right)}{\partial \varepsilon^{2}}-\omega^{2}\rho I\rho\_{\text{fl}}\left(\varepsilon\right)=0,\tag{58}$$

$$
\epsilon^{\text{idx}} \mathcal{W}\_n(-1) - \mathcal{W}\_n(1) = 0,\tag{59}
$$

$$
\epsilon^{\text{ik}} q\_{\text{\textquotedblleft}}(-1) - q\_{\text{\textquotedblleft}}(1) = 0,\tag{60}
$$

$$
\epsilon^{\text{ika}} \phi\_n'(-1) - \phi\_n'(1) = 0,\tag{61}
$$

$$\varepsilon^{\bar{j}\bar{k}\bar{n}} \left[ \varphi\_{\bar{n}}(-1) - \frac{1}{L\_{\chi}} \mathsf{W}'\_{n}(-1) \right] - \left[ \varphi\_{\bar{n}}(1) - \frac{1}{L\_{\chi}} \mathsf{W}'\_{n}(1) \right] - \varepsilon^{\bar{j}\bar{k}\bar{n}} \frac{1}{k\_{\text{s}}GA} \frac{m\_{\bar{z}}k\_{z}\alpha^{2}}{k\_{z} - m\_{\bar{z}}\alpha^{2}} \left| \mathsf{W}\_{n}(-1) = 0,\tag{62}$$

Substituting Equations (1)–(6) into Equations (57)–(62), the boundary conditions and governing equations discretized by DQM can be obtained as below:

$$\frac{k\_{\text{S}}G A}{L\_{\text{x}}} \sum\_{j=1}^{N\_{\text{x}}} A\_{ij}^{(1)} \varphi\_{j} - \frac{k\_{\text{S}}G A}{L\_{\text{x}}^{2}} \sum\_{j=1}^{N\_{\text{x}}} A\_{ij}^{(2)} W\_{j} - \omega^{2} \rho A \mathcal{W}\_{i} = 0, \; i = 2, 3, \dots, N\_{\text{x}} - 1,\tag{63}$$

$$k\_5 G A \rho\_i - \frac{k\_5 G A}{L\_x} \sum\_{j=1}^{N\_x} A\_{ij}^{(1)} \mathcal{W}\_j - \frac{E I}{L\_x^2} \sum\_{j=1}^{N\_x} A\_{ij}^{(2)} \rho\_j - \omega^2 \rho I \rho\_i = 0, \; i = 2, 3, \dots, N\_x - 1,\tag{64}$$

$$e^{i\mathbf{k}\mathbf{x}}\mathcal{W}\_1 - \mathcal{W}\_{\mathbf{N}\mathbf{x}} = \mathbf{0},\tag{65}$$

$$x^{\text{ika}}\varphi\_1 - \varphi\_{\text{Nx}} = 0,\tag{66}$$

$$\epsilon^{ikx}\sum\_{j=1}^{N\_{\mathfrak{X}}} A\_{1j}^{(1)}{}\_{j}\boldsymbol{\wp}\_{j} - \sum\_{j=1}^{N\_{\mathfrak{X}}} A\_{Nxj}^{(1)}{}\_{j}\boldsymbol{\wp}\_{j} = \boldsymbol{0},\tag{67}$$

$$\epsilon^{\mathrm{ika}} \Big[ \varphi\_1 - \frac{1}{L\_\chi} \sum\_{j=1}^{N\_\pi} A\_{1j}^{(1)} \mathcal{W}\_j \Big] - \left[ \varphi\_{\mathcal{N}\_\Gamma} - \frac{1}{L\_\chi} \sum\_{j=1}^{N\_\pi} A\_{\mathcal{N}\_\Gamma}^{(1)} \mathcal{W}\_j \right] - \epsilon^{\mathrm{ika}} \frac{1}{k\_\sigma \mathcal{G} A} \frac{m\_\pi k\_\sigma a^2}{k\_\pi - m\_\pi a^2} \mathcal{W}\_1 = 0,\tag{68}$$

Equations (63)–(68) can be expressed as a matrix equation form as shown below. Because the third term on the left side of Equation (68) is the non-linear term of ω2, it is impossible to solve the relationship between wave vector *k* and ω by using the conventional matrix-partitioning method. Next, we use the proposed unconventional matrix-partitioning method and the variable substitution method to solve this problem.

$$
\begin{bmatrix}
\mathbf{K}\_{qq}(\omega) & \mathbf{K}\_{qb} & \mathbf{K}\_{qd} \\
\mathbf{K}\_{bq} & \mathbf{K}\_{bb} & \mathbf{K}\_{bd} \\
\mathbf{K}\_{dq} & \mathbf{K}\_{db} & \mathbf{K}\_{dd}
\end{bmatrix} \cdot \begin{pmatrix}
\boldsymbol{\ell}\boldsymbol{I}\_{q} \\
\mathbf{U}\_{b} \\
\mathbf{U}\_{d}
\end{pmatrix} - \omega^{2} \begin{bmatrix}
\mathbf{0} & \mathbf{0} & \mathbf{0} \\
\mathbf{0} & \mathbf{0} & \mathbf{0} \\
\mathbf{0} & \mathbf{0} & \mathbf{M}\_{dd}
\end{bmatrix} \cdot \begin{pmatrix}
\boldsymbol{\ell}\boldsymbol{I}\_{q} \\
\mathbf{U}\_{b} \\
\mathbf{U}\_{d}
\end{pmatrix} = \mathbf{0},\tag{69}
$$

$$\mathbf{M}\_{dd} = \begin{bmatrix} \rho A \mathbf{I} & \mathbf{0} \\ \mathbf{0} & \rho I \mathbf{I} \end{bmatrix} \tag{70}$$

where *Uq* = *W*1; **U***<sup>b</sup>* = (*WN*,ϕ1,ϕ*N*) *<sup>T</sup>*; *Ud* <sup>=</sup> (*W*2, *<sup>W</sup>*3, ... , *WN*−1,ϕ2,ϕ3, ... ,ϕ*N*−1) *<sup>T</sup>*; **<sup>K</sup>***qb* is a 1 <sup>×</sup> <sup>3</sup> vector; **K***qd* is a 1 × (2*N* − 4) vector; **K***bq* is a 3 × 1 vector; **K***bb* is a 3 × 3 matrix; **K***bd* is a 3 × (2*N* − 4) matrix; **K***dq* is a (2*N* − 4) × 1 vector; **K***db* is a (2*N* − 4) × 3 matrix; **K***dd* is a (2*N* − 4) × (2*N* − 4) matrix; **M***dd* is a (2*N* − 4) × (2*N* − 4) matrix; K*qq* is a scalar as shown below:

$$K\_{qq} = \frac{c\_1}{p\_z} + c\_{2\prime} \tag{71}$$

where *pz* = *kz* <sup>−</sup> *mz*ω2; *<sup>c</sup>*<sup>1</sup> = <sup>−</sup>*eikak*<sup>2</sup> *<sup>z</sup>*/*ksGA*; *<sup>c</sup>*<sup>2</sup> <sup>=</sup> <sup>−</sup>*eikaA*(1) 1,1/*Lx* <sup>−</sup> *<sup>A</sup>*(1) *Nx*,1/*Lx* <sup>+</sup> *eikakz*/*ksGA*; Performing a partitioned matrix operation on Equation (69), the following three matrix equations can be obtained:

$$\left(\frac{c\_1}{p\_z} + c\_2\right) l I\_q + \mathbf{K}\_{qb} \cdot \mathbf{U}\_b + \mathbf{K}\_{qd} \cdot \mathbf{U}\_d = 0,\tag{72}$$

$$\mathbf{K}\_{bq}\mathbf{U}\_q + \mathbf{K}\_{bb} \cdot \mathbf{U}\_b + \mathbf{K}\_{bd} \cdot \mathbf{U}\_d = 0,\tag{73}$$

$$\mathbf{K}\_{\rm dq} \mathbf{U}\_q + \mathbf{K}\_{\rm db} \cdot \mathbf{U}\_b + \mathbf{K}\_{\rm dd} \cdot \mathbf{U}\_d - \omega^2 \mathbf{M}\_{\rm dd} \cdot \mathbf{U}\_d = 0,\tag{74}$$

The following procedure is similar to the Euler model. Substituting Equations (41)–(43) into Equations (72)–(74), one can obtain a standard quadratic eigenvalue equation after simplification.

$$r\_z^2 \mathbf{H}\_2 \cdot \mathbf{U}\_d + r\_z \mathbf{H}\_1 \cdot \mathbf{U}\_d + \mathbf{H}\_0 \cdot \mathbf{U}\_d = 0,\tag{75}$$

For any given wave vector *k* in the first Brillouin zone, *rz* can be got by calculating the eigenvalue of matrix **H***f*. One can get the circular frequency ω from the Equation (46), thus the dispersion relation and bending vibration band gaps can be plotted.

#### **4. Numerical Results and Discussions**

In this section, we first gave a study on the convergence of the proposedmethod. Next, the correctness and accuracy of the present method were verified by comparing the results with the existing literature. It is worth mentioning that, compared with the plane wave expansion method and the finite element method, the present method demonstrated high accuracy and computational efficiency. Finally, the parameter analysis was carried out to discuss the effects of shear deformation and rotary inertia, the lumped mass *mz*, the spring stiffness coefficient *kz*, and the lattice size *a* on the band gap.

Yu et al. [45] used the transfer-matrix method to calculate the bending vibration band gap of a LR beam. For the convenience of comparison and analysis, Yu et al.'s [45] geometry and material parameters of the LR beam are adopted. Figure 2 illustrates a straight beam with locally resonant oscillators. The material of the LR beam is aluminum, and the shape of the cross-section is a ring with outer radius and inner radius *<sup>r</sup>*<sup>1</sup> <sup>=</sup> <sup>1</sup> <sup>×</sup> <sup>10</sup>−<sup>2</sup> m and *<sup>r</sup>*<sup>0</sup> <sup>=</sup> <sup>7</sup> <sup>×</sup> <sup>10</sup>−<sup>3</sup> m, respectively. Each locally resonant oscillator consists of a rubber ring and a copper ring coated on the outside, and their outer radii are *<sup>r</sup>*<sup>2</sup> <sup>=</sup> 1.5 <sup>×</sup> <sup>10</sup>−<sup>2</sup> <sup>m</sup> and *<sup>r</sup>*<sup>3</sup> <sup>=</sup> 1.95 <sup>×</sup> <sup>10</sup>−<sup>2</sup> m, respectively. Both rings have the same width *lz* <sup>=</sup> <sup>1</sup> <sup>×</sup> <sup>10</sup>−<sup>2</sup> m, and the lattice constant of the LR beam is taken as *a* = 7.5 <sup>×</sup> 10−<sup>2</sup> m. Material parameters of the LR beam are listed in Table 1. Unless otherwise specified, values of the parameters in the following studies are consistent with those here.

**Figure 2.** (**a**) Illustration of a straight beam with locally resonant oscillators. (**b**) The sketch of the locally resonant oscillator.


**Table 1.** Material parameters.

For a ring rubber, the radial equivalent stiffness can be expressed as follows:

$$k\_z = \frac{\pi (3.29H\_z^2 + 5) G\_{rubber} l\_z}{\ln(r\_2/r\_1)}\tag{76}$$

where *Hz* = 1/(*r*<sup>1</sup> + *r*2)*ln*(*r*2/*r*1) is the shape coefficient.

#### *4.1. Convergence Study*

To study the convergence of the present method proposed in this work, Table 2 shows the fundamental frequencies of the locally resonant beam. We give a series of results corresponding to various numbers of discrete points *Nx*. It is found that when the number of sampling points *Nx* ≥ 8, the frequency converges rapidly. In the rest part of this paper, 15 discrete points are selected to calculate and analyze the bending vibration band gap characteristics of LR beams in order to ensure good accuracy.

**Table 2.** Fundamental frequency of a locally resonant (LR) beam.


#### *4.2. Validation*

As Figure 3 shows, the curves and scatters represent the dispersion relationship of a straight beam with LR oscillator structures. The dispersion relationship curve does not cover the full frequency range and the frequency range through which no dispersion curve covers are represented by a shadow zone, which is the band gap. The bending wave in the band gap frequency range cannot propagate through the beam. The 1st band gap locates between 308.722–478.943 *Hz*, and the range obtained by Yu et al. [45] is 309.1–479.4 *Hz*. From a qualitative point of view, the scatter points are all on the curve. From a quantitative point of view, the relative error between the present results and Yu et al.'s is within 0.13%. The above proves that our numerical solutions match well with those in the existing literature.

**Figure 3.** Dispersion relationship and band gaps for a locally resonant (LR) beam.

Consider a finite locally resonant beam consisting of eight preceding periodic cells, in which a harmonic displacement excitation in the y-direction between 0 and 800 Hz is applied at one end and the frequency response function (FRF) at the other end is shown in Figure 4a. The solid line represents the frequency response function and the dashed line represents the input spectrum (0 dB). Within the frequency range marked by a double arrow, the FRF has a maximum attenuation of more than 60 dB. For infinite structures, the imaginary part of k causes attenuating vibration. The larger the absolute value of the imaginary part, the stronger the spatial attenuation of the evanescent wave. As shown in Figure 4b, the band gap frequency range of the infinite structure is in good agreement with that of the finite structure, which also proves the correctness of the present method from the perspective of the evanescent modes.

**Figure 4.** Transmission properties of the finite locally resonant (LR) beam and band gap characteristics of the corresponding infinite beam in the range of 0–800 Hz: (**a**) Frequency response function of the finite beam. (**b**) The imaginary part of wave vectors of the infinite beam.

#### *4.3. Method Advantages*

To demonstrate the advantages of the present approach, the following case is compared with the PWE method and the FEM method. Han et al. [11] used the modified transfer matrix method (MTM) and the PWE method to calculate the band gaps of a PC beam. The frequency ranges of the first three band gaps obtained by MTM, PWE, FEM, and the present method are listed in Table 3. As with the TM method, the MTM obtains an analytical solution; therefore, the closer the results of the other three numerical methods to the results obtained by the MTM, the higher the accuracy. The relative errors between FEM and DQM are within 0.014%, and both methods show higher accuracy than the PWE method. The following compares the methods from the perspective of computational efficiency. A total of 40 wave vectors are selected at equal intervals in the first Brillouin zone, and the same computer and software are used for the calculation. Both the number of the DQM discrete points and the number of the FEM nodes are taken as 30. The relevant simulation times and computational resources (memory) are listed in Table 4. It can be seen that the present method has a shorter simulation time and requires less memory than the FEM method. Therefore, by comparing with two classical, widely-used methods, the high accuracy and computational efficiency of the present method are demonstrated.



**Table 4.** Simulation times and memory required of the present method and finite element method (FEM) (Intel (R) Xeon(R) E5620 CPU, Mathematica 11.1, simulation time is in second, memory is in KB).


#### *4.4. Parameter Studies*

In this part, we conducted a parametric analysis to investigate the effect of shear deformation and rotary inertia, the lumped mass *mz*, the spring stiffness coefficient *kz*, and the lattice size *a* on the 1st band gap of the locally resonant beam. The emphasis is on changes in the lower edge as well as the width of the 1st band gap. Since most of the vibrations presented in the project are low-frequency vibrations, for engineering purposes, the following research focuses on ways to decrease the corresponding frequency and widen the band gap.

#### 4.4.1. Effects of Shear Deformation and Rotary Inertia

In this case, the material and geometric parameters of a LR beam are as follows: *<sup>E</sup>* <sup>=</sup> 4.35 <sup>×</sup> <sup>10</sup><sup>9</sup> Pa, <sup>ρ</sup> <sup>=</sup> 1180 kg/m3, *<sup>G</sup>* <sup>=</sup> <sup>5</sup> <sup>×</sup> <sup>10</sup><sup>8</sup> Pa, *<sup>A</sup>* <sup>=</sup> <sup>1</sup> <sup>×</sup> <sup>10</sup>−<sup>4</sup> <sup>m</sup>2, *<sup>I</sup>* <sup>=</sup> 8.333 <sup>×</sup> <sup>10</sup>−<sup>10</sup> <sup>m</sup>4, *ks* <sup>=</sup> 0.8333, *<sup>a</sup>* <sup>=</sup> 0.075 m. Both the Euler-Bernoulli beam model (EB) and the Timoshenko beam model (TB) are used to calculate the band structure of the LR beam. The results of the two models are plotted in Figure 5. The 1st band gap locates between 260.301–743.608 *Hz* by EB, and 253.44–734.597 *Hz* by TB, respectively. The shear deformation will reduce the stiffness of the beam, and the rotary inertia will increase the inertia of the beam, both of which will reduce the natural frequency of the beam. Therefore, the band gap frequency range calculated by TB is lower than that calculated by EB.

**Figure 5.** Dispersion relationship and band gaps for a locally resonant (LR) beam calculated by Timoshenko beam model (TB) and Euler-Bernoulli beam model (EB). Half of the dispersion relationship is plotted because of its symmetry.

4.4.2. Effects of Lumped Mass *mz*

The mass ratio *mz*/*m* is introduced here to characterize the relative magnitude of the lumped mass, where *m* = ρ*Aa* is the mass of the LR beam per period. Figure 6 manifests the changes in the lower edge and width of the 1st band gap over the mass ratio *mz*/*m* from 0 to 500. As is demonstrated in Figure 6a that the lower edge falls considerably when the mass ratio (*mz*/*m*) is small. Since then, it remains more or less stable and finally tends to 0 Hz. According to Figure 6b, the width of the 1st band gap grows rapidly when the mass ratio (*mz*/*m*) is small. After that, it reaches a plateau and tends to be stable. According to the above, we have found that it may be possible to simultaneously decrease the corresponding frequency range and widen the band gap by increasing the lumped mass *mz* within a certain range.

**Figure 6.** Effects of lumped mass *mz* on the 1st band gap: (**a**) Lower edge of the 1st band gap. (**b**) Width of the 1st band gap.

4.4.3. Effects of the Spring Stiffness Coefficient *kz*

In order to facilitate the description of the relative magnitude of the spring stiffness coefficient *kz*, a stiffness ratio *kz*/*k* is introduced here, where *k* = *EI*/*a* is the equivalent stiffness of the LR beam in a period. Figure 7 illustrates the changes in the lower edge and width of the 1st band gap over the stiffness ratio *kz*/*k* from 0 to 250. It is apparent from Figure 7 that both the width and the lower edge experience a steady growth with the stiffness ratio *kz*/*k* increasing. Therefore, we can conclude that it is hard to simultaneously decrease the corresponding frequency range and widen the band gap by merely changing the stiffness coefficient of the LR oscillator.

**Figure 7.** Effects of the spring stiffness coefficient *kz* on the 1st band gap: (**a**) Lower edge of the 1st band gap. (**b**) Width of the 1st band gap.

#### 4.4.4. Effects of the Lattice Constant *a*

˖

When propagating in the locally resonant beam, the elastic wave is reflected at the nodes, which satisfies the Bragg band gap mechanism. Thus, there are also Bragg band gaps in the LR beam. In the case of small lattice constant, the LR band gap is separated from the Bragg band gap frequency range. Here, we mainly focus on the effects on the 1st LR band gap, so the lattice size *a* selected here is between (0, 0.075], which can separate the LR band gap from the Bragg band gap. As is shown in Figure 8, with the lattice constant *a* increasing, the band gap width decreases gradually, showing an inverse correlation (Figure 8b), but the lower edge always levels off at about 309 Hz (Figure 8a). It shows that the starting frequency of the locally resonant band gap is independent of the lattice constant because its frequency is determined by *<sup>f</sup>*<sup>0</sup> <sup>=</sup> <sup>√</sup> *kz*/*mz*/(2π) [46], which is the resonant frequency of the locally resonant oscillator. Altering the lattice constant can only affect the band gap width.

**Figure 8.** Effects of the lattice constant *a* on the 1st band gap: (**a**) Lower edge of the 1st band gap. (**b**) Width of the 1st band gap.

#### **5. Conclusions**

The differential quadrature method has been developed to calculate the elastic band gaps from the Bragg reflection mechanism in periodic structures efficiently and accurately. However, there have been no reports on the successful use of this method to calculate the band gaps of locally resonant structures. Hence, this paper proposes a numerical method to calculate and study the flexural vibration band gap of a locally resonant beam. The proposed method is based on the DQM method, an unconventional matrix-partitioning method, and a variable substitution method. According to the analysis and research in this work, the following conclusions can be obtained.


Referring to the method and solving process in this paper, future research can apply the proposed method to study the band gap characteristics of 2-D or 3-D locally resonant structures.

**Author Contributions:** Conceptualization, X.L. and T.W.; methodology, X.L. and Y.D.; software, T.W. and Z.L.; validation, Y.R.; formal analysis, T.W.; writing—original draft preparation, T.W.; writing—review and editing, T.W.; visualization, X.J.; supervision, X.L. and X.J.; project administration, X.L.; funding acquisition, X.L.

**Acknowledgments:** This work was supported by the National Natural Science Foundation of China (Grant No. 51879231, 51679214, 51338009, 51409228) and the Zhejiang Provincial Key Research and Development Program (2018C03031).

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

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
