Wave sensitivity analysis for periodic and arbitrarily complex composite structures

Purpose 
 
 
 
 
This paper aims to present the development of a numerical continuum-discrete approach for computing the sensitivity of the waves propagating in periodic composite structures. The work can be directly used for evaluating the sensitivity of the structural dynamic performance with respect to geometric and layering structural modifications. 
 
 
 
 
Design/methodology/approach 
 
 
 
 
A structure of arbitrary layering and geometric complexity is modelled using solid finite element (FE). A generic expression for computing the variation of the mass and the stiffness matrices of the structure with respect to the material and geometric characteristics is hereby given. The sensitivity of the structural wave properties can thus be numerically determined by computing the variability of the corresponding eigenvalues for the resulting eigenproblem. The exhibited approach is validated against the finite difference method as well as analytical results. 
 
 
 
 
Findings 
 
 
 
 
An intense wavenumber dependence is observed for the sensitivity results of a sandwich structure. This exhibits the importance and potential of the presented tool with regard to the optimization of layered structures for specific applications. The model can also be used for computing the effect of the inclusion of smart layers such as auxetics and piezoelectrics. 
 
 
 
 
Originality/value 
 
 
 
 
The paper presents the first continuum-discrete approach specifically developed for accurately and efficiently computing the sensitivity of the wave propagation data for periodic composite structures irrespective of their size. The considered structure can be of arbitrary layering and material characteristics as FE modelling is used.


Introduction
Layered and complex structures are nowadays widely used within the aerospace, automotive, construction and energy sectors with a general increase tendency. The wave propagation data for such structures are often employed for energy harvesting, vibration control, health monitoring and vibroacoustic transmission modelling purposes. Optimizing the layer characteristics of such structures for certain objectives is often a challenging task due to the large number of varying parameters to be considered as well as due to the lack of exact modelling approaches. The same it true about taking into account for the effect of these parametric uncertainties on the structural behaviour.
The numerical analysis of wave propagation within periodic structures was firstly considered in the pioneering work of the author of [1]. The work was extended to two dimensional media in [2]. In [3] the authors proposed a structured linearization method using state space eigenvalue problem for large matrices among other considerations for smoothening the solution of wave propagation in periodic structures. The WFE method was introduced in [4,5] in order to facilitate the post-processing of the eigenproblem solutions and further improve the computational efficiency of the method. The vibration of a uniform waveguide using the same technique was investigated by the authors of [6,7,8]. The WFE method for two dimensional structures was introduced in [9] . The same method was used in [10] in order to compute the dynamic response of two dimensional infinite structures. The WFE has recently found applications in predicting the vibroacoustic and dynamic performance of composite panels and shells [11,12,13,14,15,16,17] , with pressurized shells [18,19] and complex periodic structures [20,21,22] having been investigated. The variability of acoustic transmission through layered structures [23,24], as well as wave steering effects in anisotropic composites [25] have been modelled through the same methodology.
Structural sensitivity analysis is of great importance for understanding the overall impact of a design parameter variation to the structural performance which is to be optimised. Accurate sensitivity models are an important tool for design optimization, system identification as well as for statistical mechanics analysis. Many authors [26,27,28,29] have been focusing on the eigenvalue derivative analysis of a structural system. With regard to the variability analysis of the waves travelling within a structural medium the conducted work has been mainly focused on deriving expressions [30,31] of the stochastic wave parameters from analytical models. In [32] the authors conduct a design sensitivity analysis by a wave based approach. Considering numerical approaches the authors in [33] used Bloch's theorem in conjunction with the FE method in order to calculate the sensitivity of the acoustic waves within an auxetic honeycomb, while with regard to the computation of the variability of the propagating waves the authors in [34,35] have presented a stochastic WFE method approach for computing the stochastic wave propagation in one dimensional media.
In this paper a continuum-discrete approach for efficiently computing the sensitivity of the wave propagation data for periodic structures is presented.
The considered structure can be of arbitrary layering and material characteristics as FE modelling is employed. The effect of local parameter variation (e.g. varying the stiffener thickness or adding mass to a single location) can also be considered. The sensitivities of the mass and stiffness matrices for a solid FE are computed with respect to any structural parameter including the material characteristics and the thickness of the element. The sensitivity of the propagating structural waves can thus be numerically determined. The exhibited approach is validated against the FD method as well as analytical results.
The paper is organized as follows: In Sec.2 the formulation of the sensitivity of the waves propagating within the periodic structure is elaborated.
More precisely in Sec.2.1 a general approach for two dimensional waveguides is adopted while in Sec.2.2 a more efficient approach for one dimensional waveguides is considered. In Sec.3 the method is validated by comparison to analytical results for a metallic waveguide as well as by comparison to a FD approach for a layered sandwich waveguide. Conclusions on the presented work are given in Sec.4, while the formulation of the generic expressions for the stiffness and mass matrices of a solid FE are presented in the Appendix A.
2. Sensitivity of the wave propagation data 2.1. General formulation

Formulation of the PST for an arbitrary structural segment
A periodic segment of a panel having arbitrary layering and complexity is hereby considered (see Fig.1) with l x , l y its dimensions in the x and y directions respectively. The segment is modelled using a conventional FE software. The mass and stiffness matrices of the segment M and K are extracted and the DoF set q is reordered according to a predefined sequence such as: corresponding to the internal, the interface edge and the interface corner DoF (see Fig.1). The free harmonic vibration equation of motion for the modelled segment is written as: The analysis then follows as in [36] with the following relations being assumed for the displacement DoF under the passage of a time-harmonic wave: with ε x and ε y the propagation constants in the x and y directions related to the phase difference between the sets of DoF. The wavenumbers k x , k y are directly related to the propagation constants through the relation: with x the reduced set of DoF: The equation of free harmonic vibration of the modelled segment can now be written as: with * denoting the Hermitian transpose. The most practical procedure for extracting the wave propagation characteristics of the segment from Eq. (6) is injecting a set of assumed propagation constants ε x , ε y . The set of these constants can be chosen in relation to the direction of propagation towards which the wavenumbers are to be sought and according to the desired resolution of the wavenumber curves. Eq.(6) is then transformed into a standard eigenvalue problem and can be solved for the eigenvector x which describe the deformation of the segment under the passage of each wave type at an angular frequency equal to the square root of the corresponding eigenvalue A complete description of each passing wave including its x and y directional wavenumbers and its wave shape for a certain frequency is therefore acquired. It is noted that the periodicity condition is defined modulo 2π, therefore solving Eq.(6) with a set of ε x , ε y varying from 0 to 2π will suffice for capturing the entirety of the structural waves. Further considerations on reducing the computational expense of the problem are discussed in [36].
It should be noted that only propagating waves will be considered in the subsequent analysis. Evanescent waves can also be captured by introducing imaginary values for ε x , ε y however the precise computation of these waves would require a very fine resolution of the propagation constants, drastically increasing the computational effort.

Parametric sensitivity
It is noted that matrices K = R * KR and M = R * MR in Eq.(6) are Hermitian therefore their resulting eigenvalues are real and the eigenvectors will be orthogonal. Assuming the known eigenvalue λ 0i and the corresponding eigenvector x 0i for the problem described in Eq.(6), then the following is true Now if matrices K 0 ,M 0 are changed by a small amount, say δK,δM then the eigenvalue λ 0i and the corresponding eigenvector x 0i will also be perturbed so that A direct consequence of the orthogonality of the eigenvectors is that with δ k i the Kronecker delta function. Substituting Eq.(8) into Eq.(7), we get then expanding and using Eq.(7) we can write and removing the higher-order terms simplifies to The orthogonality properties of the unperturbed eigenvectors in Eq.(6) allow for using them as a basis for expressing the perturbed eigenvectors. That is the perturbation of the eigenvector x 0i can be expressed as with ǫ ik small unknown constants. Substituting Eq.(13) into Eq. (12) gives which can also be written as Making use of Eq.(7) it is true that Again due to the orthogonality of the eigenvectors in Eq.(9), the summations can be removed by left multiplying by x ⊤ 0i , therefore and eliminating the equal terms gives Rearranging the expression with regard to δλ i we get However, because of Eq.(9) the expression of the eigenvalue perturbation is given as When the partial derivatives of K,M with regard to a design parameter β are known then the sensitivity of an eigenvalue λ i to this design parameter will be equal to The global mass and stiffness matrices M, K of the structural segment are formed by adding the local mass and stiffness matrices of individual It is known however that Eq. (22) can therefore be written as: with ω i the angular frequency at which the set of ε x , ε y is true for the i wave type described by the x 0i deformation. The generic symbolic expressions of the m, k matrices for an orthotropic structural segment modelled with a linear solid FE is given in Appendix A. The wavenumber sensitivity ∂k i ∂β can be deduced as with c g,i = ∂ω i ∂k i the group velocity associated with the wave type i at frequency ω i which can be evaluated by avoiding any numerical differentiation as exhibited in [37].

Condensation process
When wave propagation is considered only in the x direction, the problem can be condensed using a transfer matrix approach as in [4]. The DSM of the structural segment is again partitioned, this time with regard to its left/right sides and internal DoF with q the displacement and f the force vectors. Using a Guyan type condensation for the internal DoF the problem can be expressed as Assuming that no external forces are applied on the segment the displacement continuity and force equilibrium equations at the interface of two consecutive periodic segments s and s + 1 give Using Eqs. (27),(28) the relation of the displacements and forces of the left and right sides of the segment can be written as and the expression of the symplectic transfer matrix T can be written as Free wave propagation is described by the eigenproblem whose solution γ i is related to the structural wavenumber k xi by

Parametric sensitivity
Assuming the known eigenvalue γ 0i and the corresponding eigenvector for the problem described in Eq.(32), then the following is true Now if the matrix T 0 is changed by a small amount, say δT then the eigenvalue γ 0i and the corresponding eigenvector z 0i will also be perturbed so Expanding we have However using Eq.(34) and removing the higher-order terms eventually simplifies the relation to The symplectic orthogonality properties [5] of the unperturbed eigenvectors of T 0 suggest that if γ i is an eigenvalue of T 0 with the corresponding right eigenvector z i , then 1/γ i will also be an eigenvalue of T 0 with a corresponding left eigenvector y i = z i J n , with J n a 2n × 2n matrix given as The following will also apply if the eigenvectors are appropriately scaled and with Γ 0 the diagonal matrix containing the eigenvalues of T 0 , while is the matrix containing the left eigenvectors of T 0 and the matrix containing the right eigenvectors of T 0 . The orthogonality property allows for expressing the perturbation of the eigenvector z 0i as with ǫ ik small unknown constants. Substituting Eq.(42) into Eq.(37) and making use of Eq.(34) once again gives and due to Eq. (38), the summations can be removed by left multiplying by Eliminating the equal terms gives Due to the orthogonality of the eigenvectors it is true that therefore the expression of the eigenvalue perturbation δγ i as a function of the transfer matrix perturbation δT can be written as When the partial derivative of T with regard to a design parameter β is known then the sensitivity of an eigenvalue γ i to β is equal to It is known however that Eq.(48) can therefore be written as:  Fig.2 refers to the sensitivity of ω at which a specific ε x , ε y set occurs for the first bending wave propagating in the structure. Real values are considered for ε x , ε y throughout this paper. The result is compared to the CPT theory by differentiating the expression relating the flexural wavenumber k f to the angular frequency ω f the flexural stiffness of a monolayer panel. Very good agreement is observed between the two curves with the maximum divergence being equal to 0.18%. The perturbation of ω for the first membrane wave propagating in the structure with relation to E is presented in Fig.3. The result is compared to the CPT theory by differentiating the expression relating the membrane wavenumber k m to the angular frequency ω m Again excellent agreement is observed between the analytical and numerical results. The perturbation of ω for the first shear wave propagating in the structure with relation to E is presented in Fig.4. The result is com- Good agreement is generally observed between the results, however there is a maximum discrepancy of 4.3% between the two curves. This is probably due to an insufficiency of the analytical approach for modelling the dynamic behaviour at higher frequencies. Again an increase in thickness will stiffen the structure, shifting the same wavenumber towards higher frequencies. There is no dependence of the shear and membrane waves on the thickness of the results are compared to a FD approach throughout this section. In order to implement the FD approach a perturbation of 0.001% was considered for each structural parameter. The resulting FD sensitivity can be computed by with ω p the perturbated angular frequency at which a certain wavenumber value occurs for β p and ω 0 the corresponding angular frequency for β 0 . The perturbation of ω with respect to E f l , E f u of the structure is presented in  Fig.13. In both cases it can be observed that increasing v for the facesheets will shift the flexural wavenumber curve towards higher frequencies. This phenomenon however will be much more intense for δv f u especially at higher wavenumber values. Again excellent agreement is observed between the FD results and the presented approach.
The perturbation of ω with respect to v c for the sandwich core is presented in Fig.14. The effect of δv c is softening up to a certain wavenumber value, beyond which a sudden exponential increase is observed which stiffens the flexural structural behaviour. The perturbation of ω with respect to ρ f l and ρ f u is presented in Fig.15. As expected both δρ f l and δρ f u will shift the wavenumber curve to lower frequencies. This effect will be greater for the at which the effect of δρ f l and δρ f u will be the same. Beyond this critical wavenumber the softening effect will be more intense for δρ f l than δρ f u . The perturbation of ω with respect to ρ c for the sandwich core is presented in

Conclusions
A numerical continuum-discrete approach for computing the sensitivity of the waves propagating in periodic structures was presented in this paper.
The main conclusions of the work are as follows: (i) A formulation of the two dimensional wave propagation problem was exhibited; this was followed by the derivation of the sensitivity of the waves propagating in the periodic structure. The structure can be of arbitrary layering, material characteristics and geometric complexity. The effect of local parameter variation can also be considered.
(ii) A reduced transfer matrix formulation of the wave propagation analysis in one dimensional periodic media was presented, followed by the derivation of the sensitivity of the wavenumber values. As for the first formulation, the structure can be of arbitrary layering, material characteristics and geometric complexity. This approach however tends to be less computationally efficient due to the involved inversions and multiplications of DSM partitions.
(iii) A monolayer panel was modelled by the presented approach and the resulting sensitivity values were validated through the CPT.
(iv) Moreover, an asymmetric sandwich structure was also modelled. The sensitivity of the propagating flexural wave with respect to the structural parameters of the facesheets and the core were computed and the results were successfully compared to the ones provided by the FD method.
y 1 y 2 y 3 y 4 y 5 y 6 y 7 y 8 The displacement interpolations are expressed as Linear shape functions are assumed for the element The element stiffness matrix k is formally given by the volume integral while the element mass matrix m can be determined as while ρ m is the mass density of the material. It is also noted that The Jacobian matrix of the element is while the the flexibility matrix of the element for an orthotropic material D −1 can generally be written as The assumption of the undeformed FE being a rectangular parallelepiped is hereby adopted. The coordinates x 1 , x 2 , x 3 , x 4 , x 5 , x 6 , x 7 , x 8 , y 1 , y 2 , y 3 , y 4 , y 5 , y 6 , y 7 , y 8 , and z 1 , z 2 , z 3 , z 4 , z 5 , z 6 , z 7 , z 8 , can then be replaced by L x , L y , L z in the expression of B. The generic expression for m is thus given as while the symbolic generic expression of k can be derived exactly in the same way but is hereby intentionally omitted for the sake of brevity.
The generic sensitivity expressions ∂k ∂β i , ∂m ∂β i as well as ∂ 2 k ∂β j ∂β i , ∂ 2 m ∂β j ∂β i with β i , β j being design parameters can therefore be calculated as a function of E x , E y , E z , v xy , v xz , v yz , G xy , G xz , G yz , L x , L y , L z by differentiating over the generic expressions for k, m.