Flexural Wave Band Gaps in Phononic Crystal Euler-Bernoulli Beams Using Wave Finite Element and Plane Wave Expansion Methods

We investigate theoretically and experimentally the forced response of flexural waves propagating in a 1D phononic crystal (PC) Euler-Bernoulli beam, composed by steel and polyethylene, and its band structure. The finite element, spectral element, wave finite element, wave spectral element, conventional and improved plane wave expansion methods are applied. We demonstrate that the vibration attenuation of the unit cell can be improved choosing correctly the polyethylene and steel quantities and we suggest the best percentages of these materials, considering different unit cell lengths. An experiment with a 1D PC beam is proposed and the numerical results can localize the band gap position and width close to the experimental results. A small Bragg-type band gap with low attenuation is observed between 405 Hz 720 Hz. The 1D PC beam with unit cells of steel and polyethylene presents potential application for vibration control.


Introduction
Artificial periodic composites known as phononic crystals (PCs), consisting of a periodic array of scatterers embedded in a host medium, have been quite studied [1][2][3][4][5][6][7][8][9][10][11][12][13][14] .They have received renewed attention because they exhibit band gaps where there are only mechanical (elastic or acoustic) evanescent waves.There are no mechanical propagating waves in the band gaps.This ability of creating phononic band gaps is similar to the electronic and photonic band gaps in semiconductors and photonic crystals, respectively.
The physical origin of phononic and photonic band gaps can be understood at micro-scale using the classical wave theory to describe the Bragg and Mie resonances based on the scattering of mechanical and electromagnetic waves propagating within the crystal 15 .
LRPCs, differently from the traditional PCs which create the Bragg-type band gaps, present the locally resonant (LR) band gaps.These LR band gaps can be obtained in a frequency range of orders of magnitude lower than that given by the Bragg limit.LR band gaps arise in the vicinity of the natural frequency of the local resonators while Bragg-type band gaps typically occur at wavelengths of the order of the unit cell size.The concept of an EM or LRPC generally involves the inclusion of local resonators, and the periodicity is advantageous (creates the Bragg-type band gaps), but it is not necessary in an EM 31,32 .
Most of the studies concerning PCs focused on investigation of bulk mechanical waves 1,2,[7][8][9][10]12 and its results have shown that the band gaps may appear because of the contrast between the physical properties, for instance elastic modulus and density of the inclusions and matrix. Oter important properties that influence the band gap width are the inclusion geometry, filling fraction and PC lattice.The band gaps may also be affected by the physical nature of PC, which can be: solid/solid 9,10,12 , fluid/fluid 33 and mixed solid/fluid 34 PCs.
Only few researchers used methods developed recently 46,48 .Junyi et al. 46 developed an inverse method to calculate the band structure of one dimensional periodic structures based on Bloch wave boundary conditions and wave superposition, Materials Research whereas Nobrega & Dos Santos 48 used the wave finite element (WFE) 52 and wave spectral element (WSE) methods.
In this context, the main purpose of this study is to investigate the Bragg-type band gap formation, band structure, also known as dispersion relation, and attenuation constant of a 1D PC beam using the FE, SE, WFE, WSE, conventional plane wave expansion (CPWE) and improved plane wave expansion (IPWE) methods.Furthermore, we demonstrate the accuracy and efficiency of the methods cited for modeling PCs.We also compare the analytical and numerical results with an experimental validation.A small band gap with low attenuation is observed between 405 Hz -720 Hz for the real 1D PC beam.
To the best of our knowledge, it is the first time the band structure of a PC is compared using both CPWE, IPWE, WFE and WSE methods.Only few studies 48 have used the advantages of the recent WFE and WSE methods to study PCs.In addition, the attenuation constant surface of a 1D PC unit cell is obtained from varying the percentages of polyethylene and steel.

1D Phononic Crystal Euler-Bernoulli Beam Model by FE, SE, WFE, WSE, CPWE and IPWE Methods
This section presents the formulation for a 1D PC Euler-Bernoulli beam using the FE, SE, WFE, WSE, CPWE and IPWE methods.FE and SE are methods well stablished in literature 49,50 .SE is an analytical method and has solution only for simple structures, such as rods, beams and plates and shells for specified boundary conditions, while the FE is a numerical method and can handle complex structures.WFE and WSE are methods that present some special advantages, such as that the discretization is only on the PC unit cell.
CPWE method, also known as ω(k) method, is a semianalytical method used to predict the band structure.CPWE method presents a slow convergence, mainly for systems with a large property mismatching.To solve this convergence problem, we use the IPWE method.Cao et al. 53 proposed the IPWE method to handle PCs and showed that this method provides much more accurate numerical results than the CPWE.

1D Phononic crystal beam model
Figure 1 sketches a 1D PC beam with a periodic array of unit cells containing two different materials, i.e. steel (blue) and polyethylene (white).The PC beam lattice parameter is a.Each unit cell is composed by 2/3 of steel and 1/3 of polyethylene, Figure 1 (a), and 2/3 of polyethylene and 1/3 of steel, Figure 1 (c), called as model 1 and model 2, respectively.It is important to mention that we consider the Euler-Bernoulli (EB) beam theory.

SE and FE methods for 1D PC EB beams
We calculate the dynamic stiffness matrix using two approaches, i.e. one analytical, SE method and the other numerical, FE method.The dynamic stiffness matrix of EB beam element modeled by SE method 49 is: (1) where where k is the Bloch wave vector, also known as wavenumber, E is the Young's modulus, ρ is the density, I is the second moment of area, L se is the length of the spectral element, i 1 = -, S is the cross section area of the beam and ω is the angular frequency.If the geometry of the beam is uniform and regarding just one material, the global dynamic stiffness matrix of the EB beam modeled by SE method can be considered as equal to the dynamic stiffness matrix of the EB beam element, i.e.D beamse = D e beamse .However, we discretize each unit cell in three spectral elements, one element for each part of the unit cell in Figure 1 (b) and (d).Thus, the global dynamic stiffness matrix can be obtained by the assembly of the dynamic stiffness matrices of the EB beam elements modeled by SE method.
The dynamic stiffness matrix can also be obtained by the FE method, D e beamfe = K e beamfe -ω 2 M e beamfe , where 50 : where L fe is the length of the finite element.The global dynamic stiffness matrix is given by where K beamfe and M beamfe are the global stiffness and mass

WSE and WFE methods for 1D PC EB beams
Considering a finite EB beam divided into unit cells (the unit cells are meshed with an equal number of nodes on their left-and right-hand edges), see Figure 1, one can obtain a dynamic stiffness of the unit cell, Dcell u , from FE, Dcellfe u , or from SE, Dcellse u , methods.The Dcell u can be portioned and the dynamic equation of motion can be written as: (10)   where u is the vector of the displacement degrees of freedom, p is the vector of the applied forces, i represents the interior degrees of freedom of the unit cell and l and r represent the left and right boundaries of the unit cell, respectively.The interior degrees of freedom can be eliminated using the first row of Eq. ( 10) 54 , which results in: (11)   Inserting Eq. ( 11) into Eq.( 10) leads to: (12)   where Equation (12), which relates the forces and the displacements on the two sides of the unit cell, is the starting point for the WFE 54 and WSE analysis.Note that D cell can be obtained by FE method, D cellfe , or by SE method, D cellse .Thus, this is the difference between the WFE and WSE methods.
One can reformulate Eq. ( 12) in terms of state vectors as: where T is a symplectic 55 matrix and it is known as transfer matrix, q r is the right state vector and q l is the left state vector.The continuity of displacements and equilibrium forces at the boundary between cells m and m + 1, yields, Substituting Eqs. ( 14) and (15) in Eq. ( 13), one can write: Invoking Floquet-Bloch's theorem 56,57 , say: where μ = -ika is the attenuation constant of the unit cell.Substituting Eq. ( 17) in Eq. ( 16) leads to: As pointed out by Zhong & Williams 55 , direct calculation of the eigenvalue problem in Eq. ( 18) can lead to numerical ill-conditioning.According to Mencik 52 , this ill-conditioning can be explained by the Bauer-Fike theorem 58 , which says that the problem is that the eigenvector of T can be illconditioned.This can be explained as it is partitioned into displacement and force components, thus the values can be largely disparate.To solve this issue, Zhong and Williams 55 proposed (for symmetric unit cells) a homogeneous generalized eigenvalue problem of the following form: where I n is the n-dimensional unit matrix (identity matrix) and substituting Eqs.(19) and (20) in Eq. ( 13), we can note that T = NL -1 .Thus, rewriting Eq. ( 18): The solutions of Eq. ( 21) are denoted as {(e μj , w j )} j , where {(e μj )} j are the eigenvalues and {w j } j are the eigenvectors associated with a given mode j.One can notice that each eigenvector of Eq. ( 18), and using Eq. ( 19), can be obtained by q lj = Lw j .It can also be shown that {(e -μj , (JLw) j )} j are also eigenvalues of Eq. ( 21) and left eigenvectors, respectively, because one can rewrite Eq. ( 21) as e μ N T (JLw) = L T (JLw) ∴ (JLw) T (N -e -μ L) = 0 55 , where: Therefore, the 2n eigenvalues obtained from Eq. ( 21), when ordered appropriately, can be subdivided into two groups.The first corresponds to the waves travelling to the right, e μj , j = 1,2,...,n with |e μj | < 1 and the second corresponds to the waves travelling to the left, e -μj , j = 1,2,...,n with |e -μj | > 1.Note that to consider these hypotheses the substructure needs to have damping.If the substructure does not have damping, |e ±μj | = 1.
For the eigenvalues {e μj } j and {e -μj } j , one can obtain the associated eigenvectors {q lj } j , and {q * lj } j , respectively.Each eigenvector can be split into displacement and force   , where R is the diagonal symmetry transformation matrix.Tracking the frequency evolution of each wave mode is a crucial step of the WFE and WSE methods, for more details see Mencik 52 .
After some mathematical manipulation 52 , to predict the harmonic response of a structure composed of N unit cells, we can write: (23) (24)   where q l (1) is the state vector of the structure first node, q r (N) is the state vector of the structure last node and q is the vector form of the modal amplitude.Equations (23-24)  and the calculation of the modal amplitudes are discussed by Mencik 52 .

CPWE and IPWE methods for 1D PC EB beams
From EB beam theory, the governing equation for flexural vibration of a uniform beam system can be written as: (25)   where u is the transversal displacement.
Considering the 1D PC beam illustrated in Figure 1, however, an infinite beam for this formulation, we can note that E = E(x), I = I(x), ρ = ρ(x) and S is constant.In order to eliminate the factor time in Eq. ( 25), we apply the temporal Fourier transform.
Applying the Floquet-Bloch's theorem 56,57 , expanding u as a Fourier series, and considering wave propagation on the x axis (k 2 = k 3 = 0), we can write: m ∈ ℤ.Note that g is a one-dimensional vector, because we consider one-dimensional periodicity.Furthermore, we may write: where P is one of α = ρS, β = EI and g has the same expressions of g, with m ∈ ℤ.Note that we use g to highlight the difference between the expansions of material properties and the displacement.For CPWE method, the material parameters are all directly expanded in Fourier series according to the spatial periodicity.However, in the IPWE method, the inverses of the material parameters (excluding the mass density) are expanded in Fourier series in order to get a good convergence 53 , thus P is one of the α and β for the CPWE and one of the α and 1/β for the IPWE.Substituting Eqs. ( 26) and ( 27) in Eq. ( 25), with g g g = + u r , multiplying by e -ig.x and integrating over the unit cell, one may write (for all x) for CPWE and IPWE methods, respectively: Equations ( 28) and ( 29) represent a generalized eigenvalue problem of ω 2 (k) and should be solved for each k into the FIBZ.
The Fourier coefficients are: where the indexes A and B of Eq. ( 30) are related to the inclusion (steel) and the matrix (polyethylene), respectively, considering model 1, Figure 1 (a), or can also be related to the inclusion (polyethylene) and the matrix (steel), respectively, considering model 2, Figure 1 (c).F(g) is the structure function and f a LA = r is the filling fraction, where L A is the inclusion length.The structure function F(g) is defined as: The integral in Eq. ( 31) is performed over the unit cell length and results in: (32)   where g = ||g||.

Numerical validation
The PC EB beam parameters and material properties are summarized in Table 1, where the subscripts A and B, as mentioned before, refer to steel and polyethylene, considering , for for g g g g P fP f P P P F 1 (a), and polyethylene and steel, considering model 2, Figure 1 (c), respectively.Note that the structural damping, η A , η B , also known as loss factors, are included as a complex Young's modulus, It is important to mention that we limited the most part of the results analysis until 10240 Hz.

Method comparison
The forced response of the PC beam is analyzed considering a free-free boundary condition and an excitation force as a cosine-shaped pulse only on the left side of the beam.Figure 2 shows the PC beam displacement of left (first node) and right (last node) sides for the models 1 and 2.
We can see in Figure 2 some regions where the resonances do not appear, however, it is difficult to localize exactly the band gaps.To overcome this limitation, we plot the frequency response function (FRF) in Figure 3 and the transmittance in Figure 4, defined as the division between the displacements of the last and first nodes.For the FRF, we choose the receptance, i.e. the division between the displacement of the first or the last nodes, and the force, which gives H 11 or H 21 , also known as point receptance and transfer receptance, respectively.
From Figures 3 (a-b) and 4 (a), that is to say model 1, we can see that the band gap is opened up between 2780 Hz and 5798 Hz.This band gap is known as a Bragg-type band gap, because the mechanism involved is the Bragg scattering.Thus, the frequency location is governed by the Bragg's law, a = n(λ/2), (n ∈ ℕ * >0 ), where λ is the wavelength of waves in the host material.The Bragg's law implies that it is difficult to achieve a low frequency Bragg-type band gap in PCs with small size.For model 2, Figures 3 (c-d) and 4 (b), it is not possible to identify clearly the band gap formation, because its low attenuation, as will be discussed in band structure results.

Influence of unit cell length
In Figure 5 (a-b), we compare the transmittance for different lattice parameters calculated by WSE method, considering models 1 and 2, respectively.We consider a = 0.212, 0.106, 0.0707, 0.053, which results in N = 2, 4, 6, 8, for a fixed beam length L = 0.424 in Table 1.We choose the WSE method because it is an analytical method.
For models 1 and 2, and the first case, a = 0.212, N = 2, we can see in Figure 5        In Figure 5 (c) and (d), we show the transmittance for N = 2, calculated by WSE, WFE, FE and SE methods for models 1 and 2, respectively.We can see that the WFE and FE methods do not converge with WSE and SE methods in high frequencies.This happens because the discretization of the unit cell (two finite elements for each part of the unit cell, as mentioned before) using the FE method is not enough in higher frequencies.

Complex elastic band structure
Figure 6 illustrates the real and imaginary parts of the elastic band structure, considering the data in Table 1, N = 10, for models 1 (a-b) and 2 (c-d).Figure 6 (a-c) shows the real part of the reduced Bloch wave vector (dimensionless), ka/π, using WSE, WFE and IPWE methods and Figure 6  (b-d) shows the imaginary part of the reduced Bloch wave vector using WSE and WFE methods.
The same Bragg-type band gap observed in Figures 3  (a-b) and 4 (a) for model 1, between 2780 Hz and 5798 Hz, can be observed in Figure 6 (a-b).For model 2, the Braggtype band gap is not observed in Figures 3 (c-d) and 4 (b), as mentioned before, however, it exists and it is opened up between 2643 Hz and 3175 Hz, as we can see in Figure 6 (c-d).
The unit cell attenuation constant, μ, is an important information that can be analyzed from the imaginary part of the Bloch wave vector, Figure 6 (b-d), because they are related as ℜ(μ) = ℑ(k)a for a complex k.In Figure 6 (b-d), we can see that the attenuation performance of the Bragg-type band gaps is better for the model 1, because the existence of a broader band gap with higher attenuation.
In IPWE and CPWE calculations, we consider 101 plane waves in the Fourier series expansion and we only show the FIBZ 59 , [0, π/a], in Figure 6 (a) and (c).It can be seen a good matching from WSE, WFE and IPWE methods in Figure 6 (a) and (c).
Figure 7 (a-b) shows the comparison between the IPWE and CPWE methods inside the FBZ 59 , i.e. [-π/a, π/a], for models 1 and 2, respectively.In Figure 7, only the first 10 branches are illustrated.The matching between IPWE and CPWE does not occur only for the higher branches, even considering a high number of planes waves.
In order to demonstrate the matching between WSE and IPWE methods for other values of N, we plot in Figure 8 (a) and (c) the real part of reduced Bloch wave vector for N = 2, 4, 6, 8, 10.It can be observed a good matching between WSE and IPWE methods.In Figure 8 (b) and (d), it is presented the imaginary part of reduced Bloch wave vector for N = 2, 4, 6, 8, 10 and its behaviour is similar to the discussed in Figure 5 (a-b), however, we also have the information of the unit cell attenuation performance and not only the band gap width.In Figure 8 (b-d), the best attenuation performance for N = 4, 6, 8, 10 is found for model 1, because more broad band gaps with higher attenuation are created.However, for N = 2, the best attenuation performance between models 1 and 2 is not clear, because both models present many band gaps with high attenuation.

Influence of polyethylene and steel quantities on unit cell attenuation
Figure 9 (a-j) shows the influence of the polyethylene quantity (5% to 95%), model 1 (a,c,e,g,i), and steel quantity (5% to 95%), model 2 (b,d,f,h,j), on the unit cell attenuation for N = 2 (a-b), N = 4 (c-d), N = 6 (e-f), N = 8 (g-h) and N = 10 (i-j), respectively, calculated by WSE method.It is important to mention that until Figure 8, we have only considered 1/3 (≈ 33,33%) of polyethylene for model 1, and 1/3 (≈ 33,33%) of steel for model 2, as illustrated in Figure 1.The influence of the polyethylene and steel quantities on the unit cell attenuation performance is complex and depends on the unit cell length.
For model 1, no attenuation is observed for N = 2, 4, 6, 8, 10 until 98 Hz, 406 Hz, 874 Hz, 1612 Hz and 2500 Hz, see Figure 9 (a,c,e,g,i), respectively, independently of the polyethylene quantity.Whereas for model 2, no attenuation is observed for N = 2, 4, 6, 8, 10 until 98.75 Hz, 402.5 Hz, 910.6 Hz, 1593 Hz and 2488 Hz, see Figure 9 (b,d,f,h,j), respectively, independently of the steel quantity.Thus, depending on the application, choosing the polyethylene and steel quantities, on models 1 and 2, correctly it is not an easy task, because it is related to the unit cell length and in which frequency the band gap appears.

Experimental Validation
A 1D PC beam is used to perform an experimental test.The PC beam is similar to the model 1 proposed in Figure 1 (a-b) in a free-free boundary condition, however, a B = 0.041 m and a A = 0.0325 m, with a = 2a A + a B = 0.106 m.The properties are the same described in Table 1, with N = 4.
The measurement instruments used in the experimental setup are summarized in Table 2. Figure 10 shows the experimental setup with the details of the impact hammer and accelerometer position.By using an impact force excitation applied to the right and left ends of the PC beam, acceleration measurements are taken on the right end of the PC beam.
However, we choose to plot the displacement, that is , where a c is the acceleration measured.Inertance point and transfer FRFs are measured with 5 averages, with the frequency discretization of 0.625 Hz. Figure 11 (a-d) illustrates the displacement of the last beam node (right side), the transmittance and the FRFs H 11 and H 21 , respectively.The numerical results present good agreement with the experimental results, however, FE and WFE methods do not match in higher frequencies with the analytical methods, as discussed before.Furthermore, there is some mismatch related to the experimental results.The numerical band gap widths do not match exactly with the experimental band gap widths, as expected, because our numerical model may not capture all the real aspects of the PC beam, such as the material used to glue the polymer and metal, the material properties may not be exactly the same of Table 1, among others.In addition, we use the EB beam theory, perhaps considering higher beam theories, such as Timoshenko beam theory 60 , the results may be improved in high frequencies.
Figure 12 shows the complex elastic band structure of the real PC beam, illustrated in Figure 10, with the Bragg-type band gaps.The band gap widths observed in Figure 11 can be confirmed in Figure 12.Furthermore, we may observe in Figure 12 a small Bragg-type band gap between 405 Hz -720 Hz with low attenuation.

Conclusion
We obtain the forced response and the complex elastic band structure of a 1D PC beam proposed by models 1 and 2. The forced response is obtained by the WFE, WSE, FE and SE methods and a good matching is observed, except for high frequencies, where WFE and FE do not match with the spectral analytical methods.The real part of the Bloch wave vector is calculated by the WFE, WSE and IPWE methods and it is shown a good agreement.The CPWE method is not accurate compared to the IPWE method for the higher bands.
The influence of the unit cell length is also studied and for longer unit cells, Bragg-type band gaps are opened up in low frequencies.The polyethylene or steel quantity into the unit cell, considering models 1 and 2, respectively, is an important variable and its influence in the unit cell attenuation constant it is complex and depends of the unit cell length.We show some ranges of frequency, considering fixed unit cell lengths, which no unit cell attenuation is observed independently of the polyethylene and steel quantities.We also show the regions that present higher attenuation.In a general way, the model 1 proposed presents better unit cell attenuation behaviour than model 2.
The analytical and numerical results present a good agreement with the experimental results and they can localize the band gap position and width close to the experimental.A small Bragg-type band gap with low attenuation is observed between 405 Hz -720 Hz.The 1D PC beam with unit cells of steel and polyethylene presents potential application for vibration management.
S S S S S S S S V X W W W W W W W W W W W matrices, respectively.They are obtained by assembling the stiffness and mass matrices of the EB beam elements modeled by FE method.Initially, we discretize each unit cell in six finite elements, two elements for each part of the unit cell in Figure1(b) and (d).

=
important relation between the components of each eigenvector is u * lj Ru lj and q * lj = -Rq lj 52

( 26 )
where x = xe 1 , u k (x) is the amplitude of the Bloch wave, note that u k (x) = u k (x + r) and u(x + r) = u(x)e ik.r , e ik.r is called the Floquet-Bloch periodic boundary condition, k = vb 1 , v ∈ ℚ is the symmetry point within the first irreducible Brillouin zone (FIBZ)59 in reciprocal space, , a 0 r # & , or one may write k = k 1 e 1 , k 1 ∈ ℝ is the point coordinate within the FIBZ for the reciprocal space, b 1 is the basis vector in the reciprocal space defined as b 1 = a 2r e 1 , a 1 = ae 1 is the component of the lattice vector r = pa 1 , p ∈ ℤ.The reciprocal lattice vector is defined as g = a 2r me 1 , (a) and (b) four and six Bragg-type band gaps, respectively, and the first band gap widths are 122.5 Hz -456.9Hz and 103.8 Hz -128.8Hz, respectively.For the other cases, a = 0.106, 0.0707, 0.053, N = 4, 6, 8, Figure 5 (a) and (b), the first band gap widths are 456.9Hz -1204 Hz , 1008 Hz -2325 Hz and 1782 Hz -3849, for

Figure 1 .
Figure 1.Schematic representation of a PC beam with N unit cells of steel (blue) and polyethylene (white) considering models 1 (a) and 2 (c).The PC beam unit cells are illustrated in (b) and (d).

Figure 2 .
Figure 2. PC beam displacement of the left (a-c) side, first node, and right (b-d) side, last node, calculated by WSE, WFE, FE and SE methods, considering models 1 (a-b) and 2 (c-d).

Figure 3 .
Figure 3. Point receptance (a-c) and transfer receptance (b-d) of the PC beam calculated by WSE, WFE, FE and SE methods, considering models 1 (a-b) and 2 (c-d).

Figure 4 .
Figure 4. Transmittance of the PC beam calculated by WSE, WFE, FE and SE methods, considering models 1 (a) and 2 (b).

Figure 6 .
Figure 6.Complex elastic band structure of the PC beam considering N = 10 for models 1 (a-b) and 2 (c-d).The real part of the reduced Bloch wave vector (dimensionless) (a-c) is calculated by WSE, WFE and IPWE methods and the imaginary part of the reduced Bloch wave vector is calculated by WSE and WFE methods.

Figure 7 .
Figure 7. Real part of the reduced Bloch wave vector (dimensionless) calculated by IPWE and CPWE methods for models 1 (a) and 2 (b).

Figure 8 .
Figure 8. Complex elastic band structure of the PC beam considering N = 2, 4, 6, 8, 10 for models 1 (a-b) and 2 (c-d).The real part of the reduced Bloch wave vector (dimensionless) (a-c) is calculated by WSE and IPWE methods and the imaginary part of the reduced Bloch wave vector is calculated by WSE method (b-d).

Figure 12 .
Figure 12.Complex elastic band structure of the real PC beam.The real part of the reduced Bloch wave vector (dimensionless) (a) calculated by WSE, WFE and IPWE methods and the imaginary part of reduced Bloch wave vector calculated by WSE and WFE methods.

Table 1 .
Beam geometric parameters and material properties.
Figure 10.Experimental setup of the PC beam.