An 8-Node Plane Hybrid Element for Structural Mechanics Problems Based on the Hellinger-Reissner Variational Principle

1  Introduction

The FEM is a computational tool for engineering design analysis; it fundamentally revolutionized the way we perform scientific modeling and engineering design, including automobiles, aircraft, marine structures, bridges, highways, and high-rise buildings [1–9]. Traditional isoparametric elements are often limited by various numerical problems, such as limited computational accuracy and sensitivity to mesh distortion. For instance, the accuracy of the 4-node isoparametric element (Q4) may be obviously limited for bending problems due to shear locking. To improve the performance of these finite element methods, new finite element formulations need to be derived and applied to mechanical problems.

Over the past few decades, numerous techniques have been proposed for developing new finite element models based on various generalized variational principles. Pian et al. [10,11] pioneered the hybrid element method using an assumed stress field in the element and an assumed displacement field along element boundaries. Cen et al. [12] improved on Pian’s hybrid element approach to develop 8- and 12-node plane elements. The new elements perform well and can work in highly distorted meshes. Sze et al. [13–18] developed a hybrid stress-assumed natural strain 8-node solid shell element unaffected by shear, membrane, trapezoidal, thickness, and dilatational locking. This element is computationally more efficient than conventional hybrid elements because it adopts orthogonal assumed stress modes and enforces admissible sparsity in the flexibility matrix. Moreover, this element was generalized for smart structure modeling by the piezoelectric effect. Jog et al. [19] presented 4-node and 9-node elements for the development of hybrid axisymmetric elements based on the Hellinger-Reissner principle within the context of linear elasticity. Several examples were presented to show the excellent performance of both elements in various situations. Ma et al. [20] proposed a 24 degree of freedom (DOF) quadrilateral hybrid stress element for couple stress theory. Reduced integration and a stress smoothing technique were introduced to improve the accuracy of this element. Numerical results showed that this element could pass the enhanced patch test for a convergence condition. Bussamra et al. [21] proposed hexahedral hybrid-mixed finite elements for free vibration analysis of three-dimensional solids. These elements were free of shear and volumetric locking and had low sensitivity to mesh distortion. Agrawal et al. [22] used hybrid finite element methodology to develop an efficient and robust method based on the Hellinger-Reissner variational principle towards developing an efficient and robust finite element-based contact strategy. The proposed contact formulation can effectively model the contact interaction of thin as well as thick geometries as well as contact between bodies made of almost incompressible materials. Xie et al. [23] developed an atomistic-informed crystal plasticity finite element method to capture the evolution of geometrically compatible dislocation patterns as well as crystal plasticity in body-centered cubic crystals at micron or submicron scales. Liu et al. [24] established a novel three-dimensional progressive damage model based on the generalized mixed finite element method to investigate the strength and failure behaviors of notched composite laminate plates. The stress distribution of notched isotropic plates under tension was studied, and the results verified the good accuracy of the method. Ping and co-workers [25–27] used a specific finite element eigen-analysis method to predict the singular stress field around the tips of cracks, notches, and inclusions within the framework of the Hellinger-Reissner variational principle. Based on this, they [28] proposed an adaptive FEM to examine the issues with fracture propagation in anisotropic materials, bi-materials, and isotropic materials. Ramtekkar et al. [29] developed a refined 6-node two-dimensional mixed finite element to analyze a laminated composite using the minimum potential energy principle. The results obtained through the model showed excellent agreement with the elastic solution. Additionally, smoothed FEMs can also increase numerical solution accuracy, and the core idea behind this method is to combine the strain smoothing methodology of meshfree methods into the standard finite element method [30–32]. According to [33,34], smoothed FEMs can be viewed as specific assumed strain or stress methods for generalized variational principles.

The quadrilateral isoparametric elements contain Lagrangian elements and Serendipity elements. The Lagrangian elements are more computationally expensive than the Serendipity elements since they contain internal nodes. Thus, the Serendipity elements are preferred in practical applications, especially the 8-node isoparametric element (Q8). As the biquadratic interpolation functions are adopted for the Q8 element, it has better accuracy compared with the Q4 element and the Wilson-incompatible element. However, the Q8 element is extremely sensitive to mesh distortion [35]. The design of 8-node plane elements that can both improve the accuracy and overcome the above shortcomings has been a research hotspot. Dhananjaya et al. [36] presented a new 8-node Serendipity quadrilateral plate bending element based on Mindlin-Reissner theory using the integrated force method. This new element performs excellently in both thin and moderately thick plate bending situations. Dang et al. [37] presented a hybrid element model derived from a two-field variational functional; the plausible equilibrated stress field within each element was derived using Airy’s stress functions. The results indicated that this hybrid element was less sensitive to geometric distortion and was more accurate for stress than standard finite elements. Wang et al. [38,39] derived a new 8-node assumed stress quasi-conforming plane element based on the Bernstein basis function and Airy’s stress function. This element was not affected by mesh distortions and could adjust to the mesh shape, which deformed into a triangle or concave quadrangle with curved elements. However, derivations of 8-node plane hybrid elements for solving structural mechanics problems based on the Hellinger-Reissner variational principle have rarely been reported in previous literature. In this paper, the formulation of the PH-Q8-15β element is derived based on the Hellinger-Reissner variational principle, and the new element is tested and verified by numerical examples.

The outline of this paper is as follows. In Section 2, a new 8-node plane hybrid finite element is proposed based on the Hellinger-Reissner variational principle, and subsequently, the finite element formulation is derived. Numerical examples with results and discussions are presented in Section 3. Finally, a brief conclusion is presented in Section 4.

2  Finite Element Equations of the 8-Node Hybrid Element

2.1 The Hellinger-Reissner Variational Principle

Consider a linear elastic body with volume Ω in the Cartesian coordinate system. The domain boundary is denoted by ∂Ω=∂Ωu⋃∂Ωσ, where ∂Ωu and ∂Ωσ are boundaries over which the displacements and tractions are prescribed, respectively. When the elastic body is in static equilibrium subjected to external forces, the following relationships are satisfied:

The equilibrium equation

σij,j+P¯i=0inΩ,(1)

The strain–displacement relationship

εij=12(ui,j+uj,i)inΩ,(2)

The constitutive relationship

εij=∂ψ∂σij=SijklσklinΩ,(3)

The Neumann condition

σijnj=T¯ion∂Ωσ,(4)

The Dirichlet condition

ui=u¯ion∂Ωu.(5)

where σij and εij represent the components of the stress tensor and strain tensor, respectively; ui is the component of the displacement vector; P¯i and T¯i are the known load distributions in Ω and along the boundary ∂Ωσ, respectively; nj is the component of the outwards unit normal vector; Sijkl is the flexibility tensor; and the complementary energy density ψ is expressed as:

ψ=∫0σijεijdσij,(6)

Then, Eqs. (1)–(5) can be expressed in terms of two independent fields, ui and σij, by the complementary energy principle and the Lagrangian multiplier method:

ΠHR=∫Ω[−ψ(σij)+12σij(ui,j+uj,i)−P¯iui]dΩ−∫∂ΩσT¯iuidΓ−∫∂Ωuσijnj(ui−u¯i)dΓ,(7)

The domain of the elastic body is divided into n elements. Generally, the trial solution of displacement ui in Eq. (7) is compatible within the entire domain of the elastic body, and the stress trial solution σij is allowed to be piecewise continuous. We can define the energy functional from Eq. (7) as the sum of the individual elements,

ΠHR=∑nΠHRe=∑n{∫Ωe[−ψ(σij)+12σij(ui,j+uj,i)−P¯iui]dΩe−∫∂ΩσeT¯iuidΓ},(8)

where Ωe is the domain of an element e.

The stationary condition of Eq. (8) is:

δΠHR=0.(9)

2.2 Interpolation Functions of Variable Fields in Elements

A standard 8-node element in the natural coordinate system is shown in Fig. 1a, and a curved edge element after the isoparametric transformation in the Cartesian coordinate system is shown in Fig. 1b. For isoparametric elements, the interpolation functions of the coordinates and displacements are the same, i.e.,

X(ξ,η)=NXe,(10)

U(ξ,η)=Nqe,(11)

Figure 1: Plane quadrilateral 8-node isoparametric element

where X=(xe,ye)T is the coordinate matrix; Xe=(x1e,y1e,⋯,x8e,y8e)T is the node coordinate vector; U=(ue,ve)T is the displacement matrix; qe=(u1e,v1e,⋯,u8e,v8e)T is the node displacement vector; and the shape function matrix N can be expressed as:

N(ξ,η)=[N1I2,N2I2,N3I2,N4I2,N5I2,N6I2,N7I2,N8I2],(12)

where In is the n-th order identity matrix, and the shape functions are:

{Ni=14(1+ξiξ)(1+ηiη)(ξiξ+ηiη−1)(i=1,2,3,4),Ni=12(1−ξ2)(1+ηiη)(i=5,7),Ni=12(1−η2)(1+ξiξ)(i=6,8).(13)

where (ξi,ηi) is the coordinate values of the i point.

The strain matrix can be obtained from Eq. (11):

ε=(εxxεyyγxy)=(LN)qe=Bqe,(14)

where the geometric matrix B(x,y) can be expressed as:

B(x,y)=LN =[∂∂x0∂∂y0∂∂y∂∂x]T[N1I2,N2I2,N3I2,N4I2,N5I2,N6I2,N7I2,N8I2]=[B1,B2,B3,B4,B5,B6,B7,B8],(15)

and Bi is given by:

Bi=[∂Ni∂x0∂Ni∂y0∂Ni∂y∂Ni∂x]T(i=1,2,⋯,8).(16)

The relationship between the natural coordinates (ξ,η) and the Cartesian coordinates (x,y) is as follows:

(∂∂ξ∂∂η)=J(∂∂x∂∂y),(17)

where the Jacobi matrix J is given by:

J=[xξyξxηyη]=[∂x∂ξ∂y∂ξ∂x∂η∂y∂η].(18)

In the hybrid element analysis, the performance of the elements constructed on the basis of different assumed stress fields differs. Therefore, the key to deriving a hybrid element with superior performance lies in the selection of stress functions. Because the plane stresses are required to satisfy the continuity in each element, the stress functions are derived based on Eq. (1) by using the Airy stress function. To overcome the zero energy modes and ensure the efficiency of the calculation, it is necessary to determine the number of stress parameters and express them in nβ. For an 8-node plane hybrid element, its displacement DOF and rigid body DOF are expressed as nqe=16 and n0=3, respectively. We select 15 stress parameters in accordance with the theory presented in [40], nβ≥nqe−n0. The stress matrix is expressed as follows:

σ(ξ,η)=Pβ=PLβL+PHβH,(19)

where σ=(σξξ,σηη,σξη)T is the stress vector; P(ξ,η) and β are the stress function matrix and the stress parameter vector, respectively; βL=(β1,β2,β3)T and βH=(β4,⋯,β15)T are the low-order stress parameter vector and the high-order stress parameter vector, respectively; PL=I3 represents the third-order identity matrix; and the high-order stress matrix PH can be expressed as:

PH=[ξηξηη2000000000000ξηξηξ2000000000000ξηξ2η2],(20)

For a generic element, the stress function in the Cartesian coordinate system can be converted from Eq. (19) by the transformation matrix,

σ=Pβ=PLβL+TcPHβH,(21)

where σ=(σxx,σyy,σxy) T; the transformation matrix Tc is defined as:

Tc=[xξ2xη22xξxηyξ2yη22yξyηxξyξxηyηxξyη+xηyξ]ξ,η=0.(22)

Notably, the form of Eq. (20) obtained by using the Airy stress function is usually not unique. Eq. (20) possesses the advantages of good accuracy, good stability and strong universality in numerical examples. For special problems, such as fracture problems, more complex assumed stress functions can be used to construct the optimal element.

2.3 Derivation of the Finite Element Formulation

Substituting the displacement matrix U, the strain matrix ε and the stress matrix σ into Eq. (9), we can obtain:

ΠHR=∑eΠHRe=βTGqe−12βTHβ−(qe)TQe,(23)

where G=∫AePT(DN)tdAe, H=∫Ae(PTSP)tdAe, and Qe=∫Ae(NTP¯)tdAe+∫le(NTT¯)tdle; t is the element thickness; P¯ and T¯ are the body load vector and the surface load vector, respectively. For plane stress problems, the flexibility matrix S is:

S=(1E−μE0−μE1E000(1−μ)2).(24)

where E is the Young’s modulus and μ is the Poisson’s ratio.

From the stationary condition of the energy functional, δΠHRe(β,qe)=0, the following equations can be obtained:

Gqe−Hβ=0,(25)

GTβ−Qe=0,(26)

Eq. (25) imposes a strong geometric constraint on the displacement vector qe, whereas Eq. (26) imposes a weak equilibrium constraint on the stress parameter β. As a result, the hybrid element is stiffened in two directions at the same time. In general, more parameters in Eq. (19) imply a better approximation. However, when nβ≫nqe − n0, the element stiffens to the point where locking behaviors may occur, so 15 parameters are a reasonable choice. It can be seen that the 15 stress modes chosen in Eq. (19) are independent of one another, and Zhang et al. [41] demonstrated that linear independence of assumed stress modes is a necessary and sufficient condition for the flexibility matrix H to be non-singular. As a result of Eq. (25), it follows that:

β=H−1Gqe,(27)

Substituting Eqs. (27) into (26), we obtain:

(GTH−1G)qe=Keqe=Qe,(28)

where Ke is the stiffness matrix of the element.

By using Eqs. (23), (27) and (28), the energy functional of the elastic body can be expressed as:

ΠHR=∑eΠHRe=∑e[−12(qe)TKeqe+(qe)TQe]=−12qTKq+qTQ,(29)

where K is the stiffness matrix of the elastic body; q and Q are the displacement vector and the load vector of the node, respectively. From the variation of Eq. (29), we can obtain a set of algebraic equations:

Kq=Q.(30)

After imposing displacement boundary conditions, the stiffness matrix K becomes non-singular, allowing the Eq. (30) to be solved. In general, the stiffness matrix K generated from the Hellinger-Reissner variational principle can be decomposed into a traditional stiffness with a reduced integration method and a higher-order stiffness, with the latter improving hybrid element computational efficiency and stability [42,43].

3  Numerical Examples

In this section, seven problems are presented to evaluate the performance of the PH-Q8-15β element. Several elements were taken from the literature for comparison, which are listed in Table 1. In the following, the geometric and physical parameters of each model are dimensionless, and the relative error is defined as:

Err=|γ−γrefγref|.(31)

where γ and γref are the numerical solution and the reference solution, respectively.

3.1 Patch Test

A rectangular panel (0.24 × 0.001 × 0.12) is divided into five irregular elements, as shown in Fig. 2. The material properties of the panel are Young's modulus E=1500 and Poisson’s ratio μ=0.3. The low- and high-order displacement fields given by Eqs. (32) and (33) are applied to the boundary nodes. The numerical results for the internal nodes in Fig. 2 are consistent with the analytical solutions, which demonstrates that the new element passes the patch test.

{u(x,y)=x+y2,v(x,y)=y+x2,(32)

{u(x,y)=14+x+3y−2x2−4xy+52y2,v(x,y)=1+x2+2y−23x2+175xy+32y2.(33)

Figure 2: The patch test of membrane elements

3.2 Cantilever Beam with Five Irregular Elements

A cantilever beam with five irregular elements is subjected to the concentrated load F and bending moment M in Fig. 3. The geometry and material properties of the beam are as follows: length L=10, height h=2, width b=1, Young’s modulus E=1500, and Poisson’s ratio μ=0.3. Numerical results of the vertical deflection vA at point A and the stress σBx at point B are given in Tables 2 and 3. We find that compared with the traditional isoparametric elements, PS(a) element and QACM8 element, the PH-Q8-15β element has better accuracy.

Figure 3: A cantilever beam with five irregular elements

3.3 Sensitivity Test for Mesh Distortion and Convergence Analysis

A cantilever beam is subjected to two different loads, as shown in Fig. 4. The beam is discretized using two elements, and e is the distortion parameter. The geometry and material properties of the beam are the same as those in Subsection 3.2. The normalized results of the vertical displacement vA are listed in Tables 4 and 5. The tables show that the PH-Q8-15β element provides stable numerical solutions, even on distorted mesh shapes. Its stability is slightly less than that of the QC-Q8-B element.

Figure 4: A cantilever beam divided by two elements with distortion parameter e

The displacement norm error is calculated using a typical mesh type, as illustrated in Fig. 5. Numerical results are shown in Table 6. We can observe that the results for vertical displacement attain first order convergence.

Figure 5: Typical mesh type models for convergence test

3.4 Bending for a Curving Beam with the Uniform Cross Section

A curved cantilever beam subjected to a concentrated force F and the mesh shapes of the models are shown in Figs. 6 and 7, respectively. The geometric and physical parameters of the beam were selected from the literature for comparison [39]. For a thick curved beam, the geometry and material properties are as follows: inner radius Ri=10, middle radius Rm=12.5, thickness t=5, width b=1, Young’s modulus E=1000, Poisson’s ratio μ=0.3, and concentrated force F=600. For a thin curved beam, the geometry and material properties are as follows: the inner radius Ri=4.12, the middle radius Rm=4.22, the thickness t=0.2, the width b=0.1, the Young’s modulus E=106, the Poisson’s ratio μ=0.25, and the concentrated force F=1. Based on the Castigliano energy theorem, the analytical solution of the radial displacement uA for this problem is [46]:

uA=πFRm4(Rm2EI+561Gbt+1Ebt).(35)

where Rm, t and b represent the middle radius, thickness and width, respectively; E is Young’s modulus; the moment of inertia is I=bh312; and the shear modulus is G=E2(1+µ).

Figure 6: Curved cantilever beam

Figure 7: Typical mesh type models for the curved cantilever beam

Numerical results of the radial displacement uA are listed in Tables 7 and 8. From the tables, we find that the results are more accurate for the new element than the other elements. Moreover, as the number of PH-Q8-15β elements increases, the double-layer mesh can improve the numerical accuracy of the thick curved beam. Therefore, the thin curved beam should be numerically calculated by the single-layer mesh.

The effects of the mesh type on the numerical accuracy for the PH-Q8-15β element for different t/Ri ratios are investigated in Fig. 8. The vertical coordinate in the figures represents the relative error between the numerical results obtained by different meshes and the analytical solution for given t/Ri. The geometry and material properties are as follows: the inner radius Ri=10, the width b=0.1, the Young’s modulus E=1000, the Poisson’s ratio μ=0.3, and the concentrated force F=1. The figures show that when the value of t/Ri is close to 0.2, the single-layer mesh has better numerical accuracy, and it becomes more obvious with decreasing t/Ri. Therefore, when t/Ri≤0.2, the single-layer mesh should be preferred for a curved beam.

Figure 8: Influence of mesh type and t/Ri on the numerical accuracy of uA for the curved beam

3.5 Locking Test

This test, proposed by Macneal and Harder, is a well-known benchmark for testing sensitivity to mesh distortion [47]. Three different mesh shapes in Fig. 9 are adopted: rectangular, parallelogram and trapezoidal. Two load types are considered: the bending moment and the concentrated force. The geometry and material properties are as follows: length L=6, height h=0.2, width b=0.1, Young’s modulus E=1500, and Poisson’s ratio μ=0.3. The normalized results of the vertical deflection vA are listed in Table 9. Better and more accurate solutions can be obtained by the PH-Q8-15β element than by the other 8-node elements, and the PH-Q8-15β element can overcome shear and trapezoidal lockings.

Figure 9: MacNeal thin beams

Additionally, shear locking can be alleviated by using bubble functions for a general plane stress problem; the functions are typically high-order displacement polynomials that are zero on the element boundaries [48,49]. This strategy, however, may improve the accuracy of elements with linear displacement modes (e.g., the Q4 element) but is not appropriate for the Q8 element. Because the Q8 element has a fully displacement mode that simulates the constant bending strain state.

Membrane locking can be examined by the analysis of a ring, and many arch elements exhibit membrane locking in the analysis of deep arches even though they may give good results in the analysis of shallow arches [50]. A ring is subjected to two pinching concentrated loads F in Fig. 10a. Because of the symmetry, only a quarter of the ring is analyzed with appropriate boundary conditions in Fig. 10b. The geometry and material properties are as follows: inner radius Ri=5, thickness t=0.1, width b=0.1, Young’s modulus E=1×107, Poisson’s ratio μ=0.3, and concentrated force F=1.

Figure 10: Pinched ring

According to the conclusion in Subsection 3.4, the single-layer mesh is adopted for the quarter of the thin pinched ring. Numerical results of the radial deflection vA are listed in Table 10. The result obtained by eight PH-Q8-15β elements is consistent with the reference solution and converges to the reference solution quickly. Therefore, the PH-Q8-15β element can overcome membrane locking and can be successfully used to solve the pinched ring. We find that the numerical accuracy of the PH-Q8-15β element is lower than the HR-Q9 element and the MHR-Q9 element when the number of elements is 2.

3.6 Axisymmetric Annular Plate with Uniformly Distributed Internal Pressure

An axisymmetric annular plate is subjected to a uniformly distributed internal pressure q in Fig. 11a. The annular plate modelled by a quarter plate with appropriate boundary conditions is shown in Fig. 11b. The geometry and material properties are as follows: inner radius Ri=10, outer radius Ro=15, thickness t=0.1, Young’s modulus E=105, Poisson’s ratio μ=0.3, and internal pressure q=1. The polar coordinate system is established at the central of the circular hole, and the analytical solutions of this problem for the radial stress σρ and the hoop stress σθ can be expressed in Eq. (36):

σρ=−Ro2ρ2−1Ro2Ri2−1q,σφ=Ro2ρ2+1Ro2Ri2−1q.(36)

Figure 11: Axisymmetric annular plate with uniformly distributed internal pressure

The stress solutions of σρ and σθ are functions of the variable ρ, independent of the angle θ. Fig. 12 shows the effects of the polar radius ρ and the number of elements on the radial stress σρ and the hoop stress σθ, and the numerical results obtained by the PH-Q8-15β element are compared with the analytical solutions. With increasing ρ, the absolute values of σρ and σθ gradually decrease, and the numerical results gradually converge to the analytical solutions with the increase in the number of elements. Therefore, the PH-Q8-15β element can be successfully used for this problem.

Figure 12: The variation in the radial stress σρ and the hoop stress σθ relative to the polar radius ρ and number of elements for the axisymmetric annular plate

3.7 Straight/Curved Beams with Varying Cross Sections

A symmetric tapered beam is subjected to a concentrated load F in Fig. 13. The geometry and material properties are as follows: length L=10, height of clamped end h=1, height of cantilever end h=0.5, width b=1, Young’s modulus E=105, Poisson’s ratio μ=0.25, and concentrated force F=1. Table 11 reports the numerical results of the vertical deflection vA obtained by the PH-Q8-15β element and the NPB element. The “Axis Meshes” in Table 11 represent the number of evenly divided elements in the horizontal direction of the symmetric tapered beam. The table shows that the PH-Q8-15β element has better accuracy than the NPB element.

Figure 13: A symmetric tapered beam

A cantilever curving beam with a varying cross section is subjected to a concentrated load F in Fig. 14. The geometry and material properties are as follows: the middle radius Rm=5, the thickness of the clamped end t1=0.6(1 + a), the thickness of the cantilever end t2=0.6(1 − a), where a is the taper ratio, the width b=0.4, the Young’s modulus E=3×1010, the Poisson’s ratio μ=0.3, and the concentrated force F=1.

Figure 14: A curved beam with varying cross section

The single-layer mesh is adopted, as shown in Fig. 7. The results obtained by the 80 TNCB elements in [53] are taken as reference solutions and listed in Table 12. Fig. 15 shows the effects of the number of elements on the relative error between the numerical results obtained by the PH-Q8-15β element and the reference solution under different taper ratios a. The horizontal coordinate in figures represents the number of elements. We find that the PH-Q8-15β element can approximate the radial displacement uA and the tangential displacement vA with reasonable accuracy. Moreover, with the increase of the number of elements, the solution converges to the reference solution, which reflects the superior performance of the PH-Q8-15β element.

Figure 15: The variation in the relative error relative to the number of elements N for the curved beam with varying cross section

4  Conclusion

Because of the shortcomings of the Q8 element, a new PH-Q8-15β element is proposed for structural mechanics problems in this paper. We obtain the formulation of this element based on the Hellinger-Reissner variational principle. Fifteen unknown parameters are adopted in the selection of stress modes according to the design principle of Pian. Several numerical benchmark problems are presented to verify the performance of the PH-Q8-15β element. Numerical examples demonstrate that the PH-Q8-15β element has good accuracy and is not limited by various common locking behaviors of plane elements. Furthermore, the PH-Q8-15β element has better performance than the other elements in thin/thick curved beams with uniform cross sections. When the inner radius Ri of a curved beam is constant and t/Ri≤0.2, the curved beam should be modelled by a single-layer mesh. Finally, the PH-Q8-15β element can be successfully applied to the problems of straight/curved beams with varying cross sections, and numerical results obtained with only a few elements are able to accurately approximate the reference solution.

In the future, we will extend the new PH-Q8-15β element to intelligent structures considering the flexoelectric effect and micromachinery and nanomachinery involving nonlocal strain gradient theory.

Acknowledgement: The authors would like to thank the editor and the reviewers for their constructive comments which improve the quality of the paper.

Funding Statement: This paper was supported by the National Natural Science Foundation of China (No. 11572210).

Author Contributions: The authors confirm contribution to the paper as follows: methodology and writing original draft: Haonan Li; investigation and software: Wei Wang; investigation and data curation: Quan Shen; validation and project administration: Linquan Yao. All authors reviewed the results and approved the final version of the manuscript.

Availability of Data and Materials: Data and materials will be made available by reasonable requests.

Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.

References

1. Liu, W. K., Li, S., Park, H. S. (2022). Eighty years of the finite element method: Birth, evolution, and future. Archives of Computational Methods in Engineering, 29(6), 4431–4453. https://doi.org/10.1007/s11831-022-09740-9 [Google Scholar] [CrossRef]

2. Truong, V. H., Liu, J., Meng, X. H., Jiang, C., Nguyen, T. T. (2018). Uncertainty analysis on vehicle-bridge system with correlative interval variables based on multidimensional parallelepiped model. International Journal of Computational Methods, 15(5), 1850030. https://doi.org/10.1142/S0219876218500305 [Google Scholar] [CrossRef]

3. Li, D., Tian, J., Shi, S., Wang, S., Deng, J. et al. (2023). Lightweight design of commercial vehicle cab based on fatigue durability. Computer Modeling in Engineering & Sciences, 136(1), 421–445. https://doi.org/10.32604/cmes.2023.024133 [Google Scholar] [PubMed] [CrossRef]

4. Wang, Y., Li, Z., Li, Z., Wang, Y. (2023). Finite element simulation of radial tire building and shaping processes using an elasto-viscoplastic model. Computer Modeling in Engineering & Sciences, 135(2), 1187–1208. https://doi.org/10.32604/cmes.2022.022596 [Google Scholar] [PubMed] [CrossRef]

5. Mehraban, A., Tufo, H., Sture, S., Regueiro, R. (2021). Matrix-free higher-order finite element method for parallel simulation of compressible and nearly-incompressible linear elasticity on unstructured meshes. Computer Modeling in Engineering & Sciences, 129(3), 1283–1303. https://doi.org/10.32604/cmes.2021.017476 [Google Scholar] [PubMed] [CrossRef]

6. Liu, H., Zhang, Q., Ma, J. (2021). Thermo-mechanical dynamics of two-dimensional FG microbeam subjected to a moving harmonic load. Acta Astronautica, 178, 681–692. https://doi.org/10.1016/j.actaastro.2020.09.045 [Google Scholar] [CrossRef]

7. van Vinh, P., Huy, L. Q. (2022). Finite element analysis of functionally graded sandwich plates with porosity via a new hyperbolic shear deformation theory. Defence Technology, 18(3), 490–508. https://doi.org/10.1016/j.dt.2021.03.006 [Google Scholar] [CrossRef]

8. Bach, D. P., Brancherie, D., Cauvin, L. (2022). An embedded-FEM approach accounting for the size effect in nanocomposites. Computational Mechanics, 70(4), 745–762. https://doi.org/10.1007/s00466-022-02194-7 [Google Scholar] [CrossRef]

9. Choi, M. J., Klinkel, S., Sauer, R. A. (2022). An isogeometric finite element formulation for frictionless contact of Cosserat rods with unconstrained directors. Computational Mechanics, 70(6), 71107. https://doi.org/10.1007/s00466-022-02223-5 [Google Scholar] [CrossRef]

10. Pian, T. H. H. (1964). Derivation of element stiffness matrices by assumed stress distributions. AIAA Journal, 2(7), 1333–1336. https://doi.org/10.2514/3.2546 [Google Scholar] [CrossRef]

11. Pian, T. H. H., Sumihara, K. (1984). Rational approach for assumed stress finite elements. International Journal for Numerical Methods in Engineering, 20(9), 1685–1695. https://doi.org/10.1002/nme.1620200911 [Google Scholar] [CrossRef]

12. Cen, S., Fu, X. R., Zhou, M. J. (2011). 8- and 12-node plane hybrid stress-function elements immune to severely distorted mesh containing elements with concave shapes. Computer Methods in Applied Mechanics & Engineering, 200(29–32), 2321–2336. https://doi.org/10.1016/j.cma.2011.04.014 [Google Scholar] [CrossRef]

13. Sze, K. Y., Yao, L. Q. (2000). A hybrid stress ANS solid-shell element and its generalization for smart structure modelling. Part I—Solid-shell element formulation. International Journal for Numerical Methods in Engineering, 48, 545–564. [Google Scholar]

14. Sze, K. Y., Yao, L. Q., Yi, S. (2000). A hybrid stress ANS solid-shell element and its generalization for smart structure modelling. Part II—Smart structure modelling. International Journal for Numerical Methods in Engineering, 48, 565–582. [Google Scholar]

15. Sze, K. Y., Yao, L. Q. (2000). Modelling smart structures with segmented piezoelectric sensors and actuators. Journal of Sound and Vibration, 235(3), 495–520. https://doi.org/10.1006/jsvi.2000.2944 [Google Scholar] [CrossRef]

16. Sze, K. Y., Lo, S. H., Yao, L. Q. (2002). Hybrid-stress solid elements for shell structures based upon a modified variational functional. International Journal for Numerical Methods in Engineering, 53, 2617–2642. https://doi.org/10.1002/nme.402 [Google Scholar] [CrossRef]

17. Sze, K. Y., Yao, L. Q., Pian, T. H. H. (2002). An eighteen-node hybrid-stress solid-shell element for homogenous and laminated structures. Finite Elements in Analysis and Design, 38(4), 353–374. https://doi.org/10.1016/S0168-874X(01)00089-0 [Google Scholar] [CrossRef]

18. Yao, L. Q., Lu, L. (2003). Hybrid-stabilized solid-shell model of laminated composite piezoelectric structures under non-linear distribution of electric potential through thickness. International Journal for Numerical Methods in Engineering, 58, 1499–1522. https://doi.org/10.1002/nme.823 [Google Scholar] [CrossRef]

19. Jog, C. S., Annabattula, R. (2006). The development of hybrid axisymmetric elements based on the Hellinger-Reissner variational principle. International Journal for Numerical Methods in Engineering, 65, 2279–2291. https://doi.org/10.1002/nme.1552 [Google Scholar] [CrossRef]

20. Ma, X., Chen, W. (2013). 24-DOF quadrilateral hybrid stress element for couple stress theory. Computational Mechanics, 53(1), 159–172. https://doi.org/10.1007/s00466-013-0899-7 [Google Scholar] [CrossRef]

21. Bussamra, F. L. S., Lucena Neto, E., Cardoso, F. R. (2018). Three-dimensional hybrid-mixed stress elements for free vibration analysis. Finite Elements in Analysis and Design, 140(5–6), 50–58. https://doi.org/10.1016/j.finel.2017.11.004 [Google Scholar] [CrossRef]

22. Agrawal, M., Nandy, A., Jog, C. S. (2019). A hybrid finite element formulation for large-deformation contact mechanics. Computer Methods in Applied Mechanics & Engineering, 356(6), 407–434. https://doi.org/10.1016/j.cma.2019.07.017 [Google Scholar] [CrossRef]

23. Xie, Y., Li, S. (2021). Finite temperature atomistic-informed crystal plasticity finite element modeling of single crystal tantalum (α-Ta) at micron scale. International Journal for Numerical Methods in Engineering, 122, 4660–4697. https://doi.org/10.1002/nme.6741 [Google Scholar] [CrossRef]

24. Liu, W., Yu, F., He, Z., Qing, G. (2020). A high-precision progressive damage model based on generalized mixed finite element method. Archive of Applied Mechanics, 90(3), 559–571. https://doi.org/10.1007/s00419-019-01625-x [Google Scholar] [CrossRef]

25. Chen, M. C., Ping, X. C., Xie, H. M., Liu, Z. W. (2008). Numerical and experimental analyses of singular electro-elastic fields around a V-shaped notch tip in piezoelectric materials. Engineering Fracture Mechanics, 75(18), 5029–5041. https://doi.org/10.1016/j.engfracmech.2008.05.011 [Google Scholar] [CrossRef]

26. Ping, X. C., Chen, M. C., Zheng, B. B., Xu, B. (2013). An effective numerical analysis of singular stress fields in dissimilar material wedges under thermo-mechanical loads. Engineering Fracture Mechanics, 106(17), 22–37. https://doi.org/10.1016/j.engfracmech.2013.03.028 [Google Scholar] [CrossRef]

27. Ping, X. C., Wang, C. G., Cheng, L. P., Chen, M. C., Xu, J. Q. (2018). A super crack front element for three-dimensional fracture mechanics analysis. Engineering Fracture Mechanics, 196(4), 1–27. https://doi.org/10.1016/j.engfracmech.2018.04.016 [Google Scholar] [CrossRef]

28. Wang, C. M., Ping, X. C., Wang, X. X. (2023). An adaptive finite element method for crack propagation based on a multifunctional super singular element. International Journal of Mechanical Sciences, 247(12), 108191. https://doi.org/10.1016/j.ijmecsci.2023.108191 [Google Scholar] [CrossRef]

29. Ramtekkar, G. D., Patel, K. S., Desai, Y. M. (2023). A six-node refined mixed finite element model for the analysis of fiber reinforced polymer composite beams. Composite Structures, 304(2), 116418. https://doi.org/10.1016/j.compstruct.2022.116418 [Google Scholar] [CrossRef]

30. Liu, G. R., Nguyen-Thoi, T. (2010). Smoothed finite element methods. New York: CRC Press. [Google Scholar]

31. Nguyen, T. K., Nguyen, V. H., Thanh, C. D., Vo, T. P., Nguyen-Xuan, H. (2016). Static and vibration analysis of isotropic and functionally graded sandwich plates using an edge-based MITC3 finite elements. Composites Part B: Engineering, 107(8), 162–173. https://doi.org/10.1016/j.compositesb.2016.09.058 [Google Scholar] [CrossRef]

32. Lee, C., Kim, S., Lee, P. S. (2021). The strain-smoothed 4-node quadrilateral finite element. Computer Methods in Applied Mechanics & Engineering, 373, 113481. https://doi.org/10.1016/j.cma.2020.113481 [Google Scholar] [CrossRef]

33. Liu, G. R., Nguyen-Xuan, H., Nguyen-Thoi, T. (2010). A theoretical study on the smoothed FEM (S-FEM) models: Properties, accuracy and convergence rates. International Journal for Numerical Methods in Engineering, 84(10), 1222–1256. https://doi.org/10.1002/nme.2941 [Google Scholar] [CrossRef]

34. Lee, C., Park, J. (2021). A variational framework for the strain-smoothed element method. Computers & Mathematics with Applications, 94, 76–93. https://doi.org/10.1016/j.camwa.2021.04.025 [Google Scholar] [CrossRef]

35. Lee, N. S., Bathe, K. J. (1993). Effects of element distortions on the performance of isoparametric elements. International Journal for Numerical Methods in Engineering, 36, 3553–3576. https://doi.org/10.1002/nme.1620362009 [Google Scholar] [CrossRef]

36. Dhananjaya, H. R., Pandey, P. C., Nagabhushanam, J. (2009). New eight node serendipity quadrilateral plate bending element for thin and moderately thick plates using integrated force method. Structural Engineering and Mechanics, 33(4), 485–502. https://doi.org/10.12989/sem.2009.33.4.485 [Google Scholar] [CrossRef]

37. Dang, T. D., Hung, N. D. (2013). A hybrid element model for structural mechanics problems. European Journal of Mechanics-A/Solids, 42(3), 469–479. https://doi.org/10.1016/j.euromechsol.2013.08.001 [Google Scholar] [CrossRef]

38. Wang, C., Wang, Y., Yang, C., Zhang, X., Hu, P. (2017). 8-node and 12-node plane elements based on assumed stress quasi-conforming method immune to distorted mesh. Engineering Computations, 34(8), 2731–2751. https://doi.org/10.1108/EC-11-2016-0404 [Google Scholar] [CrossRef]

39. Wang, C., Lu, X., Zhang, X., Hu, P. (2018). 8-node quasi-conforming plane element by using Bernstein basis functions. European Journal of Mechanics-A/Solids, 70(1–5), 127–140. https://doi.org/10.1016/j.euromechsol.2018.02.003 [Google Scholar] [CrossRef]

40. Pian, T. H. H., Wu, C. C. (2006). Hybrid and incompatible finite element methods. New York: Chapman & Hall/CRC. [Google Scholar]

41. Zhang, C. H., Feng, W., Huang, Q. (2002). The stress subspace of hybrid stress element and the diagonalization method for flexibility matrix H. Applied Mathematics and Mechanics, 23(11), 1263–1273. https://doi.org/10.1007/BF02439457 [Google Scholar] [CrossRef]

42. Sze, K. Y. (2003). A novel approach for devising higher-order hybrid elements. International Journal for Numerical Methods in Engineering, 36(19), 3303–3316. https://doi.org/10.1002/nme.1620361907 [Google Scholar] [CrossRef]

43. Sze, K. Y., Yang, X. M., Yao, L. Q. (2004). Stabilized plane and axisymmetric piezoelectric finite element models. Finite Elements in Analysis and Design, 40(9–10), 1105–1122. https://doi.org/10.1016/j.finel.2003.06.002 [Google Scholar] [CrossRef]

44. Cen, S., Chen, X. M., Fu, X. R. (2007). Quadrilateral membrane element family formulated by the quadrilateral area coordinate method. Computer Methods in Applied Mechanics & Engineering, 196(41–44), 4337–4353. https://doi.org/10.1016/j.cma.2007.05.004 [Google Scholar] [CrossRef]

45. Fu, X. R., Cen, S., Li, C. F. (2010). Analytical trial function method for development of new 8-node plane element based on the variational principle containing Airy stress function. Engineering Computations, 27(4), 442–463. https://doi.org/10.1108/02644401011044568 [Google Scholar] [CrossRef]

46. Tayşi, N., Göĝüş, M. T., Özakça, M. (2011). Structural analysis of arches in plane with a family of simple and accurate curved beam elements based on Mindlin-Reissner model. Journal of Mechanics, 27(1), 129–138. https://doi.org/10.1017/jmech.2011.14 [Google Scholar] [CrossRef]

47. Macneal, R. H., Harder, R. L. (1985). A proposed standard set of problems to test finite element accuracy. Finite Elements in Analysis and Design, 1(1), 3–20. https://doi.org/10.1016/0168-874X(85)90003-4 [Google Scholar] [CrossRef]

48. Kemp, B. L., Cho, C. M., Lee, S. W. (1998). A four-node solid shell element formulation with assumed strain. International Journal for Numerical Methods in Engineering, 43(5), 909–924. https://doi.org/10.1002/(SICI)1097-0207(19981115)43:5<909::AID-NME450>3.0.CO;2-X [Google Scholar] [CrossRef]

49. Nassehi, V., Parvazinia, M. (2009). A multiscale finite element space-time discretization method for transient transport phenomena using bubble functions. Finite Elements in Analysis and Design, 45(5), 315–323. https://doi.org/10.1016/j.finel.2008.10.005 [Google Scholar] [CrossRef]

50. Shi, G., Voyiadjis, G. Z. (1991). Simple and efficient shear flexible two-node arch/beam and four-node cylindrical shell/plate finite elements. International Journal for Numerical Methods in Engineering, 31, 759–776. https://doi.org/10.1002/nme.1620310408 [Google Scholar] [CrossRef]

51. Sze, K. Y., Fan, H., Chow, C. L. (2010). Elimination of spurious pressure and kinematic modes in biquadratic nine-node plane element. International Journal for Numerical Methods in Engineering, 38(23), 3911–3932. https://doi.org/10.1002/nme.1620382302 [Google Scholar] [CrossRef]

52. Auricchio, F., Balduzzi, G., Lovadina, C. (2015). The dimensional reduction approach for 2D non-prismatic beam modelling: A solution based on Hellinger-Reissner principle. International Journal of Solids and Structures, 63(12), 264–276. https://doi.org/10.1016/j.ijsolstr.2015.03.004 [Google Scholar] [CrossRef]

53. Tufekci, E., Eroglu, U., Aya, S. A. (2017). A new two-noded curved beam finite element formulation based on exact solution. Engineering with Computers, 33(2), 261–273. https://doi.org/10.1007/s00366-016-0470-1 [Google Scholar] [CrossRef]