iconOpen Access


Static Analysis of Anisotropic Doubly-Curved Shell Subjected to Concentrated Loads Employing Higher Order Layer-Wise Theories

Francesco Tornabene*, Matteo Viscoti, Rossana Dimitri

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

* Corresponding Author: Francesco Tornabene. Email: email

(This article belongs to this Special Issue: Theoretical and Computational Modeling of Advanced Materials and Structures)

Computer Modeling in Engineering & Sciences 2023, 134(2), 1393-1468. https://doi.org/10.32604/cmes.2022.022237


In the present manuscript, a Layer-Wise (LW) generalized model is proposed for the linear static analysis of doublycurved shells constrained with general boundary conditions under the influence of concentrated and surface loads. The unknown field variable is modelled employing polynomials of various orders, each of them defined within each layer of the structure. As a particular case of the LW model, an Equivalent Single Layer (ESL) formulation is derived too. Different approaches are outlined for the assessment of external forces, as well as for non-conventional constraints. The doubly-curved shell is composed by superimposed generally anisotropic laminae, each of them characterized by an arbitrary orientation. The fundamental governing equations are derived starting from an orthogonal set of principal coordinates. Furthermore, generalized blending functions account for the distortion of the physical domain. The implementation of the fundamental governing equations is performed by means of the Generalized Differential Quadrature (GDQ) method, whereas the numerical integrations are computed employing the Generalized Integral Quadrature (GIQ) method. In the post-processing phase, an effective procedure is adopted for the reconstruction of stress and strain through-the-thickness distributions based on the exact fulfillment of three-dimensional equilibrium equations. A series of systematic investigations are performed in which the static response of structures with various curvatures and lamination schemes, calculated by the present methodology, have been successfully compared to those ones obtained from refined finite element three-dimensional simulations. Even though the present LW approach accounts for a two-dimensional assessment of the structural problem, it is capable of well predicting the three-dimensional response of structures with different characteristics, taking into account a reduced computational cost and pretending to be a valid alternative to widespread numerical implementations.


1  Introduction

Laminated materials are very often required in many engineering applications. In particular, an increasing need for structures with complex geometric shapes characterized by smart non-conventional fabrics is much more evident [14]. In this way, an optimization of the material properties can be obtained from the model, since the lamination scheme can be selected according to the structural needs, especially when the geometry cannot be varied due to architectural and functional requirements [5,6]. On the other hand, the introduction of curvature provides an optimization of the stress distribution. Nevertheless, classical models can lead to erroneous predictions due to an unusual structural behaviour coming from the absence of material and geometric symmetry planes [710]. In order to prevent the invalidation of the design process and exposing the final product to some safety risks, new simple but accurate methodologies should be developed. Furthermore, a proper mathematical modelling of complex appliances, embedding all the curvature effects and constitutive couplings, can be very cumbersome in its conception as well as computationally demanding, due to the huge number of variables occurring in the structural problem [11,12]. Three-dimensional elasticity solutions are the most accurate approaches for the correct prediction of a structural response [1315]. However, in the case of anisotropic materials and complicated structural shapes, a closed-form solution can not be easily found. Therefore, a large numerical system should be developed with a significative number of Degrees of Freedom (DOFs). This is likely to come across several numerical issues like the computational stability. For this reason, simplified two-dimensional formulations have been developed throughout literature, so that the solution is found with a lower computational effort [1620]. In particular, two main approaches can be traced, namely the Equivalent Single Layer (ESL) and Layer-Wise (LW) formulations [2124]. According to ESL, the three-dimensional doubly-curved solid is reduced to a surface located in its middle thickness. In this way, a 2-manifold is derived with its geometric parameters describing the shape of the actual structure [25,26]. Moreover, the field variable, as well as the primary and secondary ones, are reduced to the surface at issue. A key aspect of this approach is the homogenization of the stacking sequence so that a set of equivalent properties can be computed. In contrast, the two-dimensional LW implementations account for a displacement field expansion within each lamina, thus providing an accurate local description of all the mechanical quantities. Furthermore, compatibility conditions are developed between two adjacent laminae so that the consistency of the whole model is assessed [27,28]. Within the LW approach, the accuracy of the solution may be increased if the order of the adopted interpolating polynomials gets higher, or equivalently if a generic lamina is divided in some virtual sub-layers following a local-global strategy, thus leading to sub-laminate formulations [2936]. As far as the displacement field assumption is concerned, classical approaches provide a linear through-the-thickness assumption of the unknown in-plane field variables, whereas a rigid behaviour is assumed in the out-of-plane directions. In this way, no stretching effects can be predicted, as well as the softcore behaviour of the lamination scheme [37]. In addition, a smooth displacement field assumption does not fit the actual interlaminar deformation during the deflection of the structure. For this reason, a higher order axiomatic expression for the out-of-plane variable is crucial in both LW and ESL approaches, even though the latters should embed in themselves with the well-known zig-zag function in order to obtain accurate results [3843]. A milestone for the assessment of the kinematic field variable is its unified description [44,45], which incorporates both higher order theories and classical approaches like the First Order Shear Deformation Theory (FSDT) and the Third Order Shear Deformation Theory (TSDT), outlined in references [46,47]. Since it allows the introduction of a generalized set of thickness functions in both in-plane and out-of-plane directions, polynomials, non-polynomials and trigonometric functions can be employed to this purpose [4850], depending on the material properties and relative thickness of the laminae embedded in the structure. Moreover, it is possible to develop an advanced set of axiomatic thickness functions, whose expressions is obtained from the mechanical shear properties of the lamination scheme, leading to the so-called refined zig-zag theories [5153].

The fundamental governing equations have been analytically solved only for simple geometries and lamination schemes, such as frames, plates, cylinders and spherical panel. On the other hand, only cross-ply orthotropic laminates can be adopted, otherwise no mathematical procedures are available at the moment for the solutions of such differential sets of equations. For this reason, in references [54,55] the ESL structural problem for structures with single and double curvatures accounting for a generally anisotropic lamination scheme is numerically developed by means of the Generalized Differential Quadrature (GDQ) method. We recall that such methodology discretizes directly the derivatives of a given function, thus allowing to solve the problem directly in the strong form, as it has been shown in references [5660]. Belonging to the class of spectral collocation methods, it embeds a series of numerical techniques like the gaussian quadrature and the classical differential quadrature [61,62], accounting for a rectangular computational grid. The accuracy of the method comes from the proper selection of a set of discrete points from the physical domain, as well as the computation of the quadrature coefficients. In references [63,64] its accuracy has been compared to other numerical techniques. Moreover, it has been shown that it is a very reliable procedure for the analysis of lattice three-dimensional cells [6567], as well as Functionally Graded Materials (FGMs) employing a reduced computational cost [6871]. Moving from the above discussed GDQ procedure, the Generalized Integral Quadrature (GIQ) method turns out to be an effective strategy for the numerical implementation of integrals with a domain collocation strategy. In references [72,73], the interested reader can find an extended treatise on the topic.

In classical variational approaches like the well-known Finite Element Method (FEM), an axiomatic set of shape functions are provided, leading to a weak form of the governing equations [7476]. However, this methodology induces some drawbacks in the solution due to the discrepancy between the geometry and the adopted interpolating function. As a matter of fact, this problem can be overcome if a domain is discretized with a fine mesh, which in turn increases the computational cost. In contrast, the Iso-Geometric Approach (IGA) adopts the actual geometry of the structure employed in the Computer Aided Design (CAD) procedure itself for the assessment of the unknown field variable of the governing equations. Starting from the pioneer works reported in reference [77], the IGA methodology has been successfully applied to several structural problems related to arbitrary shapes [78,79]. It has been shown that IGA is a very efficient methodology in the case of structures of arbitrary geometries, especially when a significative domain distortion is required. In particular, the best performances of such methodology are reached when Non-Uniform Rational Basis Spline (NURBS) curves with higher order basis functions are adopted because they show a significative computational stability. In references [80,81], the curves at issue are presented in a comprehensive way, together with an iterative procedure for their computation. Furthermore, in references [8283], a meshfree Galerkin method is applied for the buckling analysis of shallow shells with single and double curvatures accounting for geometric interpolation functions for the assessment of the unknown field variable. In reference [84], the meshfree approach has been adopted for the finite rotation analysis of structures with different curvatures.

Another topic related to the analysis of shell structures relies on the computation of concentrated loads within the structural problem. In classic domain decomposition procedures like FEM, generalized forces are applied at a specific node of the domain mesh [85,86]. If the load is not applied in a computational point, a set of equivalent forces are derived by means of an interpolating procedure employing the adopted shape functions. The same approach can be found in references [87,88] where the GDQ algorithm has been adopted within each element in which the domain has been divided. As a consequence, a concentrated force pretends to be a boundary condition within the differential model. On the other hand, in the case of a problem developed within a single domain, in references [8991], an interesting procedure based on GDQ and GIQ methods accounts for a differential-integral implementation on a rectangular plate under a concentrated load, taking into account the main features of the well-known Dirac-Delta function [92]. On the other hand, a concentrated load can be seen as a particular case of a surface pressure acting on a very small area. Nevertheless, the distribution governing parameters should be set according to the actual dimension of the structure object of analysis [93,94].

In the present manuscript, a LW formulation is derived for the static analysis of generally anisotropic shell structures with a double curvature under the action of concentrated loads. The unknown displacement field variable is described, within each layer, employing a generalized approach with different higher order thickness functions. Furthermore, the displacement compatibility conditions between two adjacent layers are fulfilled. The geometry of the shell is described with curvilinear principal coordinates and a mapping procedure is adopted for arbitrarily-shaped structures. A generally anisotropic elastic behaviour is considered within each lamina. The fundamental set of differential equations is derived following an energy approach, together with the kinematic and static boundary conditions by means of a strong formulation of the structural problem. Then, non-conventional external constraints are enforced within the two-dimensional model. The differential problem is tackled numerically with the GDQ method, accounting for the discretization of the derivatives of a generic order, whereas integrals are solved with the GIQ numerical algorithm. Since the numerical model embeds in itself a smooth variation of derivatives, both surface distributed and concentrated loads are modelled with a comprehensive set of bivariate distributions, whose governing parameters have been carefully calibrated. Furthermore, the Dirac-Delta function is adopted in its classical and generalized GDQ discrete version for the assessment of concentrated loads [9597]. The main elements of novelty of the proposed method are the efficient numerical implementation of the concentrated load, which is a singularity among surface tractions, directly in the continuum model. Furthermore, the normalization of the distribution with respect to the surface area provides an effective calibration of the shape and position parameters. In the post-processing stage, an effective reconstruction of physical quantities throughout the entire shell thickness is assessed based on the three-dimensional static balance equations for a laminated anisotropic solid applied to each layer of the stacking sequence. A significative number of numerical investigations have been attached to the manuscript. The accuracy of the numerical predictions has been checked with respect to refined three-dimensional Finite Element models developed with a commercial package, showing a very good agreement between different approaches. Moreover, the inconsistency of higher order ESL formulations has been outlined for very thick structures characterized by very complex lamination schemes with a huge number of laminae with various material syngonies. Then, the solution has been checked for both single and double curvatures. The proposed higher order LW formulation has been added to the Differential Quadrature for Mechanics of Anisotropic Shells, Plates, Arches and Beams (DiQuMASPAB) project [98], a free research software which provides the static and the dynamic response of doubly-curved shell structures with various ESL and LW theories.

2  Doubly-Curved Shell Geometry

A doubly-curved shell is a three-dimensional solid within the Euclidean space (Fig. 1). For this reason, if e1,e2,e3 are the unit vectors of a global coordinate system, the position vector R(α1,α2,α3) of an arbitrary point of the structure can be described in terms of the following relation [21]:


Figure 1: Geometric assessment of a doubly-curved shell according to a LW approach. Representation of a generic thickness function of different orders defined in each k-th layer of the stacking sequence for k=1,,l, being l the total number of laminae occurring in the lamination scheme


where f1,f2,f3 are functions of the variables α1,α2,α3. It should be said that Eq. (1) assumes a physical meaning if the variations αi0αiαi1 for i=1,2,3 are declared, being αi0,αi1 the extremes of the variation intervals. If a laminated structure composed by l laminae is considered, the overall thickness of the structure h can be computed as the sum of the widths hk of each layer, with k=1,,l:


As a matter of fact, the association α3=ζ is performed so that the axis at issue is oriented alongside the thickness direction of the shell. Moreover, a reference surface r(α1,α2) is assessed, located in the middle thickness of the structure, whose parametric directions α1,α2 are defined from the principal geometric features of the shell. Referring to a generic k-th layer of a laminated structure, a unit vector set Oα1(k)α2(k)ζ(k) is introduced for k=1,,l. On the other hand, if hk assumes a constant value throughout the entire structure, α1(k)α2(k) can be obtained from an affine transformation alongside ζ direction of α1,α2 in-plane coordinates of the shell reference surface, thus setting αi(k)=αi for i=1,2. As a consequence, a key relation for the LW geometric and mechanic computation is introduced [24], defining the differential variation of local and global thickness coordinates ζ(k) and ζ:


In this way, a reference surface r(k)(α1,α2) is defined in the middle thickness of each k-th layer of the stacking sequence starting from the global geometric quantity r(α1,α2) referred to the global curvilinear coordinate Oα1α2ζ according to the following expression [21], as shown in Fig. 1.


where ζk,ζk+1 denote the extreme locations of the k-th layer along the shell thickness direction, whereas n(α1,α2) accounts for the normal unit vector of the reference surface r(α1,α2), defined as follows:


where r,i=r/αi denotes the partial derivative of the shell reference surface with respect to the already introduced principal direction αi for i=1,2 [21]. Starting from Eq. (4), the thickness curvature parameter Hi(k)(α1,α2) referred to the αi=α1,α2 principal direction is computed:


Furthermore, the Lamé parameters Ai(k)(α1,α2) and the principal curvature radii Ri(k)(α1,α2) of the reference surface of the k-th layer can be computed as:


where n is the normal unit vector defined in Eq. (5), whereas r,i(k)=r(k)/αi and r,ii(k)=2r(k)/αi2 account for the first and the second order derivative of Eq. (4) with respect to αi=α1,α2, respectively. In addition, Ai=Ai(α1,α2) and Ri=Ri(α1,α2) for i=1,2 are referred to the surface r(α1,α2) located in the middle thickness of the entire laminated structure. From the main outcomes of the differential geometry, such quantities are calculated according to the following expressions, setting r,ii=2r/αi2 [21]:


Having in mind all these premises, the three-dimensional position vector R(k)(α1,α2,ζ) of a generic point belonging to the k-th layer can be referred to r(α1,α2) as follows (Fig. 1):


being zk=2ζ(k)/hk a dimensionless thickness coordinate belonging to the interval [1,1]. It is useful to compute the first order derivative of zk with respect to the local thickness coordinate ζ(k), so that:


2.1 Arbitrarily-Shaped Shells

When the parametrization of the reference surface does not account for a curvilinear set of principal coordinates, the two-dimensional physical domain is distorted so that a rectangular dimensionless parent element is obtained, described in terms of the natural coordinates ξ1[1,+1] and ξ2[1,+1], as it has been schematically shown in Fig. 2. If we denote with (α1(i),α2(i)) the location of the i-th corner of the distorted geometry, for i=1,,4, the blending functions presented in the following can be employed [21]:


Figure 2: Isogeometric mapping of the physical domain employing NURBS curves. Definition of the local reference system along the edges of the distorted shell for the assessment of boundary conditions. Derivation of the rectangular computational domain employing natural coordinates

α1(ξ1,ξ2)=12((1ξ2)α¯1(1)(ξ1)+(1+ξ1)α¯1(2)(ξ2)+(1+ξ2)α¯1(3)(ξ1)+(1ξ1)α¯1(4)(ξ2))+ 14((1ξ1)(1ξ2)α1(1)+(1+ξ1)(1ξ2)α1(2)+(1+ξ1)(1+ξ2)α1(3)+(1ξ1)(1+ξ2)α1(4))(11)

 α2(ξ1,ξ2) =12((1ξ2)α¯2(1)(ξ1)+(1+ξ1)α¯2(2)(ξ2)+(1+ξ2)α¯2(3)(ξ1)+(1ξ1)α¯2(4)(ξ2))+ 14((1ξ1)(1ξ2)α2(1)+(1+ξ1)(1ξ2)α2(2)+(1+ξ1)(1+ξ2)α2(3)+(1ξ1)(1+ξ2)α2(4))(12)

In the previous equation, (α¯1(j),α¯2(j)) denotes the description in terms of α1,α2 of the j=1,,4 edge of the structure. To describe a generic curve C(u) with aub alongside the physical domain, a combination of n control points Pi for i=1,,n is employed, as follows (Fig. 2):


where Ri,p(u) is a rational B-Spline of p-th order which can be computed from the following expression, being wi a proper weighting coefficient:


Setting for simplicity [a,b]=[0,1] and introducing a predefined knot vector U=[a,,ap+1,up+1,,ump1,b,,bp+1] a recursive relationship can be adopted to compute the B-Spline basis function Ni,p(u) of p-th order. Starting with p=0, it is [21]:


When arbitrarily-shaped structures are investigated, the fundamental relations, provided in terms of (α1,α2)[α10,α11]×[α20,α21], should be expressed in terms of ξ1,ξ2 coordinates within the interval [1,1]×[1,1]. If we denote with J the Jacobian matrix of the coordinate transformation α1=α1(ξ1,ξ2), α2=α2(ξ1,ξ2), it can be stated that:


Accordingly, an inversion of Eq. (16) can be performed if det(J)=ξ1,α1ξ2,α2ξ1,α2ξ2,α10, leading to the definition of the inverse Jacobian matrix J1:


Once Eqs. (11) and (12) have been assessed, it is possible to provide an expression to define the first order partial derivatives with respect to α1,α2 principal coordinates in terms of the natural ones ξ1(α1,α2) and ξ2(α1,α2), starting from the well-known derivation chain rule:


being ξi,αj=ξi/αj for i,j=1,2. From a comparison of Eqs. (17) and (18), it is possible to provide the complete expression of ξi,αj coefficients, leading to [21]:


Based on Eq. (18), it is possible to express the second order derivatives with respect to the principal coordinate α1,α2 in terms of ξ1,ξ2. One gets:


Eq. (20) assesses the dependence of the second order derivatives with respect to α1,α2 in terms of the first and second order derivatives with respect to ξ1,ξ2 natural coordinates. Accordingly, the equation at issue is not bi-linear in the present formulation. In the following, the interested reader can find the complete expression of coefficients introduced in Eq. (20):




The first order derivatives det(J)ξ1,det(J)ξ2 of the determinant det(J) of the Jacobian matrix introduced in Eq. (16) with respect to ξ1 and ξ2, respectively, read as follows:


A useful nomenclature is now introduced, so that the edges of the physical domain can be univocally identified. Referring to the dimensionless rectangular parent element described in terms of the natural coordinates ξ1,ξ2, the following definitions [21] are outlined (Fig. 2):

West edge(W)ξ2=1(α1,α2)=(α¯1(1)(ξ1),α¯2(1)(ξ1))South edge(S)ξ1=1(α1,α2)=(α¯1(2)(ξ2),α¯2(2)(ξ2))East edge(E)ξ2=1(α1,α2)=(α¯1(3)(ξ1),α¯2(3)(ξ1))North edge(N)ξ1=1(α1,α2)=(α¯1(4)(ξ2),α¯2(4)(ξ2))(25)

3  Unified Formulations for Kinematic Relations

In the present section a unified assessment of the kinematic field variable is presented following the LW methodology. Referring to a generic k-th lamina of the laminate, each component U1(k)(α1,α2,ζ(k)) for i=1,,3 of the three-dimensional displacement field column vector U(k)(α1,α2,ζ(k)) can be expanded up to an arbitrary N-th order, thus introducing for each τ=0,,N+1 a generic function Fταi(k) for i=1,,3 dependent from the local thickness coordinate ζ(k) [24]:


where u1(kτ),u2(kτ),u3(kτ) are the generalized displacement field components defined for each τ-th kinematic expansion order lying on the r(k)(α1,α2) reference surface of the k-th lamina. In Fig. 1, one can find a graphic representation of Eq. (26).

Since a laminated doubly-curved shell structure is considered with k=1,,l, where l is the total number of superimposed laminae, the displacement field variable should fulfil the interlaminar compatibility conditions. To this purpose, an interpolation methodology based on higher order polynomials is followed for the definition of Fταi(k) for τ=0,,N+1. Accordingly, the dimensionless thickness coordinate zk introduced in Eq. (9) is considered:


where Lτ(zk) for τ=1,,N is defined employing various interpolating polynomials. As can be seen, Eq. (27) refers to the three-dimensional displacement distribution throughout the thickness. As a consequence, the interlaminar compatibility conditions are directly satisfied by the axiomatic assumptions of the unknown field variable itself. The first order derivative of the generalized thickness functions introduced in Eq. (27) with respect to ζ(k) can be computed as:


As can be seen, the derivatives zk/ζ(k) occurring in Eq. (28) are calculated according to Eq. (10).

If power functions are introduced in Eqs. (27) and (28), the following expressions should be adopted for each τ=1,,N [24]:


On the other hand, a formulation based on higher order Lagrange interpolating polynomials reads as follows:


If trigonometric functions are adopted in Eq. (27), Fταi(k)=Lτ(zk) for each τ=1,,N reads as:


Furthermore, the Jacobi orthogonal polynomials Jτ(γ,δ)(ζ(k)) can be employed to define the LW thickness function Fταi(k) for τ=1,,N, leading to:


where Jτ(γ,δ)(ζ(k)) of characteristic parameters γ,δ are calculated employing a recursive procedure:


setting A=2τ+γ+δ2, B=A2, C=A1 and D=2(τ1)(τ+γ+δ1)B. Accordingly, the first order derivative of Lτ(zk) of Eq. (32) with respect to ζ(k) occurring in Eq. (28) can be computed as:


Based on a LW higher order approach (26), the kinematic relations are derived starting from those referred to the three-dimensional solid described in Eq. (9), which are briefly recalled for the sake of completeness:


where ε(k)=[ε1(k)ε2(k)γ12(k)γ13(k)γ23(k)ε3(k)]T are the three-dimensional strain vector referred to the k-th layer. Furthermore, the kinematic through-the-thickness differential operator Dζ(k) is defined as:


whereas DΩ(k)αi for i=1,2,3 are defined so that they embed all the in-plane coordinates α1,α2 derivatives:


In Eq. (37) the differential vectors D¯Ω(k)αi have been introduced, whose extended version accounts as follows:


Accordingly, coefficients (D¯Ω(k)αi)j with j=1,,9 and αi=α1,α2,α3 read as:


The kinematic relation of the three-dimensional solid reported in Eq. (35), can be arranged if the unified LW displacement field assessment of Eq. (26) is substituted, leading to the introduction of the LW generalized strain vector ε(kτ)αi(α1,α2)=[ε1(kτ)αiε2(kτ)αiγ1(kτ)αiγ2(kτ)αiγ13(kτ)αiγ23(kτ)αiω13(kτ)αiω23(kτ)αiε3(kτ)αi]T, defined for each τ=0,,N+1 and αi=α1,α2,α3 [24]:


It is useful to introduce, for each αi=α1,α2,α3, the vector Z(kτ)αi referred to the τ-th kinematic expansion order, setting τ=0,,N+1:


Eventually, Eq. (40) turns into:


being ε(kτ)αi=DΩ(k)αiu(kτ).

4  Anisotropic Constitutive LW Relations

We now focus on the elastic constitutive behaviour of a generic doubly-curved laminated structure. Thus, each k-th layer of the stacking sequence, for k=1,,l, is modelled by means of a three-dimensional relationship valid for generally anisotropic materials. In this perspective, a local reference system denoted with Oα^1(k)α^2(k)ζ^(k) is derived from the direct application of the Neumann’s Principle to the periodic unit volume of each lamina. As a matter of fact, such material axes are intended to be featured so that one axis is parallel to the shell outward normal direction, namely ζ^(k)=ζ. If we denote with Eij(k) for i,j=1,,6 the generic stiffness constant linking the i-th component of the three-dimensional stress vector σ^(k)=[σ^1(k)σ^2(k)τ^12(k)τ^13(k)τ^23(k)σ^3(k)]Treferred to the material reference system to the corresponding j-th element of the three-dimensional strain vectorε^(k)=[ε^1(k)ε^2(k)γ^12(k)γ^13(k)γ^23(k)ε^3(k)]T, the generally anisotropic behaviour of the k-th layer can be expressed as [21]:


A key aspect of the present LW formulation is the assessment of all the fundamental governing relations into the geometric reference system of each layer oriented alongside the reference surface principal directions. To this purpose, an angle ϑk is identified in each k-th layer accounting for the deviation between α^1(k) and α1 material and geometric directions, respectively. If we denote with T(k)(ϑk) the rotation orthogonal matrix referred to each lamina of the structure, the constitutive relationship of Eq. (43) can be rotated with the following linear transformation so that it is referred to the geometric coordinate axes Oα1α2ζ of the reference surface of the structure thus leading to the rotated stiffness matrix E¯(k)=T(k)E(k)(T(k))T for each k=1,,l, whose generic component is denoted with E¯ij(k) for i,j=1,..,6:


being σ(k)=[σ1(k)σ2(k)τ12(k)τ13(k)τ23(k)σ3(k)]T and ε(k)=[ε1(k)ε2(k)γ12(k)γ13(k)γ23(k)ε3(k)]T the three-dimensional stress and strain vectors, respectively, referred to the shell geometric reference system. Referring to the anisotropic stiffness matrix E(k) occurring in Eq. (43), it usually consists in the three-dimensional elastic coefficients Eij(k)=Cij(k) for i,j=1,,6. In plane stress conditions (σ^3(k)=0) within the two-dimensional model in Eq. (43), the reduced elastic coefficients Eij(k)=Qij(k) are adopted. In particular, they are derived from a correction of the three-dimensional stiffness matrix, as follows:


It should be remarked that the constitutive relationship of Eq. (44), expressed for each k-th layer, has a three-dimensional connotation. As a matter of fact, it should be reduced to the local reference surface introduced in Eq. (4). From the computation of the variation δΦk of the elastic strain energy of the doubly-curved solid for each k=1,,l, one gets [21]:


Introducing in the previous relation the unified assessment of the displacement field variable of Eq. (26) and the LW kinematic relation of Eq. (42), for each τ-th kinematic expansion order with τ=0,,N+1 the generalized stress resultant vector S(kτ)αi(α1,α2)=[N1(kτ)αiN2(kτ)αiN12(kτ)αiN21(kτ)αiT1(kτ)αiT2(kτ)αiP1(kτ)αiP2(kτ)αiS3(kτ)αi]T is introduced, with αi=α1,α2,α3:


being A(kτη)αiαj the generalized constitutive operator, computed for each τ,η=0,,N+1 and αi=α1,α2,α3 according to the following definition:


In a more expanded form, A(kτη)αiαj matrix introduced in the previous equation reads, for each k-th layer, as:

A(kτη)αiαj=[ A(kτη)[00]αiαjA(kτη)[01]αiαjA(kτη)[10]αiαjA(kτη)[11]αiαj](49)

setting sub-matrices A(kτη)[00]αiαj,A(kτη)[01]αiαj,A(kτη)[10]αiαj,A(kτη)[11]αiαj as follows:





The generic component of Eqs. (50)(53) are obtained from a through-the-thickness homogenization of the mechanical properties of each k-th lamina according to the following expression, setting 0Fταi(k)/ζ(k)0=Fταi(k) and 0Fηαj(k)/ζ(k)0=Fηαj(k) [21]:


In the previous relation, coefficient B¯nm(k) is defined so that B¯nm(k)=E¯nm(k) for each n,m=1,..,6. Accordingly, if the LW definition of the displacement field of Eq. (26) accounts for a constant out-of-plane displacement field assumption, such quantity should be corrected by means of the well-known shear correction factor κ(ζ)=5/6, namely:


Referring to a generic τ=0,,N+1 and αi=α1,α2,α3, it is possible to express the higher order LW constitutive relationship of Eq. (47) in terms of the generalized displacement field vector u(kη)=[u1(kη)u2(kη)u3(kη)]T introduced in Eq. (26) for each reference surface of the k-th layer, taking into account the kinematic relation (40):


In Appendix A, an extended version of the components of the previously-introduced matrix O(kτη)αiαj can be found. It will be seen that the present version of the generalized elastic law is a key for the definition of both static and kinematic boundary conditions within the higher order LW framework.

5  Governing Equations

In the present section the fundamental relations of the static problem for a laminated doubly-curved structure are derived in the LW framework employing a higher order displacement field assumption. In particular, an energy approach will be followed, accounting for the curvature effects of the geometry. A generalized methodology is proposed for the assessment of surface loads acting on the structure, and an effective solution is provided for the implementation of concentrated loads. A consistent form of the fundamental governing equations is provided, together with the natural and non-conventional boundary conditions, characterized by three-dimensional capabilities.

5.1 External Loads

The present LW formulation considers a two-dimensional structural assessment for each layer of the laminated structure. Accordingly, a generic k-th lamina of constant thickness hk, with k=1,,l, is intended to be loaded at its intrados (ζ(k)=hk/2) by the static loads q¯1(k),q¯2(k),q¯3(k), along α1,α2,α3 principal directions, whereas the tractions q¯1(k+),q¯2(k+),q¯3(k+) are applied at the extrados (ζ(k)=+hk/2). For the sake of conciseness, the vector q¯(k±)=[q¯1(k±)q¯2(k±)q¯3(k±)]T is introduced. Thus, its components assume the following general form so that a general load case ϕ(k±)(α1,α2) can be assigned to the structure [21]:


Accordingly, if a uniform load is applied to the arbitrary k-th lamina, Eq. (57) is computed so that:


In addition, a Gaussian function has been implemented in the two-dimensional model, setting α1m(k),α2m(k) the position parameters, Δ1(k),Δ2(k)>0 the variances of the bivariate distribution and ρ12(k)[1,1] the correlation factor:


Furthermore, a Super-Elliptic shape of surface loads is introduced so that Δ1(k),Δ2(k) are the shape factors of the distribution of power coefficient n(k), whereas α1m(k),α2m(k) assume the role of position parameters and β(k) accounts for the orientation of the principal axes with dispersion:


As a matter of fact, for n(k)=2, Eq. (60) provides the well-known elliptic distribution. A general surface loading can be modelled also by means of a bivariate Fourier series in which n1(k)=n2(k) terms are assigned to the α1 and α2 direction, so that:


If the Jacobi polynomials Jl(γ(k),δ(k)) of governing parameters γ(k),δ(k) are employed for the assessment of general loads within each k-th lamina, the bivariate load distribution of Eq. (57) assumes the following form, setting k=1,,l:


with ϕi(k±)(αi), for i=1,2, reading as:


The complete expression of p01(k),p02(k),εl1(k),εl2(k) can be computed as:



Following the approach outlined in Eq. (57), a concentrated load can be embedded within the two-dimensional LW formulation by properly setting a load distribution. Let us consider for a generic k-th layer a concentrated load vector Q(k+)=[Q1(k+)Q2(k+)Q3(k+)]T of magnitude Q(k+) applied at ζ(k)=+hk/2 and a vector Q(k)=[Q1(k)Q2(k)Q3(k)]T referred to the layer bottom surface located at ζ(k)=hk/2 of magnitude Q(k). For the sake of conciseness, we will adopt the compact notation Q(k±)=[Q1(k±)Q2(k±)Q3(k±)]T. Accordingly, the deviation of the vector at issue from α1,α2,ζ shell principal directions is identified by means of the angles φ1(k±),φ2(k±),φ3(k±), respectively. As a matter of fact, each component Qi(k±) of Q(k±), is derived as follows [21]:


Thus, an equivalent surface traction q~i(k±) is provided for each i=1,2,3, defining the effects of the concentrated load components of Eq. (66):


The surface tractions of Eq. (67) are conveniently arranged in the vector q~(k±)=[q~1(k±)q~2(k±)q~3(k±)]T. Generally speaking, the bivariate function ϕ~(k±)(α1,α2) occurring in the previous equation consists of the well-known Dirac-Delta function δ, namely:


The function at issue has a singularity for (α1m(k),α2m(k)), accounting for the following properties:


On the other hand, the numerical modelling of the δ function in the continuum smooth model can be quite difficult because no closed-form analytical expressions can be provided in Eq. (69). As a consequence, the physical meaning of concentrated load could be not properly interpreted. For this reason, the concentrated load is modelled as a particular case of surface load according to Eq. (57), as here applied to a very small area. Furthermore, the distribution is normalized with respect to the area of the extrados or intrados of the shell so that it perfectly fulfils the main properties of the Dirac-Delta function in Eq. (69).

When an arbitrary bivariate distribution ϕ(k±)(α1,α2) is selected for the description of a concentrated load among those reported in Eqs. (58)(62), the following relation should be considered:


where Iϕ(k±) denotes the surface integral of ϕ(k±)(α1,α2) performed on the top (+) or bottom () surface of the k-th layer calculated by means of the rectangular domain [α10,α11]×[α20,α21]:


In this way, the surface integral of ϕ~(k±)(α1,α2) distribution employed in Eq. (67) fits the main features of the Dirac-Delta function already outlined in Eq. (69), thus giving:


In other words, Eq. (72) requires that the resultant of the corresponding pressure associated to the concentrated load should be equal to the three-dimensional applied vector Q(k±) in all its components. A surface pressure q(k±)=[q1(k±)q2(k±)q3(k±)]T is introduced for each αi=α1,α2,α3 principal direction, consisting in a contribution referred to the generally-shaped load q¯(k±)=[q¯1(k±)q¯2(k±)q¯3(k±)]T assessed in Eq. (57) and the vector q~(k±)=[q~1(k±)q~2(k±)q~3(k±)]T embedding the effects of concentrated loads:


Since a higher order displacement field assumption has been considered according to Eq. (26) in each layer of the structure, the Static Equivalence Principle is applied so that the computation of the virtual work associated to the vectors q(k+)=[q1(k+)q2(k+)q3(k+)]T and q(k+)=[q1(k)q2(k)q3(k)]T gets into the derivation of a generalized surface load vector q(kτ)=[q1(kτ)q2(kτ)q3(kτ)]T defined within the reference surface of all the layers of the stacking sequence for k=1,,l. One gets:


In the previous equation, Fταi(k±) refers to the computation of the thickness function Fταi(k) at the top and bottom, respectively, of the k-th layer. In the same way, the main curvature parameters H1(k±) and H2(k±) are introduced within the generic lamina. As a matter of fact, in all the simulations reported in the present manuscript, a perfect bonding between two adjacent laminae is assumed, therefore the structure can be loaded, according to Eq. (74), only at its top and bottom surfaces, located at ζ=+h/2 and ζ=h/2, respectively.

5.2 Fundamental Relations

We now account for the energy procedure employing the well-known Minimum Potential Energy Principle [21] for the determination of the static response of the structure under the action of static loads. According to the LW approach, a stationary configuration of the potential energy Π of the shell is computed, taking into account the variation of the total elastic strain energy δΦ and virtual external work δLe:


Integrating by parts, the following equilibrium relations are derived in terms of S(kτ)αi and q(kτ), reading as [21]:

i=13DΩ(k)αiS(kτ)αi+q(kτ)=0for τ=0,,N+1,k=1,,l(76)

The differential operators DΩ(k)αi=DΩ(k)α1,DΩ(k)α2,DΩ(k)α3 are defined for each k-th layer as:


setting the following definitions:


An extended version of the components of the vectors introduced in Eq. (78) has been now reported:


The fundamental relations for the static assessment of an anisotropic doubly-curved shell in terms of generalized displacement components u1(kτ),u2(kτ),u3(kτ) are easily derived from Eq. (76) combined with Eq. (42) and (47), leading to an expression referred to a generic τ-th kinematic expansion order [21]:


where the fundamental operator L(kτη) is defined for each τ,η-th in a generic k-th layer with k=1,,l as follows:


The complete expression of Lij(kτη)αiαj with τ,η=0,,N+1 and i,j=1,2,3 can be found in Appendix B.

Starting from the physical interpretation of the kinematic variables introduced in Eq. (27), it should be recalled that the generalized displacement field vectors corresponding to τ=0 and τ=N+1 are defined in such a way that the compatibility conditions that should be enforced between two adjacent layers are implicitly enforced, namely:


being u(k(N+1)) and u((k+1)0) the generalized displacement field vectors associated to the extrados and intrados of the k-th and the (k+1)-th layer, respectively. Furthermore, Eq. (80) is assembled so that all the expansion orders of the kinematic assumption in Eq. (26) are considered, leading to the final fundamental relation:


Note that the interlaminar compatibility conditions of Eq. (82) are modelled in Eq. (83) by means of the identity matrix I located in the proper position of the fundamental operator. In this way, an independent equation is added for each k-th layer so that the generalized displacement field associated to τ=0 is the same to the (k+1)-th lamina, referred to a kinematic expansion order τ=N+1.

From the present energy formulation it is also possible to derive the conventional external constraints associated to the boundaries of the physical domain. If a generalized displacement field component is assigned within each k-th layer of the laminate, the following relations should be enforced at the shell edges [21]:


where u¯1(kτ),u¯2(kτ),u¯3(kτ) are the prescribed values of the generalized displacement components. On the other hand, a prescribed set of generalized stress resultants leads to the definition of the following boundary conditions within each k-th layer, for k=1,,l [21]:


with N¯1(kτ),N¯12(kτ),T¯1(kτ) enforced at α1=α1i with i=0,1, whereas N¯21(kτ),N¯2(kτ),T¯2(kτ) are referred to α2=α2i for i=0,1. As a particular case of what exerted in Eq. (84), a fully clamped (C) configuration is outlined when all the generalized displacement field components are neglected for each k=1,,l and τ=0,,N+1:


In a similar way, the free (F) edge boundary condition moves from Eq. (85) as follows:


Referring to Eq. (85), the generalized stress resultants employed for the assessment of boundary conditions are enforced if the components of a boundary stress vector σ¯(k)=[σ¯1(k)σ¯2(k)τ¯12(k)τ¯13(k)τ¯23(k)σ¯3(k)]T are applied at the edges of the structure. Accordingly, σ¯(k) is intended to be obtained from the sum of a prescribed stress vector σ~(k)=[σ~1(k)σ~2(k)τ~12(k)τ~13(k)τ~23(k)σ~3(k)]T and a vector σ(k)=[σ1(k)σ2(k)τ12(k)τ13(k)τ23(k)σ3(k)]T dependent from the three-dimensional displacement field vector U(k)=[U1(k)U2(k)U3(k)]T:


Referring to the shell sides located at α1=α1s for s=0,1 within the physical domain, the following definitions can be outlined in each k-th layer of the stacking sequence [21]:


where λ¯=λ¯(ζ) accounts for a generalized through-the-thickness normalized distribution of stress components, whereas β(α1s,α2) denotes the in-plane distribution of prescribed stresses. In the same way, the static boundary conditions are applied at α2=α2s for s=0,1 and α1[α10,α11] at each k-th lamina, with k=1,,l, according to the following expressions:


being β(α1,α2s) the axiomatic assumed in-plane stress distribution. In the present formulation, a general distribution λ¯ of applied stresses has been considered for each shell edge acting along the thickness h of the structure according to Eq. (2). For instance, a constant dispersion has been modelled so that λ¯=1. Moreover, a linear dispersion can be assessed as:


Furthermore, a parabolic stress distribution accounts for a polynomial expression of λ¯:


Starting from Eqs. (89) and (90), a generalized set of non-conventional boundary conditions and prescribed stresses is developed if general in-plane univariate distributions β(α1s,α2)=βs(α2) and β(α1,α2s)=βs(α1) with s=0,1, are associated to the components of the applied stresses vector σ¯(k) for k=1,,l. To this purpose, a dimensionless coordinate ξ¯r with r=1,2 is introduced within the closed interval [0,1] for a smart assessment of general boundary conditions, namely:


In the case of constant in-plane distribution of stresses, the relation βs(αr)=1 for r=1,2 is assumed. Furthermore, two different analytical univariate expressions for βs(αr)=βs((αr1αr0)ξ¯r)=βs(ξ¯r) have been provided. A Double–Weibull (W) distribution accounts as follows:


where ξ~r=1ξ¯r, whereas ξ¯m,ξ~m[0,1] and p refer to the position and shape parameters, respectively. In addition, a Super-Elliptic (S) dispersion has been modelled according to the following expression:


By the way, based on Eqs. (89) and (90) a set of non-conventional constraints can be enforced to the model if the stress components of σ(k) are provided in each point of the three-dimensional shell edge by a series of linear elastic springs. Referring to the regions located at α1=α1s for s=0,1, it can be said that [21]:


where k1f(k)α1s,k2f(k)α1s and k3f(k)α1s are the stiffness of the springs in the α1,α2,α3 directions, respectively. Accordingly, the corresponding relations for α2=α2s for s=0,1 read as:


being U1(k),U2(k),U3(k) the three-dimensional displacement field associated to the k-th layer with respect to the geometric principal reference system. Eqs. (96) and (97) are, thus, embedded in the LW framework if the unified assessment of the displacement field of Eq. (26) with higher order theories is taken into account. Eq. (96) turns into:


whereas Eq. (97) assumes the following form:


As a consequence, for each τ-th kinematic expansion order with τ=0,,N+1, two different sets of generalized stress resultants are derived which are associated to σ~(k) and σ(k) stress components, according to the following definitions:


The generalized boundary stress resultants acting at α2=α2s for s=0,1 read as follows:


Referring to Eq. (100), the contribution coming from the applied stresses and linear springs can be arranged in the following compact matrix form:


In the same way, the following relation for α2=α2s can be assessed starting from Eq. (101):


Fundamental coefficients Li(p)αnsfb(kτη)αi occurring in Eqs. (102) and (103) are computed for each k-th layer and τ,η=0,,N+1, according to the following effective expression:


In the case of arbitrarily-shaped domains, the application of the blending functions of Eqs. (11) and (12) requires a rearrangement of natural boundary conditions assessed in Eqs. (86) and (87). To this purpose, a right-handed reference system is introduced from the geometric properties of a generic curve lying on the r(α1,α2) reference surface. The corresponding unit vectors, denoted with nn, ns and nζ, read as [21]:


where nri for i=1,2,3 and with r=n,s,ζ are the components of the vectors at issue with respect to the shell geometric reference system Oα1α2ζ. In particular, since nr with r=n,s,ζ stands for the local principal directions of an arbitrary curve belonging to the reference surface r(α1,α2), it should be said that nn3=ns3=nζ1=nζ2=0 and nζ3=1. Referring to a particular τ-th kinematic expansion order, for τ=0,,N+1, the generalized displacement field vector u(τ) can be expressed with respect to such coordinate system according to the following transformation relation:


being un(kτ),us(kτ),uζ(kτ) the components of u(τ) referred to nn, ns and nζ directions, respectively. In the same way, static boundary conditions can be enforced on a distorted domain in terms of generalized higher order stresses Nn(kτ)α1, Nns(kτ)α2 and Tζ(kτ)α3 referred to the local coordinate system at issue:


For prescribed displacements u¯n(kτ),u¯s(kτ),u¯ζ(kτ) or stresses N¯n(kτ)α1,N¯ns(kτ)α2,T¯ζ(kτ)α3 alongside the edges of an arbitrarily-shaped shell, the mechanical and kinematic constraints reported in the following are derived from the minimum potential energy principle [21]:


6  Equivalent Single Layer Theory

In the present manuscript a LW formulation is presented for laminated anisotropic doubly-curved shells. Since a generalized approach has been followed, the structure can be geometrically described in terms of r(α1,α2) referring to the global curvilinear coordinate system Oα1,α2ζ assessed in Eq. (1). As a consequence, the three-dimensional position vector can be expressed as [21]:


where r(α1,α2) is located in the middle thickness of the entire structure. The geometric Lamè parameters A1(α1,α2),A2(α1,α2), of the shell, as well as the main curvature radii R1(α1,α2), R2(α1,α2), are thus calculated by means of Eq. (8).

As far as the unified formulation of the displacement field is concerned, a set of u(τ) generalized vectors for τ=0,,N+1 is defined on the reference surface r(α1,α2), thus turning Eq. (26) into the following one:


The thickness function matrix F(τ) is defined employing a power expansion for the displacement field. In the case of laminated structures, the Murakami’s function is adopted. For more details on the topic, the interested reader can refer to reference [24]. For concentrated and surface loads within the ESL formulation, the generalized distributions of Eqs. (58)(62) become independent from the k-th lamina.

7  Numerical Implementation with the GDQ Method

In the present section the LW model for anisotropic doubly-curved shells outlined in the manuscript is numerically tackled by means of the GDQ method. Belonging to the class of spectral collocation algorithms, the GDQ approach represents a quadrature procedure to discretize the n-th order derivatives of an arbitrary function. Referring to a generic univariate function f=f(x) with x[a,b], it gives [21]:


where N>n due to the consequences of the Weierstrass Interpolation Theorem. The weighting coefficients occurring in Eq. (111) are calculated from the following recursive procedure:


Accordingly, in the present manuscript the computational domain has been discretized in IN and IM points along α1,α2 principal directions, respectively, according to the non-uniform Chebyshev-Gauss-Lobatto (CGL) harmonic distribution [21]. Referring to a dimensionless domain [1,1], the generic point xi of the distribution with i=1,IQ is introduced so that:


Starting from the GDQ rule in Eq. (111), the well-known GIQ method is assessed so that integrations restricted to a generic interval [xi,xj] of a univariate function f=f(x) can be numerically tackled, setting xk discrete points, with k=1,,l:


GIQ weighting coefficients wkij=wjkwik are computed following the procedure reported in reference [21]. We now focus on the numerical assessment of the concentrated loads within the computational domain by means of the Dirac-Delta function according to what exerted in Eq. (69). Following the procedure suggested by references [90,91], the GDQ and GIQ rules are adopted for the implementation of the concentrated load in the present strong form problem according to Eq. (69), assuming that the vector at issue is applied in one of the selected discrete computational points. The GDQ algorithm of Eq. (111) for the discretization of the fundamental differential relations of Eq. (83) is adopted for all the internal discrete points of the computational domain, whereas the static and kinematic external constraints are numerically enforced at boundaries. Referring to the inner nodes of the IN×IM two-dimensional grid, it gives for r=2,,IN1 and s=2,,IM1:


where A1(rs)(k),A2(rs)(k) are evaluated in each point of the computational grid according to Eq. (7), whereas the integral properties of the Dirac-Delta function of Eq. (69) are then adopted for the discretization of the fundamental governing equations in the computational point corresponding to that of the physical domain, denoted with (α1m(k),α2m(k)), where the concentrated load has been applied according to Eq. (66). Some remarks are reported in references [85,90] for more details.

Furthermore, the Dirac-Delta function of Eq. (68) has been also implemented according to the generalized approach in reference [89]. Moving from the methodology presented in Eq. (115), the concentrated load application point (α1m(k),α2m(k)) can be selected regardless the nature of the computational grid. In particular, a procedure based on the Lagrange interpolating Polynomials is followed so that the applied load is transferred to the IN×IM discrete set of point starting from an arbitrary location within the physical domain. One gets for r=2,,IN1 and s=2,,IM1 [90]:




Once the fundamental differential problem outlined in Eq. (83) has been implemented by means of the GDQ method according to Eq. (111), the overcoming computational problem is efficiently solved by means of a proper condensation of the linear system. To this purpose, the unknown DOFs, embedded in the vector δ, are arranged so that δb and δd are referred to the boundary (“b” points) and the inner nodes (“d” points), respectively [21]. One gets:


where fb,fd are the external load vectors associated to the “b” and “d” points, respectively. If Eq. (118) is expressed only in terms of δd vector, the following reduced linear system is obtained:


8  Post-Processing

The present formulation accounts the static structural assessment of laminated doubly-curved shell structures employing a two-dimensional formulation by LW approach. Since the solution is located at the middle surface of each k-th layer, the higher order displacement field assumption of Eq. (26) can be adopted for the derivation of the three-dimensional response of the solid. On the other hand, the results may not fulfil the external tractions applied at the intrados and the extrados of each lamina. For this reason, a correction of stresses should be performed.

From the closed interval [hk/2,hk/2], representing the k-th lamina thickness, a set of IT points is selected, whose generic one ζ(g)(k) with g=1,,IT is derived from the CGL harmonic distribution of Eq. (113). Then, the three-dimensional displacement field vector is U(ijg)(k) evaluated in each (α1(i),α2(j),ζ(g)(k)) point of the three-dimensional solid, for i=1,,IN, j=1,,IM and g=1,,IT, employing the unified approach of Eq. (26):


In the same way, the discrete form ε(ijg)(k)=[ε1(ijg)(k)ε2(ijg)(k)γ12(ijg)(k)γ13(ijg)(k)γ23(ijg)(k)ε3(ijg)(k)]T of the three-dimensional strain vector is calculated according to Eq. (35), setting i=1,,IN, j=1,,IM and g=1,,IT:


Starting from the previous equation, it is possible to derive in the arbitrary point located at (α1(i),α2(j),ζ(g)(k)) the corresponding membrane stresses σ1(ijg)(k),σ2(ijg)(k),τ12(ijg)(k), for each k-th layer, according to the elastic constitutive law of Eq. (44), leading to:


At this point, the three-dimensional equilibrium equations of a doubly-curved solid written in curvilinear principal coordinates should be recalled, remembering that dζ=dζ(k):


with a(k),b(k),c(k) reading as:


From the three-dimensional relations reported in Eq. (123), the out-of-plane stress components τ13(ijg)(k),τ23(ijg)(k) are derived for each (α1(i),α2(j),ζ(g)(k)) point, once the in-plane stresses are computed in Eq. (122), and their first order derivatives are calculated by means of the GDQ method of Eq. (111). Furthermore, two different loading conditions have been considered for each k-th lamina, leading to prescribed values of τ13(k),τ23(k) at the top and the bottom surfaces of the layer. Nevertheless, only the tractions referred to ζ(k)=+hk/2 can be enforced. With particular reference to the first lamina k=1, if q1(ij)(1),q2(ij)(1) denote the in-plane components of the load vector already defined in Eq. (73) for each (α1(i),α2(j)), applied at ζ(1)=h1/2, one gets:


Since a perfect bonding has been considered at the interlaminar level, the following relations should be considered too:


The values of τ13(k),τ23(k) shear stresses obtained from Eq. (123) can now be corrected so that the in-plane loading conditions referred to the l-th lamina are fulfilled, namely τ13(ijIT)(l)=q1(ij)(l+) and τ23(ijIT)(l)=q2(ij)(l+). To this end, two vectors τ~13(ij),τ~23(ij) of IL=lIT components are introduced at each (α1(i),α2(j)) point with i=1,,IN and j=1,,IM, setting s=kg:


Each element of the vectors at issue can now be corrected so that out-of-plane shear stresses fulfil the in-plane loading conditions at the top surface of the shell [21]:


where τ~13(ijg)(k)=τ~13(ijs) and τ~23(ijg)(k)=τ~23(ijs). Taking into account the recovered stress of Eq. (128), from the third relation of Eq. (123) the actual dispersion of the pressure along ζ direction is provided, enforcing the compatibility conditions between two adjacent laminae σ~3(ijIT)(k)=σ~3(ij1)(k+1) for k=1,,l1, together with the loading conditions at the intrados of the structure:


In the same way of what defined in Eq. (127), the vector σ~3 of normal stresses is introduced:


If the notation σ~3(ijg)(l)=σ~3(ijs) is adopted with s=kg=1,,IL=lIT, the out-of-plane normal stress satisfies the load boundary condition enforced at the shell top surface if the following correction of σ~3 components is performed [21]:


Furthermore, the out-of-plane three-dimensional strain components profile can now be adjusted employing the recovered distributions of out-of-plane stresses τ~13(ijs),τ~23(ijs),σ~3(ijs) for s=1,,IL. The following relations are adopted for each k-th layer:


where detA(ijg)(k) reads as:


The complete procedure for the derivation of the recovered γ13(ijg)(k),γ23(ijg)(k),ε3(ijg)(k) strain profiles is explained in reference [21].

9  Applications and Results

In the present section the LW formulation presented in the manuscript is applied to some structures of different curvature and materials, subjected to various load cases. In particular, the advantages of the present formulation compared with other trustworthy models are outlined. At a first stage, a fully-clamped beam subjected to a central concentrated load is considered, and the governing parameters of the load distributions presented in Eqs. (58)(70) are calibrated with respect to a closed-form analytical solution. After that, the static deflection of some thick shells characterized by zero, single and double curvature are calculated under different kinds of loads and boundary conditions. Furthermore, different kinds of lamination schemes have been considered, accounting for various numbers of laminae with both softcore and hardcore behaviours. In this way, the inconsistency of higher order ESL theories is outlined when such structures are employed, whereas the proposed higher order LW formulation provides very good performances with respect to three-dimensional refined solutions.

9.1 Validation for Concentrated Load Distributions

Let us consider a beam of length Lx=10m characterized by a rectangular cross section of dimensions b=1m and h=0.1m made of isotropic Aluminium (E=70GPa,ν=0.3,ρ=2707kg/m3). The structure is clamped at its extremities, and it is subjected to a concentrated vertical load Q3(+)=1000N applied at its mid-span (Fig. 3). A reference value for the maximum deflection of the structure has been calculated from the well-known Euler-Bernoulli Theory (EBT) according to the following expression:


Figure 3: Calibration of the parameters of the distribution employed for the assessment of concentrated loads on the two-dimensional physical domain starting from a clamped beam subjected to a concentrated force at the middle span. A reference value of the maximum deflection has been calculated with the EBT closed-form solution. Accuracy of the numerical assessment of the Dirac-Delta function and the Generalized Dirac-Delta function (a). Accuracy of the results obtained when the Gaussian and the Super-Elliptic distribution are employed (b). Sensitivity analysis of a Jacobi polynomials-based modelling of the concentrated load (c). In (d), the precision of a two-dimensional model employing various numbers of Fourier series terms is outlined


being I=bh3/12 the moment of inertia of the rectangular cross section. The same structure has been analysed with the GDQ method by means of Eq. (83) employing in Eq. (110) the FSDT kinematic field assumptions. In this way, the structure has been investigated based on the ESL two-dimensional approach. For each investigation, the percentage error e% of the solution with respect to wmaxEBT has been plotted, computed by means of the following expression:


As a matter of fact, in this case the boundary conditions are assigned so that a FCFC configuration is obtained according to Eqs. (86) and (87), following the nomenclature of Eq. (25). Thus, a parametric study has been performed so that the sensitivity of the main governing parameters of the concentrated load distributions presented in this manuscript is shown.

In Fig. 3a, the Dirac-Delta function of Eq. (68) has been adopted for the two-dimensional simulations according to the GDQ approach presented in Eq. (115). Moreover, the generalized version of the dispersion at issue has been studied employing the reduction to the computational nodes of the applied load according to the procedure of Eq. (116). Setting IM=31 for the numerical assessment of the two-dimensional model, the number IN of points along the beam length has been varied, showing that an increased grid dimension leads to more accurate results. More specifically, when a specific point of the computational grid is located in the middle span of the structure, i.e., an even number IN is adopted, a percentage error (135) lower than 1% is obtained with a very reduced number of points. On the other hand, when a higher order interpolation procedure is required, at least IN=24 discrete points are required to obtain stable and accurate results. This is due to the fact that the procedure based on a higher order interpolation of the Dirac-Delta function performs a reduction of the applied load to the adopted grid points, therefore a small accuracy loss is noticed. On the other hand, the GDQ-based integral-differential procedure for the numerical implementation of the Dirac-Delta function do not require any interpolation since the load is applied directly at the grid points.

In Figs. 3b3d, the sensitivity of the continuous distribution parameters is checked, setting IN=IM=31. In particular, it has been shown that the Gaussian distribution of Eq. (59) provides a very good agreement with analytical solutions if the variance parameters Δ1(k)=Δ2(k) are set equal to 2%. In the case of Super-Elliptic distributions (60) of various power exponents (Fig. 3b), an increase of n(k) does not lead to an improvement of precision of the simulation. In any case, stable results are reached if Δ1(k)=Δ2(k)=3%. This means that concentrated loads can be effectively described without any loss of accuracy if an equivalent pressure is applied to an area with a radius smaller than 3% of one edge of the physical domain. The efficiency of the formulation for concentrated loads by means of the Jacobi polynomials of Eq. (63) is shown in Fig. 3c, where the percentage error of Eq. (135) is shown for different γ(k) and δ(k) so that various polynomials belonging to the class at issue are employed. A sensitivity analysis with respect to n1(k)=n2(k) is shown in Fig. 3c. In particular, it is shown that for γ(k)=δ(k)=0.5 stable results with negligible errors are reached for n1(k)=n2(k)=70, whereas other polynomials require higher order polynomials. When γ(k)=δ(k)=1, the proposed formulation is not capable of providing good results in any case. When the concentrated loads are implemented by means of a Fourier series expansion according to Eq. (61), at least 300 terms are required in the truncated series, as it has been shown in Fig. 3d, telling that the singularity can be properly modelled with the superimposition of at least 300 sinusoidal functions.

Once the distribution of the main governing parameters are checked in Fig. 3, the percentage error of Eq. (135) has been computed with respect to the discretization of the physical domain for all the continuous distributions presented in this manuscript (Fig. 4). As it has been shown in the previous simulations, the employment of the Dirac-Delta function in its discrete form for the assessment of concentrated loads provides an excellent agreement with closed-form solutions even though a significative reduced number of DOFs are employed in the model. Among continuous functions, if the Super-Elliptic distribution of Eq. (60) is adopted, a fast stabilization of results with an excellent level of accuracy is seen regardless the selection of the distribution governing parameter. For this case, a constant value δ=0.02 has been selected, whereas the order of the distribution has been varied from n(k)=2 to n(k)=10. In a similar way, the Gaussian distribution of Eq. (59), corrected according to Eq. (70), provides very good results even with IN=11 grid points. The employment of the Fourier function of Eq. (61) with n(k)=315 terms rapidly leads to reduced values of e% even though an oscillation is seen by varying the dimension IN of the discrete computational grid. Furthermore, for a Jacobi distribution, high values of IN are required to obtain a stable behaviour. The sensitivity of the computational grid has been checked for two different Jacobi polynomials of order n(k)=71, namely the Legendre polynomials (γ(k)=δ(k)=0), and the Chebyshev polynomials of I kind, characterized by γ(k)=δ(k)=0.5. As can be seen, the best agreement with respect to the closed-form solution of Eq. (134) is achieved when the Legendre polynomials are employed.


Figure 4: Parametric investigation of a clamped beam subjected to a concentrated load applied at the middle span of the structure. The discrepancy error has been calculated between the maximum deflection provided by the closed-form EBT and the numerical implementation proposed in the manuscript. The sensitivity of the computational grid has been outlined with respect to each load distribution

9.2 Validation of the LW Formulation

Once the governing parameters of the load distributions have been validated, the proposed LW methodology has been employed for the analysis of some structures of different features subjected to various loads and non-conventional boundary conditions. Accordingly, the examples of investigations have been selected so that the main advantages of the LW solutions are checked with respect to more simplified two-dimensional ESL methodologies.

In all simulations presented in this section, concentrated loads have been modelled by means of the Dirac-Delta function in Eq. (68), according to the numerical GDQ procedure. Furthermore, the three-dimensional static response of the structures at issue has been calculated from refined three-dimensional simulations employing a classical 3D FEM. In this way, a reference solution is provided for validation purposes of the proposed methodology. The results have been provided in terms of the through-the-thickness dispersion of the three-dimensional kinematic and mechanical quantities in some points of the structure.

Different laminated structures have been considered, accounting for various stacking sequences with both hardcore and softcore behaviours. To this end, different classes of materials have been employed, i.e., anisotropic, orthotropic and isotropic materials. For the first class we consider a Triclinic material [95], whose stiffness matrix is defined with respect to the material symmetry planes according to conventions Eq. (43):


For the sake of completeness, the density of the material has been taken equal to ρ(k)=7750kg/m3. In order to provide a generally anisotropic material with softcore features, the so-called Triclinic-Soft material has been derived from Eq. (136), setting each stiffness constants equal to 1/1000 of the original one.

Referring to the group of orthotropic materials, the material stiffness properties have been expressed employing the well-known engineering constants, namely the three elastic moduli E1(k),E2(k),E3(k), the shear moduli G12(k),G13(k),G23(k) and the Poisson coefficients ν12(k),ν13(k),ν23(k). The relationship between the quantities at issue and the three-dimensional stiffness constants Eij(k) for i,j=1,,6 can be found in reference [21]. In the following, the material properties of an orthotopic Carbon Fibre Reinforced Polymer (CFRP) of density ρ(k)=1824kg/m3 have been collected [95]:


The orthotropic Graphite-Carbon Epoxy (ρ(k)=1760kg/m3) reads as follows:


In addition, the following lattice material named 3D Augmented Re-entrant Cellular Structure (3D ARCS), with density ρ(k)=66.468kg/m3, has been considered:


Among the isotropic materials, a foam of density ρ(k)