Next Article in Journal
Use of Biopolymers in Mucosally-Administered Vaccinations for Respiratory Disease
Next Article in Special Issue
Rate-Dependent Cohesive Zone Model for Fracture Simulation of Soda-Lime Glass Plate
Previous Article in Journal
Two-Step Spark Plasma Sintering Process of Ultrafine Grained WC-12Co-0.2VC Cemented Carbide
Previous Article in Special Issue
Investigation of the Flow Properties of CBM Based on Stochastic Fracture Network Modeling
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

by
Michele Bacciocchi
1,2,*,
Raimondo Luciano
3,
Carmelo Majorana
4 and
Angelo Marcello Tarantino
5
1
Department of Civil, Chemical, Environmental, and Materials Engineering (DICAM), University of Bologna, Viale del Risorgimento, 40136 Bologna, Italy
2
Dipartimento di Economia, Scienze e Diritto (DESD), University of San Marino, Via Consiglio dei Sessanta, 47891 Dogana, San Marino
3
Engineering Department, University of Napoli Parthenope, Via Ammiraglio Ferdinando Acton, 80133 Napoli, Italy
4
Department of Civil, Environmental and Architectural Engineering (DICEA), University of Padova, Via F. Marzolo, 35131 Padova, Italy
5
Department of Engineering “Enzo Ferrari” (DIEF), University of Modena and Reggio Emilia, Via Vivarelli, 41125 Modena, Italy
*
Author to whom correspondence should be addressed.
Materials 2019, 12(15), 2444; https://doi.org/10.3390/ma12152444
Submission received: 27 June 2019 / Revised: 23 July 2019 / Accepted: 29 July 2019 / Published: 31 July 2019
(This article belongs to the Special Issue Advances in Structural Mechanics Modeled with FEM)

Abstract

:
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.

1. 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].
In the current paper, the FE technique is implemented in a computational code to solve the free vibration analysis of three-layered sandwich plates with a damaged soft-core and non-uniform mechanical properties. In particular, the structures are made by an isotropic core that undergoes a progressive uniform damage, which is modeled as a decay of the mechanical properties expressed in terms of engineering constants. The damage model is the one illustrated in the book by Lemaitre and Chaboche [21], and it is representative of the formation of microcracks and discontinuities within the considered medium. Since these defects are uniformly distributed and affect the central layer of the plates independently from the direction, this phenomenon is known as “isotropic damage” and it is fully described by a scalar parameter. Further details concerning the phenomenological aspects related to the damage of materials can be found in the book by Reddy and Miravete [22], whereas some structural applications which investigate the effect of damage are illustrated in [23,24,25,26,27,28,29,30,31,32,33].
The soft-core of the structures is surrounded and strengthened by two external sinks (or face-sheets), which are stiffer and thinner than the central layer. These plies are made of a polymer matrix which contains randomly oriented Carbon nanotubes (CNTs) and keeps together straight Carbon fibers. The role of the matrix in composite materials is clearly illustrated in the books by Vinson [34], Jones [35], Reddy [36], and Barbero [37]. These works should be taken into account as complete references for the analysis and modelling of composite materials. As highlighted in the papers [38,39], the matrix is reinforced at different levels. At the nano-scale, the CNTs provide the matrix with additional stiffness [40,41,42,43,44,45,46,47,48,49]. The single fiber of Carbon nanotube (CNT) is modeled as a transversely isotropic cylinder as suggested by Odegard et al. [50] and it assumed to be randomly oriented in the matrix, following the approach developed in the paper by Shi et al. [51]. The Eshelby-Mori-Tanaka scheme is applied at this level to evaluate the mechanical properties of this enriched matrix [52,53], which has isotropic features due to the random orientation of CNTs [54]. On the other hand, the reinforcing phase at the micro-scale is given by oriented straight fibers. The semi-empirical Halpin-Tsai homogenization procedure [55,56,57], which is based on the use of the Hill’s elastic moduli [58,59], is employed at this level to obtain an accurate estimation of the overall mechanical features of the composite skins. Therefore, a multi-phase approach has been developed to this aim [60,61,62,63].
The development of micromechanics theories able to evaluate the mechanical features of fiber-reinforced composites has always fascinated many researchers, as illustrated in a complete manner by Chamis and Sendeckyj [64]. For completeness purposes, it should be recalled that different methods can be found in the literature as alternatives to the Halpin-Tsai semi-empirical homogenization procedure, such as self- consistent models [65,66], variational methodologies [67,68], and techniques based on the mechanics of materials [69,70,71]. Further details can be found also in the papers [72,73].
With respect to previous contributions, the straight fibers are graded in the thickness direction. In other words, the volume fraction distribution of the reinforcing fibers is non-uniform along the thickness of the sinks. This approach is typically applied in the class of granular composites, also known as Functionally Graded (FG) materials, to characterize the gradual variation of isotropic constituents [74,75,76,77,78,79,80]. On the other hand, in this research the graded constituent is represented by the orthotropic reinforcing fibers, which can be also arbitrarily oriented in the planar direction. A power-law function is used to characterize the through-the-thickness variation of their volume fraction. FG materials are typically used as constituents in plates and shells [81,82,83,84,85,86], beams [87,88], micro- and nano-structures [89,90,91].
Thus, it should be noted that the plates are characterized by noticeable differences in terms of mechanical properties at the layer interfaces. For this reason, a suitable structural model must be introduced to describe accurately their mechanical behavior. As clearly highlighted in the papers by Carrera [92,93,94,95,96], laminated and sandwich structures are characterized by a piece-wise continuous displacement field along the thickness direction. The differences in terms of transverse deformability among the layers give rise to a change in the slope of the displacement components between two adjacent layers. This outcome is known as zig-zag effect. As specified in [97,98,99], this effect is well-captured by Layer-Wise (LW) models, in which the degrees of freedom are assumed as independent parameters 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.

2. 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 = 2 h 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 11 ( k ) , E 22 ( k ) , the shear moduli G 12 ( k ) , G 13 ( k ) , G 23 ( k ) , and the Poisson’s ratio ν 12 ( k ) . 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 11 ( k ) = E 22 ( k ) and the Poisson’s ratio ν ( k ) = ν 12 ( k ) , whereas its shear modulus G ( k ) = G 12 ( k ) = G 13 ( k ) = G 23 ( k ) is given by:
G ( k ) = E ( k ) 2 ( 1 + ν ( k ) ) .
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.

2.1. Mechanical Properties of the Face-Sheets

A multiscale approach is employed to evaluate the overall mechanical properties of the three-phase 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 C , l C , m C , n C , p C . Its density ρ C is required to compute the volume fraction of CNTs V C as follows:
V C = ( ρ C w C ρ M ρ C ρ M + 1 ) 1 ,
where ρ M represents the density of the polymer matrix, whereas w C 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 V M = 1 V C . The polymer matrix is isotropic, and it is fully characterized by its Young’s modulus E M and Poisson’s ratio ν M . In the following, its bulk modulus K M and shear modulus G M are required in order to apply the Eshelby-Mori-Tanaka approach [52,53]. The following definitions are needed for this purpose:
K M = E M 3 ( 1 2 ν M ) ,    G M = E M 2 ( 1 + ν M ) .
These quantities are noticeably affected by the presence of randomly oriented CNTs. As illustrated in the paper by Shi et al. [51], if the agglomeration of CNTs is neglected, the bulk modulus K M * and the shear modulus G M * of the enriched matrix are given by:
K M * = K M + V C ( δ C 3 K M α C ) 3 ( V M + V C α C ) ,    G M * = G M + V C ( η C 2 G M β C ) 2 ( V M + V C β C ) ,
where the following quantities, which can be computed by defining the Hill’s elastic moduli of CNTs, are introduced:
α C = 3 ( K M + G M ) + k C + l C 3 ( G M + k C ) , β C = 1 5 ( 4 G M + 2 k C + l C 3 ( G M + k C ) + 4 G M G M + p C + 2 ( G M ( 6 K M + 8 G M ) ) G M ( 3 K M + G M ) + m C ( 3 K M + 7 G M ) ) , δ C = 1 3 ( n C + 2 l C + ( 2 k C + l C ) ( 3 K M + G M l C ) G M + k C ) , η C = 1 5 ( 2 3 ( n C l C ) + 8 G M p C G M + p C + 2 ( k C l C ) ( 2 G M + l C ) 3 ( G M + k C ) + 8 m C G M ( 3 K M + 4 G M ) 3 K M ( m C + G M ) + G P M ( 7 m C + G M ) ) .
Finally, the evaluation of K M * and G M * allows to compute the Young’s modulus E M * and the Poisson’s ratio ν M * of the polymer matrix enriched by CNTs:
E M * = 9 K M * G M * 3 K M * + G M * ,        ν M * = 3 K M * 2 G M * 6 K M * + 2 G M * ,
whereas its density is given by:
ρ M * = ( ρ C ρ M ) V C + ρ M .
The mechanical properties of a single CNT fiber in terms of its Hill’s elastic moduli, as well as its density, are listed in Table 1. Such properties are valid for a single-walled Carbon nanotube with 10 as chiral index and armchair structure. Further details about the mechanical characterization of CNT can be found in [38].
In order to apply the Halpin-Tsai approach that allows to compute the overall mechanical properties of the composite given by this enriched matrix including straight Carbon fibers, the Hill’s elastic moduli of the matrix k M * , l M * , m M * , n M * , p M * are needed. For this purpose, the following definitions are introduced:
k M * = E M * 2 ( 1 + ν M * ) ( 1 2 ν M * ) ,           l M * = 2 ν M * k M * ,           m M * = ( 1 2 ν M * ) k M * , n M * = 2 ( 1 ν M * ) k M * ,           p M * = ( 1 2 ν M * ) k M * .
Analogously, the Hill’s elastic moduli of the Carbon fibers k F , l F , m F , n F , p F are required. The reinforcing fibers are assumed as transversely isotropic and their mechanical properties are given by the corresponding Young’s moduli E 11 F , E 22 F , shear modulus G 12 F and Poisson’s ratios ν 12 F , ν 23 F . Once these quantities are known, the Hill’s elastic moduli can be easily evaluated:
k F = E 22 F 2 ( 1 ν 23 F 2 ν 21 F ν 12 F ) ,           l F = 2 ν 12 F k F ,           m F = 1 ν 23 F 2 ν 21 F ν 12 F 1 + ν 23 F k F , n F = 2 ( 1 ν 23 F ) E 11 F E 22 F k F ,           p F = G 12 F ,
where ν 21 F = E 22 F ν 12 F / E 11 F . 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
V ˜ F = ( ρ F w F ρ M * ρ F ρ M * + 1 ) 1 .
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 V F = V F ( z ) = V ˜ F f ( k ) ( z ) 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:
f ( k ) = { f 1 ( k ) ( z ) = ( z z k z k + 1 z k ) α f 2 ( k ) ( z ) = ( z k + 1 z z k + 1 z k ) α
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 :
k = k M * ( k F + m M * ) V M * + k F ( k M * + m M * ) V F ( k F + m M * ) V M * + ( k M * + m M * ) V F , l = V F l F + V M * l M * + l F l M * k F k M * ( k V F k F V M * k M * ) , m = m M * 2 V F m F ( k M * + m M * ) + 2 V M * m F m M * + V M * k M * ( m F + m M * ) 2 V F m M * ( k M * + m M * ) + 2 V M * m F m M * + V M * k M * ( m F + m M * ) , n = V F n F + V M * n M * + ( l F l M * k F k M * ) 2 ( k V F k F V M * k M * ) , p = ( p F + p M * ) p M * V M * + 2 p F p M * V F ( p F + p M * ) V M * + 2 p M * V F ,
being V M * = 1 V F . The engineering constants of these layers can be evaluated according to the following definitions:
E 11 ( k ) = n l 2 k ,        E 22 ( k ) = 4 m ( k n l 2 ) k n l 2 + m n ,        ν 12 ( k ) = l 2 k ,        G 12 ( k ) = G 13 ( k ) = p ,        G 23 ( k ) = m .
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 E 11 ( k ) ( z ) ,    E 22 ( k ) ( z ) ,    ν 12 ( k ) ( z ) ,    G 12 ( k ) ( z ) ,    G 13 ( k ) ( z ) and G 23 ( k ) ( z ) , for k = 1 , 3 . Finally, the density of the composite face-sheets is given by:
ρ ( k ) = ( ρ F ρ M * ) V F + ρ M * .
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 ν 21 ( k ) = E 22 ( k ) ν 12 ( k ) / E 11 ( k ) is also required to compute the Poisson’s ratio ν 21 ( k ) . 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.

2.2. 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:
E ( k ) = ( 1 D ) E M
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.

3. 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:
U x ( e ) ( x , y , z , t ) = u x ( e ) ( x , y , t ) + z ϕ x ( e ) ( x , y , t ) + F z ( z ) ψ x ( e ) ( x , y , t ) U y ( e ) ( x , y , z , t ) = u y ( e ) ( x , y , t ) + z ϕ y ( e ) ( x , y , t ) + F z ( z ) ψ y ( e ) ( x , y , t ) U z ( e ) ( x , y , z , t ) = u z ( e ) ( x , y , t )
in which the three-dimensional displacement components are denoted by U x ( e ) , U y ( e ) , U z ( e ) . 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:
F z = ( 1 ) k 2 z z k + 1 z k z k + 1 + z k z k + 1 z k
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 ) , u y ( e ) , u z ( e ) represent the translational displacements along x , y , z , ϕ x ( e ) , ϕ y ( e ) are the rotations about the principal axes y , x respectively, whereas ψ x ( e ) , ψ y ( e ) 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.
Due to this approximation, the degrees of freedom can be written as a function of the corresponding nodal displacements u x , i ( e ) , u y , i ( e ) , u z , i ( e ) , ϕ x , i ( e ) , ϕ y , i ( e ) , ψ x , i ( e ) , ψ y , i ( e ) :
u x ( e ) ( x , y , t ) = i = 1 9 N i ( x , y ) u x , i ( e ) ( t ) = N ¯ u x ( e ) u y ( e ) ( x , y , t ) = i = 1 9 N i ( x , y ) u y , i ( e ) ( t ) = N ¯ u y ( e ) u z ( e ) ( x , y , t ) = i = 1 9 N i ( x , y ) u z , i ( e ) ( t ) = N ¯ u z ( e ) ϕ x ( e ) ( x , y , t ) = i = 1 9 N i ( x , y ) ϕ x , i ( e ) ( t ) = N ¯ ϕ x ( e ) ϕ y ( e ) ( x , y , t ) = i = 1 9 N i ( x , y ) ϕ y , i ( e ) ( t ) = N ¯ ϕ y ( e ) ψ x ( e ) ( x , y , t ) = i = 1 9 N i ( x , y ) ψ x , i ( e ) ( t ) = N ¯ ψ x ( e ) ψ y ( e ) ( x , y , t ) = i = 1 9 N i ( x , y ) ψ y , i ( e ) ( t ) = N ¯ ψ x ( e )
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:
u x ( e ) = [ u x , 1 ( e ) u x , 9 ( e ) ] T ,       u y ( e ) = [ u y , 1 ( e ) u y , 9 ( e ) ] T ,       u z ( e ) = [ u z , 1 ( e ) u z , 9 ( e ) ] T , ϕ x ( e ) = [ ϕ x , 1 ( e ) ϕ x , 9 ( e ) ] T ,       ϕ y ( e ) = [ ϕ y , 1 ( e ) ϕ y , 9 ( e ) ] T , ψ x ( e ) = [ ψ x , 1 ( e ) ψ x , 9 ( e ) ] T ,       ψ x ( e ) = [ ψ y , 1 ( e ) ψ y , 9 ( e ) ] T .
This notation is useful to define the vector u ¯ ( e ) which includes all the nodal degrees of freedom:
u ¯ ( e ) = [ u x ( e ) u y ( e ) u z ( e ) ϕ x ( e ) ϕ y ( e ) ψ x ( e ) ψ y ( e ) ] T
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:
x ( e ) = i = 1 9 N i ( ξ , η ) x i ( e ) , ​​​       y ( e ) = i = 1 9 N i ( ξ , η ) y i ( e )
in which the i -th node of the element under consideration is identified by the couple of nodal coordinates x i ( e ) , y i ( e ) . Such coordinates are included in the corresponding vectors x ( e ) , y ( e ) , which assume the following aspects:
x e = [ x 1 ( e ) x 9 ( e ) ] T ,       y e = [ y 1 ( e ) y 9 ( e ) ] T .
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:
B ξ = [ N 1 ξ N 9 ξ ] ,       B ξ = [ N 1 η N 9 η ] .
At this point, the Jacobian matrix J can be introduced as specified in [11]:
J = [ x ( e ) ξ y ( e ) ξ x ( e ) η y ( e ) η ] = [ i = 1 9 x i ( e ) N i ξ i = 1 9 y i ( e ) N i ξ i = 1 9 x i ( e ) N i η i = 1 9 y i ( e ) N i η ] = [ B ξ x ( e ) B ξ y ( e ) B η x ( e ) B η y ( e ) ] .
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:
[ B x B y ] = J 1 [ B ξ B η ]
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:
ε x x ( e ) = U x ( e ) x = u x ( e ) x + z ϕ x ( e ) x + F z ψ x ( e ) x = ε x 0 ( e ) + z k x 0 ( e ) + F z ε x 1 = B x u x ( e ) + z B x ϕ x ( e ) + F z B x ψ x ( e ) ε y y ( e ) = U y ( e ) y = u y ( e ) y + z ϕ y ( e ) y + F z ψ y ( e ) y = ε y 0 ( e ) + z k y 0 ( e ) + F z ε y 1 = B y u y ( e ) + z B y ϕ y ( e ) + F z B y ψ y ( e ) γ x y ( e ) = U y ( e ) x + U x ( e ) y = u y ( e ) x + u x ( e ) y + z ( ϕ y ( e ) x + ϕ x ( e ) y ) + F z ( ψ y ( e ) x + ψ x ( e ) y ) =           = γ x y 0 ( e ) + z k x y 0 ( e ) + F z γ x y 1 ( e ) = B x u y ( e ) + B y u x ( e ) + z ( B x ϕ y ( e ) + B y ϕ x ( e ) ) + F z ( B x ψ y ( e ) + B y ψ x ( e ) ) γ x z ( e ) = U z ( e ) x + U x ( e ) z = u z ( e ) x + ϕ x ( e ) + F z z ψ x ( e ) = γ x z 0 ( e ) + F z z γ x z 1 ( e ) = B x u z ( e ) + N ¯ ϕ x ( e ) + F z z N ¯ ψ x ( e ) γ y z ( e ) = U z ( e ) y + U y ( e ) z = u z ( e ) y + ϕ y ( e ) + F z z ψ y ( e ) = γ y z 0 ( e ) + F z z γ y z 1 ( e ) = B y u z ( e ) + N ¯ ϕ y ( e ) + F z z N ¯ ψ y ( e ) ,
where ε x x ( e ) , ε y y ( e ) , γ x y ( e ) are the membrane strains, and γ x z ( e ) , γ y z ( e ) are the transverse shear strains. It should be observed that the normal strain along the thickness direction ε z z ( e ) is omitted since the plane-strain assumption entails that ε z z ( e ) = 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 ε x 0 ( e ) , ε y 0 ( e ) , γ x y 0 ( e ) are the well-known membrane strains, k x 0 ( e ) , k y 0 ( e ) , k x y 0 ( e ) the bending and twisting curvatures, and γ x z 0 ( e ) , γ y z 0 ( e ) 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 σ z z k ( e ) = 0 is assumed by hypothesis, whereas the other stress components are given by:
σ x x k ( e ) = Q ¯ 11 ( k ) ε x x ( e ) + Q ¯ 12 ( k ) ε y y ( e ) + Q ¯ 16 ( k ) γ x y ( e ) σ y y k ( e ) = Q ¯ 12 ( k ) ε x x ( e ) + Q ¯ 22 ( k ) ε y y ( e ) + Q ¯ 26 ( k ) γ x y ( e ) τ x y k ( e ) = Q ¯ 16 ( k ) ε x x ( e ) + Q ¯ 26 ( k ) ε y y ( e ) + Q ¯ 66 ( k ) γ x y ( e ) τ x z k ( e ) = Q ¯ 44 ( k ) γ x z ( e ) + Q ¯ 45 ( k ) γ y z ( e ) τ y z k ( e ) = Q ¯ 45 ( k ) γ x z ( e ) + Q ¯ 55 ( k ) γ y z ( e )
in which Q ¯ i j ( k ) 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 i j ( k ) expressed as a function of the engineering constants of the k -layer defined below, can be found in the book by Reddy [36]:
Q 11 ( k ) = E 11 ( k ) 1 ν 12 ( k ) ν 21 ( k ) ,     Q 22 ( k ) = E 22 ( k ) 1 ν 12 ( k ) ν 21 ( k ) ,     Q 12 ( k ) = ν 12 ( k ) E 22 ( k ) 1 ν 12 ( k ) ν 21 ( k ) ,     Q 66 ( k ) = G 12 ( k ) ,     Q 44 ( k ) = G 13 ( k ) ,     Q 55 ( k ) = G 23 ( k ) .
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:
t 1 t 2 ( δ Κ ( e ) δ Φ ( e ) ) d t = 0
where δ Κ ( 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 δ Κ ( e ) is given by:
δ Κ ( e ) = k = 1 3 x y z k z k + 1 ρ ( k ) ( δ U x ( e ) U ¨ x ( e ) + δ U y ( e ) U ¨ y ( e ) + δ U y ( e ) U ¨ y ( e ) ) d x d y d z
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:
δ Φ ( e ) = k = 1 3 x y z k z k + 1 ( σ x x k ( e ) δ ε x x ( e ) + σ y y k ( e ) δ ε y y ( e ) + τ x y k ( e ) δ γ x y ( e ) + τ x z k ( e ) δ γ x z ( e ) + τ y z k ( e ) δ γ y z ( e ) ) d x d y d z
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:
[ N x 0 ( e ) N y 0 ( e ) N x y 0 ( e ) M x 0 ( e ) M y 0 ( e ) M x y 0 ( e ) T x 0 ( e ) T y 0 ( e ) N x 1 ( e ) N y 1 ( e ) N x y 1 ( e ) T x 1 ( e ) T y 1 ( e ) ] = [ A 11 A 12 A 16 B 11 B 12 B 16 0 0 F 11 F 12 F 16 0 0 A 12 A 22 A 26 B 12 B 22 B 26 0 0 F 12 F 22 F 26 0 0 A 16 A 26 A 66 B 16 B 26 B 66 0 0 F 16 F 26 F 66 0 0 B 11 B 12 B 16 D 11 D 12 D 16 0 0 G 11 G 12 G 16 0 0 B 12 B 22 B 26 D 12 D 22 D 26 0 0 G 12 G 22 G 26 0 0 B 16 B 26 B 66 D 16 D 26 D 66 0 0 G 16 G 26 G 66 0 0 0 0 0 0 0 0 κ A 44 κ A 45 0 0 0 κ H 44 κ H 45 0 0 0 0 0 0 κ A 45 κ A 55 0 0 0 κ H 45 κ H 55 F 11 F 12 F 16 G 11 G 12 G 16 0 0 L 11 L 12 L 16 0 0 F 12 F 22 F 26 G 12 G 22 G 26 0 0 L 12 L 22 L 26 0 0 F 16 F 26 F 66 G 16 G 26 G 66 0 0 L 16 L 26 L 66 0 0 0 0 0 0 0 0 κ H 44 κ H 45 0 0 0 κ P 44 κ P 45 0 0 0 0 0 0 κ H 45 κ H 55 0 0 0 κ P 45 κ P 55 ] [ ε x 0 ( e ) ε y 0 ( e ) γ x y 0 ( e ) k x 0 ( e ) k y 0 ( e ) k x y 0 ( e ) γ x z 0 ( e ) γ y z 0 ( e ) ε x 1 ( e ) ε y 1 ( e ) γ x y 1 ( e ) γ x z 1 ( e ) γ y z 1 ( e ) ] .
The superscripts “ 0 ” and “ 1 ” have the same meaning discussed previously. In particular, N x 0 ( e ) , N y 0 ( e ) , N x y 0 ( e ) are the membrane forces, M x 0 ( e ) , M y 0 ( e ) , M x y 0 ( e ) the bending and twisting moments, and T x 0 ( e ) , T y 0 ( e ) 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]:
A i j = k = 1 3 z k z k + 1 Q ¯ i j ( k ) d z ,        B i j = k = 1 3 z k z k + 1 z Q ¯ i j ( k ) d z ,        D i j = k = 1 3 z k z k + 1 z 2 Q ¯ i j ( k ) d z .
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:
F i j = k = 1 3 z k z k + 1 F z Q ¯ i j ( k ) d z ,        G i j = k = 1 3 z k z k + 1 z F z Q ¯ i j ( k ) d z ,        H i j = k = 1 3 z k z k + 1 F z z Q ¯ i j ( k ) d z , L i j = k = 1 3 z k z k + 1 F z 2 Q ¯ i j ( k ) d z ,        P i j = k = 1 3 z k z k + 1 ( F z z ) 2 Q ¯ i j ( k ) d z .
It should be specified that the integrals in (33)–(34) are computed numerically since the stiffnesses Q ¯ i j ( k ) 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:
K ( e ) u ¯ ( e ) + M ( e ) u ¯ ¨ ( e ) = 0
where K ( e ) is the element stiffness matrix, M ( e ) the element mass matrix, and u ¯ ¨ ( e ) 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:
K ( e ) = [ K 11 K 17 K 71 K 77 ] .
On the other hand, the mass matrix M ( e ) is given by:
M ( e ) = [ M 11 M 17 M 71 M 77 ] .
The matrices K i j and M i j , for i , j = 1 , , 7 are defined in the Appendix A. In particular, the following inertia terms are required to compute the mass matrix:
I 0 = k = 1 3 z k z k + 1 ρ ( k ) d z ,        I 1 = k = 1 3 z k z k + 1 z ρ ( k ) d z ,        I 2 = k = 1 3 z k z k + 1 z 2 ρ ( k ) d z , I 3 = k = 1 3 z k z k + 1 F z ρ ( k ) d z ,        I 4 = k = 1 3 z k z k + 1 z F z ρ ( k ) d z ,        I 5 = k = 1 3 z k z k + 1 F z 2 ρ ( k ) d z .
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:
K u + M u ¨ = 0
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.
If N P denotes the number of nodes, the model is characterized by N d o f s = 7 N 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 2 N P u z , 1       u z , N P 2 N P + 1 3 N P                            ϕ x , 1       ϕ x , N P 3 N P + 1 4 N P ϕ y , 1       ϕ y , N P 4 N P + 1 5 N P ψ x , 1       ψ x , N P 5 N P + 1 6 N P ψ y , 1       ψ y , N P 6 N P + 1 7 N P ] T .
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 K , M is given by N d o f s × N d o 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.

3.1. 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:
x y G ( x , y ) d x d y = 1 1 1 1 G ( ξ , η ) det J d ξ d η
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 ξ , η . From the numerical point of view, this integral can be converted into the following weighted linear sum:
1 1 1 1 G ( ξ , η ) det J d ξ d η I = 1 M J = 1 N G ( ξ I , η J ) det J | ξ I , η J W I W J
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].

3.2. 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:
( K ω 2 M ) d = 0
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 π .

4. 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.

4.1. 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.

4.2. 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 11 F = 5.6466 TPa , E 22 F = 7.0800 TPa , G 12 F = 1.9445 TPa , ν 12 F = 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 1 ( 1 ) , f 2 ( 2 ) with α = 1 , whereas “Case 2” is given by f 2 ( 1 ) , f 1 ( 2 ) with α = 1 . The frequencies are presented in Table 4 in dimensionless form as
ϖ = ω L x 2 h ρ M E M
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.

4.3. 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 2 ( 1 ) = f 2 ( 3 ) with α = 1 ; “Scheme 3” is accomplished by using f 2 ( 1 ) , f 1 ( 3 ) 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 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.

4.4. 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.

4.5. 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.

4.6. 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.

4.7. 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.

5. 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.

Author Contributions

Conceptualization, M.B., R.L., C.M. and A.M.T.; methodology, M.B., R.L., C.M. and A.M.T.; software, M.B., R.L., C.M. and A.M.T.; validation, M.B., R.L., C.M. and A.M.T.; formal analysis, M.B., R.L., C.M. and A.M.T.; investigation, M.B., R.L., C.M. and A.M.T.; resources, M.B., R.L., C.M. and A.M.T.; data curation, M.B., R.L., C.M. and A.M.T.; writing—original draft preparation, M.B., R.L., C.M. and A.M.T.; writing—review and editing, M.B., R.L., C.M. and A.M.T.; visualization, M.B., R.L., C.M. and A.M.T.; supervision, M.B., R.L., C.M. and A.M.T.; project administration, M.B., R.L., C.M. and A.M.T.

Funding

This research received no external funding.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

The submatrices K i j ( e ) , 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
K 11 = x y ( B x T ( A 11 B x + A 16 B y ) + B y T ( A 16 B x + A 66 B y ) ) d x d y K 12 = x y ( B x T ( A 12 B y + A 16 B x ) + B y T ( A 26 B y + A 66 B x ) ) d x d y K 13 = 0 K 14 = x y ( B x T ( B 11 B x + B 16 B y ) + B y T ( B 16 B x + B 66 B y ) ) d x d y K 15 = x y ( B x T ( B 12 B y + B 16 B x ) + B y T ( B 26 B y + B 66 B x ) ) d x d y K 16 = x y ( B x T ( F 11 B x + F 16 B y ) + B y T ( F 16 B x + F 66 B y ) ) d x d y K 17 = x y ( B x T ( F 12 B y + F 16 B x ) + B y T ( F 26 B y + F 66 B x ) ) d x d y
K 21 = K 12 T K 22 = x y ( B y T ( A 22 B y + A 26 B x ) + B x T ( A 26 B y + A 66 B x ) ) d x d y K 23 = 0 K 24 = x y ( B y T ( B 12 B x + B 26 B y ) + B x T ( B 16 B x + B 66 B y ) ) d x d y K 25 = x y ( B y T ( B 22 B y + B 26 B x ) + B x T ( B 26 B y + B 66 B x ) ) d x d y K 26 = x y ( B y T ( F 12 B x + F 26 B y ) + B x T ( F 16 B x + F 66 B y ) ) d x d y K 27 = x y ( B y T ( F 22 B y + F 26 B x ) + B x T ( F 26 B y + F 66 B x ) ) d x d y
K 31 = K 13 T ,         K 32 = K 23 T K 33 = x y ( B x T ( κ A 44 B x + κ A 45 B y ) + B y T ( κ A 45 B x + κ A 55 B y ) ) d x d y K 34 = x y ( B x T ( κ A 44 N ¯ ) + B y T ( κ A 45 N ¯ ) ) d x d y K 35 = x y ( B x T ( κ A 45 N ¯ ) + B y T ( κ A 55 N ¯ ) ) d x d y K 36 = x y ( B x T ( κ H 44 N ¯ ) + B y T ( κ H 45 N ¯ ) ) d x d y K 37 = x y ( B x T ( κ H 45 N ¯ ) + B y T ( κ H 55 N ¯ ) ) d x d y
K 41 = K 14 T ,         K 42 = K 24 T ,         K 43 = K 34 T K 44 = x y ( B x T ( D 11 B x e + D 16 B y ) + B y T ( D 16 B x + D 66 B y ) ) d x d y + x y N ¯ T κ A 44 N ¯ d x d y K 45 = x y ( B x T ( D 12 B y + D 16 B x ) + B y T ( D 26 B y + D 66 B x ) ) d x d y + x y N ¯ T κ A 45 N ¯ d x d y K 46 = x y ( B x T ( G 11 B x + G 16 B y ) + B y T ( G 16 B x + G 66 B y ) ) d x d y + x y N ¯ T κ H 44 N ¯ d x d y K 47 = x y ( B x T ( G 12 B y + G 16 B x ) + B y T ( G 26 B y + G 66 B x ) ) d x d y + x y N ¯ T κ H 45 N ¯ d x d y
K 51 = K 15 T ,         K 52 = K 25 T ,         K 53 = K 35 T ,         K 54 = K 45 T K 55 = x y ( B y T ( D 22 B y + D 26 B x ) + B x T ( D 26 B y + D 66 B x ) ) d x d y + x y N ¯ T κ A 55 N ¯ d x d y K 56 = x y ( B y T ( G 12 B x + G 26 B y ) + B x T ( G 16 B x + G 66 B y ) ) d x d y + x y N ¯ T κ H 45 N ¯ d x d y K 57 = x y ( B y T ( G 22 B y + G 26 B x ) + B x T ( G 26 B y + G 66 B x ) ) d x d y + x y N ¯ T κ H 55 N ¯ d x d y
K 61 = K 16 T ,         K 62 = K 26 T ,         K 63 = K 36 T ,         K 64 = K 46 T ,         K 65 = K 56 T K 66 = x y ( B x T ( L 11 B x + L 16 B y ) + B y T ( L 16 B x + L 66 B y ) ) d x d y + x y N ¯ T κ P 44 N ¯ d x d y K 67 = x y ( B x T ( L 12 B y + L 16 B x ) + B y T ( L 26 B y + L 66 B x ) ) d x d y + x y N ¯ T κ P 45 N ¯ d x d y
K 71 = K 17 T ,         K 72 = K 27 T ,         K 73 = K 37 T ,         K 74 = K 47 T ,         K 75 = K 57 T ,         K 76 = K 67 T K 77 = x y ( B y T ( L 22 B y + L 26 B x ) + B x T ( L 26 B y + L 66 B x ) ) d x d y + x y N ¯ T κ P 55 N ¯ d x d y
Analogously, the submatrices M i j ( e ) , for i , j = 1 , , 7 , which appear in the element mass matrix M ( e ) introduced in (37), are defined below. The submatrices are presented row-by-row as in the previous case
M 11 = x y N ¯ T I 0 N ¯ d x d y ,         M 14 = x y N ¯ T I 1 N ¯ d x d y ,         M 16 = x y N ¯ T I 3 N ¯ d x d y
M 22 = x y N ¯ T I 0 N ¯ d x d y ,         M 25 = x y N ¯ T I 1 N ¯ d x d y ,         M 27 = x y N ¯ T I 3 N ¯ d x d y
M 33 = x y N ¯ T I 0 N ¯ d x d y
M 41 = M 14 T ,         M 44 = x y N ¯ T I 2 N ¯ d x d y ,         M 46 = x y N ¯ T I 4 N ¯ d x d y
M 52 = M 25 T ,         M 55 = x y N ¯ T I 2 N ¯ d x d y ,         M 57 = x y N ¯ T I 4 N ¯ d x d y
M 61 = M 16 T ,         M 64 = M 46 T ,         M 66 = x y N ¯ T I 5 N ¯ d x d y
M 72 = M 27 T ,         M 75 = M 57 T ,         M 77 = x y N ¯ T I 5 N ¯ d x d y
The null submatrices ( M i j ( e ) = 0 ) are omitted in definitions (A8)–(A14). The full integration scheme is always used to compute the elements of the mass matrix.
It can be easily noted from definitions (A1)–(A14) that the stiffness and mass matrices are both symmetrical operators, thus K i j ( e ) = K j i ( e ) ,    M i j ( e ) = M j i ( e ) . It should be recalled that the same definitions K i j ( e ) ,    M i j ( e ) are valid also for the well-known RM theory, for i , j = 1 , , 5 .

References

  1. Duncan, W.J.; Collar, A.R. A method for the solution of oscillations problems by matrices. Phil. Mag. 1934, 17, 865–909. [Google Scholar] [CrossRef]
  2. Duncan, W.J.; Collar, A.R. Matrices applied to the motions of damped systems. Phil. Mag. 1935, 19, 197–219. [Google Scholar] [CrossRef]
  3. Hrennikoff, A. Solution of Problems of Elasticity by the Frame–Work Method. ASME J. Appl. Mech. 1941, 8, A619–A715. [Google Scholar]
  4. Courant, R. Variational methods for the solution of problems of equilibrium and vibration. B. Am. Math. Soc. 1943, 49, 1–23. [Google Scholar] [CrossRef]
  5. Clough, R.W. The finite element method in plane stress analysis. In Proceedings of the 2nd ASCE conference in electronics computation, Pittsburgh, PA, USA, 8–9 September 1960. [Google Scholar]
  6. Melosh, R.J. Basis for derivation of matrices for the direct stiffness method. AIAA J. 1963, 1, 1631–1637. [Google Scholar] [CrossRef]
  7. Oden, J.T. Finite Elements of Nonlinear Continua; McGraw-Hill: New York, NY, USA, 1972. [Google Scholar]
  8. Oden, J.T.; Reddy, J.N. An Introduction to the Mathematical Theory of Finite Elements; John Wiley: New York, NY, USA, 1976. [Google Scholar]
  9. Hinton, E. Numerical Methods and Software for Dynamic Analysis of Plates and Shells; Pineridge Press: Swansea, UK, 1988. [Google Scholar]
  10. Zienkiewicz, O.C. The Finite Element Method; McGraw–Hill: New York, NY, USA, 1991. [Google Scholar]
  11. Reddy, J.N. An Introduction to the Finite Element Method; McGraw–Hill: New York, NY, USA, 1993. [Google Scholar]
  12. Hughes, T.J.R. The Finite Element Method-Linear Static and Dynamic Finite Element Analysis; Dover Publications: New York, NY, USA, 2000. [Google Scholar]
  13. Ferreira, A.J.M. MATLAB Codes for Finite Element Analysis; Springer: New York, NY, USA, 2008. [Google Scholar]
  14. Martínez-Pañeda, E. On the finite element implementation of functionally graded materials. Materials 2019, 12, 287. [Google Scholar] [CrossRef]
  15. Nguyen, H.N.; Nguyen, T.Y.; Tran, K.V.; Tran, T.T.; Nguyen, T.T.; Phan, V.D.; Do, T.V. A finite element model for dynamic analysis of triple-layer composite plates with layers connected by shear connectors subjected to moving load. Materials 2019, 12, 598. [Google Scholar] [CrossRef] [PubMed]
  16. Leonetti, L.; Fantuzzi, N.; Trovalusci, P.; Tornabene, F. Scale effects in orthotropic composite assemblies as micropolar continua: A comparison between weak-and strong-form finite element solutions. Materials 2019, 12, 758. [Google Scholar] [CrossRef]
  17. Liu, P.; Bui, T.Q.; Zhu, D.; Yu, T.T.; Wang, J.W.; Yin, S.H.; Hirose, S. Buckling failure analysis of cracked functionally graded plates by a stabilized discrete shear gap extended 3-node triangular plate element. Compos. Part B Eng. 2015, 77, 179–193. [Google Scholar] [CrossRef]
  18. Hosseini, S.S.; Bayesteh, H.; Mohammadi, S. Thermo-mechanical XFEM crack propagation analysis of functionally graded materials. Mat. Sci. Eng. A 2013, 561, 285–302. [Google Scholar] [CrossRef]
  19. Yin, S.; Yu, T.; Bui, T.Q.; Liu, P.; Hirose, S. Buckling and vibration extended isogeometric analysis of imperfect graded Reissner-Mindlin plates with internal defects using NURBS and level sets. Comput. Struct. 2016, 177, 23–38. [Google Scholar] [CrossRef]
  20. Singh, S.K.; Singh, I.V.; Mishra, B.K.; Bhardwaj, G.; Singh, S.K. Analysis of cracked plate using higher-order shear deformation theory: Asymptotic crack-tip fields and XIGA implementation. Comput. Method. Appl. M. 2018, 336, 594–639. [Google Scholar] [CrossRef]
  21. Lemaitre, J.; Chaboche, J.L. Mechanics of Solid Materials; Cambridge University Press: New York, NY, USA, 1990. [Google Scholar]
  22. Reddy, J.N.; Miravete, A. Practical Analysis of Composite Laminates; CRC Press: Boca Raton, FL, USA, 1995. [Google Scholar]
  23. Tarantino, A.M. Equilibrium paths of a hyperelastic body under progressive damage. J. Elast. 2014, 114, 225–250. [Google Scholar] [CrossRef]
  24. Lanzoni, L.; Tarantino, A.M. Damaged hyperelastic membranes. Int. J. Nonlinear Mech. 2014, 60, 9–22. [Google Scholar] [CrossRef]
  25. Lanzoni, L.; Tarantino, A.M. Equilibrium configurations and stability of a damaged body under uniaxial tractions. Z. Angew. Math. Phys. 2015, 66, 171–190. [Google Scholar] [CrossRef]
  26. Lanzoni, L.; Tarantino, A.M. A simple nonlinear model to simulate the localized necking and neck propagation. Int. J. Nonlinear Mech. 2016, 84, 94–104. [Google Scholar] [CrossRef]
  27. Savino, V.; Lanzoni, L.; Tarantino, A.M.; Viviani, M. Simple and effective models to predict the compressive and tensile strength of HPFRC as the steel fiber content and type changes. Compos. Part B Eng. 2018, 137, 153–162. [Google Scholar] [CrossRef] [Green Version]
  28. Falope, F.O.; Lanzoni, L.; Tarantino, A.M. Modified hinged beam test on steel fabric reinforced cementitious matrix (SFRCM). Compos. Part B Eng. 2018, 146, 232–243. [Google Scholar] [CrossRef]
  29. Dezi, L.; Menditto, G.; Tarantino, A.M. Homogeneous structures subjected to successive structural system changes. J. Eng. Mech. ASCE 1990, 116, 1723–1732. [Google Scholar] [CrossRef]
  30. Dezi, L.; Tarantino, A.M. Time dependent analysis of concrete structures with variable structural system. ACI Mater. J. 1991, 88, 320–324. [Google Scholar]
  31. Dezi, L.; Menditto, G.; Tarantino, A.M. Viscoelastic heterogeneous structures with variable structural system. J. Eng. Mech. ASCE 1992, 119, 238–250. [Google Scholar] [CrossRef]
  32. Dezi, L.; Tarantino, A.M. Creep in continuous composite beams. Part II: Parametric study. J. Eng. Mech. ASCE 1993, 119, 2112–2133. [Google Scholar]
  33. Tarantino, A.M. Homogeneous equilibrium configurations of a hyperelastic compressible cube under equitriaxial dead-load tractions. J. Elast. 2008, 92, 227–254. [Google Scholar] [CrossRef]
  34. Vinson, J.R. The Behavior of Shells Composed of Isotropic and Composite Materials; Springer: New York, NY, USA, 1993. [Google Scholar]
  35. Jones, R.M. Mechanics of Composite Materials, 2nd ed.; Taylor & Francis: Philadelphia, PA, USA, 1999. [Google Scholar]
  36. Reddy, J.N. Mechanics of Laminated Composite Plates and Shells-Theory and Analysis, 2nd ed.; CRC Press: Boca Raton, FL, USA, 2004. [Google Scholar]
  37. Barbero, E.J. Introduction to Composite Materials Design; CRC Press: Boca Raton, FL, USA, 2011. [Google Scholar]
  38. Tornabene, F.; Bacciocchi, M.; Fantuzzi, N.; Reddy, J.N. Multiscale approach for three-phase cnt/polymer/ fiber laminated nanocomposite structures. Polym. Compos. 2019, 40, E102–E126. [Google Scholar] [CrossRef]
  39. Bacciocchi, M.; Tarantino, A.M. Time-dependent behavior of viscoelastic three-phase composite plates reinforced by carbon nanotubes. Compos. Struct. 2019, 216, 20–31. [Google Scholar] [CrossRef]
  40. Popov, V.N.; Van Doren, V.E. Elastic properties of single-walled carbon nanotubes. Phys. Rev. B 2000, 61, 3078–3084. [Google Scholar] [CrossRef]
  41. Qian, D.; Wagner, G.J.; Liu, W.K.; Yu, M.F.; Ruoff, R.S. Mechanics of carbon nanotubes. Appl. Mech. Rev. 2002, 55, 495–533. [Google Scholar] [CrossRef]
  42. Fidelus, J.D.; Wiesel, E.; Gojny, F.H.; Schulte, K.; Wagner, H.D. Thermo–mechanical properties of randomly oriented carbon/epoxy nanocomposites. Compos. Part A Appl. S. 2005, 36, 1555–1561. [Google Scholar] [CrossRef]
  43. Ray, M.C.; Batra, R.C. Effective Properties of Carbon Nanotube and Piezoelectric Fiber Reinforced Hybrid Smart Composites. J. App. Mech. T. ASME 2009, 76, 034503. [Google Scholar] [CrossRef]
  44. Song, Y.S.; Youn, J.R. Modeling of effective elastic properties for polymer based carbon nanotube composites. Polymer 2006, 47, 1741–1748. [Google Scholar] [CrossRef]
  45. Coiai, S.; Passaglia, E.; Pucci, A.; Ruggeri, G. Nanocomposites based on thermoplastic polymers and functional nanofiller for sensor applications. Materials 2015, 8, 3377–3427. [Google Scholar] [CrossRef]
  46. Bhattacharya, M. Polymer Nanocomposites-a comparison between carbon nanotubes, graphene, and clay as nanofillers. Materials 2016, 9, 262. [Google Scholar] [CrossRef] [PubMed]
  47. Acierno, S.; Barretta, R.; Luciano, R.; Marotti de Sciarra, F.; Russo, P. Experimental evaluations and modeling of the tensile behavior of polypropylene/single-walled carbon nanotubes fibers. Compos. Struct. 2017, 174, 12–18. [Google Scholar] [CrossRef]
  48. Wang, G.; Wang, Y.; Luo, Y.; Luo, S. Carbon nanomaterials based smart fabrics with selectable characteristics for in–line monitoring of high-performance composites. Materials 2018, 11, 1677. [Google Scholar] [CrossRef] [PubMed]
  49. Arena, M.; Viscardi, M.; Barra, G.; Vertuccio, L.; Guadagno, L. Multifunctional performance of a nano–modified fiber reinforced composite aeronautical panel. Materials 2019, 12, 869. [Google Scholar] [CrossRef] [PubMed]
  50. Odegard, G.M.; Gates, T.S.; Wise, K.E.; Park, C.; Siochi, E.J. Constitutive modeling of nanotube–reinforced polymer composites. Compos. Sci. Technol. 2003, 63, 1671–1687. [Google Scholar] [CrossRef]
  51. Shi, D.L.; Huang, Y.Y.; Hwang, K.C.; Gao, H. The effect of nanotube waviness and agglomeration on the elastic property of carbon nanotube–reinforced composites. J. Eng. Mater. T. ASME 2004, 126, 250–257. [Google Scholar] [CrossRef]
  52. Eshelby, J.D. The determination of the elastic field of an ellipsoidal inclusion, and related problems. P. Roy. Soc. Lond. A Mat. 1957, 241, 376–396. [Google Scholar]
  53. Mori, T.; Tanaka, K. Average stress in matrix and average elastic energy of materials with misfitting inclusions. Acta Metall. 1973, 21, 571–574. [Google Scholar] [CrossRef]
  54. Safaei, B.; Moradi-Dastjerdi, R.; Qin, Z.; Behdinan, K.; Chu, F. Determination of thermoelastic stress wave propagation in nanocomposite sandwich plates reinforced by clusters of carbon nanotubes. J. Sandw. Struct. Mater. 2019. [Google Scholar] [CrossRef]
  55. Halpin, J.C. Effects of Environmental Factors on Composite Materials; Technical Report AFML-TR-67-423; Air Force Materials Lab Wright-Patterson AFB: Greene, OH, USA, 1969. [Google Scholar]
  56. Tsai, S.W. Structural Behavior of Composite Materials; Philco Corporation: Newport Beach, CA, USA, 1964. [Google Scholar]
  57. Tsai, S.W. Strength Characteristics of Composite Materials; Philco Corporation: Newport Beach, CA, USA, 1965. [Google Scholar]
  58. Hill, R. Theory of mechanical properties of fibre-strengthened materials: I. Elastic behavior. J. Mech. Phys. Solids 1964, 12, 199–212. [Google Scholar] [CrossRef]
  59. Hill, R. Theory of mechanical properties of fibre–strengthened materials: II. Inelastic behavior. J. Mech. Phys. Solids 1964, 12, 213–218. [Google Scholar] [CrossRef]
  60. Thostenson, E.T.; Li, W.Z.; Wang, D.Z.; Ren, Z.F.; Chou, T.W. Carbon nanotube/carbon fiber hybrid multiscale composites. J. Appl. Phys. 2002, 91, 6034–6037. [Google Scholar] [CrossRef]
  61. Bekyarova, E.; Thostenson, E.T.; Yu, A.; Kim, H.; Gao, J.; Tang, J.; Hahn, H.T.; Chou, T.W.; Itkis, M.E.; Haddon, R.C. Multiscale carbon nanotube-carbon fiber reinforcement for advanced epoxy composites. Langmuir 2007, 23, 3970–3974. [Google Scholar] [CrossRef] [PubMed]
  62. Kim, M.; Park, Y.B.; Okoli, O.I.; Zhang, C. Processing, characterization, and modeling of carbon nanotube-reinforced multiscale composites. Compos. Sci. Technol. 2009, 69, 335–342. [Google Scholar] [CrossRef]
  63. Rafiee, M.; He, X.Q.; Mareishi, S.; Liew, K.M. Modeling and stress analysis of smart CNTs/fiber/polymer multiscale composite plates. Int. J. Appl. Mech. 2014, 6, 1450025. [Google Scholar] [CrossRef]
  64. Chamis, C.C.; Sendeckyj, G.P. Critique on theories predicting thermoelastic properties of fibrous composites. J. Compos. Mater. 1968, 2, 332–358. [Google Scholar] [CrossRef]
  65. Hill, R. Theory of mechanical properties of fibre–strengthened materials: III. Self–consistent model. J. Mech. Phys. Solids 1964, 13, 189–198. [Google Scholar] [CrossRef]
  66. Chou, T.W. A self-consistent approach to the elastic stiffness of short–fiber composites. J. Compos. Mater. 1980, 14, 178–188. [Google Scholar] [CrossRef]
  67. Hashin, Z. The elastic moduli of heterogeneous materials. J. Appl. Mech. T. ASME 1962, 29, 143–150. [Google Scholar] [CrossRef]
  68. Hashin, Z.; Rosen, B.W. The elastic moduli of fiber-reinforced materials. J. Appl. Mech. T. ASME 1964, 31, 223–232. [Google Scholar] [CrossRef]
  69. Ekvall, J.C. Elastic Properties of Orthotropic Monofilament Laminates; Lockheed Aircraft Corporation: Burbank, CA, USA, 1961; 61-AV-56. [Google Scholar]
  70. Chen, C.H.; Cheng, S. Mechanical properties of fiber reinforced composites. J. Compos. Mater. 1967, 1, 30–41. [Google Scholar] [CrossRef]
  71. Lei, Z.X.; Liew, K.M.; Yu, J.L. Free vibration analysis of functionally graded carbon nanotube-reinforced composite plates using the element-free kp-Ritz method in thermal environment. Compos Struct. 2013, 106, 128–138. [Google Scholar] [CrossRef]
  72. Gusella, F.; Cluni, F.; Gusella, V. Homogenization of dynamic behaviour of heterogeneous beams with random Young’s modulus. Eur. J. Mech. A Solid. 2019, 73, 260–267. [Google Scholar] [CrossRef]
  73. Gusella, F.; Cluni, F.; Gusella, V. Homogenization of the heterogeneous beam dynamics: The influence of the random Young’s modulus mixing law. Compos. Part B Eng 2019, 167, 608–614. [Google Scholar] [CrossRef]
  74. Reddy, J.N.; Chin, C.D. Thermomechanical analysis of functionally graded cylinders and plates. J. Therm. Stresses 1998, 21, 593–626. [Google Scholar] [CrossRef]
  75. Reddy, J.N. Analysis of functionally graded plates. Int. J. Numer. Meth. Engng. 2000, 47, 663–684. [Google Scholar] [CrossRef]
  76. Vel, S.S.; Batra, R.C. Three-dimensional exact solution for the vibration of functionally graded rectangular plates. J. Sound Vib. 2004, 272, 703–730. [Google Scholar] [CrossRef]
  77. Batra, R.C.; Jin, J. Natural frequencies of a functionally graded rectangular plate. J. Sound Vib. 2005, 282, 509–516. [Google Scholar] [CrossRef]
  78. Kim, J.; Reddy, J.N. A general third-order theory of functionally graded plates with modified couple stress effect and the von Kármán nonlinearity: Theory and finite element analysis. Acta Mech. 2015, 226, 2973–2998. [Google Scholar] [CrossRef]
  79. Kim, J.; Reddy, J.N. Modeling of functionally graded smart plates with gradient elasticity effects. Mech. Adv. Mater. Struct. 2017, 24, 437–447. [Google Scholar] [CrossRef]
  80. Alexandrov, S.; Wang, Y.C.; Lang, L. A theory of elastic/plastic plane strain pure bending of FGM sheets at large strain. Materials 2019, 12, 456. [Google Scholar] [CrossRef] [PubMed]
  81. Tornabene, F. Free vibration analysis of functionally graded conical, cylindrical shell and annular plate structures with a four-parameter power-law distribution. Comput. Method. Appl. Mech. Eng. 2009, 198, 2911–2935. [Google Scholar] [CrossRef]
  82. Tornabene, F.; Viola, E. Free vibration analysis of functionally graded panels and shells of revolution. Meccanica 2009, 44, 255–281. [Google Scholar] [CrossRef]
  83. Sofiyev, A.H.; Kuruoglu, N. Dynamic instability of three-layered cylindrical shells containing an FGM interlayer. Thin Wall. Struct. 2015, 93, 10–21. [Google Scholar] [CrossRef]
  84. Alibeigloo, A. Thermo elasticity solution of sandwich circular plate with functionally graded core using generalized differential quadrature method. Compos. Struct. 2016, 136, 229–240. [Google Scholar] [CrossRef]
  85. Civalek, Ö.; Baltacıoglu, A.K. Free vibration analysis of laminated and FGM composite annular sector plates. Compos. Part B Eng. 2019, 157, 182–194. [Google Scholar] [CrossRef]
  86. Nguyen, H.N.; Tan, T.C.; Luat, D.T.; Phan, V.D.; Thom, D.V.; Minh, P.V. Research on the buckling behavior of functionally graded plates with stiffeners based on the third–order shear deformation theory. Materials 2019, 12, 1262. [Google Scholar] [CrossRef] [PubMed]
  87. Lanc, D.; Vo, T.P.; Turkalj, G.; Lee, J. Buckling analysis of thin–walled functionally graded sandwich box beams. Thin Wall. Struct. 2015, 86, 148–156. [Google Scholar] [CrossRef]
  88. Lanc, D.; Turkalj, G.; Vo, T.; Brnic, J. Nonlinear buckling behaviours of thin-walled functionally graded open section beams. Compos. Struct. 2016, 152, 829–839. [Google Scholar] [CrossRef] [Green Version]
  89. Kim, J.; Zur, K.K.; Reddy, J.N. Bending, free vibration, and buckling of modified couples stress-based functionally graded porous micro-plates. Compos. Struct. 2019, 209, 879–888. [Google Scholar] [CrossRef]
  90. Barretta, R.; Feo, L.; Luciano, R.; Marotti de Sciarra, F.; Penna, R. Functionally graded Timoshenko nanobeams: A novel nonlocal gradient formulation. Compos. Part B Eng. 2016, 100, 208–219. [Google Scholar] [CrossRef]
  91. Apuzzo, A.; Barretta, R.; Faghidian, S.A.; Luciano, R.; Marotti de Sciarra, F. Nonlocal strain gradient exact solutions for functionally graded inflected nano-beams. Compos. Part B Eng. 2019, 164, 667–674. [Google Scholar] [CrossRef]
  92. Carrera, E. Co reissner-mindlin multilayered plate elements including Zig-Zag and interlaminar stress continuity. Int. J. Numer. Meth. Eng. 1996, 39, 1797–1820. [Google Scholar] [CrossRef]
  93. Carrera, E. Developments, ideas and evaluations based upon the Reissner’s mixed theorem in the modeling of multilayered plates and shells. Appl. Mech. Rev. 2001, 54, 301–329. [Google Scholar] [CrossRef]
  94. Carrera, E. Historical review of Zig-Zag theories for multilayered plates and shells. Appl. Mech. Rev. 2003, 56, 287–308. [Google Scholar] [CrossRef]
  95. Carrera, E. Theories and finite elements for layered plates and shells: A unified compact formulation with numerical assessment and benchmarking. Arch. Comput. Meth. Eng. 2003, 10, 215–296. [Google Scholar] [CrossRef]
  96. Carrera, E. On the use of the Murakami’s Zig-Zag function in the modeling of layered plates and shells. Comput. Struct. 2004, 82, 541–554. [Google Scholar] [CrossRef]
  97. Maturi, D.A.; Ferreira, A.J.M.; Zenkour, A.M.; Mashat, D.S. Analysis of laminated shells by murakami’s Zig–Zag theory and radial basis functions collocation. J. Appl. Math. 2013, 2013, 14. [Google Scholar] [CrossRef]
  98. Brischetto, S.; Carrera, E.; Demasi, L. Improved bending analysis of sandwich plates using a Zig–Zag function. Compos. Struct. 2009, 89, 408–415. [Google Scholar] [CrossRef]
  99. Hu, H.; Belouettar, S.; Daya, E.M.; Potier-Ferry, M. Evaluation of kinematic formulations for viscoelastically damped sandwich beam modeling. J. Sandw. Struct. Mater. 2006, 8, 477–495. [Google Scholar] [CrossRef]
  100. Murakami, H. Laminated composite plate theory with improved in-plane responses. J. Appl. Mech. 1986, 53, 661–666. [Google Scholar] [CrossRef]
  101. Lenci, S.; Tarantino, A.M. Chaotic dynamics of an elastic beam resting on a Winkler–type soil. Chaos Soliton. Fract. 1996, 7, 1601–1614. [Google Scholar] [CrossRef]
Figure 1. Geometric features of a laminated sandwich plate with an inner damaged soft-core.
Figure 1. Geometric features of a laminated sandwich plate with an inner damaged soft-core.
Materials 12 02444 g001
Figure 2. Through-the-thickness representation of the function f ( k ) employed to characterize the volume fraction distribution of the fibers in the k -th layer of the plate, for increasing values of the exponent α : (a) f 1 ( k ) ; (b) f 2 ( k ) .
Figure 2. Through-the-thickness representation of the function f ( k ) employed to characterize the volume fraction distribution of the fibers in the k -th layer of the plate, for increasing values of the exponent α : (a) f 1 ( k ) ; (b) f 2 ( k ) .
Materials 12 02444 g002
Figure 3. Discrete element: Node numbering and natural coordinates.
Figure 3. Discrete element: Node numbering and natural coordinates.
Materials 12 02444 g003
Figure 4. Discrete domain, identification of the elements and node numbering.
Figure 4. Discrete domain, identification of the elements and node numbering.
Materials 12 02444 g004
Figure 5. Through-the-thickness representation of the volume fraction distribution V F for several schemes: (a) Scheme 1; (b) Scheme 2; (c) Scheme 3; (d) Scheme 4.
Figure 5. Through-the-thickness representation of the volume fraction distribution V F for several schemes: (a) Scheme 1; (b) Scheme 2; (c) Scheme 3; (d) Scheme 4.
Materials 12 02444 g005
Figure 6. Variation of the first natural frequency f 1 [Hz] for sandwich plates with a damaged soft-core due to an increasing damage D , for three different lamination schemes: (a) (0°/core/0°); (b) (30° /core/45°); (c) (−45° /core/45°).
Figure 6. Variation of the first natural frequency f 1 [Hz] for sandwich plates with a damaged soft-core due to an increasing damage D , for three different lamination schemes: (a) (0°/core/0°); (b) (30° /core/45°); (c) (−45° /core/45°).
Materials 12 02444 g006
Figure 7. Variation of the first three natural frequencies f 1 , f 2 , f 3 [Hz] for sandwich plates with a damaged soft-core due to an increasing value of the exponent α of the through-the-thickness distribution of the fiber volume fraction, for two different schemes: (a) Scheme 2; (b) Scheme 3.
Figure 7. Variation of the first three natural frequencies f 1 , f 2 , f 3 [Hz] for sandwich plates with a damaged soft-core due to an increasing value of the exponent α of the through-the-thickness distribution of the fiber volume fraction, for two different schemes: (a) Scheme 2; (b) Scheme 3.
Materials 12 02444 g007
Figure 8. Variation of the first three natural frequencies f 1 , f 2 , f 3 [Hz] for a sandwich plates with a damaged soft-core for several lamination schemes depending on the angle parameter θ which defines the in-plane fiber orientation.
Figure 8. Variation of the first three natural frequencies f 1 , f 2 , f 3 [Hz] for a sandwich plates with a damaged soft-core for several lamination schemes depending on the angle parameter θ which defines the in-plane fiber orientation.
Materials 12 02444 g008
Figure 9. Variation of the first five natural frequencies f n [Hz], for n = 1 , , 5 , of sandwich plates with a damaged soft-core for varying mechanical properties of the constituents: (a) Variation of CNT mass fraction; (b) Variation of fiber mass fraction.
Figure 9. Variation of the first five natural frequencies f n [Hz], for n = 1 , , 5 , of sandwich plates with a damaged soft-core for varying mechanical properties of the constituents: (a) Variation of CNT mass fraction; (b) Variation of fiber mass fraction.
Materials 12 02444 g009
Figure 10. First five mode shapes of the mechanical configurations defined in Table 4: (a) Case 1; (b) Case 2; (c) Case 3; (d) Case 4; (e) Case 5; (f) Case 6; (g) Case 7.
Figure 10. First five mode shapes of the mechanical configurations defined in Table 4: (a) Case 1; (b) Case 2; (c) Case 3; (d) Case 4; (e) Case 5; (f) Case 6; (g) Case 7.
Materials 12 02444 g010
Table 1. Mechanical properties of a single CNT fiber.
Table 1. Mechanical properties of a single CNT fiber.
Hill’s Elastic ModuliDensity
k C = 271 GPa ρ C = 1400 kg / m 3
l C = 88 GPa
m C = 17 GPa
n C = 1089 GPa
p C = 442 GPa
Table 2. Mechanical properties of the layer constituents.
Table 2. Mechanical properties of the layer constituents.
ConstituentYoung’s ModuliShear ModuliPoisson’s RatiosDensity
Carbon fibers E 11 F = 230 GPa G 12 F = 50 GPa ν 12 F = 0.20 ρ F = 1800 kg / m 3
E 22 F = 15 GPa ν 23 F = 0.25
Epoxy resin E M = 3.27 GPa ν M = 0.38 ρ M = 1200 kg / m 3
Table 3. First ten natural frequencies of a sandwich plate with an inner soft-core for three different lamination schemes: comparison with the 3D-FE solution.
Table 3. First ten natural frequencies of a sandwich plate with an inner soft-core for three different lamination schemes: comparison with the 3D-FE solution.
Mode3D-FE
N d o f s = 29673
RMZ
N d o f s = 3087
RM
N d o f s = 2205
%diff
(RMZ)
%diff
(RM)
Lamination scheme: (0°/core/0°)
1142.126137.126169.1273.52%19.00%
2195.625189.516224.8193.12%14.92%
3289.303280.546324.7083.03%12.24%
4303.824288.429408.3905.07%34.42%
5346.073329.618453.3754.75%31.01%
6414.415401.529463.4763.11%11.84%
7421.211402.429533.7744.46%26.72%
8498.657469.706635.3515.81%27.41%
9527.513504.420650.8394.38%23.38%
10534.573505.028678.3795.53%26.90%
Lamination scheme: (30°/core/45°)
1129.776125.889147.1342.99%13.38%
2211.140203.739244.1333.51%15.63%
3265.283253.703327.3024.37%23.38%
4299.738288.099354.3093.88%18.21%
5371.903354.527462.2244.67%24.29%
6398.033381.561476.9604.14%19.83%
7426.714404.619562.8405.18%31.90%
8483.191459.287616.6014.95%27.61%
9506.505484.458616.7244.35%21.76%
10541.210512.228704.8095.36%30.23%
Lamination scheme: (−45°/core/45°)
1107.180105.297112.9301.76%5.36%
2206.687201.053225.2342.73%8.97%
3206.687201.516225.2342.50%8.97%
4290.914281.838322.7093.12%10.93%
5346.351334.314392.8463.48%13.42%
6348.741337.043395.2403.35%13.33%
7419.280403.125479.2723.85%14.31%
8419.280404.002479.2723.64%14.31%
9515.903494.071608.8114.23%18.01%
10515.903496.715608.8113.72%18.01%
Table 4. Comparison of the frequency parameter ϖ for a two phase fully clamped plate.
Table 4. Comparison of the frequency parameter ϖ for a two phase fully clamped plate.
ModeLei et al. [71]
(kp-Ritz)
Lei et al. [71]
(Commercial FE)
RMZ
(MIX)
RMZ
(HT)
Case 1
116.66716.70716.67117.383
222.13822.25322.09823.381
332.23732.37832.21133.630
432.42432.85732.24234.397
535.67435.80935.65237.437
637.36737.44737.40340.615
Case 2
118.04518.08318.05519.071
223.49823.60623.44825.152
333.91534.33833.71236.116
434.36134.46734.32436.449
537.36737.44737.40339.892
637.69337.78637.63740.616
Table 5. Mechanical configurations investigated to show the variation of the mode shapes due to the stacking sequence, the damage parameter and the exponent of the reinforcing fiber distribution.
Table 5. Mechanical configurations investigated to show the variation of the mode shapes due to the stacking sequence, the damage parameter and the exponent of the reinforcing fiber distribution.
CaseStacking
Sequence
Damage
D
Exponent
α
1(0°/core/0°)0.001
2(0°/core/0°)0.501
3(0°/core/0°)0.5012
4(30°/core/45°)0.001
5(30°/core/45°)0.501
6(30°/core/45°)0.504
7(30°/core/45°)0.5012

Share and Cite

MDPI and ACS Style

Bacciocchi, M.; Luciano, R.; Majorana, C.; Tarantino, A.M. Free Vibrations of Sandwich Plates with Damaged Soft-Core and Non-Uniform Mechanical Properties: Modeling and Finite Element Analysis. Materials 2019, 12, 2444. https://doi.org/10.3390/ma12152444

AMA Style

Bacciocchi M, Luciano R, Majorana C, Tarantino AM. Free Vibrations of Sandwich Plates with Damaged Soft-Core and Non-Uniform Mechanical Properties: Modeling and Finite Element Analysis. Materials. 2019; 12(15):2444. https://doi.org/10.3390/ma12152444

Chicago/Turabian Style

Bacciocchi, Michele, Raimondo Luciano, Carmelo Majorana, and Angelo Marcello Tarantino. 2019. "Free Vibrations of Sandwich Plates with Damaged Soft-Core and Non-Uniform Mechanical Properties: Modeling and Finite Element Analysis" Materials 12, no. 15: 2444. https://doi.org/10.3390/ma12152444

APA Style

Bacciocchi, M., Luciano, R., Majorana, C., & Tarantino, A. M. (2019). Free Vibrations of Sandwich Plates with Damaged Soft-Core and Non-Uniform Mechanical Properties: Modeling and Finite Element Analysis. Materials, 12(15), 2444. https://doi.org/10.3390/ma12152444

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop