Analysis of Bending Waves in Phononic Crystal Beams with Defects

Existing investigations on imperfect phononic crystal beams mainly concern periodic multi-span beams carrying either one or two channel waves with random or deterministic disorder in span-length. This paper studies the two channel bending waves in phononic crystal beams consisting of many phases of materials with defects introduced as one structural segment having different cross-sectional dimensions or material parameters. The method of reverberation-ray matrix (MRRM) based on the Timoshenko beam theory, which can conduct high-frequency analysis, is extended for the theoretical analysis of dispersion and transmission of bending waves. The supercell technique and the Floquet–Bloch theorem are adopted for modeling the dispersion characteristics, and the whole finite structural model is used to calculate the transmission spectra. Experimental measurements and numerical calculations are provided to validate the displacement transmission obtained by the proposed MRRM, with the effect of damping on transmission spectra being concerned. The high-frequency calculation applicability of the proposed MRRM is also confirmed by comparing the present results with the corresponding ones either using the transfer matrix method (TMM) or MRRM based on Euler—Bernoulli beam theory. The influences of defect size, defect form, and unit-cell number on the transmission spectra and the band structures are discussed. The drawn conclusions may be useful for designing or evaluating the defected phononic crystal beams in bending wave control. In addition, our conclusions are especially potential for identifying the defect location through bending wave signals.


Method of Reverberation-Ray Matrix
Consider the bending waves in a defected bi-coupled phononic crystal beam consisting of n rigidly connected unit cells, as depicted in Figure 1, whose defect is introduced in the {j}th unit cell.All unit cells are comprised of rigidly connected m beams denoted from left to right as (l) to (m), as shown in Figure 2a,b for perfect and defected unit cells, respectively.These beams are long in comparison with the cross-sectional dimensions and carry the bending waves including flexural and thickness-shear modes.They are assumed to be single-phase and to have straight axes as well as uniform cross-section.The joints in any unit cell-i.e., the endpoints and the points where adjacent beams connect-are denoted in order from left to right as 1 throughout M + 1 (M = m).The difference between the perfect and the defected unit cells lies in the cross-sectional dimensions and the material parameters of the (j)th structural beam, as indicated from the comparison between Figure 2a,b.Such described point defect can represent the in-service degradation, the intentionally introduced disorder, or the manufacturing error.The parameters related to the structural beam where the defect is introduced are represented by an apostrophe.As described in Figure 2b, A (j) , I (j) and κ (j) are the area, the second moment of inertia, and the shear correction factor of cross-section of the beam where the defect is introduced; while ρ (j) , E (j) , and G (j) are the density, the Young's modulus, and the shear modulus of the beam material.η v(j) and η ϕ(j) are the damping coefficients associated with the transversal and rotational velocities of the beam in the defected unit cell, respectively.The counterpart parameters without apostrophes are the corresponding parameters of the perfect beams, as described both in Figure 2a

Analysis Model
According to the Floquet-Bloch theorem [4,6,15], if the dispersion characteristics of the perfect phononic crystal beam corresponding to Figure 1 is to be analyzed, the model in Figure 2a can be employed.However, for analyzing the dispersion characteristics or the transmission spectra of the defected phononic crystal beam, it is appropriate to employ the system model in Figure 1 by virtue of the supercell technique [5,65] and the Floquet-Bloch theorem [4,6,15].To differentiate the beams and joints in different unit cells, the cell number is attached afterwards as necessary.For example,

Analysis Model
According to the Floquet-Bloch theorem [4,6,15], if the dispersion characteristics of the perfect phononic crystal beam corresponding to Figure 1 is to be analyzed, the model in Figure 2a can be employed.However, for analyzing the dispersion characteristics or the transmission spectra of the defected phononic crystal beam, it is appropriate to employ the system model in Figure 1 by virtue of the supercell technique [5,65] and the Floquet-Bloch theorem [4,6,15].To differentiate the beams and joints in different unit cells, the cell number is attached afterwards as necessary.For example,

Analysis Model
According to the Floquet-Bloch theorem [4,6,15], if the dispersion characteristics of the perfect phononic crystal beam corresponding to Figure 1 is to be analyzed, the model in Figure 2a can be employed.However, for analyzing the dispersion characteristics or the transmission spectra of the defected phononic crystal beam, it is appropriate to employ the system model in Figure 1 by virtue of the supercell technique [5,65] and the Floquet-Bloch theorem [4,6,15].To differentiate the beams and joints in different unit cells, the cell number is attached afterwards as necessary.For example, (j) − {i} and J − {i} denote the (j) beam and J joint in the {i} unit cell, respectively.It should be emphasized that henceforth we use the {lowercase letter}, (lowercase letter), and uppercase letter to represent the unit-cell, the member, and the joint labels, respectively.For the convenience of describing the system Crystals 2018, 8, 21 5 of 33 and the joints, a global coordinates (XYZ) with its origin at the very left end of unit {1} is established following the right-handed screw rule as in Figure 1.The constituent beams, however, should be described in the local coordinates.
As a unique feature of the method of reverberation-ray matrix (MRRM), a pair of dual local coordinates is set up for any constituent beam, as depicted by (xyz) JK and (xyz) K J in Figure 3a for a typical beam (j) that is sometimes also called as member JK or K J according to the pertinent coordinates.In the constituent beam, the generalized displacements including the deflection v and rotation ϕ as well as the generalized forces including the shear force Q and the bending moment M, can be described in either local coordinate system as represented in Figure 3b by the variables with superscripts of double letters denoting the pertinent coordinates.These quantities are the variables describing the dynamic states of the constituent beam.The bending waves propagating within the constituent beam are also described in dual local coordinates, as indicated in Figure 3c  represent the unit-cell, the member, and the joint labels, respectively.For the convenience of describing the system and the joints, a global coordinates ( ) X YZ with its origin at the very left end of unit {1} is established following the right-handed screw rule as in Figure 1.The constituent beams, however, should be described in the local coordinates.
As a unique feature of the method of reverberation-ray matrix (MRRM), a pair of dual local coordinates is set up for any constituent beam, as depicted by ( ) JK xyz and ( ) KJ xyz in Figure 3a for a typical beam ( ) j that is sometimes also called as member JK or KJ according to the pertinent coordinates.In the constituent beam, the generalized displacements including the deflection v and rotation ϕ as well as the generalized forces including the shear force Q and the bending moment M , can be described in either local coordinate system as represented in Figure 3b by the variables with superscripts of double letters denoting the pertinent coordinates.These quantities are the variables describing the dynamic states of the constituent beam.The bending waves propagating within the constituent beam are also described in dual local coordinates, as indicated in Figure 3c   The physical variables of a joint are described in the global coordinates if they are coordinate dependent.Figure 4a ) are exclusively and respectively denoted by YL u , ZL in order to indicate clearly that they are the coupling variables between adjacent supercells.Note that the translational displacements/forces of joints and of members are deemed positive in the forward coordinate direction.In addition, the rotational displacements/moments of joints and of members are deemed positive around the counter-clockwise.These quantities are clearly indicated in Figures 3b and 4b.The physical variables of a joint are described in the global coordinates if they are coordinate dependent.Figure 4a depicts the global coordinates (XYZ) on the background with the left exterior joint 1 − {1}, any interior joint J, and the right exterior joint [M + 1] − {n}.The joint label will usually be attached to joint variables as superscript to denote their pertinence.As shown in Figure 4b, u J Y and θ J Z are the translational and rotational displacements of the typical interior joint J, respectively, which are combined to form the joint generalized displacements.p J Y and m J Z are the corresponding transversal force and rotational moment, respectively, which are combined to form the joint generalized forces.Nevertheless, the generalized displacements/forces of the left (right) exterior joint 1 − {1} ([M + 1] − {n}) are exclusively and respectively denoted by u YL , θ ZL /p YL , m ZL (u YR , θ ZR /p YR , m ZR ) in order to indicate clearly that they are the coupling variables between adjacent supercells.
Note that the translational displacements/forces of joints and of members are deemed positive in the forward coordinate direction.In addition, the rotational displacements/moments of joints and of members are deemed positive around the counter-clockwise.These quantities are clearly indicated in Figures 3b and 4b.

Analysis of Members and Phase Relations
This subsection analyzes the equations related to the bending wave propagation within the constituent beams.

Governing Equations and Their Wave Solutions
Considering the dynamics of any constituent beam, the equations governing the bending waves can be established based on the Timoshenko beam theory.These equations apply from null up to relatively high frequency ranges since the shear deformation and the rotary inertia are further involved in the Timoshenko beam theory comparing to the Euler-Bernoulli beam theory.In either local coordinates of a structural member, these equations with concrete superscripts being omitted are written as where the involved quantities have been explained previously.Note that as shown in Figure 2, the geometrical and material parameters are all constant for each constituent beam.The damping By virtue of Fourier transform or separation of variables, Equation (1) can be reduced to the wave equations as follows

Analysis of Members and Phase Relations
This subsection analyzes the equations related to the bending wave propagation within the constituent beams.

Governing Equations and Their Wave Solutions
Considering the dynamics of any constituent beam, the equations governing the bending waves can be established based on the Timoshenko beam theory.These equations apply from null up to relatively high frequency ranges since the shear deformation and the rotary inertia are further involved in the Timoshenko beam theory comparing to the Euler-Bernoulli beam theory.In either local coordinates of a structural member, these equations with concrete superscripts being omitted are written as where the involved quantities have been explained previously.Note that as shown in Figure 2, the geometrical and material parameters are all constant for each constituent beam.The damping coefficients η v and η ϕ are usually obtained from η v = 2ξ v ρω and η ϕ = 2ξ ϕ ρω by specifying the damping ratios, ξ v and ξ ϕ , respectively.By virtue of Fourier transform or separation of variables, Equation (1) can be reduced to the wave equations as follows where ω is the circular frequency; c l = E/ρ and c t = κG/ρ are the speeds of quasi-longitudinal and quasi-transversal waves in the beam, respectively; and the over caret henceforth signifies the quantities in the frequency domain (the Fourier transformed quantities).
According to the theory of ordinary differential equations, the solutions to the spectral generalized displacements can be expressed in the form of traveling waves as follows v(x) = (a 2 e ik 2 x + d 2 e −ik 2 x ) + (a 3 e ik 3 x + d 3 e −ik 3 x ), where i = √ −1 is the imaginary unit, and are the wavenumbers of flexural and thickness-shear modes, respectively.Since the square root operation can lead to two values with opposite sign, the wavenumbers k 2,3 in Equation ( 4) can always be chosen to satisfy Im(k 2,3 ) < 0 or Im(k 2,3 ) = 0, Re(k 2,3 ) ≥ 0. Thus, the first and the second terms in the former (latter) parenthesis of Equation ( 3), when combined with the time variation e iωt in the inverse Fourier transform, signify backward and forward waves of flexural (thickness-shear) mode.In MRRM, these waves are also referred to as the arriving and the departing waves, respectively, with a p and d p being the corresponding undetermined amplitudes of flexural (p = 2) or thickness-shear (p = 3) mode.The velocities c 2,3 = ω/k 2,3 of flexural and thickness-shear waves are frequency dependent, which indicates that these waves are dispersive.It is also indicated from Equation (3) that the two waves are coupled between each other, since the solutions to the generalized displacements are the combination of contributions from both waves.
is the influence coefficient of rotation displacement corresponding to p wave mode, as the influence coefficients of transverse displacement corresponding to the two modes are both unity.
Substituting Equation (3) into the frequency-domain constitutive relations of the structural member as follows one obtains the solutions to the generalized forces expressed in form of flexural and thickness-shear waves as where ζ p = (ik p − β p )κGA and µ p = ik p β p EI are the influence coefficients of shear force and of bending moment corresponding to p wave mode, respectively.The solutions in Equations ( 3) and ( 6) based on Timoshenko beam theory will be reduced to those based on Euler-Bernoulli beam theory, as the wavenumbers are expressed ) and the influence coefficients are given by β p = ik p , ζ p = −(ik p ) 3 EI, and µ p = (ik p ) 2 EI (p = 2, 3).

Relations between Wave Solutions in Dual Coordinates
Consider the compatibility conditions between the generalized displacements (generalized forces) in the pair of dual local coordinates of any structural member.Taking the typical member (j) (also JK or K J in the pertinent local coordinates) as shown in Figure 3 for example, at the cross-section x JK in (xyz) JK , i.e., the cross-section l (j) − x JK in (xyz) K J , we have vJK (x JK ) = − vKJ (l (j) − x JK ), φJK (x JK ) = φKJ (l (j) − x JK ), in view of the sign conventions of the generalized displacements and forces.

Phase Relations of a Member
Substituting the solutions to the generalized displacements (forces) provided in Equation (3) (Equation ( 6)) into the former (latter) two formulas of Equation (7), noticing that the wavenumber, the length, and the influence coefficients of member (j) are irrelevant to the local coordinates, i.e., k JK p = k KJ p = k p(j) , l JK = l KJ = l (j) , and , and observing the functions e ik p x and e −ik p x cannot be zero, one obtains the local phase relations of the typical member (j) where T are the arriving and departing wave vectors of typical member (j), respectively, with T being the corresponding vectors in respectively pertinent local coordinates as shown in Figure 3.
T is a variant of the arriving wave vector of member (j), which has the same components but arranged in different sequence as a (j) .U (j) = [0, I 2 ; I 2 , 0] = [U (j) ] −1 is an orthogonal permutation matrix used to transform a (j) into a (j) , with I 2 the second order identity matrix.P (j) =< P JK , P K J > is the phase matrix of member (j) with P JK = P K J =< −e −ik 2(j) l (j) , −e −ik 3(j) l (j) > the phase matrices of member JK and K J, respectively.Here and after, the superscript T indicates the transpose of a vector or matrix, and < • > denotes the (block) diagonal matrix with elements (or sub-matrices) only on the main diagonal.It should be noted that the chosen wavenumbers k p(j) in Equation ( 4) can always guarantee Re[−ik p(j) l (j) ] ≤ 0. Thus, no exponentially growing function appears in Equation (8).

Global Phase Relations of the System
Grouping the phase relations of all members in sequence of the unit-cell number and the beam number, the global phase relations are obtained as follows where a G and d G are respectively the global arriving and departing wave vectors composed of a G is a variant of the global arriving wave vector a G , which is composed of and related to a G by the 4mn × 4mn orthogonal global permutation matrix U G formed as Crystals 2018, 8, 21 9 of 33 P G is the 4mn × 4mn diagonal global phase matrix composed of Note that the phase matrices of the counterpart members in all unit cell, except the one where the defect is introduced, are actually identical.Therefore, it will not cost too much effort to form P G .
It should also be pointed out that the exponentially growing function is avoided in the phase relations given by Equations ( 8) and ( 9).This leads to the unconditionally numerical stability of our MRRM formulation from one aspect.
The global phase relation gives altogether 4mn equations for the 4mn unknown arriving and the 4mn unknown departing wave amplitudes in respective vectors a G and d G .Obviously, another 4mn equations should be supplemented to determine all the unknowns.

Analysis of Joints and the Scattering Relations
This subsection analyzes the equations related to the bending wave scattering at the structural joints.

Continuous Equations and Scattering Relations at an Interior Joint
Refer to Figure 4, the spectral equations of compatibility of displacements and equilibrium of forces at a typical interior joint J should be which serve as the governing equations of joint J. Substitution of the solutions to the generalized displacements and forces (in Equation ( 3) and Equation (6), respectively) into Equation ( 14) leads to the local scattering relation of interior joint J as follows where T are the arriving and the departing wave vectors of the typical interior joint J, respectively.
J K are the scattering matrix and the source vector of joint J, respectively, with D J , A J and v J K being Note that ) are the influence coefficients (of rotation, shear force, and bending moment) for members (j − 1) and (j), respectively.If the dispersion characteristics are analyzed, the applied generalized forces should not be considered, and the source vector s J in this case is a vector of zeros.

Coupling Equations and Scattering Relations of the Pair of Exterior Joints
Also refer to Figure 4, the spectral equations describing the compatibility of displacements and the equilibrium of forces at the left exterior joint 1 − {1} and the right exterior joint [M + 1] − {n} are expressed respectively as For the dispersion characteristic analysis and for the transmission spectra analysis, these two equations are subsequently treated in different ways, so split discussion is given as follows.

Dispersion characteristic analysis:
Due to Floquet-Bloch theorem [4,6,15], periodic conditions exist between the counterpart generalized displacements/forces of the exterior joint pair 1 − {1} and where µ and q wavenumber of characteristic wave in the supercell respectively, and L is the length of the supercell.The imaginary and the real parts of µ denoting by µ I and µ R , i.e., the real and the imaginary parts of qL denoting by q R L and q I L, are known as the phase constant and the attenuation constant, respectively.It should be stated that both µ and q can be used to represent the feature of the characteristic wave, the wavenumber q is preferred in our study for its closer connection to the wavelength and phase velocity.Substitution of Equations ( 17) and ( 18) into Equation (19) gives the coupling equations of exterior joint pair from which the scattering relations are further obtained as where and ) are the arriving and the departing wave vectors of exterior joint 1 − {1} ([M + 1] − {n}), respectively.The scattering matrix can be computed from where Γ

Transmission spectra analysis:
To compute the transmission, the type of the left (right) exterior joint 1 − {1} ([M + 1] − {n}) should first be specified.The left and the right exterior joints can be any one among the four possible joint types shown in Tables 1 and 2, respectively, where the corresponding known quantities forming the excitation vector and the unknown quantities vector are also provided.For the left (right) exterior joint, two formulas associated with the known quantities in Equation ( 17) (Equation ( 18)) are used to determine the wave scattering at the joint.Taking the free end of left exterior joint as an example, the two equilibrium equations of forces in Equation ( 17) after the solutions to the generalized forces of member 12 − {1} are introduced, lead to the scattering relation of left exterior joint 1 − {1} where the wave vectors are identical to those given before, the scattering matrix can be computed from the matrices and vector provided in Table 1.For other types of left exterior joint, the scattering relations can be expressed in the same form as Equation ( 24), but with the scattering matrix and the source vector obtained from the other corresponding form of D 1−{1} , A 1−{1} , and v Similarly, the scattering relations of the right exterior joint are written consistently as where the wave vectors are the same as before, and the scattering matrix −{n} and the source vector can be obtained from the matrices and vectors given in Table 2.
Table 1.The known and unknown quantities of four types of left exterior joint 1 − {1} as well as the matrices and vector for computing the scattering matrix and the source vector.− {n} as well as the matrices and vector for computing the scattering matrix and the source vector.

. Global Scattering Relations of the System
Grouping the scattering relations of all joints in order of the unit-cell number and the joint number, one obtains the global scattering relation where the global arriving and departing wave vectors in form of are exactly the same as those given in Equation (10).The 4mn × 4mn global scattering matrix S G and the 4mn-dimensional source vector s G are formed as where s G = 0 4mn×1 for the dispersion characteristics analysis and for the transmission spectra analysis.Note that the scattering matrices of the counterpart joints in all unit cells, except those of the exterior joints and those of the joints adjacent to the defected location, are actually identical.Therefore, it will not take too much efforts to form S G .
It should also be emphasized that the exponential functions in the solutions to generalized displacements and forces of structural members disappear in the scattering relations, as indicated in Equations ( 16) and ( 22) and Tables 1 and 2, because the axial coordinates of the joints in the adopted local coordinate systems are always zero for all correlated structural members.This is the main advantage of introducing the dual local coordinates in the MRRM.The uniquely appeared exponential function e iqL in the scattering relations (in Equation (22) as the dispersion characteristics are analyzed), can be easily specified to avoid exponentially growing case.In this way, numerical stability can be unconditionally realized in the MRRM from the other aspect.
The global scattering relation in Equation (26) supplements the global phase relation with another 4mn equations for the 4mn unknown arriving wave amplitudes in a G and the 4mn unknown departing wave amplitudes in d G , thus the wave vectors a G and d G can be uniquely determined by combining the two kinds of global relations.

Analysis of Bending Wave Dispersion and Transmission from the System Equation
Substitution of Equation ( 9) into Equation (26) gives the system equation where the 4mn × 4mn R(ω, q) = S G (ω, q)U G P G (ω) is referred to as reverberation-ray matrix.
The coefficient matrix [I − R] in this system equation is called the system matrix.
If the dispersion characteristics are analyzed, s G = 0 4mn×1 satisfies and the system equation becomes homogeneous.For existing non-trivial solutions, the determinant of the system matrix should vanish to give det[I − R(ω, q)] = 0, (30) which provides the dispersion relation and determines the feature of the characteristic bending waves in the defected bi-coupled phononic crystal beam.If any one of three in frequency ω, wavenumber q (or wavelength λ = 2π/q) and phase velocity c = ω/q is specified, the other two can be determined by properly solving Equation ( 30) through root searching technique, and then the dispersion spectra representing the propagation feature of characteristic bending waves can be plotted.The most commonly used dispersion curves are the ω − q spectra, which are also called the band structures.
If the transmission spectra are analyzed, S G should be ω−dependent but q−irrelevant (inferred from Tables 1 and 2) and should be blocked-diagonal because of S (1−{1})([M+1]−{n}) = 0 2×2 and S ([M+1]−{n})(1−{1}) = 0 2×2 .Thus, the system equation is in-homogeneous and gives the solution to the global departing wave vector d G as Back substituting the solution to d G into the global phase relation in Equation (9) gives the solution to a G .By further employing Equations ( 3) and ( 6), the solutions of the generalized displacements and forces at any cross-section of whichever beam member in whatever unit cell can be determined, respectively.After conducting the abovementioned calculations circularly for each frequency ω in a specified range with self-given increment, the transmission spectra can then be plotted.About this process, two points should be emphasized.The first point is that the generalized inversion (Moore-Penrose pseudo inversion) may be used to solve Equation (31) for avoiding program interruption as the system matrix becomes ill-conditioned.The second point is that when the generalized forces of all interior joints vanish, the commonly required transmission that relates the known quantity of one exterior joint to the unknown quantity of the other exterior joint can be given.

Validation of Theoretical Method by Experimental Measurement and Numerical Simulation
To verify the above-derived MRRM, we consider a defected phononic crystal beam with eight binary unit cells, as depicted in Figure 5a.The corresponding specimen actually made for the experiment is shown in Figure 5b.All unit cells in the considered structure consist of a PMMA (polymethyl methacrylate) beam and an aluminum beam with the same length 80 mm.The defect is in the cross-sectional height h of the PMMA constituent beam of the fifth unit cell counted from left.All the PMMA and aluminum beams have the cross-section of b × h = 15 mm × 15 mm (b is the cross-sectional width), except that the PMMA beam in the defected unit cell has the cross-section of 15 mm × 8 mm.The difference between the PMMA beams in the perfect and the defected unit cells can be clearly seen from the comparison in Figure 5a and from the contrast between Figure 5c,d.
The geometrical and material parameters of these constituent beams are listed in Table 3 for the convenience of subsequent calculations.Among the geometrical parameters, the length and the cross-sectional dimensions are measured.Considering the cross-section is rectangular, the area and the second moment of inertia are calculated and the shear coefficient is directly specified.Among the material parameters, the Young's modulus of the PMMA is experimentally obtained beforehand by matching the steel ball impact-induced frequency responses of a PMMA plate to those predicted by the finite element simulations [65].All the other material parameters are given by referring to a materials handbook [71].Note that, except the calculations in the following first part, the computations thenceforth do not consider the damping of the phononic crystal beam by setting all the damping ratios as zero.Also, for the convenience of engineering applications, the engineering frequency f = ω/(2π), the dimensionless wavenumber qL/π, and the decibel (dB) transmission are used hereafter.Table 3.The geometrical and material parameters of the constituent beams of the defected phononic crystal beam in Figure 5. Firstly, using the MRRM proposed in Section 2 on the basis of Timoshenko beam theory, the transmission coefficients relating the deflection at the left exterior joint (output end) to that at the right exterior joint (input end) are calculated within frequency range f ≤ 8 kHz.Five cases are specified separately for the damping ratios ξ v = ξ ϕ of all PMMA and aluminum beams, i.e., 0, 0.005, 0.01, 0.015, and 0.02, so that the effect of damping loss on the transmission spectra can be inspected.The resulting transmission spectra are provided and compared with the counterpart results by experiment as well as by FEM with various element sizes in Figure 6.The experiment is based on the fiber Bragg grating (FBG) displacement sensing system proposed in [65], where the other setups and the experiment process are also given.The FEM simulation is performed by the COMSOL Multiphysics ® software (Comsol Inc., Stockholm, Sweden), and the modeling process can be referred to literature [72].The three-dimensional (3D) tetrahedral element is utilized in the modeling.We separately adopt four cases of meshing schemes-i.e., the fine, extra fine, extremely fine, and extremely finer situations predefined in COMSOL-in order to perform the mesh refinement study.The total numbers of elements corresponding to the four situations are 1499, 2068, 4324, and 6227, respectively.Figure 6a compares the result by the MRRM theoretical analysis with the damping ratios being 0.01, the result by the FBG experiment, and the result by the FEM with the isotropic structural loss factors being 0.01 [65] and with 6227 elements.The FEM result corresponding to the most refined mesh (i.e., extremely finer mesh) is selected in Figure 6a for comparison, because it should be the most accurate FEM result among the four cases of element sizes according to the convergence criterion of FEM analysis [68]. Figure 6a indicates that the transmission spectra obtained by MRRM theory, by FBG experiment, and by FEM simulation are agreed in general within the considered frequency range.In particular, excellent agreement is discerned below 3.5 kHz.The agreement of these results in Figure 6a can validate simultaneously the proposed MRRM, the FBG experiment, and the FEM simulation.However, the discrepancy between these transmission spectra above 3.5 kHz, as shown in Figure 6a, does exist and becomes more distinguishable with the increasing frequency.The difference between the result by MRRM together with that by FEM and the experimental measurement may be due to the fact that the ambient noise in the experiment is not considered in the computational models.That the material constants used in the computational models may not exactly identical to the actual values in the experimental specimen can also cause the difference.For instance, the comparison between the MRRM results associated with various damping ratios and the experimental result in Figure 6b demonstrates that the damping parameter has very distinct influence on the computational transmission spectra.Nevertheless, it is pretty difficult to determine the damping parameter exactly since quite a lot of mechanisms are related to the energy loss in a practical structure [73].Figure 6b also illustrates that although the resonant amplitudes decrease obviously with the increasing of damping ratio, the resonant frequencies are nearly unchanged with the damping ratio.Meanwhile, the defect modes become harder and harder to be discerned with the increasing of damping ratio.The discrepancy between the MRRM result and the other two results may be caused by three reasons.The first reason is that the MRRM uses the one-dimensional (1D) Timoshenko beam theory to model the structural beams.Nevertheless, the experiment specimen is actually 3D solid, and the FEM simulation adopts 3D tetrahedral element to assemble the 3D solid structure.The second reason is that the MRRM uses viscous damping to model the energy loss mechanism.Nevertheless, the experiment specimen actually possesses multifarious energy loss mechanisms, and the FEM simulation adopts the hysteretic damping model [73].To avoid introducing this damping factor that causes the difference, in the following analysis zero damping ratios or loss factors are specified hereafter.The third reason is that the FEM simulation using 6227 elements may have not yet provided convergent results at high frequency ranges (e.g., above 4 kHz).To explore the convergence of the FEM simulation and to know more about the discrepancy between the FEM and the MRRM results, Figure 6c compares the undamped transmission spectra obtained by the FEM with four meshing schemes and also the result by the MRRM.It is seen from Figure 6c that with the increasing of the total elements in the FEM model, the FEM results shift to the lower frequency side and move closer to the MRRM result.This fact is in accordance with the convergence process of the FEM analysis [68].The discrete FEM model generally overvalues the stiffness of the real structure for the following reason.The FEM model has finite degree of freedom (DOF), while the real structure has infinite DOF.Thus, modeling the real structure by a FEM model is equivalent to adding extra constraints to the real structure.Since adding constraints to a structure will increase its stiffness.Consequently, the stiffness of a FEM model is always bigger than that of the associated real structure.Thus, the FEM analysis gives higher resonant frequencies than the exact ones.The fact seen from Figure 6c also implies that the MRRM result is very close to the exact solution of the 3D solid structure.Consequently, the proposed MRRM works very well for the analysis of transmission coefficients of defected phononic crystal beams.In addition, the excellent agreement of the result corresponding to 6227 elements and the result associated with 4324 elements below 3 kHz confirms the convergence of both results within 3 kHz.However, the distinguishable difference between the two results above 3 kHz demonstrates that the result associated with 4324 elements has not converged yet and the convergence of the FEM results becomes slower for higher frequency analysis.From the above comparison between the results of FEM with various elements and the result of MRRM, it is known that for providing satisfactory result especially at high frequency range, the FEM needs a tremendous number of elements to model each constituent beam, while the MRRM needs only one member to model each constituent beam.So the main advantage of MRRM versus FEM should be that the MRRM needs less computational effort.To explore the convergence of the FEM simulation and to know more about the discrepancy between the FEM and the MRRM results, Figure 6c compares the undamped transmission spectra obtained by the FEM with four meshing schemes and also the result by the MRRM.It is seen from Figure 6c that with the increasing of the total elements in the FEM model, the FEM results shift to the lower frequency side and move closer to the MRRM result.This fact is in accordance with the convergence process of the FEM analysis [68].The discrete FEM model generally overvalues the stiffness of the real structure for the following reason.The FEM model has finite degree of freedom (DOF), while the real structure has infinite DOF.Thus, modeling the real structure by a FEM model is equivalent to adding extra constraints to the real structure.Since adding constraints to a structure will increase its stiffness.Consequently, the stiffness of a FEM model is always bigger than that of the associated real structure.Thus, the FEM analysis gives higher resonant frequencies than the exact ones.The fact seen from Figure 6c also implies that the MRRM result is very close to the exact solution of the 3D solid structure.Consequently, the proposed MRRM works very well for the analysis of transmission coefficients of defected phononic crystal beams.In addition, the excellent agreement of the result corresponding to 6227 elements and the result associated with 4324 elements below 3 kH z confirms the convergence of both results within 3 kH z .However, the distinguishable difference between the two results above 3 kH z demonstrates that the result associated with 4324 elements has not converged yet and the convergence of the FEM results becomes slower for higher frequency analysis.

Cross-Sectional Area
From the above comparison between the results of FEM with various elements and the result of MRRM, it is known that for providing satisfactory result especially at high frequency range, the FEM needs a tremendous number of elements to model each constituent beam, while the MRRM needs only one member to model each constituent beam.So the main advantage of MRRM versus FEM should be that the MRRM needs less computational effort.
(a) Secondly, in the low frequency range ( ), the above undamped transmission spectra of the defected phononic crystal beam by the MRRM are compared in Figure 7a with those of the corresponding perfect structure computed also by the MRRM.For the verifying purpose, the measured transmission spectra by the FBG experiment is also depicted in Figure 7a.In addition, the band structures ( q ω − dispersion curves) of the defected phononic crystal beam as a supercell are computed in the same frequency range by the TMM and by the MRRM with the help of the supercell technique [5,65].The resulting band structures are depicted in Figure 7b,c, respectively.To simultaneously examine the effect of the supercell technique and the influence of the defect, the band structures of the corresponding perfect phononic crystal beam within Secondly, in the low frequency range ( f ≤ 2.5 kHz), the above undamped transmission spectra of the defected phononic crystal beam by the MRRM are compared in Figure 7a with those of the corresponding perfect structure computed also by the MRRM.For the verifying purpose, the measured transmission spectra by the FBG experiment is also depicted in Figure 7a.In addition, the band structures (ω − q dispersion curves) of the defected phononic crystal beam as a supercell are computed in the same frequency range by the TMM and by the MRRM with the help of the supercell technique [5,65].The resulting band structures are depicted in Figure 7b,c, respectively.
To simultaneously examine the effect of the supercell technique and the influence of the defect, the band structures of the corresponding perfect phononic crystal beam within f ≤ 2.5 kHz are also computed by the MRRM, with the results based on an eight unit-cells model and one unit-cell model provided in Figure 7d,e, respectively.It needs to point out that here and after the attenuation constant spectra corresponding to the evanescent modes below cut-off frequencies of structural members are not entirely depicted within the considered frequency range, since it is nonsignificant to the frequency bands.Also note that the shadowed frequency range here in Figure 7 and after is the band gaps of the perfect phononic crystal beam.
Figure 7a demonstrates that introducing a point defect into the phononic crystal beam alters the resonant frequencies in the sense of finite structure, which can be seen from the peak shifting of the transmission spectra.Generally speaking, the transmission spectra at higher frequency ranges are more sensitively influenced by the defect.The resonant frequencies shift to the lower frequency side, because the defect introduced by reducing the cross-sectional height of the PMMA constituent beam in the fifth unit cell, actually causes a decrease of each order of natural frequencies of the whole finite structure.The shifting of the resonant frequency that is very close to the low bounding frequency of the passband into the adjacent lower-frequency stopband of the original perfect phononic crystal beam forms the localized (defect) mode, as depicted in Figure 7a.This process of forming the localized (defect) mode due to the above-introduced defect is very similar to the formation of the localized (defect) mode in a diatomic crystal by introducing a defected heavy atom that has been well-discussed in the solid state theory [74].Note that the reduction of cross-sectional height simultaneously decreases the stiffness (EI, κGA) and the inertia (ρI, ρA) of the structural beam, but the decreasing of the stiffness is much pronounced than the decreasing of the inertia since the moduli of PMMA are about six orders bigger than its density.Thus, the general effect of the defect introduced by reducing the cross-sectional dimensions of one constituent beam in the phononic crystal beam is consistent with the influence of the defect introduced by replacing an original mass by a heavier mass in the crystal lattice [74].Figure 7b shows that the TMM becomes numerically instable for analyzing dispersion curves over 250 Hz, while Figure 7c indicates that the MRRM is always numerically stable in the whole considered frequency range.The comparison of the band structures in Figure 7d,e indicates that, for the perfect phononic crystal beam, the supercell technique does not change the frequency distribution of passbands and stopbands, but does alter the phase constant spectra in the passbands and the attenuation constant spectra in the stopbands.The phase constant spectra in each passband associated with supercell containing n unit cells can be obtained by folding n times the phase constant spectra associated with one unit cell within the passband frequency range.This fact has already been noticed in [65] and should be attributed to the zone folding effect [70].The attenuation constant of the supercell containing n unit cells is the n time magnification of that of one unit cell in respective stopband.The comparison of the band structures in Figure 7c,d illustrates in the sense of infinite structure that introducing defect opens new stopbands within each of the original passbands and may form new propagating modes (defect modes) within the original stopbands of the perfect phononic crystal beam.The newly appeared stopbands and defect passbands caused by introducing a single defect in the original perfect phononic crystal beam have also been noticed in [65].They are probably formed due to more diversified wave interference phenomena induced by adding new types of reflected and transmitted waves at the interfaces next to the defected PMMA constituent beam, which are the same as the reason that the band structures of the perfect phononic crystal beam are formed [70].The comparison of Figure 7a with Figure 7c,d further illustrates the coincidences of the wave properties predicted by the transmission and by the dispersion for respective type of phononic crystal beam.These coincidences on one hand validate that the proposed MRRM works very well for the analyses of transmission and dispersion of the characteristic bending waves within the considered frequency range.On the other hand, they confirm that when defects are introduced the relations between the natural frequencies of the finite phononic crystal beam and the band structures of the counterpart infinite phononic crystal beam are identical to the previously found relations for the perfect periodic beams [75].Consequently, the effects of the point defect represented by the transmission spectra in the sense of finite structure are inherently related to its effects expressed by band structures in the sense of an infinite structure.Thirdly, to verify that the proposed MRRM can be applied for high-frequency analysis, the undamped transmission spectra of both the defected and the perfect phononic crystal beams obtained by the MRRM, the band structures computed by the TMM and the band structures calculated by the MRRM of the defected phononic crystal beam within 29 kHz ≤ f ≤ 38 kHz are provided in Figure 8a-c, respectively.The counterpart band structures of the perfect phononic crystal beam based on supercells containing eight unit cells and one unit cell are shown in Figure 8d,e, respectively.This frequency range contains the ninth and tenth passbands as well as the ninth stopband of the perfect phononic crystal beam, as can be judged from Figure 8f that provides the band structures of the perfect phononic crystal beam within f ≤ 47 kHz based on one unit cell.Therefore, the frequency range 29 kHz ≤ f ≤ 38 kHz is appropriate to represent the high frequency range.In Figure 8, the coincidence of the frequency bands predicted by the transmission and by the dispersion of respective type of phononic crystal beam is illustrated within 29 kHz ≤ f ≤ 38 kHz. Figure 8a demonstrates that the above-found effect of the defect on the shifting of transmission spectra in the sense of finite structure also holds in the high frequency range.Figure 8b indicates that the TMM is infeasible for the dispersion characteristic analysis in the considered high frequency range, while Figure 8c denotes that the MRRM works very well for the analysis of high-frequency dispersion.The comparisons of Figure 8c-e illustrate that the above-found effect of the supercell technique and influence of the defect on the band structures in the sense of infinite structure also hold in the high frequency range.
Fourthly, the validated MRRM is used to demonstrate the dependence of the computed band structures on the beam theories, i.e., the Timoshenko and Euler-Bernoulli beam theories.The computed band structures within f ≤ 2.5 kHz of the defected phononic crystal beam on basis of the two beam theories are compared in Figure 9.The contrast shows that at the low frequency range f ≤ 0.5 kHz, the Timoshenko and Euler-Bernoulli beam theories lead to equivalent results.Above 0.5 kHz, the discrepancy becomes observable and increases with the frequency.Note that the Euler-Bernoulli beam theory always leads to band structures with exaggerated frequency, because this beam theory, which ignores the shear deformation and rotary inertia, actually overvalues the stiffness of the structural member.Therefore, the analysis based on the Euler-Bernoulli beam theory ceases to be effective in relatively high frequency ranges, while the analysis based on the Timoshenko beam theory is applicable for high-frequency analysis since the shear deformation and rotary inertia are included.theory is applicable for high-frequency analysis since the shear deformation and rotary inertia are included.

Results and Discussion
Using the validated MRRM based on the Timoshenko beam theory and taking the abovementioned perfect phononic crystal beam as the benchmark, this section studies the influences of factors including the defect size, the defect form, and the unit-cell number on the transmission spectra and the band structures.Since it is noticed from the above discussions that the transmission spectra and the band structures at the low and high frequency ranges exhibit a similar pattern, here we focus the frequency within . The influence rule can be summarized by comparing the transmission spectra and band structures of various defected phononic crystal beam with those of the benchmark perfect phononic crystal beam.

Effect of Defect Size and Defect Form
Two forms of defects-i.e., the defect in the cross-sectional dimension (height) and the defect in the material parameter (moduli) of the left beam of the fifth unit cell-are considered here to inspect the effect of defect size and defect form on the transmission spectra and band structures.Note that the weakened defect may represent the in-service degradation, and both the weakened and strengthened defects can describe the intentionally introduced disorder or manufacturing error.In the first case of the two defect forms, the cross-sectional height is specified as 12 mm and 18 mm in addition to the already considered 8 mm and 15 mm , while the cross-sectional width remains unchanged.If the cross-sectional height remains unchanged while the width varies or both the height and width vary, a similar influence rule could be anticipated since actually the cross-sectional area and the second moment of inertia are utilized in the calculation.Hence, these situations are not discussed here to save space.For phononic crystal beams with four kinds of cross-sectional height, the transmission spectra are compared in Figure 10a within are given in Figure 10b and 10c, respectively, because the first several bands and the first defect mode can be clearly illustrated at respective frequency ranges.In the second case of the two defect forms, the Young's modulus is specified as 1.0 GPa , 3.0 GPa , and 6.0 GPa in addition to the already considered 4.5 GPa , while the Poisson's ratio remains the same as

Results and Discussion
Using the validated MRRM based on the Timoshenko beam theory and taking the above-mentioned perfect phononic crystal beam as the benchmark, this section studies the influences of factors including the defect size, the defect form, and the unit-cell number on the transmission spectra and the band structures.Since it is noticed from the above discussions that the transmission spectra and the band structures at the low and high frequency ranges exhibit a similar pattern, here we focus the frequency within f ≤ 1500 Hz.The influence rule can be summarized by comparing the transmission spectra and band structures of various defected phononic crystal beam with those of the benchmark perfect phononic crystal beam.

Effect of Defect Size and Defect Form
Two forms of defects-i.e., the defect in the cross-sectional dimension (height) and the defect in the material parameter (moduli) of the left beam of the fifth unit cell-are considered here to inspect the effect of defect size and defect form on the transmission spectra and band structures.Note that the weakened defect may represent the in-service degradation, and both the weakened and strengthened defects can describe the intentionally introduced disorder or manufacturing error.In the first case of the two defect forms, the cross-sectional height is specified as 12 mm and 18 mm in addition to the already considered 8 mm and 15 mm, while the cross-sectional width remains unchanged.If the cross-sectional height remains unchanged while the width varies or both the height and width vary, a similar influence rule could be anticipated since actually the cross-sectional area and the second moment of inertia are utilized in the calculation.Hence, these situations are not discussed here to save space.For phononic crystal beams with four kinds of cross-sectional height, the transmission spectra are compared in Figure 10a within f ≤ 1500 Hz.The band structures within f ≤ 150 Hz and 300 Hz ≤ f ≤ 1100 Hz are given in Figure 10b,c, respectively, because the first several bands and the first defect mode can be clearly illustrated at respective frequency ranges.In the second case of the two defect forms, the Young's modulus is specified as 1.0 GPa, 3.0 GPa, and 6.0 GPa in addition to the already considered 4.5 GPa, while the Poisson's ratio remains the same as before.The shear modulus G = E/[2(1 + ν)] is computed accordingly.The transmission spectra within f ≤ 1500 Hz and the band structures within f ≤ 150 Hz and within 300 Hz ≤ f ≤ 1100 Hz of the phononic crystal beams with four kinds of material modulus are compared in Figure 11a-c, respectively.Figure 10a illustrates that, referring to the transmission spectra of the benchmark perfect phononic crystal beam, the transmission spectra of the defected phononic crystal beam move toward higher and lower frequency sides as the defects are introduced by respectively increasing and decreasing the cross-sectional height of the left beam in the fifth unit cell.This is due to the increasing and decreasing of the resonant frequencies of the finite structure in the two defected cases, whose reason has been discussed after Figure 6.It is interesting to note that defect mode occurs as the height reduces to as low as 12 mm and 8 mm.However, defect mode does not appear as the height enlarges to 18 mm, because the resonant frequency that is very close to the high bounding frequency of the passband has not shifted into the adjacent higher-frequency stopband of the original perfect phononic crystal beam.Figure 10b,c shows that the defect introduces new stopbands within the original passband of the perfect phononic crystal beam by opening the phase constant spectra at q R L/π = N (N is an integer) and by slightly changing the shape in-between adjacent q R L/π = N.As the defect is formed by reducing the cross-sectional dimension, the phase constant spectra generally move toward the lower frequency side during the process of opening at q R L/π = N and of slightly changing the shape in between adjacent q R L/π = N.This is because that the flexural and shear stiffness of the structural beam decrease much more obviously than the decreasing of rotary and transversal inertias as the cross-sectional dimension reduces.Thus, the characteristic frequencies of the whole system decrease.The greater the reduction, the lower the characteristic frequencies of the whole system are, and the farther is the movement of the phase constant spectra.However, a different order of phase constant spectrum moves at different rate, thus new stopbands with various bandgaps forms between adjacent phase constant spectra.For the newly introduced stopbands, the band gaps are generally wider and the stopband central frequencies are generally lower as the defect degree is bigger.Forming the defect by increasing the cross-sectional dimension shows opposite influence rule except that 'the bigger the defect degree, the wider the stopbands' still holds.At the same time, a new passband may occur within the original stopband of the perfect phononic crystal beam, as seen from the spectra associated with 12 mm or 8 mm height.This is due to the fact that introducing a defect causes the original different orders of phase constant spectra to move at different degrees, and coincidently some phase constant spectra may move into one original stopband to form new defect modes.Consequently, there is also a possibility that no phase constant spectra move into the original stopband, as can be seen from the spectra associated with a height of 18 mm.
The comparisons of counterpart subfigures in Figures 10 and 11 indicate that introducing defect by reducing (increasing) the material moduli of one constituent beam has the same effect on the transmission spectra and band structures as bringing defect by reducing (increasing) the cross-sectional dimension.Because reducing (increasing) the material moduli of one constituent beam has actually led to the reduction (increasing) of the flexural/shear stiffness without changing the inertias of that beam and further the reduction (increasing) of the characteristic frequency of the whole system.This effect plays the same role as reducing (increasing) the cross-sectional dimension of that constituent beam.Consequently, inversion of the defect by only observing the variations of the transmission spectra and the band structures cannot give out a unique solution.

The Effect of Unit-Cell Number
As all the geometrical and material parameters remain the same as those in Table 3, the unit-cell number n of the defected phononic crystal beam varies from 8 into 6, 4 and 2 to inspect the effect of unit-cell number on the transmission spectra and the band structures.The defect PMMA beam in all the four cases is located next to and on the right side of the finite structure center.The transmission spectra and the band structures within f ≤ 1500 Hz of the four defected phononic crystal beams (with two, four, six, and eight unit cells) are given and compared in Figure 12a,b, respectively.

Conclusions
Based on the Timoshenko beam theory, the method of reverberation-ray matrix (MRRM) is proposed in this paper for the theoretical analysis of dispersion and transmission of bending waves in bi-coupled phononic crystal beams with defects consisting of many phases of materials.The defects can be either in one or many parameters of one or many structural members, corresponding to single (point) or multiple defects.This work only considers the influence of point defects.The defects appearing in the cross-sectional dimensions or in the material parameters can be taken into consideration to model disorders like in-service degradation, intentional design, or manufacturing error.The proposed MRRM are validated by comparing the theoretically computed transmission spectra with the experimentally measured counterparts (based on the FBG displacement sensing system) and with the FEM simulation (based on the COMSOL Multiphysics ® software) results.Numerical examples are provided to study the effects of beam theories, defect occurrence, defect size, defect type, and unit-cell number on the transmission spectra and the band structures of bending Figure 12a illustrates the zone folding effect [70] of the transmission spectra as the unit-cell number increases.Figure 12b demonstrates the folding of band structures by increasing the unit-cell number.Nevertheless, the folding does not occur in exactly corresponding frequency ranges, since the defect also influences the evolution of transmission spectra and band structures of the defect phononic crystal beam.In a fixed frequency range, the transmission spectra of the defect phononic crystal beam with more units of cells show more resonant peaks due to more occurrence of wave scattering.In general, the more are the units cells in the finite phononic crystal beam with defect, the more and narrower are the frequency bands in a specified frequency range.

Conclusions
Based on the Timoshenko beam theory, the method of reverberation-ray matrix (MRRM) is proposed in this paper for the theoretical analysis of dispersion and transmission of bending waves in bi-coupled phononic crystal beams with defects consisting of many phases of materials.The defects can be either in one or many parameters of one or many structural members, corresponding to single (point) or multiple defects.This work only considers the influence of point defects.The defects appearing in the cross-sectional dimensions or in the material parameters can be taken into consideration to model disorders like in-service degradation, intentional design, or manufacturing error.The proposed MRRM are validated by comparing the theoretically computed transmission spectra with the experimentally measured counterparts (based on the FBG displacement sensing system) and with the FEM simulation (based on the COMSOL Multiphysics ® software) results.Numerical examples are provided to study the effects of beam theories, defect occurrence, defect size, defect type, and unit-cell number on the transmission spectra and the band structures of bending waves in defected phononic crystal beams.The influence of damping on the transmission spectra has also been discussed.From these analyses, we can draw some conclusions as follows: (1) From the comparisons of computational results and element requirements of the MRRM and those of the FEM, it is found that the main advantage of MRRM versus FEM should be that the MRRM needs less computational effort.From the comparison of the results by the MRRM and by the TMM (transfer matrix method) both at low and high frequency ranges, it has been proven that the proposed MRRM indeed eliminates the numerical instability problem in the TMM.Based on the Timoshenko or the Euler-Bernoulli beam theory, the proposed MRRM can be applied to calculate the dynamics of phononic crystal beams from zero up to high frequency as long as the respective beam theory is still valid to model the constituent beams.(2) According to the fact that defects actually increase (decrease) the natural/characteristic frequencies, the transmission/phase-constant spectra shift to the higher (lower) frequency side with respect to those of the perfect counterpart.A localized mode forms as a resonant frequency that is very close to the high (low) bounding frequency of the passband shifts into the adjacent higher (lower) stopband.The defect modes become harder and harder to be discerned from the transmission spectra as the damping in the defected phononic crystal beam increases, because the resonant amplitudes decrease more obviously while the resonant frequencies are nearly unchanged with the increasing of damping.In addition, the effect of the point defect seen from the band structures for infinite structures can be reflected by the transmission spectra for practical finite structures.The transmission spectra and the band structures at higher frequency ranges are generally more sensitive to the point defect.
(3) The defect in the cross-sectional dimensions has identical effect on the transmission spectra and the band structures as the defect in the material moduli does.The defects introduced by predominantly reducing the flexural/shear stiffness of the structural beam, like reducing the cross-sectional dimension or the material moduli, move the transmission spectra and the phase constant spectra in each passband towards the lower frequency direction, and vice versa.The bigger the defect degree is, generally the wider the bandgaps are.(4) Increasing the unit-cell number by n times will folding n times the transmission spectra and the phase constant spectra in each passband, and will magnify the attenuation constants n times in the stopbands.The larger is the unit-cell number, the more and narrower are the frequency bands in a specified frequency range.

Figure 1 .
Figure 1.A bi-coupled periodic beam consisting of n unit cells with a defect in the { } j th unit cell.

Figure 2 .
Figure 2. The configuration of: (a) The perfect unit cell { } i ; (b) The defected unit cell { } j .
denote the ( ) j beam and J joint in the { } i unit cell, respectively.It should be emphasized that henceforth we use the { } lo w e r c a s e le tte r , ( ) lo w e r c a s e le tte r , and u p p e r c a s e le tte r to

Figure 1 .
Figure 1.A bi-coupled periodic beam consisting of n unit cells with a defect in the {j}th unit cell.
in Figure2a,b.Other parameters include the length of the unit cell the finite phononic crystal beam L n l = × , with ( ) l α the length of structural beam ( )

Figure 1 .
Figure 1.A bi-coupled periodic beam consisting of n unit cells with a defect in the { } j th unit cell.

Figure 2 .
Figure 2. The configuration of: (a) The perfect unit cell { } i ; (b) The defected unit cell { } j .
denote the ( ) j beam and J joint in the { } i unit cell, respectively.It should be emphasized that henceforth we use the { } lo w e r c a s e le tte r , ( ) lo w e r c a s e le tte r , and u p p e r c a s e le tte r to

Figure 2 .
Figure 2. The configuration of: (a) The perfect unit cell {i}; (b) The defected unit cell {j}.
by the wave amplitudes with the subscripts 2 and 3 denoting the flexural and thickness-shear modes, respectively.For examples, a JK 2 and d JK 2 are the wave amplitudes of the backward (arriving) and the forward (departing) flexural modes in the (xyz) JK coordinate system, respectively, while a JK 3 and d JK 3 are those amplitudes of the thickness-shear modes.

2 JK a and 2 JKd 3 JK a and 3 JKd
by the wave amplitudes with the subscripts 2 and 3 denoting the flexural and thickness-shear modes, respectively.For examples, are the wave amplitudes of the backward (arriving) and the forward (departing) flexural modes in the ( ) JK xyz coordinate system, respectively, while are those amplitudes of the thickness-shear modes.

Figure 3 .
Figure 3.The description of a typical structural member ( ) j : (a) The dual local coordinates; (b) The generalized displacements and forces; (c) The wave amplitudes.

.θ
depicts the global coordinates ( ) X YZ on the background with the left exterior joint 1 {1} −, any interior joint J , and the right exterior joint [ The joint label will usually be attached to joint variables as superscript to denote their pertinence.As shown in Figure4b, are the translational and rotational displacements of the typical interior joint J , respectively, which are combined to form the joint generalized displacements.corresponding transversal force and rotational moment, respectively, which are combined to form the joint generalized forces.Nevertheless, the generalized displacements/forces of the left (right) exterior joint

Figure 3 .
Figure 3.The description of a typical structural member (j): (a) The dual local coordinates; (b) The generalized displacements and forces; (c) The wave amplitudes.

Figure 4 .
Figure 4.The description of the left exterior joint 1 {1} − , a typical interior joint J and the right exterior joint [ 1] { } M n + − : (a) The global coordinates and the dual local coordinates of the relative members; (b) The generalized displacements and forces of these joints and of their relative members; (c) The wave amplitudes of the relative members.

Figure 4 .
Figure 4.The description of the left exterior joint 1 − {1}, a typical interior joint J and the right exterior joint [M + 1] − {n} : (a) The global coordinates and the dual local coordinates of the relative members; (b) The generalized displacements and forces of these joints and of their relative members; (c) The wave amplitudes of the relative members.

Figure 5 .
Figure 5.The defected phononic crystal beam with eight binary unit cells, each of which consists of a PMMA beam and an aluminum beam having identical length 80 mm : (a) The schematic of the structure; (b) The experiment specimen; (c) The PMMA beam in a perfect unit cell of the experiment specimen; (d) The PMMA beam in the defected unit cell of the experiment specimen.

Figure 5 .
Figure 5.The defected phononic crystal beam with eight binary unit cells, each of which consists of a PMMA beam and an aluminum beam having identical length 80 mm: (a) The schematic of the structure; (b) The experiment specimen; (c) The PMMA beam in a perfect unit cell of the experiment specimen; (d) The PMMA beam in the defected unit cell of the experiment specimen.

Figure 6 .
Figure 6.The transmission (dB) spectra of deflections at two ends of the defected phononic crystal beam in Figure 5: (a) Results obtained by the MRRM with 1% damping ratios, by the experiment and by the FEM with 1% isotropic structural loss factors and with 6227 elements; (b) Results obtained by the MRRM with five cases of damping ratios and their comparisons with the experiment result; (c) Undamped results obtained by the FEM with four meshing schemes and their comparisons with the MRRM result.

Crystals 2018, 8 , 21 4 of 9 Figure 7 .Figure 8 .Figure 8 .Figure 8 .
Figure 7.The transmissions and the band structures within low frequency range 2.5 kHz f ≤ : (a) Theoretical and experimental transmissions of the defected phononic crystal beam and the computed transmission of the counterpart perfect phononic crystal beam by the MRRM; (b) Band structures of the defected phononic crystal beam as a supercell by the TMM; (c) Band structures of the defected phononic crystal beam as a supercell by the MRRM; (d) Band structures of the perfect phononic crystal beam based on supercell of eight unit cells by the MRRM; (e) Band structures of the perfect phononic crystal beam based on one unit cell by the MRRM.

Figure 9 .
Figure 9.The comparison of band structures of the defected phononic crystal beam as a supercell within 2.5 kHz f ≤ by the MRRM based respectively on the Timoshenko and Euler-Bernoulli beam

Figure 9 .
Figure 9.The comparison of band structures of the defected phononic crystal beam as a supercell within f ≤ 2.5 kHz by the MRRM based respectively on the Timoshenko and Euler-Bernoulli beam theories.

Figure 10 .Figure 10 .
Figure 10.The transmission spectra and the band structures of the four phononic crystal beams with the PMMA beam in the fifth unit cell having cross-sectional height 8 mm , 12 mm , 15 mm , and 18 mm , respectively: (a) Transmission spectra within 1500 Hz f ≤ ; (b) Phase and attenuation constant spectra within 150 Hz f ≤ ; (c) Phase and attenuation constant spectra within 300 Hz 1100 Hz f ≤ ≤ .

Figure 10 .Figure 11 .Figure 11 .
Figure 10.The transmission spectra and the band structures of the four phononic crystal beams with the PMMA beam in the fifth unit cell having cross-sectional height 8 mm , 12 mm , 15 mm , and 18 mm , respectively: (a) Transmission spectra within 1500 Hz f ≤ ; (b) Phase and attenuation constant spectra within 150 Hz f ≤ ; (c) Phase and attenuation constant spectra within 300 Hz 1100 Hz f ≤ ≤ .

Figure 12 .
Figure 12.The computed transmission spectra and band structures within 1500 Hz f ≤ of the four defected phononic crystal beams (with two, four, six and eight unit cells): (a) The transmission coefficients; (b) The phase and attenuation constants spectra.

Figure 12 .
Figure 12.The computed transmission spectra and band structures within f ≤ 1500 Hz of the four defected phononic crystal beams (with two, four, six and eight unit cells): (a) The transmission coefficients; (b) The phase and attenuation constants spectra.

Table 2 .
The known and unknown quantities of four types of right exterior joint[M + 1]