Study on the Explicit Formula of the Triangular Flat Shell Element Based on the Analytical Trial Functions for Anisotropy Material

Abstract

This paper presents a novel way to formulate the triangular flat shell element. The basic analytical solutions of membrane and bending plate problem for anisotropy material are studied separately. Combining with the conforming displacement along the sides and hybrid element strategy, the triangular flat shell elements based on the analytical trial functions (ATF) for anisotropy material are formulated. By using the explicit integral formulae of the triangular element, the matrices used in proposed shell element are calculated efficiently. The benchmark examples showed the high accuracy and high efficiency.

1. Introduction

The triangular flat shell element is included widely in finite element software packages. In order to improve the performances of the triangular flat shell element, most of studies focused on two challenges [1–15].

One challenge is the drilling degree of freedom in the membrane element. Introducing the in-plane drilling degree of freedom, the singularity of global stiffness can be avoided in the analysis of shell structures when the neighboring elements around a common node are close to coplanar. Olson introduced the drilling degree of freedom to describe the displacement of the flat shell element [4]. Mohr constructed the hybrid membrane element with the drilling degree of freedom [5]. Allman presented a rational interpolation function to construct the displacement fields associated with the drilling degree of freedom in the triangular membrane elements [6]. Many researchers made different possibilities to improve accuracy of the membrane element by introducing the drilling degree of freedom to the corner nodes of triangular or quadrilateral elements [7–9].

Another challenge is the shear locking of some bending plate elements based on the Mindlin/Reissner thick plate theory [10–15]. The pioneering work of Adini, Clough, and Melosh studied the famous thin plate element ACM, which showed the difficulties in constructing the 𝐶₁ conforming displacement along the boundaries of the thin plate elements [10]. The thick plate element based on the Mindlin/Reissner theory can avoid the strict 𝐶₁ conforming condition [11–15], but the shear locking problem in these elements became the main challenge, which raised the interest of many well-known researchers, including Zienkiewicz et al. [11], Hughes et al. [12], Hinton and Huang [13], Belytschko et al. [14], Bathe and Dvorkin [15].

Long gave a general strategy to deal with the nonconforming displacement along the boundary between the neighboring elements, which was named as the generalized conforming element [16]. Taking advantage of the conforming boundary displacement, the nonconforming displacement model can be used freely to describe the inner field of the element. The generalized conforming equations which are utilized to determine the parameters of the proposed displacement model were studied systematically in the past decades [17–20]. Cen and Long gave a rational way to define the displacement field in the thick plate elements [19]. Fu et al. proved that the hybrid element using the balanced stress fields and the generalized conforming element introducing the associated generalized conforming equations can derive the same transform matrix between the inner parameters and the nodal displacements in the element [20].

Fu et al. and Cen studied the hybrid element based on the variational principle containing the stress function [21]. A series of works showed that elements based on the analytical trial functions (ATF) would be insensitive to the distortion of the element shape [22–25]. The hybrid element strategy was improved through turning to the stress trial functions which are derived from the stress functions (e.g., the Airy stress functions for plane problem). The idea of ATF can be traced back to the first element proposed by Clough which is named as CST (constant stress triangular) element [26]. Some works showed that the element CST can be derived from three rigid body displacements and three displacements of constant stress [1]. On the other hand, the first hybrid element also employed five analytical trial functions of stress [27].

This paper studies the triangular flat shell element based on the analytical trial functions of anisotropy. There are mainly four parts as follows.(a)The first part extracts the general analytical solutions of plane and bending plate problem in anisotropy elasticity (in Section 2).(b)The second part gives the hybrid element strategy based on the analytical trial functions of stress derived in the first part (in Section 3).(c)The third part studies the explicit formulae of the hybrid element based on explicit integral formulae of the triangular element (in Section 4).(d)The fourth part shows some numerical examples as benchmark to study the accuracy and the efficiency of the proposed element model (in Section 5).

The proposed shell element of anisotropy has high accuracy and computational efficiency.

2. The Analytical Trial Functions (ATF)

2.1. The Analytical Trial Functions of the Plane Problem of Anisotropy

The general analytical solutions play very important roles in many methods of computational mechanics [21–27]. Article [28] provided an efficient approach named characteristic differential equation method (CDEM) to derive the general analytical solutions from the governing differential equations directly.

Following the strategy of CDEM to study the plane problem of anisotropy, stress-displacement relationship can be described as 𝜎_𝑥𝜎_𝑦𝜏_𝑥𝑦=𝐶₁₁𝐶₁₂𝐶₁₆𝐶₁₂𝐶₂₂𝐶₂₆𝐶₁₆𝐶₂₆𝐶₆₆𝜕𝑢𝜕𝑥𝜕𝑣𝜕𝑦𝜕𝑣+𝜕𝑥𝜕𝑢𝜕𝑦. (1) And the equilibrium equations in terms of displacement can be written as 𝐶₁₁𝜕²𝑢𝜕𝑥²+𝐶₆₆𝜕²𝑢𝜕𝑦²+2𝐶₁₆𝜕²𝑢𝜕𝑥𝜕𝑦+𝐶₁₆𝜕²𝑣𝜕𝑥²+𝐶₂₆𝜕²𝑣𝜕𝑦²+𝐶₆₆+𝐶₁₂𝜕²𝑣𝐶𝜕𝑥𝜕𝑦=0,₁₆𝜕²𝑢𝜕𝑥²+𝐶₂₆𝜕²𝑢𝜕𝑦²+𝐶₆₆+𝐶₁₂𝜕²𝑢𝜕𝑥𝜕𝑦+𝐶₆₆𝜕²𝑣𝜕𝑥²+𝐶₂₂𝜕²𝑣𝜕𝑦²+2𝐶₂₆𝜕²𝑣𝜕𝑥𝜕𝑦=0. (2)(3) Utilizing the determinant expansion equation of the differential operator matrix in (2), the characteristic differential equation can be given as 𝐶₁₁𝐶₆₆−𝐶²₆₁𝜕⁴𝜑𝜕𝑥⁴+𝐶₂₂𝐶₆₆−𝐶²₆₂𝜕⁴𝜑𝜕𝑦⁴+2𝐶₁₁𝐶₆₂−2𝐶₆₁𝐶₂₁𝜕⁴𝜑𝜕𝑥³+𝜕𝑦2𝐶₆₁𝐶₂₂−2𝐶₆₂𝐶₂₁𝜕⁴𝜑𝜕𝑥𝜕𝑦³+𝐶₁₁𝐶₂₂−2𝐶₆₆𝐶₂₁−𝐶²₁₂+2𝐶₆₁𝐶₆₂×𝜕⁴𝜑𝜕𝑥²𝜕𝑦²=0. (4)

Then the general analytical solutions of the displacement can be derived from the characteristic solution 𝜑 of (4): 𝑢𝑣=𝐶₆₆𝜕²𝜕𝑥²+𝐶₂₂𝜕²𝜕𝑦²+2𝐶₂₆𝜕²−𝐶𝜕𝑥𝜕𝑦₁₆𝜕²𝜕𝑥²+𝐶₂₆𝜕²𝜕𝑦²+𝐶₆₆+𝐶₁₂𝜕²𝜕𝑥𝜕𝑦𝜑. (5) Substituting (5) into (1), the general stress solutions can be derived from the characteristic solutions 𝜑. Table 1 gives seven fundamental analytical stress solutions of the anisotropy plane problem, which can be selected as trial functions of the following numerical method.

Table 1 

The fundamental stress solutions of the anisotropic plane problem.

2.2. The Analytical Trial Functions of the Bending Plate Problem of Anisotropy

The characteristic differential equation of the displacement solutions can be written as [29] −𝑡³𝐶12¹¹𝜕⁴𝑤𝜕𝑥⁴+2𝐶₁₂+4𝐶₆₆𝜕⁴𝑤𝜕𝑥²𝜕𝑦²+4𝐶₆₁𝜕⁴𝑤𝜕𝑥³𝜕𝑦+4𝐶₆₂𝜕⁴𝑤𝜕𝑥𝜕𝑦³+𝐶₂₂𝜕⁴𝑤𝜕𝑦⁴+𝑞𝑥,𝑦=0. (6)

Deriving from the general analytical displacement solutions of (6), Table 2 presents seven general analytical stress solutions of the anisotropy plate problem, which will be selected as trial functions of the following hybrid element in Section 3.

Table 2 

The fundamental stress solutions of anisotropic material used in bending plate problem.

3. The Hybrid Element Strategy Based on the Analytical Trial Functions

3.1. The Hybrid Membrane Element Based on the Analytical Trial Functions of Anisotropy

In the sense of novel hybrid element strategy [21], the basic stress solutions of the anisotropic plane problem in Table 1 can be employed as the trial function of the inner stress field in the element, which is denoted as {𝝈}=𝜎_𝑥𝜎_𝑦𝜏_𝑥𝑦=𝜎_𝑥1𝜎_𝑥2⋯𝜎_𝑥𝑁𝜎_𝑦1𝜎_𝑦2⋯𝜎_𝑦𝑁𝜏_𝑥𝑦1𝜏_𝑥𝑦2⋯𝜏_𝑥𝑦𝑁𝛽₁𝛽₂⋮𝛽_𝑁=𝐒𝜷, (7) where [𝐒] is matrix of the analytical trial functions, whose dimension is 3×𝑁. 𝑁 will be determined by the number of trial functions involved in the element. {𝜷} is the vector of the unknowns of the inner parameters. The modified complementary energy can be presented as [21] Π=Π_𝐶−𝐻_PC, (8) where the complementary energy Π_𝐶 is given by Π_𝐶=12_Ω{𝝈}^𝑇𝐃⁻¹{𝝈}𝑡𝑑𝐴. (9) The constitutive matrix [𝐃] is 𝐃=𝐶₁₁𝐶₁₂𝐶₁₆𝐶₁₂𝐶₂₂𝐶₂₆𝐶₁₆𝐶₂₆𝐶₆₆. (10) Substituting (7) into (9) yields Π_𝐶=12_Ω𝜷^𝑇𝐒^𝑇𝐃⁻¹𝐒𝜷1𝑡𝑑𝐴=2𝜷^𝑇𝐕𝜷, (11) where 𝐕=_Ω𝐒^𝑇𝐃⁻¹𝐒𝑡𝑑𝐴. (12) The additional complementary energy 𝐻_PC is given by 𝐻_PC=_Γ𝑇_𝑥𝑢+𝑇_𝑦𝑣̃𝑡𝑑𝑠, (13) in which 𝑇_𝑥 and 𝑇_𝑦 are the boundary forces along the side Γ around the element, determined by the inner stress field in (7): 𝑇_𝑥𝑇_𝑦=𝜎𝑙0𝑚0𝑚𝑙^𝑥𝜎_𝑦𝜏_𝑥𝑦=𝐋{𝝈}=𝐋𝐒𝜷, (14) where 𝐋=𝑙0𝑚0𝑚𝑙, (15) in which (𝜃)𝑦𝑙=sin=−^𝑖−𝑦_𝑗𝑙_𝑖𝑗(𝜃)=𝑥,𝑚=−cos^𝑖−𝑥_𝑗𝑙_𝑖𝑗. (16)

As showed in Figure 1, 𝜃 is the angle between the side 𝑙_𝑖−𝑗 and the axis 𝑥 and 𝑙_𝑖𝑗 is the length of the side 𝑙_𝑖−𝑗.

Figure 1 

Normal direction and tangential direction.

The conforming displacements 𝑢 and 𝑣̃ along the side 𝑙_𝑖−𝑗 can be obtained from the displacements of 𝑢_𝑛 and 𝑢_𝑠 in Figure 1, which can be assumed as 𝑢_𝑛(𝑠)=𝐻₁(𝑠)𝑢_𝑛𝑖+𝐻₂(𝑠)𝜃_𝑖+𝐻₃(𝑠)𝑢_𝑛𝑗+𝐻₄(𝑠)𝜃_𝑗,𝑢_𝑠(𝑠)=𝐺₁(𝑠)𝑢_𝑠𝑖+𝐺₂(𝑠)𝑢_𝑠𝑗, (17) where 𝑠 is the natural coordinate 𝐿_𝑗 along the side 𝑙_𝑖−𝑗 in the triangular element; its value is defined as 𝑠_𝑖=0 and 𝑠_𝑗=1.

𝐻_𝑖(𝑠) are the Hermite functions, which have the expressions as follows: 𝐻₁(𝑠)=(1−𝑠)²(),𝐻2𝑠+1₂(𝑠)=𝑙_𝑖𝑗(1−𝑠)²𝐻𝑠,₃(𝑠)=𝑠²(),𝐻3−2𝑠₄(𝑠)=𝑙_𝑖𝑗𝑠²(𝑠−1). (18)𝐺_𝑖(𝑠) is the linear shape function 𝐺₁(𝑠)=1−𝑠,𝐺₂(𝑠)=𝑠. (19) Involving the relationship between displacement 𝑢_𝑛/𝑢_𝑠 and displacement 𝑢/𝑣̃, the conforming displacement along the side 𝑙_𝑖−𝑗 can be defined as 𝑢(𝑠)=𝑢_𝑖𝑙²𝐻₁(𝑠)+𝑚²𝐺₁(𝑠)−𝑣_𝑖̃𝑙𝑚𝐺₁(𝑠)−𝐻₁(𝑠+𝜃_𝑖)𝑙𝐻₂(𝑠)+𝑢_𝑗𝑙²𝐻₃(𝑠)+𝑚²𝐺₂(𝑠)−𝑣_𝑗̃𝑙𝑚𝐺₂(𝑠)−𝐻₃(𝑠)+𝜃_𝑗𝑙𝐻₄(𝑠),𝑣(̃𝑠)=−𝑢_𝑖𝑙𝑚𝐺₁(𝑠)−𝐻₁(𝑠)+𝑣_𝑖̃𝑙²𝐺₁(𝑠)+𝑚²𝐻₁(𝑠)+𝜃_𝑖𝑚𝐻₂(𝑠)−𝑢_𝑗𝑙𝑚𝐺₂(𝑠)−𝐻₃(𝑠)+𝑣_𝑗̃𝑙²𝐺₂(𝑠)+𝑚²𝐻₃(𝑠)+𝜃_𝑗𝑚𝐻₄(𝑠). (20)(21) It can also be denoted as 𝑢𝑣̃=𝐍𝐪_𝑖𝑗, (22) where [𝐍] is the matrix of the trial functions of the displacements defined in side 𝑙_𝑖−𝑗 as 𝐍=𝑙²𝐻₁+𝑚²𝐺₁−𝑙𝑚𝐺₁−𝐻₁𝑙𝐻₂𝑙²𝐻₃+𝑚²𝐺₂−𝑙𝑚𝐺₂−𝐻₃𝑙𝐻₄−𝑙𝑚𝐺₁−𝐻₁𝑙²𝐺₁+𝑚²𝐻₁𝑚𝐻₂−𝑙𝑚𝐺₂−𝐻₃𝑙²𝐺₂+𝑚²𝐻₃𝑚𝐻₄. (23)

And {𝐪_𝑖𝑗} is part of the unknown nodal displacement vector, 𝐪_𝑖𝑗=𝑢_𝑖𝑣_𝑖̃𝜃_𝑖𝑢_𝑗𝑣_𝑗̃𝜃_𝑗^𝑇. (24) Substituting (14) and (22) into (13), we have 𝐻_PC=_𝑇{𝜷}^𝑇𝐒^𝑇𝐋^𝑇𝐍{𝐪}𝑡𝑑𝑠={𝜷}^𝑇𝐇{𝐪}, (25) where 𝐇=_𝑇𝐒^𝑇𝐋^𝑇𝐍𝑡𝑑𝑠. (26)

Thus the modified complementary energy Π can be presented as Π=Π_𝐶−𝐻_PC=12𝜷^𝑇𝐕𝜷−𝜷^𝑇𝐇{𝐪}. (27)

With the principle of the modified complementary energy, 𝛿Π=0, 𝜕Π𝜕𝜷={𝟎},𝜕Π𝜕{𝐪}={𝟎}. (28)

From the first part of (28), the relationship between the inner parameters {𝜷} and the nodal displacement {𝐪} can be determined as 𝜷=𝐕⁻¹𝐇{𝐪}. (29)

Substituting (29) into (27), 1Π=−2{𝐪}^𝑇𝐇^𝑇𝐕⁻¹𝐇{𝐪}. (30)

The stiffness matrix can be calculated by 𝐊=𝐇^𝑇𝐕⁻¹𝐇. (31)

Two triangular membrane elements based on the analytical trial functions of anisotropy are studied in this paper. The one containing five items of basic stress solutions in Table 1 is named as ATF-TR5; the other containing seven items is named as ATF-TR7.

3.2. The Hybrid Plate Element Based on the Analytical Trial Functions of Anisotropy

The basic stress solutions of the anisotropic plate problem in Table 2 can be employed as the trial function of the inner stress field in the plate element, which is denoted as {𝐅}=𝑀_𝑥𝑀_𝑦𝑀_𝑥𝑦𝑄_𝑥𝑄_𝑦=𝑀_𝑥1𝑀_𝑥2⋯𝑀_𝑥𝑁𝑀_𝑦1𝑀_𝑦2⋯𝑀_𝑦𝑁𝑀_𝑥𝑦1𝑀_𝑥𝑦2⋯𝑀_𝑥𝑦𝑁𝑄_𝑥1𝑄_𝑥2⋯𝑄_𝑥𝑁𝑄_𝑦1𝑄_𝑦2⋯𝑄_𝑦𝑁𝛽₁𝛽₂⋮𝛽_𝑁=𝐒𝜷, (32) where [𝐒] is matrix of the analytical trial functions, whose dimension is 5×𝑁. 𝑁 will be determined by the number of trial functions involved in the element. {𝜷} is the vector of the unknowns of the inner parameters.

The modified complementary energy of the plate problem can be presented as Π=Π_𝐶−𝐻_PC=_𝐴𝑒12{𝐅}^𝑇𝐀⁻¹{𝐅}𝑑𝑥𝑑𝑦−_𝑆𝑒𝜓_𝑛𝑀_𝑛+𝜓_𝑠𝑀_𝑛𝑠−𝑤𝑄_𝑛=𝑑𝑠_𝐴𝑒12{𝐅}^𝑇𝐀⁻¹{𝐅}𝑑𝑥𝑑𝑦−_𝑆𝑒𝐅^𝑇𝐝𝑑𝑠, (33) where [𝐀] is the constitutive matrix of the anisotropy plate problem and {𝐅} is the vector of the boundary forces along the side of the element; it can be denoted as 𝐅=𝑀_𝑛𝑀_𝑛𝑠−𝑄_𝑛=𝐋𝐅}=𝑙{²𝑚²2𝑙𝑚00−𝑙𝑚𝑙𝑚𝑙²−𝑚²×𝑀00000−𝑙−𝑚_𝑥𝑀_𝑦𝑀_𝑥𝑦𝑄_𝑥𝑄_𝑦=𝐋𝐒𝜷. (34){𝐝} is the conforming displacement; it can be denoted as 𝐝=𝜓_𝑛𝜓_𝑠𝑤. (35)

As showed in Figure 2, the conforming displacement along the side 𝑙_𝑖−𝑗 can be defined as 𝑤_𝑖𝑗=𝐻₁(𝑠)𝑤_𝑖+𝐻₃(𝑠)𝑤_𝑗−𝐻₂(𝑠)𝑥_𝑖𝑗𝜓_𝑥𝑖−𝑦_𝑖𝑗𝜓_𝑦𝑖𝑙_𝑖𝑗−𝐻₄(𝑠)𝑥_𝑖𝑗𝜓_𝑥𝑗−𝑦_𝑖𝑗𝜓_𝑦𝑗𝑙_𝑖𝑗,𝜓_𝑠𝑖𝑗=𝑍₂(𝑠)𝑤_𝑗−𝑤_𝑖𝑙_𝑖𝑗−𝑍₃(𝑠)𝑥_𝑖𝑗𝜓_𝑥𝑖+𝑦_𝑖𝑗𝜓_𝑦𝑖𝑙_𝑖𝑗−𝑍₄(𝑠)𝑥_𝑖𝑗𝜓_𝑥𝑗+𝑦_𝑖𝑗𝜓_𝑦𝑗𝑙_𝑖𝑗,𝜓_𝑛𝑖𝑗=𝐺₁(𝑠)−𝑦_𝑖𝑗𝜓_𝑥𝑖+𝑥_𝑖𝑗𝜓_𝑦𝑖𝑙_𝑖𝑗+𝐺₂(𝑠)−𝑦_𝑖𝑗𝜓_𝑥𝑗+𝑥_𝑖𝑗𝜓_𝑦𝑗𝑙_𝑖𝑗, (36) where 𝑍₂=6𝑠(1−𝑠),𝑍₃=1−4𝑠+3𝑠²,𝑍₄=3𝑠²−2𝑠. (37) And it is defined as 𝑦_𝑖𝑗=𝑦_𝑖−𝑦_𝑗,𝑥_𝑖𝑗=𝑥_𝑖−𝑥_𝑗. (38) Equation (36) can also be denoted as 𝐝_𝑖𝑗=𝜓_𝑛𝑖𝑗𝜓_𝑠𝑖𝑗𝑤_𝑖𝑗=𝐍𝐪_𝑖𝑗, (39) where [𝐍] is the matrix of the trial functions of the displacements defined in side 𝑙_𝑖−𝑗 as 𝐍=𝑦0−^𝑖𝑗𝐺₁𝑙_𝑖𝑗𝑥_𝑖𝑗𝐺₁𝑙_𝑖𝑗𝑦0−^𝑖𝑗𝐺₂𝑙_𝑖𝑗𝑥_𝑖𝑗𝐺₂𝑙_𝑖𝑗−𝑍₂𝑙_𝑖𝑗−𝑥_𝑖𝑗𝑍₃𝑙_𝑖𝑗−𝑦_𝑖𝑗𝑍₃𝑙_𝑖𝑗𝑍₂𝑙_𝑖𝑗−𝑥_𝑖𝑗𝑍₄𝑙_𝑖𝑗−𝑦_𝑖𝑗𝑍₄𝑙_𝑖𝑗𝐻₁−𝑥_𝑖𝑗𝐻₂𝑙_𝑖𝑗𝑦_𝑖𝑗𝐻₂𝑙_𝑖𝑗𝐻₃−𝑥_𝑖𝑗𝐻₄𝑙_𝑖𝑗𝑦_𝑖𝑗𝐻₄𝑙_𝑖𝑗. (40) And {𝐪_𝑖𝑗} is part of the unknown nodal displacement vector: 𝐪_𝑖𝑗=𝑤_𝑖𝜓_𝑥𝑖𝜓_𝑦𝑖𝑤_𝑗𝜓_𝑥𝑗𝜓_𝑦𝑗^𝑇. (41)

Figure 2 

The conforming displacement along the side of the element.

Following the hybrid element strategy proposed in Section 3.1, the modified complementary energy Π can be presented as 1Π=2𝜷^𝑇𝐕𝜷−𝜷^𝑇𝐇{𝐪}. (42)

Similarly, the stiffness matrix can be calculated by 𝐊=𝐇^𝑇𝐕⁻¹𝐇. (43)

Two triangular plate elements based on the analytical trial functions of anisotropy are studied in this paper. The one containing five items of basic stress solutions in Table 2 is named as ATF-TP5; the other containing seven items is named as ATF-TP7.

4. The Explicit Formulae of the Hybrid Element Based on the Analytical Trial Functions

The triangular flat shell elements based on the analytical trial functions of anisotropy can be presented in the natural coordinates of the triangular element. For example, the matrix [𝐒] and, [𝐒] of the stress trial functions contain the variables 𝑥 and 𝑦, which can be written as 𝑥=𝑥_𝑖𝐿_𝑖+𝑥_𝑗𝐿_𝑗+𝑥_𝑚𝐿_𝑚,𝑦=𝑦_𝑖𝐿_𝑖+𝑦_𝑗𝐿_𝑗+𝑦_𝑚𝐿_𝑚, (44) while the matrix [𝐍] and matrix [𝐍] of the conforming displacements contain the variables 𝑠, which can be denoted as the natural coordinate too, that is in the side 𝑙_𝑖−𝑗𝑠=𝐿_𝑗. (45)

According to (16), the matrix [𝐋] and [𝐋] can be expressed in the coordinates of the nodes in the element. For example, in the side 𝑙_𝑖−𝑗, 𝐋=−𝑦_𝑖𝑗𝑙_𝑖𝑗0𝑥_𝑖𝑗𝑙_𝑖𝑗0𝑥_𝑖𝑗𝑙_𝑖𝑗−𝑦_𝑖𝑗𝑙_𝑖𝑗[𝐋]=𝑦_𝑖𝑗,𝑙_𝑖𝑗²𝑥_𝑖𝑗𝑙_𝑖𝑗²𝑥−2^𝑖𝑗𝑦_𝑖𝑗𝑙²_𝑖𝑗𝑥00_𝑖𝑗𝑦_𝑖𝑗𝑙²_𝑖𝑗−𝑥_𝑖𝑗𝑦_𝑖𝑗𝑙²_𝑖𝑗𝑦_𝑖𝑗𝑙_𝑖𝑗²−𝑥_𝑖𝑗𝑙_𝑖𝑗²𝑦00000^𝑖𝑗𝑙_𝑖𝑗−𝑥_𝑖𝑗𝑙_𝑖𝑗. (46)(47)

In utilized equations (44)~(47), all variables in the matrices [𝐇], [𝐕], [𝐊], [𝐇], [𝐕], and [𝐊] can be expressed in the natural coordinates 𝐿_𝑖, 𝐿_𝑗, and 𝐿_𝑘. To the integral formula of the natural coordinates in the triangular element, we have _𝐴𝐿^𝑎_𝑖𝐿^𝑏_𝑗𝐿^𝑐_𝑚𝑑𝑥𝑑𝑦=𝑎!𝑏!𝑐!(𝑎+𝑏+𝑐+2)!2𝐴,_𝑙𝐿^𝑎_𝑖𝐿^𝑏_𝑗𝑑𝑠=𝑎!𝑏!()!𝑙𝑎+𝑏+1^𝑖𝑗, (48) where 𝐴 is the area of the triangular element. Utilizing the explicit integral formulae of (48) in the triangular element, the explicit formula of the proposed elements in this paper can be obtained.

For example, in the plane elements ATF-TR5 and ATF-TR7, we have 𝑉()=𝐶1,1𝐴𝑡₂₂𝐶₆₆−𝐶²₂₆Δ_𝐷,𝑉()𝐶1,2=−𝐴𝑡₁₂𝐶₆₆−𝐶₁₆𝐶₂₆Δ_𝐷,𝑉()=𝐶2,2𝐴𝑡₁₁𝐶₆₆−𝐶²₁₆Δ_𝐷,𝑉()=𝐶1,3𝐴𝑡₁₂𝐶₂₆−𝐶₁₆𝐶₂₂Δ_𝐷,𝑉()𝐶2,3=−𝐴𝑡₁₁𝐶₂₆−𝐶₁₂𝐶₁₆Δ_𝐷,𝑉()=𝐶3,3𝐴𝑡₁₁𝐶₂₂−𝐶²₁₂Δ_𝐷,𝑉()=𝑡𝑥1,4^𝑗+𝑥_𝑚𝐶₂₂𝐶₆₆−𝐶²₂₆𝐶₆₆3,𝑉()=𝑦2,4𝐴𝑡_𝑗+𝑦_𝑚𝐶₁₁𝐶₆₆−𝐶²₁₆𝐶₂₆3,Δ_𝐷=𝐃=𝐶₁₁𝐶₂₂𝐶₆₆−𝐶₁₁𝐶²₂₆−𝐶²₁₂𝐶₆₆+2𝐶₁₂𝐶₁₆𝐶₂₆−𝐶²₁₆𝐶₂₂. (49)

In the plate elements ATF-TP5 and ATF-TP7, we have 𝑉(1,1)=12𝐴_𝑡𝐷𝐶₂₂𝐶₆₆−𝐶²₂₆,𝑉(1,2)=−12𝐴_𝑡𝐷𝐶₁₂𝐶₆₆−𝐶₁₆𝐶₂₆,𝑉(2,2)=12𝐴_𝑡𝐷𝐶₁₁𝐶₆₆−𝐶²₁₆,𝑉(1,3)=12𝐴_𝑡𝐷𝐶₁₂𝐶₂₆−𝐶₁₆𝐶₂₂,𝑉(2,3)=−12𝐴_𝑡𝐷𝐶₁₁𝐶₂₆−𝐶₁₂𝐶₁₆,𝑉(3,3)=12𝐴_𝑡𝐷𝐶₁₁𝐶₂₂−𝐶²₁₂,𝑉()=𝑥1,44𝐴_𝑗+𝑥_𝑚𝑡³,𝑉()𝐴2,4=0,_𝑡𝐷=𝐴/𝑡³Δ_𝐷. (50)

5. Numerical Examples

5.1. Cantilever Shell

As shown in Figure 3, a slender cantilever shell has the length 𝑏 of 32 m, the height ℎ of 2 m, and the thickness 𝑡 of 0.01 m. It is made of the anisotropy T700, which has the material parameters 𝐶₁₁=137.54 GPa, 𝐶₁₂=4.14 GPa, 𝐶₂₂=9.79 GPa, and 𝐶₆₆=4.47 GPa, 𝐶₄₄=4.47 GPa, 𝐶₅₅=4.47 GPa. There is an angle 𝜃 between the principle axis 1 and the axis 𝑦, designated as 0°, 45°, and 90°. Figure 4 shows four schemes of the mesh to the proposed cantilever shell.

Figure 3 

Cantilever beam subjected to force/moment.

(a)

(b)

(c)

(d)

Figure 4 

Schemes of the mesh.

Two load conditions are studied in the example. One loading is the unit in-plane moment 𝑀_𝑥=1 N·m on the end side 𝐶. Table 3 gives the results of the displacement 𝑉_𝐶𝑧 obtained from the proposed elements ATF-TR5 and ATF-TR7, respectively. Comparing with the element S3, which is employed in the software ABAQUS, ATF-TR5 and ATF-TR7 are more accurate.

Table 3 

Deflection 𝑉_𝐶𝑧/10⁻⁶ m of the cantilever shell subjected to 𝑀_𝑥.

The other loading is the unit out-plane load 𝑃_𝑥=1 N on the end point 𝐶; Table 4 gives the results of the displacement 𝑉_𝐶𝑥 obtained from the proposed element, ATF-TP5 and ATF-TP7. Comparing with the element S3, ATF-TP5 and ATF-TP7 are more accurate.

Table 4 

Deflection 𝑉_𝐶𝑥/m of the cantilever shell subjected to 𝑃_𝑥.

5.2. Clamped Square Plate

As shown in Figure 5, a clamped square plate has the length of 1 m in each side and the thicknesses 𝑡 of 0.001 m, 0.01 m, and 0.1 m. It is also made of the anisotropy T700 and is subjected to the uniformly distributed loading 𝑞=1 Pa. Table 5 gives the deflection 𝑤_𝐶 of the central point 𝐶 in different mesh schemes.

Table 5 

Deflections 𝑤_𝐶 of the central point C.

Figure 5 

Schemes of the mesh.

5.3. Cook’s Skew Beam

As showed in Figure 6, Cook’s skew flat shell [30] is studied. The thickness 𝑡 of the shell is 1 m. It is also made of the anisotropy T700 and is subjected to the uniformly distributed load which has the summation of 𝑃=1 N. Table 6 gives the deflection of the central point 𝐶 in different mesh schemes.

Table 6 

Deflections (10⁻⁹ m) of the central point C in different mesh scheme.

Figure 6 

Cook’s skew shell.

6. Conclusions

This paper presents a novel way to formulate the flat shell element of anisotropy. The basic analytical solutions and the proposed conforming displacement are employed to formulate the triangular elements ATF-TR5, ATF-TR7, ATF-TP5, and ATF-TP7. The high accuracy and computational efficiency of these elements are proved by the benchmark.

Acknowledgments

This project is supported by the National Natural Science Foundation of China (no. 11272340) and the National Basic Research Programs of China (no. 2010CB731503).

Supplementary Materials

The results of ATF-TR5 and ATF-TR7 in Table 3 are recalculated using the analytical stress solutions of isotropy in Table 3(b). The results in Table 3(a) show that the element ATF-TR5A and ATF-TR7A, which are based on the analytical stress solutions of isotropy, can perform as well as the element ATF-TR5 and ATF-TR7. All of the elements based on the ATF perform better than S3.

1. Supplementary Table