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

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.

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,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,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,21,22,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,39,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:

(1)

Based on the differential quadrature method, we propose a numerical method for calculating the flexural vibration band gaps of a locally resonant beam. In this paper, the differential quadrature method is applied to calculate the band gaps of locally resonant structures for the first time.

(2)

By using an unconventional matrix-partitioning method and a variable substitution method, we can transform the problem of solving the dispersion relation into a standard quadratic eigenvalue problem, and the problem of the non-linear term after using the DQM method can be solved easily.

(3)

Compared with the plane wave expansion method and the finite element method, the high accuracy and computational efficiency of the present method are demonstrated.

(4)

The proposed method has high precision and a rapid speed of convergence.

2. Method

2.1. Differential 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(\epsilon )={\displaystyle \sum _{j=1}^{N_x}{p}_j(\epsilon )}{f}_j,$$

(1)

$$\frac{{\partial }^{r}f({\epsilon }_i)}{\partial x^r}={\displaystyle \sum _{i=1}^{N_x}{A}_{ij}^{(r)}\cdot f({\epsilon }_i)},$$

(2)

where $N_x$ is the total number of discrete points in the calculation domain; $i,j=1,2,\dots ,N_x$; $f({\epsilon }_j)=f_j$; $p_j(\epsilon )$ are the Lagrange interpolation polynomials, and $A_{ij}^r$ is the weight coefficient of the r^th derivative defined by [41]:

$$A_{ij}^{(1)}=\frac{{\displaystyle \prod _{r=1,r\ne i}^{N_x}({\epsilon }_i-{\epsilon }_r)}}{({\epsilon }_i-{\epsilon }_j){\displaystyle \prod _{r=1,r\ne j}^{N_x}({\epsilon }_j-{\epsilon }_r)}},$$

(3)

$$A_{ij}^{(r)}=r\left(A_{ij}^{(r-1)}{A}_{ij}^1-\frac{A_{ij}^{(r-1)}}{x_i-x_j}\right),$$

(4)

where $r$ = 1, 2, …, N_x−1; i,j = 1, 2, …, N_x, (i ≠ j); and $A_{ii}^{(r)}$ are defined as:

$$A_{ii}^{(r)}=-{\displaystyle \sum _{j=1,j\ne i}^{N_x}{A}_{ij}^{(r)}},$$

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

$${\epsilon }_i=-\cos \left(\frac{\pi (i-1)}{N_x-1}\right), i=1,2,\dots ,N_x,$$

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

$$\left[\begin{array}{ccc}{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{array}\right]\cdot \left(\begin{array}{c}{U}_q\\ {\mathbf{U}}_b\\ {\mathbf{U}}_d\end{array}\right)-{\omega }^2\left[\begin{array}{ccc}0& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& {\mathbf{M}}_{dd}\end{array}\right]\cdot \left(\begin{array}{c}{U}_q\\ {\mathbf{U}}_b\\ {\mathbf{U}}_d\end{array}\right)=\mathbf{0},$$

(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,$$

(8)

$${\mathbf{H}}_f=-{\left[\begin{array}{cc}\mathbf{0}& {\mathbf{H}}_2\\ -\mathbf{I}& \mathbf{0}\end{array}\right]}^{-1}\cdot \left[\begin{array}{cc}{\mathbf{H}}_0& {\mathbf{H}}_1\\ \mathbf{0}& \mathbf{I}\end{array}\right],$$

(9)

$${\omega }^2=\left(-1/m_z{c}_4\right)r_z+(k_z/m_z+c_1/m_z{c}_4),$$

(10)

For any given wave vector k in the first Brillouin zone, r_z can be obtained by calculating the eigenvalue of matrix H_f. 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 m_z and a spring with stiffness k_z, 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.

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(x,t)}{\partial x^4}+\rho A\frac{{\partial }^2{w}_n(x,t)}{\partial t^2}=0,$$

(11)

where $\rho$ 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\omega t)$, where $W(x)$ is the vibration amplitude of the beam at $x$, and $\omega$ is the circular frequency, Equation (11) can be rewritten as follows:

$$EI\frac{{\partial }^4{W}_n(x)}{\partial x^4}-{\omega }^2\rho A{W}_n(x)=0,$$

(12)

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

$$f_{n+1}(t)-m_z{\ddot{Z}}_{n+1}(t)=0,$$

(13)

where $f_{n+1}(t)$ is the interaction force between the beam and the oscillator at node $x_{n+1}$, $Z_{n+1}(t)=V_{n+1}\exp (-i\omega 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 $V_{n+1}$.

According to Hooke's law, $f_{n+1}(t)$ can be expressed as follows:

$$f_{n+1}(t)=k_z[w(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),$$

(14)

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

$$V_{n+1}=\frac{k_z}{k_z-m_z{\omega }^2}{W}_{n+1}(0),$$

(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 $x_{n+1}$.

$$W_{n+1}(0)=W_n(a),$$

(16)

$$W_{n+1}^{\prime }(0)=W_n^{\prime }(a),$$

(17)

$$EI{W}_{n+1}^{″}(0)=EI{W}_n^{″}(a),$$

(18)

$$EI{W}_{n+1}^{‴}(0)-F_{n+1}=EI{W}_n^{‴}(a),$$

(19)

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

$$e^{ika}{W}_n(0)=W_n(a),$$

(20)

$$e^{ika}{W}_n^{\prime }(0)=W_n^{\prime }(a),$$

(21)

$$e^{ika}EI{W}_n^{″}(0)=EI{W}_n^{″}(a),$$

(22)

$$e^{ika}EI{W}_n^{‴}(0)-e^{ika}{F}_n=EI{W}_n^{‴}(a),$$

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

$$\epsilon =\frac{x-x_n^l}{L_x}-1,$$

(24)

where $n$ is the number of the cell; $x_n^l$ represents the coordinates of the left side of the n^th beam unit; $L_x$ is equal to $0.5a$, and $\epsilon$ is the local coordinate in the normalized computational domain.

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}_n(\epsilon )}{\partial {\epsilon }^4}-{\omega }^2\rho A{W}_n(\epsilon )=0,$$

(25)

$$e^{ika}{W}_n(-1)-W_n(1)=0,$$

(26)

$$e^{ika}{W}_n^{\prime }(-1)-W_n^{\prime }(1)=0,$$

(27)

$$e^{ika}{W}_n^{″}(-1)-W_n^{″}(1)=0,$$

(28)

$$e^{ika}{W}_n^{‴}(-1)-W_n^{‴}(1)-e^{ika}\frac{L_x^3}{EI}\cdot \frac{m_z{k}_z{\omega }^2}{k_z-m_z{\omega }^2}{W}_n(-1)=0,$$

(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}{\displaystyle \sum _{j=1}^{N_x}{A}_{ij}^{(4)}{W}_j}-{\omega }^2\rho A{W}_i=0, i=3,4,\dots ,N_x-2,$$

(30)

$$e^{ika}{W}_1-W_{Nx}=0,$$

(31)

$$e^{ika}{\displaystyle \sum _{j=1}^{N_x}{A}_{1j}^{(1)}{W}_j}-{\displaystyle \sum _{j=1}^{N_x}{A}_{Nxj}^{(1)}{W}_j}=0,$$

(32)

$$e^{ika}{\displaystyle \sum _{j=1}^{N_x}{A}_{1j}^{(2)}{W}_j}-{\displaystyle \sum _{j=1}^{N_x}{A}_{Nxj}^{(2)}{W}_j}=0,$$

(33)

$$e^{ika}{\displaystyle \sum _{j=1}^{N_x}{A}_{1j}^{(3)}{W}_j}-{\displaystyle \sum _{j=1}^{N_x}{A}_{Nxj}^{(3)}{W}_j}-e^{ika}\frac{L_x^3}{EI}\cdot \frac{m_z{k}_z{\omega }^2}{k_z-m_z{\omega }^2}{W}_1=0,$$

(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 ω², 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.

$$\left[\begin{array}{ccc}{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{array}\right]\cdot \left(\begin{array}{c}{U}_q\\ {\mathbf{U}}_b\\ {\mathbf{U}}_d\end{array}\right)-{\omega }^2\left[\begin{array}{ccc}0& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& {\mathbf{M}}_{dd}\end{array}\right]\cdot \left(\begin{array}{c}{U}_q\\ {\mathbf{U}}_b\\ {\mathbf{U}}_d\end{array}\right)=\mathbf{0},$$

(35)

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

(36)

where $U_q=W_1$; ${\mathbf{U}}_b={\left(W_2,W_{N-1},W_N\right)}^T$; $U_d={\left(W_3,W_4,\dots ,W_{N-2}\right)}^T$; K_qb is a 1 × 3 vector; K_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,$$

(37)

where $p_z=k_z-m_z{\omega }^2$; $c_1=-e^{ika}{k}_z^2{L}_x^3/EI$; $c_2=e^{ika}{A}_{1,1}^{(3)}-A_{N,1}^{(3)}+e^{ika}{k}_z{L}_x^3/EI$; Performing a partitioned matrix operation on Equation (35), the following three matrix equations can be obtained:

$$\left(\frac{c_1}{p_z}+c_2\right)U_q+{\mathbf{K}}_{qb}\cdot {\mathbf{U}}_b+{\mathbf{K}}_{qd}\cdot {\mathbf{U}}_d=0,$$

(38)

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

(39)

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

(40)

Performing matrix operation and simplification on Equation (39), and using U_d and U_q to represent U_b.

$${\mathbf{U}}_b={\mathbf{S}}_1\cdot U_q+{\mathbf{S}}_2\cdot {\mathbf{U}}_d,$$

(41)

where ${\mathbf{S}}_1=-{\mathbf{K}}_{bb}^{-1}\cdot {\mathbf{K}}_{bq}$; ${\mathbf{S}}_2=-{\mathbf{K}}_{bb}^{-1}\cdot {\mathbf{K}}_{bd}$, substituting Equation (41) into Equation (38) for calculation and simplification, the following equation can be obtained:

$$U_q=\left(\frac{1}{r_z}{\mathbf{S}}_5+{\mathbf{S}}_6\right)\cdot {\mathbf{U}}_d,$$

(42)

where ${\mathbf{S}}_5=\left(c_1/c_4\right){\mathbf{S}}_4$; ${\mathbf{S}}_6=\left(-1/c_4\right){\mathbf{S}}_4$; $r_z=c_4{p}_z+c_1$; $c_4=c_2+{\mathbf{K}}_{qb}\cdot {\mathbf{S}}_1$; and ${\mathbf{S}}_4={\mathbf{K}}_{qd}+{\mathbf{K}}_{qb}\cdot {\mathbf{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,$$

(43)

where ${\mathbf{S}}_7={\mathbf{S}}_1\cdot {\mathbf{S}}_5$; and ${\mathbf{S}}_8={\mathbf{S}}_1\cdot {\mathbf{S}}_6+{\mathbf{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,$$

(44)

$${\mathbf{H}}_f=-{\left[\begin{array}{cc}\mathbf{0}& {\mathbf{H}}_2\\ -\mathbf{I}& \mathbf{0}\end{array}\right]}^{-1}\cdot \left[\begin{array}{cc}{\mathbf{H}}_0& {\mathbf{H}}_1\\ \mathbf{0}& \mathbf{I}\end{array}\right],$$

(45)

$${\omega }^2=\left(-1/m_z{c}_4\right)r_z+(k_z/m_z+c_1/m_z{c}_4),$$

(46)

where ${\mathbf{H}}_2=\left(\rho A/m_z{c}_4\right)\mathbf{I}$; ${\mathbf{H}}_1={\mathbf{K}}_{dq}\cdot {\mathbf{S}}_6+{\mathbf{K}}_{db}\cdot {\mathbf{S}}_8+{\mathbf{K}}_{dd}-\rho A\left(k_z/m_z+c_1/m_z{c}_4\right)\mathbf{I}$; and ${\mathbf{H}}_0={\mathbf{K}}_{dq}\cdot {\mathbf{S}}_5+{\mathbf{K}}_{db}\cdot {\mathbf{S}}_7$. For any given wave vector k in the first Brillouin zone, r_z 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 n^th cell is shown below [8]:

$$k_{s}GA\left(\frac{\partial {\phi }_n(x)}{\partial x}-\frac{{\partial }^2{W}_n(x)}{\partial x^2}\right)-{\omega }^2\rho A{W}_n(x)=0,$$

(47)

$$k_{s}GA\left({\phi }_n(x)-\frac{\partial W_n(x)}{\partial x}\right)-EI\frac{{\partial }^2{\phi }_n(x)}{\partial x^2}-{\omega }^2\rho I{\phi }_n(x)=0,$$

(48)

where $k_s$ is the shear coefficient; $G$ is the shear modulus, and $\phi$ 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 $x_{n+1}$.

$$W_{n+1}(0)=W_n(a),$$

(49)

$${\phi }_{n+1}(0)={\phi }_n(a),$$

(50)

$$EI{\phi }_{n+1}^{\prime }(0)=EI{\phi }_n^{\prime }(a),$$

(51)

$$k_{s}GA\left[{\phi }_{n+1}(0)-W_{n+1}^{\prime }(0)\right]-F_{n+1}=k_{s}GA\left[{\phi }_n(a)-W_n^{\prime }(a)\right],$$

(52)

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

$$e^{ika}{W}_n(0)=W_n(a),$$

(53)

$$e^{ika}{\phi }_n(0)={\phi }_n(a),$$

(54)

$$e^{ika}EI{\phi }_n^{\prime }(0)=EI{\phi }_n^{\prime }(a),$$

(55)

$$e^{ika}{k}_{s}GA\left[{\phi }_n(0)-W_n^{\prime }(0)\right]-e^{ika}{F}_n=k_{s}GA\left[{\phi }_n(a)-W_n^{\prime }(a)\right],$$

(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_x}\frac{\partial {\phi }_n(\epsilon )}{\partial \epsilon }-\frac{1}{L_x^2}\frac{{\partial }^2{W}_n(\epsilon )}{\partial {\epsilon }^2}\right)-{\omega }^2\rho A{W}_n(\epsilon )=0,$$

(57)

$$k_{s}GA\left({\phi }_n(\epsilon )-\frac{1}{L_x}\frac{\partial W_n(\epsilon )}{\partial \epsilon }\right)-\frac{EI}{L_x^2}\frac{{\partial }^2{\phi }_n(\epsilon )}{\partial {\epsilon }^2}-{\omega }^2\rho I{\phi }_n(\epsilon )=0,$$

(58)

$$e^{ika}{W}_n(-1)-W_n(1)=0,$$

(59)

$$e^{ika}{\phi }_n(-1)-{\phi }_n(1)=0,$$

(60)

$$e^{ika}{\phi }_n^{\prime }(-1)-{\phi }_n^{\prime }(1)=0,$$

(61)

$$e^{ika}\left[{\phi }_n(-1)-\frac{1}{L_x}{W}_n^{\prime }(-1)\right]-\left[{\phi }_n(1)-\frac{1}{L_x}{W}_n^{\prime }(1)\right]-e^{ika}\frac{1}{k_{s}GA}\frac{m_z{k}_z{\omega }^2}{k_z-m_z{\omega }^2}{W}_n(-1)=0,$$

(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_{s}GA}{L_x}{\displaystyle \sum _{j=1}^{N_x}{A}_{ij}^{(1)}{\phi }_j-}\frac{k_{s}GA}{L_x^2}{\displaystyle \sum _{j=1}^{N_x}{A}_{ij}^{(2)}{W}_j}-{\omega }^2\rho A{W}_i=0, i=2,3,\dots ,N_x-1,$$

(63)

$$k_{s}GA{\phi }_i-\frac{k_{s}GA}{L_x}{\displaystyle \sum _{j=1}^{N_x}{A}_{ij}^{(1)}{W}_j-}\frac{EI}{L_x^2}{\displaystyle \sum _{j=1}^{N_x}{A}_{ij}^{(2)}{\phi }_j}-{\omega }^2\rho I{\phi }_i=0, i=2,3,\dots ,N_x-1,$$

(64)

$$e^{ika}{W}_1-W_{Nx}=0,$$

(65)

$$e^{ika}{\phi }_1-{\phi }_{Nx}=0,$$

(66)

$$e^{ika}{\displaystyle \sum _{j=1}^{N_x}{A}_{1j}^{(1)}{\phi }_j}-{\displaystyle \sum _{j=1}^{N_x}{A}_{Nxj}^{(1)}{\phi }_j}=0,$$

(67)

$$e^{ika}\left[{\phi }_1-\frac{1}{L_x}{\displaystyle \sum _{j=1}^{N_x}{A}_{1j}^{(1)}}{W}_j\right]-\left[{\phi }_{N_x}-\frac{1}{L_x}{\displaystyle \sum _{j=1}^{N_x}{A}_{N_{x}j}^{(1)}}{W}_j\right]-e^{ika}\frac{1}{k_{s}GA}\frac{m_z{k}_z{\omega }^2}{k_z-m_z{\omega }^2}{W}_1=0,$$

(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 ω², 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.

$$\left[\begin{array}{ccc}{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{array}\right]\cdot \left(\begin{array}{c}{U}_q\\ {\mathbf{U}}_b\\ {\mathbf{U}}_d\end{array}\right)-{\omega }^2\left[\begin{array}{ccc}0& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& {\mathbf{M}}_{dd}\end{array}\right]\cdot \left(\begin{array}{c}{U}_q\\ {\mathbf{U}}_b\\ {\mathbf{U}}_d\end{array}\right)=\mathbf{0},$$

(69)

$${\mathbf{M}}_{dd}=\left[\begin{array}{cc}\rho A\mathbf{I}& \mathbf{0}\\ \mathbf{0}& \rho I\mathbf{I}\end{array}\right],$$

(70)

where $U_q=W_1$; ${\mathbf{U}}_b={\left(W_N,{\phi }_1,{\phi }_N\right)}^T$; $U_d={\left(W_2,W_3,\dots ,W_{N-1},{\phi }_2,{\phi }_3,\dots ,{\phi }_{N-1}\right)}^T$; K_qb is a 1 × 3 vector; K_qd is a 1 × (2N − 4) vector; K_bq is a 3 × 1 vector; K_bb is a 3 × 3 matrix; K_bd is a 3 × (2N − 4) matrix; K_dq is a (2N − 4) × 1 vector; K_db is a (2N − 4) × 3 matrix; K_dd is a (2N − 4) × (2N − 4) matrix; M_dd is a (2N − 4) × (2N − 4) matrix; K_qq is a scalar as shown below:

$$K_{qq}=\frac{c_1}{p_z}+c_2,$$

(71)

where $p_z=k_z-m_z{\omega }^2$; $c_1=-e^{ika}{k}_z^2/k_{s}GA$; $c_2=-e^{ika}{A}_{1,1}^{(1)}/L_x-A_{N_x,1}^{(1)}/L_x+e^{ika}{k}_z/k_{s}GA$; 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)U_q+{\mathbf{K}}_{qb}\cdot {\mathbf{U}}_b+{\mathbf{K}}_{qd}\cdot {\mathbf{U}}_d=0,$$

(72)

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

(73)

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

(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,$$

(75)

For any given wave vector k in the first Brillouin zone, r_z 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 proposed method. 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 m_z, the spring stiffness coefficient k_z, 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 r₁ = 1 × 10⁻² m and r₀ = 7 × 10⁻³ m, respectively. Each locally resonant oscillator consists of a rubber ring and a copper ring coated on the outside, and their outer radii are r₂ = 1.5 × 10⁻² m and r₃ = 1.95 × 10⁻² m, respectively. Both rings have the same width l_z = 1 × 10⁻² m, and the lattice constant of the LR beam is taken as a = 7.5 × 10⁻² m. Material parameters of the LR beam are listed in

Table 1. Material parameters.

Parameters	Items	Values
ρAl	Aluminum density (kg/m3)	2600
EAl	Elastic modulus of aluminum (Pa)	7 × 1010
ρrubber	Rubber density (kg/m3)	1300
Erubber	Elastic modulus of rubber (Pa)	7.7 × 105
Grubber	Shear modulus of rubber (Pa)	2.6 × 105
ρCu	Copper density (kg/m3)	8950

. Unless otherwise specified, values of the parameters in the following studies are consistent with those here.

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)},$$

(76)

where H_z = 1/(r₁ + r₂)ln(r₂/r₁) is the shape coefficient.

4.1. Convergence Study

To study the convergence of the present method proposed in this work,

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

Wave Vector k	Present Frequency (Hz)
	Nx = 6	Nx = 7	Nx = 8	Nx = 9	Nx = 15
0.0	0.000	0.000	0.000	0.000	0.000
0.3	277.547	277.635	277.638	277.638	277.638
0.5	304.968	304.972	304.973	304.973	304.973
0.7	308.056	308.057	308.057	308.057	308.057
1.0	308.722	308.722	308.722	308.722	308.722

shows the fundamental frequencies of the locally resonant beam. We give a series of results corresponding to various numbers of discrete points N_x. It is found that when the number of sampling points N_x ≥ 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.

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 H_z, and the range obtained by Yu et al. [45] is 309.1–479.4 H_z. 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.

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.

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. The first three band gaps in a phononic crystals (PC) beam.

Method	Band Gaps (Hz)
	First		Second		Third
	Lower	Upper	Lower	Upper	Lower	Upper
Present	616.378	1097.33	3038.92	6334.27	9601.44	11929.7
FEM	616.379	1097.33	3038.94	6334.57	9602.36	11931.3
MTM	616	1098	3038	6335	9601	11930
PWE	617	1103	3053	6358	9662	11931

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

Method	Present	FEM
Simulation time	8.8438	17.0781
Memory	154584	160380

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

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 m_z, the spring stiffness coefficient k_z, 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: E = 4.35 × 10⁹ Pa, ρ = 1180 kg/m³, G = 5 × 10⁸ Pa, A = 1 × 10⁻⁴ m², I = 8.333 × 10⁻¹⁰ m⁴, k_s = 0.8333, a = 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 H_z by EB, and 253.44–734.597 H_z 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.

4.4.2. Effects of Lumped Mass m_z

The mass ratio m_z/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 m_z/m from 0 to 500. As is demonstrated in Figure 6a that the lower edge falls considerably when the mass ratio (m_z/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 (m_z/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 m_z within a certain range.

4.4.3. Effects of the Spring Stiffness Coefficient k_z

In order to facilitate the description of the relative magnitude of the spring stiffness coefficient k_z, a stiffness ratio k_z/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 k_z/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 k_z/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.

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 $f_0=\sqrt{k_z/m_z}/(2\pi )$ [46], which is the resonant frequency of the locally resonant oscillator. Altering the lattice constant can only affect the band gap width.

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.

(1)

The governing equation and periodic boundary conditions are discretized by the DQM method. The matrix equation is transformed into a standard quadratic eigenvalue problem by the partitioned matrix operation and variable substitution. Thus, the dispersion relationship and the band gap characteristics of a locally resonant beam can be solved. By extending the proposed method further, it can also be suitable for 2-dimensional or 3-dimensional LR structures.

(2)

By comparing with the results of the existing literature, the validity of the proposed method is developed from both propagation modes and evanescent modes. Convergence studies indicate that accurate enough results could be got when the number of discrete points N_x ≥ 8.

(3)

By comparing with the plane wave expansion method and the finite element method, the high accuracy and computational efficiency of the present method are demonstrated.

(4)

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. In addition, one can decrease the lower edge of the 1st band gap by increasing the mass ratio or decreasing the stiffness ratio. The starting frequency of the LR band gap has nothing to do with the lattice constant because it is defined by the resonant frequency of the LR oscillator. The band gap frequency range calculated by the Timoshenko beam theory is lower than that calculated by the Euler–Bernoulli beam theory.

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.