Computer Modeling in Engineering & Sciences

---

DOI: 10.32604/cmes.2022.022210

ARTICLE

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

Francesco Tornabene*, Matteo Viscoti and Rossana Dimitri

Department of Innovation Engineering, School of Engineering, University of Salento, Lecce, 73100, Italy
*Corresponding Author: Francesco Tornabene. Email: francesco.tornabene@unisalento.it
Received: 28 February 2022; Accepted: 11 May 2022

1  Introduction

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

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

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

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

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

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

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

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

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

Once the 2D solution is provided on the reference surface of a doubly-curved shell, the 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) [116]. As far as the through-the thickness equilibrium-based recovery procedure is concerned, the generalized constitutive relationship can be adopted as usual, but a correction of shear and membrane stresses is essential so the boundary conditions related to applied loads can be fulfilled [117–120].

In the present work, an ESL theoretical formulation is proposed for generally anisotropic 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 [121]. A general distribution of surface loads has been considered in the model, setting a Gaussian and a Super-Elliptic distribution for the application of external loads along each principal direction of the shell. In particular, the calibration of distribution parameters leads to model different load cases. The fundamental governing equations, derived from the Minimum Potential Energy Principle, are numerically tackled by the GDQ method. The Generalized Integral Quadrature (GIQ) method has been adopted for the computation of integrals [16]. The GIQ moves from the above-discussed GDQ procedure, and it allows performing numerical integrations following an effective quadrature approach.

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

2  Geometrical Description of a Doubly-Curved Shell

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

R(α1,α2,α3)=f1(α1,α2,α3)e1+f2(α1,α2,α3)e2+f3(α1,α2,α3)e3(1)

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

R(α1,α2,ζ)=r(α1,α2)+ζn(α1,α2)(2)

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

n=r,1×r,2|r,1×r,2|(3)

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 R(α1,α2,ζ)∈[α10,α10]×[α20,α21]×[−h/h22,h/h22]. Accordingly, shell thickness h(α1,α2) can be defined starting from Eq. (2), according to the following expression taken from reference [16]:

h(α1,α2)=∑k=1lhk(α1,α2)=∑k=1l(ζk+1(α1,α2)−ζk(α1,α2))(4)

where hk(α1,α2)=ζk+1(α1,α2)−ζk(α1,α2) is the width of the k-th lamina of the stacking sequence composed of l∈N layers located within the interval [ζk,ζk+1]⊆[−h/h22,h/h22]. If the second order derivatives r,ij=∂2r/∂2r(∂αi∂αj)(∂αi∂αj) of the reference surface with αi,αj=α1,α2 are calculated, it is possible to compute the principal radii of curvature Ri=R1,R2 in each point of the physical domain:

Ri=−r,i⋅r,ir,ii⋅nfori=1,2(5)

In addition, the well-known Lamè parameters Ai=A1,A2 are defined as follows:

Ai=r,i⋅r,ifori=1,2(6)

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

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

Hi(α1,α2,ζ)=1+ζRi(α1,α2)fori=1,2(7)

Accordingly, in the case of straight shells along the αi principal direction, Eq. (7) turns into Hi(α1,α2,ζ)=1.

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

hk(α1,α2)=hk0h~(α1,α2)=hk0(1+∑j4δjϕj(α1,α2)+δ¯)(8)

being hk0 the reference height of the k-th layer and δj for j=1,…,4 a scaling parameter, whereas δ¯ is a constant dimensionless shift number. In Eq. (8), thickness variation h~(α1,α2) has been described in terms of four analytical expressions ϕj(α1,α2) employing a dimensionless coordinate defined along each αi=α1,α2 principal direction of the shell, according to the following equation:

α¯i=αi−αi0αi1−αi0fori=1,2(9)

In the present formulation, the expressions of ϕj(α1,α2) for j=1,…,4 have been implemented, being αjm∈[0,1], nj∈N and pj∈R [16]:

ϕj(α1,α2)={α¯ipje−(α¯iαjm)pj(sin⁡(π(njα¯i+αjm)))p1fori=1,2j=1,2ϕj(α1,α2)={(1−α¯i)pje−(1−α¯iαjm)pj(sin⁡(π(nj(1−α¯i)+αjm)))pjfori=1,2j=3,4(10)

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

n n = [ n n1 ​​  n n2   n n3 ] T n s = [ n s1   n s2  n s3 ] T n ζ = [ n ζ1   n ζ1   n ζ2  ] T (11)

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

3  NURBS Blending of the Physical Domain

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

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))(12)

In Eq. (12), the symbols α¯i(j) with i=1,2 and j=1,…,4 account for the parametric representation of the outer sides of the reference surface alongside the physical domain, collected in the vectors α¯i=[α¯i(1)(ξ1)α¯i(2)(ξ2)α¯i(3)(ξ1)α¯i(4)(ξ2)]T for i=1,2. In the same way, the four corners of r(α1,α2) denoted with (α1(j),α2(j)) for j=1,…,4 are arranged so that αi=[αi(1)αi(2)αi(3)αi(4)]T with i=1,2. Based on these definitions, it is possible to assess the blending expressions αi(ξ1,ξ2) of the physical domain according to the following matrix-form expressions:

[α1(ξ1,ξ2)α2(ξ1,ξ2)]=[b¯00b¯][B¯00B¯][α~1α~2](13)

being α~i=[α¯iαi]T for i=1,2. Accordingly, the following definitions for b¯ and B¯ should be assessed:

b¯=[bb],B¯=[12I00−14B](14)

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

b=[(1−ξ2)(1+ξ1)(1+ξ2)(1−ξ1)]B=[(1−ξ1)0000(1−ξ2)0000(1+ξ1)0000(1+ξ2)](15)

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

[∂∂α1∂∂α2]=[∂ξ1∂α1∂ξ2∂α1∂ξ1∂α2∂ξ2∂α2][∂∂ξ1∂∂ξ2]=[ξ1,α1ξ2,α1ξ1,α2ξ2,α2][∂∂ξ1∂∂ξ2](16)

where the definitions ξi,αj=∂ξi/∂ξi∂αj∂αj for i,j=1,2 are introduced. In the case of arbitrarily-shaped physical domains, it is very likely that the partial derivatives with respect to ξ1,ξ2 should be referred to the principal coordinate system as well. For this reason, the following transformation should be declared, setting J the Jacobian matrix of Eq. (13):

[∂∂ξ1∂∂ξ2]=[∂α1∂ξ1∂α2∂ξ1∂α1∂ξ2∂α2∂ξ2][∂∂α1∂∂α2]=J[∂∂α1∂∂α2](17)

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

J−1=1det(J)[∂α2∂ξ2−∂α2∂ξ1−∂α1∂ξ2∂α1∂ξ1](18)

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

[∂∂α1∂∂α2]=J−1[∂∂ξ1∂∂ξ2]=[∂ξ1∂α1∂ξ2∂α1∂ξ1∂α2∂ξ2∂α2][∂∂ξ1∂∂ξ2](19)

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

ξ1,α1=∂ξ1∂α1=1det(J)∂α2∂ξ2,ξ1,α2=∂ξ1∂α2=−1det(J)∂α1∂ξ2ξ2,α1=∂ξ2∂α1=−1det(J)∂α2∂ξ1,ξ2,α2=∂ξ2∂α2=1det(J)∂α1∂ξ1(20)

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

det(Jξ1)=∂α1∂ξ1∂2α2∂ξ1∂ξ2−∂α2∂ξ1∂2α1∂ξ1∂ξ2+∂α2∂ξ2∂2α1∂ξ12−∂α1∂ξ2∂2α2∂ξ12det(Jξ2)=−∂α1∂ξ2∂2α2∂ξ1∂ξ2+∂α2∂ξ2∂2α1∂ξ1∂ξ2−∂α2∂ξ1∂2α1∂ξ22+∂α1∂ξ1∂2α2∂ξ22(21)

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

∂2∂α1∂α2=ξ1,α1ξ1,α2∂2∂ξ12+ξ2,α1ξ2,α2∂2∂ξ22+(ξ1,α1ξ2,α2+ξ1,α2ξ2,α1)∂2∂ξ1∂ξ2+ξ1,α1α2∂∂ξ1+ξ2,α1α2∂∂ξ2(22)

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

ξ1,α1α2=1det(J)2(−∂α2∂ξ2∂2α1∂ξ1∂ξ2+∂α2∂ξ2∂α1∂ξ2det(Jξ1)det(J)+∂α2∂ξ1∂2α1∂ξ22−∂α2∂ξ1∂α1∂ξ2det(Jξ2)det(J))ξ2,α1α2=1det(J)2(−∂α2∂ξ1∂2α1∂ξ1∂ξ2+∂α1∂ξ2∂α1∂ξ1det(Jξ1)det(J)+∂α2∂ξ2∂2α1∂ξ12−∂α2∂ξ1∂α1∂ξ1det(Jξ2)det(J))(23)

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

∂2∂αi2=ξ1,αi2∂2∂ξ12+ξ2,αi2∂2∂ξ22+2ξ1,αiξ2,αi∂2∂ξ1∂ξ2+ξ1,αiαi∂∂ξ1+ξ2,αiαi∂∂ξ2fori=1,2(24)

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

ξ1,α1α1=1det(J)2(∂α2∂ξ2∂2α2∂ξ1∂ξ2−(∂α2∂ξ2)2det(Jξ1)det(J)−∂α2∂ξ1∂2α2∂ξ22+∂α2∂ξ1∂α2∂ξ2det(Jξ2)det(J))ξ1,α2α2=1det(J)2(∂α1∂ξ2∂2α1∂ξ1∂ξ2−(∂α1∂ξ2)2det(Jξ1)det(J)−∂α1∂ξ1∂2α1∂ξ22+∂α1∂ξ1∂α1∂ξ2det(Jξ2)det(J))ξ2,α1α1=1det(J)2(−∂α2∂ξ2∂2α2∂ξ12+∂α2∂ξ2∂α2∂ξ1det(Jξ1)det(J)+∂α2∂ξ1∂2α2∂ξ1∂ξ2−(∂α2∂ξ1)2det(Jξ2)det(J))ξ2,α2α2=1det(J)2(−∂α1∂ξ2∂2α1∂ξ12+∂α1∂ξ2∂α1∂ξ1det(Jξ1)det(J)+∂α1∂ξ1∂2α1∂ξ1∂ξ2−(∂α1∂ξ1)2det(Jξ2)det(J))(25)

In the present work, the β-th edge (α¯1(β),α¯2(β)) of the mapped physical domain for β=1,…,4, whose components α¯i(β)(ξj) along αi direction have been collected in the vector α¯i for i=1,2, has been geometrically parametrized with respect to ξj for j=1,2 in terms of NURBS curves [16].

4  ESL Kinematic Equations Employing Generalized Shape Functions

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

U(α1,α2,ζ)=∑τ=0N+1Fτ(ζ)u(τ)(α1,α2)⇔[U1U2U3]=∑τ=0N+1[Fτα1000Fτα2000Fτα3][u1(τ)u2(τ)u3(τ)](26)

Employing the unified notation of Eq. (26), a generalized displacement field component vector u(τ)(α1,α2)=[u1(τ)u2(τ)u3(τ)]T is assessed for each τ-th order of the kinematic expansion, defined in each point of the reference surface r(α1,α2). The dependence of U from the out-of-plane coordinate ζ is taken from the definition, for each τ=0,…,N+1, of a series of thickness functions Fταi along each αi=α1,α2,α3 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 Fταi are properly selected, a wide range of effects can be effectively predicted, thus awarding the model of 3D capabilities. In the present work, Fταi(ζ) have been chosen defined from power polynomials [16]:

Fταi(ζ)={ζτforτ=0,…,N(−1)kzkforτ=N+1(27)

Referring to τ=N+1, a dimensionless coordinate zk has been defined in each k-th layer of the stacking sequence starting from the through-the-thickness coordinate of the top (ζk+1) and the bottom (ζk) surface delimiting the k-th sheet, respectively. If k=1,…,l with l denoting the total number of laminae, one gets:

zk=(2ζk+1−ζkζ−ζk+1+ζkζk+1−ζk)(28)

As a matter of fact, the employment of the unified notation of Eq. (26) lets us to perform a systematic analysis with different ESL theories, since the theoretical formulation is independent from the actual selection of the thickness function analytical expression. Nevertheless, since various polynomials orders can be selected in Eq. (27), a nomenclature is adopted for a smarter identification of the axiomatic through-the-thickness assumptions of the simulation. In particular, the acronym ED(Z)−N 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 u(N+1)(α1,α2) is introduced according to what exerted in Eq. (27), the letter “Z” is included. In this way, the formulation accounts for the zigzag function in each αi=α1,α2,α3 component [16].

Once the ESL assessment of the displacement field component vector U(α1,α2,ζ)=[U1U2U3]T 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 IN×IM points, whose generic element can be identified with (α1f,α2g) for f=1,…,IN and g=1,…,IM. For this purpose, an interpolation algorithm is introduced for u(τ)(α1,α2) at each τ=0,…,N+1 kinematic expansion order. A very useful tool is the vectorization operator [16], accounting for a rearrangement with a 1D vector of a 2D matrix. If we denote with A a matrix of dimensions IN×IM of components Aij with i=1,…,IN and j=1,…,IM, a column vector A→=Vec(A) of dimensions INIM×1 is defined, whose arbitrary component Ak for k=1,…,INIM reads as:

A→=Vec(A)   ⇔   Ak=(Aij)k for k=i+(j−1)INi=1,…,IN    j=1,…,IM(29)

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

u→i(τ)=[ui(11)(τ)⋯ui(IN1)(τ)ui(12)(τ)⋯ui(IN2)(τ)ui(1IM)(τ)⋯ui(INIM)(τ)]Tfori=1,2,3(30)

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:

u¯(τ)=[u→1(τ)Tu→2(τ)Tu→3(τ)T]T(31)

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

u(τ)=NTu¯(τ)(32)

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

[u1(τ)(α1,α2)u2(τ)(α1,α2)u3(τ)(α1,α2)]=∑f=1IN∑g=1IM[lf(α1)lg(α2)000lf(α1)lg(α2)000lf(α1)lg(α2)][u1(τ)(α1f,α2g)u2(τ)(α1f,α2g)u3(τ)(α1f,α2g)](33)

where ui(τ)(α1,α2) for i=1,2,3 in a generic point is obtained from the interpolation of all the values ui(τ)(α1f,α2g) assumed by the generalized displacement field component in each (α1f,α2g) point of the above introduced discrete grid. To this end, the well-known univariate Lagrange interpolating polynomials lη(αi) of (IP−1)-th order defined along the αi=α1,α2 principal direction of the physical domain are employed, with η=f,g and IP=IN,IM. According to reference [16], lη(αi) can be computed as follows:

lη(α i)=∏j=1I P(α i−α ij)(α i−α iη)∏j=1,  j≠ηI P(α iη−α ij)             I P=I N,I Mfor   i=1,2                      η=f,g         (34)

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

NT(α1,α2)=[lα2⊗lα1000lα2⊗lα1000lα2⊗lα1](35)

In Eq. (35), the Lagrange interpolating polynomials introduced in Eq. (34) are collected in a 1×IP vector with IP=IN,IM denoted with lαi for each αi=α1,α2, according to the following expanded form:

l α i=[ l 1(α i)   ⋯   l η(α i)   ⋯   l IP(α i) ]             I P=I N,I Mfor   i=1,2                      η=f,g         (36)

In the same way, vector lαi(1) accounting for the first order derivative lη(1)(αi)=∂lη/∂lη∂αi∂αi of the Lagrange interpolating polynomials collected according to Eq. (36) with respect to αi=α1,α2 is defined as follows, setting η=f,g:

N¯T(α 1,α 2)=lα 2⊗lα 1 ⇔ NT=[ N¯T000N¯T000N¯T ](37)

Since the Kronecker tensorial product, denoted with ⊗, has been adopted in Eq. (35), it is useful to introduce a N¯T row vector of dimensions 1×INIM so that NT matrix can be expressed in a compact form as follows [16]:

N¯T(α1,α2)=lα2⊗lα1⇔NT=[N¯T000N¯T000N¯T](38)

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

U(α1,α2,ζ)=∑τ=0N+1Fτu(τ)=∑τ=0N+1FτNTu¯(τ)(39)

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

ε=DU=Dζ(∑i=13DΩαi)U(40)

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

Dζ=[1H10000000001H20000000001H11H20000000001H10∂∂ζ00000001H20∂∂ζ000000000∂∂ζ](41)

In the same way, DΩαi for αi=α1,α2,α3 read as:

DΩα1=[ D¯Ωα1  0    0 ]DΩα2=[ 0    D¯Ωα2  0 ]DΩα3=[ 0    0    D¯Ωα3 ](42)

being

D¯Ωα1=[1A1∂∂α11A1A2∂A2∂α1−1A1A2∂A1∂α21A2∂∂α2−1R10100]TD¯Ωα2=[1A1A2∂A1∂α21A2∂∂α21A1∂∂α1−1A1A2∂A2∂α10−1R2010]TD¯Ωα3=[1R11R2001A1∂∂α11A2∂∂α2001]T(43)

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

ε=∑τ=0N+1∑i=13DζDΩαiFτu(τ)=∑τ=0N+1∑i=13Z(τ)αiDΩαiu(τ)=∑τ=0N+1∑i=13Z(τ)αiε(τ)αi(44)

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

Z(τ)αi=[FταiH1000000000FταiH2000000000FταiH1FταiH2000000000FταiH10∂Fταi∂ζ0000000FταiH20∂Fταi∂ζ000000000∂Fταi∂ζ](45)

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

ε(τ)αi=DΩαiu(τ)=DΩαiNTu¯(τ)=Bαiu¯(τ)(46)

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

Bα19× (I​N×I​M)=[ 1A 1lα2⊗lα1(1)1A 1A 2∂A 2∂α1lα2⊗lα1−1A 1A 2∂A 1∂α2lα2⊗lα11A 2lα2(1)⊗lα1−1R 1lα2⊗lα10lα2⊗lα100 ],      Bα29× (I​N×I​M)=[ 1A 1A 2∂A 1∂α2lα2⊗lα11A 2lα2(1)⊗lα11A 1lα2⊗lα1(1)−1A 1A 2∂A 2∂α1lα2⊗lα10−1R 2lα2⊗lα10lα2⊗lα10 ],      Bα39× (I​N×I​M)=[ 1R 1lα2⊗lα11R 2lα2⊗lα1001A 1lα2⊗lα1(1)1A 2lα2(1)⊗lα100lα2⊗lα1 ](47)

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

5  Generally Anisotropic Lamination Scheme Assessment

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 k-th layer with k=1,…,l is characterized by a thickness hk according to the conventions of Eq. (8) and a reference system O′α^1(k)α^2(k)ζ^(k) orientated along the material symmetry axes. Accordingly, the formulation has been developed by assuming ζ^(k)=ζ. In addition, a characteristic angle θk describing the orientation of α^1(k) axis with respect to α1 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 σ^(k)(α^1(k),α^2(k),ζ)=[σ^1(k)σ^2(k)τ^12(k)τ^13(k)τ^23(k)σ^3(k)]T and ε^(k)(α^1(k),α^2(k),ζ)=[ε^1(k)ε^2(k)γ^12(k)γ^13(k)γ^23(k)ε^3(k)]T from the introduction of a generally anisotropic stiffness matrix E(k) for each k-th layer. The generalized constitutive relationship can be thus defined as [16]:

σ^(k)=E(k)ε^(k)fork=1,…,l(48)

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

E(k)=[E11(k)E12(k)E16(k)E14(k)E15(k)E13(k)E12(k)E22(k)E26(k)E24(k)E25(k)E23(k)E16(k)E26(k)E66(k)E46(k)E56(k)E36(k)E14(k)E24(k)E46(k)E44(k)E45(k)E34(k)E15(k)E25(k)E56(k)E45(k)E55(k)E35(k)E13(k)E23(k)E36(k)E34(k)E35(k)E33(k)]fork=1,…,l(49)

where Eij(k) with i,j=1,…,6 relates the i-th component of σ^(k) vector to the j-th strain component of the ε^(k) vector. Accordingly, the coefficients employed in Eq. (49) account for the 3D elastic stiffnesses of the material, namely Eij(k)=Cij(k). When E(k) is referred to the plane stress assumption (σ3(k)=0), the reduced stiffness coefficients Eij(k)=Qij(k) can be employed, defined from the expression reported hereafter:

Qij(k)=Cij(k)−Cj3(k)Ci3(k)C33(k)fori,j=1,…,6,k=1,…,l(50)

Nevertheless, Eq. (48) should be referred to the curvilinear principal reference system of the physical domain, accounting for the 3D stress and strain vectors σ(k)(α1,α2,ζ)=[σ1(k)σ2(k)τ12(k)τ13(k)τ23(k)σ3(k)]T and ε(k)(α1,α2,ζ)=[ε1(k)ε2(k)γ12(k)γ13(k)γ23(k)ε3(k)]T. To this end, a generalized rotation matrix T(k) is defined according to reference [30] within each layer starting from the previously discussed angle θk employed to compute the rotated stiffness matrix E¯(k)=T(k)E(k)(T(k))T with components E¯ij(k) for i,j=1,…,6 referred to the shell geometric reference system O′α1α2ζ. More in detail, Eq. (48) takes the following form:

σ(k)=E¯(k)ε(k)⇔[σ1(k)σ2(k)τ12(k)τ13(k)τ23(k)σ3(k)]=[E¯11(k)E¯12(k)E¯16(k)E¯14(k)E¯15(k)E¯13(k)E¯12(k)E¯22(k)E¯26(k)E¯24(k)E¯25(k)E¯23(k)E¯16(k)E¯25(k)E¯66(k)E¯46(k)E¯56(k)E¯36(k)E¯14(k)E¯24(k)E¯46(k)E¯44(k)E¯45(k)E¯34(k)E¯15(k)E¯25(k)E¯56(k)E¯45(k)E¯55(k)E¯35(k)E¯13(k)E¯23(k)E¯36(k)E¯34(k)E¯35(k)E¯33(k)][ε1(k)ε2(k)γ12(k)γ13(k)γ23(k)ε3(k)]fork=1,…,l(51)

On the other hand Eq. (51), defined for the k-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 S(τ)αi(α1,α2)=[N1(τ)αiN2(τ)αiN12(τ)αiN21(τ)αiT1(τ)αiT2(τ)αiP1(τ)αiP2(τ)αiS3(τ)αi]T, associated to each τ-th order of the kinematic expansion of Eq. (39). In this way, Eq. (51) can be expressed in terms of S(τ)αi and ε(η)αi vectors, according to the following equations directly derived from the static equivalence principle:

S(τ)αi=∑η=0N+1∑j=13A(τη)αiαjε(η)αi    for    τ=0,...,N+1,   α i=α 1,α 2,α 3 (52)

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

whose generic component A nm  (pq)(τη)[ fg ]α iα j is defined as:

Anm  (pq)(τη)[ fg ]α iα j=∑k=1l∫ζ kζ k+1B¯ nm(k)∂fF ηα j∂ζf∂gF τα i∂ζgH 1H 2H 1pH 2qdζ           τ,η=0,...,N+1           n,m=1,...,6for     p,q=0,1,2           α i,α j=α 1,α 2,α 3           f,g=0,1(54)

being Fηαj=∂0Fηαj/∂0Fηαj∂ζ0∂ζ0 and Fταi=∂0Fταi/∂0Fταi∂ζ0∂ζ0. As can be seen, four blocks can be traced in Eq. (53) accounting for the presence of the derivatives of the f,g-th order within the expression reported in Eq. (54). Furthermore, in Eq. (54) the term B¯nm(k) refers to the 3D rotated stiffness coefficient E¯ij(k) introduced in Eq. (51), possibly adjusted by the shear correction factor κ(ζ) [16] in the case of displacement thickness functions employed in Eq. (26) neglecting the out-of-plane stretching effects:

B¯nm(k)={E¯nm(k)forn,m=1,2,3,6κ(ζ)E¯nm(k)forn,m=4,5(55)

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

S(τ)αi=∑η=0N+1∑j=13A(τη)αiαjBαju¯(η)    for    τ=0,...,N+1,   α i=α 1,α 2,α 3 (56)

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

S(τ)αi=∑η=0N+1∑j=13A(τη)αiαjDΩαjNTu¯(η)=∑η=0N+1∑j=13O(τη)αiαjNTu¯(η)for    τ=0,…,N+1,αi=α1,α2,α3(57)

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

6  Equilibrium Relations in the Variational Form

In the present section, the equilibrium equations for an arbitrary doubly-curved shell are derived by employing an energy approach. The stationary configuration is specifically enforced to the total potential energy virtual variation δΠ. The contribution of external loads should be inserted within the ESL framework, accounting for various pre-determined distributions along the top (ζ=+h/h22) and the bottom (ζ=−h/h22) surfaces of the shell. Moreover, the formulation embeds the contribution of linear elastic springs distributed along the shell edges. Suppose we denote with Les,Leb 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 [16]:

δΠ=δΦ−δLes−δLeb=0(58)

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

6.1 Elastic Strain Energy

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

δΦ=∑k=1l∫α1∫α2∫ζkζk+1(δε(k))Tσ(k)A1A2H1H2dα1dα2dζ(59)

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

δΦ=∑τ=0N+1∑i=13∫α1∫α2(δε(τ)αi)TS(τ)αiA1A2dα1dα2(60)

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:

δ Φ=∑τ=0N+1∑η=0N+1∑i=13∑j=13∫α 1∫α 2(Bα iδu¯(τ))TA(τη)α iα j(Bα ju¯(η))A 1A 2dα 1dα 2          =∑τ=0N+1∑η=0N+1(δu¯(τ))T(∑i=13∑j=13∫α 1∫α 2(Bα i)TA(τη)α iα j(Bα j)A 1A 2dα 1dα 2)u¯(η)(61)

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

KΦ(τη)=∫α1∫α2[(Bα1)TA(τη)α1α1Bα1(Bα1)TA(τη)α1α2Bα2(Bα1)TA(τη)α1α3Bα3(Bα2)TA(τη)α2α1Bα1(Bα2)TA(τη)α2α2Bα2(Bα2)TA(τη)α2α3Bα3(Bα3)TA(τη)α3α1Bα1(Bα3)TA(τη)α3α2Bα2(Bα3)TA(τη)α3α3Bα3]A1A2dα1dα2=∫α1∫α2[KΦ(τη)α1α1KΦ(τη)α1α2KΦ(τη)α1α3KΦ(τη)α2α1KΦ(τη)α2α2KΦ(τη)α2α3KΦ(τη)α3α1KΦ(τη)α3α2KΦ(τη)α3α3]A1A2dα1dα2(62)

In Appendix B, an extended expression of each stiffness matrix coefficient KΦ(τη)αiαj is reported, setting τ,η=0,…,N+1 and αi,αj=α1,α2,α3.

6.2 External Loads General Distributions

Let us consider two generic vectors qa(1) and qa(2) of external loads applied to the doubly-curved shell structure at issue, whose components are referred in the above discussed geometric reference system O′α1α2ζ as shown in Fig. 2:

q a(1)=[ q 1a(1)   q 2a(1)   q 3a(1) ]T=[ q 1a(+)   q 2a(+)   q 3a(+) ]Tq a(2)=[ q 1a(2)   q 2a(2)   q 3a(2) ]T=[ q 1a(−)   q 2a(−)   q 3a(−) ]T(63)

where symbols (+) and (−) account for quantities referred to ζ=+h/h22 and ζ=−h/h22, respectively. Accordingly, a general variation q¯(α1,α2) dependent on in-plane principal coordinates can be assigned to each component qia(j) for i=1,2,3 and j=1,2 of external load vectors introduced in Eq. (63). If we denote with qi(j) a reference scaling value of loads applied along the αi principal direction, the following relation can be assessed (Fig. 2):

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

qia(j)(α1,α2,±h2)=qi(j)|ζ=±h/2q¯(α1,α2)       for  i=1,2,3j=1,2    (64)

For a constant loading distribution, it is:

q¯(α1,α2)=1(65)

In the present manuscript, two different distributions for qa(j) components with j=1,2 are declared. A bivariate Gaussian distribution is now reported:

q¯(α1,α2)=e−12(1−ρ122)((α1−α1mΔ1)2+(α2−α2mΔ2)2−2ρ12α1−α1mΔ1α2−α2mΔ2)(66)

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

q¯(α1,α2)=e−(|(α1−α1m)cos⁡β+(α2−α2m)sin⁡βΔ1|n+ |−(α1−α1m)sin⁡β+(α2−α2m)cos⁡βΔ2|n)(67)

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

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

δLes(3D)=∫α1∫α2(∑j=12(∑i=13qi(j)δUi(j))H1(j)H2(j))A1A2dα1dα2(68)

being Hi(j) for j=1,2 the through-the-thickness curvature parameters calculated by means of Eq. (7) for ζ=+h/h22 and ζ=−h/h22. On the other hand, if we consider a virtual variation δu(τ) of the generalized displacement field component vector introduced in Eq. (26), a new generalized vector of external actions q(τ)(α1,α2)=[q1(τ)q2(τ)q3(τ)]T lying on the reference surface can be introduced for each τ=0,…,N+1. Thus, the virtual work δLes reads as follows:

δLes=∑τ=0N+1∫α1∫α2∑i=13qi(τ)δui(τ)A1A2dα1dα2=∑τ=0N+1∫α1∫α2(δu(τ))Tq(τ)A1A2dα1dα2(69)

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

δLes(3D)=δLes(70)

Starting from Eq. (70), the extended expressions of the q(τ) components are obtained for the τ-th kinematic expansion order, being Fταi(j) for j=1,2 the thickness function employed for the assessment of the displacement field components according to Eq. (26) calculated at ζ=+h/h22 and ζ=−h/h22, respectively:

qi(τ)=∑j=12qi(j)Fταi(j)H1(j)H2(j)fori=1,2,3(71)

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

δLes=∑τ=0N+1∫α1∫α2(NTδu¯(τ))Tq(τ)A1A2dα1dα2=∑τ=0N+1(δu¯(τ))TQ(τ)(72)

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

Q(τ)=∫α1∫α2[N¯q1(τ)N¯q2(τ)N¯q3(τ)]A1A2dα1dα2(73)

6.3 Boundary Springs

The last contribution to the total potential energy occurring in Eq. (58) accounts for the influence of the elastic constraints distributed along the edges of the structure. In particular, a set of linear elastic springs of stiffness kif(k)αjm orientated along the shell side located at αj=αjm with αj=α1,α2,α3 and m=0,1 distributed along the αi=α1,α2,α3 principal direction has been considered, according to the theoretical development reported in reference [16]. The superscript k stands for the stacking sequence layer. Accordingly, they induce some boundary stresses after applying the corresponding 3D displacement component Ui=U1,U2,U3. If the physical domain consists of a rectangular 2D interval [α10,α11]×[α20,α21], different stresses are induced on each boundary. Referring to α1=α1m for m=0,1, the following definitions can be stated at the k-th lamina level, setting α2∈[α20,α21]:

σ¯1(k)(α1m,α2,ζ)=−k1f(k)α1mf(α1m,α2)U1(α1m,α2,ζ)τ¯12(k)(α1m,α2,ζ)=−k2f(k)α1mf(α1m,α2)U2(α1m,α2,ζ)for m=0,1k=1,…,lτ¯13(k)(α1m,α2,ζ)=−k3f(k)α1mf(α1m,α2)U3(α1m,α2,ζ)(74)

where f(α1m,α2) accounts for an arbitrary univariate function in the α2 coordinate located at α1=α1m. On the other hand, stresses τ¯12(k),σ¯2(k),τ¯23(k) can be found at α2=α2m according to the following definitions:

τ¯12(k)(α1,α2m,ζ)=−k1f(k)α2mf(α1,α2m)U1(α1,α2m,ζ)σ¯2(k)(α1,α2m,ζ)=−k2f(k)α2mf(α1,α2m)U2(α1,α2m,ζ)form=0,1k=1,…,lτ¯23(k)(α1,α2m,ζ)=−k3f(k)α2mf(α1,α2m)U3(α1,α2m,ζ)(75)

being α1∈[α10,α11]. In the same way, Eq. (75) contains a univariate dispersion of linear elastic springs denoted with f(α1,α2m) for α2=α2m. In order to provide an effective expression of both f(α1m,α2) and f(α1,α2m) regardless of the actual dimension of the physical domain, a new dimensionless coordinate ξ¯∈[0,1] is introduced, according to the following definition:

ξ¯=αr−αr0αr1−αr0for r=1,2(76)

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

f(ξ)=e−(ξ¯−ξ¯mξ~m)p(77)

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

f(ξ)=1−e−(ξ¯ξ¯m)p+e−(1−ξ¯ξ~m)p(78)

where ξ¯m,ξ~m∈[0,1] 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 α1=α1m with m=0,1:

[N¯1(τ)α1N¯12(τ)α2T¯1(τ)α3]=−∑η=0N+1[L1(2)α1mfb(τη)α1000L2(2)α1mfb(τη)α2000L3(2)α1mfb(τη)α3][u1(η)u2(η)u3(η)]=−∑η=0N+1L(1)α1mfb(τη)u(η)for m=0,1forα2∈[α20,α21](79)

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

[ N¯ 21(τ)α1N¯ 2(τ)α2T¯ 2(τ)α3 ]=−∑η=0N+1[ L 1(1)α 2mfb(τη)α 1000L 2(1)α 2mfb(τη)α 2000L 3(1)α 2mfb(τη)α 3 ][ u 1(η)u 2(η)u 3(η) ]=−∑η=0N+1L (1)α 2mfb(τη)u (η)       for   m=0,1               for   α 1∈[ α 20,α 21 ](80)

Fundamental coefficients Li(p)αnmfb(τη)αi occurring in Eqs. (79) and (80) can be computed for each τ,η=0,…,N+1 according to the following relation, setting n,p=1,2, i=1,2,3 and n,p=1,2:

Li(p)αnmfb(τη)αi=∑k=1l∫ζkζk+1kif(k)αnmFηαiFταiHpdζ(81)

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

δLeb=δLeb1+δLeb2=∑τ=0N+1(∮α2(δu(τ))TS¯α1(τ)A2α2+∮α1(δu(τ))TS¯α2(τ)A1α1)(82)

where S¯α1(τ) components are defined according to what exerted in Eq. (80). S¯α2(τ) 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 δLeb comes out:

δLeb=δLeb1+δLeb2=∑τ=0N+1((δu¯(τ))T(Q¯α1(τ)+Q¯α2(τ)))(83)

being Q¯α1(τ),Q¯α2(τ) the ESL boundary load vectors of Eqs. (79) and (80) expressed in terms of the generalized interpolation of the displacement field, as follows [16]:

Q¯α1(τ)=∮α2[N¯N¯1(τ)α1N¯N¯12(τ)α2N¯T¯1(τ)α3]A2α2,Q¯α2(τ)=∮α1[N¯N¯21(τ)α1N¯N¯2(τ)α2N¯T¯2(τ)α3]A1α1(84)

Thus, the generalized displacement vector u¯(τ) 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:

δLeb=∑τ=0N+1(δu¯(τ))T∑η=0N+1(∑m=01(∮α2NL(1)α2mfb(τη)NTA2α2+∮α1NL(1)α1mfb(τη)NTA1α1))u¯(η)=∑τ=0N+1(δu¯(τ))T∑η=0N+1K¯b(τη)u¯(η)(85)

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

6.4 Fundamental Governing Equations

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

∑η=0N+1(K¯Φ(τη)−K¯b(τη))u¯(η)=Q(τ)⇔∑η=0N+1K¯(τη)u¯(η)=Q(τ)forτ=0,…,N+1(86)

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

[K¯(00)K¯(01)⋯K¯(0(N))K¯(0(N+1))K¯(10)K¯(11)⋯K¯(1(N))K¯(1(N+1))⋮⋮⋱⋮⋮K¯((N)0)K¯((N)1)⋯K¯((N)(N))K¯((N)(N+1))K¯((N+1)0)K¯((N+1)1)⋯K¯((N+1)(N))K¯((N+1)(N+1))][U¯(0)U¯(1)⋮U¯(N)U¯(N+1)]=[Q(0)Q(1)⋮Q(N)Q(N+1)]⇔K¯u¯=Q(87)

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

KΦ(τη)=∫−11∫−11[KΦ(τη)α1α1KΦ(τη)α1α2KΦ(τη)α1α3KΦ(τη)α2α1KΦ(τη)α2α2KΦ(τη)α2α3KΦ(τη)α3α1KΦ(τη)α3α2KΦ(τη)α3α3]A1A2det(J)dξ1dξ2(88)

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

Q(τ)=∫−11∫−11[N¯q1(τ)N¯q2(τ)N¯q3(τ)]A1A2det(J)dξ1dξ2(89)

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

Q¯ n(i)(τ)=∫0L(i)[ N¯N¯ n(τ)α 1N¯N¯ ns(τ)α 2N¯T¯ ζ(τ)α 3 ] dsn(i)      for    i=1,...,4          τ=0,...,N+1(90)

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

L(i)=∫−11A12(dα1dξj)2+A22(dα2dξj)2dξjfori=1,…,4(91)

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

Q¯n(i)(τ)=∫−11[N¯N¯n(τ)α1N¯N¯ns(τ)α2N¯T¯ζ(τ)α3]L(i)2dη(i)fori=1,..,4(92)

Once the quantity Q¯n(i)(τ) 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 K¯ 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 Q is employed for the external surface loads.

In order to compute the generalized stress resultants N¯n(τ)α1,N¯ns(τ)α2,T¯ζ(τ)α3 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 un(τ),us(τ),uζ(τ) referred to nn,ns,nζ local reference system in terms of u1(τ),u2(τ),u3(τ) for each τ-th order of the kinematic expansion [16]:

[un(τ)us(τ)uζ(τ)]=[nn1nn20ns1ns20001][u1(τ)u2(τ)u3(τ)](93)

In the same way, the generalized stress resultants N¯n(τ)α1,N¯ns(τ)α2,T¯ζ(τ)α3 should be expressed in terms of the S(τ)αi vector for αi=α1,α2,α3 according to the following relations [16]:

[Nn(τ)α1Nns(τ)α2Tζ(τ)α3]=[nn12nn22nn1nn2nn1nn200nn1ns1nn2ns2nn1ns2nn2ns1000000nn1nn2][N1(τ)α1N2(τ)α1N12(τ)α1N21(τ)α1T1(τ)α3T2(τ)α3](94)

7  Numerical Implementation with the GDQ Method

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

∂n+mf(ξ 1,ξ 2)∂ξ 1n∂ξ 2m|ξ 1=ξ1i, ξ 2=ξ2j≅∑k=1I Nςikξ 1(n)(∑l=1IMςjlξ 2(m)f(ξ 1k,ξ 2l))   i = 1, 2,..., INj = 1, 2,..., IM(95)

where IN,IM denote the number of the selected discrete points along ξ1,ξ2 directions, respectively. As can be seen, Eq. (95) accounts for a quadrature rule for derivative purposes. The weighting coefficients ςpvξr(q) for r=1,2, p=i,j, q=n,m and v=k,l are computed utilizing the recursive expression valid for each q≥1 [16]:

ςpvξ r(1)=L(1)(ξ rp)(ξ rp−ξ rv)L(1)(ξ rv),          ςpvξ r(q)=q(ςpvξ r(1)ςppξ r(q−1)−ςpvξ r(q−1)ξ rp−ξ rv)p≠vςppξ r(q)=−∑v=1 v≠pNςpvξ r(q)p=v(96)

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

ςpvξr(0)=δpvforr=1,2(97)

being δpv the Kronecker delta operator. For a generic non-uniform grid of IN×IM points, the GDQ weighting coefficients for the derivative of the (n+m)-th order are calculated by means of Eqs. (96) and (97) are collected in a INIM×INIM squared matrix denoted with ςξ1(n)ξ2(m). Starting from Eq. (96), the matrix at issue can be obtained from a proper expansion of matrix ςξr(q) for q=n,m and r=1,2 of dimensions IQ×IQ with IQ=IN,IM, respectively, accounting for the derivatives of n-th and m-th order along ξ1,ξ2:

ςξ1(n)ξ2(m)=ςξ2(m)⊗ςξ1(n)(98)

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

f→(n+m)=ςξ1(n)ξ2(m)f→(99)

where f is a IN×IM matrix whose general component is f(ξ1i,ξ1j) for i=1,…,IN and j=1,…,IM. In the same way, the generic element of f(n+m) is ∂n+mf(ξ1i,ξ2j)/∂n+mf(ξ1i,ξ2j)∂ξ1n∂ξ2m∂ξ1n∂ξ2m.

A 2D Legendre-Gauss-Lobatto (LGL) grid of dimensions IN×IM referred to the domain [−1,1]×[−1,1] has been here adopted, such that:

(1−ξ r2)A I​P−1(ξ r)=0      for    I P=I N,I Mr=1,2           (100)

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

A I​P−1(ξ r)=ddξ r(L I​P(ξ r))      for    I P=I N,I M        r=1,2                    ξ r∈[ −1,1 ]             (101)

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

ςα 1(1)α 2(0)=P→ 11∘ςξ 1(1)ξ 2(0)+P→ 21∘ςξ 1(0)ξ 2(1)ςα 1(0)α 2(1)=P→ 12∘ςξ 1(1)ξ 2(0)+P→ 22∘ςξ 1(0)ξ 2(1)(102)

being P→ij with i,j=1,2 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 Prij with r,i,j=1,2:

ςα1(2)α2(0)=P→11o2∘ςξ1(2)ξ2(0)+P→21o2∘ςξ1(0)ξ2(2)+2P→11∘P→21∘ςξ1(1)ξ2(1)+P→111∘ςξ1(1)ξ2(0)+P→211∘ςξ1(0)ξ2(1)ςα1(0)α2(2)=P→12o2∘ςξ1(2)ξ2(0)+P→22o2∘ςξ1(0)ξ2(2)+2P→12∘P→22∘ςξ1(1)ξ2(1)+P→122∘ςξ1(1)ξ2(0)+P→222∘ςξ1(0)ξ2(1)ςα1(1)α2(1)=P→11∘P→12∘ςξ1(2)ξ2(0)+P→21∘P→22∘ςξ1(0)ξ2(2)+(P→11∘P→22+P→12∘P→21)∘ςξ1(1)ξ2(1)+P→112∘ςξ1(1)ξ2(0)+P→212∘ςξ1(0)ξ2(1)(103)

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 IN×IM 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 δb, whereas those referring to the former are arranged within the vector δd. Thus, the fundamental system accounts as follows [16]:

[K¯bbK¯bdK¯dbK¯dd][δbδd]=[Q¯bQ¯d](104)

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

(K¯dd−K¯db(K¯bb)−1K¯bd)δd=Q¯d−K¯db(K¯bb)−1Q¯b(105)

The numerical integrations that are involved in the formulation of the present manuscript are performed by means of the GIQ procedure [16]. Referring to a generic univariate function f=f(x) defined in the closed interval [a,b] with a<b∈R, a discrete grid is selected from the 1D domain, namely xk∈[a,b] with k=1,2,…,IN. Accordingly, the integral of f restricted to the sub-interval [xi,xj]⊆[a,b] with xi,xj belonging to the grid at issue can be computed as a weighted sum of the values f(xk) for k=1,2,…,IN assumed by the function at each point of the computational grid of [a,b]:

∫xixjf(x)dx=∑k=1INwkijf(xk)(106)

where wkij=wjk−wik are referred to the interval [xi,xj]. The coefficients wjk,wik are computed following the extensive procedure of reference [16], starting from the shifted GDQ weighting coefficients reported hereafter, setting ε=1×10−10:

ς¯ ijζ(1)=ζ i−εζ j−ες ijζ(1)            for   i≠jς¯ ijζ(1)=ς iiζ(1)+1ζ i−ε       for   i=j(107)

Then, they are collected in the matrix ς¯ζ(1) of dimensions IN×IN. Thus, the coefficients wik,wjk referred to the integral of f restricted to the interval [ε,xr]⊆[a,b] are the elements of the GIQ matrix, reading as:

W=(ς¯ζ(1))−1(108)

All the terms wkij introduced in Eq. (106) can be collected in a vector of dimensions 1×IN denoted with w(x), setting x the integration variable. In the same way, the 1×IM vector w(y) can be introduced when the integration of the univariate function f=f(y) is required. For a bivariate function f=f(x,y), Eq. (106) turns into the following expression, setting a generic 2D IN×IM grid from a rectangular domain [a,b]×[c,d] of extremes a<b∈R and c<d∈R [16]:

∫cd∫abf(x,y)dxdy=∑i=1IN∑j=1IMwi1INwj1IMf(xi,yj)=∑k=1INMwk1INMfk(109)

where fk values for k=1,…,INM=INIM are arranged employing the above introduced by-column vectorization operation of Eq. (29). On the other hand, wk1INIM GIQ coefficients are computed using the previously introduced vectors w(x),w(y), as reported:

w(xy)=w(y)⊗w(x)=[w11IM⋯wIM1IM]⊗[w11IN⋯wIN1IN](110)

Referring to the bivariate function f=f(x,y), 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, w(x0) and w(0y) are defined starting from Eq. (110):

w(0y)=w(y)⊗o(x)w(x0)=o(y)⊗w(x)(111)

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

∫abf(x,y)dx=∑i=1IN∑j=1IMwi1INwj(0)f(xi,yj)∫cdf(x,y)dy=∑i=1IN∑j=1IMwi(0)wj1IMf(xi,yj)(112)

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

P=(w(ξ1ξ2)∘(A→1∘A→2∘D→J)T)⊗O(113)

being O a ones squared matrix of dimensions 3INIM. Accordingly, w(ξ1ξ2)=w(ξ2)⊗w(ξ1) following Eq. (110), whereas DJ, A1, and A2 accounts for the determinant of the Jacobian matrix J, as well as the Lamè parameters A1,A2 of the shell. If the index k=n+(m−1)IN with n=1,…,IN and m=1,…,IM is assessed, we obtain the vectorized form of KΦ(τη), accounting for the IN×IM grid of the computational domain:

K~Φ(τη)=[KΦ(1)(τη)αiαjKΦ(2)(τη)αiαj⋯KΦ(k)(τη)αiαj⋯KΦ(INIM)(τη)αiαj]T(114)

Employing Eqs. (113) and (114), the stiffness matrix K^Φ(nm)(τη) for an arbitrarily-shaped shell structure can be computed:

K^Φ(nm)(τη)=PK~Φ(τη)(115)

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

N~1(τ)α1=w(0α2)T∘A→2∘N¯1(τ)α1N~12(τ)α2=w(0α2)T∘A→2∘N¯12(τ)α2T~1(τ)α3=w(0α2)T∘A→2∘T¯1(τ)α3N~21(τ)α1=w(α10)T∘A→1∘N¯21(τ)α1N~2(τ)α2=w(α10)T∘A→1∘N¯2(τ)α2T~2(τ)α3=w(α10)T∘A→1∘T¯2(τ)α3(116)

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

N˜n(τ)α1=L(1)2w˜(α10)T∘N¯n(τ)α1N˜ns(τ)α2=L(1)2w˜(α10)T∘N¯ns(τ)α2T˜ζ(τ)α3=L(1)2w˜(α10)T∘T¯ζ(τ)α3(117)

The edge characterized by ξ1=1 behaves as:

N˜n(τ)α1=L(2)2w˜(0α2)T∘N¯n(τ)α1N˜ns(τ)α2=L(2)2w˜(0α2)T∘N¯ns(τ)α2T˜ζ(τ)α3=L(2)2w˜(0α2)T∘T¯ζ(τ)α3(118)

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

N˜n(τ)α1=L(3)2w˜(α10)T∘N¯n(τ)α1N˜ns(τ)α2=L(3)2w˜(α10)T∘N¯ns(τ)α2T˜ζ(τ)α3=L(3)2w˜(α10)T∘T¯ζ(τ)α3(119)

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

N˜n(τ)α1=L(4)2w˜(0α2)T∘N¯n(τ)α1N˜ns(τ)α2=L(4)2w˜(0α2)T∘N¯ns(τ)α2T˜ζ(τ)α3=L(4)2w˜(0α2)T∘T¯ζ(τ)α3(120)

8  Post-Processing Analysis

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

ξ g=−cos (g−1I T−1π)      for   ξ g∈[ −1,1 ]g=1,...,I T(121)

being IT the number of the selected discrete points. Referring to the interval [ζk,ζk+1] for k=1,…,l representing the k-th layer of the staking sequence within the 3D physical domain, Eq. (121) should be transformed according to the following relation:

ζg=ζ k+1−ζ k2(ξg+1)+ζ k      for   ξ g∈[ −1,1 ]g=1,...,I Tk=1,...,l(122)

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 U(ijg)(k) with i=1,…,IN, j=1,…,IM and g=1,…,IT for each k-th layer of the stacking sequence [16]:

U(ijg)(k)=∑τ=0N+1Fτ(g)(k)u(ij)(τ)(α1i,α2j)⇔[U1(ijg)(k)(α1i,α2j,ζg(k))U2(ijg)(k)(α1i,α2j,ζg(k))U3(ijg)(k)(α1i,α2j,ζg(k))]=∑τ=0N+1[Fτ(g)α1(k)(ζg(k))000Fτ(g)α2(k)(ζg(k))000Fτ(g)α3(k)(ζg(k))][u1(ij)(τ)(α1i,α2j)u2(ij)(τ)(α1i,α2j)u3(ij)(τ)(α1i,α2j)](123)

being Fτ(g)α1(k)(ζg(k)),Fτ(g)α2(k)(ζg(k)),Fτ(g)α3(k)(ζg(k)) the generalized thickness functions set calculated at the (α1i,α2j,ζg(k)) point of the 3D solid, whereas u1(ij)(τ),u2(ij)(τ),u3(ij)(τ) are the generalized displacement field components for the τ-th kinematic expansion order referred to the point (α1i,α2j) of the physical domain, for i=1,…,IN, j=1,…,IM. In the same way, the generalized strain component vector ε(ij)(τ)αq for q=1,2,3 and τ=0,…,N+1 are computed in terms of the discrete 3D strain component vector ε(ijg)(k)=[ε1(ijg)(k)ε2(ijg)(k)γ12(ijg)(k)γ13(ijg)(k)γ23(ijg)(k)ε3(ijg)(k)]T within each k-th lamina, leading to the following expression:

ε(ijg)(k)=∑τ=0N+1∑q=13Z(ijg)(kτ)αiε(ij)(τ)αq(124)

Referring to the generic (α1i,α2j,ζg(k)) point of the shell structure, membrane stresses can be calculated according to the generally anisotropic law of Eq. (51), leading to the following relation, setting E¯ij(kg) for i,j=1,…,6 the elastic stiffness coefficients of the k-th layer referred to the geometric reference system employed for the computation of σ1(k),σ2(k),τ12(k) stresses:

[σ1(ijg)(k)σ2(ijg)(k)τ12(ijg)(k)]=[E¯11(ijg)(k)E¯12(ijg)(k)E¯16(ijg)(k)E¯14(ijg)(k)E¯15(ijg)(k)E¯13(ijg)(k)E¯12(ijg)(k)E¯22(ijg)(k)E¯26(ijg)(k)E¯24(ijg)(k)E¯25(ijg)(k)E¯23(ijg)(k)E¯16(ijg)(k)E¯26(ijg)(k)E¯66(ijg)(k)E¯46(ijg)(k)E¯56(ijg)(k)E¯36(ijg)(k)][ε1(ijg)(k)ε2(ijg)(k)γ12(ijg)(k)γ13(ijg)(k)γ23(ijg)(k)ε3(ijg)(k)](125)

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

[∂∂ζ+1R1+ζ+1R2+ζ000∂∂ζ+1R1+ζ+1R2+ζ000∂∂ζ+1R1+ζ+1R2+ζ][τ13(k)τ23(k)σ3(k)]=[a(k)b(k)c(k)](126)

where R1,R2 are the principal curvature radii of the reference surface r(α1,α2), whereas the terms a(k),b(k),c(k) depend on the quantities calculated in Eq. (125), reading as:

a(k)=−1A 1(1+ζ/R 1)∂σ 1(k) ∂α 1+σ 2(k)−σ 1(k) A 1A 2(1+ζ/R 2)∂A 2∂α 1−1A 2(1+ζ/R 2)∂τ 12(k) ∂α 2−2τ 12(k)A 1A 2(1+ζ/R 1)∂A 1∂α 2b(k)=−1A 2(1+ζ/R 2)∂σ 2(k)∂α 2+σ 1(k)−σ 2(k)A 1A 2(1+ζ/R 1)∂A 1∂α 2−1A 1(1+ζ/R 1)∂τ 12(k) ∂α 1−2τ 12(k)A 1A 2(1+ζ/R 2)∂A 2∂α 1c(k)=−1A 1(1+ζ/R 1)∂τ 13(k)∂α 1−τ 13(k)A 1A 2(1+ζ/R 2)∂A 2∂α 1−1A 2(1+ζ/R 2)∂τ 23(k) ∂α 2−τ 23(k)A 1A 2(1+ζ/R 1)∂A 1∂α 2+σ 1(k)R 1+ζ+σ 2(k)R 2+ζ(127)

It should be noted that the expression of c(k) contains the terms τ13(k),τ23(k). 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 σ1(k),σ2(k),τ12(k) with respect to α1,α2 occurring in the first two expressions of Eq. (127) should be performed employing the GDQ method for each k-th layer of the lamination scheme, setting k=1,…,l:

 σ1,1(ijg)(k)=∂σ 1(k)∂α 1|(ijg)≅∑q=1I Nς iqα 1(1)σ 1(qjg)(k),        σ2,2(ijg)(k)=∂σ 2(k)∂α 2|(ijg)≅∑q=1I Mς jqα 2(1)σ 2(qjg)(k)τ12,1(ijg)(k)=∂τ 12(k)∂α 1|(ijg)≅∑q=1I Nς iqα 1(1)τ 12(qjg)(k),        τ12,2(ijg)(k)=∂τ 12(k)∂α 2|(ijg)≅∑q=1I Mς jqα 2(1)τ 12(qjg)(k)(128)

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

τ~13(ij1)(1)=q1a(ij)(−)τ~23(ij1)(1)=q2a(ij)(−)(129)

being q1a(ij)(−),q2a(ij)(−) the in-plane static loads referred to the (α1i,α2j) point of the bottom skin of the structure following the convention in Eq. (63). On the other hand, for an arbitrary k-th layer, the first two relations of Eq. (126) can be solved at each discrete grid point (α1i,α2j,ζg(k)) by means of the GDQ computations of Eq. (128) enforcing the unknown stress compatibility conditions between the k-th and the (k+1)-th laminae as follows:

τ˜ 13(ijI T)(k)=τ˜ 13(ij1)(k+1)τ˜ 23(ijI T)(k)=τ˜ 23(ij1)(k+1)     for    i=1,...,I N  j=1,...,I M  k=1,...,l−1(130)

For the surface loads q1a(ij)(+) and q2a(ij)(+) 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 [16]:

τ13(ijs)=τ˜13(ijs)+(ζ s+h(ij)2)q 1a(ij)(+)−τ˜13(ijIL)h(ij)τ23(ijs)=τ˜23(ijs)+(ζ s+h(ij)2)q 2a(ij)(+)−τ˜23(ijIL)h(ij)     for    i=1,...,I Nj=1,...,I Ms=1,...,I L=l⋅I T(131)

In the previous equation an index s=kg has been introduced, accounting for g=1,…,IT points defined according to Eq. (122) for the k-th lamina alongside the entire shell thickness h(ij)=h(αi,αj) for i=1,…,IN and j=1,…,IM. As a matter of fact, stresses τ~13(ijIL),τ~23(ijIL) come from the direct solution of the balance Eq. (126). The introduction of the index s leads to vectors τ~13,τ~23 employed for the numerical implementation of Eq. (131), accounting for an assembly of out-of-plane shear stresses alongside the entire lamination scheme:

τ~13=[τ~13(ij1)(1)⋯τ~13(ijIT)(1)τ~13(ij1)(2)⋯τ~13(ijIT)(2)⋯τ~13(ij1)(k)⋯τ~13(ijIT)(k)⋯τ~13(ij1)(l)⋯τ~13(ijIT)(l)]Tτ~23=[τ~23(ij1)(1)⋯τ~23(ijIT)(1)τ~23(ij1)(2)⋯τ~23(ijIT)(2)⋯τ~23(ij1)(k)⋯τ~23(ijIT)(k)⋯τ~23(ij1)(l)⋯τ~23(ijIT)(l)]T(132)

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

∂τ13(k)∂α1|(ijm)≅∑q=1INςiqα1(1)τ13(qjm)(k),∂τ23(k)∂α2|(ijm)≅∑q=1IMςjqα2(1)τ23(qjm)(k)(133)

Now, it is possible to solve the third differential relation in (126), while enforcing the stress compatibility condition for the surface external loads q3a(−) at the boundary of the shell. Referring to the discrete point located at (α1i,α2j), one gets:

σ~3(ij1)(1)=q3a(ij)(−)(134)

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

σ~3(ijIT)(k)=σ~3(ij1)(k+1)(135)

Following a similar procedure to Eq. (132) for τ13(k),τ23(k) shear stresses, the association σ~3(ijm)(k)=σ~3(ijs) is performed, thus leading to the definition of the following vector:

σ~3=[σ~3(ij1)(1)⋯σ~3(ijIT)(1)σ~3(ij1)(2)⋯σ~3(ijIT)(2)⋯σ~3(ij1)(k)⋯σ~3(ijIT)(k)⋯σ~3(ij1)(l)⋯σ~3(ijIT)(l)]T(136)

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

σ3(ijs)=σ~3(ijm)+(ζs+h(ij)2)q3a(ij)(+)−σ~3(ijIL)h(ij)(137)

being q3a(ij)(+) the component of external loads applied at ζ=+h/h22 at each (α1i,α2j) point.

Starting from the main outcomes of Eqs. (131) and (137), the corrected values of the out-of-plane 3D deformations collected γ13(k),γ23(k),ε13(k) of the k-th layer can be calculated as well, according to what has been stated in reference [16]. To this end, the constitutive law for generally anisotropic materials introduced in Eq. (51) is employed. Performing a computation in each discrete point (α1i,α2j,ζg(k)) with i=1,…,IN, j=1,…,IM and g=1,…,IT, it gives the discrete vector x(ijg)(k)=[γ13(ijg)(k)γ23(ijg)(k)ε3(ijg)(k)]T:

A(ijg)(k)x(ijg)(k)=b(ijg)(k)(138)

where

A(ijg)(k)=[E¯44(ijg)(k)E¯45(ijg)(k)E¯34(ijg)(k)E¯45(ijg)(k)E¯55(ijg)(k)E¯35(ijg)(k)E¯34(ijg)(k)E¯35(ijg)(k)E¯33(ijg)(k)]b(ijg)(k)=[τ13(ijg)(k)τ23(ijg)(k)σ3(ijg)(k)]−[E¯14(ijg)(k)E¯24(ijg)(k)E¯46(ijg)(k)E¯15(ijg)(k)E¯25(ijg)(k)E¯56(ijg)(k)E¯13(ijg)(k)E¯23(ijg)(k)E¯36(ijg)(k)][ε1(ijg)(k)ε2(ijg)(k)γ12(ijg)(k)](139)

From the inversion of the fundamental linear system of Eq. (138), the corrected values of γ13(ijg)(k),γ23(ijg)(k),ε3(ijg)(k) 3D strains are provided at each discrete point of the shell (α1i,α2j,ζg(k)) with i=1,…,IN, j=1,…,IM and g=1,…,IT:

γ 13(ijg)(k)=1detA(ijg)(k)(E¯ 33(ijg)(k)E¯ 55(ijg)(k)−(E¯ 35(ijg)(k))2)(τ 13(ijg)(k)−E¯ 14(ijg)(k)ε 1(ijg)(k)−E¯ 24(ijg)(k)ε 2(ijg)(k)−E¯ 46(ijg)(k)γ 12(ijg)(k))              +1detA(ijg)(k)(E¯ 34(ijg)(k)E¯ 35(ijg)(k)−E¯ 33(ijg)(k)E¯ 45(ijg)(k))(τ 23(ijg)(k)−E¯ 15(ijg)(k)ε 1(ijg)(k)−E¯ 25(ijg)(k)ε 2(ijg)(k)−E¯ 56(ijg)(k)γ 12(ijg)(k))              +1detA(ijg)(k)(E¯ 35(ijg)(k)E¯ 45(ijg)(k)−E¯ 34(ijg)(k)E¯ 55(ijg)(k))(σ 3(ijg)(k)−E¯ 13(ijg)(k)ε 1(ijg)(k)−E¯ 23(ijg)(k)ε 2(ijg)(k)−E¯ 36(ijg)(k)γ 12(ijg)(k))γ 23(ijg)(k)=1detA(ijg)(k)(E¯ 34(ijg)(k)E¯ 35(ijg)(k)−E¯ 33(ijg)(k)E¯ 45(ijg)(k))(τ 13(ijg)(k)−E¯ 14(ijg)(k)ε 1(ijg)(k)−E¯ 24(ijg)(k)ε 2(ijg)(k)−E¯ 46(ijg)(k)γ 12(ijg)(k))              +1detA(ijg)(k)(E¯ 33(ijg)(k)E¯ 44(ijg)(k)−(E¯ 34(ijg)(k))2)(τ 23(ijg)(k)−E¯ 15(ijg)(k)ε 1(ijg)(k)−E¯ 25(ijg)(k)ε 2(ijg)(k)−E¯ 56(ijg)(k)γ 12(ijg)(k))              +1detA(ijg)(k)(E¯ 34(ijg)(k)E¯ 45(ijg)(k)−E¯ 35(ijg)(k)E¯ 44(ijg)(k))(σ 3(ijg)(k)−E¯ 13(ijg)(k)ε 1(ijg)(k)−E¯ 23(ijg)(k)ε 2(ijg)(k)−E¯ 36(ijg)(k)γ 12(ijg)(k))ε 3(ijg)(k)=1detA(ijg)(k)(E¯ 35(ijg)(k)E¯ 45(ijg)(k)−E¯ 34(ijg)(k)E¯ 55(ijg)(k))(τ 13(ijg)(k)−E¯ 14(ijg)(k)ε 1(ijg)(k)−E¯ 24(ijg)(k)ε 2(ijg)(k)−E¯ 46(ijg)(k)γ 12(ijg)(k))              +1detA(ijg)(k)(E¯ 34(ijg)(k)E¯ 45(ijg)(k)−E¯ 35(ijg)(k)E¯ 44(ijg)(k))(τ 23(ijg)(k)−E¯ 15(ijg)(k)ε 1(ijg)(k)−E¯ 25(ijg)(k)ε 2(ijg)(k)−E¯ 56(ijg)(k)γ 12(ijg)(k))              +1detA(ijg)(k)(E¯ 44(ijg)(k)E¯ 55(ijg)(k)−(E¯ 45(ijg)(k))2)(σ 3(ijg)(k)−E¯ 13(ijg)(k)ε 1(ijg)(k)−E¯ 23(ijg)(k)ε 2(ijg)(k)−E¯ 36(ijg)(k)γ 12(ijg)(k))(140)

In the last Eq. (140), detA(ijg)(k) refers the determinant of the matrix A(ijg)(k) of Eq. (138):

detA(ijg)(k)=E¯33(ijg)(k)E¯44(ijg)(k)E¯55(ijg)(k)+2E¯34(ijg)(k)E¯35(ijg)(k)E¯45(ijg)(k)−E¯44(ijg)(k)(E¯35(ijg)(k))2−E¯33(ijg)(k)(E¯45(ijg)(k))2−E¯55(ijg)(k)(E¯34(ijg)(k))2(141)

The corrected strain values can be employed for the correction of membrane stresses σ1(ijg)(k),σ2(ijg)(k),τ12(ijg)(k) at each point of the discrete computational domain associated to the structure.

9  Numerical Investigations

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

9.1 Materials Employed in the Analyses

In the present section we describe the mechanical input properties of materials for numerical analyses. For each component, the stiffness matrix E(k) is provided according to conventions from Eq. (49), referring to the material reference system O′α^1(k)α^2(k)ζ^(k). 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 (ρ(k)=7750kg/kgm3m3) has been reported, setting Eij(k) with i,j=1,…,6 the generalized stiffness constants:

E(k)=[E11(k)E12(k)E16(k)E14(k)E15(k)E13(k)E12(k)E22(k)E26(k)E24(k)E25(k)E23(k)E16(k)E26(k)E66(k)E46(k)E56(k)E36(k)E14(k)E24(k)E46(k)E44(k)E45(k)E34(k)E15(k)E25(k)E56(k)E45(k)E55(k)E35(k)E13(k)E23(k)E36(k)E34(k)E35(k)E33(k)]=[98.8453.920.031.05−0.150.7853.9299.190.030.55−0.1850.870.030.0322.55−0.040.250.021.050.55−0.0421.10.071.03−0.1−0.180.250.0721.14−0.1850.7850.870.021.03−0.1887.23]GPa(142)

The trigonal material (ρ(k)=2649kg/kgm3m3) accounts for the following stiffness matrix E(k):

E(k)=[E11(k)E12(k)E