Open Access

ARTICLE

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

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

* Corresponding Author: Francesco Tornabene. 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

**Received** 28 February 2022; **Accepted** 26 April 2022; **Issue published** 31 August 2022

## Abstract

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.## Keywords

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 [1–4]. 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 [7–10]. 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 [13–15]. 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 [16–20]. In particular, two main approaches can be traced, namely the Equivalent Single Layer (ESL) and Layer-Wise (LW) formulations [21–24]. 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 [29–36]. 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 [38–43]. 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 [48–50], 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 [51–53].

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 [56–60]. 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 [65–67], as well as Functionally Graded Materials (FGMs) employing a reduced computational cost [68–71]. 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 [74–76]. 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 [82–83], 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 [89–91], 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 [95–97]. 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

where

As a matter of fact, the association

In this way, a reference surface

where

where

Furthermore, the Lamé parameters

where

Having in mind all these premises, the three-dimensional position vector

being

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

In the previous equation,

where

Setting for simplicity

When arbitrarily-shaped structures are investigated, the fundamental relations, provided in terms of

Accordingly, an inversion of Eq. (16) can be performed if

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

being

Based on Eq. (18), it is possible to express the second order derivatives with respect to the principal coordinate

Eq. (20) assesses the dependence of the second order derivatives with respect to

The first order derivatives

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

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

where

Since a laminated doubly-curved shell structure is considered with

where

As can be seen, the derivatives

If power functions are introduced in Eqs. (27) and (28), the following expressions should be adopted for each

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

If trigonometric functions are adopted in Eq. (27),

Furthermore, the Jacobi orthogonal polynomials

where

setting

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

whereas

In Eq. (37) the differential vectors

Accordingly, coefficients

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

It is useful to introduce, for each

Eventually, Eq. (40) turns into:

being

4 Anisotropic Constitutive LW Relations

We now focus on the elastic constitutive behaviour of a generic doubly-curved laminated structure. Thus, each

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

being

It should be remarked that the constitutive relationship of Eq. (44), expressed for each

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

being

In a more expanded form,

setting sub-matrices

The generic component of Eqs. (50)–(53) are obtained from a through-the-thickness homogenization of the mechanical properties of each

In the previous relation, coefficient

Referring to a generic

In Appendix A, an extended version of the components of the previously-introduced matrix

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.

The present LW formulation considers a two-dimensional structural assessment for each layer of the laminated structure. Accordingly, a generic

Accordingly, if a uniform load is applied to the arbitrary

In addition, a Gaussian function has been implemented in the two-dimensional model, setting

Furthermore, a Super-Elliptic shape of surface loads is introduced so that

As a matter of fact, for

If the Jacobi polynomials

with

The complete expression of

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

Thus, an equivalent surface traction

The surface tractions of Eq. (67) are conveniently arranged in the vector

The function at issue has a singularity for

On the other hand, the numerical modelling of the

When an arbitrary bivariate distribution

where

In this way, the surface integral of

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

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

In the previous equation,

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

Integrating by parts, the following equilibrium relations are derived in terms of

The differential operators

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

where the fundamental operator

The complete expression of

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

being

Note that the interlaminar compatibility conditions of Eq. (82) are modelled in Eq. (83) by means of the identity matrix

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

where

with

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

Referring to the shell sides located at

where

being

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

In the case of constant in-plane distribution of stresses, the relation

where

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

where

being

whereas Eq. (97) assumes the following form:

As a consequence, for each

The generalized boundary stress resultants acting at

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

Fundamental coefficients

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

where

being

For prescribed displacements

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

where

As far as the unified formulation of the displacement field is concerned, a set of

The thickness function matrix

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

where

Accordingly, in the present manuscript the computational domain has been discretized in

Starting from the GDQ rule in Eq. (111), the well-known GIQ method is assessed so that integrations restricted to a generic interval

GIQ weighting coefficients

where

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

being

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

where

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

From the closed interval

In the same way, the discrete form

Starting from the previous equation, it is possible to derive in the arbitrary point located at

At this point, the three-dimensional equilibrium equations of a doubly-curved solid written in curvilinear principal coordinates should be recalled, remembering that

with

From the three-dimensional relations reported in Eq. (123), the out-of-plane stress components

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

The values of

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

In the same way of what defined in Eq. (127), the vector

If the notation

Furthermore, the out-of-plane three-dimensional strain components profile can now be adjusted employing the recovered distributions of out-of-plane stresses

where

The complete procedure for the derivation of the recovered

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

being

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

In Figs. 3b–3d, the sensitivity of the continuous distribution parameters is checked, setting

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

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

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

The orthotropic Graphite-Carbon Epoxy

In addition, the following lattice material named 3D Augmented Re-entrant Cellular Structure (3D ARCS), with density

Among the isotropic materials, a foam of density