Linear Static Behavior of Damaged Laminated Composite Plates and Shells

A mathematical scheme is proposed here to model a damaged mechanical configuration for laminated and sandwich structures. In particular, two kinds of functions defined in the reference domain of plates and shells are introduced to weaken their mechanical properties in terms of engineering constants: a two-dimensional Gaussian function and an ellipse shaped function. By varying the geometric parameters of these distributions, several damaged configurations are analyzed and investigated through a set of parametric studies. The effect of a progressive damage is studied in terms of displacement profiles and through-the-thickness variations of stress, strain, and displacement components. To this end, a posteriori recovery procedure based on the three-dimensional equilibrium equations for shell structures in orthogonal curvilinear coordinates is introduced. The theoretical framework for the two-dimensional shell model is based on a unified formulation able to study and compare several Higher-order Shear Deformation Theories (HSDTs), including Murakami’s function for the so-called zig-zag effect. Thus, various higher-order models are used and compared also to investigate the differences which can arise from the choice of the order of the kinematic expansion. Their ability to deal with several damaged configurations is analyzed as well. The paper can be placed also in the field of numerical analysis, since the solution to the static problem at issue is achieved by means of the Generalized Differential Quadrature (GDQ) method, whose accuracy and stability are proven by a set of convergence analyses and by the comparison with the results obtained through a commercial finite element software.


Introduction
Shell structures are becoming very popular due to their typical curved shapes that can characterize the structural behavior of these elements. For instance, the capability of transferring and holding external loads is substantially different from the one of a flat panel. Analogously, the free vibrations, stresses and strains, as well as the buckling loads, are highly affected by the shell curvature [1,2]. It should be pointed out that the geometry of a shell structure is completely defined when the corresponding middle surface is accurately described [1][2][3]. The characterization of these shapes is often linked to the greatest issues that could be encountered during the mechanical analysis of shells. These difficulties are even bigger if a doubly-curved surface must be studied, since each point of the domain is defined by two different radii of curvature. To the best of the authors' knowledge, this obstacle can be overcome if the doubly-curved surface at issue is described analytically by means of the differential geometry principles, as illustrated in the book by Kraus [3].
The mechanical behavior of these structures is even more intriguing to analyze if composite materials are introduced [1,4]. Indeed, the structural response of a composite shell is completely different from the corresponding structure made of conventional materials, such as isotropic mediums. The aim of these innovative materials is to enhance the mechanical behavior of the structure by combining two or more constituents, which could not provide the same response if taken separately. That is the case of a fiber-reinforced material. In this example, high-strength fibers are combined with an isotropic matrix. In general, the fibers (such as carbon or glass filaments) carry the external stresses, whereas the matrix keeps them together and acts as a shield to protect the reinforcing phase from the environmental factors, which could cause a deterioration of mechanical properties [4]. A stack of these fiber-reinforced mediums can be made to obtain the so-called laminated composites. In this circumstance, each layer (or ply) represents one of the various constituent of the composite. As a consequence, the structural response is affected by the number of layers, the orientation of the fibers, and the way they are combined. In general, each fiber-reinforced lamina is orthotropic, thus the corresponding constitutive elastic equations are required to describe properly the relation between stresses and strains. The mechanical analysis of laminated composite plates and shells is currently a recurring topic in the pertinent literature. In particular, the effect of the fiber orientation and the stacking sequence on the structural response has been hugely investigated in many papers to analyze the static [5][6][7][8][9][10][11][12][13][14][15][16][17] and dynamic  behavior of such structures.
In recent years, the idea of variable mechanical properties for a composite medium has begun to spread in order to optimize the structural performances and reach the best mechanical responses towards the environmental demands. Firstly, the well-known class of functionally graded materials (FGMs) should be mentioned [44][45][46][47][48][49][50][51][52][53][54][55][56][57][58][59][60], where a continuous variation of the material properties along the thickness of the structure is introduced to reduce the stress peaks at the layer interfaces of a laminate. Indeed, this gradual variation of the mechanical properties is applied to overcome some of the limitations of laminated composites, such as delamination and interlaminar cracks. Starting from the same ideas of FGMs, the class of functionally-graded Carbon Nanotube-reinforced composites is becoming very popular due to the introduction of nanoparticles as the reinforcing phase in the composite medium [61][62][63][64][65][66][67][68][69]. Even in this circumstance, the variation of the mechanical properties is obtained by defining the volume fraction distribution of the constituents by means of different through-the-thickness laws. It is clear that the mechanical properties of such composites turn out to be functions of the thickness coordinate. Alternatively, a different kind of composite material has been developed to obtain a variation of its properties within the domain defined by the middle surface of the structure. Such variability can be achieved by changing the orientation of the fiber in each point of the domain. This approach is known in the literature as variable angle tow (VAT) concept [70][71][72][73][74][75][76][77][78][79][80][81][82][83]. Similar variations could be obtained by varying the spatial distribution of the reinforcing fibers by increasing (or decreasing) the space among them [84], or adding a fiber-reinforced ply only in some areas of the domain [85]. This last example can be generalized by applying a smooth thickness profile to the composite structure at issue. As illustrated in [86][87][88][89][90][91][92][93][94], the thickness variation can affect the structural response through an optimal distribution of the materials.
The present paper also falls within the topic of variable stiffness structures, since it aims to investigate the static behavior of composite plates and shells with variable mechanical properties. Nevertheless, this variation is concentrated in some delimited areas of the composite to model a sudden deterioration of the material properties. The weakening at issue can occur when damage spreads within the medium that composes the structure. Consequently, the elastic properties of the structure are considerably reduced only in limited zones, whereas the undamaged areas of the structure present unchanged features. This configuration can be mathematically achieved by assigning particular smooth functions, such as the well-known two-dimensional Gaussian distribution or an ellipse shaped law, to the elastic parameters that define the mechanical properties of the structure. As highlighted in the paper by Ladevèze and Le Dantec [95], a generic damage can be described as a reduction of the material stiffness originating from microcracking and debonding, for instance. Different kinds of damage are also possible, as presented by Daudeville and Ladevèze [96], such as transverse matrix cracking and fiber ruptures. A complete review of the failures that can occur in a laminated composite medium is presented in the books by Reddy and Miravete [97] and by Murakami [98], in which it is specified also that the main issues in developing a mathematical model for damage are related to the various geometric scales involved in the failure progression. Independently from the scale of the damage, a failure is always associated with a global stiffness reduction of the structure. The stiffness changes at issue affect deflections, vibration characteristics, and the stress and strains distribution, as illustrated clearly in the work by Highsmith and Reifsnider [99]. In all the mentioned works, the reasons that cause damage are diverse. Since many variables are involved in the growth and progression of damage, this topic is often investigated from both the numerical [100][101][102][103][104][105][106][107] and experimental [108][109][110][111][112][113] points of view. Analogously, several mechanical approaches at a different scale were also proposed. Typically, the most exploited ones are variational methods [97], continuum damage models [97,[114][115][116][117], or approaches that take into account a plastic behavior of the structure [95,97,118,119]. The more appropriate approach should be chosen according to the failure to investigate. It should be mentioned that the causes that give rise to a damage have not been investigated in the present paper. Analogously, the point of the domain in which this damage arises is assumed a priori. Bearing in mind these hypotheses, several investigations are presented in this research to show the effect of the damage parameters (point of application, intensity, and width) on the linear static behavior of laminated shell structures.
Once the constitutive laws are specified, the governing equations are obtained in the framework of higher-order shear deformation theories (HSDTs), since these peculiar mechanical configurations could be ineffectively studied through the well-known first-order shear deformation theory (FSDT), as highlighted in many works during the last decades [61,[78][79][80][92][93][94][120][121][122][123][124][125][126]. Thus, an analytical model able to deal with several enriched displacement fields is employed to describe the mechanical behavior of these structures. The fundamental assumptions of this theoretical formulation can be found in [1,5,[18][19][20][127][128][129]. As far as the achievement of the solutions is concerned, a computational method is introduced. In particular, the strong form of the governing equation is solved numerically by means of the Generalized Differential Quadrature (GDQ) method [130,131]. The same numerical scheme is employed to evaluate both the geometric parameters of the doubly-curved surfaces used as reference domains and the through-the-thickness variations of stress, strain, and displacement components. For this purpose, a posteriori recovery procedure based on the three-dimensional equilibrium equations is proposed.

Shell Geometry
In this paper, differential geometry is used to evaluate those geometric parameters needed for the description of the shell middle surface, which are required in the fundamental operators of the governing equations [1,3]. Let us consider a generic shell element in a global reference system O x 1 x 2 x 3 as shown in Figure 1. Each point P within the three-dimensional medium of thickness h is identified by the vector R(α 1 , α 2 , ζ), which takes the following aspect R(α 1 , α 2 , ζ) = r(α 1 , α 2 ) + ζ n(α 1 The coordinate ζ specifies the normal direction along the shell thickness, whereas α 1 , α 2 are the orthogonal and principal curvilinear coordinates of the shell middle surface. By hypothesis, the coordinates α 1 , α 2 coincide with the parametric lines of curvature of the middle surface. It should be specified that their meaning is different according to the surface to describe [1,[92][93][94]. As shown in Figure 1, O α 1 α 2 ζ denotes the local reference system of the reference surface of the shell, which corresponds to its middle surface.
The normal direction  is identified by the following outward unit normal vector   Definitions (4) are based on differential geometry principles. It should be noted that the scalar parameters in (4) are clearly physical quantities whose magnitude depends on the coordinates  Figure 1. Doubly-curved shell element, edge identification and lamination scheme. The notation represents the stacking sequence of the composite, where θ (k) is the orientation of the k -th layer.
Each point P' of the shell middle surface is located by the position vector r(α 1 , α 2 ). Once this position vector is defined, its derivatives with respect to α 1 , α 2 are required. For the first two orders, one gets The normal direction ζ is identified by the following outward unit normal vector n(α 1 , α 2 ), described by the following vector product Once the position vector r(α 1 , α 2 ) is provided, by means of quantities in Equation (2), the well-known Lamè parameters A 1 (α 1 , α 2 ), A 2 (α 1 , α 2 ) are computed as scalar products Definitions (4) are based on differential geometry principles. It should be noted that the scalar parameters in (4) are clearly physical quantities whose magnitude depends on the coordinates α 1 , α 2 defined on the reference domain (middle surface). Further details concerning the Lamè parameters are illustrated in the book by Tornabene et al. [1].
For a doubly-curved surface, the principal radii of curvature R 1 (α 1 , α 2 ), R 2 (α 1 , α 2 ) are also required. The definitions below are valid if the coordinates α 1 , α 2 are orthogonal and principal In case of a singly-curved shell or a degenerate shell (which corresponds to a plate), one radius of curvature (or both of them for a flat structure) tends to infinite value. Finally, the geometric quantities H 1 (α 1 , α 2 , ζ), H 2 (α 1 , α 2 , ζ) are also needed. They are defined as follows and they are introduced to take into account the three-dimensional size effect related to the shell curvature. It should be stated that a finite domain is specified by setting limited values along each principal coordinate. In particular, the limitations at issue can be expressed analytically as follows in which α 0 i , α 1 i , for i = 1, 2, represent the minimum and the maximum boundary values along the coordinates α 1 , α 2 , respectively. When a laminated composite shell is considered, the overall thickness of the structure is given by where l is the total number of layers, whereas the index k stands for the geometric and mechanical parameters of the k-th ply. As it can be noted from Figure 1, h k represents the thickness of the k-th ply and it can be evaluated as follows in which ζ k and ζ k+1 are the lower and upper coordinates of the k-th layer, respectively.

Shell Structural Model
The displacement field of a generic laminated composite shell is defined by the three-dimensional displacements U 1 (α 1 , α 2 , ζ), U 2 (α 1 , α 2 , ζ), U 3 (α 1 , α 2 , ζ), which assume the following aspect (10) in which the order of kinematic expansion τ can be chosen arbitrarily. For conciseness purposes, the displacement components in Equation (10) can be collected into the vector U = U(α 1 , α 2 , ζ). According to the present formulation, which is able to deal with both the shear deformations and the stretching effect along the thickness of the structure, several HSDTs can be obtained by setting the maximum order of expansion N. The degrees of freedom of the model are given by the generalized displacements of the shell middle surface u (τ) for τ = 0, 1, 2, . . . , N, N + 1. It should be specified that the present theory belongs to the class of Equivalent Single Layer (ESL) approaches, since all the geometric and mechanical parameters are evaluated on the middle surface of the structure. It is also important to specify that only the orders τ = 0, 1 have a physical meaning. In particular, the translational displacement along α 1 , α 2 , ζ are obtained for τ = 0, whereas τ = 1 provides the corresponding rotational components. Further details about these aspects can be found in the book by Tornabene et al. [1]. The kinematic model in (10) is well-defined once the shear functions (or thickness functions) F τ (ζ) are specified for each order of kinematic expansion. A complete list of functions that can be chosen for this purpose is shown in [5]. In this paper, the power-law function ζ τ is used to define the displacement field up to the N-th order of expansion. As far as the (N + 1)-th expansion order is concerned, the corresponding thickness function F N+1 (ζ) coincides with the well-known Murakami's function Z = Z(ζ), whose analytical expression is given by where the dimensionless parameter z k (ζ) ∈ [−1, 1] is defined as This function could be required to capture the so-called zig-zag effect that could happen when peculiar lamination schemes are considered. Further details concerning the Murakami's function can be found in [92][93][94]129,130]. For the sake of completeness, it should be mentioned that this function allows the description of continuous three-dimensional displacements characterized by discontinuities in their derivatives at the interface between two adjacent layers. To sum up, the thickness functions are chosen as follows, for the corresponding orders of kinematic expansion It is clear that the maximum order of expansion N defines the structural model. Consequently, the acronyms ED N and EDZ N are used to specify the theory. As illustrated in the previous papers [93][94][95], the letter "E" specifies that the theory is based on an ESL approach, whereas the letter "D" states that the governing equations are deduced in terms of the generalized displacements. The letter "Z" is added to the notation only when the Murakami's function is embedded in the model.
The generalized strains of the middle surface ε (τ) (α 1 , α 2 ) are evaluated for each order τ as a function of the generalized displacements as follows ε (τ) = D Ω u (τ) (15) in which the kinematic operator D Ω is given by For the sake of completeness, the quantities collected in the generalized strain vector are shown below The complete treatise concerning the generalized strains and their definitions in extended notation can be found in the book by Tornabene et al. [1]. The meaning of higher-order terms is also explained. Quantities in (17) allow the evaluation also of the three-dimensional strains of the structure ε(α 1 , α 2 , ζ) where the strain vector is given by whereas the matrix Z (τ) (ζ) assumes the following form As far as the linear elastic constitutive relations are concerned, the stress components σ (k) (α 1 , α 2 , ζ) for the generic k-th orthotropic layer can be written as follows The elements of the constitutive matrix E (k) nm (α 1 , α 2 ) represent the material constants. They can be assumed equal to the reduced elastic coefficients if the hypothesis of plane stress is introduced Conversely, they are set equal to the non-reduced coefficients It should be noted that the normal stress σ n and the normal strain ε n are neglected if the plane stress hypothesis is introduced. Independently from this assumption, the elastic coefficients E (k) nm must be evaluated in the geometric local reference system O α 1 α 2 ζ. This aspect is extremely important when each orthotropic layer of the laminate has a different orientation. By means of the proper relations that allow the orientation in hand to be taken into account [1,4], the elastic coefficients can be evaluated as a function of the corresponding elastic coefficients related to the local reference system of the oriented ply C   23 ). The plane stress-reduced elastic coefficients written in the material reference system are defined as follows (27) in which the quantity ∆ (k) is given by With reference to definitions (26)- (28), all the remaining engineering constants can be deduced by means of the following relations In addition, the orientation of the material properties θ (k) must be specified for each layer for a complete mechanical characterization of the laminate. The notation θ (1) /θ (2) / . . . /θ (k) / . . . /θ (l) is used to specify the stacking sequence of the composite structure ( Figure 1). For completeness purposes, it should be recalled that an isotropic medium requires only two independent engineering constants and its properties are evaluated independently from the orientation of the material reference system. In the present paper, the aspect of damaged structures is introduced. A generic damage can be seen as a relatively concentrated deterioration of the mechanical properties of the elastic medium. Thus, peculiar functions can be introduced to model this rapid variation of the mechanical features of a layer. In particular, each engineering constant is affected by this sudden decay. Analytically speaking, all the engineering constants of the medium are multiplied by the factor Ψ (k) , which assumes the aspect below where δ (k) ∈ [0, 1] denotes the intensity of the damage. On the other hand, the functions ψ (k) E (α 1 , α 2 ) represent the two-dimensional normalized Gaussian function and an elliptic variation, respectively. The variation at issue are defined as follows ψ (k) It should be noted that such variations are applied at the point α On the other hand, the width of the damaged area is controlled by the size parameters Λ (k) 2 , which can be evaluated as a function of the corresponding quantities ∆ (k) ∈ [−1, 1] and β (k) ∈ [0, 2π] denote respectively the correlation parameter of the Gaussian distribution and the rotation of the ellipse, which is evaluated counterclockwise from the axis α 1 . For the sake of clarity, two examples of the application of the functions at issue are depicted in Figure 2 to understand the meaning of the various geometric quantities just introduced. The superscript k which characterizes each parameter means that the variations can assume different shapes through the various layers of the laminate. As a result, the engineering constants of each layer depend on the local coordinate α 1 , α 2 of the reference surface as can be observed from the following relations Once the generalized strain components (15) are evaluated, it is possible to compute also the stress resultants S (τ) (α 1 , α 2 ) for each order of kinematic expansion. For this purpose, the corresponding vector can be conveniently introduced for τ = 0, 1, 2, . . . , N, N + 1. The stress resultant in hand is defined by the following matrix notation where the stiffness matrix of the laminate A (τη) is given by for τ, η = 0, 1, 2, . . . , N, N + 1. As mentioned above, the definitions and the meaning of higher-order stress resultants are illustrated in the book by Tornabene et al. [1].   The constitutive operator A (τη) is a 9 × 9 matrix that assumes the following aspect In extended notation, relation (36) can be written as follows for a generic HSDT whose maximum order of expansion is given by N (embedded with the Murakami's function) (1) . . .
With reference to the well-known FSDT, the following correspondences can be deduced in which A, D, B are classically known as membrane stiffness matrix, bending stiffness matrix, and coupling stiffness matrix [1]. In fact, τ = 0 and τ = 1 are related respectively to the membrane stresses and bending moments (analogously for the strains), denoted by S (0) and S (1) (ε (0) and ε (1) as far as the generalized strains are concerned). On the other hand, for τ ≥ 2, relations among higher-order stresses and strains exist. Thus, higher-order terms could be coupled or uncoupled according to the mechanical properties of the structure (lamination scheme), following the same principles of the FSDT. For instance, it is well-known that the elements in B assume zero values for symmetric laminates. The same feature is repeated also for higher-order coupling terms. Each element inside the matrix (38) can be evaluated as follows for τ, η = 0, 1, 2, . . . , N, N + 1, n, m = 1, 2, 3, 4, 5, 6, and p, q = 0, 1, 2. It should be highlighted that the considerations just illustrated, in particular the correspondences in (40) nm (α 1 , α 2 ) in (41) stand for the elastic coefficients of the k-th layer. They can be computed as shown below for n, m = 4, 5 (42) in which the shear correction factor is identified by χ. Such distinction is essential to introduce the possibility to correct the distributions of the shear stresses. In fact, HSDTs in general do not require this coefficient. Nevertheless, in the present paper this correction is applied to the structural models up to the second order of kinematic expansion. Further details concerning the use of the shear correction factor in higher-order theories can be found in the works [93][94][95]. Some comments related to this correction are presented also in the next sections. Finally, it should be specified that integrals in (41) must be solved numerically. As illustrated in the previous paper by the authors [95], the Generalized Integral Quadrature (GIQ) method represents an accurate numerical tool to deal with this kind of issue. Classic integration schemes could be used as well [2].
By means of Hamilton's variational principle [1,4], a set of three equilibrium equations are obtained for each order of kinematic expansion τ. The structures are loaded only on their external surfaces by normal or shear pressures. These applied forces per unit surface are denoted by 3 . It should be noted that the superscript (+) specifies that the top surface is loaded, whereas the superscript (−) is used to identify a load applied on the other external surface. On the other hand, the subscript represents the coordinate along which the force is applied. The static equivalence principle is employed to evaluate the generalized load components acting on the shell middle surface, given by q 3 , which can be collected in the corresponding vector The following expressions can be used to define each quantity in (43) for each order of kinematic expansion represent the values that both the geometric quantities in (6) and the thickness functions assumes on the shell external surfaces (for ζ = ±h/2). At this point, the equilibrium equations are carried out for τ = 0, 1, 2, . . . , N, N + 1 where the operator D * Ω is given by By means of relations (36) and (15), the equilibrium equations can be written as a function of the generalized displacements of the middle surface where the fundamental operator L (τη) is given by the following 3 × 3 matrix for τ, η = 0, 1, 2, . . . , N, N + 1. A proper set of boundary conditions must be introduced to solve this system of partial differential equations. For the sake of conciseness, the boundary conditions are summarized in Table 1 for the external edges of the shell element depicted in Figure 1. In the following sections, the boundary conditions are specified by following the order "WSEN", where the four letters stand for the lateral edges of the shell ( Figure 1). For instance, "SSSS" and "CCCC" mean that all the four edges are simply-supported and clamped, respectively. On the other hand, "FCFC" is used to define a structure with the western and eastern edges free, whereas the others are clamped. A shell with only the northern edge free is specified as "CCCF". Finally, due to the considerable number of variables and definitions, a nomenclature section is included in Appendix A for the sake of clarity.

Numerical Solution
The GDQ technique is used to obtain and solve the discrete form of the governing equations. Indeed, the present numerical method provides the solution of the strong formulation of the fundamental Equation (47) by approximating the partial derivatives. For completeness purposes, the main aspects of this technique are presented briefly [131]. Let us consider a two-dimensional domain in which I N , I M denote respectively the total number of discrete points along x and y, where x, y represents the principal coordinates of the domain itself. The n-th order derivative with respect to x, the m-th order derivative with respect to y, and the n + m order mixed derivative of a generic function f (x, y) evaluated at the discrete point x i , y j are given by where the notation f x i , y j = f (ij) is employed. The weighting coefficients are specified by ς (n) x(ik) and ς (m) y(jl) . They can be evaluated by using the Lagrange polynomials. For conciseness purpose, only the definition of the coefficients related to the first coordinate x is given here. For the first-order derivatives, one gets whereas the weighting coefficients for higher-order derivatives can be computed recursively as assuming i, k = 1, 2, . . . , I N . The Lagrange polynomials are denoted by L(x) and they are defined as follows from which the corresponding first-order derivatives L (1) (x i ) can be deducted and evaluated at the point x i . The weighting coefficients ς (m) y(jl) can be computed following the me procedure. The weighting coefficients ς (n) x(ik) and ς (m) y(jl) can be conveniently collected in the corresponding matrices ς (n) x and ς (n) y , whose size is given by I N × I N and I M × I M , respectively. A simplified approach for the implementation of this technique can be used by collecting the grid points of the domain as specified by the blue arrow in Figure 3. Thus, the order of the point is given by   , f x y assume in each point are also collected by following the same order, so that the following algebraic vector f is obtained  I  I   I  I  I   T   I  I  I  I I  I  I The size of f is given by The weighting coefficients can be computed as where I is the identity matrix, whereas the symbol  stands for the Kronecker product. The size of each operator in (57)-(59) is specified under the various matrices.
are the weighting coefficients matrices that allow the corresponding derivatives in each point of the domain to be computed as simple matrix products Thus, the order of the point is given by ( The values that the generic function f (x, y) assume in each point are also collected by following the same order, so that the following algebraic vector f is obtained where f k = f x i , y j k , for i = 1, 2, . . . , I N , j = 1, 2, . . . , I M and k = i + (j − 1)I N . The size of f is given by (I N · I M ) × 1. The weighting coefficients can be computed as where I is the identity matrix, whereas the symbol ⊗ stands for the Kronecker product. The size of each operator in (57)-(59) is specified under the various matrices. C (n) x , C (m) y , and C (n+m) xy are the weighting coefficients matrices that allow the corresponding derivatives in each point of the domain to be computed as simple matrix products where f y , and f (n+m) xy represent the vectors of size (I N · I M ) × 1 which collect the derivatives of the function f (x, y) within the whole domain, following the same order employed in (56). It should be highlighted that this approach does not set any restriction on the choice of the grid distribution. The points can be placed within the domain following different schemes. In the present paper, the discrete points of the shell reference surface are given by where the meaning of r i and r j , for i = 1, 2, . . . , I N and j = 1, 2, . . . , I M , is specified in Table 2 for various distributions. Table 2. Different kinds of discrete grid distributions for the Generalized Differential Quadrature (GDQ) method [130]. Symbols k and N stand for i and I N respectively, along the α 1 direction. On the other hand, k and N specify j and I M respectively, along the α 2 direction. Notations L N , L N+1 are introduced to denote the corresponding Legendre polynomials.

Solution of the Static Problem
The static solution of the problem is obtained applying the GDQ technique, which allows the following discrete form of the governing equations to be written where K is the stiffness matrix of size (3 ), δ the displacement vector, and f the external load vector. Both δ and f are given by a vector characterized by the same size (3 × (N + 2) × (I N × I M )) × 1. It should be noted that the degrees of freedom of the problem are ordered following the scheme shown in Figure 3 for each order of kinematic expansion. By means of the so-called static condensation procedure, the size of the problem can be reduced by separating the degrees of freedom linked to the inner points of the domain (d) from the ones related to the boundary nodes (b). As a consequence, all the quantities in (65) are evaluated to consider this classification as follows The generalized displacements of the boundary points δ b are easily obtained as whereas the generalized displacements collected in δ d can be deduced from the following relation It can be observed that the size of the problem turns out to be 3 × (N + 2) × ((I N − 2) × (I M − 2)), which represents a reduced value if compared to the original size of the problem shown in (65). Finally, it should be specified that the derivatives needed to compute the stress resultants (36) involved in the boundary conditions for free or simply-supported edges (Table 1) are computed through the GDQ method.

Strain and Stress Recovery Procedure
Once the static solution is obtained, the three-dimensional elasticity equations in terms of stresses can be used to evaluate the through-the-thickness variation of all these quantities. In fact, the three indefinite equilibrium equations of elasticity written in an orthogonal curvilinear coordinate system allows the three-dimensional profiles of stresses and strains to be obtained, consequently, even if the higher-order structural model illustrated above is two-dimensional. The following treatise presents the main aspects of the recovery procedure at issue. Once the generalized displacements of the middle surface are computed, the kinematic model (10) is employed to obtain the three-displacements of the solid medium in each point α 1i , α 2 j of the reference surface for i = 1, 2, . . . , I N , j = 1, 2, . . . , I M , and m = 1, 2, . . . , I T , where I T denotes the number of discrete points along the thickness for each ply of the laminate. The Cheb-Gau-Lob distribution with I T = 25 is used to obtain a discrete distribution of points ζ m along the thickness of each layer for m = 1, 2, . . . , I T . It should be specified that F τ(m) represents the value that the corresponding thickness function assumes at the point ζ m .
By means of the GDQ technique, the generalized strains (15) can be computed in each grid point. As a consequence, the three-dimensional strain components related to each discrete point of the middle surface can be evaluated by using expression (18). The strains at issue assume the following aspect In particular, quantities ε 1(ijm) , ε 2(ijm) , γ 12(ijm) collected in (71) are required to obtain the membrane stresses σ 1(ijm) , σ 2(ijm) , and τ 12(ijm) in each point of the domain, by means of the constitutive relations (21). The three-dimensional equilibrium equations along the principal curvilinear coordinates are presented below The GDQ method is employed to approximate the partial derivatives of the membrane stresses as follows where ς represents the weighting coefficients for the derivatives within the shell middle surface. The discrete form of the equilibrium Equations (72)-(74) is obtained numerically through the GDQ technique as well in which ς ζ (1) mk are the weighting coefficients for the derivatives along ζ. The Lamè parameters A 1(ij) , A 2(ij) and the radii of curvature R 1(ij) , R 2(ij) are also computed in each discrete point of the reference surface. Relations (76) and (77), which must be written for each layer of the laminate, allow the shear stresses τ 1n , τ 2n to be obtained by imposing the boundary conditions on the bottom surface of the structure τ 1n(ij1) = q where q Conditions (80) are imposed through the following expressions, which represent the effective through-the-thickness variations of the shear stresses in each point of the domain for m = 1, 2, . . . , I T . It is evident that the symbol ζ m used in (81) and (82) takes into account the discrete grid distribution applied along each ply, since the correction of the stress components at issue affects the whole laminate thickness. Relations (79) must be seen as interlaminar compatibility conditions when the shear stresses of inner layers are computed. The GDQ method can be applied also to evaluate the derivatives of (81) and (82) with respect to α 1 , α 2 , respectively Equation (78) can be now employed to obtain the normal stress by applying the corresponding boundary condition at the bottom surface Analogously, another boundary condition must be introduced at the external surface of the top layer where q 3(ij) stand for the applied normal loads on the external surfaces of the shell. Finally, the effective through-the-thickness profile of the normal stress at issue is achieved for m = 1, 2, . . . , I T . The same considerations introduced for the previous profiles (81) and (82) are valid also in this circumstance. It should be mentioned that the stresses τ 1n , τ 2n , σ n just computed are continuous at the interfaces among the various layers of the structure. For any point of the three-dimensional solid, the shear strains γ 1n , γ 2n can be evaluated through the constitutive laws (21) as follows where C (m) pq represent the values that the elastic coefficients assume at the m-th discrete point of the layer. The normal strain ε n is computed by means of the constitutive relations (21) as well It should be noted that the quantity in (89) is evaluated neglecting the hypothesis of plane stress, thus the non-reduced elastic coefficients are employed. Moreover, the strain components (87)- (89) could be discontinuous at the layer interfaces since no compatibility interlaminar conditions are enforced.
Finally, the actual through-the-thickness variations of the membrane stresses σ 1 , σ 2 , τ 12 is computed through the constitutive laws (21) by using the effective value of the normal strain ε n . Consequently, the membrane stresses are computed without taking into account the hypothesis of plane stress, even for those structural theories that require that assumption.

Evaluation of the Elastic Coefficients
As mentioned in the previous sections, the integrals which appear in the definitions of the elastic nm (pq) shown in (41) are computed numerically by means of the GIQ method, whose key aspects are summarized here for a one-dimensional domain in which the only variable x is defined in the interval [a x , b x ]. Let us assume that the reference domain is discretized by N P points. One of the grid distribution of Table 2 could be equally chosen. The integral of a generic function f (x) within the interval x i , x j , with x i , x j ∈ [a x , b x ] is given by The procedure to compute the weighting coefficients for the integration w ij k are briefly presented in this paragraph. The following weighting coefficients are required for i, j = 1, 2, . . . , N P , where ς (1) x(ij) represents the corresponding weighting coefficients for the first-order derivatives that can be easily deducted from Equations (53) and (54). On the other hand, c denotes an arbitrary constant which can be assumed equal to c = b x + 10 −10 for accuracy purposes [131]. The weighting coefficients can be conveniently collected in the corresponding matrix ς (1) x , whose size is given by N P × N P . Its inverse must be computed as A generic element of the matrix W is denoted by w ij for i, j = 1, 2, . . . , N P . These quantities are used to evaluate the weighting coefficients w ij k for the integral (90) as shown below w ij k = w jk − w ik (93) assuming k = 1, 2, . . . , N P , where the indices i, j are related to the boundary points of the interval chosen for the integration. In this paper, the Cheb-Gau-Lob grid distribution is employed in each layer to obtain the quantities in in (41) as presented in (70), with m = 1, 2, . . . , N P . The total number of points is set equal to N P = 51 for each ply.

Applications
The present GDQ based approach is implemented in MATLAB code developed by the authors [132], which is employed to achieve the linear static solution of several damaged laminated plates and shells. All the numerical results are shown in this section, which is organized as follows: firstly, a set of convergence analyses in terms of displacements is carried out to prove the stability features of the numerical method for several lamination schemes and HSDTs; then, the recovery procedure is used to obtain the through-the-thickness variations of strain, stress, and displacement components of different structures. For the sake of conciseness, the structural elements are presented in Figure 4, where the employed discrete grid distribution is also specified. In particular, a square plate (Figure 4a), a singly-curved cylindrical panel with parabolic profile (Figure 4b), and a doubly-curved panel of revolution with parabolic meridian (Figure 4c), are considered.
All the mechanical and geometric characteristics of these structures are summarized in Table 3 for brevity purposes. On the other hand, Table 3 provides position vectors, boundary conditions, applied loads, lamination schemes, as well as the points of the domain for the evaluation of the through-the-thickness quantities, for each structural element. On the other hand, the mechanical properties of each constituent are listed in terms of engineering constants in Table 4. Finally, all the functions used for the damage modeling are specified in Table 5, together with all the parameters needed to represent correctly the Gaussian (or elliptic) distributions. Several HSDTs are employed in the following examples. In particular, it should be remarked that the theories up to the second order without the Murakami's function (FSDT χ=1.2 RS and ED2 χ=1.2 ) are always taken with the shear correction factor χ = 1.2, whereas the corresponding zig-zag models can be used also without correction for peculiar lamination schemes, as specified in [12,81,[93][94][95]. On the other hand, the first-order theories FSDT All the mechanical and geometric characteristics of these structures are summarized in Table 3 for brevity purposes. On the other hand, Table 3 provides position vectors, boundary conditions, applied loads, lamination schemes, as well as the points of the domain for the evaluation of the through-the-thickness quantities, for each structural element. On the other hand, the mechanical properties of each constituent are listed in terms of engineering constants in Table 4. Finally, all the functions used for the damage modeling are specified in Table 5, together with all the parameters needed to represent correctly the Gaussian (or elliptic) distributions. Several HSDTs are employed in the following examples. In particular, it should be remarked that the theories up to the second order without the Murakami's function ( , whereas the corresponding zig-zag models can be used also without correction for peculiar lamination schemes, as specified in [12,81,[93][94][95]. On the other hand, the first-order theories        Young's moduli: E 1 = 137.9 GPa, E 2 = E 3 = 8.96 GPa Shear moduli: G 12 = G 13 = 7.1 GPa, G 23 = 6.21 GPa Poisson's ratios: ν 12 = ν 13 = 0.3, ν 23 = 0.49 Table 5. Analytical expressions of the functions used to model the damage.  12 = 0, for k = 1, 2, 3 Case (a) − δ (n) = 0.00 for n = 1, 2, 3 Case (b) − δ (n) = 0.99, ∆ (n) = 5 for n = 1 Case (c) − δ (n) = 0.99, ∆ (1) = 10, ∆ (2) = 5, for n = 1, 2 Case (d) − δ (n) = 0.99, ∆ (1) = 15, ∆ (2) = 10, ∆ (3) = 5 for n = 1, 2, 3 (f) Function 6 (Ell-1) is taken into account for both these two cases, whereas the corresponding theory embedded with the Murakami's function (FSDTZ χ=1.2 RS ) is also used for the laminated composite, whose stacking sequence is given by (−45/45). Table 5 is employed to define the variable mechanical properties of each layer. It should be noted that the convergence analyses are carried out for an extremely concentrated damage characterized by a high deterioration of the mechanical properties. All the grid distributions shown in Table 2 are considered assuming I N = I M and increasing the number of points up to 51. The convergence graphs are presented in Figure 5 in terms of central deflection u (0) 3 0.5L x , 0.5L y . Analogously, the numerical values of the same quantity are shown in Table 6 for all the grid distributions.

The function Gau-1 shown in
From the various plots depicted in Figure 5, it can be noted that both the Legendre-Gauss (Leg-Gau) and Chebyshev-Gauss (Cheb-Gau) grid distributions provide a convergent behavior when the static analysis of damaged structures has to be performed. Nevertheless, the Leg-Gau distribution turned out to be the best discrete grid also when the GDQ method is employed to obtain the static response of plates and shells subjected to point or line loads [12]. Since both damages and concentrated forces represent mechanical discontinuities which this approach is able to deal with, the Leg-Gau distribution is used also in this paper.
Finally, it should be highlighted that a good convergence is reached by using many sampling points. Nevertheless, the test at issue is characterized by a noticeable mechanical discontinuity. For this reason, the value I N = I M = 51 is set in all the following examples independently from the damage features.

Isotropic Square Plate
The same flat structure of the convergence analysis is also used in this application, assuming completely clamped boundary conditions (CCCC). Two damage configurations are considered: Gau-2 and Gau-3. The former provides the increase of damage intensity δ, whereas the latter deals with the extension increase (∆ 1 = ∆ 2 = ∆). The damage is applied in the center of the plate in both these cases (α 1m = α 2m = 0.5 m). The vertical displacement profiles evaluated along the x direction for y = 0.5L y are shown in Figure 6a,b, respectively.
As expected, the displacements increase when the damage grows or when it involves a bigger area of the domain. The results are compared with the solution obtained with a finite element commercial software (Strand 7) for the undamaged case, by using a three-dimensional model (3D-FEM). The profiles in hand are shown in the graphs through a dotted notation and match perfectly with the corresponding GDQ curves. The through-the-thickness variations of strain, stress, and displacement components are presented in Figures 7-9 for the Gau-2 damage model, whereas Figures 10-12 deal with the results related to the Gau-3 damage. Even in these circumstances, the 3D-FEM for the undamaged plate is depicted by black dots and it is in good agreement with the present solutions. It can be noted that the damage features can modify the structural response. In particular, displacements are mostly affected by the damage properties. It should be specified that all the HSDTs considered in this application do not require Murakami's function since the structure is made of a sole isotropic layer. By comparing also the strain profiles related to these cases (shown respectively in Figures 7 and 10), it is possible to observe that a greater extension of damage (keeping the magnitude constant) has more influence on the strain patter than an extremely high decay of mechanical properties although concentrated in a smaller area of the domain. These two circumstances are modeled respectively through Gau-3 and Gau-2 functions.

Isotropic Square Plate
The same flat structure of the convergence analysis is also used in this application, assuming completely clamped boundary conditions (CCCC). Two damage configurations are considered: Gau-2 and Gau-3. The former provides the increase of damage intensity  , whereas the latter deals with the extension increase (  As expected, the displacements increase when the damage grows or when it involves a bigger area of the domain. The results are compared with the solution obtained with a finite element commercial software (Strand 7) for the undamaged case, by using a three-dimensional model (3D-FEM). The profiles in hand are shown in the graphs through a dotted notation and match perfectly with the corresponding GDQ curves. The through-the-thickness variations of strain, stress, and displacement components are presented in Figures 7-9 for the Gau-2 damage model, whereas Even in these circumstances, the 3D-FEM for the undamaged plate is depicted by black dots and it is in good agreement with the present solutions. It can be noted that the damage features can modify the structural response. In particular, displacements are mostly affected by the damage properties. It should be specified that all the HSDTs considered in this application do not require Murakami's function since the structure is made of a sole isotropic layer. By comparing also the strain profiles related to these cases (shown respectively in Figures 7 and 10), it is possible to observe that a greater extension of damage (keeping the magnitude constant) has more influence on the strain patter than an extremely high decay of mechanical properties although concentrated in a smaller area of the domain. These two circumstances are modeled respectively through Gau-3 and Gau-2 functions.

Laminated Square Plate
A completely clamped (CCCC) laminated plate is considered in this paragraph. The geometric features are the same as of the previous test, but the lamination scheme is given by   which each layer is made of the same material (see Table 3). Two damage types are considered and they are mathematically described by functions Gau-4 and Gau-5 of Table 4. As in the previous application, the damage is located in the central node of the structure (

Laminated Square Plate
A completely clamped (CCCC) laminated plate is considered in this paragraph. The geometric features are the same as of the previous test, but the lamination scheme is given by (−45/0/45), in which each layer is made of the same material (see Table 3). Two damage types are considered and they are mathematically described by functions Gau-4 and Gau-5 of Table 4. As in the previous application, the damage is located in the central node of the structure (α 1m = α 2m = 0.5 m). Both the intensity and the size of the deterioration of the mechanical properties are increased in the first case (Gau-4), involving each lamina in the same manner. On the other hand, the damage described by Gau-5 affects the three layers progressively, starting from the bottom ply. In addition, when the damage expands across the upper layers the damaged areas in the lower ones are increased. The graphs of the vertical displacement evaluated along coordinate x with y = 0.5L y are depicted in Figure 13a,b for the two damage models, respectively. The 3D-FEM profiles are also shown for comparison purposes. On the other hand, the through-the-thickness variations of strains, stresses, and displacements related to the function Gau-4 are presented in Figures 14-16. Figures 17-19 show the same quantities for the function Gau-5. The same considerations of the previous example are still valid in this circumstance. It should be noted that similar profiles are obtained by means of the various HSDTs. For completeness purposes, the third-order shear deformation theory is taken with and without the Murakami's function. No significant differences are observable between the ED3 and the EDZ3.
comparison purposes. On the other hand, the through-the-thickness variations of strains, stresses, and displacements related to the function Gau-4 are presented in Figures 14-16. Figures 17-19 show the same quantities for the function Gau-5. The same considerations of the previous example are still valid in this circumstance. It should be noted that similar profiles are obtained by means of the various HSDTs. For completeness purposes, the third-order shear deformation theory is taken with and without the Murakami's function. No significant differences are observable between the ED3 and the EDZ3 .

Cylindrical Panel
A FCFC sandwich structure with a central soft-core is investigated in this paragraph. The shell geometry is given by the singly-curved cylindrical surface with parabolic profile depicted in Figure 4b. The geometric parameters required to represent the surface at issue are listed in Table 3. For the sake of clarity, Figure 4d is added to explain the meaning of the parameters related to the parabolic arch which represents the profile of this cylinder. The same figure can be used as a reference to describe the parabolic meridian of the shell investigated in the following paragraph. A sandwich lamination scheme is considered here. In particular, the mechanical properties of the isotropic soft-core are considerably lower in comparison with the ones of the external orthotropic sheets (Tables 3 and 4). The damage is modeled through the elliptic function defined as Ell-1 in Table 5.

Cylindrical Panel
A FCFC sandwich structure with a central soft-core is investigated in this paragraph. The shell geometry is given by the singly-curved cylindrical surface with parabolic profile depicted in Figure 4b. The geometric parameters required to represent the surface at issue are listed in Table 3. For the sake of clarity, Figure 4d is added to explain the meaning of the parameters related to the parabolic arch which represents the profile of this cylinder. The same figure can be used as a reference to describe the parabolic meridian of the shell investigated in the following paragraph. A sandwich lamination scheme is considered here. In particular, the mechanical properties of the isotropic soft-core are considerably lower in comparison with the ones of the external orthotropic sheets (Tables 3 and 4). The damage is modeled through the elliptic function defined as Ell-1 in Table 5. The deterioration of the mechanical properties involves only the external layers and the elliptic shape describes a damage that affects the shell along all its length along the coordinate y . The parametric study aims to investigate the effect of a damage that gradually expands itself along the first coordinate  , as can be noticed from the parameters of the elliptic function Ell-1 in Table 5 ( On the other hand, the intensity of the damage is kept constant. As highlighted in many works [93][94][95], Murakami's function is needed to deal with this kind of

Cylindrical Panel
A FCFC sandwich structure with a central soft-core is investigated in this paragraph. The shell geometry is given by the singly-curved cylindrical surface with parabolic profile depicted in Figure 4b. The geometric parameters required to represent the surface at issue are listed in Table 3. For the sake of clarity, Figure 4d is added to explain the meaning of the parameters related to the parabolic arch which represents the profile of this cylinder. The same figure can be used as a reference to describe the parabolic meridian of the shell investigated in the following paragraph. A sandwich lamination scheme is considered here. In particular, the mechanical properties of the isotropic soft-core are considerably lower in comparison with the ones of the external orthotropic sheets (Tables 3 and 4). The damage is modeled through the elliptic function defined as Ell-1 in Table 5.
The deterioration of the mechanical properties involves only the external layers and the elliptic shape describes a damage that affects the shell along all its length along the coordinate y. The parametric study aims to investigate the effect of a damage that gradually expands itself along the first coordinate ϕ, as can be noticed from the parameters of the elliptic function Ell-1 in Table 5 . On the other hand, the intensity of the damage is kept constant. As highlighted in many works [93][94][95], Murakami's function is needed to deal with this kind of mechanical configuration. Thus, only the corresponding structural theories are considered up to the third-order of expansion (FSDTZ χ=1 RS , EDZ2 χ=1 , and EDZ3). The first-order and second-order models are taken without the shear correction factor χ = 1, as explained in [93][94][95]. The results of the displacement analysis are shown in Figure 20. In particular, both the tangential component u   (Figure 20b) increase when the damage affects wider areas, as expected. Similar graphs are obtained by the various HSDTs. The through-the-thickness variations of strains, stresses, and three-dimensional displacements are depicted respectively in Figures 21-23, where the so-called zig-zag effect is clearly observable. In general, comparable profiles are obtained through all the structural models, even though the FSDTZ χ=1 RS provides curves that are slightly detached, especially when the plotted quantities present a lower order of magnitude. The deterioration of the mechanical properties involves only the external layers and the elliptic shape describes a damage that affects the shell along all its length along the coordinate y . The parametric study aims to investigate the effect of a damage that gradually expands itself along the first coordinate  , as can be noticed from the parameters of the elliptic function Ell-1 in EDZ2   , and EDZ3 ). The first-order and second-order models are taken without the shear correction factor 1   , as explained in [93][94][95]. The results of the displacement analysis are shown in Figure 20. In particular, both the tangential component  

Doubly-Curved Panel of Revolution
A branch of a parabolic arch is employed to describe the meridian curve of a CCCF doubly-curved panel of revolution ( Figure 4c). As illustrated in the previous application, the meaning of the geometric parameters of the parabola listed in Table 3 is shown in Figure 4d for the sake of clarity. The damage is modeled using the mathematical expression denoted as Gau-6. In this circumstance, the origin of the Gaussian variation is located at   , where n represents the layers progressively involved by the damage (the count goes from the bottom to the upper layer of the laminated structure). Both the intensity and the size of the damage are kept constant in each ply. The displacement analysis is presented in Figure 24 for the tangential   The effect of the damage parameters on the strain, stress, and displacement profile along the shell thickness is investigated in Figures 25-27. By comparing all the curves in these figures, it is possible to note that different through-the-thickness profiles are obtained varying the order of kinematic expansion of the structural theory. Such differences are more evident in the normal strains and stresses ( n  and n  , respectively). Finally, it should be noted that the boundary conditions at the external surface are always well-enforced for all the static analyses presented in this section.

Doubly-Curved Panel of Revolution
A branch of a parabolic arch is employed to describe the meridian curve of a CCCF doubly-curved panel of revolution ( Figure 4c). As illustrated in the previous application, the meaning of the geometric parameters of the parabola listed in Table 3 is shown in Figure 4d for the sake of clarity. The damage is modeled using the mathematical expression denoted as Gau-6. In this circumstance, the origin of the Gaussian variation is located at α (k) 1m = 0.75(ϕ 1 − ϕ 0 ) + ϕ 0 , α (k) 2m = 0.50(ϑ 1 − ϑ 0 ) + ϑ 0 , for k = 1, 2, . . . , n, where n represents the layers progressively involved by the damage (the count goes from the bottom to the upper layer of the laminated structure). Both the intensity and the size of the damage are kept constant in each ply. The displacement analysis is presented in Figure 24 for the tangential u (0) 2 and vertical u (0) 3 components, from which it is possible to observe that the plotted quantities, measured along the first coordinate ϕ for ϑ = 0, increase by expanding the damage to the upper layers.
The effect of the damage parameters on the strain, stress, and displacement profile along the shell thickness is investigated in Figures 25-27. By comparing all the curves in these figures, it is possible to note that different through-the-thickness profiles are obtained varying the order of kinematic expansion of the structural theory. Such differences are more evident in the normal strains and stresses (ε n and σ n , respectively). Finally, it should be noted that the boundary conditions at the external surface are always well-enforced for all the static analyses presented in this section.

Conclusions and Remarks
A numerical analysis was proposed in the paper to deal with the linear static behavior of damaged plates and shells. The solutions were obtained through the application of the GDQ method, which was also employed to evaluate the through-the-thickness variations of stress and strain components by means of a recovery procedure based on the three-dimensional equilibrium equations for a shell structure. Several boundary conditions, geometric shapes, stacking sequences, and load combinations were considered and investigated through a variety of higher-order structural models, including the Murakami's function for the zig-zag effect when needed. The mechanical properties of the shells were taken as variable to model a damaged configuration. In particular, the two-dimensional Gaussian function and an elliptic expression were introduced to model a damage which affects all the engineering constants of the corresponding medium in concentrated areas of the shell domain. Several parametric investigations were carried out to study the effect of the various parameters that describe those functions. The following observations are noted:

•
The results were presented in terms of displacement profiles related to the shell middle surface and through-the-thickness variations of strain, stress, and displacement components at specific points of the domain. It was proven that the damage affects all these static quantities.

•
The numerical approach shows a convergent behavior when applied to this kind of structural problems. • It was proven that the damage modifies the structural response. In general, the displacement increases when the damage intensity grows or when it involves a bigger area of the domain; analogously, the overall displacement is greater when the damage spreads to the adjacent layers. • As far as stresses and strains are concerned, the effects of the damage on their through-the-thickness profiles are lower than the effects caused on the corresponding displacements; nevertheless, it is possible to observe noticeable differences between the profiles which are related respectively to the minimum and maximum values of the parameters used in the parametric analyses to model the deterioration of the mechanical properties. • Similar profiles in terms of strains, stresses, and displacements are obtained by means of different HSDTs. Consequently, all these enriched kinematic models are able to deal with damaged composite plates and shells. • No significant differences are observable between a higher-order model and the corresponding one embedded with the Murakami's function when a laminated composite is analyzed. On the other hand, the Murakami's function is required to deal with sandwich structures with an inner soft-core.
As a future development, the authors believe that the present damage model should be also employed in the dynamic field to investigate how the deterioration of the mechanical properties of a generic medium affects the corresponding structural response. In addition, the authors will attempt to compare the present approach with some experimental and numerical tests available in the literature. notations used in the paper are presented. In particular, the nomenclature employed in the previous sections is schematically presented in Table A1.