A comparison between different finite elements for elastic and aero-elastic analyses.

In the present paper, a comparison between five different shell finite elements, including the Linear Triangular Element, Linear Quadrilateral Element, Linear Quadrilateral Element based on deformation modes, 8-node Quadrilateral Element, and 9-Node Quadrilateral Element was presented. The shape functions and the element equations related to each element were presented through a detailed mathematical formulation. Additionally, the Jacobian matrix for the second order derivatives was simplified and used to derive each element's strain-displacement matrix in bending. The elements were compared using carefully selected elastic and aero-elastic bench mark problems, regarding the number of elements needed to reach convergence, the resulting accuracy, and the needed computation time. The best suitable element for elastic free vibration analysis was found to be the Linear Quadrilateral Element with deformation-based shape functions, whereas the most suitable element for stress analysis was the 8-Node Quadrilateral Element, and the most suitable element for aero-elastic analysis was the 9-Node Quadrilateral Element. Although the linear triangular element was the last choice for modal and stress analyses, it establishes more accurate results in aero-elastic analyses, however, with much longer computation time. Additionally, the nine-node quadrilateral element was found to be the best choice for laminated composite plates analysis.

coefficient matrix for in-plane action

steady aerodynamic coefficient matrix in structural coordinates

unsteady aerodynamic coefficient matrix in structural coordinates

steady aerodynamic coefficient matrix

strain-displacement matrix

displacement in global coordinates

stress-strain matrix (isotropic material properties matrix)

Jacobian Matrix for first order derivatives

stiffness matrix

mass matrix

shape function matrix

structural bending displacement field

structural bending nodal displacements

element local coordinates

displacement vector in local coordinates

structural global coordinates

volume

elasticity modulus of the wing material

flow speed

dynamic pressure

stress

reference length (half the wing root chord)

Jacobian matrix for second order derivatives

reduced frequency

plate wing thickness

wing damping ratio

flutter frequency

strain vector

reference element coordinates

air density

Introduction

Numerical methods are usually the first choice for many researchers and engineers to analyze complicated systems because of their accessibility, flexibility and ability to solve complex systems. The Finite Element Method (FEM) as one of the powerful numerical methods for structural analysis comes at the top of the list of all numerical methods. As introduced in many Refs. [1], [2], [3], [4], [5], [6], the method is mainly based on dividing the whole structure into a finite number of elements connected at nodes. The properties of the whole structure such as mass and stiffness, which are continuous in nature, are discretized over the elements and approximate solutions are obtained for the governing equations. The elements equations are assembled together to reach a global system of algebraic equations, which can be solved for the unknown solution variables of the structure. The accuracy of the FEM solution depends on many factors, such as the interpolation polynomials and subsequently the element shape functions, the number of degrees of freedoms selected for each element, the mesh size, and the type of element used. The model accuracy is a result of the deep understanding of the effect of each factor on the final results.

The selection of the element interpolation functions is a key factor in the accuracy of the FEM solution. For this reason, intensive researches have been made to develop new finite elements having different shapes and interpolation functions. There are numerous types of elements for different structural problems. In this paper, the main focus is on two-dimensional shell elements. Finite shell elements such as triangular elements [7], [8], [9], quadrilateral elements [10], [11], higher order elements [12], [13], [14], [15], [16], [17], and improved elements [18] are all tested and approved to achieve an acceptable level of accuracy. Although a vast number of elements are available in literature, researchers cannot easily figure out which element is the most suitable to select for their particular problem. The selection problem is even more difficult for engineers who are mainly interested in the application rather than the theoretical background. Additionally, the detailed mathematical formulation of some thin shell bending elements, especially the higher order ones, cannot be easily found in the literature.

Considering aero-elasticity in which the structural model is coupled to an aerodynamic model adds more complications to the problem, and makes the choice of the suitable element more challenging. Aero-elasticity is crucial for structures such as aircraft, wind turbines, and several other applications in which divergence and flutter phenomena may occur leading to catastrophic failures of the structure. Therefore, designers of these structures are constrained by the design limits and definitely need accurate FEM without being computationally expensive.

Therefore, the aim of the present work is to present a detailed mathematical formulation for different thin shell finite elements along with a complete comparison between them for specific problems in structures and aero-elasticity. The results of the selected elements are compared based on (1) solution accuracy of each element, (2) number of elements needed to achieve convergence, and (3) computational time. The comparison is for free vibration analysis, stress analysis, aero-elastic analysis, and laminated composite analysis. Five different elements are selected for the present comparison with different nature. These finite elements are

Three-node linear triangular element [1] denoted as LINTRI.

Four-node linear quadrilateral element [1] denoted as LINQUAD.

Four-node linear quadrilateral element based on deformation modes (MKQ12 [18]).

Eight-node quadrilateral element denoted as QUAD8NOD.

Nine-node quadrilateral element denoted as QUAD9NOD.

These elements are selected with different nature ranging from linear to higher order, triangular to quadrilateral, and improved to regular elements to provide wide range of variety to the present comparison. All these elements are tested using bench mark problems from the literature [19], [20] for elastic and aero-elastic analyses with analytical results and/or experimental measurements. The element shape functions are derived using MATHEMATICA [21] software and then implemented into MATLAB [22] codes to solve the selected problems.

The finite elements’ formulation

The present finite element model is based on either the classical plate theory for metallic materials or laminated plate theory for composite materials. Both are based on the Kirchhoff assumptions which neglect the transverse shear and transverse normal effects [2].

To formulate a finite shell element there is a standard procedure that is usually followed.

Start from the weak (integral) form of the governing equation.

Assume suitable interpolation polynomials for both the in-plane P and bending P displacement fields.

Calculate the coefficients of these polynomials by applying the nodal movement conditions.

Determine the shape functions for both the in-plane N and bending N actions.

Derive the strain-displacement matrix B from the shape functions’ derivatives.

Integrate to obtain the element stiffness matrix K knowing the material elasticity matrix and the strain displacement relationships.

Calculate the element mass matrix M from the element shape functions and the material density ρ, and finally,

The structural matrices K and M can be obtained by assembling the element matrices obtained in steps 6 and 7.

All these steps were followed for each element considered in the current study to derive the element shape functions, strain-displacement relationships, stiffness, and mass matrices for both the in-plane and bending actions. The element shape functions, presented in this section, are derived using MATHEMATICA software. All the strain-displacement matrices, stiffness, and mass matrices are numerically integrated using MATLAB software.

General formulation

In the present section, general formulation of the element shape functions and strain displacement matrices is developed. Based on this formulation all the shell elements’ shape functions and subsequently the elements’ equations are derived

In-plane action

An interpolation function is chosen from Pascal’s Triangle [1], [2],where u is the in-plane displacement field at any point through the element and a is a vector of constants to be determined from the nodal in-plane displacements U.nnod represents the total number of nodes in an element. The size of U equals the total element in-plane degrees of freedoms.

Finally, the in-plane shape functions can be obtained,

The in-plane strain displacement matrix can be obtained from the shape functions’ derivatives by using the Jacobian matrix definitionwhere are the elements of the Jacobian matrix and is the Jacobian determinant

Bending action

An interpolation function is chosen either from Pascal’s Triangle or based on the displacements modes,where w is the bending displacement at any point on the element, from which we can obtain the rotation around the x-axis (θ) and the y-axis (θ) using the Jacobian matrix.

Note that the Jacobian matrix elements are rearranged in so that the displacement rotations are defined as,

Then, the three bending displacements can be calculated from the equation,a is the coefficients vector to be determined from the out of plane nodal displacements W. The bending shape functions will have the formwhere the C matrix is calculated from the c matrix after applying the nodal boundary conditions, and

Then the strain-displacement matrix can be derived and simplified from the Jacobian definition for second order derivatives. They were derived and simplified by the authors to have the formwhere is the Jacobian Matrix for the second order derivatives, which can be calculated using the elements of the regular Jacobian matrix to have the form,

The inverse of the Jacobian Matrix for the second order derivatives is

Based on this simple and detailed mathematical implementation, the considered elements’ equations can be derived. All the shape functions for those elements are presented in the following sections.

Notice that and are represented as row vectors all over the present paper.

The linear triangular element (LINTRI)

The LINTRI thin-shell element has three nodes. The element has six degrees of freedom per node with a total of 18 degrees of freedom. Fig. 1a shows a schematic of the element with the element global, local, and reference coordinates. The element interpolation and shape functions are derived in the following.

Fig. 1

Finite elements local and reference coordinates.

For in-plane action

The interpolation polynomial for in-plane action has the formand subsequently the in-plane shape functions have the form

For bending action

The interpolation polynomial for bending action based on the element area coordinates has the formand subsequently the bending shape functions are

The linear quadrilateral element (LINQUAD)

The LINQUAD element consists of four nodes. It has six degrees of freedom per node with a total of 24 degrees of freedom. Fig. 1b shows a schematic of the element with the global, local, and reference coordinates. The element interpolation and shape functions are derived to be as follows.

The interpolation polynomial for in-plane action selected from Pascal’s Triangle has the formand subsequently the in-plane shape functions have the form

The interpolation basis functions for bending action selected from Pascal’s Triangle has the formand subsequently the bending shape functions are

The linear quadrilateral element based on deformation modes (MKQ12)

The MKQ12 element has four nodes. It has six degrees of freedom per node with a total of 24 degrees of freedom. It has the global, local, and reference coordinates shown in Fig. 1b. This element was introduced by Karkon and Rezaiee-Pajand [18]. It has the same in-plane shape functions of the LINQUAD element but with improved bending shape functions based on the deformation modes.

The shape functions are then

The eight-node quadrilateral element (QUAD8NOD)

The QUAD8NOD element has eight nodes. It has six degrees of freedom per node with a total of 48 degrees of freedom. Fig. 1c shows a schematic of the element with the global, local, and reference coordinates. The element interpolation and shape functions were derived in the following.

The interpolation polynomial for in-plane action selected from Pascal’s Triangle has the formand subsequently the in-plane shape functions have the form

The interpolation polynomial for bending action was selected carefully from the well-known Pascal Triangle. Initially, the following basis functions were selected;

Using the above basis functions yields a singular C matrix. The rank of the matrix turns out to be 22 instead of 24, which indicates that two terms result in repeated equations. Different terms have been replaced with higher order terms to detect the reason for the singularity. The analysis revealed that the bilinear term , the biquadratic term , and the bicubic term all yield similar equations. Therefore, the biquadratic and bicubic terms were replaced with and . Finally, the bending interpolation function for the QUAD8NOD element is;

This choice eliminates all the singularities, and subsequently the bending shape functions are

The nine-node quadrilateral element (QUAD9NOD)

The QUAD9NOD element has nine nodes. It has six degrees of freedom per node with a total of 54 degrees of freedom. It has the global, local, and reference coordinates shown in Fig. 1d. The element interpolation and shape functions were derived in the following.

The interpolation polynomial for in-plane action selected from Pascal’s Triangle has the formand subsequently the in-plane shape functions have the form

The interpolation polynomial for bending action was selected carefully from the well-known Pascal Triangle. Initially, the following basis functions were selected;

Using the above basis functions yielded singular C matrix. The rank of the matrix turned out to be 23 instead of 24, which indicated that two terms result in repeated equations. Different terms have been replaced with higher order terms to detect the origin of the singularity. Finally, the bending interpolation function for the QUAD9NOD element is;and subsequently the bending shape functions are

The test problems

It has been mentioned earlier that the aim of the present work was to compare between different shell finite elements with different behavior for elastic and aero-elastic analyses. This will enable any researcher to select the shell element which best suits his/her specific application. Different test benchmark problems are considered for elastic and aero-elastic analyses. These problems are described in detail in this section in addition to their mathematical models. These models are implemented in the next section into MATLAB codes, which are carefully constructed and validated. For each problem, a suitable number of elements was selected based on convergence analysis. The number of elements was increased till the response converged to a certain value. Then the results were compared with published experimental or analytical solutions.

Study of the natural frequencies of free vibration of an elastic square plate

Problem formulation

This problem was presented by Safizadeh et al. [20] in which an analytical solution was provided. It deals with a square plate (1 m × 1 m) with thickness t equal 0.003 m and the material properties, elastic modulus, mass density, and Poisson’s ratio

The plate is fixed along all sides.

The mathematical model

The plate natural frequencies are calculated by solving the eigenvalue problemwhere the subscript b refers to the bending action and the mass and stiffness matrices are calculated from the strain displacement matrix and the shape functions

The superscript e means that these matrices are calculated over each element and then assembled in the global coordinates.

D is the isotropic material stiffness matrix

Stress and deformation analysis of metallic plate wing

The problem formulation

The elastic stress and displacement analysis performance of the considered elements was tested by analyzing a plate-like straight-rectangular wing under aerodynamic load. The wing geometry was shown in Fig. 2. The aerodynamic analysis was performed by using the Doublet lattice method [23]. A convergence analysis was performed to select the suitable number of elements for both the aerodynamic and finite element analyses.

Fig. 2

Plate wing plane form geometry.

The material properties are

Young’s Modulus = 98E9 Pa, Poison’s ratio = 0.28, plate thickness = 0.001 m

The flow properties are

Speed = 30 m/s, density = 1.225 kg/m3, AOA = 3°

AOA is the flow angle of attack.

First the stiffness matrix for both the in-plane and bending actions are derived. The stiffness matrix of the bending action is given in Eq. (40) and the in-plane stiffness matrix has the form

The elastic problem was solved to find the displacements, then the stresses are calculated using the equationwhere is the local in-plane displacements vector and is the local bending displacements vector.

Wing aero-elastic analysis

The finite element selection in aero-elastic analysis is always a problem. The point is to select the suitable element for accurate aero-elastic calculations and load transformation between the aerodynamic model and the structural model and vice versa. For this purpose, a comparison was made between different finite elements for aero-elastic analysis.

The results are compared with published experimental results. For all the finite elements, the shape functions are used in the aero-elastic coupling, rather than the conventional spline interpolation, to make the model more accurate and consistent [24].

A straight-rectangular plate wing model made of laminated composite materials was analyzed with different laminate configurations. The wing has the same plane form shown in Fig. 2. The suitable number of elements was selected for both the aerodynamic and finite element analyses throughout convergence analyses. The lamina material properties were EL = 98 Gpa, ET = 7.9 Gpa, GLT = 5.6 Gpa,  = 0.28, ρ = 1520 kg/m3, t = 0.134e−3m

Mathematical model

The need to decrease the aircraft structural weight for economic purposes leads to an increase in the aircraft flexibility, and subsequently the tendency for aero-elastic instability. The wing aero-elastic instability was categorized into divergence and flutter analyses [25], [26], [27], [28], [29]. In the former, analysts were interested in determining the minimum speed at which wing static torsional instability takes place. In the latter, analysts were interested in determining the minimum speed at which wing dynamic instability flutter takes place. For both analyses, the doublet lattice method was used for the steady and unsteady aerodynamic analyses [23] for all elements. An exception was the linear triangular element where the vortex lattice method [30] was used in the steady aerodynamic analysis because it produces more accurate results, although it needs more elements and subsequently longer computation time. For this reason, it was not considered for the rest of elements, as the doublet lattice method was enough and smaller computational time.

The problem was solved by developing two models; one for the structural analysis and the other for the aerodynamic analysis. Then, the aerodynamic coefficient or stiffness matrix was transformed into the structural nodes by means of either spline interpolation or by the same shape functions of the finite element. The use of the shape functions of the finite element in the connection between the finite element model and the aerodynamic model was found to be more accurate and consistent than the spline method [24]. The mathematical models for both the divergence and flutter analyses are presented in this section.

The divergence analysis

Divergence can be regarded as a static torsional instability that occurs for aircraft wings at a certain flight speed. The divergence speed can be calculated by solving the eigenvalue problemwhere A is the aerodynamic stiffness transformed from the aerodynamic control points into the finite element nodes by using the element shape functions [31], and represents the dynamic pressure at which divergence takes place.where is the velocity at which the divergence occurs.

The flutter analysis

Flutter can be regarded as a dynamic instability that occurs to aircraft wings at a certain flight speed. The flutter speed was determined by solving the eigenvalue problemwhere b is a reference length (chosen to be half the wing root chord), A is the unsteady aerodynamic stiffness matrix transformed from the aerodynamic control points to the finite elements nodes by using the element shape functions [31], while k is known as the reduced frequency which is defined as

Eq. (46) comprises two unknowns; the speed V and the frequency ω, and both can be obtained by iteration where the flutter occurs at zero damping coefficient ς.

It is worth noting that Eq. (46) is nested and solved using the k-method.

Aerodynamic analysis and aero-elastic coupling

The aerodynamic coefficient matrix A for steady aerodynamic analysis was calculated using either the Vortex Lattice Method (VLM) [30] or the Doublet Lattice Method (DLM) [23]. Then, they were transformed to the structural coordinates as followingwhere Asteady is the steady aerodynamic coefficient matrix at the aerodynamic control points. GN1 and GN2d are transformation matrices calculated from the element bending shape functions.T and are geometric transformation matrices that connect between the global structural coordinates and the element local coordinates.n represents the number of aerodynamic control points in each element.

In flutter analysis, the unsteady aerodynamic coefficient matrix A was calculated using the DLM.Aunsteady is the unsteady aerodynamic coefficient matrix at the aerodynamic control points. GN2f is calculated as GN2d, by considering the lateral displacement . Bc is a boundary conditions matrix calculated as

Laminated plate elastic analysis

To present a complete picture regarding the differences between the considered finite elements, a comparison between them for the analysis of a composite laminated plate was established. A square plate was considered with 25 cm side length and 1 cm total thickness [32]. The lamina material properties were EL = 52.5 MPa, ET = 2.1 MPa, GLT = 1.05 MPa,  = 0.25. Reddy [32] presents an analytical solution for the maximum displacement of the square plate subject to a distributed pressure of 1 N/cm2. Two laminate configurations were considered as well as two different boundary conditions. The analytical results are represented in the following section.

The mathematical model for this problem is exactly the same as that of the elastic deformation problem, but with imposing the effect of composite material on the stress-strain relation. More details can be found in Reddy [32].

Results and discussion

The selected finite elements are tested by the four problems described in the previous section, and the results are demonstrated in this section. All the analyses have been performed on a personal computer with an i7 processor CPU @ 3.6 GHz, intel core and 16 GB RAM. The results for each analysis are shown in the next subsection.

The dynamic elastic analysis

The problem of the square plate presented in the previous section was implemented on a computer code using MATLAB software. Table 1 shows the predicted natural frequencies. The first column contains the analytical solution [20] for the first five natural frequencies, while the finite element results are listed in the rest of the columns. Results have shown that the natural frequencies predicted by the different types of finite elements are in general less than those predicted by the analytical model. The frequency values obtained by using the Linear Triangular Element are farthest from the analytical ones and the number of elements (N_elem) needed to reach convergence is a maximum. On the other hand, the frequencies obtained by using the Linear Quadrilateral Element based on deformation modes are closest to the analytical solution. The MKQ12 performed even better than the higher order elements, which were expected to produce accurate results in bending analysis. The number of elements needed to reach convergence in the linear quadrilateral element based on deformation modes is also the same as those of higher order elements.

Table 1

The natural frequencies of clamped square plate [Hz].

Mode	Analytical (20)	LINTRI	LINQUAD	MKQ12	QUAD8NO	QUAD9NOD
1	28.47	26.2	26.54	26.81	26.68	26.61
2	58.06	53.35	54.06	54.78	54.39	54.27
3	58.06	53.35	54.06	54.78	54.39	54.25
4	85.66	78.35	79.09	81.49	80.3	80.01
5	104.09	95.54	96.91	98.29	97.42	97.24
Nelem	–	800	225	100	100	100

The average error percentage and processing time of each element are demonstrated in Fig. 3. The product of the average error percent and processing time can be used as a measure of excellence of the finite element, where the best element has the minimum value for the product. This product is given in Table 2, which shows that the linear quadrilateral element based on deformation modes is the best element to use in this type of problems.

Fig. 3

The average error and execution time of each finite element in the frequency analysis.

Table 2

The product of the average error percentage and processing time for various finite elements.

Mode	LINTRI	LINQUAD	MKQ12	QUAD8NO	QUAD9NOD
∊avg × Time	125.5	23.8	7.7	25.8	101

Plate wing stress and displacement analysis

The elastic performance of a plate wing under steady aerodynamic load was studied using the five finite elements under consideration. The values of the maximum Von Mises stresses and maximum displacements are tabulated in Table 3 and plotted in Fig. 4 together with the execution time. Fig. 5 shows the distribution of the Von Mises stress in the wing for each element type. The stresses are calculated in the case of triangular elements over the mid-side points and then averaged over the element. In case of the quadrilateral elements, the stresses are determined at the element integration points, and then averaged over the element.

Table 3

Max. displacement and stress over the plate wing under aerodynamic load.

Element	LINTRI	LINQUAD	MKQ12	QUAD8NO	QUAD9NOD
dmax [mm]	19.1	20.3	20.3	20.4	19.9
σmax [MPa]	28.9	32.9	33.2	34.4	34.1
Nelem	192	120	120	60	60
Time [s]	20.4	5.2	8.1	6	7.9

Fig. 4

Max. displacement and stress of the plate wing in addition to the executing time.

Fig. 5

The Von Mises stresses for each element model.

Since there was neither analytical nor experimental data available for this model, the error percentage cannot be computed for this particular problem. However, it is clear from Fig. 5 that the stress distribution resulting from using higher order elements (QUAD8NOD and QUAD9NOD) is the smoothest and most realistic. This can be attributed to the higher order interpolation functions for displacements, which render the stress distribution (derived from the displacement derivatives) continuous. Hence, if all the results are considered together, the best performing element can be considered to be the QUAD8NOD element, which results in accurate displacement and stress distributions, together with a reasonable computational time. Following this element comes the QUAD9NOD in the second place. On the other hand, the LINTRI element comes as the worst element for wing stress analysis from the point of view of computation time and stress distribution as seen in Fig. 5.

Composite plate wing aero-elastic analysis

The static and dynamic aero-elastic analysis of a composite plate wing was carried out using the five elements. The results are listed in Table 4 for the smaller of the divergence and Flutter speeds. The wing aero-elastic analysis was performed for different laminate configurations. The subscript D refers to the Divergence speed while subscript F refers to the Flutter speed. The error percent, the average error, and the computation time are listed in Table 5 and plotted in Fig. 6. Considering the minimum value of the product of the average error and execution time to be the sign of excellence, we find that the QUAD9NOD element was the best choice for wing aero-elastic analysis, followed by the QUAD8NOD element.

Table 4

Divergence speed of the plate wing related to different laminate configurations.

Laminate configuration	Divergence/Flutter Speed [m/s]
	Exp [19]	LINTRI [31]	LINQUAD	MKQ12	QUAD8NOD	QUAD9NOD
[0 0 90 90 0 0]	25F	25.4/26.4	25.4/24.47	25.5/25.9	25.4F	52.6/25.3
[45 −45 0 0 −45 45]	>32	47.5F	43.8F	46.8F	45.6F	46.6F
[45 45 0 0 45 45]	28F	27.8F	26.1F	27.6F	29.2F	29F
[−45 −45 0 0 −45 −45]	12.5D	12.7/29.1	11.4/26.9	11.5/29	12.86/23.5	12.88/32.4
[30 30 0 0 30 30]	27F	27.4F	26.1F	27.16F	28.1F	27.8F
[−30 −30 0 0 −30 −30]	11.7D	12.8/48.1	11.58/33.7	11.67/35.09	13/31.4	13/30.6
Nelem	–	48	96	96	12	12
FE mesh	–	3 × 8	6 × 16	6 × 16	2 × 6	2 × 6
DLM mesh	–	8 × 12	6 × 12	6 × 12	7 × 12	6 × 12

Table 5

The error % in each analysis and the computation time.

Laminate configuration	Error %
	LINTRI	LINQUAD	MKQ12	QUAD8NOD	QUAD9NOD
[0 0 90 90 0 0]	5.6	2.1	3.7	1.6	1.3
[45 −45 0 0 −45 45]	–	–	–	–	–
[45 45 0 0 45 45]	0.7	6.8	1.4	4.4	3.6
[−45 −45 0 0 −45 −45]	1.6	8.8	7.7	2.9	3.1
[30 30 0 0 30 30]	1.5	3.3	0.6	4.2	3.2
[−30 −30 0 0 −30 −30]	5.6	2.1	3.7	1.6	1.3
∊avg %	2.3	5.2	3.4	3.3	2.8
Time [s]	120	20.3	22.8	17.6	18.6
∊avg × Time	276	106	77.5	58.1	52.1

Fig. 6

The average error and computation time for each element.

Laminated plate elastic analysis

The square laminated plate was analyzed for [0, 90]o and [−45, 45]° laminate configurations. Two boundary conditions were considered; in the first all the plate sides were simply supported, and in the second all the plate sides were clamped [32]. The analytical results are listed in Table 6, and the average error and computation time for each element are depicted in Fig. 7. The analyses are obtained for the maximum normalized bending displacement,where is the plate maximum thickness, is the total laminate thickness, L is the side length, and P is the applied pressure load.

Table 6

Laminated elastic plate maximum normalized displacement.

BC's	Laminate config	Analytical	LINTRI	LINQUAD	MKQ12	QUAD8NOD	QUAD9NOD
Simply supported	[0, 90]	1.6955	1.719	1.61	1.606	1.6958	1.6996
	[−45, 45]	0.6773	0.6902	0.8554	0.8721	0.7198	0.6925
Fixed	[0, 90]	0.3814	0.3952	0.3806	0.371	0.4096	0.3968
	[−45, 45]	0.3891	0.3901	0.294	0.2865	0.4229	0.4078
	Nelem	–	200	144	144	36	25

Fig. 7

The average error and computational time of laminated plate elastic analysis.

The results have shown that the linear triangular and the 9-node quadrilateral elements are the best elements for laminated composite analysis from the accuracy point of view. However, the linear triangular element (LINTRI) needs higher number of elements, and subsequently, longer computational time. On the other hand, the worst elements for laminated composite analysis were the LINQUAD and the MKQ12 elements as they have the maximum relative average error.

Conclusions

In the present paper, five different thin shell finite elements were considered. The five elements were the Linear Triangular Element (LINTRI), the Linear Quadrilateral Element (LINQUAD), the Linear Quadrilateral Element Based on Deformation Modes (MKQ12), the 8-Node Quadrilateral Element (QUAD8NOD), and the 9-Node Quadrilateral Element (QUAD9NOD). A simple and detailed mathematical model to derive the interpolation functions and the stiffness matrix of each element was presented. The basis functions were selected from the well-known Pascal Triangle to minimize the order of the interpolation functions. However, singularities existed and specific terms had to be removed and replaced with other terms to eliminate the source of singularities. The five elements were tested using several elastic and aero-elastic analyses through three carefully selected bench mark problems with analytical or experimental results available in the literature. In order to have a fair comparison, a convergence analysis was conducted for each element and the minimum number of elements needed for convergence was used in the comparison.

From the present investigation, it was found that the MKQ12 element was the best choice for elastic free vibration analysis of a plate, since it yields the most accurate results with the minimum execution time. The second choice is the QUAD8NOD element, and the worst results are produced by the LINTRI element. In case of elastic thin shells, and if the stress analysis was sought, the most accurate elements are naturally the higher order elements. From the point of view of time and accuracy, the best element was found to be the QUAD8NOD element, and the QUAD9NOD element comes out second. The worst element for this kind of analysis was the LINTRI element, which requires longer computational times and produces discontinuous stress distributions.

For aero-elastic analysis, the most accurate results are obtained by using the LINTRI element, however, it requires the longest computational time. The best element for this type of analysis was found to be the QUAD9NOD element considering its accuracy and computational time. The second choice was the QUAD8NOD element. It is worth noting that in spite of its bad performance in elastic analysis, the LINTRI element was found to be more consistent to use with the Vortex Lattice Method in the aero-elastic analysis, but it requires long computation time. In laminated composite plate analysis, the best recommended elements are either the LINTRI or QUAD9NOD. However, the LINTRI element requires dense mesh, and subsequently longer computational time. The present results can serve many researchers and engineers interested in elastic and aero-elastic analyses, especially those who find difficulties in finding the detailed formulation of the finite elements, and those who are confused in selecting specific elements for specific application.

Conflict of Interest

The authors have declared no conflict of interest.

Compliance with Ethics requirements

This article does not contain any studies with human or animal subjects.