|Computer Modeling in Engineering & Sciences|
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
Department of Innovation Engineering, School of Engineering, University of Salento, Lecce, 73100, Italy
*Corresponding Author: Francesco Tornabene. Email: firstname.lastname@example.org
Received: 28 February 2022; Accepted: 11 May 2022
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 soft-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 non-uniform 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
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 . 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 . 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 . An important reason comes from the endeavor to obtain the best form under fixed load cases and external constraints layup . 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 .
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 . 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  an extensive study with higher order theories for doubly-curved shell structures is provided, and in reference  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  the dynamic behaviour of anisotropic doubly-curved shells is investigated with a higher order LW approach, whereas in reference  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 . 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 .
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  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 , 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  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 . 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  it is shown that the out-of-plane stretching plays an important role in the deformation response, whereas in reference  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  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  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  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  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  moving concentrated loads have been applied to 1D structures starting from the definition of the Dirac-Delta function. Then, in reference  the formulation has been extended to 2D structural problems. Moving load problems have been treated in reference  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  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  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  FGM structures have been investigated by employing a semi-analytical model. In reference  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 , 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 post-processing 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) . 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 doubly-curved 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 . 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 . 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 , 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.
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 system of unit vectors in terms of three parameters :
being the three-dimensional position vector, and for a continuous function. If the variables are taken as a set of curvilinear coordinates, it is possible to introduce the above discussed reference surface starting from Eq. (1) as follows :
where the has been selected alongside the normal direction of . 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 with respect to , namely , are 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 that . Accordingly, shell thickness can be defined starting from Eq. (2), according to the following expression taken from reference :
where is the width of the -th lamina of the stacking sequence composed of layers located within the interval . If the second order derivatives of the reference surface with are calculated, it is possible to compute the principal radii of curvature in each point of the physical domain:
In addition, the well-known Lamè parameters are defined as follows:
Accordingly, in the case of straight shells along the principal direction, Eq. (7) turns into .
In the present manuscript, doubly-curved shells of variable thickness  have been considered; therefore, a generalized analytical expression of is provided hereafter:
being the reference height of the -th layer and for a scaling parameter, whereas is a constant dimensionless shift number. In Eq. (8), thickness variation has been described in terms of four analytical expressions employing a dimensionless coordinate defined along each principal direction of the shell, according to the following equation:
In the present formulation, the expressions of for have been implemented, being , and :
When an arbitrary curve is considered lying on the reference surface described with principal coordinates, a local reference system composed of three unit vectors , and should be defined in each point of the curve at issue. More specifically, accounts for the tangential unit vector of the curve, whereas is the normal unit vector of the curve. denotes the bi-normal direction. Referring to middle surface reference system, , and components can be computed in each point of the curve, leading to:
Since the curve at issue lies on , it should be noted that and .
In the previous paragraph, the geometry of a generic doubly-curved shell has been described starting from the reference surface , according to the ESL framework. Nevertheless, it is possible to perform a distortion of the original physical domain in order to take into account structures of arbitrary shape. To this purpose, a set of blending coordinates 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 :
In Eq. (12), the symbols with and account for the parametric representation of the outer sides of the reference surface alongside the physical domain, collected in the vectors for . In the same way, the four corners of denoted with for are arranged so that with . Based on these definitions, it is possible to assess the blending expressions of the physical domain according to the following matrix-form expressions:
being for . Accordingly, the following definitions for and should be assessed:
where is the identity matrix. On the other hand, vectors and matrix read as follows:
where the definitions for are introduced. In the case of arbitrarily-shaped physical domains, it is very likely that the partial derivatives with respect to should be referred to the principal coordinate system as well. For this reason, the following transformation should be declared, setting the Jacobian matrix of Eq. (13):
If , the inverse form of the Jacobian matrix occurring in Eq. (17) can be calculated, thus obtaining:
Accordingly, the inverse relation of Eq. (17) accounts as follows:
For the sake of completeness, first order partial derivatives with of with respect to natural coordinates can 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 principal coordinates in terms of and . The final expression for second order mixed derivatives looks as follows :
being with a transformation coefficient reading as:
In addition, one gets the following expression for pure second order derivatives with respect to direction, setting :
where the coefficients for complete expression are thus reported:
In the present work, the -th edge of the mapped physical domain for , whose components along direction have been collected in the vector for , has been geometrically parametrized with respect to for in terms of NURBS curves .
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 associated with each point of the solid of the Euclidean space, defined in Eq. (2) with respect to the principal set of coordinates , can be expressed according to the following equation (Fig. 1):
Employing the unified notation of Eq. (26), a generalized displacement field component vector is assessed for each -th order of the kinematic expansion, defined in each point of the reference surface . The dependence of from the out-of-plane coordinate is taken from the definition, for each , of a series of thickness functions along each principal 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, if are 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 :
Referring to , a dimensionless coordinate has been defined in each -th layer of the stacking sequence starting from the through-the-thickness coordinate of the top and the bottom surface delimiting the -th sheet, respectively. If with l denoting 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 is adopted. Letter “E” tells that an ESL formulation is considered. “D” stands for the displacement field as an unknown variable, whereas N denotes the maximum order of the kinematic expansion. When the additional generalized displacement field component vector 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 component .
Once the ESL assessment of the displacement field component vector has 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 of points, whose generic element can be identified with for and . For this purpose, an interpolation algorithm is introduced for at each kinematic expansion order. A very useful tool is the vectorization operator , accounting for a rearrangement with a 1D vector of a 2D matrix. If we denote with a matrix of dimensions of components with and , a column vector of dimensions is defined, whose arbitrary component for reads as:
Since all quantities in Eq. (26) are evaluated in each point of the discrete computational grid, quantities are introduced with 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 employing a matrix of generalized shape functions:
In a more expanded form, one gets :
where for in a generic point is obtained from the interpolation of all the values assumed by the generalized displacement field component in each point of the above introduced discrete grid. To this end, the well-known univariate Lagrange interpolating polynomials of -th order defined along the principal direction of the physical domain are employed, with and . According to reference , can be computed as follows:
Starting from Eq. (33) it is possible to define the generalized shape function matrix alongside the physical domain, whose dimensions depend on the selected computational grid. Eventually, one gets:
In the same way, vector accounting for the first order derivative of the Lagrange interpolating polynomials collected according to Eq. (36) with respect to is defined as follows, setting :
Since the Kronecker tensorial product, denoted with , has been adopted in Eq. (35), it is useful to introduce a row vector of dimensions so that matrix can be expressed in a compact form as follows :
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 (with ) read as :
being the 3D strain component vector and a differential operator relating to the displacement field component vector . Accordingly, the operator accounts for the partial derivatives with respect to the shell principal direction , as well as the curvature thickness parameters :
In the same way, for read as:
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 -th order :
In the previous equation, the strain vector 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, is introduced for each and for according to the following extended form :
In Eq. (44) it has also been shown that the generalized strain component vector embeds the ESL assessment of the displacement field of Eq. (26). Based on the generalized interpolation procedure performed in Eq. (32), can be related to the discrete vector of the generalized displacement field of Eq. (31) for each , i.e.,
where the kinematic operators with are defined from the previously discussed higher order interpolation algorithm performed on the discrete grid by means of the tensorial Kronecker product :
We recall that in Eq. (47) the Lagrange polynomials employed for the interpolation of the field variable along the in plane direction, as well as their first order derivative, have been collected in the vectors and , respectively.
The present manuscript investigates the problem of the static structural response of a doubly-curved 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 of l laminae has been contemplated, whose -th layer with is characterized by a thickness according to the conventions of Eq. (8) and a reference system orientated along the material symmetry axes. Accordingly, the formulation has been developed by assuming . In addition, a characteristic angle describing the orientation of axis with respect to is introduced. Each layer of the stacking sequence is intended to be generally anisotropic. The mechanical behaviour is described with respect to the 3D stress and strain component vectors and from the introduction of a generally anisotropic stiffness matrix for each -th layer. The generalized constitutive relationship can be thus defined as :
Since a linear elastic theoretical assessment of the structural problem is performed, the matrix reads as follows:
where with relates the -th component of vector to the -th strain component of the vector. Accordingly, the coefficients employed in Eq. (49) account for the 3D elastic stiffnesses of the material, namely . When is referred to the plane stress assumption , the reduced stiffness coefficients can 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 vectors and . To this end, a generalized rotation matrix is defined according to reference  within each layer starting from the previously discussed angle employed to compute the rotated stiffness matrix with components for referred to the shell geometric reference system . More in detail, Eq. (48) takes the following form:
On the other hand Eq. (51), defined for the -th layer of the 3D solid, should be assessed within the ESL framework accounting for all the l laminae 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 vector , associated to each -th order of the kinematic expansion of Eq. (39). In this way, Eq. (51) can be expressed in terms of and vectors, according to the following equations directly derived from the static equivalence principle:
where is the generalized ESL stiffness matrix :
whose generic component is defined as:
being and . As can be seen, four blocks can be traced in Eq. (53) accounting for the presence of the derivatives of the -th order within the expression reported in Eq. (54). Furthermore, in Eq. (54) the term refers to the 3D rotated stiffness coefficient introduced in Eq. (51), possibly adjusted by the shear correction factor  in the case of displacement thickness functions employed in Eq. (26) neglecting the out-of-plane stretching effects:
Thus, the generalized ESL assessment of the constitutive relationship performed in Eq. (56) can be expressed so that each component of can be expressed in terms of the generalized displacement field component vector for 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 .
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 top and the bottom surfaces of the shell. Moreover, the formulation embeds the contribution of linear elastic springs distributed along the shell edges. Suppose we denote with the 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 :
being the total elastic strain energy of the 3D solid, expressed via the ESL definition of the displacement field, as reported in Eq. (39).
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 and defined in each -th lamina of the stacking sequence concerning the geometrical principal reference system , according to the following relation:
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 Appendix B, an extended expression of each stiffness matrix coefficient is reported, setting and .
Let us consider two generic vectors and of external loads applied to the doubly-curved shell structure at issue, whose components are referred in the above discussed geometric reference system as shown in Fig. 2:
where symbols and account for quantities referred to and , respectively. Accordingly, a general variation dependent on in-plane principal coordinates can be assigned to each component for and of external load vectors introduced in Eq. (63). If we denote with a reference scaling value of loads applied along the principal direction, the following relation can be assessed (Fig. 2):
For a constant loading distribution, it is:
In the present manuscript, two different distributions for components with are declared. A bivariate Gaussian distribution is now reported:
being for the position parameter acting on the in-plane principal direction with a variance scaling factor equal to and a correlation factor between denoted with . 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 for with respect to principal shell coordinates. In addition, is the location in the physical domain of the centre of the ellipse, whereas n is a power exponent. For , 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 of the 3D displacement vector referred to the top and the bottom surface is considered, the virtual work of surface loads of Eq. (63) can be computed as:
being for the through-the-thickness curvature parameters calculated by means of Eq. (7) for and . On the other hand, if we consider a virtual variation of the generalized displacement field component vector introduced in Eq. (26), a new generalized vector of external actions lying on the reference surface can be introduced for each . Thus, the virtual work reads as follows:
Starting from Eq. (70), the extended expressions of the components are obtained for the -th kinematic expansion order, being for the thickness function employed for the assessment of the displacement field components according to Eq. (26) calculated at and , respectively:
Since the generalized displacement field component vectors with have been interpolated with a generalized shape functions set in Eq. (39), the virtual work of generalized external actions already computed in Eq. (69) should be expressed in terms of vector taking into account Eq. (32), leading to the following expression :
In Eq. (72), the generalized external actions referred to each -th order have been collected in the column vector , defined as:
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 stiffness orientated along the shell side located at with and distributed along the principal direction has been considered, according to the theoretical development reported in reference . The superscript k stands for the stacking sequence layer. Accordingly, they induce some boundary stresses after applying the corresponding 3D displacement component If the physical domain consists of a rectangular 2D interval , different stresses are induced on each boundary. Referring to for , the following definitions can be stated at the -th lamina level, setting :
where accounts for an arbitrary univariate function in the coordinate located at . On the other hand, stresses can be found at according to the following definitions:
being . In the same way, Eq. (75) contains a univariate dispersion of linear elastic springs denoted with for . In order to provide an effective expression of both and regardless of the actual dimension of the physical domain, a new dimensionless coordinate is introduced, according to the following definition:
As a matter of fact, the Super Elliptic distribution (S) of the -th order has been employed, according to the following analytical expression, setting and the position and the scaling parameter, respectively:
In the following, the expression of a Double–Weibull distribution (D) has been reported in terms of :
where are two position parameters and p is the scaling value.
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 with :
Referring to the edges of the structure located at for , Eq. (75) turns into:
Having in mind the ESL assessment of general distributions of boundary springs, the virtual work contribution to Eq. (58) can be computed as the sum of and , referred to the edges of the physical domain characterized by and , respectively, for :
where components are defined according to what exerted in Eq. (80). accounts for the generalized boundary stressed reported in Eq. (80). Introducing in Eq. (82) the weak form assessment of the displacement field of Eq. (32), the final form of comes out:
Thus, the generalized displacement vector of 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 matrix has been introduced, accounting for the ESL generalized stresses induced by the springs distributed along the edges of the physical domain, at each kinematic expansion order, as well as the generalized shape functions arranged in , based on the Lagrange interpolation procedure of Eq. (32).
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 :
For shells of arbitrary shapes, the employment of the generalized blending functions assessed in Eq. (13) with the Jacobian matrix 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 matrix and the external surface loads should be adjusted appropriately for each since integrals with respect to principal directions should be assessed in terms of natural coordinates . Referring to a specific th kinematic expansion order, stiffness matrix defined in Eq. (62) turns into the following relation:
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 directions employing the symbols . The generalized boundary springs along the -th edge of the structure with accounts as:
Since the natural coordinate set has been defined in the dimensionless rectangular interval , Eq. (90) can be written as:
Once the quantity has been successfully computed in the previous equation for each -th order, a kinematic expansion is performed such that Eq. (92) is embedded in the global variational matrices of Eq. (87). In this way, the global stiffness matrix accounts for the 3D assessment of the boundary linear elastic springs introduced for the definition of natural and non-conventional boundary conditions, whereas the term is employed for the external surface loads.
In order to compute the generalized stress resultants induced by linear elastic springs for a mapped geometry employing the approach of Eqs. (74) and (75), it is useful to express the generalized displacement field referred to local reference system in terms of for each -th order of the kinematic expansion :
In the same way, the generalized stress resultants should be expressed in terms of the vector for according to the following relations :
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 belonging to the dimensionless squared computational domain, the -th order derivative of a bivariate function concerning natural coordinates can be expressed as :
where denote the number of the selected discrete points along directions, respectively. As can be seen, Eq. (95) accounts for a quadrature rule for derivative purposes. The weighting coefficients for , , and are computed utilizing the recursive expression valid for each :
Eq. (96) is based on the interpolation of the unknown function employing the Lagrange interpolating polynomials . In the case of pure derivatives with respect to , it gives and , respectively, and GDQ coefficients come into the following relation:
being the Kronecker delta operator. For a generic non-uniform grid of points, the GDQ weighting coefficients for the derivative of the -th order are calculated by means of Eqs. (96) and (97) are collected in a squared matrix denoted with . Starting from Eq. (96), the matrix at issue can be obtained from a proper expansion of matrix for and of dimensions with , respectively, accounting for the derivatives of -th and -th order along :
Based on Eq. (97), it is . In this way, the numerical derivation employing the GDQ method with respect to can be performed as:
where is a matrix whose general component is for and . In the same way, the generic element of is .
A 2D Legendre-Gauss-Lobatto (LGL) grid of dimensions referred to the domain has been here adopted, such that:
Since the fundamental governing relations (86) account for derivations along principal directions, we consider the GDQ rule of Eq. (95) referred to natural coordinates . The matrix approach is then followed, starting from Eq. (98). Employing the generalized blending functions reported in Eq. (13), the matrices and for the first order derivatives along of a generic bivariate function can be computed as:
being with the vectorized form by means of Eq. (29) of the coordinate transformation computed by Eq. (19). In the same way, starting from Eq. (22) and Eq. (24) the discrete form of second order derivatives can be assessed introducing the quantities with :
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 entire computational 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 , whereas those referring to the former are arranged within the vector . Thus, the fundamental system accounts as follows :
Performing a static condensation of Eq. (104) with respect to according to the procedure reported in Eq. (16), the final form of the discrete system of dimensions comes out, being the identity matrix:
The numerical integrations that are involved in the formulation of the present manuscript are performed by means of the GIQ procedure . Referring to a generic univariate function defined in the closed interval with , a discrete grid is selected from the 1D domain, namely with . Accordingly, the integral of f restricted to the sub-interval with belonging to the grid at issue can be computed as a weighted sum of the values for assumed by the function at each point of the computational grid of :
where are referred to the interval . The coefficients are computed following the extensive procedure of reference , starting from the shifted GDQ weighting coefficients reported hereafter, setting :
Then, they are collected in the matrix of dimensions . Thus, the coefficients referred to the integral of f restricted to the interval are the elements of the GIQ matrix, reading as:
All the terms introduced in Eq. (106) can be collected in a vector of dimensions denoted with , setting x the integration variable. In the same way, the vector can be introduced when the integration of the univariate function is required. For a bivariate function , Eq. (106) turns into the following expression, setting a generic 2D grid from a rectangular domain of extremes and :
where values for are 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 , as reported:
Referring to the bivariate function , it is also possible performing a univariate integration along a parametric line of the rectangular computational domain parallel to x and y, respectively, employing the GIQ rule. For this purpose, and are defined starting from Eq. (110):
For the sake of completeness, symbol are row vectors of ones of dimensions and , 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 employing principal coordinates. As it has been shown, for arbitrarily-shaped structures for is described by means of natural coordinates in Eq. (88). For this purpose, it is useful to define a matrix of dimensions for the transformation at issue:
being a ones squared matrix of dimensions . Accordingly, following Eq. (110), whereas , , and accounts for the determinant of the Jacobian matrix , as well as the Lamè parameters of the shell. If the index with and is assessed, we obtain the vectorized form of , accounting for the grid of the computational domain:
The edge characterized by behaves as:
On the other hand, for it is:
Finally, the edge identified with reads as:
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 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 of the computational points is provided for the dimensionless interval :
being the number of the selected discrete points. Referring to the interval for representing the -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 vector with , and for each -th layer of the stacking sequence :
being the generalized thickness functions set calculated at the point of the 3D solid, whereas are the generalized displacement field components for the -th kinematic expansion order referred to the point of the physical domain, for , . In the same way, the generalized strain component vector for and are computed in terms of the discrete 3D strain component vector within each -th lamina, leading to the following expression:
Referring to the generic point of the shell structure, membrane stresses can be calculated according to the generally anisotropic law of Eq. (51), leading to the following relation, setting for the elastic stiffness coefficients of the -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  expressed in curvilinear principal coordinates are employed, instead, as:
where are the principal curvature radii of the reference surface , whereas the terms depend on the quantities calculated in Eq. (125), reading as:
It should be noted that the expression of 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 in-plane stresses with respect to occurring in the first two expressions of Eq. (127) should be performed employing the GDQ method for each -th layer of the lamination scheme, setting :
The first order differential relations (126) should fulfill a single boundary condition. Referring to the first lamina of the stacking sequence , the external loads applied at the bottom side of the structure satisfy the following boundary conditions:
being the in-plane static loads referred to the point of the bottom skin of the structure following the convention in Eq. (63). On the other hand, for an arbitrary -th layer, the first two relations of Eq. (126) can be solved at each discrete grid point by means of the GDQ computations of Eq. (128) enforcing the unknown stress compatibility conditions between the -th and the -th laminae as follows:
For the surface loads and at 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 :
In the previous equation an index has been introduced, accounting for points defined according to Eq. (122) for the -th lamina alongside the entire shell thickness for and . As a matter of fact, stresses come from the direct solution of the balance Eq. (126). The introduction of the index s leads 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 of stresses are found employing Eq. (131), the index association and is performed. In this way, the first order derivatives of with respect to can be performed within each -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 loads at 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 the -th and -th lamina, being :
Following a similar procedure to Eq. (132) for shear stresses, the association is performed, thus leading to the definition of the following vector:
Therefore, the load boundary conditions at the top surface of the shell read as :
being the component of external loads applied at at each point.
Starting from the main outcomes of Eqs. (131) and (137), the corrected values of the out-of-plane 3D deformations collected of the -th layer can be calculated as well, according to what has been stated in reference . To this end, the constitutive law for generally anisotropic materials introduced in Eq. (51) is employed. Performing a computation in each discrete point with , and , it gives the discrete vector :
From the inversion of the fundamental linear system of Eq. (138), the corrected values of 3D strains are provided at each discrete point of the shell with , and :
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.
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.
In the present section we describe the mechanical input properties of materials for numerical analyses. For each component, the stiffness matrix is provided according to conventions from Eq. (49), referring to the material reference system . In 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 material has been reported, setting with the generalized stiffness constants:
The trigonal material accounts for the following stiffness matrix :