Static Analysis of Doubly-Curved Shell Structures of Smart Materials and Arbitrary Shape Subjected to General Loads Employing Higher Order Theories and Generalized Differential Quadrature Method

Francesco Tornabene,Matteo Viscoti and Rossana Dimitri

Department of Innovation Engineering,School of Engineering,University of Salento,Lecce,73100,Italy

ABSTRACT The article proposes an Equivalent Single Layer(ESL)formulation for the linear static analysis of arbitrarily-shaped shell structures subjected to general surface loads and boundary conditions.A parametrization of the physical domain is provided by employing a set of curvilinear principal coordinates.The generalized blending methodology accounts for a distortion of the structure so that disparate geometries can be considered.Each layer of the stacking sequence has an arbitrary orientation and is modelled as a generally anisotropic continuum.In addition,re-entrant auxetic three-dimensional honeycomb cells with soた-core behaviour are considered in the model.The unknown variables are described employing a generalized displacement field and pre-determined through-the-thickness functions assessed in a unified formulation.Then,a weak assessment of the structural problem accounts for shape functions defined with an isogeometric approach starting from the computational grid.A generalized methodology has been proposed to define two-dimensional distributions of static surface loads.In the same way, boundary conditions with three-dimensional features are implemented along the shell edges employing linear springs.The fundamental relations are obtained from the stationary configuration of the total potential energy, and they are numerically tackled by employing the Generalized Differential Quadrature(GDQ)method,accounting for nonuniform computational grids.In the post-processing stage,an equilibrium-based recovery procedure allows the determination of the three-dimensional dispersion of the kinematic and static quantities.Some case studies have been presented,and a successful benchmark of different structural responses has been performed with respect to various refined theories.

KEYWORDS 3D honeycomb;anisotropic materials;differential quadrature method;general loads and constraints;higher order theories;linear static analysis;weak formulation

1 Introduction

New advances in many engineering applications are based on the employment of smart innovative materials and appliances characterized by non-conventional shapes [1,2].Very complex mechanical properties are very frequently required, together with an increased need for lightweight materials[3].As a matter of fact, when unconventional materials and geometries are adopted, the structural behaviour is not easy to be predicted a priori with simple but effective models due to the absence of symmetry planes[4,5].

In this perspective, laminated structures are widespread solutions for achieving unconventional and optimized structural performances.As it is well-known, classical composite materials usually induce an orthotropic behaviour of the overall structure[6].Nevertheless,new solutions can be found in literature where very complex deformation mechanisms can be induced, together with coupling effects.Moreover, several engineering manufacturing processes employ architectured geometries that have been conceived for the fulfillment of different aesthetic, ergonomic, and aerodynamic requirements[7].An important reason comes from the endeavor to obtain the best form under fixed load cases and external constraints layup[8].In fact,an effective reduction of the material is sought in this way.All shell structures are spatially distributed volumes [9,10].For this reason, advanced design skills are required for a correct mathematical modelling of their shape and material properties.Moreover,the accuracy of a correct modelling is crucial for its reliability since the presence of curvature significantly orients the structural behaviour,as well as the constitutive relationships[11].

As a matter of fact, three-dimensional (3D) solid simulations based on the widespread Finite Element Method (FEM) are nowadays the most adopted strategy for the structural assessment of curved and layered shells[12,13].Accordingly,such models lead to a very demanding computational cost for the simulation[14].For this reason,two-dimensional(2D)alternative approaches have been derived.The most famous strategy is Equivalent Single Layer(ESL)[15-17],which assesses the entire structural problem on a reference surface whose generic point is located in the middle thickness of the solid.In particular,in reference[16]an extensive study with higher order theories for doubly-curved shell structures is provided,and in reference[17]an ESL higher order model is developed for laminated composite curved structures.In the same way,the Layer-Wise(LW)approach[18-22]introduces a 2-manifold for each layer of the stacking sequence.In reference[19]the dynamic behaviour of anisotropic doubly-curved shells is investigated with a higher order LW approach, whereas in reference [20] an equivalent LW model is developed, accounting for an efficient description of the kinematic field,referring to the intrados and extrados of the structure.In the so-called Sub-Laminate Theory(SLT)a series of reference surfaces are defined for a specific group of the entire lamination scheme[23-25].For the LW and the SLT,the compatibility conditions between the various reference surfaces are taken into account during the global matrices assembly.Generally speaking, LW provides more accurate results than the ESL approach,but it can be very high computationally demanding since the overall number of Degrees of Freedom (DOFs) comes to be proportional to the layers within the laminate[26].Within the ESL framework, the number of unknown variables is not dependent on the actual stacking sequence.Moreover,suppose a higher order set of axiomatic assumptions is included in the model.In that case,the accuracy of ESL 2D formulations can be compared to that of the 3D FEM and LW solutions even for generally anisotropic lamination schemes[27].

As it is well known,classical finite element simulations are characterized by a local pre-determined polynomial approximation of the unknown field variable between two adjacent points.On the other hand, when shell elements with curvatures are considered, some issues like the shear locking can emerge,thus invalidating the simulation[28,29].For this reason,very refined meshes are implemented so that very accurate results are requested.In this way, the geometric and simulation discretization error is reduced,but the computational cost increases a lot.

In contrast, the IsoGeometric Approach (IGA) accounts for the employment of the actual geometry of the structure directly taken from the Computer Aided Design (CAD) process in the theoretical modelling [30-34].As a consequence, arbitrary interpolating functions are developed.Throughout literature, IGA has been employed for both the strong and weak implementations of disparate structural problems with curved edges and shapes.In particular, in reference [34] the IGA approach has been adopted for the buckling analysis of 3D plates, taking into account 2D computation of Non-Uniform Rational B-Spline (NURBS) curves, as well as the enforcement of boundary conditions.Furthermore,the hybrid meshfree collocation method is employed in reference[35],and a computationally efficient formulation is obtained,which is characterized by the advantages of IGA and meshfree methods.

When laminated doubly-curved shells are investigated by means of 2D models,a key aspect is the parametrization of the reference surfaces.It is feasible that a set of principal curvilinear coordinates are provided starting from the geometric properties of the structure.In the case of shells of arbitrary shape,a distortion algorithm is proposed so that the physical domain is described with a rectangular computational interval.To this end a set of generalized blending functions are developed[36,37].They are based on the geometric description of the shell edges and the location of the corners of the structure within the physical domain.When IGA methodology is followed, NURBS curves are employed for this task.In references[38,39],an interesting insight into the computation of such curves is provided.

Within the ESL 2D simulations,since the 3D shell structure is reduced to its equivalent reference surface as described in references [40,41] the three-dimensional unknown field variable should be expressed in terms of a set of predetermined through-the-thickness shape functions [42-45], each of them arranged employing a generalized matrix formulation [46,47].Accordingly, reference [42]employs power polynomials for the assessment of all thickness functions set and derives a closed-form analytical solution, whereas in references [43,44] a higher order power assumption is adopted only for in-plane displacement components.As far as ESL theories is concerned, both polynomial and trigonometric functions can be adopted along each shell principal direction[48-50].As the order of the kinematic expansion increases,the so-called Higher Order Shear Deformation Theories(HSDTs)come out[16].The above-mentioned generalized approach also embeds classical theories like the First Order Shear Deformation Theory (FSDT) [51,52] and the Third Order Shear Deformation Theory(TSDT) [53-55] which considers a first and a third-order polynomial expansion for the in-plane kinematic field,respectively.The out-of-plane direction is characterized by a constant distribution of the displacement field instead,thus neglecting any kind of stretching and shear effect.The employment of non-uniform through-the-thickness assumptions is crucial for modelling the actual deflection of the structure.Actually, in reference [56] it is shown that the out-of-plane stretching plays an important role in the deformation response, whereas in reference [57] the shear contribution in the total deflection of composite laminated structure is outlined.Even though an axiomatic selection of thickness functions can be assessed,some refined theories have been derived for orthotropic laminated structures in which the through-the-thickness hypotheses are based on the mechanical characteristics of the constituent materials[58-61].In particular,in reference[58]the classical FSDT theory is refined with the introduction of a shear thickness function derived from the actual lamination scheme of the selected plate.Furthermore,in reference[59]the same approach has been applied to mono-dimensional(1D) beams.In the case of laminated shells the unknown field variable dispersion should fulfil the interlaminar continuity since a piecewise transverse displacement field and shear stress distribution should be modelled.Moreover, it is feasible that the shear effect is embedded in the formulation,since experimental tests show that its contribution turns out to be prominent in laminated structures[62,63].For this reason,the LW approach is the most suitable strategy.On the other hand,the ESL approach can provide good results too if the so-called zigzag functions are adopted for the field variable description.In reference [64] an interesting review for multilayered theories is reported, where the selection of thickness functions accounts for continuous transverse stresses and a discontinuous first order derivative of the same quantities for both in-plane and out-of-plane components.In references[65,66] instead, a discontinuous first order derivative of the in-plane displacement components is provided in an effective way from the introduction of a piecewise function at the first order of the kinematic expansion.

Some considerations should be made on the treatment of external surface loads distributions within a 2D formulation.Classical solutions for plate problems [67,68] account in general for a mere computation of normal external loads within the equilibrium equations.In the case of higher order assumptions on the displacement field variable, the external loads should be inserted within the 2D approach.An effective procedure for the assessment of external loads is based on the Static Equivalence Principle[16,69].The same methodology can be applied to classical refined models when a higher order displacement field assumption is taken along each in-plane field variable.Throughout literature, a huge variety of works can be found where the static performance of plates and shells is evaluated under a uniform distribution of normal and tangential loads [70-75].In reference [75] a consistent method is proposed for concentrated loads within an ESL model directly in a strong form.In particular,a smooth continuous function is taken into account,and the calibration of the shape and position parameters can lead to a proper modelling of the singularity induced by point and line loads.On the other hand, there are only few works from literature dealing with a general distribution of surface pressure.A different group of studies investigates the problem of concentrated and line loads,starting from the numerical approximation methods of the Delta-Dirac function.Among others,the interested reader can refer to the articles by Eftekhari et al.[76-78],where the Dirac-Delta function is discretized with success taking into account the employed numerical method.In particular,in reference[76] moving concentrated loads have been applied to 1D structures starting from the definition of the Dirac-Delta function.Then,in reference[77]the formulation has been extended to 2D structural problems.Moving load problems have been treated in reference [78] with the well-known Heaviside function.Generally speaking,the approximating function of a singularity problem should fulfill some characteristic properties derived from the theoretical definition of concentrated loads; otherwise, it should be properly normalized.These kinds of pressures are singularities in a differential continuum model, such that various numerical techniques have been proposed for their approximation with smooth functions [79-81].According to the author’s knowledge, there are some works accounting for a Gaussian distribution (see reference [82] among others).An interesting problem related to concentrated loads on curved structures is the so-called pinched cylinder.In references[83-85]a series of theoretical developments have been provided for a circular cylinder subjected to a point load applied at the diameter of the surface.

As far as the numerical implementation is concerned, several articles can be found in which spectral collocation methods are adopted within IGA framework.Among them, the Generalized Differential Quadrature (GDQ) method [86-90] has demonstrated to be a reliable tool since it overcomes many computational drawbacks.Based on polynomial functional approximation methods,the GDQ has been successfully applied to a series of structural problems[91,92].In particular,some works can be found where GDQ has been adopted for the dynamic analysis of smart structures,as well as multifield problems [93,94].It has been shown that it is an extension of pseudospectral collocation methods since the similar levels of accuracy are obtained when the same computational grid is adopted[95-97].On the other hand,the generalized procedure for the calculation of weighting coefficients based on Lagrange polynomials for the global interpolation is not based on the actual selection of discrete points and the derivation order.In particular,the GDQ method provides the best performance among all the numerical techniques belonging to the weighted residual class.In reference[98]the outcomes of the free vibration analysis of thin-walled beams calculated by means of the GDQ method are compared to those coming from a 1D Ritz formulation.Other works[99-102]provide the theoretical assessment together with some parametric static and dynamic analyses of shells of different curvatures characterized by innovative material like Carbon Nanotubes (CNTs) and Functionally Graded Materials(FGMs).In reference[103]FGM structures have been investigated by employing a semi-analytical model.In reference[104]an effective homogenization algorithm is developed for the mechanical assessment of sandwich shells with both honeycomb and re-entrant patterns.Referring to pantographic lattice structures, the widespread method for the computation of equivalent elastic properties is the well-known Neumann’s principle[105],accounting for the computation of generalized stiffnesses starting from simple independent load cases[106-108].The accuracy of the solution lies in the theoretical hypothesis considered for the numerical modelling.In particular,some considerations should be made on the nodal area within the honeycomb cell,as well as the bending and shear effects of vertical and inclined struts[109-111].

Once the 2D solution is provided on the reference surface of a doubly-curved shell, the postprocessing algorithm is crucial for the reconstruction of both static and kinematic quantities along the three-dimensional solid.Actually,both the ESL and the LW formulation account for the assessment of the generalized unknown field variables lying on the reference surface.Accordingly, in the case of simulations performed with the FEM, the solution is calculated at the discrete set of nodes of the computational grid.Then, a reconstruction along the continuum domain should be assessed.Referring to the latter, there are three main classes of methodologies for the solution extrapolation,namely the Superconvergent Patch Recovery(SPR)[112,113],the Averaging Method(AVG)[114,115]and the Projection Method (PROJ) [116].As far as the through-the thickness equilibrium-based recovery procedure is concerned, the generalized constitutive relationship can be adopted as usual,but a correction of shear and membrane stresses is essential so the boundary conditions related to applied loads can be fulfilled[117-120].

In the present work,an ESL theoretical formulation is proposed for generally anisotropic doublycurved shells of arbitrary shapes by means of the ESL method.The reference surface of the structure is described employing the main outcomes of differential geometry,setting an orthogonal curvilinear set of principal coordinates.A NURBS based generalized set of blending functions is employed for the distortion of the physical domain in the case of arbitrarily-shaped structures.The displacement field component vector is intended to be the unknown variable of the problem,and it is expressed via the employment of a generalized set of thickness functions.Actually,a higher order kinematic expansion is performed together with a consistent zigzag function to properly catch all the possible stretching and warping effects within each layer, as well as coupling issues at the interlaminar stage.A generalized weak expression of the field variable is introduced,accounting for higher order shape functions for the interpolation alongside the reference surface [121].A general distribution of surface loads has been considered in the model, setting a Gaussian and a Super-Elliptic distribution for the application of external loads along each principal direction of the shell.In particular,the calibration of distribution parameters leads to model different load cases.The fundamental governing equations,derived from the Minimum Potential Energy Principle,are numerically tackled by the GDQ method.The Generalized Integral Quadrature(GIQ)method has been adopted for the computation of integrals[16].The GIQ moves from the above-discussed GDQ procedure, and it allows performing numerical integrations following an effective quadrature approach.

The theoretical formulation proposed in this manuscript has been implemented in the DiQuMASPAB software [122], a research code based on the GDQ method accounting for static and dynamic simulations on doubly-curved shell structures.All the numerical examples of mechanical properties have been included in the built-in material library.The graphic user interface provided by the software allows to implement the lamination scheme, the geometry of the structure, as well as external load and boundary conditions features.The structures that are considered for the numerical validation account for different mechanical issues, namely the presence of a single and double curvature, the number of layers and the presence of the soft core.Different loading conditions have been adopted.The equilibrium-based 3D reconstruction of mechanical quantities has been demonstrated to be very efficient in accordance with the predictions of refined 2D and 3D solutions.The same level of accuracy has been outlined in the case of lamination schemes with a lattice layer characterized by a macro mechanical orthotropic behaviour.Any performance reduction has been observed when a mapping of the physical domain has been required.

2 Geometrical Description of a Doubly-Curved Shell

According to the ESL formulation,a doubly-curved shell should be referred to its reference surface whose points are located at the middle thickness of the structure.It should be remarked that a generic three-dimensional solid can be expressed with respect to a Cartesian global reference systemO x1x2x3of unit vectors e1,e2,e3in terms of three parametersα1,α2,α3:

being R(α1,α2,α3)the three-dimensional position vector,andfi(α1,α2,α3)fori=1,2,3 a continuous function.If the variablesα1,α2,α3are taken as a set of curvilinear coordinates,it is possible to introduce the above discussed reference surface r(α1,α2)starting from Eq.(1)as follows[16]:

where theα3=ζhas been selected alongside the normal direction of r.A schematic representation of the ESL assessment of the shell geometry has been reported in Fig.1.The normal unit vector can be computed if the partial derivatives of r(α1,α2)with respect toαi=α1,α2,namelyare known:

More specifically,Eq.(2)describes a doubly-curved shell when all the variables vary in a closed interval.In particular, it should be stated thatAccordingly,shell thicknessh(α1,α2)can be defined starting from Eq.(2),according to the following expression taken from reference[16]:

wherehk(α1,α2)=ζk+1(α1,α2)-ζk(α1,α2)is the width of thek-th lamina of the stacking sequence composed ofl∈N layers located within the intervalIf the second order derivativesof the reference surface withαi,αj=α1,α2are calculated,it is possible to compute the principal radii of curvatureRi=R1,R2in each point of the physical domain:

In addition,the well-known Lamè parametersAi=A1,A2are defined as follows:

Figure 1: Graphic representation of the ESL methodology for the structural assessment of doublycurved shells of arbitrary shape.NURBS mapping of the structure and definition of a 2D rectangular computational domain

Starting from Eqs.(1) and (2), the parametersHi(α1,α2,ζ)are defined, accounting for the curvature effect alongside the thickness directionζ:

Accordingly, in the case of straight shells along theαiprincipal direction, Eq.(7) turns intoHi(α1,α2,ζ)=1.

In the present manuscript,doubly-curved shells of variable thickness[16]have been considered;therefore,a generalized analytical expression ofh=h(α1,α2)is provided hereafter:

beingh0kthe reference height of thek-th layer andδjforj= 1,...,4 a scaling parameter, whereasis a constant dimensionless shift number.In Eq.(8),thickness variation(α1,α2)has been described in terms of four analytical expressionsφj(α1,α2)employing a dimensionless coordinate defined along eachαi=α1,α2principal direction of the shell,according to the following equation:

In the present formulation,the expressions ofφj(α1,α2)forj= 1,...,4 have been implemented,beingαjm∈[0,1],nj∈N andpj∈R[16]:

When an arbitrary curve is considered lying on the reference surface r(α1,α2)described with principal coordinates,a local reference system composed of three unit vectors nn,nsand nζshould be defined in each point of the curve at issue.More specifically,nsaccounts for the tangential unit vector of the curve, whereas nnis the normal unit vector of the curve.nζdenotes the bi-normal direction.Referring toO＇α1α2ζmiddle surface reference system,nn,nsand nζcomponents can be computed in each point of the curve,leading to:

Since the curve at issue lies on r(α1,α2), it should be noted thatnn3=ns3=nζ1=nζ2= 0 andnζ3=1.

3 NURBS Blending of the Physical Domain

In the previous paragraph, the geometry of a generic doubly-curved shell has been described starting from the reference surface r(α1,α2), according to the ESL framework.Nevertheless, it is possible to perform a distortion of the original physical domainin order to take into account structures of arbitrary shape.To this purpose, a set of blending coordinatesξ1,ξ2∈[-1,1]are introduced along the parametric lines of the mapped shell, as it has been shown in Fig.1.In particular, the four edges of the physical are identified on the physical domain according to the following nomenclature[16]:

where I is the identity matrix.On the other hand,vectors b and matrix B read as follows:

Starting from Eq.(13),it is possible to derive an expression for partial derivatives with respect to in-plane directionsα1,α2in terms of natural coordinatesξ1,ξ2thanks to the well-known chain rule[16]:

If det(J)/=0,the inverse form J-1of the Jacobian matrix occurring in Eq.(17)can be calculated,thus obtaining:

Accordingly,the inverse relation of Eq.(17)accounts as follows:

From a direct comparison between Eqs.(16) and (19), the definition of coefficientsξi,αjcan be assessed[16]:

For the sake of completeness,first order partial derivatives detwithi= 1,2 of det(J)with respect to natural coordinatesξ1,ξ2can be expressed as:

Following a similar procedure of that adopted in Eqs.(16)-(20), the chain rule can also be employed for the computation of second order derivatives with respect toα1,α2principal coordinates in terms ofThe final expression for second order mixed derivatives looks as follows[16]:

beingξj,α1α2withj=1,2 a transformation coefficient reading as:

In addition,one gets the following expression for pure second order derivatives with respect toαidirection,settingi=1,2[16]:

where the coefficientsξj,αiαifori,j=1,2 complete expression are thus reported:

In the present work,theβ-th edgeof the mapped physical domain forβ= 1,...,4,whose componentsalongαidirection have been collected in the vectorfori=1,2,has been geometrically parametrized with respect toξjforj=1,2 in terms of NURBS curves[16].

4 ESL Kinematic Equations Employing Generalized Shape Functions

We now deal with the definition of the ESL assessment of the displacement field variable employing a generalized set of shape function for the weak formulation of the structural problem.The 3D displacement field component vector U(α1,α2,ζ)= [U1U2U3]Tassociated with each point of the solid of the Euclidean space,defined in Eq.(2)with respect to the principal set of coordinatesα1,α2,ζ,can be expressed according to the following equation(Fig.1):

Employing the unified notation of Eq.(26), a generalized displacement field component vectoris assessed for eachτ-th order of the kinematic expansion,defined in each point of the reference surface r(α1,α2).The dependence of U from the out-of-plane coordinateζis taken from the definition,for eachτ=0,...,N+1,of a series of thickness functionsalong eachαi=α1,α2,α3principal direction.It should be noted that the unknown field variable arrangement introduced in Eq.(26)comes into the definition of a 2D model for an arbitrary shell.Moreover,ifare properly selected,a wide range of effects can be effectively predicted,thus awarding the model of 3D capabilities.In the present work,have been chosen defined from power polynomials[16]:

Referring toτ=N+1,a dimensionless coordinatezkhas been defined in eachk-th layer of the stacking sequence starting from the through-the-thickness coordinate of the top(ζk+1)and the bottom(ζk)surface delimiting thek-th sheet,respectively.Ifk= 1,...,lwithldenoting the total number of laminae,one gets:

As a matter of fact, the employment of the unified notation of Eq.(26) lets us to perform a systematic analysis with different ESL theories, since the theoretical formulation is independent from the actual selection of the thickness function analytical expression.Nevertheless,since various polynomials orders can be selected in Eq.(27),a nomenclature is adopted for a smarter identification of the axiomatic through-the-thickness assumptions of the simulation.In particular, the acronym ED(Z)-Nis adopted.Letter “E”tells that an ESL formulation is considered.“D”stands for the displacement field as an unknown variable,whereasNdenotes the maximum order of the kinematic expansion.When the additional generalized displacement field component vector u(N+1)(α1,α2)is introduced according to what exerted in Eq.(27),the letter“Z”is included.In this way,the formulation accounts for the zigzag function in eachαi=α1,α2,α3component[16].

Once the ESL assessment of the displacement field component vector U(α1,α2,ζ)=[U1U2U3]Thas been performed,Eq.(26), defined in each point of the physical domain, should be referred to a 2D computational grid composed of a generic discrete set ofIN×IMpoints, whose generic element can be identified with（α1f,α2g）forf=1,...,INandg=1,...,IM.For this purpose,an interpolation algorithm is introduced for u(τ)(α1,α2)at eachτ= 0,...,N+1 kinematic expansion order.A very useful tool is the vectorization operator [16], accounting for a rearrangement with a 1D vector of a 2D matrix.If we denote with A a matrix of dimensionsIN×IMof componentsAijwithi= 1,...,INandj= 1,...,IM,a column vector= Vec(A)of dimensionsINIM×1 is defined,whose arbitrary componentAkfork=1,...,INIMreads as:

Since all quantities in Eq.(26) are evaluated in each point of the discrete computational grid,quantities u(τ)iare introduced withi= 1,2,3 for eachτ-th kinematic expansion order.Taking into account Eq.(29),they are arranged as follows:

Once column vectors of Eq.(30)have been correctly defined the following definition can be stated,accounting for all the values assumed by the generalized displacement field component vector in each point of the computational grid:

Thus, it is possible to introduce the discrete form of Eq.(26) employing a global interpolation algorithm.Accordingly, the following relation can be assessed for eachτ-th order generalized displacement field component vector u(τ)employing a matrix NTof generalized shape functions:

In a more expanded form,one gets[16]:

Starting from Eq.(33) it is possible to define the generalized shape function matrix NT(α1,α2)alongside the physical domain, whose dimensions 3×3INIMdepend on the selected computational grid.Eventually,one gets:

In Eq.(35),the Lagrange interpolating polynomials introduced in Eq.(34)are collected in a 1×IPvector withIP=IN,IMdenoted withlαifor eachαi=α1,α2,according to the following expanded form:

In the same way, vectoraccounting for the first order derivativeof the Lagrange interpolating polynomials collected according to Eq.(36) with respect toαi=α1,α2is defined as follows,settingη=f,g:

Since the Kronecker tensorial product,denoted with ⊗,has been adopted in Eq.(35),it is useful to introduce arow vector of dimensions 1×INIMso that NTmatrix can be expressed in a compact form as follows[16]:

To sum up, the ESL assessment of the kinematic field introduced in Eq.(26) is rearranged accounting for the interpolation procedure of Eq.(32),leading to the following expression:

Now, the generalized form of the displacement field of Eq.(39) is adopted for the definition of the kinematic relations for a doubly-curved shell structure, according to the ESL approach.More specifically, the kinematic relations for a 3D structure in principal coordinatesαi=α1,α2,α3(withα3=ζ)read as[16]:

beingε(α1,α2,α3)= [ε1ε2γ12γ13γ23ε3]Tthe 3D strain component vector and D a differential operator relatingεto the displacement field component vector U(α1,α2,α3).Accordingly, the Dζoperator accounts for the partial derivatives with respect to the shell principal directionα3=ζ, as well as the curvature thickness parametersH1,H2:

Introducing in Eq.(40) the ESL assessment of the displacement field settled in Eq.(26), the kinematic relations turn into the following equation,accounting for an expansion up to the(N+1)-th order[16]:

In the previous equation,the strain vectorε(α1,α2,α3)referred to a 3D solid has been expressed in terms of the ESL generalized strain component vector,defined for an arbitrary order of the kinematic expansion.As can be seen,Z(τ)αiis introduced for eachαi=α1,α2,α3and forτ=0,...,N+1 according to the following extended form[16]:

In Eq.(44)it has also been shown that the generalized strain component vectorε(τ)αiembeds the ESL assessment of the displacement field of Eq.(26).Based on the generalized interpolation procedure performed in Eq.(32),ε(τ)αican be related to the discrete vectorof the generalized displacement field of Eq.(31)for eachτ=0,...,N+1,i.e.,

where the kinematic operators Bαiwithαi=α1,α2,α3are defined from the previously discussed higher order interpolation algorithm performed on theIN×IMdiscrete grid by means of the tensorial Kronecker product ⊗[16]:

We recall that in Eq.(47) the Lagrange polynomials employed for the interpolation of the field variable along theαi=α1,α2in plane direction, as well as their first order derivative, have been collected in the vectorslα iand,respectively.

5 Generally Anisotropic Lamination Scheme Assessment

The present manuscript investigates the problem of t he static structural response of a doublycurved shell with generally anisotropic laminated structures characterized by general orientations and material syngonies.It should be remarked for the sake of completeness that a sequence ofllaminae has been contemplated, whosek-th layer withk= 1,...,lis characterized by a thicknesshkaccording to the conventions of Eq.(8) and a reference systemorientated along the material symmetry axes.Accordingly, the formulation has been developed by assuming=ζ.In addition, a characteristic angleθkdescribing the orientation ofaxis with respect toα1is introduced.Each layer of the stacking sequence is intended to be generally anisotropic.The mechanical behaviour is described with respect t o the 3D stress and strain component vectorsfrom the introduction of a generally anisotropic stiffness matrix E(k)for eachk-th layer.The generalized constitutive relationship can be thus defined as[16]:

Since a linear elastic theoretical assessment of the structural problem is performed,the matrix E(k)reads as follows:

wherewithi,j=1,...,6 relates thei-th component ofvector to thej-th strain component of thevector.Accordingly,the coefficients employed in Eq.(49)account for the 3D elastic stiffnesses of the material, namelyWhen E(k)is referred to the plane stress assumptionthe reduced stiffness coefficientscan be employed, defined from the expression reported hereafter:

Nevertheless,Eq.(48)should be referred to the curvilinear principal reference system of the physical domain,accounting for the 3D stress and strain vectorsandTo this end, a generalized rotation matrix T(k)is defined according to reference [30] within each layer starting from the previously discussed angleθkemployed to compute the rotated stiffness matrix= T(k)E(k)（T(k)）Twith componentsfori,j= 1,...,6 referred to the shell geometric reference systemO＇α1α2ζ.More in detail,Eq.(48)takes the following form:

On the other hand Eq.(51),defined for thek-th layer of the 3D solid,should be assessed within the ESL framework accounting for all thellaminae occurring in the lamination scheme.To this purpose,the actual dispersion of stress components are employed for the definition of the generalized stress resultant component vectorassociated to eachτ-th order of the kinematic expansion of Eq.(39).In this way, Eq.(51) can be expressed in terms of S(τ)αiandε(η)αivectors, according to the following equations directly derived from the static equivalence principle:

where A(τη)αiαjis the generalized ESL stiffness matrix[16]:

whose generic componentis defined as:

Employing the kinematic relations performed in Eq.(46)utilizing the generalized interpolation,Eq.(52)gets into:

Thus,the generalized ESL assessment of the constitutive relationship performed in Eq.(56)can be expressed so that each component of S(τ)αi(α1,α2)can be expressed in terms of the generalized displacement field component vectorforη= 0,...,N+1 at each point of the 2D discrete grid alongside the physical domain.Eventually,one gets:

In Appendix A, the interested reader can find the complete expressions for each component of matrix O(τη)αiαj.

6 Equilibrium Relations in the Variational Form

In the present section,the equilibrium equations for an arbitrary doubly-curved shell are derived by employing an energy approach.The stationary configuration is specifically enforced to the total potential energy virtual variationδΠ.The contribution of external loads should be inserted within the ESL framework, accounting for various pre-determined distributions along the topand the bottomsurfaces of the shell.Moreover,the formulation embeds the contribution of linear elastic springs distributed along the shell edges.Suppose we denote withLes,Lebthe virtual work of surface loads and boundary springs,respectively,the equilibrium configuration satisfies the following equation,accounting for the stationary assessment of the virtual variationδof the computed quantities[16]:

beingΦthe total elastic strain energy of the 3D solid, expressed via the ESL definition of the displacement field,as reported in Eq.(39).

6.1 Elastic Strain Energy

We focus on the total elastic strain energyΦof the 3D shell described according to the ESL approach of Eq.(2).The virtual variation ofΦcan be computed in terms of the 3D stress and strain vectorsσ(k)andε(k)defined in eachk-th lamina of the stacking sequence concerning the geometrical principal reference systemO＇α1α2ζ,according to the following relation:

Employing the ESL assessment of the displacement field of Eq.(26), the through-the-thickness integration performed in Eq.(59) is avoided.Thus,δΦcan be computed in terms of S(τ)αiandε(τ)αi,leading to[16]:

Taking into account the weak formulation expression of generalized strain resultant vector of Eq.(46)in the previous equation,the virtual variation of the total elastic strain energy reads as:

In Eq.(61)the global stiffness matrixof the doubly-curved shell is introduced.Referring to the genericτ,η-th kinematic orders,takes the following expanded form[16]:

In Appendix B, an extended expression of each stiffness matrix coefficientis reported,settingτ,η=0,...,N+1 andαi,αj=α1,α2,α3.

6.2 External Loads General Distributions

Let us consider two generic vectorsandof external loads applied to the doubly-curved shell structure at issue,whose components are referred in the above discussed geometric reference systemO＇α1α2ζas shown in Fig.2:

where symbols(+)and(-)account for quantities referred torespectively.Accordingly,a general variation1,α2)dependent on in-plane principal coordinates can be assigned to each componentfori=1,2,3 andj=1,2 of external load vectors introduced in Eq.(63).If we denote witha reference scaling value of loads applied along theαiprincipal direction,the following relation can be assessed(Fig.2):

Figure 2:Generalized 2D distributions(α1,α2)along shell curvilinear principal coordinates α1,α2

For a constant loading distribution,it is:

In the present manuscript, two different distributions forcomponents withj= 1,2 are declared.A bivariate Gaussian distribution is now reported:

beingαimfori=1,2 the position parameter acting on theαiin-plane principal direction with a variance scaling factor equal toΔiand a correlation factor betweenα1,α2denoted withρ12.In the same way,a Super Elliptic distribution has been employed,according to the following expression:

whereβis the angular orientation of the bivariate dispersion principal axes of lengthΔifori= 1,2 with respect toα1,α2principal shell coordinates.In addition,(α1m,α2m)is the location in the physical domain of the centre of the ellipse,whereasnis a power exponent.Forn= 2,Eq.(67)turns into the well-known Elliptic distribution.In Fig.2 we plot the distributions presented in Eqs.(66)and(67),for fixed parameters configurations.

If a virtual variationδU(j)of the 3D displacement vector referred to the top(j=1)and the bottom(j=2)surface is considered,the virtual workof surface loads of Eq.(63)can be computed as:

beingH(j)iforj= 1,2 the through-the-thickness curvature parameters calculated by means of Eq.(7)forOn the other hand, if we consider a virtual variationδu(τ)of the generalized displacement field component vector introduced in Eq.(26),a new generalized vector of external actionslying on the reference surface can be introduced for eachτ=0,...,N+1.Thus,the virtual workδLesreads as follows:

The static equivalence principle employs both Eqs.(68)and(69),so that[16]:

Starting from Eq.(70),the extended expressions of the q(τ)components are obtained for theτ-th kinematic expansion order,beingforj=1,2 the thickness function employed for the assessment of the displacement field components according to Eq.(26) calculated atrespectively:

Since the generalized displacement field component vectors u(τ)withτ= 0,...,N+1 have been interpolated with a generalized shape functions set N in Eq.(39),the virtual workδLesof generalized external actions already computed in Eq.(69)should be expressed in terms ofvector taking into account Eq.(32),leading to the following expression[16]:

In Eq.(72),the generalized external actions referred to eachτ-th order have been collected in the 3INIM×1 column vector Q(τ),defined as:

6.3 Boundary Springs

The last contribution to the total potential energy occurring in Eq.(58)accounts for the influence of the elastic constraints distributed along the edges of the structure.In particular, a set of linear elastic springs of stiffnessorientated along the shell side located atαj=αm jwithαj=α1,α2,α3andm= 0,1 distributed along theαi=α1,α2,α3principal direction has been considered,according to the theoretical development reported in reference [16].The superscriptkstands for the stacking sequence layer.Accordingly,they induce some boundary stresses after applying the corresponding 3D displacement componentUi=U1,U2,U3.If the physical domain consists of a rectangular 2D intervaldifferent stresses are induced on each boundary.Referring toα1=αm1form=0,1,the following definitions can be stated at thek-th lamina level,setting

As a matter of fact, the Super Elliptic distribution (S) of thep-th order has been employed,according to the following analytical expression, setting∈[0,1] and∈]0,1] the position and the scaling parameter,respectively:

In the following,the expression of a Double-Weibull distribution(D)has been reported in terms of:

Starting from Eqs.(74) and (75) which have been expressed for a 3D structure, the stresses coming from the linear elastic springs distributions should be computed in terms of generalized stress resultants,accounting for the ESL assessment of the displacement field.Substituting Eq.(26)in the previously cited relations,one gets forα1=α1mwithm=0,1:

Referring to the edges of the structure located atα2=α2mform=0,1,Eq.(75)turns into:

Fundamental coefficientsoccurring in Eqs.(79)and(80)can be computed for eachτ,η=0,...,N+1 according to the following relation,settingn,p=1,2,i=1,2,3 andn,p=1,2:

Having in mind the ESL assessment of general distributions of boundary springs,the virtual work contributionδLebto Eq.(58)can be computed as the sum ofδLeb1andδLeb2,referred to the edges of the physical domain characterized byα1=α1mandα2=α2m,respectively,form=0,1[16]:

Thus,the generalized displacement vectorof Eq.(32)can be outlined in Eq.(83),accounting for the definitions of Eq.(84), together with what has been exerted in Eqs.(79) and (80).The final form of the virtual work associated with boundary springs reads as:

Note that the fundamental stiffness matrixhas been introduced, accounting for the ESL generalized stresses induced by the springs distributed along the edges of the physical domain,at eachτ,η=0,...,N+1 kinematic expansion order,as well as the generalized shape functions arranged in N,based on the Lagrange interpolation procedure of Eq.(32).

6.4 Fundamental Governing Equations

The final form of the fundamental governing equations of the static problem can be outlined from Eq.(58),taking into account the contributions of the elastic strain energy of Eq.(61),the virtual work of surface loads of Eq.(72)and that of external boundary springs of Eq.(85).Referring to theτ-th order of the kinematic expansion,one gets[16]:

In extended form, Eq.(86) can be assembled with respect to the displacement field ESL higher order model of Eq.(26),leading to the following expression:

For shells of arbitrary shapes,the employment of the generalized blending functions assessed in Eq.(13) with the Jacobian matrix J of the transformation defined in Eq.(17) should be taken into account in the computation of all the contributions ofΠoccurring in Eq.(58).Accordingly, the generalized stiffness matrixand the external surface loads Q(τ)should be adjusted appropriately for eachτ=0,...,N+1 since integrals with respect to principal directionsα1,α2should be assessed in terms of natural coordinatesξ1,ξ2∈[-1,1].Referring to a specificτth kinematic expansion order,stiffness matrix defined in Eq.(62)turns into the following relation:

The ESL assessment of external surface load vector Q(τ)reported in Eq.(73)reads as follows[16]:

As far as the definition of generalized boundary springs for arbitrarily-shaped structures is concerned,it is useful to refer to the local coordinate system introduced in Eq.(11)along each edge of the shell.As a consequence,boundary stresses can be expressed with respect to nn,ns,nζdirections employing the symbolsThe generalized boundary springs along thei-th edge of the structure withi=1,...,4 accounts as:

Starting from the blending coordinates transformation of Eq.(13),the lengthL(i)of the shell side can be computed according to the following relation[16],settingj=1,2:

Since the natural coordinate setξ1,ξ2has been defined in the dimensionless rectangular interval[-1,1]×[-1,1],Eq.(90)can be written as:

In order to compute the generalized stress resultantsinduced by linear elastic springs for a mapped geometry employing the approach of Eqs.(74)and(75),it is useful to express the generalized displacement fieldu(nτ),u(sτ),u(ζτ)referred to nn,ns,nζlocal reference system in terms ofu(1τ),u(2τ),u(3τ)for eachτ-th order of the kinematic expansion[16]:

In the same way,the generalized stress resultantsshould be expressed in terms of the S(τ)αivector forαi=α1,α2,α3according to the following relations[16]:

7 Numerical Implementation with the GDQ Method

In the present section we deal with the numerical assessment of the fundamental governing equations for the static problem derived in Eq.(87)in their assembled form.The employed numerical technique that has been adopted is the well-known GDQ method.According to this methodology,the discrete form of the derivatives is directly provided.Referring to a generic point（ξ1,ξ2）∈[-1,1]×[-1,1]belonging to the dimensionless squared computational domain,the(n+m)-th order derivative of a bivariate functionconcerning natural coordinatesξ1,ξ2can be expressed as[16]:

whereIN,IMdenote the number of the selected discrete points alongξ1,ξ2directions,respectively.As can be seen,Eq.(95)accounts for a quadrature rule for derivative purposes.The weighting coefficientsforr= 1,2,p=i,j,q=n,mandv=k,lare computed utilizing the recursive expression valid for eachq≥1[16]:

Eq.(96)is based on the interpolation of the unknown function employing the Lagrange interpolating polynomialsL.In the case of pure derivatives with respect toξ1,ξ2,it givesm= 0 andn= 0,respectively,and GDQ coefficients come into the following relation:

beingδpvthe Kronecker delta operator.For a generic non-uniform grid ofIN×IMpoints,the GDQ weighting coefficients for the derivative of the(n+m)-th order are calculated by means of Eqs.(96)and(97)are collected in aINIM×INIMsquared matrix denoted withςξ1(n)ξ2(m).Starting from Eq.(96),the matrix at issue can be obtained from a proper expansion of matrixςξr(q)forq=n,mandr= 1,2 of dimensionsIQ×IQwithIQ=IN,IM, respectively, accounting for the derivatives ofn-th andm-th order alongξ1,ξ2:

Based on Eq.(97),it isςξr(0)=I.In this way,the numerical derivation employing the GDQ method with respect toξ1,ξ2can be performed as:

where f is aIN×IMmatrix whose general component isf（ξ1i,ξ1j）fori= 1,...,INandj= 1,...,IM.In the same way,the generic element of f(n+m)is

A 2D Legendre-Gauss-Lobatto(LGL)grid of dimensionsIN×IMreferred to the domain[-1,1]×[-1,1]has been here adopted,such that:

In Eq.(100),AIP-1(ξr)withr=1,2 denotes the Lobatto polynomials of theIP-th order,which can be calculated from theIP-th derivation of the Legendre polynomialsLIP-1(ξr)with respect toξr[16]:

Since the fundamental governing relations (86) account for derivations alongα1,α2principal directions, we consider the GDQ rule of Eq.(95) referred to natural coordinatesξ1,ξ2.The matrix approach is then followed, starting from Eq.(98).Employing the generalized blending functions reported in Eq.(13), the matricesςα1(1)α2(0)andςα1(0)α2(1)for the first order derivatives alongα1,α2of a generic bivariate functionf (α1,α2)can be computed as:

The employment of the GDQ method allows to obtain the discrete form of the governing equations.When the assembly of the fundamental relations is performed alongside the entireIN×IMcomputational domain,it is useful to distinguish the inner domain nodes“d”from those located along the boundaries(“b”nodes).As a matter of fact,the DOFs associated with the latter are collected in a column vector denoted withδb,whereas those referring to the former are arranged within the vectorδd.Thus,the fundamental system accounts as follows[16]:

Performing a static condensation of Eq.(104) with respect toδdaccording to the procedure reported in Eq.(16),the final form of the discrete system of dimensions 3(N+2)INIM×3(N+2)INIMcomes out,being I the identity matrix:

The numerical integrations that are involved in the formulation of the present manuscript are performed by means of the GIQ procedure[16].Referring to a generic univariate functionf=f (x)defined in the closed interval [a,b] witha ＜b∈R, a discrete grid is selected from the 1D domain,namelyxk∈[a,b] withk= 1,2,...,IN.Accordingly, the integral offrestricted to the sub-intervalwithxi,xjbelonging to the grid at issue can be computed as a weighted sum of the valuesf (xk)fork= 1,2,...,INassumed by the function at each point of the computational grid of[a,b]:

where=wjk-wikare referred to the intervalThe coefficientswjk,wikare computed following the extensive procedure of reference [16], starting from the shifted GDQ weighting coefficients reported hereafter,settingε=1×10-10:

Then, they are collected in the matrixof dimensionsIN×IN.Thus, the coefficientswik,wjkreferred to the integral offrestricted to the interval[ε,xr]⊆[a,b]are the elements of the GIQ matrix,reading as:

All the termsintroduced in Eq.(106) can be collected in a vector of dimensions 1 ×INdenoted with w(x), settingxthe integration variable.In the same way, the 1×IMvector w(y)can be introduced when the integration of the univariate functionf=f (y)is required.For a bivariate functionf=f (x,y),Eq.(106)turns into the following expression,setting a generic 2DIN×IMgrid from a rectangular domain[a,b]×[c,d]of extremesa ＜b∈R andc ＜d∈R[16]:

wherefkvalues fork= 1,...,INM=INIMare arranged employing the above introduced by-column vectorization operation of Eq.(29).On the other hand,GIQ coefficients are computed using the previously introduced vectors w(x),w(y),as reported:

Referring to the bivariate functionf=f (x,y), it is also possible performing a univariate integration along a parametric line of the rectangular computational domain parallel toxandy,respectively, employing the GIQ rule.For this purpose, w(x0)and w(0y)are defined starting from Eq.(110):

For the sake of completeness, symbol o(x),o(y)are row vectors of ones of dimensions 1×INand 1×IM,respectively.Thus,the following relations can be derived:

We recall that in Eq.(87)the global stiffness matrix of the shell has been provided,starting from a parametrization of the reference surface r(α1,α2)employing principal coordinates.As it has been shown,for arbitrarily-shaped structuresforτ,η=0,...,N+1 is described by means of natural coordinatesξ1,ξ2in Eq.(88).For this purpose,it is useful to define a matrix of dimensions 3INIM×3I2NI2Mfor the transformation at issue:

being O a ones squared matrix of dimensions 3INIM.Accordingly, w(ξ1ξ2)=w(ξ2)⊗w(ξ1)following Eq.(110), whereas DJ, A1, and A2accounts for the determinant of the Jacobian matrix J, as well as the Lamè parametersA1,A2of the shell.If the indexk=n+(m-1)INwithn= 1,...,INandm=1,...,IMis assessed,we obtain the vectorized form of,accounting for theIN×IMgrid of the computational domain:

Employing Eqs.(113)and(114),the stiffness matrixfor an arbitrarily-shaped shell structure can be computed:

In the following, the discrete form of the generalized external loads of Eq.(84) is provided by means of the GIQ method presented in Eq.(106):

In the same way,Eq.(116)becomes suitable for an arbitrary mapped domain.According to the nomenclature introduced in Eq.(12),it gives forξ2=-1:

The edge characterized byξ1=1 behaves as:

On the other hand,forξ2=1 it is:

Finally,the edge identified withξ1=-1 reads as:

8 Post-Processing Analysis

Once the fundamental governing relations reported in Eq.(87) have been solved by means of the GDQ method, the solution in terms of the generalized displacement field component vector u(τ)is obtained.Actually, the 3D distributions of both kinematic and mechanical quantities along the shell thickness should be recovered starting from the outcomes of the 2D model.To this purpose, the Chebyshev-Gauss-Lobatto (CGL) non-uniform grid is introduced.The locationξgof the computational points is provided for the[-1,1]dimensionless interval[16]:

beingITthe number of the selected discrete points.Referring to the interval[ζk,ζk+1]fork= 1,...,lrepresenting thek-th layer of the staking sequence within the 3D physical domain,Eq.(121)should be transformed according to the following relation:

Employing the unified ESL assessment of the displacement field introduced in Eq.(26) it is possible to derive the through-the-thickness distribution of the displacement field component vectorwithi=1,...,IN,j=1,...,IMandg=1,...,ITfor eachk-th layer of the stacking sequence[16]:

Referring to the generic（α1i,α2j,ζg(k)）point of the shell structure, membrane stresses can be calculated according to the generally anisotropic law of Eq.(51), leading to the following relation,settingfori,j=1,...,6 the elastic stiffness coefficients of thek-th layer referred to the geometric reference system employed for the computation of,,stresses:

Nevertheless, the remaining stress components,,are not calculated starting from the constitutive relation (51).The 3D equilibrium equations [16] expressed in curvilinear principal coordinates are employed,instead,as:

whereR1,R2are the principal curvature radii of the reference surface r(α1,α2), whereas the termsa(k),b(k),c(k)depend on the quantities calculated in Eq.(125),reading as:

It should be noted that the expression ofc(k)contains the terms,.For this reason, the first two relations of Eq.(126) should be solved independently in the first step.Then, the solution of the third equation can be easily found.To this end, the discrete form of the derivatives of inplane stresses,,with respect toα1,α2occurring in the first two expressions of Eq.(127)should be performed employing the GDQ method for eachk-th layer of the lamination scheme,settingk=1,...,l:

The first order differential relations(126)should fulfill a single boundary condition.Referring to the first lamina of the stacking sequence(k=1),the external loads applied at the bottom side of the structuresatisfy the following boundary conditions:

For the surface loadsandat the top side,an adjustment algorithm is provided starting from the solution of Eq.(126)by means of the boundary conditions(129)and(130).Accordingly,the exact fulfilment of the compatibility conditions between stresses and in-plane external loads is ensured at the external shell skins[16]:

In the previous equation an indexs=kghas been introduced,accounting forg=1,...,ITpoints defined according to Eq.(122)for thek-th lamina alongside the entire shell thicknessfori= 1,...,INandj= 1,...,IM.As a matter of fact,stressescome from the direct solution of the balance Eq.(126).The introduction of the indexsleads to vectors,employed for the numerical implementation of Eq.(131),accounting for an assembly of out-of-plane shear stresses alongside the entire lamination scheme:

Once the corrected valuesτ13(ijs),τ23(ijs)of stresses are found employing Eq.(131),the index association=τ13(ijs)and=τ23(ijs)is performed.In this way,the first order derivatives of,with respect toα1,α2can be performed within eachk-th lamina based on the GDQ algorithm:

Now, it is possible to solve the third differential relation in (126), while enforcing the stress compatibility condition for the surface external loadsat the boundary of the shell.Referring to the discrete point located at,one gets:

In the same way,the following relation should be fulfilled at the interface level between thek-th and(k+1)-th lamina,beingk=1,...,l-1:

Following a similar procedure to Eq.(132)for,shear stresses,the associationis performed,thus leading to the definition of the following vector:

Therefore,the load boundary conditions at the top surface of the shell read as[16]:

beingthe component of external loads applied atpoint.

Starting from the main outcomes of Eqs.(131)and(137),the corrected values of the out-of-plane 3D deformations collectedof thek-th layer can be calculated as well,according to what has been stated in reference[16].To this end,the constitutive law for generally anisotropic materials introduced in Eq.(51)is employed.Performing a computation in each discrete pointwithi=1,...,IN,j=1,...,IMandg=1,...,IT,it gives the discrete vector

where

From the inversion of the fundamental linear system of Eq.(138), the corrected values of3D strains are provided at each discrete point of the shellwithi=1,...,IN,j=1,...,IMandg=1,...,IT:

The corrected strain values can be employed for the correction of membrane stresses,,at each point of the discrete computational domain associated to the structure.

9 Numerical Investigations

We now present a series of case studies to validate the proposed methodology.Different structures have been considered, characterized by a variety of lamination schemes and curvatures.More specifically, the linear static response of zero-, singly- and doubly-curved panels have been investigated.Geometries of arbitrary shape have been considered,accounting for the generalized mapping technique of Eq.(13).The results provided by the ESL methodology presented in the manuscript have been compared to predictions obtained from refined 3D FEM simulations with parabolic elements as provided by a commercial software.The numerical investigations check for the accuracy of the model for different displacement field axiomatic assumptions according to Eq.(26), governing parameters of the employed distributions of external loads,and boundary linear springs distributions.

9.1 Materials Employed in the Analyses

In the present section we describe the mechanical input properties of materials for numerical analyses.For each component, the stiffness matrix E(k)is provided according to conventions from Eq.(49),referring to the material reference systemIn particular,a triclinic and a trigonal material have been employed belonging to the class of generally anisotropic materials.Referring to the class of orthotropic medium,a classic composite graphite-epoxy has been considered.Moreover,some considerations are provided for a 3D honeycomb lattice cell,characterized by an orthotropic softcore behaviour.

In the following, the 3D stiffness matrix of the triclinic materialhas been reported,settingwithi,j=1,...,6 the generalized stiffness constants:

Referring to the orthotropic graphite-epoxythe mechanical properties have been provided in terms of the nine independent engineering constants, namely the Young modulithe shear moduliand the Poisson’s coefficients

The correlations between quantities reported in Eq.(144) and anisotropic elastic constantsfori,j=1,...,6 have been reported in Appendix C.

A 3D lattice structure is now presented, called 3D Augmented Re-entrant Cellular Structure(3D ARCS), which has been extensively outlined in references [109,110].In Fig.3 a geometric representation of a 3D ARCS unit cell,together with all the geometric features,is presented.As can be seen,the main characteristics of the cell can be referred to a 2D pattern,which has been properly highlighted.When referring to the principal reference system of the shellO＇α1α2ζ,the heights of both vertical ribs are denoted withand,whereas the length of the inclined rib of the angleϑis identified withl.All the cell frames have a squared cross section of widtht.Due to the geometric and material symmetry of the unit cell,the equivalence of in-plane directionsα1,α2can be declared according to the Neumann’s principle.The constituent material is intended to be isotropic,characterized by an elastic modulusEsand a densityρs.It is useful to introduce for the unit cell,the dimensionless cell slenderness denoted withα,β,γ,η[109]:

Figure 3: Geometric representation of a 3D Augmented Re-entrant Cellular Structure (3D ARCS)unit cell.The configuration of the 2D unit pattern provides an auxetic behaviour of the equivalent continuum along in-plane directions α1,α2

Starting from reference [109], the classical beam theory is adopted to model each frame, here assumed as thin frame with rigid node area.A preliminary geometric overlapping verification should be performed.The vertical struts do not superimpose to each other if the following geometric relation is respected:

The following stiffnesses are conveniently introduced, accounting for the extensional behaviour and the bending deflections of the frame[109]:

Moreover,the generalized stiffness parameterφis defined as:

Based on the Eqs.(148)and(149),the quantitiesA,Bare defined as:

Following the procedure in reference [109], the 3D ARCSk-th layer can be homogenized as an orthotropic medium.The correspondent Elastic moduliE(k)1,E(k)2,E(k)3are expressed as:

In the following,the expressions for Poisson’s orthotropic coefficientscan be found:

In addition, the homogenized density of the 3D ARCS layer is computed according to the following expression,applying the procedure of void fraction introduced in reference[110]:

9.2 Straight Panels

The first set of numerical investigations considers some rectangular plates of constant thickness with general boundary conditions and subjected to a static uniform load.All the mechanic and geometric properties of the structure are reported in Fig.4.A convergence study starts considering a fully clamped rectangular plate of dimensionsLx= 1.5 m andLx= 1 m subjected to uniform pressure=-1.0×104Pa.The lamination scheme accounts for two external triclinic layers and a central orthotropic lamina made of graphite-epoxy.A reference 3D FEM solution has been provided,and the vertical deflection of the central part of the structure has been computed.The results are summarized in Table 1,with a clear fast convergence for all the displacement field assumptions.

Figure 4: Geometric and mechanical properties for an anisotropic rectangular plate enforced with non-conventional boundary conditions and subjected to a general loading

Table 1:Convergence analysis of a fully-clamped rectangular plated under uniform surface pressure.The vertical deflection at the centre of the structure has been computed by setting a different number of grid points,as well as various kinematic expansion orders.A reference 3D FEM solution has been provided for the comparison.The results are reported in 10-6 m

In particular, two different configurations have been considered.For the Case 01, four layers with general orientations(0/30/45/70)have been employed in the lamination scheme,embedding the triclinic material of Eq.(142)and the trigonal material introduced in Eq.(143).As far as the external constraints is concerned, the super elliptic distribution of Eq.(77) has been employed, leading to a half edge clamped condition.A general uniform pressure（α1,α2）= 1 has been applied on the top surface of the shell,characterized by both normaland tangential components,with respect toα1,α2,ζprincipal directions.A reference solution has been computed with a 3D FEM model,employing 371685 DOFs.The GDQ simulations, instead, were performed on a computational grid based on the LGL distribution of Eq.(100)characterized byIN= 37 andIM= 31.Different higher order theories have been considered according to Eq.(27),together with the zigzag effects.The static structural response has been reported in Figs.5-7.The 3D displacement field componentsU1,U2,U3through-the-thickness distributions of Fig.5 taken from the point located in the centre of the structure show that an almost linear distribution of in-plane components is predicted with all the assumptions,whereas a constant dispersion ofU3is assumed.

Figure 5: Through-the-thickness distributions atof the threedimensional displacement field vector U(α1,α2,ζ)of a rectangular plate composed of four generally anisotropic layers under static loads（=-10000 Pa）enforced with general boundary conditions（FC）

Figure 6: Through-the-thickness distributions atof the threedimensional strain vector ε(α1,α2,ζ) of a rectangular plate composed of four generally anisotropic layers under static loads（=-10000 Pa）enforced with general boundary conditions（）

On the other hand, the best agreement with the 3D FEM outcomes is provided by the EDZ4 theory for all the components.As far as the 3D strain vector dispersions are concerned (Fig.6),the GDQ approach predicts well the kinematic variables.In this case, the simulation with a fourth kinematic expansion order without the zigzag effect (ED4) seems to be more accurate, even though all the approaches provide very good results.The employment of a higher order assumption for the displacement field is crucial for theγ13distortion component, since in the third layer the 3D FEM predicts a slight nonlinear dispersion.The through-the-thickness distributions of stresses reported in Fig.7 employing the EDZ4 higher order assumption of the displacement field are perfectly in line with those belonging to the 3D reference solution.It can be seen that lower order ED1 and EDZ1 shell theories cannot give results of the same accuracy,even when embedding a zigzag function within the formulation.

Figure 7: Through-the-thickness distributions atof the threedimensional stress vector σ(α1,α2,ζ) of a rectangular plate composed of four generally anisotropic layers under static loads（=-10000 Pa）enforced with general boundary conditions（FC）

Figure 8: Through-the-thickness distributions atof the threedimensional displacement field vector U(α1,α2,ζ) of a fully-clamped (CCCC) rectangular plate composed of three generally anisotropic layers under static loads（=-5000 Pa）enforced with general boundary conditions（FC）

Figure 9: Through-the-thickness distributions atof the threedimensional strain vector ε(α1,α2,ζ)of a fully-clamped(CCCC)rectangular plate composed of three generally anisotropic layers under static loads（=-5000 Pa）enforced with general boundary conditions（FC）

Figure 10: Through-the-thickness distributions atof the threedimensional stress vector σ(α1,α2,ζ)of a fully-clamped(CCCC)rectangular plate composed of three generally anisotropic layers under static loads（=-5000 Pa）enforced with general boundary conditions（FC）

The second investigation was performed on a fully-clamped (CCCC) rectangular plate under a uniform static pressure along the ζ coordinate.The structure is characterized by two external layers of the triclinic material assessed in Eq.(142) and a central core made of an orthotropic 3D ARCS, whose mechanical properties are computed employing the homogenization provided by Eqs.(151)-(154).In particular,the geometric input quantities have been selected so that=0.0024 m,= l = 0.0012 m andThe tips and struts thicknesses are chosen so that t = 0.05 l.The constitutive material is intended to be isotropic with an elastic modulus Es= 196 GPa and a density equal toThe central isotropic layer is obtained from a superimposition of five cells.Thickness plots have been derived atand they have been reported in Figs.8-10.In this case, since the central layer behaves like a softcore, two different stable reference solutions have been computed,namely the 3D FEM and a fourth-order Layer-Wise(LD4)calculated by means of the GDQ method[19].For the sake of completeness,the simulation has been performed also with the EDZ4 displacement field theory,employing the strong formulation[27].Fig.8 shows that the employment of higher order theories,together with the zigzag assumption,is required for seeking results in line with other theories for both in-plane and out-of-plane displacement field components.This is due to the fact that the central softcore layer of 3D ARCS induces an unconventional in-plane deflection of the structure due to the static pressure.Similar considerations can be made for the 3D strain components of Fig.9.Note that both the 3D FEM and LD4 predict a behaviour typical of lamination schemes embedding softcore layers.Also for this case, the zigzag effect is very evident,since common ED-Ntheories are not adequate.Referring to out-of-plane distortionsγ13andγ23,the squeezing effect of the central core of the plate is outlined.Actually,the ESL simulations not including the zigzag axiomatic assumption of Eq.(27)are not capable of predicting the right structural response.The EDZ3 solution is in line with LD4 and 3D FEM predictions for all strain components.Also for the stress component vector through-the-thickness dispersions reported in Fig.10,the employment of the zigzag function leads to good predictions in line with reference simulations,regardless the kinematic expansion order.In the case ofσ(k)3three-dimensional stress component,higher order theories provide results in line with the 3D FEM outcomes.Lower order theories account only for the fulfilment of the external loads,because of the recovery algorithm assessed in Eq.(137).

9.3 Panels with Curvatures

In the present section the linear static analysis is performed for laminated structures characterized by the presence of curvature.Accordingly, a singly-curved shell and two doubly-curved panels have been considered,each of them enforced with generalized external constraints.The lamination schemes account for layers characterized by a general orientation with respect to the curvilinear principal reference system.As far as surface pressures are concerned, various loading conditions have been considered.They have been modelled within the ESL framework employing the distributions1,α2)outlined in Eqs.(65)-(67).

In Fig.11,there are all the information regarding the analysed cylindrical panel.According to the ESL approach, the shape of the structure is described starting from the reference surface, as shown in Eq.(2).In the following, the expression of the reference surface is reported employing principal curvilinear coordinatesα1,α2,settingRbthe radius of the circle[16]:

The stacking sequence accounts for three layers.Namely, a central triclinic material has been covered by two external sheets of graphite-epoxy of the same thickness.A uniform surface load=-1.0·104Pa has been applied to the top surfaceof the shell.Boundary conditions account for the clamping of the South (S) edge, whereas the East (E) side is constrained in a portion of its length by means of the super elliptic distribution described in Eq.(77).The structural response has been outlined in Figs.12-14 which refer to the dispersions of kinematic and static quantities along with the shell thickness in the region located at the centre of the physical domain.3D displacement field componentsU1,U2,U3have been reported in Fig.12.It is shown that the 3D FEM solution is properly predicted by the higher order models,especially for the out-of-plane component which cannot be predicted by lower order assumptions like the ED1 displacement field.The in-plane components are characterized by a linear distribution instead.On the other hand,an excellent agreement between different ESL approaches can be found in Fig.13,regardless of the employment of the zigzag function.Also,for this case,ED1 predictions are not in line with the outcomes from refined 3D FEM.As can be seen from Fig.14,also the 3D stress component vectors coming from different ESL approaches are in line with those obtained from the commercial software.Referring to theτ13,τ23andσ3dispersions,the exact fulfilment of the external load conditions can be easily seen.The worst accuracy is obtained whenN=1.

Figure 11: Geometric and mechanical properties for an anisotropic cylindrical panel enforced with non-conventional boundary conditions and subjected to general loads

Figure 12: Through-the-thickness distributions atof the threedimensional displacement field vector U(α1,α2,ζ)of a cylindrical panel composed of three generally anisotropic layers under static loads（=-10000 Pa）enforced with general boundary conditions（F）

Figure 13: Through-the-thickness distributions atof the threedimensional strain vector ε(α1,α2,ζ) of a cylindrical panel composed of three generally anisotropic layers under static loads（=-10000 Pa）enforced with general boundary conditions（F）

Next numerical investigation is performed on a doubly-curved panel of translation of mapped geometry.In particular,a hyperbolic paraboloid has been considered,whose reference surface equation in principal coordinates is reported[16]:

Figure 14: Through-the-thickness distributions atof the threedimensional stress vector σ(α1,α2,ζ)of a cylindrical panel composed of three generally anisotropic layers under static loads（=-10000 Pa）enforced with general boundary conditions（F）

beingkα1,kα2two characteristic parameters depending on the curvature of the principal lines.A geometric representation of the structure at issue is reported in Fig.15,together with all the lamination scheme information and the boundary conditions.The physical domain mapping has been assessed by means of the NURBS-based procedure of Eq.(13)with all the information related to each surface edge reported in Fig.15.External constraints have been modelled on the structure setting the West(W)side fully clamped.The employment of the Double-Weibull distribution defined in Eq.(78)accounts for the perfect clamping of the two external points of the East(E)edge.A uniform load has been applied in a quarter of the structure on the top surface.Within the GDQ simulation,such discrete variation of the load has been modelled by means of the Super Elliptic distribution reported in Eq.(67).In the out-of-plane direction,linear elastic springs follows a uniform dispersion（(ζ)=1）.In Fig.16,the dispersion of the displacement field componentsU1,U2,U3atare collected coming from various higher order theories.The refined 3D FEM outcomes are taken as references.It is shown that theU1deflection is well predicted by the proposed formulation with a higher order assumption of the displacement field.On the other hand,other in-plane field variable is not well predicted.As far as the vertical deflection is concerned, the structure at issue is characterized by a certain stretching along with the shell thickness.The shape of the distribution is depicted as well.The three-dimensional strain components are collected in Fig.17.It is shown that the 3D FEM solution can be fit by higher order assumptions of the displacement field together with the zigzag function.In particular, best performances are provided by the EDZ4 solution.The same considerations can be repeated for the three-dimensional stress vectorσ(α1,α2,ζ)of Fig.18.In fact, the shape of the distributions is well predicted by all the proposed theories.It is interesting to note that the implemented recovery algorithm leads to the perfect fulfilment of the equilibrium boundary conditions coming from external surface loads.

Figure 15: Geometric and mechanical properties of an anisotropic mapped Hyperbolic Paraboloid enforced with non-conventional boundary conditions and subjected to general loads

Figure 16: Through-the-thickness distributions atof the threedimensional displacement field vector U(α1,α2,ζ) of a hyperbolic paraboloid of arbitrary shape composed of three generally anisotropic layers under static loads（=-5000 Pa）enforced with general boundary conditions（F）

Figure 17: Through-the-thickness distributions atof the threedimensional strain vector ε(α1,α2,ζ) of a hyperbolic paraboloid of arbitrary shape composed of three generally anisotropic layers under static loads（=-5000 Pa）enforced with general boundary conditions（F）

The last numerical investigation has been performed on a super elliptic panel of revolution enforced with general external constraints along the West(W)mapped edge.Since the shell at issue is a doubly-curved revolution panel,the reference surface equation is computed starting from the general relation based on spherical coordinates[16]:

whereR0(α1)denotes the radius of the parallel,whereasx3(α2)is valuated in terms ofα2:

In Eq.(158),a1,a2are two geometric characteristic parameters,whereasndenotes the order of the meridian curve.In Fig.19 a representation of the domain distortion has been provided,together with the knot vector,the weights and the control points coordinate alongside the physical domain.

Figure 18: Through-the-thickness distributions atof the threedimensional stress vector σ(α1,α2,ζ) of a hyperbolic paraboloid of arbitrary shape composed o three generally anisotropic layers under static loads（=-5000 Pa）enforced with general boundary conditions（F）f

Figure 19: Geometric and mechanical properties of a super elliptic panel of revolution of mapped geometry enforced with non-conventional boundary conditions and subjected to general loads

It is interesting to note that the employed NURBS procedure allows tackling straight curves as a particular case.A general surface load has been applied to the structure, following a smooth distribution according to the Gaussian bivariate function of Eq.(67).The central core is orthotropic,in line with the 3D ARCS employed for the rectangular plate; since it has already been shown that the LD4 simulation predicts well the refined 3D FEM solutions, this time the LW model has been adopted as a reference structural response.Moreover, since the lamination scheme is characterized by a central softcore and two external layers of triclinic material (142), the ESL solution has been calculated only with the employment of the zigzag hypothesis on the displacement field and variousNth kinematic expansion orders,according to what exerted in Eq.(27).The static deflections,computed athave been collected in Fig.20.The in-plane displacement field components show that the EDZ4 assumption in both the strong and weak form best fits the results provided by the LD4 reference simulation in both triclinic layers and in the central 3D ARCS softcore.A great accordance is noticed between different approaches involving the strong and weak form of the EDZ4 structural theory.The shape of the distribution coming from the LD4 model is also well predicted by the ESL theory with zigzag functions characterized byN= 4.On the other hand, the ESL solution provides a linear distribution in the displacement field in-plane components, whereas the LD4 simulation accounts for a nonlinear dispersion due to the softcore behaviour of the lattice layer.Referring to the through-the-thickness dispersions of the three-dimensional strains reported in Fig.21,the squeezing phenomenon due to the presence of the 3D ARCS is evident.The EDZ4 theory in both its strong and weak formulation well predicts the deformations of the two external triclinic layers.Referring to the central core,the structural response is strongly dependent on the choice of the axiomatic displacement field thickness function.The employment of a higher order for the kinematic expansion provides results close to the LW outcomes,especially for the case of out-of-plane distortionsγ13,γ23and stretching componentε3.

Figure 20: Through-the-thickness distributions atof the threedimensional displacement field vector U(α1,α2,ζ)of a super elliptic panel of revolution of arbitrary shape composed of two generally anisotropic layers under static loads（=-3000 Pa）enforced with general boundary conditions（FCF）

Figure 21: Through-the-thickness distributions atof the three dimensional strain vector ε(α1,α2,ζ) of a super elliptic panel of revolution of arbitrary shap composed of two generally anisotropic layers under static loads（=-3000 Pa）enforced wit general boundary conditions（FCF）-eh

Some concluding remarks are made on the three-dimensional stress componentsσ(α1,α2,ζ),as visible in Fig.22.Since the equivalent elastic stiffnesses provided via the implementation of the engineering constants from Eqs.(151)and(152)are lower than those of a triclinic material provided in Eq.(142),in the central area of the structure very low in-plane normal and shear stresses are obtained.On the contrary,the two outer layers provide higher values of stresses,following a linear dispersion.A parabolic distribution is traced in theτ13,τ23andσ3representation, due to the enforcement of the equilibrium boundary conditions.A great accordance can be found between the LD4 reference simulation and the results obtained from the ESL solution regardless the maximum order of the kinematic expansion.

Figure 22: Through-the-thickness distributions atof the threedimensional stress vector σ(α1,α2,ζ) of a super elliptic panel of revolution of arbitrary shape composed of two generally anisotropic layers under static loads（=-3000 Pa）enforced with general boundary conditions（FCF）

Above all, it is shown that the proposed formulation is very reliable for the static study of arbitrarily-shaped laminated structures under static loads.Moreover, the accuracy of the model is not affected by the presence of a layer with lattice features.

10 Conclusions

In the present work,an ESL formulation based on higher order theories has been proposed for the linear static analysis of structures with double curvature and arbitrary shape accounting for a general assessment of external boundary conditions and lamination schemes,involving a generalized in-plane and out-of-plane dispersion of linear elastic springs.Furthermore, a general shape of the external loads has been provided based on various bivariate distributions.The formulation has been developed following a generalized variational approach, accounting for an interpolation of the field variable of higher order based on Lagrange polynomials and a generic number of grid points.Differently from classical variational approaches, the structural problem is now characterized by refined shape functions,with a significant reduction in the computational cost.As a matter of fact,the computation of the fundamental matrices is performed by means of a numerical integration throughout the entire physical domain.

In the post-processing phase, an equilibrium-based recovery procedure has allowed for an exact fulfilment of boundary conditions coming from the externally applied surface loads.A series of validating examples has been proposed, where the static response of structures with different curvatures has been successfully checked.The accuracy of the solution has been evaluated from a direct comparison of the results obtained from refined 3D and LW theories.The proposed formulation has been demonstrated to be very reliable for disparate situations in terms of materials,lamination schemes and load conditions.

Funding Statement:The authors received no specific funding for this study.

Conflicts of Interest:The authors declare that they have no conflicts of interest to report regarding the present study.

Appendix A.Higher order coefficients for boundary conditions

In the previous paragraphs, referring to a genericτ-th kinematic expansion order withτ=0,...,N+ 1, a relation has been provided in Eq.(57) in a compact matrix form that computes the generalized stress resultant vector S(τ)αiof Eq.(52) forαi=α1,α2,α3in terms of generalized displacement field component vector u(τ)introduced in Eq.(26).In the following, we report an extended version of coefficientsof the matrix O(τη)αiαj, beingi,j,r= 1,2,3 andg= 1,...,9.Since the fundamental governing equations have been assessed in aIN×IMdiscrete grid,fundamental coefficients have been provided in the correspondent discrete formdefined as follows:

whereςα1(n)α2(m)contains the GDQ weighting coefficients for the derivative of the(n+m)-th order referred to the principal directions of the shellα1,α2.For the sake of completeness,it should be recalled that the symbol°stands for the Hadamard product.In the following,the interested reader can find the extended version of,settingg=1,...,9,s=1,2,3 andi,j,r=1,2,3.

First column of the generalized operator

Second column of the generalized operator

Third column of the generalized operator

Appendix B.Higher order fundamental stiffness matrix coefficients

In the present section we report the complete expression of the fundamental coefficientsof the stiffness matrix of Eq.(62),defined for eachτ,η=0,...,N+1 andαi,αj=α1,α2,α3principal direction of the shell.As can be seen,a formulation employing the Kronecker tensorial product ⊗is provided,as well as vectorslα1,lα2of the Lagrange interpolating polynomials of Eq.(34),as well as their first order derivatives,For the sake of completeness,the terms of the generalized stiffness matrix A(τη)αiαjare computed according to the expression reported in Eq.(54)for eachτ,η=0,...,N+1 andαi,αj=α1,α2,α3.All the terms at issue are collected by row.

First row of the fundamental stiffness matrix

Second row of the fundamental stiffness matrix

Third row of the fundamental stiffness matrix

Appendix C.Constitutive relationship for an orthotropic medium

With particular references of numerical investigations accounting for orthotropic materials, the expressions ofwithi,j= 1,...,6 andk= 1,...,lof the three-dimensional constitutive relationship reported in Eq.(49)should be considered.Actually,the mechanical properties of thekth lamina are provided in terms of the well-known engineering constants,namely the elastic moduli,,,the shear moduli,,and the Poisson’s coefficients,,[16]:

In the previous equation,quantityΔ(k)reads as follows:

It should be recalled that the Poisson’s coefficients satisfy the well-known reciprocity relation[16]: