Modeling Bloch Waves in Prestressed Phononic Crystal Plates

The study aims at investigating the effect of a generic state of prestress on the passbands and bandgaps of a phononic crystal plate. To this end, an Updated Lagrangian scheme is developed, consisting in a two-step procedure: ﬁrst, a static geometrically nonlinear analysis of a representative unit cell undergoing the action of an applied external load is conducted and then the Floquet-Bloch decomposition is applied to the linearized equations of the acousto-elasticity for the unit cell in the deformed conﬁguration. In addition, a formula for the calculation of the energy velocity is proposed. In the case of an epoxy plate with cylindrical steel inclusions, it is shown that, even in the presence of prestress inducing full reversible deformation state, the bandgap experiences a shift towards higher frequencies when the cell is subjected to a compressive prestress, whereas a frequency downshift is registered when the cell is subjected to traction. In particular, it is demonstrated that the frequency downshift of the bandgap for the phononic plate undergoing a tensile prestress is approximately 3.5% with respect to the case of the phononic plate under compression. The results presented herein provide insights in the behavior of phononic crystal plates with tunable dispersive properties, and suggest new leverages for wave manipulation valuable in many application ﬁelds such as wave ﬁlters, waveguiding and beam splitting, sensing devices, and vibration shielding.

However, most of the configurations proposed so far are limited by the fact that their unconventional dynamic properties are accessible to specific frequencies only, once the unit cell design is conceived, being its geometrical and mechanical properties constant with respect to time. To overcome this limitation, periodic systems with adaptive elastic properties have received great attention in the recent years. Indeed, introducing an additional degree of freedom, enables the reversible tailoring of the dispersive diagram of the structure, considerably enhancing their functionality and technological applications. To achieve this goal, several approaches have been proposed and explored so far. They include active systems, such as piezoelectric materials (Casadei et al., 2010;Bergamini et al., 2014;Kherraz et al., 2016), temperature (Jim et al., 2009;Airoldi et al., 2011;Cheng et al., 2011;Wu et al., 2018) and magneto-based techniques (Robillard et al., 2009;Matar et al., 2012;Guo and Wei, 2016;, actuated polymers via electromagnetic waves (Walker et al., 2014(Walker et al., , 2017 as well as the application of an external mechanical load. Among them, the latter has shown to be particularly promising to achieve tunable dispersive properties. Bigoni et al. (2008) formulated a theoretical model for an orthotropic, prestressed (compressible) elastic layer vibrating on an elastic half space assuming longwave asymptotics for the solution. It was found that the influence of the prestress over bandgaps and passbands, so that it could be exploited as a tuning parameter to shift the dispersion curves of the system. A modeling tool for the prediction of controlled bandgap structures responding to flexural vibrations was provided by Gei et al. (2009), also relaxing the hypothesis of perfect periodicity, allowing thus to consider quasi-periodic structures (Gei, 2010). Wang and Bertoldi (2012) showed how mechanical deformations enabled the tuning of the phononic bandgaps in 3D periodic elastomeric structures. The ability of the elastomers to undergo small as well as large strain deformations guaranteed the reversibility and repeatability of the process. Both the linear and nonlinear regimes of elastic deformation were explored, including different geometrical topology for triggering mechanical instability. Galich et al. (2016) showed how local deformations experienced by a composite structure can bring to local stiffening influencing the instability-induced interfaces on elastic wave propagation in finitely deformed layered materials. Finally, Norris and Parnell (2012) tried to exploit the effects of the prestress to theoretically show the feasibility of a hyperelastic cloaking via transformation elasticity.
In this paper, the effects of an applied mechanical load on the dispersive diagram for phononic crystal plates is investigated in the case of elastic deformations, so to have a complete reversibility of the phenomena. The paper is organized as follows: section 2 provides the description of an Updated Lagrangian scheme, in which a representative unit cell of the phononic crystal plate is studied in its static and dynamic deformed configurations. In section 3, the procedure to extract the band diagrams and in particular the energy velocity using the equations described in section 2 is detailed. Finally, section 4 investigates the effect of a reversible compressive and tensile state of prestress on an epoxy plate with cylindrical inclusions.

Static Analysis
In what follows, the equations of static and dynamic equilibrium are derived with reference to the generic infinite phononic crystal (PC) shown in Figure 1, where C 0 indicates the static undeformed configuration, C the static configuration resulting from the application of an external load, and C ′ the dynamic configuration undergoing a harmonic motion. The so called Updated Lagrangian (UL) scheme is employed to analyze the PC in the three different equilibrium configurations. The UL approach consists of two steps: in the first step, C 0 is used as the reference configuration to calculate the displacement and stress fields relative to the configuration C, while in the second step, C is used as the new reference configuration to evaluate the dynamic equilibrium of the unit cell in the configuration C ′ .
Due to the periodicity of the system under consideration, in the undeformed configuration C 0 any scalar, vector or tensor function φ satisfies the condition φ(x 0 + r 0 m) = φ(x 0 ), where m ∈ Z 2 and r 0 = [r 01 , r 02 ] is the matrix of lattice vectors. This allows to restrict the computational domain to a unit cell of domain 0 and boundary ∂ 0 (Collet et al., 2011;Mazzotti et al., 2017). Similarly, the unit cell in the configuration C satisfies the periodic condition φ(x + rm) = φ(x), with r = [r 1 , r 2 ]. For this cell, the mapping between the coordinates x 0 = {x 01 , x 02 , x 03 } T of a material point in C 0 and the coordinates x = {x 1 , x 2 , x 3 } T of the same point in C is established by the deformation gradient F(x 0 ) = ∇ 0 x = I + ∇ 0 u 0 , where ∇ 0 denotes the gradient operator defined with respect to C 0 , u 0 is the pseudo-static displacement resulting from the application of external, 0 -periodic, volume (f V ) and surface (t S ) loads. Following Zhang and Parnell (2017), the deformation gradient can be decomposed as F(x 0 ) = F L (x 0 )F P (x 0 ), where F L indicates the affine deformation gradient of the lattice points such that r = F L r 0 , while F P denotes a periodic non-affine deformation.
The specific material density of the unit cell in C 0 is ρ 0 (x 0 ), while, assuming a hyperelastic material behavior described by the Murnaghan's model (Murnaghan et al., 1937;Pau and Lanza di Scalea, 2015;Dubuc et al., 2017Dubuc et al., , 2018, the tensor of tangential elastic moduli with respect to C 0 is expressed by D 0 = 4∂ 2 /(∂C∂C), where the elastic energy density has the form in which λ and µ denote the first and second Lamé parameters, respectively, (l, m, n) the third order Murnaghan parameters, C = F T F the right Cauchy-Green deformation tensor, E = 1 2 (F T F − I) the Green-Lagrange strain tensor and I 1 (E), I 2 (E), and I 3 (E) its first, second and third invariants, respectively.
Finding the coordinates of the deformed configuration C and the associated stress fields requires the solution of a variational functional of the form (Bonet and Wood, 2008) (2) subjected to the Dirichlet boundary conditions in which δu 0 is an arbitrary admissible kinematic variation of u 0 (x 0 ), S = D 0 : E the second Piola-Kirchhoff stress tensor and n(x 0 ) the outward pointing normal at x 0 ∈ ∂ 0 . The application of a standard Galerkin approach to Equation (2) results in the generalized system of equations where K(Q 0 ) is a static stiffness matrix, P 0 represents the global vector of nodal forces, Q 0 denotes a global vector of independent nodal displacements and Ŵ 0 is a mapping operator resulting from Equation (3) and realizing the condition U 0 = Ŵ 0 Q 0 , in which U 0 indicates the full vector of nodal displacements. In this work, the solution of Equation (4) is carried out using the weak form module implemented in Comsol Multiphysics 5.3 (Comsol, 2017).

Geometry and Mechanical Properties Updating
Once the set of independent static displacements Q 0 is obtained, the reference configuration is updated from C 0 to C by calculating the corresponding nodal coordinates x = x 0 + Ŵ 0 Q 0 (x 0 ). The updated material properties in C are given by ρ = ρ 0 (det F) −1 and D ijkl = (detF) −1 F iI F jJ F kK F lL (D 0 ) IJKL , while the Cauchy stress tensor is obtained from the relation σ = (det F) −1 FSF T . The geometry of the unit cell in the configuration C is then remeshed and used as the basis for the linear dynamic analysis presented in the next section.

Dynamic Analysis Using the Floquet-Bloch Decomposition
Assuming C as the new reference configuration, the position vector for the unit cell in the dynamic deformed configuration C ′ is given by is a time-harmonic perturbation superimposed on C, being t the time and ω the angular frequency. In the rest of the paper, the timedependence exp(−iωt) is dropped for conciseness. According to the small-on-large displacement hypotheses usually assumed in acousto-elasticity (Mazzotti et al., 2012;Pau and Lanza di Scalea, 2015;Shim et al., 2015), only small perturbations u(x, t) are considered. In this case, the approximations x ′ ≈ x and C ′ ≈ C hold and any displacement-dependent vector and tensor field can be obtained in a linearized incremental form. Following this procedure, the stress tensor is linearized in the direction of u(x, t) by applying a first order Taylor series expansion of σ (x) about x, resulting in (Mazzotti et al., 2012) where e(x) = 1 2 [∇u(x) + (∇u(x)) T ] denotes the linearized Green-Lagrange strain tensor.
From the application of the Floquet-Bloch theorem, any small harmonic perturbation u(x) can be expressed in the form (Collet et al., 2011) in whichũ(x) is a -periodic displacement amplitude and k ∈ is the Bloch wavenumber vector, being the reciprocal unit cell defined in C by the reciprocal lattice vector basis g j , j = 1, 2, satisfying the condition r i · g j = 2πδ ij , where δ ij is the Kronecker delta. These are related to the reciprocal lattice vector basis g 0j in C 0 by g j = F −T L g 0j , j = 1, 2 (Zhang and Parnell, 2017). In the reciprocal lattice domain, the Bloch wavenumber vector is described in terms of its orientation angle ϑ as k(ϑ) = k (ϑ), where ϑ is defined with respect to g 1 , k = k(ϑ) 2 and (ϑ) = {cosϑ, sinϑ, 0} T . The relation between the orientation vector in C 0 and C writes = F −T L 0 . From Equation (6), and by defining the k-shifted gradient of a generic -periodic vector fieldφ(x) ∈ C 3 as the solution of the elastodynamic problem for free vibrations of the unit cell in C subjected to an initial stress σ 0 can be obtained from the variational statement subjected to the Dirichlet boundary conditioñ in which (·) * stands for the conjugate of a complex vector or tensor field,ẽ k (x, ϑ) = 1 2 [∇ kũ (x) + (∇ kũ (x)) T ] follows from Equation (7) and n(x) represents the outpointing surface normal at x ∈ ∂ . After the application of a Galerkin discretization scheme, Equations (8) and (9) lead to the following generalized linear eigenvalue problem which forms the basis of the band structure analysis for the prestressed PC. In Equation (10), Ŵ is a mapping operator implementing the Dirichlet boundary condition in Equation (9) such thatŨ(ϑ) = ŴQ(ϑ) is verified, whereŨ(ϑ) is the global vector of nodal displacement amplitudes andQ(ϑ) a subvector ofŨ(ϑ) collecting only its independent components. The mass operator M and the stiffness operators K 1 , K 2 , and K 3 are given by where (e) denotes the domain of the e-th finite element of the mesh, e (·) stands for the standard direct stiffness assembling procedure, N(x) is a matrix of shape functions for the e-th element, 0 (x) is a block-diagonal matrix of the form while the different compatibility operators are expressed as being z i a unit vector identifying the i-th coordinate in the Cartesian frame of reference. The operators L i and L 0i are given in the Appendix. It should be noted that, since the mass and stiffness matrices in Equations (11)-(14) are evaluated with respect to the deformed configuration C, a new mesh for the deformed geometry of the unit cell needs to be generated after the static analysis described in section 2.1 has been completed and before the eigenvalue analysis is carried out.

BAND STRUCTURE ANALYSIS AND ENERGY VELOCITY EXTRACTION
The homogeneous problem in Equation (10) can be solved in the wavenumber k(ϑ, ω) and corresponding Floquet eigenvectors Q(ϑ, ω) for any fixed direction ϑ and circular frequency ω ∈ R, ω > 0, from which the band structures of the propagation and attenuation constants of a specific mode m are obtained by taking the real and imaginary components of the corresponding wavenumber k m (ϑ, ω), respectively. For the particular case of a lossless structure (Im(D) = 0) immersed in vacuum, the modes supported by the crystal can be either propagative (Im(k) = 0) or evanescent (Re(k) = 0), the latter belonging to the so called deaf frequency range (bandgap). In addition to the Bloch wavenumber k(ϑ, ω), from the computed set of solutions (k m (ϑ, ω), Q m (ϑ, ω)) it is possible to extract the band structure of the energy velocity, which corresponds to the velocity of propagation of packets of waves having close central frequency (Brillouin, 1953;Mazzotti et al., 2012). The energy velocity of a specific Bloch mode can be found as the ratio between the time-averaged energy flux in the direction ϑ over one period T = 2π/ω and the time-averaged total mechanical energy over the same period, i.e.
where φ = ω 2π t+2π/ω t φdt denotes the time-averaging operation, I (ϑ, ω) denotes the total energy flux along the orientation ϑ of the Bloch wavenumber k, K (ϑ, ω) indicates the total kinetic energy while W (ϑ, ω) and W σ 0 (ϑ, ω) represent the total stored elastic energy related to the harmonic motion and prestress, respectively. The time-averaged expressions for these quantities are given by Frontiers in Materials | www.frontiersin.org Equation (18) can be evaluated at any given solution k m (ϑ, ω) by means of a Gauss quadrature scheme over the finite element mesh of the unit cell in the configuration C, in which the displacement, strain and stress fields can be post-processed from Q m (ϑ, ω) using nodal interpolations.

NUMERICAL APPLICATIONS: EPOXY PLATE WITH CYLINDRICAL STEEL INCLUSIONS
The numerical method presented in the previous section is here applied to study a PC plate made of steel cylinders embedded in an epoxy matrix. In the undeformed configuration, the plate is 2.5 mm thick, while the cylindrical inclusions have a radius of 3.0 mm and are arranged in a square lattice of 10.0 mm side length, as shown in Figure 2A. The material properties for the steel and the epoxy are reported in Table 1.
To study how an initial state of stress affects the passbands and bandgaps of the PC plate, two different deformed configuration are considered. In the first configuration, the plate is subjected to a state of compressive stress representative of a normal displacement u 0 (x 0 ) · n 0 (x 0 ) = −0.025 mm applied at each point x 0 belonging to a lateral face of the unit cell. In the second case, a generic tensile state is applied by assuming a normal displacement u 0 (x 0 ) · n 0 (x 0 ) = 0.025 mm for each point x 0 belonging to a lateral face of the unit cell. Figures 2B,C report the von Mises stress consequent to the compression and tension prestress condition, respectively. It is worth noticing that the maximum magnitude of the applied prestretch has been limited to the typical strength at yield of both epoxy and steel in order to ensure the elastic behavior of the materials. This choice enables a full reversibility of the undeformed configuration, and thus of the dispersive behavior of the PC plate, once the load is removed, which is of more practical interest with respect to the case of a permanently induced deformation.
The band diagrams in terms of wavenumber k and energy velocity c e versus frequency are shown in Figure 3. In these diagrams, the blue and red dots denote the dispersion curves relative to the unit cell subjected to the compressive stress and tensile stress shown in Figures 2B,C, respectively. It should be noted that, since the deformation due to prestress applied to each face of the cell is isotropic in the xy-plane in both cases, the deformation gradient F and its affine component F L are diagonal and, as a consequence, the orientation of the reciprocal lattice vectors g j in the deformed configurations does not change with respect to that in the deformed configuration. This implies that the orientation angle θ of the Bloch wavevector also remains unchanged between the undeformed and deformed configurations. The results presented in Figure 3 have been obtained for an angle θ = 0 measured with respect to g 01 and g 1 . In classical Bloch analysis, this angle corresponds to the Ŵ − X direction of the irreducible Brillouin zone.
From the dispersion analysis of the undeformed unit cell, it can be found that, at frequencies below 110 kHz, only one bandgap is present along the direction ϑ = 0, which is located in the [68.3-88.6] kHz frequency range (see Figure A1 in Appendix B). This result, which has not been reported in the Figure 3 for the sake of clarity, can be readily obtained by carrying out the analysis of sections 2.3 and 3 on the undeformed geometry and by setting σ 0 = 0 in Equation (15). However, when the cell is subjected to a compressive stress as shown in Figure 2B, the bandgap experiences a shift towards higher frequencies, being it located, in this case, in the [70.7 − 87.3] kHz frequency range. On the other hand, when the cell is subjected to traction, the lower bound of the bandgap is observed at 68.3 kHz while its upper bound at 84.3 kHz, which corresponds to a frequency downshift of approximately 3.5% of the bandgap with respect to the case of the PC plate under compression. These results suggest that, even in the elastic regime, a deformation due to prestress can lead to significant changes in the passband and bandgap behavior.
Furthermore, looking at the behavior of the energy velocity for the two deformed configurations, it is possible to infer that, similarly to the case of homogeneous plates (Pau and Lanza di Scalea, 2015), a generic state of compression leads to an increase of the energy velocity, whereas a state of traction leads, in general, to its decrease. This situation can be readily deduced from the inspection of the energy velocity band diagrams in Figure 3, from which it can be noted that the largest variations take place in proximity of the cutoff frequencies. Moreover, in the case of prestressed homogeneous plates (Dubuc et al., 2017), characteristic spikes in the energy velocity dispersion curves appear due to mode coupling. This behavior can also be observed for the third fundamental mode at around 36 kHz when the PC plate is subjected to a compression.
Finally, in the case of the PC plate under compression, it is possible to observe the presence of a cutoff wavenumber for the first fundamental mode at about 108 rad/m and null frequency. This value can be derived by solving the eigenvalue problem reported in Equation (10) in k(ϑ) by imposing ω = 0. In this case, the eigenvalue problem corresponds to that of the buckling load and, consequently, the lowest eigenvalue (the cutoff wavenumber) corresponds to a bifurcation point that indicates the onset of buckling for the infinite PC plate. It should be noted that the specific value of the cutoff wavenumber depends on the magnitude and distribution of the compressive stress in the PC.

CONCLUSIONS
An Updated Lagrangian computational scheme has been presented for the calculation of the band diagrams of phononic crystal plates subjected to a generic state of prestress. The scheme involves the solution of a geometrically nonlinear static problem for a representative unit cell and the application of the Floquet-Bloch theorem to the linearized equations of acousto-elasticity on the deformed configuration of the unit cell undergoing the static load. For the case of an epoxy plate with cylindrical steel   inclusions, it has been demonstrated that the existence of a prestress state of compression or tension can lead to significant changes in the passbands and bandgaps of a phononic crystal plate, even in the case of prestress inducing full reversible deformation state, which is of more practical interest with respect to the case of a permanently induced deformation. It has been observed that when the cell is subjected to a compressive stress, the bandgap experiences a shift toward higher frequencies, whereas when the cell is subjected to traction, we observe a frequency downshift of approximately the 3.5% of the bandgap with respect to the case of the phononic plate under compression.
Similarly to the case of homogeneous plates, a generic state of compression leads to an increase of the energy velocity, whereas a state of traction tends to lower its value. The largest variations have been observed in proximity of the cutoff frequencies.
Characteristic spikes in the energy velocity dispersion curves appear due to mode coupling for the third fundamental mode at around 36 kHz when the phononic plate is subjected to a compression.
The results presented herein provide insights in the behavior of phononic crystal plates with tunable dispersive properties, and suggest new leverages for wave manipulation valuable in many application fields such as wave filters, waveguiding and beam splitting, sensing devices, and vibration shielding.

AUTHOR CONTRIBUTIONS
MMa and MMi conceived the idea and wrote the first draft of the paper. MMa performed numerical calculations. All the authors discussed the manuscript.

APPENDIX A
The operators L i and L 0i are expressed as

APPENDIX B
The band diagram in terms of wavenumber k and energy velocity c e for the undeformed unit cell of Figure 2 at ϑ = 0 is reported in Figure A1.  Figure 2A.