Sage Journals: Discover world-class research

Abstract

Two novel plate-bending elements are developed and investigated in this study. Elements with 13 and 15 degree-of-freedoms are named AT13 and AT15, respectively. These triangular elements are formulated in a semi-analytic way. For this aim, the basic elasticity function is employed with unknown parameters. Subsequently, the trial-and-error procedure is used to determine the unidentified constants. Besides, the achieved results are compared with those obtained by displacement-based triangular elements with the same degrees-of-freedom (TUBA13 and TUBA15). In this research, both stress and displacement responses of diverse structures are assessed. After performing extensive numerical studies, the findings clearly demonstrate the superiorities of the proposed elements.

Keywords

Triangular element complementary energy method plate problem finite element method variational techniques

Introduction

Conventional elasticity theories employ strain as the only deformation measure. So, these theories do not take the internal length scale of the material into account. Conversely, numerical methods are usually used to solve the edge effect problem, and among them, the FEM is generally regarded as the most efficient and popular one.^1–4 Nevertheless, due to the inherent theoretical deficiencies, only poor performances can be expected when the usual finite elements are used.⁵ Among these, Remešíková et al.⁶ enhanced the performance of the triangular elements using the discrete Lagrangian schemes. The obtained results were also compared with those obtained experimentally. Results demonstrated that the performance of triangular elements was extensively improved because of using the grasshopper component for Rhinoceros. In addition to Lagrangian schemes, the zigzag theory was also developed to assess the structural performance of elements accurately.^7–9 Recently, Yekkalam Tash and Navayi Neya¹⁰ developed a new element to assess the performance of thick bending plates. They showed the accuracy of the proposed element, and plates with different thicknesses were evaluated. For this aim, the displacement of potential techniques was employed. They found that the ultimate deformation of the plates with different thicknesses was moved from the center to the thinner edges. In another investigation, Magnucka-Blandzi et al.¹¹ evaluated the flexural behavior and buckling potential of a circular plate. By using two elements, the three-layers plate was measured. The achieved results illustrated the accuracy of the proposed elements to assess the performance of bending plates. In 2018, Zhang et al.¹² utilized the finite difference method to improve meshes’ performance in order to measure the behavior of circular plates. Elements were established using energy formulations. Based on the obtained results, using energy schemes led to improving the performance of bending plate elements considerably.¹³

In 2013, Li et al.¹⁴ utilized a wavelet way to modify plate behavior prediction. During the integration procedure, multiple boundary circumstances could be applied forthrightly. The obtained outcomes demonstrated the positive role of wavelet theories in order to improve the accuracy of bending plate elements. Levyakov and Kuznetsov¹⁵ presented a new triangular element in order to solve composite plates. The obtained formulas permit one to freely measure the strain energy of the element. Hence, the stress and strain components of the element edges could be measured accurately. Numerical consequences were given to illustrate the nonlinear capacity of the proposed element. Han et al.¹⁶ used wavelet theories to assess the behavior of bending plates. For this purpose, a rectangular element was used. This element was established by the combination of both wavelet and Monte Carlo schemes. Their conclusion showed that using wavelet theories played an effective role in obtaining high numerical accuracy and converging fast. In 2010, Choo et al.¹⁷ studied the performance of the triangular and rectangular plate elements. To develop new elements, the hybrid-Trefftz model was utilized. Their achieved outcomes demonstrated that their proposed elements were robust, accurate, and free of shear locking in the thin plates.

Dey et al.¹⁸ developed a new element in order to evaluate the behavior of composite shells under vibration. Their element could predict the displacement of a member with low degree-of-freedoms. To show the accuracy of the proposed element, different numerical examples were utilized. According to the conclusion of their study, free vibration examination of plates using the presented element could reach highly accurate outcomes and reduce the calculation cost. In 2012, Dehghan and Sabouri¹⁹ used triangular and quadrilateral elements in order to solve the Pennes bioheat transfer equation. For the triangular elements, the error was achieved when quadrature points coincided with the nodal points. This scheme was also utilized to solve the equations in order to calculate the temperature influence on the thin bending plates. Zhang et al.²⁰ evaluated the performance of triangular elements under uniform load in a 3D piezoelectric medium. The proposed models were presented simply by a linear combination of three kinds of elementary functions, linear, trigonometric and logarithm ones. Their elements were useful to assess the behavior of the bending plate with highly accurate values. Recently, Geng et al.²¹ used a B-spline wavelet to improve the accuracy of finite element schemes in order to anticipate the dynamic behavior of thin plates. The numerical consequences demonstrated that both the computing efficiency and stability of the proposed schemes were higher than those of the conventional finite element methods. In 2018, Shirmohammadi and Bahrami²² assessed the dynamic performance of circular plates using a spectral element scheme. A novel formulation was projected in the frequency field for conducting a spectral element matrix. The obtained results showed the high accuracy of their proposed formulation for the members under dynamic and impact loads. To remedy this defect, the Airy stress function was used in this study.

For researchers, the FEM is a frequently utilized numerical approach in existing designs and simulations.²³ Nonetheless, a strong process is still plagued by several numerical issues that have yet to be resolved. One of these is the mesh distortion sensitivity problem.²⁴ The convergence of stress solutions is frequently not as good as that of displacement solutions due to the inherent theoretical flaws of displacement-based element models. The discordant displacement modes,^25–27 the improved strain method,²⁸ the selectively reduced integration system,²⁹ the new spline FEM,³⁰ the unsymmetric interpolation technique,³¹ and natural coordinate ways^32–36 have all been proposed in the last 50 years for developing strong finite element models that are insensitive to mesh distortion. Although the preceding methods can help to improve the FEM’s robustness, it should be noted that most models will fail once the element form is sufficiently distorted. In the case of a rounded quadrangle degenerating into a triangle or a concave quadrangle, for example. Many scholars investigated alternative meshless approaches^34–42 to overcome this challenge, which has substantially greater processing costs.

This paper is dedicated to the behavior study of four novel plate-bending elements. All of them have triangular shapes, but two of these are stress-based, and the others are displacement-based. Besides presenting element formulations, a lot of problems are solved to show the merit of the new elements. In order to find a semi-analytic solution, the errors between the exact displacement-based finite element and recommended scheme are also minimized. All the obtained results show the superior of the presented new elements.

Elements formulation

The deformed shape of a plate is defined using the transverse displacement at the midplane (w). Having this field variable, the rotation and the corresponding displacements at any point could be obtained using the following relationship, as demonstrated in Figure 1. Where, $θ_{x}$ and $θ_{y}$ indicate the rotation angle between the transverse displacement and the direction of the y- and x-axis, respectively.

θ = {\begin{matrix} θ_{x} \\ θ_{y} \end{matrix}} = - [\begin{matrix} \frac{\partial}{\partial x} \\ \frac{\partial}{\partial y} \end{matrix}] w = - D_{1} w

(1)

u = {\begin{matrix} u_{x} \\ u_{y} \\ u_{z} \end{matrix}} = {\begin{matrix} z θ_{x} \\ z θ_{y} \\ w \end{matrix}}

(2)

In which, $u_{x}$ , $u_{y}$ , and $u_{z}$ indicate the displacement components in the x, y, and z directions. Besides, the rotational components could be achieved by $θ_{x}$ and $θ_{y}$ , relative to that of y- and x-direction, respectively. The curvatures are introduced from the definition of the corresponding non-zero strains:

\begin{matrix} ε = {\begin{matrix} ε_{xx} \\ ε_{yy} \\ 2 ε_{xy} \end{matrix}} = [\begin{matrix} \frac{\partial}{\partial x} \\ 0 \\ \frac{\partial}{\partial y} \end{matrix} \begin{matrix} 0 \\ \frac{\partial}{\partial y} \\ \frac{\partial}{\partial x} \end{matrix}] {\begin{matrix} u_{x} \\ u_{y} \end{matrix}} = z D_{2} {\begin{matrix} θ_{x} \\ θ_{y} \end{matrix}} \\ = - z [\begin{matrix} \frac{\partial^{2}}{\partial x^{2}} \\ \frac{\partial^{2}}{\partial y^{2}} \\ 2 \frac{\partial^{2}}{\partial x \partial y} \end{matrix}] w = z ℵ \end{matrix}

(3)

ℵ = {\begin{matrix} ℵ_{xx} \\ ℵ_{yy} \\ 2 ℵ_{xy} \end{matrix}} = - [\begin{matrix} \frac{\partial^{2}}{\partial x^{2}} \\ \frac{\partial^{2}}{\partial y^{2}} \\ 2 \frac{\partial^{2}}{\partial x \partial y} \end{matrix}] w = - D_{2} D_{1} w = - D w

(4)

Figure 1.

Displacement and rotation components of an infinitesimal plate element.⁴³

Furthermore, the moments at any point are defined as the integral through the thickness of the plate, h, of the stresses times the distance to the midplane. Besides, the shear forces correspond to the integrals of the out-of-plane shear stresses and could be calculated as follows.⁴³

n = {\begin{matrix} m_{xx} \\ m_{yy} \\ m_{xy} \end{matrix}} = \int_{- h / 2}^{h / 2} z {\begin{matrix} σ_{xx} \\ σ_{yy} \\ σ_{xy} \end{matrix}} dz

(5)

q = {\begin{matrix} q_{x} \\ q_{y} \end{matrix}} = \int_{- h / 2}^{h / 2} {\begin{matrix} σ_{xz} \\ σ_{yz} \end{matrix}} dz

(6)

Where, $m_{ij}$ and $q_{i}$ are the moment and shear forces. If $\bar{p}$ is a distributed transverse load, the equilibrium of an infinitesimal element of the plate could be written in the following form:

\frac{\partial q_{x}}{\partial x} + \frac{\partial q_{y}}{\partial y} + \bar{p} = 0

(7)

\frac{\partial m_{xx}}{\partial x} + \frac{\partial m_{xy}}{\partial y} = q_{x}

(8)

\frac{\partial m_{xy}}{\partial x} + \frac{\partial m_{yy}}{\partial y} = q_{y}

(9)

D_{1}^{T} q + \bar{p} = 0

(10)

D_{2}^{T} m = q

(11)

These formulas can be combined as:

\frac{\partial^{2} m_{xx}}{\partial x^{2}} + \frac{\partial^{2} m_{yy}}{\partial y^{2}} + 2 \frac{\partial^{2} m_{yy}}{\partial x \partial y} + \bar{p} = 0

(12)

D^{T} m + \bar{p} = 0

(13)

The plane stress linear elastic constitutive law leads to the following relations between moments and curvatures:

m = \frac{E h^{3}}{12 (1 - v^{2})} [\begin{matrix} 1 & v & 0 \\ v & 0 & 0 \\ 0 & 0 & \frac{1 - v}{2} \end{matrix}] ℵ

(14)

ℵ = \frac{12}{E h^{3}} [\begin{matrix} 1 & - v & 0 \\ - v & 1 & 0 \\ 0 & 0 & 2 (1 + v) \end{matrix}] m

(15)

$E$ , h, and v are the elastic moduli, thickness and poison ratio, respectively. These relations have the compact notation as $m = k ℵ$ and $ℵ = fm$ , with $f = k^{- 1}$ . For simplicity, initial or thermal strains are not considered. Compared to the classical elasticity, micro-rotation $ϕ_{i}$ is introduced as a kinematic variable in addition to displacement $u_{i}$ in Cosserat elasticity. The strain $ε_{ij}$ depends on both displacement and the micro-rotation. Hence, the latter gives rise to the curvature $k_{ij}$ , can be written as follows⁴⁴:

ε_{ij} = u_{i, j} + ϵ_{ijk} ϕ_{k}

(16)

k_{ij} = ϕ_{j, i}

(17)

Where, $ϵ_{ijk}$ is the permutation symbol. It can be easily noticed that the above equations are reduced to the classical counterpart when the micro-rotation is neglected. The constitutive relations between the stress and the strain, the couple-stress, and the curvature are defined as:

σ_{ij} = C_{ijkl} ε_{kl}

(18)

m_{ij} = {\hat{C}}_{ijkl} k_{kl}

(19)

In which:

C_{ijkl} = μ [(1 + a) g_{ik} g_{jl} + (1 - a) g_{il} g_{jk} + \frac{λ}{μ} g_{ij} g_{kl}]

(20)

{\hat{C}}_{ijkl} = 4 μ l^{2} (g_{ik} g_{jl} + b g_{il} g_{jk} + c g_{ij} g_{kl})

(21)

$g^{ij}$ is the metric tensor, λ and $μ$ are Lamé constants, $l$ is the characteristic length of the material, while a, b, and c are the micropolar material constants. The equilibrium equations take the forms:

σ_{ij, i} + {\bar{p}}_{j} = 0, m_{ij, i} + ϵ_{ijk} σ_{ik} + {\bar{q}}_{j} = 0

(22)

Here, the ${\bar{p}}_{j}$ and ${\bar{q}}_{j}$ are the prescribed body force and body couple per unit volume. Without losing generality, the following forms of essential and natural boundary conditions are considered:

u_{i} + {\bar{u}}_{j} = 0, ϕ_{i} - {\bar{ϕ}}_{j} = 0 on Γ_{D}

(23)

σ_{ij} n_{i} - {\bar{σ}}_{j} = 0, m_{ij} n_{i} - {\bar{m}}_{j} = 0 on Γ_{N}

(24)

Γ_{D} \cap Γ_{N} = null and Γ_{D} \cup Γ_{N} = \partial Ω

(25)

$Ω$ is the problem domain and $\partial Ω$ is the boundary of $Ω$ . Also, $n_{i}$ indicates the unit outward normal vector component of $\partial Ω$ and $Γ_{D}$ and $Γ_{N}$ are the portion of $\partial Ω$ prescribed with the essential and natural boundary conditions, respectively. Therefore, the total potential energy functional for a Cosserat continuum is expressed as

\begin{matrix} Π = \int_{Ω} (\frac{1}{2} C^{ijkl} ε_{ij} ε_{kl} + \frac{1}{2} \hat{C}^{ijkl} k_{il} k_{kl} - {\bar{p}}_{j} u_{j} - {\bar{q}}_{j} ϕ_{j}) d Ω \\ - \int_{Γ_{N}} ({\bar{σ}}_{j} u_{j} + {\bar{m}}_{j} ϕ_{j}) d Γ \end{matrix}

(26)

It should be stated that the boundary conditions have to be satisfied with equations (7)–(9). Moreover, the stress and strain expressions can be simplified by introducing the symmetric “S” and anti-symmetric “A” stress and strain tensors, namely:

σ_{ij} = σ_{ij}^{S} + σ_{ij}^{A} and ε_{ij} = ε_{ij}^{S} + ε_{ij}^{A}

(27)

Where

σ_{ij}^{S} = \frac{1}{2} (σ_{ij} + σ_{ji})

(28)

σ_{ij}^{A} = \frac{1}{2} (σ_{ij} - σ_{ji})

(29)

ε_{ij}^{S} = \frac{1}{2} (ε_{ij} + ε_{ji}) = \frac{1}{2} (u_{i, j} + u_{j, i})

(30)

ε_{ij}^{A} = \frac{1}{2} u_{j, i} - u_{i, j} - ϵ_{ijk} ϕ_{k}

(31)

Then, the constitutive relations could be rewritten as:

σ = C_{s} ε, τ = C_{A} γ, m = C_{m} k

(32)

In which

ε = [ε_{11} ε_{22} ε_{33} 2 ε_{23}^{S} 2 ε_{31}^{S} 2 ε_{12}^{S}]

(33)

σ = [σ_{11} σ_{22} σ_{33} σ_{23}^{S} σ_{31}^{S} σ_{12}^{S}]

(34)

γ = {[2 ε_{23}^{A} 2 ε_{31}^{A} 2 ε_{12}^{A}]}^{T}

(35)

τ = {[σ_{23}^{A} σ_{31}^{A} σ_{12}^{A}]}^{T}

(36)

$ε$ , $σ$ , $γ$ , and $τ$ are vectors of symmetric strain, symmetric stress, anti-symmetric strain, and anti-symmetric stress components, respectively. Moreover, m and k are the vectors of couple-stress and curvature components, correspondingly. $C_{s}$ , $C_{A}$ , and $C_{m}$ are the corresponding material constitutive matrices. By virtue of the above-induced notations, the function in equation (26) for finite element formulation can be further expressed as:

Π_{p} = \sum_{e} Π_{p}^{e} and Π_{p}^{e} = \frac{1}{2} \int_{Ω^{e}} (ε^{T} C_{S} ε + γ^{T} C_{A} γ + k^{T} C_{m} k) d Ω - f^{e}

(37)

Where

f^{e} = \int_{Ω^{e}} ({\bar{p}}_{j} u_{j} + {\bar{q}}_{j} ϕ_{j}) d Ω + \int_{Γ_{N}^{e}} ({\bar{σ}}_{j} u_{j} + {\bar{m}}_{j} ϕ_{j}) d Γ

(38)

$f^{e}$ is the elemental load potential, $Ω^{e}$ denotes the element domain and $Γ_{N}^{e} = \partial Ω^{e} \cap Γ_{N}$ . By relaxing the first two constitutive relations in the potential energy functional, the following partial Hellinger-Reissner functional can be written for Cosserat continua:

Π_{HR} = \sum_{e} Π_{HR}^{e}

(39)

\begin{matrix} Π_{HR}^{e} = \\ \int_{Ω^{e}} (- \frac{1}{2} σ^{T} C_{S}^{- 1} σ + σ^{T} ε - \frac{1}{2} τ^{T} C_{A}^{- 1} τ + τ^{T} γ + \frac{1}{2} k^{T} C_{m} k) d Ω - f^{e} \end{matrix}

(40)

In this function, $σ$ and $τ$ need neither be continuous across the element boundary nor satisfy the equilibrium condition. Within the element, the stress and strain fields can be expressed symbolically as:

σ = P β, ε = B_{S} d, τ = Q α + \bar{τ}, γ = B_{A} d, k = B_{m} d

(41)

Where, P and Q are the shape function matrices for symmetric stress, r, and anti-symmetric stress, $τ$ , respectively. Also, $β$ and $α$ are the vectors of coefficients and d is the vector comprising the nodal displacements and micro-rotations. $B_{S}$ , $B_{A}$ , and $B_{m}$ are the relationship matrices between the symmetric strain, anti-symmetric strain, curvature, and the elementwise nodal displacement and micro-rotation vector, respectively. Moreover, $\bar{τ}$ is a prescribed quantity to be discussed later. Substitution of equation (26) into the second expression of equation (25) yields:

\begin{matrix} Π_{HR}^{e} = - \frac{1}{2} β^{T} H_{S} β + β^{T} G_{S} d - \frac{1}{2} α^{T} H_{A} α \\ + α^{T} G_{S} d - α^{T} M + Ld + \frac{1}{2} d^{T} K_{m} d - f^{e} + c_{1} \end{matrix}

(42)

In which

\begin{matrix} {H_{S}, H_{A}, G_{S}, G_{A}, M, L, K_{m}} = \\ \int_{Ω^{e}} {P^{T} C_{S}^{- 1} P, Q^{T} C_{A}^{- 1} Q, P^{T} B_{S}, Q^{T} B_{A}, Q^{T} C_{A}^{- 1} \bar{τ}, {\bar{τ}}^{T} B_{A}, B_{m}^{T} C_{m} B_{m}} d Ω \end{matrix}

(43)

Here, $c_{1}$ is a constant scalar arising from $\bar{τ}$ . Variation of the functional $w . r . t . β$ and $α$ lead to the next relations:

H_{S} β = G_{S} d

(44)

H_{A} α = G_{A} d - M

(45)

Finally, the following equation could be concluded:

\begin{matrix} Π_{HR}^{e} = \frac{1}{2} d^{T} (G_{S}^{T} H_{S}^{- 1} G_{S} + G_{A}^{T} H_{A}^{- 1} G_{A} + K_{m}) d \\ - d^{T} (G_{A}^{T} H_{A}^{- 1} M - L^{T}) - f^{e} + c_{1} \end{matrix}

(46)

In this study, two well-known elements with 13 and 15 nodes, named TUBA13 and TUBA15, are utilized. These elements are formulated by using the compulsory energy scheme. Besides, the newly presented stress-based elements are named AT13 and AT15. Figure 2 shows the proposed elements. To increase the degree of approximation and look for the adequate distribution of DOFs in the element, both $w$ and $w_{n}$ along each side should be uniquely defined by the DOFs presented on that side. It should be noted that the purpose of nodes 13, 14, 15 in TA15 is to increase the accuracy of the element to find responses inside the element. Also, shear locking can occur in first-order, fully-integrated elements (C3D8) that are subjected to bending. This occurs when artificial shear strain develops due to an inability of the element edges to bend. So, to solve this issue, in TUBA and TA elements (TA elements were established based on well-known TUBA displacement-based elements), the mid-side nodes will follow the bending curve, avoiding the shear strain in the horizontal plane and shear locking, as also mentioned by Argyris et al.⁴⁵

Figure 2.

Proposed triangular elements TA: (a) TA13 and (b) TA15.

In the finite element scheme, the complementary energy function has the next form:

Π_{C} = Π_{C}^{*} + V_{C}^{*} = \frac{1}{2} \int \int_{A^{e}}^{1} σ^{T} C σ tdA - \int_{Γ}^{1} T^{T} \bar{U} tds

(47)

In which:

Π_{C}^{*} = \frac{1}{2} \int \int_{A^{e}}^{1} σ^{T} C σ tdA

(48)

V_{C}^{*} = - \int_{Γ}^{1} T^{T} \bar{U} tds

(49)

\begin{matrix} σ = {\begin{matrix} σ_{11} \\ σ_{22} \\ σ_{12}^{S} \end{matrix}} = [\begin{matrix} λ + 2 μ & λ & 0 \\ λ & λ + 2 μ & 0 \\ 0 & 0 & μ \end{matrix}] \\ {\begin{matrix} ε_{11} \\ ε_{22} \\ 2 ε_{12}^{S} \end{matrix}}, σ_{12}^{A} = μ a, 2 ε_{12}^{A} \end{matrix}

(50)

T = {\begin{matrix} T_{11} \\ \begin{matrix} T_{22} \\ T_{12} \end{matrix} \end{matrix}} = 4 μ l^{2} [\begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{matrix}] {\begin{matrix} k_{11} \\ \begin{matrix} k_{22} \\ k_{12} \end{matrix} \end{matrix}}

(51)

\bar{U} = {\begin{matrix} u \\ v \\ θ \end{matrix}}

(52)

Where, $Π_{C}^{*}$ and $V_{C}^{*}$ indicate the complementary energy within the element and along the element boundaries, respectively. In addition, t, $σ$ , T, U are the thickness of the element, the stress vector, the surface traction force vector through the element boundaries, the displacement vector along element boundaries and the elastic flexibility matrix, respectively. Furthermore, for the plate, $μ = \frac{E}{2 (1 + ν)}$ , $λ = \frac{ν E}{(1 + ν) (1 - 2 ν)}$ . It should be reminded that $E$ and $μ$ are Young’s modulus and Poisson’s ratio, respectively. Therefore, the stress vector, $σ$ , can be expressed using the Airy stress function, $φ$ , as follows:

σ = {\begin{matrix} σ_{11} \\ σ_{22} \\ σ_{12}^{S} \end{matrix}} = {\begin{matrix} \frac{\partial^{2} φ}{\partial y^{2}} \\ \frac{\partial^{2} φ}{\partial x^{2}} \\ - \frac{\partial^{2} φ}{\partial x \partial y} \end{matrix}} = \tilde{R} (φ)

(53)

At this stage, the complementary energy function can be written in the following form:

Π_{C} = Π_{C}^{*} + V_{C}^{*} = \frac{1}{2} \int \int_{A^{e}}^{1} \tilde{R} {(φ)}^{T} C \tilde{R} (φ) tdA - \int_{Γ}^{1} (L \tilde{R} (φ))^{T} \bar{U} tds

(54)

Different parts of the last equation can be expressed by:

Π_{C}^{*} = \frac{1}{2} \int \int_{A^{e}}^{1} \tilde{R} {(φ)}^{T} C \tilde{R} (φ) dA

(55)

V_{C}^{*} = - \int_{Γ}^{1} (L \tilde{R} (φ))^{T} tds

(56)

Up to now, the element complementary energy functional containing the Airy stress function is established based on the fundamental relationships of the finite element method. In the plate bending problems without the body forces, the Airy stress function satisfies the following equation:

\nabla^{4} φ = \frac{\partial^{4} φ}{\partial y^{4}} + 2 \frac{\partial^{4} φ}{\partial x^{2} \partial y^{2}} + \frac{\partial^{4} φ}{\partial x^{4}}

(57)

In the first stage, the Airy stress function can be defined in terms of unknown parameters. A general form of this function is as follows:

φ = \sum_{i = 1}^{n} φ_{i} β_{i}

(58)

Where, n is the number of analytical solutions, and the other parts are given below:

φ = [φ_{1} φ_{2} \dots \dots φ_{n}] and β = [β_{1} β_{2} \dots \dots β_{n}]

(59)

$φ_{i}$ (i = 1 − n) are the number of analytical solutions for stress function and $β_{1}$ (i = 1 − n) (i ¼ 1–15) are the number of unknown constants. Upon substitution of equation (58) into equation (54), the subsequent equation will be achieved:

Π_{C}^{*} = \frac{1}{2} β^{T} M β

(60)

V_{C}^{*} = - β M β^{T}

(61)

M = \int \int_{A^{e}}^{1} S^{T} CStdA

(62)

Here, S and M are the matrix expressions. Having these matrices, the plate bending analysis will be achieved. Therefore, the matrix M can be written as follows:

M = \int_{- 1}^{1} \int_{- 1}^{1} S {(ξ_{1}, ξ_{2})}^{T} CS (ξ_{1}, ξ_{2}) t | j | d ξ_{1} d ξ_{2}

(63)

Where, $| j |$ is the Jacobian determinant. Besides, the following equation could be established:

V_{C}^{*} = - β^{T} H q^{e}

(64)

H = \int_{Γ}^{1} S^{T} L^{T} \bar{N} tdS

(65)

Matrices, $\bar{N}$ and H can be obtained according to the below relationships for TA13 and TA15:

\begin{matrix} \bar{N} 13 = [\begin{matrix} N_{1}^{0} & 0 & 0 \\ 0 & N_{1}^{0} & 0 \\ 0 & 0 & N_{1}^{0} \end{matrix} \begin{matrix} N_{2}^{0} & 0 & 0 \\ 0 & N_{2}^{0} & 0 \\ 0 & 0 & N_{2}^{0} \end{matrix} \begin{matrix} N_{3}^{0} & 0 & 0 \\ 0 & N_{3}^{0} & 0 \\ 0 & 0 & N_{3}^{0} \end{matrix} \begin{matrix} N_{4}^{0} & 0 & 0 \\ 0 & N_{4}^{0} & 0 \\ 0 & 0 & N_{4}^{0} \end{matrix} \\ \begin{matrix} N_{5}^{0} & 0 & 0 \\ 0 & N_{5}^{0} & 0 \\ 0 & 0 & N_{5}^{0} \end{matrix} \begin{matrix} N_{6}^{0} & 0 & 0 \\ 0 & N_{6}^{0} & 0 \\ 0 & 0 & N_{6}^{0} \end{matrix} \begin{matrix} N_{7}^{0} & 0 & 0 \\ 0 & N_{7}^{0} & 0 \\ 0 & 0 & N_{7}^{0} \end{matrix} \begin{matrix} N_{8}^{0} & 0 & 0 \\ 0 & N_{8}^{0} & 0 \\ 0 & 0 & N_{8}^{0} \end{matrix} \begin{matrix} N_{9}^{0} & 0 & 0 \\ 0 & N_{9}^{0} & 0 \\ 0 & 0 & N_{9}^{0} \end{matrix} \\ \begin{matrix} N_{10}^{0} & 0 & 0 \\ 0 & N_{10}^{0} & 0 \\ 0 & 0 & N_{10}^{0} \end{matrix} \begin{matrix} N_{11}^{0} & 0 & 0 \\ 0 & N_{11}^{0} & 0 \\ 0 & 0 & N_{11}^{0} \end{matrix} \begin{matrix} N_{12}^{0} & 0 & 0 \\ 0 & N_{12}^{0} & 0 \\ 0 & 0 & N_{12}^{0} \end{matrix} \begin{matrix} N_{13}^{0} & 0 & 0 \\ 0 & N_{13}^{0} & 0 \\ 0 & 0 & N_{13}^{0} \end{matrix} \end{matrix}

(66)

\begin{matrix} \bar{N} 15 = [\begin{matrix} N_{1}^{0} & 0 & 0 \\ 0 & N_{1}^{0} & 0 \\ 0 & 0 & N_{1}^{0} \end{matrix} \begin{matrix} N_{2}^{0} & 0 & 0 \\ 0 & N_{2}^{0} & 0 \\ 0 & 0 & N_{2}^{0} \end{matrix} \begin{matrix} N_{3}^{0} & 0 & 0 \\ 0 & N_{3}^{0} & 0 \\ 0 & 0 & N_{3}^{0} \end{matrix} \begin{matrix} N_{4}^{0} & 0 & 0 \\ 0 & N_{4}^{0} & 0 \\ 0 & 0 & N_{4}^{0} \end{matrix} \\ \begin{matrix} N_{5}^{0} & 0 & 0 \\ 0 & N_{5}^{0} & 0 \\ 0 & 0 & N_{5}^{0} \end{matrix} \begin{matrix} N_{6}^{0} & 0 & 0 \\ 0 & N_{6}^{0} & 0 \\ 0 & 0 & N_{6}^{0} \end{matrix} \begin{matrix} N_{7}^{0} & 0 & 0 \\ 0 & N_{7}^{0} & 0 \\ 0 & 0 & N_{7}^{0} \end{matrix} \begin{matrix} N_{8}^{0} & 0 & 0 \\ 0 & N_{8}^{0} & 0 \\ 0 & 0 & N_{8}^{0} \end{matrix} \begin{matrix} N_{9}^{0} & 0 & 0 \\ 0 & N_{9}^{0} & 0 \\ 0 & 0 & N_{9}^{0} \end{matrix} \\ \begin{matrix} N_{10}^{0} & 0 & 0 \\ 0 & N_{10}^{0} & 0 \\ 0 & 0 & N_{10}^{0} \end{matrix} \begin{matrix} N_{11}^{0} & 0 & 0 \\ 0 & N_{11}^{0} & 0 \\ 0 & 0 & N_{11}^{0} \end{matrix} \begin{matrix} N_{12}^{0} & 0 & 0 \\ 0 & N_{12}^{0} & 0 \\ 0 & 0 & N_{12}^{0} \end{matrix} \begin{matrix} N_{13}^{0} & 0 & 0 \\ 0 & N_{13}^{0} & 0 \\ 0 & 0 & N_{13}^{0} \end{matrix} \\ \begin{matrix} N_{14}^{0} & 0 & 0 \\ 0 & N_{14}^{0} & 0 \\ 0 & 0 & N_{14}^{0} \end{matrix} \begin{matrix} N_{15}^{0} & 0 & 0 \\ 0 & N_{15}^{0} & 0 \\ 0 & 0 & N_{15}^{0} \end{matrix}] \end{matrix}

(67)

These functions are required to develop new TUBA elements. According to the number of nodes and degrees of freedom in these elements, stress function can be determined. The newly developed stress-based elements are called TA13 and TA15. After searching, examining and matching, all stress function parameters are found and listed in Tables 1 and 2.

Table 1.

Basic analytical solutions of stress function and stresses for the plane problem using TA13.

TA13	$i$	1	2	3	4	5	6	7	8	9	10	11
	$φ_{i}$	$x^{2}$	$xy$	$y^{2}$	$x^{3}$	$x^{2} y$	$y^{2} x$	$y^{3}$	$x^{4}$	$x^{3} y$	$y^{2} x^{2}$	$y^{3} x$
	$σ_{11}$	0	0	2	0	0	2x	6 years	0	0	$2 x^{2}$	6xy
	$σ_{22}$	2	0	0	6x	2 years	0	0	$12 x^{2}$	6xy	$2 y^{2}$	0
	$σ_{12}^{S}$	0	−1	0	0	−2x	−2y	0	0	$- 3 x^{2}$	−4xy	$- 3 y^{2}$
	$i$	12	13	14	15	16	17	18	19	20	21	22
	$φ_{i}$	$y^{4}$	$x^{5}$	$x^{4} y$	$x^{3} y^{2}$	$y^{3} x^{2}$	$y^{4} x$	$x^{5}$	$x^{6}$	$x^{5} y$	$x^{4} y^{2} - y^{4} x^{2}$	$y^{5} x$
	$σ_{11}$	$12 x^{2}$	0	0	$4 x^{3}$	$6 y x^{2}$	12 $y^{2} x$	$20 y^{3}$	0	0	$2 x^{4} - 12 y^{2} x^{2}$	$20 y^{3} x$
	$σ_{22}$	0	$20 x^{3}$	12 $x^{2}$ y	$6 x y^{2}$	$4 y^{3}$	0	0	$30 x^{4}$	$20 x^{3}$ y	$12 x^{2} y^{2} - 2 y^{4}$	0
	$σ_{12}^{S}$	0	0	$- 4 x^{3}$	−6 $x^{2}$ y	−6 $y^{2} x$	$- 4 y^{3}$	0	0	$- 5 x^{4} y$	$8 x^{3} y - 8 y^{3} x$	$- 5 y^{4} x$
	$i$	23			24		25		26		27
		$φ_{i}$	$y^{6}$			$x^{4} y$		$x^{3} y^{2} - y^{3} x^{2}$		$y^{4} x$		$x^{5}$
		$σ_{11}$	$30 y^{4}$			0		$2 x^{3} - 6 y x^{2}$		12 $y^{2} x$		$20 y^{3}$
		$σ_{22}$	$0$			12 $x^{2}$ y		$6 x y^{2} - 2 y^{3}$		0		0
		$σ_{12}^{S}$	0			$- 4 x^{3}$		$6 x^{2} y - 6 y^{2} x$		$- 4 y^{3}$		0

Table 2.

Basic analytical solutions of stress function and stresses for the plane problem using TA15.

TA15

i

φ_{i}

x^{2}

xy

y^{2}

x^{3}

x^{2} y

y^{2} x

y^{3}

x^{4}

x^{3} y

y^{3} x

y^{4}

σ_{11}

6 years

6xy

12 x^{2}

σ_{22}

2 years

12 y^{2}

6xy

σ_{12}^{S}

−1

−2x

−2y

- 3 x^{2}

- 3 y^{2}

i

φ_{i}

y^{2} x^{2}

y^{3} x

y^{4}

x^{5}

x^{2} y + y^{2} x

x^{3} y + y^{3} x

x^{5}

x^{3} y^{2} + y^{3} x^{2}

x^{4} y

y^{5}

x^{6}

σ_{11}

2 x^{2}

6xy

12 x^{2}

2 x

6 yx

2 x^{3} + 6 y x^{2}

20 y^{3}

σ_{22}

2 y^{2}

20 x^{3}

x^{2} y + y^{2} x

6 xy

20 x^{3}

6 x y^{2} + 2 y^{3}

12 x^{2} y

30 x^{4}

σ_{12}^{S}

−4xy

- 3 y^{2}

−2x − 2y

- 3 x^{2} - 3 y^{2}

- 6 x^{2} y - 6 y^{2} x

- 4 x^{3}

i

φ_{i}

x^{5} y

x^{4} y^{2} + y^{4} x^{2}

x^{3} y^{3}

y^{5} x

y^{6}

σ_{11}

2 x^{4} + 12 y^{2} x^{2}

6 x^{3} y

20 y^{3} x

30 y^{3}

σ_{22}

20 x^{3} y

12 x^{2} y^{2} + 2 y^{4}

6 x y^{3}

σ_{12}^{S}

- 5 x^{4}

- 8 x^{3} y - 8 y^{3} x

- 9 x^{2} y^{2}

- 5 y^{4}

i

φ_{i}

7 x^{4} y - 5 y^{4} x + 3 xy

−15

x^{3} y^{3} + 7 x^{2} y^{2}

−2

x^{4} y^{3} + 7 x^{3} y^{4}

−5

x^{2} y^{3} + 4 x^{3} y^{2}

\frac{13}{2} x^{5} y - \frac{15}{4} x y^{5}

σ_{11}

−

60 y^{2} x

−90

x^{3} y + 14 x^{2}

−12

x^{4} y + 84 x^{3} y^{2}

−30

x^{2} y + 8 x^{3}

- 75 x y^{3}

σ_{22}

84 x^{4} y

−90

x y^{3} + 14 y^{2}

−24

x^{2} y^{3} + 42 x y^{4}

−10

y^{3} + 24 x y^{2}

130 x^{3} y

σ_{12}^{S}

- 28 x^{3} + 20 y^{3} - 3

135

x^{2} y^{2} - 28 xy

x^{3} y^{2} - 84 x^{2} y^{3}

x y^{2} - 24 x^{2} y

- \frac{65}{2} x^{4} + \frac{75}{4} y^{4}

i

φ_{i}

\frac{11}{3} x^{6} y - \frac{15}{4} x y^{6}

- 15 x^{7} + 12 y^{7}

- 7 x^{3} y^{4} + 5 y^{3} x^{3} - 8 y^{2} x^{2}

σ_{11}

- \frac{225}{2} x y^{4}

504 y^{5}

- 84 x^{3} y^{2} + 30 y x^{3} - 16 x^{2}

σ_{22}

110 x^{4} y

- 504 x^{5}

- 42 x y^{4} + 30 y^{3} x - 16 y^{2}

σ_{12}^{S}

- 22 x^{5} + \frac{45}{2} y^{5}

84 x^{2} y^{3} - 45 y^{2} x^{2} + 32 yx

At this stage, the degrees of freedom for the new elements are defined. As it is shown in Figure 2, the element nodal displacement vector, $q^{e}$ , for TA13 has the next form:

\begin{matrix} {q^{e}}_{13} = [w_{1} w_{x 1} w_{y 1} w_{xx 1} w_{xy 1} w_{yy 1} w_{2} w_{x 2} w_{y 2} w_{xx 2} \\ w_{xy 2} w_{yy 2} w_{3} w_{x 3} w_{y 3} w_{xx 3} w_{xy 3} w_{yy 3} w_{n 4} w_{n 5} w_{n 6} \\ w_{7} w_{8} w_{9} w_{n 10} w_{n 11} w_{n 12} w_{13}] \end{matrix}

(68)

TA15 has the following degrees-of-freedom:

\begin{matrix} {q^{e}}_{15} = [w_{1} w_{x 1} w_{y 1} w_{xx 1} w_{xy 1} w_{yy 1} w_{2} w_{x 2} w_{y 2} w_{xx 2} \\ w_{xy 2} w_{yy 2} w_{3} w_{x 3} w_{y 3} w_{xx 3} w_{xy 3} w_{yy 3} w_{x 4} w_{y 4} w_{x 5} w_{y 5} \\ w_{x 6} w_{y 6} w_{n 7} w_{n 8} w_{n 9} w_{x 10} w_{y 10} w_{x 11} w_{y 11} w_{x 12} w_{y 12} w_{13} \\ w_{14} w_{15}] \end{matrix}

(69)

Having these matrices, plate bending analyses will be obtained. After performing the required calculations, the next result for TA13 will be found:

\begin{matrix} S_{13} = [\begin{matrix} 0 & 0 & 2 \\ 2 & 0 & 0 \\ 0 & - 1 & 0 \end{matrix} \begin{matrix} 0 & 0 & 2 x \\ 6 x & 2 y & 0 \\ 0 & - 2 x & - 2 y \end{matrix} \begin{matrix} 6 y & 0 & 0 \\ 0 & 12 x^{2} & 6 xy \\ 0 & 0 & - 3 x^{2} \end{matrix} \begin{matrix} 2 x^{2} & 6 xy & 12 x^{2} \\ 2 y^{2} & 0 & 0 \\ - 4 xy & - 3 y^{2} & 0 \end{matrix} \\ \begin{matrix} 0 & 0 & 4 x^{3} \\ 20 x^{3} & 12 x^{2} y & 6 x y^{2} \\ 0 & - 4 x^{3} & - 6 x^{2} y \end{matrix} \begin{matrix} 6 y x^{2} & 12 y^{2} x & 20 y^{3} \\ 4 y^{3} & 0 & 0 \\ - 6 y^{2} x & - 4 y^{3} & 0 \end{matrix} \begin{matrix} 0 \\ 30 x^{4} \\ 0 \end{matrix} \begin{matrix} 0 \\ 20 x^{3} y \\ - 5 x^{4} y \end{matrix} \begin{matrix} 2 x^{4} - 12 y^{2} x^{2} \\ 12 x^{2} y^{2} - 2 y^{4} \\ 8 x^{3} y - 8 y^{3} x \end{matrix} \\ \begin{matrix} 20 y^{3} x \\ 0 \\ - 5 y^{4} x \end{matrix} \begin{matrix} 30 y^{4} \\ 0 \\ 0 \end{matrix} \begin{matrix} 0 \\ 12 x^{2} y \\ - 4 x^{3} \end{matrix} \begin{matrix} 2 x^{3} - 6 y x^{2} \\ 6 x y^{2} - 2 y^{3} \\ 6 x^{2} y - 6 y^{2} x \end{matrix} \begin{matrix} 12 y^{2} x \\ 0 \\ - 4 y^{3} \end{matrix} \begin{matrix} 20 y^{3} \\ 0 \\ 0 \end{matrix}] \end{matrix}

(70)

For TA15, the subsequent result is obtained:

\begin{matrix} S_{15} = [\begin{matrix} 0 & 0 & 2 \\ 2 & 0 & 0 \\ 0 & - 1 & 0 \end{matrix} \begin{matrix} 0 & 0 & 2 x \\ 6 x & 2 y & 0 \\ 0 & - 2 x & - 2 y \end{matrix} \begin{matrix} 6 y & 0 & 0 \\ 0 & 12 x^{2} & 6 xy \\ 0 & 0 & - 3 x^{2} \end{matrix} \begin{matrix} 6 xy & 12 x^{2} & 2 x^{2} \\ 0 & 0 & 2 y^{2} \\ - 3 y^{2} & 0 & - 4 xy \end{matrix} \\ \begin{matrix} 6 xy & 12 x^{2} & 0 \\ 0 & 0 & 20 x^{3} \\ - 3 y^{2} & 0 & 0 \end{matrix} \begin{matrix} 2 x & 6 yx & 0 \\ x^{2} y + y^{2} x & 6 xy & 20 x^{3} \\ - 2 x - 2 y & - 3 x^{2} - 3 y^{2} & 0 \end{matrix} \begin{matrix} 2 x^{3} + 6 y x^{2} \\ 6 x y^{2} + 2 y^{3} \\ - 6 x^{2} y - 6 y^{2} x \end{matrix} \begin{matrix} 0 \\ 12 x^{2} y \\ - 4 x^{3} \end{matrix} \\ \begin{matrix} 20 y^{3} & 0 & 0 \\ 0 & 30 x^{4} & 20 x^{3} y \\ 0 & 0 & - 5 x^{4} \end{matrix} \begin{matrix} 2 x^{4} + 12 y^{2} x^{2} & 6 x^{3} y & 20 y^{3} x \\ 12 x^{2} y^{2} + 2 y^{4} & 6 x y^{3} & 0 \\ - 8 x^{3} y - 8 y^{3} x & - 9 x^{2} y^{2} & - 5 y^{4} \end{matrix} \begin{matrix} 30 y^{3} \\ 0 \\ 0 \end{matrix} \begin{matrix} 60 y^{2} x \\ 84 x^{4} y \\ - 28 x^{3} + 20 y^{3} - 3 \end{matrix} \\ \begin{matrix} - 90 x^{3} y + 14 x^{2} & - 12 x^{4} y + 84 x^{3} y^{2} & - 30 x^{2} y + 8 x^{3} \\ - 90 x y^{3} + 14 y^{2} & - 24 x^{2} y^{3} + 42 x y^{4} & - 10 y^{3} + 24 x y^{2} \\ 135 x^{2} y^{2} - 28 xy & 24 x^{3} y^{2} - 84 x^{2} y^{3} & 30 x y^{2} - 24 x^{2} y \end{matrix} \begin{matrix} - 75 x y^{3} \\ 130 x^{3} y \\ - \frac{65}{2} x^{4} + \frac{75}{4} y^{4} \end{matrix} \\ \begin{matrix} - \frac{225}{2} x y^{4} & 504 y^{5} & - 84 x^{3} y^{2} + 30 y x^{3} - 16 x^{2} \\ 110 x^{4} y & - 504 x^{5} & - 42 x y^{4} + 30 y^{3} x - 16 y^{2} \\ - 22 x^{5} + \frac{45}{2} y^{5} & 0 & 84 x^{2} y^{3} - 45 y^{2} x^{2} + 32 yx \end{matrix}] \end{matrix}

(71)

Matrices of $\bar{N}$ can be obtained in the below form:

\begin{matrix} \bar{N} 13 = [\begin{matrix} N_{1}^{0} & 0 & 0 \\ 0 & N_{1}^{0} & 0 \\ 0 & 0 & N_{1}^{0} \end{matrix} \begin{matrix} N_{2}^{0} & 0 & 0 \\ 0 & N_{2}^{0} & 0 \\ 0 & 0 & N_{2}^{0} \end{matrix} \begin{matrix} N_{3}^{0} & 0 & 0 \\ 0 & N_{3}^{0} & 0 \\ 0 & 0 & N_{3}^{0} \end{matrix} \begin{matrix} N_{4}^{0} & 0 & 0 \\ 0 & N_{4}^{0} & 0 \\ 0 & 0 & N_{4}^{0} \end{matrix} \\ \begin{matrix} N_{5}^{0} & 0 & 0 \\ 0 & N_{5}^{0} & 0 \\ 0 & 0 & N_{5}^{0} \end{matrix} \begin{matrix} N_{6}^{0} & 0 & 0 \\ 0 & N_{6}^{0} & 0 \\ 0 & 0 & N_{6}^{0} \end{matrix} \begin{matrix} N_{7}^{0} & 0 & 0 \\ 0 & N_{7}^{0} & 0 \\ 0 & 0 & N_{7}^{0} \end{matrix} \begin{matrix} N_{8}^{0} & 0 & 0 \\ 0 & N_{8}^{0} & 0 \\ 0 & 0 & N_{8}^{0} \end{matrix} \begin{matrix} N_{9}^{0} & 0 & 0 \\ 0 & N_{9}^{0} & 0 \\ 0 & 0 & N_{9}^{0} \end{matrix} \\ \begin{matrix} N_{10}^{0} & 0 & 0 \\ 0 & N_{10}^{0} & 0 \\ 0 & 0 & N_{10}^{0} \end{matrix} \begin{matrix} N_{11}^{0} & 0 & 0 \\ 0 & N_{11}^{0} & 0 \\ 0 & 0 & N_{11}^{0} \end{matrix} \begin{matrix} N_{12}^{0} & 0 & 0 \\ 0 & N_{12}^{0} & 0 \\ 0 & 0 & N_{12}^{0} \end{matrix} \begin{matrix} N_{13}^{0} & 0 & 0 \\ 0 & N_{13}^{0} & 0 \\ 0 & 0 & N_{13}^{0} \end{matrix}] \end{matrix}

(72)

\begin{matrix} \bar{N} 15 = [\begin{matrix} N_{1}^{0} & 0 & 0 \\ 0 & N_{1}^{0} & 0 \\ 0 & 0 & N_{1}^{0} \end{matrix} \begin{matrix} N_{2}^{0} & 0 & 0 \\ 0 & N_{2}^{0} & 0 \\ 0 & 0 & N_{2}^{0} \end{matrix} \begin{matrix} N_{3}^{0} & 0 & 0 \\ 0 & N_{3}^{0} & 0 \\ 0 & 0 & N_{3}^{0} \end{matrix} \begin{matrix} N_{4}^{0} & 0 & 0 \\ 0 & N_{4}^{0} & 0 \\ 0 & 0 & N_{4}^{0} \end{matrix} \\ \begin{matrix} N_{5}^{0} & 0 & 0 \\ 0 & N_{5}^{0} & 0 \\ 0 & 0 & N_{5}^{0} \end{matrix} \begin{matrix} N_{6}^{0} & 0 & 0 \\ 0 & N_{6}^{0} & 0 \\ 0 & 0 & N_{6}^{0} \end{matrix} \begin{matrix} N_{7}^{0} & 0 & 0 \\ 0 & N_{7}^{0} & 0 \\ 0 & 0 & N_{7}^{0} \end{matrix} \begin{matrix} N_{8}^{0} & 0 & 0 \\ 0 & N_{8}^{0} & 0 \\ 0 & 0 & N_{8}^{0} \end{matrix} \begin{matrix} N_{9}^{0} & 0 & 0 \\ 0 & N_{9}^{0} & 0 \\ 0 & 0 & N_{9}^{0} \end{matrix} \\ \begin{matrix} N_{10}^{0} & 0 & 0 \\ 0 & N_{10}^{0} & 0 \\ 0 & 0 & N_{10}^{0} \end{matrix} \begin{matrix} N_{11}^{0} & 0 & 0 \\ 0 & N_{11}^{0} & 0 \\ 0 & 0 & N_{11}^{0} \end{matrix} \begin{matrix} N_{12}^{0} & 0 & 0 \\ 0 & N_{12}^{0} & 0 \\ 0 & 0 & N_{12}^{0} \end{matrix} \begin{matrix} N_{13}^{0} & 0 & 0 \\ 0 & N_{13}^{0} & 0 \\ 0 & 0 & N_{13}^{0} \end{matrix} \\ \begin{matrix} N_{14}^{0} & 0 & 0 \\ 0 & N_{14}^{0} & 0 \\ 0 & 0 & N_{14}^{0} \end{matrix} \begin{matrix} N_{15}^{0} & 0 & 0 \\ 0 & N_{15}^{0} & 0 \\ 0 & 0 & N_{15}^{0} \end{matrix}] \end{matrix}

(73)

Moreover, H has the next formula:

H = \int_{Γ_{ij}}^{1} S^{T} L^{T} \bar{N} tds + \int_{Γ_{jk}}^{1} S^{T} L^{T} \bar{N} tds + \int_{Γ_{kl}}^{1} S^{T} L^{T} \bar{N} tds

(74)

Where, $Γ_{ij}$ , $Γ_{jk}$ , and $Γ_{kl}$ indicate the element edges. The direction cosines of the outer normal of each element edge, $l,$ and $m$ , could be written as follows:

l = \frac{dy}{ds} . m = - \frac{dx}{ds}

(75)

By inserting equations (74) and (75) into equation (54), the subsequent element complementary energy function can be found:

Π_{C}^{*} = \frac{1}{2} β^{T} M β - β^{T} H q^{e}

(76)

To find the solution by using the principle of minimum complementary energy, $Π_{C}$ should be minimized:

\frac{\partial Π_{C}}{\partial β} = 0

(77)

After calculating the nodal displacement vector, $q^{e}$ , the unknown constant vector, $β$ can be achieved by the next relation:

β = M^{- 1} H q^{e}

(78)

Substitution of equation (78) into equation (74) yields:

Π_{C}^{*} = \frac{1}{2} {q^{e}}^{T} K^{*} q^{e}

(79)

K^{*} = (M^{- 1} H)^{T} H

(80)

In the last equation, matrix K* can be considered as the equivalent stiffness matrix. This matrix can be used in the way as the conventional finite element technique. After finding the element nodal displacement vector, $q^{e}$ , the element stresses can be written as:

σ = S M^{- 1} H q^{e}

(81)

Having the stress function for each element, stresses of all points will be in hand. In fact, the stress value at any point within the element can be determined by inserting the Cartesian coordinates of that point into equation (81).

Numerical examples

To show the accuracy of the proposed new elements, TA13 and 1315, nine different benchmark problems are analyzed for both stresses and displacements. Moreover, the obtained results are compared with those achieved using well-known displacement-based elements (TUBA13 and TUBA15). The findings will clearly demonstrate the qualities of the new suggested elements in analyses of plate bending structures. All of the used units are consistent.

Benchmark 1: In this structure, an L-shaped plate containing a hole is assessed.^46,47 This plate is shaped by one-quarter of a double-symmetric frame. Load is applied on the left side, as represented in Figure 3. Dirichlet boundary conditions associated with symmetry planes are imposed on the simply supported boundary at the right-hand side. Therefore, on this side, the rotation is free to occur while the transitional degree of freedom in two directions in the plane is fixed. The achieved numerical results are drawn in Table 3. Regarding this table, the obtained outcomes are compared with those presented by TUBA13 and TUBA15 .^46,47 These outcomes reveal that the new elements have an excellent performance in terms of both displacement and stress responses. This could be associated with the use of Airy’s stress function. Because, both traditional elements, TUBA13 and TUBA 15 were displacement-based establishment elements while TA13 and TA15 are stress-based elements. Additionally, the trial-and-error methodology was used in this study which leads to reducing the error between the exact values and those obtained using the new elements, TA13 and TA15. Besides, all findings demonstrate that increasing the number of degrees-of-freedom leads to raising the accuracy of elements. For both displacements and stresses, TA15 reaches near-exact solutions, relative to the other elements.

Table 3.

Stress and displacement outcomes of a square panel with a central circular cutout.

Element	Displacement	Stress	Displacement error (%)	Stress error (%)
TUBA13	0.9388	1.23747	1.53	0.73
TA13	0.9405	1.23981	1.33	0.54
TUBA15	0.9484	1.24631	0.52	0.03
TA15	0.9509	1.24650	0.26	0.01
Exact	0.9534	1.24657	-	-

Figure 3.

Geometries of the left: plate with hole plate.^46,47

Benchmark 2: In this section, another L-shaped plate is considered, as it is demonstrated in Figure 4.^46,47 This panel is considered one-quarter of a double-symmetric frame. For this structure, Dirichlet boundary conditions associated with symmetry planes are imposed on the lower boundary and the right-hand side. As a result, in two simple supports, the rotation is permitted at these two sides while the transitional degrees of freedom in the two plane directions are fixed. A constant load in the horizontal direction ${\bar{t}}_{1}^{P} = - L^{- 1}$ kN/mm is applied to the body. The obtained results are represented in Table 4. Besides, the relative error of the center L-shaped domain pate is represented in Figure 5. Regarding Table 4 and Figure 5, the new TA13 and TA15 have more accurate stress and displacement responses. It should be highlighted that the influences of utilizing a stress function scheme to enhance the performance of TA elements on displacement outcomes are more significant than the stress consequences. This might be connected to the application of Airy’s stress function. Due to the fact that TA13 and TA15 are stress-based elements while both traditional elements, TUBA13 and TUBA15, were establishment elements based on displacement. The error between the exact values and those produced utilizing the additional elements, TA13 and TA15, was also reduced in this study by applying the trial-and-error process.

Table 4.

Maximum value of stress and deformation.

	Mesh $N \times N$	$16 \times 16$	$32 \times 32$	$100 \times 100$	Exact	Error (%) ( $16 \times 16$ )	Error (%) ( $32 \times 32$ )	Error (%) ( $100 \times 100$ )
Stress	TUBA13	0.0155	0.014	0.011	0.010	55	40	10
	TA13	0.0141	0.0122	0.010		41	22	0
	TUBA15	0.0135	0.012	0.010		35	20	0
	TA15	0.0122	0.011	0.010		22	10	0
Displacement	TUBA13	0.131	0.128	0.121	0.120	9.2	6.6	0.8
	TA13	0.129	0.125	0.120		7.5	4.1	0
	TUBA15	0.128	0.126	0.120		6.6	5.0	0
	TA15	0.125	0.122	0.120		4.2	1.7	0

Figure 4.

Geometries of the center L-shaped domain pate.^46,47

Figure 5.

Relative error of the center L-shaped domain pate.

Benchmark 3: Figure 6 shows a quarter of a 4 * 4 panel with a circular cutout. Along x₁ = 4 and x₂ = 4, traction and double-traction, derived from the following axial symmetric displacement field are prescribed⁴⁶:

{\begin{matrix} u_{1} \\ u_{2} \end{matrix}} = \frac{P}{2 μ r f_{3}} [f_{1} - \frac{a}{1 - 2 v f_{2}}] {\begin{matrix} \cos θ \\ \sin θ \end{matrix}}

(82)

\begin{matrix} f_{1} = (v - 1) a^{4} + (2 v - 1) [2 l^{2} - r^{2} (v - 1)] a^{2} \\ - 4 (v - 1) r^{2} l^{2} (2 v - 1) k_{1} (\frac{a}{l}) \end{matrix}

(83)

f_{2} = [a^{2} - r^{2} (v - 1)] K_{0} (\frac{a}{l}) - 2 lr K_{1} (\frac{r}{l})

(84)

f_{3} = a^{2} + 4 l^{2} (v - 1) K_{1} (\frac{a}{l}) + 2 al (v - \frac{1}{2}) K_{0} (\frac{a}{l})

(85)

Figure 6.

A quarter of a square panel with a central circular cutout.

$(r, θ)$ Here, are the polar coordinates, $a = 1$ is the radius of the cutout, $ρ$ is radial stress at $r = \infty$ , and $K_{n}$ is the n-th order modified Bessel functions of the second kind. The obtained outcomes are reported in Table 5. According to the achieved consequences, the new triangular element, TA15, gives accurate results. The application of Airy’s stress function may be related to these accurate outcomes. Because TA13 and TA15 are stress-based elements while TUBA13 and TUBA15 were displacement-based establishment elements. Furthermore, the study’s use of the trial-and-error method helped to diminish the discrepancy between exact values and results achieved by employing the additional elements, TA13 and TA15.

Table 5.

Stress and displacement outcomes of a square panel with a central circular cutout.

Mesh	TUBA13	TA13	TUBA15	TA15
Normalized stress
$4 \times 4$	1.21	1.18	1.11	1.09
Error (%) ( $4 \times 4)$	21.0	18.0	11.0	9.0
$8 \times 8$	1.19	1.15	1.05	1.03
Error (%) ( $8 \times 8)$	19.0	15.0	5.0	3.0
$16 \times 16$	1.17	1.12	1.01	1.00
Error (%) ( $16 \times 16)$	17.0	12.0	1.0	0.0
Normalized displacement
$4 \times 4$	1.15	1.11	1.09	1.05
Error (%) ( $4 \times 4)$	15.0	11.0	9.0	5.0
$8 \times 8$	1.11	1.08	1.05	1.03
Error (%) ( $8 \times 8)$	11.0	8.0	5.0	3.0
$16 \times 16$	1.07	1.04	1.03	1.01
Error (%) ( $16 \times 16)$	7.0	4.0	3.0	1.0

Benchmark 4:Figure 7 illustrates a 60° rhombic plate, which is subjected to a uniformly distributed transverse loading q. The thickness (a), span width (a), and Poisson’s ratio (µ) are 0.1, 5, and 0.3, respectively. The plate edges, AB and DC, are hard and simply supported, and the other two sides are both free. Since all elements are skewed because of the plate’s rhombic shape, it is a suitable problem to assess the performances of the newly suggested elements in the oblique meshes. The achieved numerical results are represented in Table 6. Regarding this table, the authors’ elements, TA13, and TA15 give very good results against a variety of the elements in terms of different stress and deformation responses. According to this table, it can be seen that the suggested elements have high accuracy because of using the presented functions.

Table 6.

The dimensionless deflections and stress resultants.

	Mesh	$16 \times 16$	$32 \times 32$	$100 \times 100$	Exact	Error (%) ( $16 \times 16$ )	Error (%) ( $32 \times 32$ )	Error (%) ( $100 \times 100$ )
$\frac{w_{D} D}{q a^{4}}$	TUBA13	0.0062	0.0056	0.0047	0.0041	51.4	36.5	14.9
	TA13	0.0051	0.0048	0.0045		24.6	17.1	10.5
	TUBA15	0.0054	0.0050	0.0042		31.8	22.1	2.6
	TA15	0.0048	0.0042	0.0040		17.0	3.1	1.4
$\frac{M_{xC}}{q a^{2}}$	TUBA13	0.0524	0.0512	0.0487	0.0481	9.8	6.2	1.6
	TA13	0.0514	0.0501	0.0485		6.0	4.7	1.4
	TUBA15	0.0513	0.0504	0.0482		6.1	5.1	1.2
	TA15	0.0501	0.0486	0.0480		4.6	1.0	1.1
$\frac{M_{yC}}{q a^{2}}$	TUBA13	0.0538	0.0521	0.0487	0.0482	12.1	8.6	1.0
	TA13	0.0511	0.0489	0.0485		6.3	1.1	0.6
	TUBA15	0.0513	0.0502	0.0483		6.1	4.3	0.2
	TA15	0.0498	0.0486	0.0481		3.4	0.8	0.2
$\frac{T_{yB}}{qa}$	TUBA13	0.4267	0.4248	0.4221	0.4190	1.8	1.4	0.7
	TA13	0.4228	0.4202	0.4180		0.9	0.3	0.2
	TUBA15	0.4251	0.4239	0.4210		1.4	1.1	0.4
	TA15	0.4234	0.4221	0.4189		1.0	0.7	0.1

Figure 7.

The 60° skew plate with two opposite edges is hard simply-supported.²⁵

Benchmark 5: As shown in Figure 8, a square plate, with the span-thickness ratio a/h = 50 and Poisson’s ratio µ = 0.3 is evaluated. This plate is subjected to a uniformly distributed transverse load with the intensity of q. Its top and bottom edges are hard simply-supported, and the other two are both free or soft simply supported. For this structure, the obtained outcomes are compared with those presented by Kant and Gadgil.⁴⁶ Their results are employed as reference solutions here. In Table 7, the dimensionless displacements and resultants are calculated at selected listed points. Regarding this table, the proposed elements lead to high accuracy for a triangular element with 15 nodes, especially for a square plate. Once again, the answers demonstrate that the suggested elements are able to furnish the exact solutions for the plate-bending structures. Moreover, the stresses at the different points are represented in Figures 9 and 10. Regarding these figures, the authors’ new triangular elements, TA13 and TA15, anticipate the stress over the plate more accurately than those obtained using the previous elements. Utilizing Airy’s stress function may be related to this. Due to the fact that TA13 and TA15 are stress-based elements, TUBA13 and TUBA15 were both conventional elements that were established based on displacement. This study also employed the trial-and-error methodology, which helped to reduce the discrepancy between the exact values and those produced utilizing the additional elements, TA13 and TA15.

Table 7.

The dimensionless deflections and stress.

	Mesh $N \times N$	$100 \times 100$	Exact	Error (%)
$\frac{w_{D} D}{q a^{4}}$	TUBA13	0.1057	0.1050	0.7
	TA13	0.1053		0.3
	TUBA15	0.1051		0.1
	TA15	0.1049		0.1
$\frac{M_{xC}}{q a^{2}}$	TUBA13	0.0273	0.0268	1.9
	TA13	0.0270		0.7
	TUBA15	0.0269		0.4
	TA15	0.0267		0.3
$\frac{M_{yC}}{q a^{2}}$	TUBA13	0.1303	0.1220	6.6
	TA13	0.1298		6.3
	TUBA15	0.1240		1.6
	TA15	0.1228		0.6
$\frac{M_{yD}}{q a^{2}}$	TUBA13	0.1361	0.1300	4.6
	TA13	0.1335		2.7
	TUBA15	0.1320		1.5
	TA15	0.1308		0.6

Figure 8.

Square plate with two opposite edges hard simply-supported.

Figure 9.

Shear force distribution lengthways the edge DC.

Figure 10.

Twisting moment distribution lengthways of the lowest border AB.

Benchmark 6: In this section, a skew plate is evaluated. This plate is simply supported on all edges. As it is demonstrated in Figure 11, this plate is assessed under a uniform distributed load. The skew angle results in interior vertex angles of 30″ and 150″. The plate properties are given by E = 10 × 106, v = 0.3, side length a = 100, plate thickness t = 1.0 and uniform loading q = 1.0. In Table 8, the obtained results of the displacement at different locations of a plate are represented. High accuracy for the proposed elements is obtained from this study. It should be stated that the achieved outcomes are also compared with those obtained by conventional displacement-based elements for the “thin” plate.^46–48 Regarding Table 8, the proposed elements provide excellent results with the least cost of all the elements considered in this paper.

Table 8.

Skew cantilever plate.

	Deflection at location
	1	2	3	4	5	6
TUBA13	0.282	0.179	0.102	0.107	0.042	0.015
Error (%)	5.05	12.25	15.70	17.06	23.00	31.80
TA13	0.287	0.186	0.114	0.115	0.049	0.018
Error (%)	3.36	8.82	5.78	10.82	12.40	18.18
TUBA15	0.292	0.198	0.115	0.112	0.047	0.019
Error (%)	1.68	2.94	4.93	13.1	16.07	13.63
TA15	0.295	0.200	0.119	0.121	0.051	0.020
Error (%)	0.67	1.96	1.64	6.20	8.90	7.91
Exact	0.297	0.204	0.121	0.129	0.056	0.022

Figure 11.

Typical mesh for Morley skew plate example.

Benchmark 7: Figure 12 shows one quadrant of a circular plate. Outcomes for the clamped and simply supported boundary conditions are computed for thick (R/t = 5/l) and relatively thin (R/t = 5/0. l) plates. Besides, simply supported boundary conditions are considered. The material characteristics of the plate are E = 10.92, and v = 0.3, which gives a plate stiffness, D, of $t^{3}$ , and the uniform loading q is 1.0. Each $Δ θ^{k}$ is also restrained on a boundary segment.⁴⁶ The obtained results are compared with those presented by TUBA13 and TUBA15, and they are shown in Table 9. According to these consequences, both new elements, TA13 and TA15, can improve the answers considerably. This might be related to the application of Airy’s stress scheme. Because TA13 and TA15 are stress-based elements, while TUBA13 and TUBA15 were displacement-based elements. Moreover, the trial-and-error process was applied in this investigation, which helped to reduce the error between the precise values and those obtained utilizing the new elements, TA13 and TA15.

Table 9.

Responses of a skew plate.

Element	$w (\times 10^{2} q l^{4} / D)$	$M_{y} (\times 10 q l^{2})$
TUBA13	0.7824	0.9519
Error (%)	1.52	0.73
TA13	0.7897	0.9521
Error (%)	0.60	0.70
TUBA15	0.7889	0.9575
Error (%)	0.70	0.14
TA15	0.7921	0.9581
Error (%)	0.30	0.08
Exact	0.7945	0.9589

Figure 12.

Typical mesh f2or quadrant of a circular plate.

Benchmark 8: In this example, the action of a sinusoidal load on two parallel sides of the original rectangular plate is considered.⁴⁷ Specifically, it is considered the scaled moments of the “symmetrical case” corresponding to the term for m = 7 in the definition of the transverse displacement following formula:

w = \frac{(1 - v^{2})}{E h^{3}} (\frac{π bm}{2 a} \tanh (\frac{π bm}{2 a}) \cosh (\frac{π m \bar{y}}{a}) - \frac{π m \bar{y}}{a} \sinh (\frac{π m \bar{y}}{a})) \sin (\frac{π mx}{a})

(86)

where $\bar{y} = y - b / 2$ to adjust this solution to the frame in Figure 13, a = 4 and b = 2. This corresponds to having normal moments equal to

\begin{matrix} m_{yy} = \frac{π^{2} m^{2}}{6 a^{3}} \cosh (\frac{π bm}{2 a}) \sin (\frac{π mx}{a}) \\ \to 614.979305228016 \sin (5.49778714378214 x) \end{matrix}

(87)

Figure 13.

Double symmetry simplification of the rectangular plate with a sinusoidal action.

At the simply supported sides with uniform y. Therefore, at these edges, the turning is free to occur while the transitional degree of freedom in two directions in the plane is fixed. The distribution of applied moments/imposed rotation is illustrated in Figure 13. When this problem is solved with the suggested elements, the corresponding solutions are not strictly equilibrated, because the polynomial moment field cannot be exactly represented for the sinusoidal load. Nevertheless, it is possible to obtain locally equilibrated solutions by considering the problem with imposed displacements, which has the same exact solution, i.e., the sides with uniform y become clamped. In fact, the normal rotation corresponding to the exact displacements is here imposed as previously given by:

\begin{matrix} θ_{y} = \\ \mp \frac{(1 - v^{2})}{E h^{3}} \frac{π m}{2 a} sech (\frac{π bm}{2 a}) (\frac{π bm}{a} + \sinh (\frac{π bm}{a})) \sin (\frac{π mx}{a}) \\ \to \mp 610.957276165364 \sin (5.49778714378214 x) \end{matrix}

(88)

The values obtained for the different elements are presented in Table 10. It is obvious that TA13 and TA15 have better performance than TUBA13 and TUBA15 elements.

Table 10.

Rectangular plate. Computed moments at the center of the original plate.

		$m_{xx}$	$m_{yy}$	$m_{xy}$
UDL	TUBA13	2.23187E-01	4.86284E-01	1.97245E-03
Error (%)		3.91	1.51	9.82
Line		2.72453E-01	5.98578E-01	3.78954E-02
Error (%)		8.26	2.81	1.63
Point		1.25687E-01	1.94578E-01	1.27875E-02
Error (%)		10.59	8.90	11.42
UDL	TA13	2.28121E-01	4.89579E-01	1.99854E-02
Error (%)		1.79	0.84	7.14
Line		2.89142E-01	6.05748E-01	3.80257E-02
Error (%)		2.64	1.64	1.30
Point		1.34258E-01	2.00687E-01	1.30304E-02
Error (%)		4.50	6.04	9.74
UDL	TUBA15	2.31456E-01	4.92531E-01	2.14571E-02
Error (%)		0.35	0.25	0.30
Line		2.96654E-01	6.02568E-01	3.84139E-02
Error (%)		0.11	2.16	0.29
Point		1.37895E-01	2.12575E-01	1.43278E-02
Error (%)		1.91	0.47	0.75
UDL	TA15	2.32156E-01	4.93547E-01	2.15174E-02
Error (%)		0.05	0.04	0.02
Line		2.96921E-01	6.14875E-01	3.84980E-02
Error (%)		0.02	0.16	0.07
Point		1.40153E-01	2.13239E-01	1.44270E-02
Error (%)		0.30	0.16	0.06
UDL	Exact	2.32280E-01	4.93748E-01	2.15210E-02
Line		2.96987E-01	6.15874E-01	3.85248E-02
Point		1.40581E-01	2.13587E-01	1.44358E-02

Benchmark 9: In this example, a simply supported rectangular plate is evaluated, as illustrated in Figure 14.⁴⁶ It should be stated that the moments under the point load are infinite. Table 11 shows the obtained outcomes using TA13, TA15, TUBA13, and TUBA15 elements. The following stress resultants are directly obtained from the compatible curvatures, using the constitutive relations, and their derivatives for the shear force. Although, standard stress smoothing/recovery procedures could be applied, in this study the raw results are presented, which are those that must be used for the computation of error bounds. The fact is that the twisting moment at the intersection of the two perpendicular sliding supports is zero because those supports imply that the twisting curvature is zero. Hence, the fundamental condition is of a kinematic nature, and it is strongly respected by any compatible solution. According to Table 11, the improvement influences of the new proposed element based on the complementary energy scheme on the displacement and stress components outcomes are more significant, relative to those of displacement-based elements, TUBA13 and TUBA15.

Table 11.

Stress and displacement outcomes of a rectangular plate.

Element	Displacement	Stress	Displacement error (%)	Stress error (%)
TUBA13	1.22054	1.85620	1.52	0.73
TA13	1.23015	1.86152	0.74	0.44
TUBA15	1.22794	1.85328	0.92	0.88
TA15	1.23901	1.86811	0.03	0.09
Exact	1.23942	1.86985	-	-

Figure 14.

Rectangular plate: symmetry simplification, geometry, mechanical characteristics, boundary conditions, and initial finite element mesh.

Conclusion

In this study, two new elements were proposed for analyzing the plate bending problems. One of them had 13, and the other one had 15 nodes. Both elements were formulated based on Airy’s functions. To build these elements, a complementary energy function was employed within the element. In this energy expression, the Airy stress function was applied as a functional variable. In addition, some basic analytical solutions were assigned for the stress functions. These trial functions were matched with each element’s number of DOF. After some mathematical operations, the element equation was established. To demonstrate the accuracy of the presented elements, several benchmarks were solved. Comparison studies were also carried out with two related compatible and displacement-based elements of TUBA13 and TUBA15. All numerical solutions indicated that accurate results could be obtained for the displacements, as well as; the stress by using these new elements.

In this study, the influence of shear stress is not considered which is one of the limitations of this study. Additionally, the proposed methods are proper for thin plates. So, it is recommended for future studies to measure and develop the Airy stress function for thick plates considering the shear stress influence.

Footnotes

Acknowledgements

The authors are thankful to Professor Mohammad Rezaiee-Pajand, a Faculty member of the Ferdowsi University of Mashhad, Iran, for his unstinting help, constructive suggestions, and scientific advice.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

ORCID iD

Ehsan Noroozinejad Farsangi

References

Milroy

Hinduja

Davey

The elastostatic three-dimensional boundary element method: analytical integration for linear isoparametric triangular elements. Appl Math Model 1997; 21(12): 763–782.

Srivastava

Contractor

DN.

Efficient evaluation of integrals in three-dimensional boundary element method using linear shape functions over plane triangular elements. Appl Math Model 1992; 16(6): 282–290.

Kumar

Chakrabarti

Bhargava

, et al. Probabilistic failure analysis of laminated sandwich shells based on higher order zigzag theory. J Sandwich Struct Mater 2015; 17(5): 546–561.

Kumar

Chakrabarti

Bhargava

Accurate dynamic response of laminated composites and sandwich shells using higher order zigzag theory. Thin-Walled Struct 2014; 77: 174–186.

Tan

Bettess

JA.

Elastic buckling of isotropic triangular flat plates by finite elements. Appl Math Model 1983; 7(5): 311–316.

Remešíková

Šagát

Novysedlák

Discrete Lagrangian algorithm for finding geodesics on triangular meshes. Appl Math Model 2019; 76: 396–427.

Kumar

Chakrabarti

Bhargava

Vibration of laminated composites and sandwich shells based on higher order zigzag theory. Eng Struct 2013; 56: 880–888.

Kumar

Chakrabarti

Bhargava

Finite element analysis of laminated composite and sandwich shells using higher order zigzag theory. Compos Struct 2013; 106: 270–281.

Kumar

Chakrabarti

Ketkar

Analysis of laminated composite skew shells using higher order shear deformation theory. Lat Am J Solids Struct 2013; 10: 391–919.

10.

Yekkalam Tash

Navayi Neya

. An analytical solution for bending of transversely isotropic thick rectangular plates with variable thickness. Appl Math Model 2020; 77(2): 1582–1602.

11.

Magnucka-Blandzi

Magnucki

Stawecki

Bending and buckling of a circular plate with symmetrically varying mechanical properties. Appl Math Model 2021; 89: 1198–1205.

12.

Zhang

Wang

Challamel

, et al. Uncovering the finite difference model equivalent to Hencky bar-net model for axisymmetric bending of circular and annular plates. Appl Math Model 2018; 61: 300–315.

13.

Karimipour

Ghalehnovi

de Brito

New model for the lap-splice length of tensile reinforcement in concrete elements. J Struct Eng 2021; 147.

14.

Ouyang

Jiang

, et al. Integration modified wavelet neural networks for solving thin plate bending problem. Appl Math Model 2013; 37(5): 2983–2994.

15.

Levyakov

Kuznetsov

VV.

Complete system of two-dimensional invariants for formulation of triangular finite element of composite shells. Appl Math Model 2012; 36(11): 5183–5198.

16.

Han

Ren

Huang

A wavelet-based stochastic finite element method of thin plate bending. Appl Math Model 2007; 31(2): 181–193.

17.

Choo

Choi

Lee

BC.

A new hybrid-Trefftz triangular and quadrilateral plate elements. Appl Math Model 2010; 34(1): 14–23.

18.

Dey

Haldar

Sengupta

, et al. An efficient plate element for the vibration of composite plates. Appl Math Model 2016; 40(9–10): 5589–5604.

19.

Dehghan

Sabouri

A spectral element method for solving the Pennes bioheat transfer equation by using triangular and quadrilateral elements. Appl Math Model 2012; 36(12): 6031–6049.

20.

Zhang

Fan

Zhao

, et al. Fundamental solutions of uniform loads over triangular elements in a three-dimensional piezoelectric medium. Appl Math Model 2014; 38(19–20): 4778–4795.

21.

Geng

Zhang

Chen

, et al. Mid-frequency dynamic characteristics prediction of thin plate based on B-spline wavelet on interval finite element method. Appl Math Model 2018; 62: 526–541.

22.

Shirmohammadi

Bahrami

Dynamic response of circular and annular circular plates using spectral element method. Appl Math Model 2018; 53: 156–166.

23.

Demasi

Hong

Santarpia

Reissner’s mixed variational theorem and triangular finite element discretizations: an energetic interpretation. Thin-Walled Struct 2021; 166: 107994.

24.

Kamiński

Uncertainty analysis in solid mechanics with uniform and triangular distributions using stochastic perturbation-based finite element method. Finite Elem Anal Des 2022; 200: 103648.

25.

Zhang

Convergence and supercloseness of a finite element method for a two-parameter singularly perturbed problem on Shishkin triangular mesh. Appl Math Comput 2022; 416: 126753.

26.

Zhou

A family of quadratic finite volume element schemes for anisotropic diffusion problems on triangular meshes. J Comput Appl Math 2022; 402: 113794.

27.

Nie

A monotone finite volume element scheme for diffusion equations on triangular grids. Comput Math Appl 2022; 105: 1–12.

28.

Liu

Yang

Error estimations in the balanced norm of finite element method on Bakhvalov–Shishkin triangular mesh for reaction–diffusion problems. Appl Math Lett 2022; 123: 107523.

29.

Selzer

Allgeier

Therrien

, et al. Finite-volume flux reconstruction and semi-analytical particle tracking on triangular prisms for finite-element-type models of variably-saturated flow. Adv Water Resour 2021; 154: 103944.

30.

Wang

Song

Development and implementation of new triangular finite element based on MGE theory for bi-material analysis. Results Phys 2019; 13: 102231.

31.

Zhang

A numerical scheme with divergence free H-div triangular finite element for the Stokes equations. Appl Numer Math 2021; 167: 211–217.

32.

Cui

Duan

Huo

, et al. A high order cell-based smoothed finite element method using triangular and quadrilateral elements. Eng Anal Bound Elem 2021; 128: 133–148.

33.

Nguyen

Ibrahimbegovic

Hybrid-stress triangular finite element with enhanced performance for statics and dynamics. Comput Methods Appl Mech Eng 2020; 372: 113381.

34.

Ebrahimi

Edalati

Valizadeh

, et al. Increase the effectiveness of AMTMDs and PMTMDs on the seismic behaviour of structures case study: ten-stories short period concrete building. Eng Struct 2021; 237: 112122.

35.

Karimipour

Ghalehnovi

Edalati

, et al. Investigation responses of the diagrid structural system of high-rise buildings equipped with tuned mass damper using new dynamic method. J Build Mater Sci 2021; 1(2).

36.

Garg

Donough

Song

, et al. Three-dimensional scaled boundary finite element method to study interfacial imperfections in thick laminated composite plates undergoing bi-axial bending. Eng Anal Bound Elem 2023; 146: 34–43.

37.

Yin

Mei

Finite element implementation of general triangular mesh for Riesz derivative. Partial Differ Equat Appl Math 2021; 4: 100188.

38.

Zhou

A family of quadratic finite volume element schemes over triangular meshes for elliptic equations. Comput Math Appl 2020; 79(9): 2473–2491.

39.

Ghalehnovi

Yousefi

Karimipour

, et al. Investigation of the behaviour of steel-concrete-steel sandwich slabs with Bi-directional corrugated-strip connectors. Appl Sci 2020; 10: 8647.

40.

Liu

Zhang

, et al. A new semi-analytical solution of bending, buckling and free vibration of functionally graded plates using scaled boundary finite element method. Thin-Walled Struct 2021; 163: 107776.

41.

Garg

Prusty

Song

, et al. Cylindrical bending of thick laminated composite plates using scaled boundary finite element method. Eng Anal Bound Elem 2020; 120: 73–81.

42.

Wang

Ren

, et al. A scaled boundary finite element method for bending analysis of fiber-reinforced piezoelectric laminated composite plates. Int J Mech Sci 2019; 161–162: 105011.

43.

Moitinho

Almeida

Maunder

Tiago

A general degree semi-hybrid triangular compatible finite element formulation for Kirchhoff plates. Int J Numer Methods Eng 2019; 120: 56–85.

44.

Fleck

Muller

Ashby

, et al. Strain gradient plasticity: theory and experiment. Acta Metallurgica Mater 1994; 42(2): 475–487.

45.

Argyris

Fried

Scharpf

DW.

The TUBA family of plate elements for the matrix displacement method. Aeronaut J 1968; 72(692): 701–709.

46.

Haitao

Chen

ZC.

Euler–Bernoulli beam under arbitrary dynamic loads. J Mech Phys Solids 2012; 57: 1835–1867.

47.

Chen

Yang

Guo

, et al. An energy finite element method for high frequency vibration analysis of beams with axial force. Appl Math Model 2018; 61: 521–539.

48.

Stricklin

Richardson

, et al. On isoparametricvs linear strain triangular elements. Int J Numer Methods Eng 1977; 11(6): 1041–1043.

Representing capabilities of novel semi-analytical triangular plate elements

Abstract

Keywords

Introduction

Elements formulation

Numerical examples

Conclusion

Footnotes

Acknowledgements

Declaration of conflicting interests

Funding

ORCID iD

References