Effect of Curvilinear Reinforcing Fibers on the Linear Static Behavior of Soft-Core Sandwich Structures

The present research deals with the linear static behavior of soft-core sandwich plates and shells. The external skins are reinforced by curvilinear fibers. Their curved paths are described by a general mathematical law that allows the definition of arbitrary placements. The mechanical behavior of these structures is modeled through several Higher-order Shear Deformation Theories (HSDTs) including the zig-zag effect, based on an Equivalent Single Layer (ESL) approach. The solution of the governing equations is achieved numerically by means of the Generalized Differential Quadrature (GDQ) method. A huge number of parametric investigations is proposed in graphical and tabular forms to highlight the influence of the fiber orientation on the static response. The results prove that the structural behavior is affected by such parameters. Thus, the desired structural behavior can be modified by means of a proper choice of the fiber orientation.

A brief literature review is presented below to introduce this topic.For this purpose, the authors would like to apologize for the possible omissions of those contributions that could be considered as important steps in the development of this subject.A more complete literature review can be found in the PhD Thesis by Groh [57].A pioneering application of curvilinear fibers in composite structures was presented by Hyer and Charette [58] to deal with a rectangular plate with centrally located circular hole.Hyer et al. [59] continued the same topic by studying also the manufacturing problem.A plane elasticity problem of a symmetrically laminated composite with curvilinear fibers was investigated by Gürdal and Olmedo [60] to evaluate the stiffness variation effect on the elastic response of the structure.In their work, the orientation of the fibers was defined by two parameters, which consisted of the fiber angle at the center of the laminate and the fiber angle at a specified distance from the center.The manufacturing process of the proposed approach was progressively studied by Waldarth et al. [61], Setoodeh and Gürdal [62], and Setoodeh et al. [63,64].Jegley et al. [65] proved that the load carrying capability of panels reinforced by curvilinear fibers was superior to the one that characterized panels with straight fibers.Analogously, the stress concentrations around the hole were reduced due to the curvilinear paths of the fibers.Abdalla et al. [66] developed a maximization procedure for the natural frequency of composite panels reinforced by curved fibers.The post-buckling progressive damage behavior and the structural failure of variable-stiffness composite panels were studied by Lopes et al. [67] taking into account the residual thermal stresses caused by the curing process.Gürdal et al. [68] investigated the in-plane and buckling responses of flat rectangular composite laminates with variable stiffness.Their research proved that it was possible to vary either the buckling load or the in-plane stiffness of the structure.The variable-stiffness concept and the two-parameter law for the fiber orientation were extended also to conical and cylindrical shells by Blom et al. [69][70][71].Akhavan and Ribeiro [72] investigated the natural frequencies and vibrational modes of variable-stiffness laminated composite plates reinforced by curvilinear fibers by means of a p-version finite element approach.The same topic was analyzed also by Honda and Narita [73], who proposed an analytical method for this purpose.Díaz et al. [74] presented a numerical method for obtaining the interlaminar stresses in variable stiffness composite panels.Coburn et al. [75] highlighted the improved buckling performances of variable-stiffness laminates compared to the ones reinforced by straight fibers.They employed the Rayleigh-Ritz procedure in the theoretical framework of the First-order Shear Deformation Theory (FSDT).Raju et al. [76] performed a geometrically nonlinear analysis of symmetric laminated composite plates with curvilinear fibers subjected to in-plane shear loads.In particular, the effect of fiber orientation on the buckling and post-buckling behavior was investigated.The same kind of analysis was performed by White et al. [77], who extended the research to variable-stiffness cylinders under axial compression forces.Yazdani and Ribeiro [78] developed a layer-wise p-version finite element approach able to capture accurately the free vibrations of thick composite laminates reinforced by curvilinear fibers.Ribeiro and Stoykov [79] and Ribeiro [80,81] investigated respectively the free and forced geometrically nonlinear vibrations of variable-stiffness cylindrical shells.Analogously, Akhavan and Ribeiro [82] developed a third-order shear deformation theory to study the geometrically nonlinear periodic forced vibrations of imperfect laminates reinforced by curvilinear fibers.To this aim, a p-version finite element method was used.An accurate variable-kinematic approach was developed instead by Vescovini and Dozio [83] for the vibration and buckling analysis of moderately thick variable-stiffness laminated and sandwich plates.Groh and Weaver [84] proposed a pseudo-spectral differential quadrature method to analyze the mechanical behavior of sandwich panels with variable-stiffness face sheets, whose governing equations are derived by using the Hellinger-Reissner model.On the other hand, the Koiter's approach and the finite element method were employed in the paper by Madeo et al. [85] to analyze the post-buckling behavior of variable-stiffness composite plates.Khani et al. [86] investigated failure loads, failure modes and weights of variable-stiffness panels with cut-out by means of several experimental tests.Their work proved that variable-stiffness laminates could sustain larger loads before failure.Demasi et al. [87] proposed a multi-theory approach based on the Generalized Unified Formulation (GUF), including Equivalent Single Layer, zig-zag, and Layer-Wise models, for the static analysis of variable-stiffness composites.The Rayleigh-Ritz method was employed by Oliveri et al. [88] to investigate the thermo-mechanical post-buckling behavior of variable-stiffness multilayered plates, taking into account geometrical nonlinearities through the Von Kármán's hypotheses.A theoretical framework based on higher-order theories was developed by Tornabene et al. [89,90] for the free vibration and static analyses of doubly-curved laminated composite shells reinforced by curvilinear fibers placed according to classical schemes.A numerical solution was proposed in these circumstances.
Finally, it should be mentioned that the presence of defects or imperfections, such as gaps and overlaps, characterizes laminated composites reinforced by curvilinear fibers, due to the peculiar manufacturing process, as highlighted in the papers [91][92][93][94].This aspect could be even more evident in shell structures due to their curvature.For this purpose, an innovative fiber placement method named as Continuous Tow Shearing (CTS) was recently developed by Kim et al. [95] to minimize such defects in the manufacture of variable-stiffness composites.
A structural model based on an ESL approach is developed in this paper to investigate the linear static behavior of doubly-curved shells reinforced by curvilinear fibers placed according to a general scheme.In particular, several Higher-order Shear Deformation Theories (HSDTs) are developed to this aim [96,97].Since soft-core sandwich structures are analyzed, the Murakami's function is embedded in the model to capture the effective zig-zag effect that occurs in this peculiar circumstance [98][99][100].This ESL approach can reach the same accuracy of a layer-wise theory without increasing the degrees of freedom of the problem and the computational cost, consequently.Further details concerning layer-wise theories could be found in the paper by Tornabene [12] and by Naumenko and Eremeyev [101].In addition, as highlighted in the papers [102][103][104][105][106][107][108][109][110][111][112][113][114], lower-order theories such as the well-known FSDT could be inadequate in these cases.Once the theoretical framework is introduced, a massive set of parametric investigations is presented to show the effect of the fiber orientation on the static response of these structural elements.The Generalized Differential Quadrature (GDQ) method is employed to get the numerical solution of the governing equations [115][116][117][118].

Shell Geometry
Let us introduce a global reference coordinate system Ox 1 x 2 x 3 in the three-dimensional space, as shown in Figure 1.Each direction x i is clearly identified by the corresponding unit vector e 1 .A generic shell element is defined in such space.Its shape can be described through the position vector of the middle surface r(α 1 , α 2 ), where α 1 , α 2 are the principal and orthogonal curvilinear coordinates of the surface in hand (Figure 1), for reinforced by curvilinear fibers placed according to classical schemes.A numerical solution was proposed in these circumstances.Finally, it should be mentioned that the presence of defects or imperfections, such as gaps and overlaps, characterizes laminated composites reinforced by curvilinear fibers, due to the peculiar manufacturing process, as highlighted in the papers [91][92][93][94].This aspect could be even more evident in shell structures due to their curvature.For this purpose, an innovative fiber placement method named as Continuous Tow Shearing (CTS) was recently developed by Kim et al. [95] to minimize such defects in the manufacture of variable-stiffness composites.
A structural model based on an ESL approach is developed in this paper to investigate the linear static behavior of doubly-curved shells reinforced by curvilinear fibers placed according to a general scheme.In particular, several Higher-order Shear Deformation Theories (HSDTs) are developed to this aim [96,97].Since soft-core sandwich structures are analyzed, the Murakami's function is embedded in the model to capture the effective zig-zag effect that occurs in this peculiar circumstance [98][99][100].This ESL approach can reach the same accuracy of a layer-wise theory without increasing the degrees of freedom of the problem and the computational cost, consequently.Further details concerning layer-wise theories could be found in the paper by Tornabene [12] and by Naumenko and Eremeyev [101].In addition, as highlighted in the papers [102][103][104][105][106][107][108][109][110][111][112][113][114], lower-order theories such as the well-known FSDT could be inadequate in these cases.Once the theoretical framework is introduced, a massive set of parametric investigations is presented to show the effect of the fiber orientation on the static response of these structural elements.The Generalized Differential Quadrature (GDQ) method is employed to get the numerical solution of the governing equations [115][116][117][118].

Shell Geometry
Let us introduce a global reference coordinate system Ox x x in the three-dimensional space, as shown in Figure 1.Each direction i x is clearly identified by the corresponding unit vector 1 e .A generic shell element is defined in such space.Its shape can be described through the position vector of the middle surface ( ) , where 1 2 , α α are the principal and orthogonal curvilinear coordinates of the surface in hand (Figure 1), for  The coordinate ζ denotes the thickness direction.Thus, O α 1 α 2 ζ stands for the local reference system of the shell middle surface.The shell volume is identified by its constant thickness h.It should be recalled that the structure is made of l layers, therefore the overall thickness is given by where h k stands for the thickness of the k-th layer, which can be computed as follows Since only sandwich structures are analyzed in the present paper, the shells are made of three layers: two external thinner plies of equal thickness h s , and a thicker core of thickness h c .Thus, the total thickness is given by Each point of this three-dimensional structure is analytically described through the following position vector where n(α 1 , α 2 ) represents the outward unit normal vector.This quantity can be computed by means of the cross product (denoted by "×") as follows in which r ,i = ∂r/∂α i .Once the position vector r(α 1 , α 2 ) of the middle surface is given, the Lamè parameters A 1 , A 2 of the doubly-curved surface which define the shape of the shell can be easily evaluated by means of the dot product (denoted by "•") as shown below Analogously, the radii of curvature R 1 , R 2 of the shell middle surface can be calculated through the following relations

Higher-Order ESL Model
The present higher-order structural model allows the study of several displacement fields by varying the maximum order of kinematic expansion N. The three-dimensional displacements ζ) can be written as follows in which F τ (ζ) are the thickness-functions that define the kinematic expansion.The generalized displacements u (τ) 3 (α 1 , α 2 ), which are defined on the shell middle surface, are the degrees of freedom.The thickness functions are assumed as power law functions.In other words, one gets for τ = 0, 1, 2, . . ., N. The (N + 1)-th function, instead, denotes the Murakami's function Z(ζ) defined below This peculiar function must be inserted in the structural model to capture the effective behavior of soft-core sandwich structures.Otherwise, the so-called zig-zag effect cannot be properly modeled [90].The notation EDZN is used to identify the shear-deformation theory including the Murakami's function.
The generalized strain components, which are defined on the shell middle surface as well, can be included in the corresponding vector ε (τ) (α 1 , α 2 ) shown below for τ = 0, 1, 2, . . ., N, N + 1.These quantities can be evaluated as a function of the generalized displacements as follows in which D Ω is the kinematic differential operator defined below The three-dimensional strain components can be included in the vector ε(α 1 , α 2 , ζ), which assumes the aspect below These quantities can be computed once the generalized strains are evaluated.The following definition can be used for this purpose where the matrix Z (τ) (α 1 , α 2 , ζ) has the meaning shown below for each order τ of kinematic expansion The geometric parameters H 1 , H 2 can be evaluated as a function of the principal radii of curvature Consequentially, the three-dimensional stress components can be evaluated as a function of the strain components by using the elastic constitutive laws, assuming a linear behavior of the composite materials.The stress components at issue can be collected in the vector σ(α 1 , α 2 , ζ) defined below The following relation written in compact matrix form allows a definition of such quantities as a function of the corresponding stress components where C (k) is the constitutive operator for the k-th orthotropic layer defined in the local reference system of the structure.It assumes the following aspect Quantities included in C (k) take into account the orientation of the reinforcing fibers θ (k) for the k-th ply.The following relations can be used to evaluate them where c = cos θ (k) and s = sin θ (k) .The elastic constants in ( 23)-( 35) can be computed as a function of the engineering constants of the materials of the k-th layer as shown below (36) where the quantity ∆ (k) has the following meaning in which 23 are the shear moduli, and ν are the Poisson's ratios.The following relations are required for a complete mechanical characterization of the orthotropic layer, for i, j = 1, 2, 3 The evaluation of the elastic coefficients is simplified if an isotropic layer is considered, since 23 , and the mechanical properties are not affected by the orientation of the lamina.In this circumstance, the shear modulus is given by In this paper, the core is always assumed as isotropic; therefore, its orientation is meaningless.On the other hand, the two external skins are orthotropic and characterized by a variable orientation of the fibers θ (k) .In other words, the reinforcing fibers are placed according to curvilinear paths.The value of θ (k) can be defined through the following general law (40) in which the constants θ i is given below in terms of the following dimensionless coordinates, for i = 1, 2 The functions ψ 3 denote a variation along the principal direction α 1 and can assume different meanings (power law, sine-wave, or exponential), as shown below where n 1 , . On the other hand, the functions ψ 4 define a variation along the principal direction α 2 and have the same meaning of the previous ones where 6 involve both the principal directions α 1 , α 2 and can be expressed in the following form for i = 5, 6.The Gaussian distribution can be used for this purpose and one gets where Analogously, an ellipse-shaped function could be chosen and one gets in which a, b define the semi-axes of the ellipse, which could be also rotated by a constant angle β measured counter-clockwise from the direction α 1 .This function is centered instead in the point (α 7m , α 8m ) within the reference domain, with i are activated by setting the corresponding value θ (k) i different from zero.It should be noted that several variations can be also combined to obtain a more complex path of the reinforcing fibers.
The present formulation does not consider any defect or imperfection that could occur in the manufacturing process.Such defects are unavoidable when the fibers are placed according to curvilinear paths, as specified in the papers by Blom et al. [91,92], Nik et al. [93] and Akbarzadeh et al. [94].Nevertheless, several production techniques are developed to reduce these imperfections [95].Thus, in this research paper it is assumed by hypothesis that such defects that occur during the manufacturing process are negligible.Further assumptions that involve the present formulation can be found in the recent paper by Tornabene et al. [56].
Once the constitutive laws are introduced, it is more advantageous to define the internal actions in terms of generalized stress resultants S (τ) (α 1 , α 2 ) as shown below for τ = 0, 1, 2, . . ., N, N + 1 The vector that collects the generalized stress resultants, which are defined on the shell middle surface, assumes the following aspect for each order of kinematic expansion τ On the other hand, the matrix A (τη) in (47) represents the constitutive operator of the present higher-order approach which can be evaluated as follows for τ, η = 0, 1, 2, . . ., N, N + 1.In extended notation, the following matrix is accomplished 11( 20 where the following definitions are required to compute its elements for τ, η = 0, 1, 2, . . ., N, N + 1, n, m = 1, 2, 3, 4, 5, 6 and p, q = 0, 1, 2 Finally, the static equilibrium of the shell can be expressed in terms of generalized displacement components by means of the following governing equation in vector form for τ = 0, 1, 2, . . ., N, N + 1.This system of partial derivative equations defines three equilibrium relations for each order of kinematic expansion.The fundamental operator L (τη) can be evaluated according to the definition shown below for τ, η = 0, 1, 2, . . ., N, N + 1.The equilibrium operator D * Ω , instead, assumes the following aspect On the other hand, the vector q (τ) introduced for each order of kinematic equation in ( 52) collects three load components for τ = 0, 1, 2, . . ., N, N + 1 as shown below where q (τ) n stand for the external load components applied within the shell middle surface along the principal directions α 1 , α 2 , ζ, respectively.Their value is computed once the surface loads applied on the shell external surfaces q (±) are defined.The superscript (+) specifies that the load component is acting on the upper surface, whereas (−) is used to characterize the force component applied on the lower surface.By means of the static equivalence principle, one gets where the symbol F (±) τ denotes the value that a generic thickness function assumes on the external surfaces of the shell.Analogously, are the values that the geometric parameters (19) assume on the same surfaces (for ζ = ±h/2).
The governing equations in (52) can be solved once the natural or essential boundary conditions are properly enforced.The compatibility requirement of the present formulation is C 1 , since both the displacements and stress resultant components are involved in the boundary conditions.The essential boundary conditions involve the generalized displacements for each order of kinematic expansion τ.The following relations are required for clamped edges On the other hand, a free edge condition is expressed analytically through the natural boundary conditions.Thus, the stress resultants are involved.In particular, the following conditions are required for α 1 = α 0 1 (or On the other hand, the natural boundary conditions for With reference to the shell element depicted in Figure 1, it can be observed that the four external edges are identified by a cardinal direction: north (N), east (E), south (S), and west (W).To simplify the definition of external restraints, the boundary conditions are specified following the sequence WSEN.The letters "C" and "F" will replace the corresponding cardinal direction notation to define a clamped or free edge, respectively.For instance, CCCC stands for a completely clamped shell.On the other hand, when the notation CCCF is used, only the northern edge is free, whereas the others are clamped.
When a structure with a closing meridian such as a complete shell of revolution is analyzed, the structural compatibility must be enforced along the two coincident edges.If the closing meridian is identified by α 1 = α 0 1 and α 1 = α 1 1 , the structural compatibility is given by On the other hand, the following relations are required if the closing edge is characterized by

Numerical Technique
The governing equations are solved numerically by means of the GDQ method.In this section, only the main features of this numerical approach are presented, since a complete and accurate treatise is presented in the recent review paper by Tornabene et al. [117].By definition, the GDQ method allows evaluation of the n-th order derivative of a smooth function f (x) at a generic discrete point x i of the reference domain as a weighted linear sum of the values that the function assumes in each point of the domain where T stands for the overall number of discrete points defined within the domain, which are placed according to the Chebyshev-Gauss-Lobatto grid distribution Consequently, a set of T discrete nodes is introduced as follows where x 1 , x T are the boundary points of the domain.The weighting coefficients ς (n) ij are computed instead through the following recursive expressions provided by Shu [115], for i, j = 1, 2, . . ., T and n = 1, 2, . . ., T − 1 where ς (1) ij stands for the weighting coefficients for the first-order derivatives defined below As far as two-dimensional domains are concerned, the present approach can be easily generalized as shown in [23,117].In these circumstances, the grid distribution ( 63) is applied along the main directions α 1 , α 2 by specifying the number of discrete points as I N , I M , respectively.
The same numerical approach can be considered as the starting point for the evaluation of integrals.The methodology in hand is known as Generalized Integral Quadrature (GIQ) method.For the sake of completeness, only the main definition is presented in the paper.Further details concerning this technique can be found in [117].The integral of a smooth function f (x) within the interval x i , x j included in the reference domain, with x i , x j ∈ [a, b] and x i > a, x j < b, can be computed as follows where w ij k stands for the weighting coefficients for the integration at issue, which can be easily evaluated starting from the weighting coefficients for the first-order derivatives (66).This numerical method is used in the present research to compute the stiffness terms (51), once the thickness direction is discretized by applying the same grid shown above.

Solution of the Static Problem
By means of the GDQ method, the discrete form of the fundamental system of Equation ( 53) is achieved, which becomes where K is the global discrete stiffness matrix, δ the vector of generalized displacements related to the discrete points of the domain, and f is the vector that collects the values of the external loads.It should be noted that the unknowns of the static problem at issue are collected in the corresponding algebraic vector δ.It is convenient at this point to collect such displacements so that the degrees of freedom linked to the inner points of the domain are separated from the ones related to the boundary points.The sub-vectors δ d , δ b are introduced, respectively.Analogously, the same classification is specified for the external load vector (f d , f b ).Thus, the system of linear Equation (68) becomes in which the stiffness matrix is also partitioned, consequently.Finally, the proper manipulations lead to the following result, which allows the obtaining of the static solution δ d in terms of generalized nodal displacements The procedure just illustrated represents the so-called static condensation procedure and allows the reduction of the initial size of the problem [4].
Once the static solution in terms of generalized displacements is accomplished, the through-thethickness profiles of stress and strain components can be evaluated by using the three-dimensional elasticity equations written in an orthogonal curvilinear coordinate system.This approach is known as recovery procedure and can be employed to evaluate the three-dimensional quantities at issue, starting from the solutions obtained by a two-dimensional shell model.The reader can find a complete treatise about this approach in the papers [21].
Finally, it should be mentioned that the present formulation is implemented in a computational code developed in MATLAB by the authors [119].

Applications
The numerical results are presented in this section.For this purpose, four structures reinforced by curvilinear fibers are analyzed.Their geometries are depicted in Figure 2. In particular, a square plate, a conical shell, a doubly-curved panel of revolution with catenary meridian and a doubly-curved panel of translation (obtained by moving a circumference over a parabola) are investigated.Their shape is described through the position vectors listed in Table 1.Analogously, the lamination schemes, as well as the mathematical laws that describe the curvilinear paths of the reinforcing fibers, are specified in Table 1.For the sake of clarity, the curvilinear paths are also presented in graphical form in Figure 3 for the lower and upper sheets of each structure, considering the maximum value of the fiber orientation.For conciseness purposes, the mechanical properties of the orthotropic skins and the isotropic core are listed in Table 2. On the other hand, the mechanical composition of the four structures, which are all loaded by an external pressure q (+) n = −10 kPa, is summarized in Table 3.The results are presented in terms of through-the-thickness profiles of strain, stress, and displacement components.Then, the profiles of the generalized displacements u 3 along two couples of orthogonal directions are also provided for each structure.
A set of parametric investigations is presented to study the effect of the fiber orientation represented by the symbol φ ∈ [0 • , 60 • ] on the linear static response.Each group of analyses is repeated to investigate also the influence of the maximum order of kinematic expansion.To this aim, three higher-order theories are considered, which are EDZ2, EDZ3 and EDZ4.
Finally, it should be mentioned that the same structures have been analyzed from the dynamic point of view (free vibration analysis) in the previous paper by Tornabene et al. [56].
The results are presented in terms of through-the-thickness profiles of strain, stress, and displacement components.Then, the profiles of the generalized displacements , , u u u along two couples of orthogonal directions are also provided for each structure.
A set of parametric investigations is presented to study the effect of the fiber orientation represented by the symbol [ ] 0 ,60 φ ∈ °° on the linear static response.Each group of analyses is repeated to investigate also the influence of the maximum order of kinematic expansion.To this aim, three higher-order theories are considered, which are EDZ2 , EDZ3 and EDZ4 .
Finally, it should be mentioned that the same structures have been analyzed from the dynamic point of view (free vibration analysis) in the previous paper by Tornabene et al. [56].The results are presented in terms of through-the-thickness profiles of strain, stress, and displacement components.Then, the profiles of the generalized displacements ( ) ( ) ( ) , , u u u along two couples of orthogonal directions are also provided for each structure.
A set of parametric investigations is presented to study the effect of the fiber orientation represented by the symbol [ ] 0 ,60 φ ∈ °° on the linear static response.Each group of analyses is repeated to investigate also the influence of the maximum order of kinematic expansion.To this aim, three higher-order theories are considered, which are EDZ2 , EDZ3 and EDZ4 .Finally, it should be mentioned that the same structures have been analyzed from the dynamic point of view (free vibration analysis) in the previous paper by Tornabene et al. [56].The discrete model of a completely clamped (CCCC) square plate is obtained by setting I N = I M = 25 grid points.In this circumstance, a three-dimensional finite element model (3D-FEM) is also developed for comparison purposes, only for straight reinforcing fibers (φ = 0).The commercial software Strand7 is employed for this reason, by using the Hexa20 brick elements.
The through-the-thickness variations of the strain, stress and displacement components are evaluated in 0.25L x , 0.25L y and are shown in Figures 4-6.On the other hand, the profiles of the generalized displacements are evaluated for constant values of the main coordinates, varying the values of the fiber orientation φ and the maximum order of kinematic expansion.These results are presented in Figures 7 and 8.It can be noted that the variation of the fiber orientation mainly affects the generalized displacements defined on the middle surface than the three-dimensional profiles of strains, stresses and displacements, since the soft-core effect is prevalent.Nevertheless, also the stress and strain components are influenced by the placement of the fibers, especially for higher values of [ ] 0 ,60 φ ∈ °°.This example proves also the validity of the current formulation, both in terms of two-dimensional and three-dimensional solutions, through the great agreement with the 3D-FEM results.These results are presented in Figures 7 and 8.It can be noted that the variation of the fiber orientation mainly affects the generalized displacements defined on the middle surface than the three-dimensional profiles of strains, stresses and displacements, since the soft-core effect is prevalent.Nevertheless, also the stress and strain components are influenced by the placement of the fibers, especially for higher values of φ ∈ [0 • , 60 • ].This example proves also the validity of the current formulation, both in terms of two-dimensional and three-dimensional solutions, through the great agreement with the 3D-FEM results.These results are presented in Figures 7 and 8.It can be noted that the variation of the fiber orientation mainly affects the generalized displacements defined on the middle surface than the three-dimensional profiles of strains, stresses and displacements, since the soft-core effect is prevalent.Nevertheless, also the stress and strain components are influenced by the placement of the fibers, especially for higher values of  3 for constant values of the main coordinates can be obtained by varying the fiber orientation, as it can be noted in Figures 12 and 13.The influence of the curvilinear path is especially evident along the circumferential lines obtained for constant values of y coordinate (Figure 12), where the corresponding displacement profiles are represented by straight segments.
, , u u u of a completely clamped square plate evaluated for constant value of x coordinate, varying the values of the fiber orientation φ and the maximum order of kinematic expansion N : (a)

Conical Shell
A completely clamped conical shell is considered in the next application.The discrete model is obtained by setting , , u u u for constant values of the main coordinates can be obtained by varying the fiber orientation, as it can be noted in Figures 12 and 13.The influence of the curvilinear path is especially evident along the circumferential lines obtained for constant values of y coordinate (Figure 12), where the corresponding displacement profiles are represented by straight segments.( )

Doubly-Curved Panel of Revolution
The doubly-curved panel investigated in this paragraph is obtained by the rotation of a catenary shaped profile about the axis of revolution.The corresponding discrete model is obtained by setting φ ∈ °° and the maximum order of kinematic expansion.Even in this circumstance, the static response is affected by the fiber orientation.
On the other hand, the profiles of the generalized displacements within the shell middle surface are collected in Figures 17 and 18 for constant values of the principal coordinates.
3 of a completely clamped conical shell evaluated for constant value of α 1 coordinate, varying the values of the fiber orientation φ and the maximum order of kinematic expansion N:

Doubly-Curved Panel of Revolution
The doubly-curved panel investigated in this paragraph is obtained by the rotation of a catenary shaped profile about the axis of revolution.The corresponding discrete model is obtained by setting I N = 31 and I M = 25 as in the previous example.The boundary conditions are specified by the notation CCCF; in other words, the northern edge is free, whereas the others are clamped.
The through-the-thickness variations of three-dimensional strains, stresses and displacements are computed in the point (0.25(ϕ 1 − ϕ 0 ) + ϕ 0 , 0.25(ϑ 1 − ϑ 0 ) + ϑ 0 ).The corresponding graphs are illustrated in Figures 14-16, varying the values of the fiber orientation φ ∈ [0 • , 60 • ] and the maximum order of kinematic expansion.Even in this circumstance, the static response is affected by the fiber orientation.illustrated in Figures 14-16, varying the values of the fiber orientation [ ] 0 ,60 φ ∈ °° and the maximum order of kinematic expansion.Even in this circumstance, the static response is affected by the fiber orientation.
On the other hand, the profiles of the generalized displacements within the shell middle surface are collected in Figures 17 and 18 for constant values of the principal coordinates.

Doubly-Curved Panel of Translation
The last structure is a completely clamped doubly-curved panel of translation.The middle surface of this peculiar geometry can be described by sliding a circumference over a parabolic arch.As far as the discrete model is concerned, the grid distribution is applied by setting Finally, the deformed shapes of the four shell structures analyzed in this section are depicted in Figure 24 for three different values of the fiber orientation ( 0 , 30 , 60 2 , u 3 of a doubly-curved panel of revolution evaluated for constant value of ϕ coordinate, varying the values of the fiber orientation φ and the maximum order of kinematic expansion N: (a) ϕ = 0.25(ϕ 1 − ϕ 0 ) + ϕ 0 ; (b) ϕ = 0.5(ϕ 1 − ϕ 0 ) + ϕ 0 .

Doubly-Curved Panel of Translation
The last structure is a completely clamped doubly-curved panel of translation.The middle surface of this peculiar geometry can be described by sliding a circumference over a parabolic arch.As far as the discrete model is concerned, the grid distribution is applied by setting I N = I M = 31.
The through-the-thickness profiles of strain, stress and displacement components are evaluated respectively in 0.25 , varying the values of the fiber orientation φ ∈ [0 • , 60 • ] and the maximum order of kinematic expansion N. Their graphs are shown in Figures 19-21.The variations of the generalized displacement components computed on the shell middle surface are shown instead in Figures 22 and 23, by selecting constant values of the principal coordinates.
Finally, the deformed shapes of the four shell structures analyzed in this section are depicted in Figure 24 for three different values of the fiber orientation (φ = 0 • , φ = 30 • , φ = 60 • ).
2 , u 3 of a doubly-curved panel of translation evaluated for constant value of α 2 coordinate, varying the values of the fiber orientation φ and the maximum order of kinematic expansion N: (a) (b) , , u u u of a completely-clamped doubly-curved panel of translation evaluated for constant value of 1 α coordinate, varying the values of the fiber orientation φ and the maximum order of kinematic expansion N : (a) ( ) 3 of a completely-clamped doubly-curved panel of translation evaluated for constant value of α 1 coordinate, varying the values of the fiber orientation φ and the maximum order of kinematic expansion N: , , u u u of a completely-clamped doubly-curved panel of translation evaluated for constant value of 1 α coordinate, varying the values of the fiber orientation φ and the maximum order of kinematic expansion N : (a) ( )

Remarks on the Use of the Murakami's Function
It can be noted that the present approach is used to perform the static analysis of sandwich structures made of three layers (two external stiffer skins and a soft-core).Murakami's function can capture their effective mechanical behavior in this circumstance.Nevertheless, sandwich structures could be characterized by more than three layers.For instance, the inner soft-core could be enclosed by couples of orthotropic stiffer layers.Consequently, a lamination scheme could be given by work by Groh and Weaver [121], the Murakami's function embedded in an ESL model could be inadequate to deal with these kinds of structural problems.A final application is presented to illustrate this aspect.For this purpose, a completely clamped laminated square plate with 1.5 m

Remarks on the Use of the Murakami's Function
It can be noted that the present approach is used to perform the static analysis of sandwich structures made of three layers (two external stiffer skins and a soft-core).Murakami's function can capture their effective mechanical behavior in this circumstance.Nevertheless, sandwich structures could be characterized by more than three layers.For instance, the inner soft-core could be enclosed by couples of orthotropic stiffer layers.Consequently, a lamination scheme could be given by θ (1) /θ (2) /Soft − core/θ (4) /θ (5) .As highlighted first in the paper by Gherlone [120], then in the work by Groh and Weaver [121], the Murakami's function embedded in an ESL model could be inadequate to deal with these kinds of structural problems.A final application is presented to illustrate this aspect.For this purpose, a completely clamped laminated square plate with L x = L y = 1.5 m and h = 0.1 m, loaded by an external pressure q (+) n = −10 kPa, is considered.Its lamination scheme is given by (0/90/Soft − core/90/0), assuming h 1 = h 2 = 0.01 m, h 4 = h 5 = 0.01 m and h 3 = 0.06 m.
The four external layers are made of orthotropic material (Graphite-epoxy), whereas the soft-core is made of Foam.Their properties are listed in Table 2.In this circumstance, the reinforcing fibers are placed according to straight paths to perform the comparison also with the three-dimensional finite element solution (3D FEM).As far as the present ESL approach is concerned, the EDZ4 is considered to be structural theory.To prove the inadequacy of the Murakami's function to capture the effective static behavior of this five-layered plate, a layer-wise theory is used for the sake of comparison.In particular, the LD4 model is used.The reader can find the meaning and the features of this theory in the paper by Tornabene [12].The through-the-thickness profiles of strains, stresses and displacements are shown in Figures 25-27.The following results prove the same aspects highlighted in the reference papers [120,121].In particular, the LD4 is in excellent agreement with the reference solution (3D FEM), whereas the EDZ4 is not able to capture the effective through-the-thickness profiles, especially in terms of strains and displacements.The discrete model is obtained by setting I N = I M = 25 grid points for the GDQ solutions.On the other hand, the commercial software Strand7 is employed for the 3D FEM solution, by using the Hexa20 brick elements.This final example justifies the use of the ESL zig-zag theories only for three-layered sandwich structures with an inner soft-core.and 3 0.06 m h = .The four external layers are made of orthotropic material (Graphite-epoxy), whereas the soft-core is made of Foam.Their properties are listed in Table 2.In this circumstance, the reinforcing fibers are placed according to straight paths to perform the comparison also with the three-dimensional finite element solution (3D FEM).As far as the present ESL approach is concerned, the EDZ4 is considered to be structural theory.To prove the inadequacy of the Murakami's function to capture the effective static behavior of this five-layered plate, a layer-wise theory is used for the sake of comparison.In particular, the LD4 model is used.The reader can find the meaning and the features of this theory in the paper by Tornabene [12].The through-the-thickness profiles of strains, stresses and displacements are shown in Figures 25-27.The following results prove the same aspects highlighted in the reference papers [120,121].In particular, the LD4 is in excellent agreement with the reference solution (3D FEM), whereas the EDZ4 is not able to capture the effective through-the-thickness profiles, especially in terms of strains and displacements.The discrete model is obtained by setting 25 grid points for the GDQ solutions.On the other hand, the commercial software Strand7 is employed for the 3D FEM solution, by using the Hexa20 brick elements.This final example justifies the use of the ESL zig-zag theories only for three-layered sandwich structures with an inner soft-core.

Discussion
A massive set of parametric investigations has been performed and presented in the paper to show the effect of the curvilinear paths of the reinforcing fibers on the linear static response of soft-core sandwich structures.A theoretical framework based on higher-order ESL theories including the Murakami's function has been employed.As far as the values of the fiber orientation are concerned, an innovative and general mathematical formulation has been developed, to describe more complex curvilinear paths (such as power law, sine-wave, exponential, Gaussian and ellipse shaped).The following observations can be gathered:

•
The soft-core effect is well-captured by means of the Murakami's function.As highlighted in previous work [56], this function is required to model the so-called zig-zag effect.These aspects are extremely clear in the through-the-thickness profiles of strain, stress and displacement components presented in the paper.

•
The higher-order theories employed in this paper provide comparable results.Thus, all these models can be used to deal with similar structural problems.With respect to first-order models, the shear correction factor and the plane stress hypothesis can be neglected, and the mechanical behavior is closer to the three-dimensional one.

•
The linear static behavior is affected by the value of the fiber orientation.Therefore, the angle of the fiber orientation represents a design parameter that can be modified and optimized during the manufacturing process to obtain the desired structural behavior in terms of static response.

•
The boundary conditions in terms of generalized displacements on the shell middle surface, as well as stress profiles, are well-enforced.In particular, null displacements are obtained for clamped edges, whereas the normal stress n σ coincides with the value of the corresponding applied load.

•
The numerical method employed in the paper represents an accurate tool to deal with similar structural problems.

27.
Through-the-thickness variation of the displacement components [m] for a completely clamped square plate evaluated in 0.25L x , 0.25L y .Comparison among the finite element method (3D FEM), layer-wise (LD4) and ESL approaches (EDZ4).

Discussion
A massive set of parametric investigations has been performed and presented in the paper to show the effect of the curvilinear paths of the reinforcing fibers on the linear static response of soft-core sandwich structures.A theoretical framework based on higher-order ESL theories including the Murakami's function has been employed.As far as the values of the fiber orientation are concerned, an innovative and general mathematical formulation has been developed, to describe more complex curvilinear paths (such as power law, sine-wave, exponential, Gaussian and ellipse shaped).The following observations can be gathered:

•
The soft-core effect is well-captured by means of the Murakami's function.As highlighted in previous work [56], this function is required to model the so-called zig-zag effect.These aspects are extremely clear in the through-the-thickness profiles of strain, stress and displacement components presented in the paper.

•
The higher-order theories employed in this paper provide comparable results.Thus, all these models can be used to deal with similar structural problems.With respect to first-order models, the shear correction factor and the plane stress hypothesis can be neglected, and the mechanical behavior is closer to the three-dimensional one.

•
The linear static behavior is affected by the value of the fiber orientation.Therefore, the angle of the fiber orientation represents a design parameter that can be modified and optimized during the manufacturing process to obtain the desired structural behavior in terms of static response.

•
The boundary conditions in terms of generalized displacements on the shell middle surface, as well as stress profiles, are well-enforced.In particular, null displacements are obtained for clamped edges, whereas the normal stress σ n coincides with the value of the corresponding applied load.

•
The numerical method employed in the paper represents an accurate tool to deal with similar structural problems.

i
, for i = 1, 2, . . ., 6, represent the maximum values of the corresponding fiber variations denoted by the functions ψ (k) i .On the other hand, θ (k) is the constant orientation that define the initial straight placement of the reinforcing fibers.The meaning of the six general functions ψ (k)

Figure 3 .
Figure 3. Curvilinear fiber paths in the external sheets of the structures: (a) Square plate; (b) Conical shell; (c) Doubly-curved panel of revolution; (d) Doubly-curved panel of translation.

Figure 3 .
Figure 3. Curvilinear fiber paths in the external sheets of the structures: (a) Square plate; (b) Conical shell; (c) Doubly-curved panel of revolution; (d) Doubly-curved panel of translation.
The discrete model of a completely clamped (CCCC) square plate is obtained by setting 25 points.In this circumstance, a three-dimensional finite element model (3D-FEM) is also developed for comparison purposes, only for straight reinforcing fibers ( 0 φ = ).The commercial software Strand7 is employed for this reason, by using the Hexa20 brick elements.The through-the-thickness variations of the strain, stress and displacement components are evaluated in ( ) shown in Figures 4-6.On the other hand, the profiles of the generalized displacements are evaluated for constant values of the main coordinates, varying the values of the fiber orientation φ and the maximum order of kinematic expansion.

Figure 4 .
Figure 4. Through-the-thickness variation of the strain components for a completely clamped square plate evaluated in ( ) 0.25 ,0.25 x y L L , varying the values of the fiber orientation φ and the maximum order of kinematic expansion N .Comparison with the finite element method (3D FEM).

Figure 4 .
Figure 4. Through-the-thickness variation of the strain components for a completely clamped square plate evaluated in 0.25L x , 0.25L y , varying the values of the fiber orientation φ and the maximum order of kinematic expansion N. Comparison with the finite element method (3D FEM).

Figure 5 .
Figure 5. Through-the-thickness variation of the stress components [Pa] for a completely clamped square plate evaluated in ( ) 0.25 ,0.25 x y L L , varying the values of the fiber orientation φ and the maximum order of kinematic expansion N .Comparison with the finite element method (3D FEM).

Figure 5 . 43 Figure 5 .
Figure 5. Through-the-thickness variation of the stress components [Pa] for a completely clamped square plate evaluated in 0.25L x , 0.25L y , varying the values of the fiber orientation φ and the maximum order of kinematic expansion N. Comparison with the finite element method (3D FEM).

Figure 6 .
Figure 6.Through-the-thickness variation of the displacement components [m] for a completely clamped square plate evaluated in ( ) 0.25 ,0.25 x y L L , varying the values of the fiber orientation φ and the maximum order of kinematic expansion N .Comparison with the finite element method (3D FEM).

Figure 6 .
Figure 6.Through-the-thickness variation of the displacement components [m] for a completely clamped square plate evaluated in 0.25L x , 0.25L y , varying the values of the fiber orientation φ and the maximum order of kinematic expansion N. Comparison with the finite element method (3D FEM).

J 43 Figure 6 .
Figure 6.Through-the-thickness variation of the displacement components [m] for a completely clamped square plate evaluated in ( ) 0.25 ,0.25 x y L L , varying the values of the fiber orientation φ and the maximum order of kinematic expansion N .Comparison with the finite element method (3D FEM).
. This example proves also the validity of the current formulation, both in terms of two-dimensional and three-dimensional solutions, through the great agreement with the 3D-FEM results.

31 N I = and 25 MI
= .The through-the-thickness profiles of strain, stress and displacement components are evaluated in the point 9-11.As it can be noted from these figures, the fiber orientation affects the static response of the sandwich structure at issue.Analogously, several profiles of the generalized displacements

Figure 9 .
Figure 9. Through-the-thickness variation of the strain components for a completely clamped conical shell evaluated in

Figure 9 .
Figure 9. Through-the-thickness variation of the strain components for a completely clamped conical shell evaluated in 0.25 α 1 1 − α 0 1 + α 0 1 , 0.25(y 1 − y 0 ) + y 0 , varying the values of the fiber orientation φ and the maximum order of kinematic expansion N.

Figure 9 .
Figure 9. Through-the-thickness variation of the strain components for a completely clamped conical shell evaluated in

Figure 10 .
Figure 10.Through-the-thickness variation of the stress components [Pa] for a completely clamped conical shell evaluated in

Figure 11 .
Figure 11.Through-the-thickness variation of the displacement components [m] for a completely clamped conical shell evaluated in

Figure 11 .Figure 12 .
Figure 11.Through-the-thickness variation of the displacement components [m] for a completely clamped conical shell evaluated in

Figure 12 .Figure 13 .
Figure 12.Profiles of generalized displacements ( ) ( ) ( )0 0 0 1 2 3, , u u u of a completely clamped conical shell evaluated for constant value of y coordinate, varying the values of the fiber orientation φ and the maximum order of kinematic expansion N : (a)

31 N I = and 25 MI
= as in the previous example.The boundary conditions are specified by the notation CCCF; in other words, the northern edge is free, whereas the others are clamped.The through-the-thickness variations of three-dimensional strains, stresses and displacements are computed in the point corresponding graphs are illustrated in Figures 14-16, varying the values of the fiber orientation [ ] 0 ,60

Figure 13 .
Figure 13.Profiles of generalized displacements u

Figure 14 .
Figure 14.Through-the-thickness variation of the strain components for a doubly-curved panel of revolution evaluated in

Figure 14 .
Figure 14.Through-the-thickness variation of the strain components for a doubly-curved panel of revolution evaluated in

Figure 15 .
Figure 15.Through-the-thickness variation of the stress components [Pa] for a doubly-curved panel of revolution evaluated in

Figure 16 .
Figure 16.Through-the-thickness variation of the displacements components [m] for a doubly-curved panel of revolution evaluated in

Figure 15 .
Figure 15.Through-the-thickness variation of the stress components [Pa] for a doubly-curved panel of revolution evaluated in

Figure 16 .
Figure 16.Through-the-thickness variation of the displacements components [m] for a doubly-curved panel of revolution evaluated in

Figure 17 .
Figure 17.Profiles of generalized displacements ( ) ( ) ( )0 0 0 1 2 3, , u u u of a doubly-curved panel of revolution evaluated for constant value of ϑ coordinate, varying the values of the fiber orientation φ and the maximum order of kinematic expansion N : (a)
the-thickness profiles of strain, stress and displacement components are evaluated respectively in φ ∈ °° and the maximum order of kinematic expansion N .Their graphs are shown in Figures19-21.The variations of the generalized displacement components computed on the shell middle surface are shown instead in Figures22 and 23, by selecting constant values of the principal coordinates.

Figure 19 .
Figure 19.Through-the-thickness variation of the strain components for a doubly-curved panel of translation evaluated in

Figure 19 .
Figure 19.Through-the-thickness variation of the strain components for a doubly-curved panel of translation evaluated in

Figure 20 .
Figure 20.Through-the-thickness variation of the stress components [Pa] for a doubly-curved panel of translation evaluated in

Figure 21 .Figure 22 .
Figure 21.Through-the-thickness variation of the displacement components [m] for a doubly-curved panel of translation evaluated in

2 0Figure 24 .
Figure 24.Deformed shapes of shell structures varying the values of the fiber orientation φ: (a) Square plate; (b) Conical shell; (c) Doubly-curved panel of revolution; (d) Doubly-curved panel of translation.

Figure 25 .
Figure 25.Through-the-thickness variation of the strain components for a completely clamped square plate evaluated in 0.25L x , 0.25L y .Comparison among the finite element method (3D FEM), layer-wise (LD4) and ESL approaches (EDZ4).

Figure 25 .
Figure 25.Through-the-thickness variation of the strain components for a completely clamped square plate evaluated in ( ) 0.25 ,0.25 x y L L .Comparison among the finite element method (3D FEM), layer-wise (LD4) and ESL approaches (EDZ4).

Figure 26 .
Figure 26.Through-the-thickness variation of the stress components [Pa] for a completely clamped square plate evaluated in ( ) 0.25 ,0.25 x y L L .Comparison among the finite element method (3D FEM), layer-wise (LD4) and ESL approaches (EDZ4).

Figure 26 .
Figure 26.Through-the-thickness variation of the stress components [Pa] for a completely clamped square plate evaluated in 0.25L x , 0.25L y .Comparison among the finite element method (3D FEM), layer-wise (LD4) and ESL approaches (EDZ4).

Figure 27 .
Figure 27.Through-the-thickness variation of the displacement components [m] for a completely clamped square plate evaluated in ( ) 0.25 ,0.25 x y L L .Comparison among the finite element method (3D FEM), layer-wise (LD4) and ESL approaches (EDZ4).
) A

Table 1 .
Definition of structural geometries and corresponding lamination schemes.

Table 2 .
Mechanical properties of the materials used in the computations.

Table 3 .
Mechanical composition of the four structures.