Free Vibrations of Sandwich Plates with Damaged Soft-Core and Non-Uniform Mechanical Properties: Modeling and Finite Element Analysis

The paper aims to investigate the natural frequencies of sandwich plates by means of a Finite Element (FE) formulation based on the Reissner-Mindlin Zig-zag (RMZ) theory. The structures are made of a damaged isotropic soft-core and two external stiffer orthotropic face-sheets. These skins are strengthened at the nanoscale level by randomly oriented Carbon nanotubes (CNTs) and are reinforced at the microscale stage by oriented straight fibers. These reinforcing phases are included in a polymer matrix and a three-phase approach based on the Eshelby-Mori-Tanaka scheme and on the Halpin-Tsai approach, which is developed to compute the overall mechanical properties of the composite material. A non-uniform distribution of the reinforcing fibers is assumed along the thickness of the skin and is modeled analytically by means of peculiar expressions given as a function of the thickness coordinate. Several parametric analyses are carried out to investigate the mechanical behavior of these multi-layered structures depending on the damage features, through-the-thickness distribution of the straight fibers, stacking sequence, and mass fraction of the constituents. Some final remarks are presented to provide useful observations and design criteria.


Introduction
Since its early development and the publication of the first research papers [1][2][3][4][5][6], the Finite Element (FE) method has shown its potentiality in solving easily and accurately many structural problems which could not be solved analytically. Nowadays, this feature is even more emphasized by the great and continuous technological advancements reached in computer sciences in terms of available computational resources. As highlighted in the books which can be certainly considered as milestones in the development of the FE method [7][8][9][10][11][12][13], the easy implementation and the possibility to reduce complex continuous problems into simpler discrete ones has encouraged its rapid spread among many researchers and engineers [14][15][16][17][18][19][20]. along each layer. Nevertheless, this approach is quite onerous in terms of computational resources. Effective and accurate solutions can be obtained more easily by introducing the Murakami's function in the displacement fields of Equivalent Single Layers (ESL) models, in which all the degrees of freedom are defined within the reference surface of the structure independently from the number of layers. The Murakami's function, in fact, is able to capture the zig-zag effect in an accurate manner without increasing excessively the computations [100]. In general, classical models for laminated structures, such as the well-known Reissner-Mindlin (RM) theory, are not able to take into account this effect. Nevertheless, these theories could be also enhanced by adding the Murakami's function in their kinematic models. The results that can be obtained are accurate and the computational cost is reduced, if compared to the one that characterizes higher-order ESL or LW models [92]. Therefore, in the present paper the in-plane displacement field of the RM model is enriched by the Murakami's function in order to capture the effective behavior of sandwich plates with damaged soft-core and non-uniform mechanical properties. It should be mentioned that further details concerning Zig-zag theories can be found in the review paper by Carrera [94].
Finally, a brief description of the research outline is presented. The geometric and mechanical characterization of the plates are illustrated after this introduction in Section 2. In particular, the multi-phase approach including the Eshelby-Mori-Tanaka scheme and the Halpin-Tsai homogenization procedure is described to provide the mechanical properties of the skin. In the same Section, the damage model is introduced for the isotropic soft-core. On the other hand, Section 3 is focused on the theoretical aspects of the Reissner-Mindlin Zig-zag (RMZ) theory. The corresponding FE model is developed, and the fundamental system of equations is deducted by means of the Hamilton's principle [36,101]. The results of the numerical applications are presented in Section 4. Here, the effect of the Murakami's function is discussed by means of the comparison between the RM and RMZ approaches, and the proposed model is validated as well. Then, several applications are illustrated to investigate the effects of the progressive damage, the non-uniform distribution of the fiber volume fraction, the in-plane fiber orientation, and the material properties on the natural frequencies of the structures and on their corresponding mode shapes. Finally, the matrix form of the fundamental operators required in the governing equations are defined in Appendix A.

Geometric and Mechanical Characterization
The paper is focused on the vibrational behavior of laminated sandwich plates with an inner damaged soft-core. The plates under consideration are characterized by a planar size given by the lengths of their sides L x , L y , in which x, y denote the principal directions of the local reference system. The extension along the coordinate z is specified by the overall thickness h of the composite structures, which is given by h = 2h s + h c , where h s stands for the thickness of the external face-sheets and h c is the thickness of the core. Thus, the plates consist in three layers, in which the external ones are made of orthotropic materials and have the same thickness whereas the central one is isotropic. The upper and lower thickness coordinates of the generic k-th ply are denoted by z k+1 , z k , respectively. The geometric features of a plate element are shown in Figure 1.
The mechanical characterization of the k-th layer is carried out in terms of the engineering constants of the orthotropic material, which are the Young's moduli E (k) 11 , E (k) 22 , the shear moduli G (k) 12 , G (k) 13 , G (k) 23 , and the Poisson's ratio ν (k) 12 . As far as the isotropic core is concerned (k = 2), only two independents parameters are needed, which are the Young's modulus E (k) = E .
(1) In the following sections, the evaluation of these engineering constants is discussed for both the external layers and the soft-core. It should be specified that a perfect bonding is assumed between two adjacent layers in the proposed model. In the following sections, the evaluation of these engineering constants is discussed for both the external layers and the soft-core. It should be specified that a perfect bonding is assumed between two adjacent layers in the proposed model.

Mechanical Properties of the Face-Sheets
A multiscale approach is employed to evaluate the overall mechanical properties of the threephase composite face-sheets [38,39]. These layers are made of fiber-reinforced materials, which means that a polymer matrix (epoxy resin) is strengthened by oriented straight fibers (Carbon fibers). The third constituent is given by randomly oriented Carbon nanotubes (CNTs), which are inserted in the matrix and are used to further increase its mechanical features. The overall properties are computed by means of a two-step approach. Firstly, the Eshelby-Mori-Tanaka scheme is applied to obtain the properties of the epoxy resin including CNTs [52,53]. At this stage, the composite turns out to be isotropic due to the fact that the nanoparticles are randomly oriented, as illustrated in the paper by Shi et al. [51]. Then, the Halpin-Tsai approach is applied to combine the features of the enriched matrix with the properties of the reinforcing straight fibers [55][56][57]. At the nano-scale level, the single fiber of CNT is modeled as a linear-elastic, transversely isotropic and homogeneous cylindrical solid, as proposed in the paper by Odegard et al. [50]. Its mechanical properties are described by five parameters, which are the Hill's elastic moduli [58,59], defined as , , , , k l m n p . Its density C ρ is required to compute the volume fraction of CNTs C V as follows: where M ρ represents the density of the polymer matrix, whereas C w stands for the mass fraction of CNTs. A uniform distribution of CNTs is assumed in each layer along the thickness direction. It should be recalled that the matrix volume fraction is given by 1 where ν F 21 = E F 22 ν F 12 /E F 11 . As shown in the previous step, the density of the fibers ρ F is required to compute the reference value of the corresponding volume fraction V F , once their mass fraction w F is defined Such constant quantity can be multiplied by a peculiar function f (k) (z) which depends on the thickness coordinate z. Consequently, the volume fraction distribution of the fibers V F is given by and it is clearly non-uniform along the thickness of the face-sheets.
In the current paper, two different functions f (k) (z) are introduced to specify the volume fraction distribution V F . It should be noted that they can be chosen arbitrarily in the two face-sheets of the sandwich structures under consideration. Power-law functions are used to this aim and the definitions of f (k) (z) are specified below: where α represents the arbitrary exponent of the distributions. As shown in Figure 2, several configurations can be obtained according to the value given to α ∈ [0, ∞]. It should be noted that the extreme values of α, which correspond to the constant values f (k) = 0, f (k) = 1, are able to characterize a uniform distributions of the fiber along the thickness or their absence (as a consequence, the polymer matrix is the only constituent of the layer).  At this point, the Halpin-Tsai approach can be applied to obtain the Hill's elastic moduli of the composite layers, which are denoted by , , , , k l m n p : The introduction of the function f (k) (z) in the definition of V F causes the dependency on the thickness coordinate z of each engineering constant specified in Equation (13). Therefore, one gets Finally, the density of the composite face-sheets is given by: In the following, the same constituents are used in the external layers assuming also the same values of the mass fractions of both CNTs and fibers. The relation ν 11 is also required to compute the Poisson's ratio ν (k) 21 . As emphasized in the introduction and illustrated in the paper [38], different approaches and homogenization techniques could be used to the same aim. For instance, the rule of the mixture represents the most exploited methodology.

Mechanical Properties of the Damaged Matrix
The core of the sandwich structures considered in the paper is made of the same polymer matrix used in the face-sheets. Nevertheless, a damage model is introduced to provide an analytical description of an irreversible rheological process that causes the decay of the mechanical properties, in terms of engineering constants. An isotropic damage is considered in the following, which is fully characterized by a scalar D as illustrated in the book by Lemaitre and Chaboche [21]. The elastic modulus of the damaged material is given by: for k = 2, in which E M is the original value of matrix Young's modulus, for 0 ≤ D < 1 It is clear that D = 0 identifies a virgin material, whereas a fully damaged material is characterized by D = 1. Having in mind relation (1), the shear modulus is subjected to the same damage. On the contrary, the Poisson's ratio ν (k) = ν M and the density ρ (k) = ρ M of the core (k = 2) are kept constant. Finally, it should be specified that the damage does not depend on the spatial coordinates and affects uniformly the core. For conciseness purposes, the mechanical properties of the undamaged epoxy resin and the Carbon fibers are summarized in Table 2.

Finite Element Model Based on A First-Order Zig-Zag Plate Theory
A linear theory is used to model the mechanical behavior of sandwich plates with an inner soft-core. With respect to the well-known Reissner-Mindlin (RM) theory, the in-plane expansion is enriched by two more degrees of freedom, which are able to capture the zig-zag effect [92][93][94][95][96]. Here, the corresponding (FE) formulation is presented. The displacement field for the generic e-th element is given by: z (x, y, t) (16) in which the three-dimensional displacement components are denoted by U z . The spatial coordinates of the plate are given by x, y, z, as shown in Figure 1, whereas t is the time variable. The zig-zag effect is modeled by means of the Murakami's function F z (z) defined below for a multilayered structure: where k identifies the generic layer. This function allows to introduce a discontinuity in the slope of the three-dimensional displacements U x (e) , U y (e) along the thickness direction at each layer interface. Its meaning is well-described in the papers by Carrera [94]. It should be noted that the current model is characterized by seven degrees of freedom per node, two more than the classical RM approach.
In particular, u x , ψ (e) y denote the magnitude of the zig-zag effect. A nine-node quadratic rectangular element is used to develop the FE formulation and each degrees of freedom is approximated by means of Lagrange interpolating functions N i , for i = 1, . . . , 9. The node numbering is illustrated in Figure 3. The node numbering is illustrated in Figure 3.  Due to this approximation, the degrees of freedom can be written as a function of the corresponding where N = N 1 · · · N 9 is the vector of the Lagrange interpolating functions. Analogously, the nodal displacements are defined in vector form as follows: This notation is useful to define the vector u (e) which includes all the nodal degrees of freedom: The interpolating functions N i assume the well-known definitions presented in the book by Reddy [11]. These polynomials are conveniently expressed as functions of the natural coordinates ξ, η, with ξ, η ∈ [−1, 1], introduced in the so-called master element (or parent element) depicted in Figure 3. Thus, the interpolating functions become N i = N i (ξ, η).
The Lagrange functions N i are also employed to define the geometric shape of each element. According to the principles of an isoparametric formulation, the coordinates x (e) , y (e) within the generic e-th element can be defined as follows: in which the i-th node of the element under consideration is identified by the couple of nodal coordinates i . Such coordinates are included in the corresponding vectors x (e) , y (e) , which assume the following aspects: The isoparametric formulation allows to move easily all the computations in the parent space. To this aim, the Jacobian matrix J is required to perform the coordinate change. In order to define this matrix, the derivates of the interpolating functions with respect to the natural coordinates ξ, η are needed and are collected in the corresponding vectors B ξ , B η defined below: At this point, the Jacobian matrix J can be introduced as specified in [11]: Assuming that the determinant of the Jacobian matrix is positive, the matrix J can be inverted. Since in the following only regular rectangular elements are considered, this assumption is always satisfied and the matrix J −1 is admissible and can be used to compute the derivatives of the interpolating functions in the physical domain defined by the coordinates x, y. The following relation is needed for this purpose: where B x and B x collect the derivatives of the shape functions with respect to x and y, respectively. These operators are required to define the compatibility equations of the RMZ model. In particular, the three-dimensional strain components for a rectangular plate can be obtained by means of the elasticity equations applied to the displacement fields (16) and assume the following definitions for the e-th element: where ε xy are the membrane strains, and γ yz are the transverse shear strains. It should be observed that the normal strain along the thickness direction ε (e) zz is omitted since the plane-strain assumption entails that ε (e) zz = 0. The generalized strains related to the plate middle surface can be easily defined from relations (26). The superscript "0" denotes those quantities that are included also in the well-known RM theory, whereas the terms related to the zig-zag effect are identified by the superscript "1". In particular, it should be noted that ε 0(e) yz the shear strains, which are also defined in the RM theory [36]. The constitutive equations are now used to characterize the stress components in the k-th layer of the laminate. The following definitions imply that the plane-stress assumption σ k(e) zz = 0 is assumed by hypothesis, whereas the other stress components are given by: ij denotes the stiffnesses of the k-th orthotropic layer evaluated in the geometric reference system. These parameters have the same meaning in each finite element, since the mechanical properties of the structure do not vary in the plate middle surface. Their well-known definition, which depends on the parameters Q (k) ij expressed as a function of the engineering constants of the k-layer defined below, can be found in the book by Reddy [36]: Quantities in (28) for k = 1, 3 depend on the thickness coordinate z, since the reinforcing fibers of the face-sheets are characterized by a non-uniform distribution in this direction. Each orthotropic layer can be also characterized by an arbitrary orientation θ (k) . The notation θ (1) /core/θ (3) is used in the next Sections to specify the orientations of the reinforcing fibers in the face-sheets and the consequent lamination scheme. Conventionally, in the above notation the layer numbering always starts from the bottom surface of the plate.
At this point, the Hamilton's variational principle should be applied to obtain the governing equations for the dynamic problem under consideration [36,101]. If t 1 , t 2 specify the boundary values of the considered time interval, the variational principle assumes the following aspect within the e-th discrete element: where δK (e) is the variation of the kinetic energy, whereas δΦ (e) represents the variation of the elastic strain energy. The kinetic energy of the sandwich structure δK (e) is given by: where the double-dot notation specifies the second-order derivatives with respect to the time variable.
On the other hand, the elastic strain energy δΦ (e) is defined as follows: The proper mathematical manipulations of the elastic strain energy provide the definitions of the stress resultants as the through-the-thickness integrals of the stress components. The stress resultants can be written in matrix form by means of the following relations, which provide their definitions as a function of the generalized strain components introduced before: The superscripts "0" and "1" have the same meaning discussed previously. In particular, xy the bending and twisting moments, and T 0(e) x , T 0(e) y the shear forces, which are included also in the RM theory [11]. The other terms are related to the zig-zag effect. The elements of the constitutive operator in (32) are now discussed. Firstly, the following terms appear also in the RM model and have the same meaning [11]: On the other hand, the stress resultants related to the zig-zag effect require the following definitions, which include the Murakami's function and its derivative with respect to the thickness coordinate: It should be specified that the integrals in (33)- (34) are computed numerically since the stiffnesses ij can be arbitrary functions of z, due to the dependency on the thickness coordinate introduced by relation (11). Finally, it is important to specify that the shear forces need the shear correction factor κ. The value of 5/6 is used to this aim.
The system of dynamic equations for the problem under consideration represents the main result of the application of the Hamilton's principle. The fundamental equation for the generic element e is given by: the vector of the second-order time derivatives of the element degrees of freedom included in u (e) . The fundamental operators K (e) assumes the following definitions: On the other hand, the mass matrix M (e) is given by: The matrices K ij and M ij , for i, j = 1, . . . , 7 are defined in the Appendix A. In particular, the following inertia terms are required to compute the mass matrix: It can be easily observed that the terms I 0 , I 1 , I 2 are included also in the RM theory. On the other hand, the parameters I 3 , I 4 , I 5 are linked to the zig-zag effect and include the Murakami's function. It should be specified that these integrals are computed numerically, since the density ρ (k) , for k = 1, 3, is an arbitrary function of the thickness coordinate.
The fundamental Equation (35) is valid for each discrete subdomain at the element level. The well-known assembly procedure is required to obtain the corresponding global system of equations. This approach allows to enforce automatically the displacement continuity at the element interfaces. Therefore, a C 0 compatibility requirement is satisfied in the current approach. Once the global matrices are obtained, the fundamental system of equations assumes the following aspect: in which K, M are the global stiffness and mass matrices, respectively. The degrees of freedom of the whole structure are collected in the vector u. The nodal displacements are listed following the scheme specified by the dashed line in Figure 4, where an example of a discrete plate domain is also depicted.
in which , K M are the global stiffness and mass matrices, respectively. The degrees of freedom of the whole structure are collected in the vector u . The nodal displacements are listed following the scheme specified by the dashed line in Figure 4, where an example of a discrete plate domain is also depicted.
The second-order time derivatives of these quantities are collected in the vector u  following the same scheme. Finally, it should be recalled that the size of the fundamental operators ,

Numerical Computation of the Fundamental Matrices
The Gauss-Legendre quadrature rule is employed to compute the fundamental matrices , K M .
By definition, the integral of a generic function ( ) , G x y defined in a two-dimensional domain can be evaluated as follows: in which the determinant of the Jacobian matrix det J is introduced. Therefore, the integral is computed in the parent space in which the reference system is given by natural coordinates ,  If N P denotes the number of nodes, the model is characterized by N do f s = 7N P as number of degrees of freedom. Consequently, the vector u can be written as follows: u = u x,1 · · · u x,N P 1→N P u y,1 · · · u y,N P N P +1→2N P u z,1 · · · u z,N P The second-order time derivatives of these quantities are collected in the vectorü following the same scheme. Finally, it should be recalled that the size of the fundamental operators K, M is given by N do f s × N do f s . The discrete system can be solved once the proper boundary conditions along the edges of the domain are enforced. In the present paper, since only fully clamped plates are considered, all the nodal displacements related to the boundary edges are all equal to zero.

Numerical Computation of the Fundamental Matrices
The Gauss-Legendre quadrature rule is employed to compute the fundamental matrices K, M. By definition, the integral of a generic function G(x, y) defined in a two-dimensional domain can be evaluated as follows: in which the determinant of the Jacobian matrix detJ is introduced. Therefore, the integral is computed in the parent space in which the reference system is given by natural coordinates ξ, η. From the numerical point of view, this integral can be converted into the following weighted linear sum: where W I , W J are the weighting coefficients, whereas ξ I , η J are the points in which the integral is computed. These nodes are the roots of Legendre polynomials [11]. The full integration is performed considering 9 evaluation points, whereas the reduced one is carried out in four points only. The position of these nodes is shown in Figure 3. It should be specified that the reduced integration is employed only to compute the elements of the stiffness matrix related to the shear forces in order to avoid the shear locking issue [11]. The analytical values of the weighting coefficients for each root of the Legendre polynomials can be found in the book by Reddy [11].

Evaluation of the Natural Frequencies
The free vibration analysis is based on a generalized eigenvalue problem from the analytical point of view. In particular, the following relation allows to compute the circular frequencies ω of the structures under consideration: where the modal amplitudes are collected in the vector d. Once relation (43) is solved, the natural frequencies f measured in Hz can be easily computed as f = ω/2π.

Numerical Applications
The numerical applications presented in this Section aim to evaluate the natural frequencies of several fully clamped sandwich plates. The geometric features are the same in each computation. In particular, a square domain defined by L x = L y = 2.5 m is considered. The thickness of the external layers is h s = 0.02 m, whereas the soft-core is defined by h c = 0.06 m. Each structure is subdivided into 100 finite elements as far as the discrete domain is concerned.
Firstly, the validity of the current approach based on the RMZ model is proved and compared with the results that could be obtained by means of the well-known RM theory. To this aim, a three-dimensional FE model is built through a commercial code. Secondly, the model is also validated with respect to the application of non-uniform distributions of the reinforcing fibers along the plate thickness. For this purpose, the results are compared with the ones available in the literature. Then, several parametric investigations are presented to discuss the effects of the damage, the through-the-thickness distributions of the reinforcing fibers, the lamination scheme and the in-plane orientation of the fibers, the mechanical properties of fibers and CNTs on the vibrational response.

Influence of the Murakami's Function and Validation of the RMZ
The first test aims to prove the need of the Murakami's function when the mechanical behavior of sandwich structures with an inner soft-core must be analyzed. In this application, the core is made by a virgin material (D = 0) and the fibers are uniformly placed along the thickness of the external face-sheets. The mechanical characterization is fully accomplished by setting w C = 0.05 and w F = 0.80. The natural frequencies are obtained by using the RM and the RMZ models, for three different lamination schemes. The same structures are investigated by means of a three-dimensional FE commercial code (twenty-node brick elements), denoted by 3D-FE in the following. The software Strand7 is employed for this purpose. The first ten natural frequencies for the sandwich plate with an inner soft-core under consideration are shown in Table 3, where the percentage differences (%diff) of the RMZ and RM solutions with respect to the 3D-FE results are also highlighted. The following aspects can be observed:

•
The RMZ theory provides natural frequencies that are close to the results given by the reference solution (3D-FE). In fact, the maximum percentage difference is about 5% for higher modes. This difference is satisfactory having in mind the approximation introduced by a two-dimensional ESL theory; • The computational cost is very different. In particular, the number of degrees of freedom in the 3D-FE model is ten times the one needed by the RMZ theory to obtain similar values; • The RM model is not adequate to evaluate the natural frequencies of a sandwich soft-core structure, as it can be observed by the percentage differences with respect to the reference solution.
The number of degrees of freedom in this circumstance is even lower if compared to the other models, but the computational saving cannot justify the poor approximation of the solution.

Validation of the Model with Respect to Non-Uniform Distributions of the Reinforcing Fibers
The current approach is validated also with respect to the application of non-uniform distribution of the reinforcing fibers in the thickness direction. In the paper by Lei et al. [71], a fully clamped square plate, characterized by the aspect ratio L x /h = 10, is analyzed by means of the kp-Ritz method. The structure is made of an epoxy resin (E M = 2.1 GPa, ν M = 0.34, ρ M = 1150 kg/m 3 ) reinforced by aligned fibers of CNTs. This configuration can be modeled as a two-phase composite, in which the Carbon fibers are characterized by the following mechanical properties E F 11 = 5.6466 TPa, E F 22 = 7.0800 TPa, G F 12 = 1.9445 TPa, ν F 12 = 0.175, ρ F = 1400 kg/m 3 . It should be noted that the approach presented in this paper can be used to evaluate the overall mechanical properties of this structure by neglecting the effect of the randomly oriented CNT particles scattered in the matrix. In other words, one gets E * M = E M , ν * M = ν M and ρ * M = ρ M . On the other hand, the aligned CNTs assume the same role of the straight fibers. Equation (10) is employed to obtain the value of w F which provides the constant V F = 0.11 specified in the reference paper.
In order to validate the current methodology, the plate is made of two orthotropic layers of equal thickness (no soft-core is included), characterized by non-uniform distributions of the CNT fibers. "Case 1" is obtained by f 2 with α = 1, whereas "Case 2" is given by f The frequencies are presented in Table 4 in dimensionless form as (44) in which ω denotes the circular frequencies. As specified in the previous sections, the Halpin-Tsai (HT) model is applied for the mechanical characterization of the structure. Nevertheless, the reference solutions are obtained by means of the rule of the mixture (MIX). The same approach has been used in the paper by Bacciocchi and Tarantino [39] for similar purposes. Therefore, only in the next application these two models (MIX and HT) are considered for the sake of comparison. Small differences are obtained depending on the homogenization procedure used in the computation. It should be emphasized that the current application does not aim to investigate the effect of the homogenization method but only the comparison with the reference solution of the frequency parameter. As it can be noted from the results shown in Table 4, a good agreement is observed with the reference solution, especially if the rule of the mixture is employed as it could be expected. In the same table, the FE results provided by a commercial code are also presented. Finally, it should be specified that the RMZ theory is considered as far as the theoretical model of the present solutions is concerned.

Effect of Damage
The same geometric and mechanical configurations are analyzed also in this paragraph. Nevertheless, the soft-core of the structures is affected by a decay of the mechanical properties and the effect of an increasing damage parameter D ∈ [0, 1] on the natural frequencies is discussed. For each lamination scheme, four different through-the-thickness distributions of the reinforcing fibers in the face-sheets are investigated, including the uniform one. "Scheme 1" denotes the uniform distribution of the fibers; "Scheme 2" is obtained by setting f 1 with α = 1; finally, "Scheme 4" has the same functions of the previous one, assuming α = 2. These configurations are graphically depicted in Figure 5. The volume fraction distribution of the fibers in the core is clearly equal to zero. The results are presented in Figure 6, where it can be observed that the first natural frequency for each configuration depends noticeably on the parameter D . The following aspects should be noted, as well: • The decrease of the frequency is clearly caused by the corresponding stiffness reduction of the structures. This expected tendency models accurately the physical behavior of structures with a lower value of stiffness. In fact, by increasing the value of D up to the unity (fully damaged core), the frequency would tend to zero; • The same behavior is obtained for each volume fraction distribution, but the maximum value of the first frequency that can be reached depends on the through-the-thickness distributions of F V ( Figure 5); • These aspects can be noted for each lamination scheme. Nevertheless, depending on the in-plane fiber orientation, the value of the first frequency could change. In addition, a peculiar choice of lamination scheme could reduce the influence of the through-the-thickness distributions of the volume fraction F V , since the curves related to the various schemes are less detached; • Finally, similar graphs could be obtained also for higher frequencies. The results are presented in Figure 6, where it can be observed that the first natural frequency for each configuration depends noticeably on the parameter D. The following aspects should be noted, as well: • The decrease of the frequency is clearly caused by the corresponding stiffness reduction of the structures. This expected tendency models accurately the physical behavior of structures with a lower value of stiffness. In fact, by increasing the value of D up to the unity (fully damaged core), the frequency would tend to zero; • The same behavior is obtained for each volume fraction distribution, but the maximum value of the first frequency that can be reached depends on the through-the-thickness distributions of V F ( Figure 5); • These aspects can be noted for each lamination scheme. Nevertheless, depending on the in-plane fiber orientation, the value of the first frequency could change. In addition, a peculiar choice of lamination scheme could reduce the influence of the through-the-thickness distributions of the volume fraction V F , since the curves related to the various schemes are less detached; • Finally, similar graphs could be obtained also for higher frequencies.

Influence of the Exponent of the Through-the-Thickness Distribution of the Fiber Volume Fraction
The current application deals with the effect of the exponent α that characterizes the through-the-thickness distribution of V F . For this purpose, the Scheme 2 and Scheme 3 of the previous test are considered. Nevertheless, different configurations are obtained since α ∈ [0, ∞], as it can be seen from Figure 2.
The geometry of the sandwich plate is kept constant and three laminations schemes are considered: (0 • /core/0 • ), (30 • /core/45 • ) and (−45 • /core/45 • ). A damaged core is modeled by setting D = 0.5. It should be recalled that the same values of the mass fraction of both CNTs and fibers (respectively w C = 0.05 and w F = 0.80) are employed. The proper choice of the exponent α allows to obtain also the following extreme cases. For α = 0, the reinforcing fibers are uniformly distributed, and the Scheme 1 of the previous application is accomplished. On the other hand, if α tends to infinity, it is easy to verify that V F = 0 and the face-sheets are made of an undamaged polymer matrix enforced only by CNTs. In this circumstance, the stiffness of the structure reaches its minimum value. In terms of natural frequencies, the results are included between these two boundary cases. The variation of the first three natural frequencies is depicted in Figure 7. The following features can be observed: • Similar behaviors are obtained for the three lamination schemes under investigation. For lower values of α, the corresponding curves are detached and the natural frequencies that can be obtained assume different values depending on the fiber orientation. By increasing the exponent α, the effect of the fibers decreases since V F draws near zero and the frequencies tends asymptotically to the same value; • The initial choice of the through-the-thickness distribution of V F (Scheme 2 or Scheme 3) affects the variation of the natural frequencies. In particular, this variation is faster for Scheme 2. In fact, the slopes of the related curves are steeper, whereas the frequency variation for Scheme 3 is a little bit more gradual; • The biggest variation of frequencies is reached for lower values of α. The decrease of the value of natural frequencies for α > 20 is negligible.

Influence of the Exponent of the Through-the-Thickness Distribution of the Fiber Volume Fraction
The current application deals with the effect of the exponent α that characterizes the through- The geometry of the sandwich plate is kept constant and three laminations schemes are considered: (0°/core/0°), (30°/core/45°) and (−45°/core/45°). A damaged core is modeled by setting

Effect of the In-Plane Fiber Orientation
A damaged sandwich plate with D = 0.25 is considered in the following application. The geometric features are kept constant with respect to the previous tests, whereas the engineering constants of the face-sheets can be computed assuming w C = 0.05 and w F = 0.80. A non-uniform distribution of the reinforcing fiber is defined according to the functions that describe Scheme 3, with α = 0.5. The aim of this numerical application is to show the dependency of the natural frequencies on the in-plane fiber orientation. Therefore, several configurations are analyzed according to the value of an angular parameter θ. The variation of the first three natural frequencies is depicted in Figure 8 for various lamination schemes depending on θ. The graphs in Figure 8 prove the following results: • As expected, the orientation of the fibers affects the dynamic response of the composite structures under consideration; • If symmetric angle-ply or cross-ply laminates, as well as antisymmetric configurations, are considered, which are denoted by (θ/core/θ) and (−θ/core/θ), the extreme values of frequencies can be obtained for θ = 0 • , 45 • , 90 • . In addition, a symmetrical behavior is obtained after reaching the value of θ = 45 • ; • This regular and symmetric behavior is lost if a laminate with a general stacking sequence, such as the last two lamination schemes, is analyzed. It should be recalled that the polymer matrix is strengthened by only straight fibers if 0 C w = .
Then, the opposite situation is taken into account. The variation of 0.40, 0.80 The results are shown in Figure 9 for the first five natural frequencies. The following observations can be deduced:

Influence of the Material Properties
Finally, this last application aims to discuss the effects of the CNT mass fraction w C and the fiber mass fraction w F on the structural response. The same geometric features are considered in this circumstance, assuming D = 0.75 as damage parameters. The through-the-thickness distribution of the reinforcing fibers V F are defined by the functions used in Scheme 3 ( Figure 5) by setting α = 4. The lamination scheme is given by (30 • /core/45 • ). Firstly, the effect of w C is investigated, by keeping constant the value of w F = 0.80, in the interval w C ∈ [0, 0.40].
It should be recalled that the polymer matrix is strengthened by only straight fibers if w C = 0. Then, the opposite situation is taken into account. The variation of w F ∈ [0.40, 0.80] is studied for a fixed value of w C = 0.05. The results are shown in Figure 9 for the first five natural frequencies.
The following observations can be deduced: • The influence of the CNT mass fraction w C is greater than the corresponding variation of the fiber mass fraction w F ; • For a small increase of w C next to zero, the variation in terms of natural frequencies that can be obtained is relevant and the behavior is non-linear; • On the other hand, bigger increases of w F do not produce the same variation of natural frequencies.
The behavior is linear in this case. It should be recalled that the polymer matrix is strengthened by only straight fibers if 0 Then, the opposite situation is taken into account. The variation of 0.40, 0.80 The results are shown in Figure 9 for the first five natural frequencies. The following observations can be deduced:

Discussion on the Mode Shapes
Finally, a brief discussion on the dependency of the previous parameters on the mode shapes is presented. For this purpose, the same geometric features of the previous tests are considered. As far as the mechanical properties of the constituents are concerned, the plate is characterized by w C = 0.05 and w F = 0.80. Several configurations are analyzed in order to investigate the effect of the stacking sequence, of the damage and of the exponent of the through-the-thickness distribution of the fibers. The Scheme 3 depicted in Figure 5 is considered here, but similar results in terms of mode shapes could be obtained also with the other schemes. The cases under considerations are summarized in Table 5 for conciseness purposes. The contour plots of the first five mode shapes related to the corresponding natural frequencies are shown in Figure 10 for the seven cases introduced in Table 5. The examples analyzed in this paragraph allow to deduce the following considerations:

•
In general, the increase of the damage parameter D does not cause any variation of the mode shapes; • The mode shapes are highly affected by the orientation of the reinforcing fibers and by the stacking sequence. This aspect can be noted by comparing the same configurations in terms of D and α, but characterized by different lamination schemes (Case 1 and Case 4, for instance); • As stated in the previous paragraphs, the increase of the exponent α reduces the influence of the reinforcing straight fibers. Therefore, the anisotropic behavior of a laminate with a general stacking sequence can be decreased. For example, the mode shapes of Case 7 tend to the ones related to Case 3 for α = 12, even if they are characterized by different fiber orientations; • The presence of a thicker isotropic core is predominant in the modal amplitudes and only noticeably variations of these mechanical parameters can define some changes in the mode shapes.

Conclusions
A set of numerical investigations has been presented to describe the mechanical behavior of laminated sandwich plates with a damaged soft-core. The FE formulation has been developed by using a nine-node quadratic rectangular Lagrange element, whereas the theoretical model for laminated plates based on the Reissner-Mindlin model has been enriched by introducing the Murakami's function. As a consequence, a FE plate theory with 7 degrees of freedom per node has been presented. The external face-sheets are made of composite materials: The polymer matrix has been strengthened by randomly oriented CNTs and oriented straight fibers. Their mechanical characterization has been carried out by using a three-phase model, which include the Eshelby-Mori-Tanaka scheme and the Halpin-Tsai approach. Several parametric tests have been performed to analyze the effect of the damage, the through-the-thickness distribution of the reinforcing fibers, the orientation of the fibers, the stacking sequence, and the mechanical features. The main achievements of the paper are summarized below: • The use of the Murakami's function is required to capture the effective mechanical behavior of sandwich structures with an inner soft-core. This aspect is very important especially if a FE commercial code is employed. In fact, it should be recalled that this function is not embedded in plate/shell formulations. Therefore, the results that can be obtained in these circumstances could be inaccurate, unless a 3D-FE modelling is pursued. Nevertheless, this approach is onerous in terms of computational time and resources; • A non-uniform distribution of the fibers along the thickness of the face-sheets could be employed to model the effective distribution of the reinforcing phase that could occur during the manufacturing process or during the structural life. This research prove that the mechanical response is affected by this parameter; • A progressive damage in the core causes a corresponding decrease of the natural frequencies, which becomes faster and faster for higher values of damages. The reinforcing layers could recover this situation. If a three-phase composite material is employed to this aim, the design of such layers could be carried out taking into account two parameters, which are the mass fractions of both CNTs and fibers. Nevertheless, a small increase of the CNT mass fraction can cause a quicker and more remarkable variation of the fundamental frequency with respect to the one that could be obtained by controlling the mass fraction of the straight fibers; • The optimal structural response can be also obtained by choosing accurately the in-plane orientation of the straight fibers. The stacking sequence, in fact, affects the value of the natural frequencies, as well as of the mode shapes.
These comments should be taken into account during the analysis of the mechanical behavior of sandwich structures subjected to a progressive damage, as well as during the process manufacturing if an optimal design has to be pursued.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
The submatrices K (e) ij , for i, j = 1, . . . , 7, which appear in the element stiffness matrix K (e) defined in (36), can be computed using the following definitions. The submatrices are presented below row-by-row  ij are valid also for the well-known RM theory, for i, j = 1, . . . , 5.