# Asymptotically accurate and locking-free finite element implementation of
first order shear deformation theory for plates

K. C. Lea,b111Corresponding author. Phone: +84 93 4152458, email:
<EMAIL_ADDRESS>H. G. Buic aInstitute for Advanced Study in
Technology, Ton Duc Thang University, Ho Chi Minh City, Vietnam
bFaculty of Civil Engineering, Ton Duc Thang University, Ho Chi Minh City,
Vietnam
cHelmholtz-Zentrum Hereon, 21502 Geesthacht, Germany

###### Abstract

A formulation of the asymptotically exact first-order shear deformation theory
for linear-elastic homogeneous plates in the rescaled coordinates and angles
of rotation is considered. This allows the development of its asymptotically
accurate and shear-locking-free finite element implementation. As
applications, numerical simulations are performed for circular and rectangular
plates, showing complete agreement between the analytical solution and the
numerical solutions based on two-dimensional theory and three-dimensional
elasticity theory.

###### keywords:

first-order shear deformation theory, plates, finite element, asymptotic
accuracy, shear-locking.

††journal: Comput. Methods Appl. Mech. Eng.

##  1 Introduction

The first-order shear deformation theory (FSDT) for plates, originally
proposed by Reissner [1],222See also Mindlin [2], whose motivation for
including transverse shear and rotatory inertia was somewhat different: his
main goal was to capture both low- and high-frequency thickness vibrations of
the plates. has since attracted the attention of both theorists and
practitioners alike, mainly because of its applicability to moderately thick
plates and the development of numerical methods [3]-[10], but also because of
the logic behind its derivation, which can be applied to other problems
[11]-[13]. The first asymptotically exact version of FSDT for plates has been
derived by Berdichevsky [14] using the variational-asymptotic method he
himself developed. The extension of his result to laminated plates using the
same method has been considered by Sutyrin [15] and Yu [16]. However, since in
the general case of laminated plates, the dimension reduction does not lead to
a FSDT, these authors have tried to optimize the parameters so that a derived
theory is as close as possible to asymptotic correctness while being a FSDT.
Le [17] has recently shown that the construction of asymptotically exact FSDT
for functionally graded (FG) plates is possible when the mass density and
elastic moduli vary across thickness such that their distributions are even
functions of the transverse coordinate. Similar to Berdichevsky’s FSDT for
homogeneous plates, his theory for FG-plates is asymptotically exact up to the
order of $h^{2}/l^{2}$, where $h$ is the plate thickness and $l$ is the
characteristic scale of change of the deformation state in the longitudinal
directions. Also worth mentioning are some recent applications of the
variational-asymptotic method to dimension reduction, homogenization,
nonlinear vibrations, and wave propagation in [18]-[23].

In numerical simulations of bending deformation of plates under a transverse
load within FSDT based on the finite element method, the so-called shear-
locking (SL) effect often occurs, especially when using low-order finite plate
elements (see the above cited papers [3]-[10]). Physically, this effect is due
to the fact that the shear stiffness, in terms of the small plate thickness
$h$, is two orders of magnitude larger than the bending stiffness. On the
other hand, the rotation angles, which are obviously of the order of the
characteristic strain, the angles caused by pure shear being even smaller, are
much less than the bending measures (or curvature changes) of the plate, since
the latter are of the order of the characteristic strain divided by the
thickness. This is also evident from the variational-asymptotic analysis of
the energy functional [24, 25]. Therefore, the numerical instability with
respect to the shear energy and also to the constitutive equations must occur
in the limit $h\to 0$ due to the multiplication of the extremely small and
large numbers when the standard low-order finite element calculation is used.
There are several sophisticated methods that make it possible to alleviate
this shear-locking effect but at the expense of computational efficiency. The
reduced and selective integration method, using two different integration
rules for bending and shear energies [26]-[29], has long been the preferred
technique for the numerical treatment of SL. The mathematical justification of
this method, based on the equivalence between the reduced integration approach
and certain mixed models, was later given by Malkus and Hughes [30].
Unfortunately, the reduced integration technique often leads to instability
due to rank deficiency and to zero-energy modes [31]-[33]. Therefore, several
alternative formulations and numerical techniques have been developed to
mitigate SL and increase the accuracy and stability of the solution. These
include the modified shear strain method [34, 35], the hybrid and mixed method
[36]-[39], the extended assumed strain method [40, 41], the assumed natural
strain method [34, 42], and the shear gap method [5] (see also the recent
paper [6] discussing relevant publications).

As far as the authors of this paper are aware, none of the existing studies on
this topic has found a simple formulation of the FSDT that is inherently free
of the shear-locking effect, regardless of the discretization scheme and
integration technique used. The goal of this paper is therefore twofold.
First, we give the formulation of the FSDT for plates in the rescaled
coordinates and rotation angles. This formulation occurs naturally when the
coordinates in the mid-plane are scaled by the plate thickness $h$ (thus
becoming dimensionless), while the rotation angles are multiplied by $h$,
resulting in equal and finite orders of the bending and shear stiffnesses as
well as the scaled rotation angles and bending measures.333Compare with [43],
where the intrinsically shear-locking-free formulation was obtained by
reparametrizing the kinematic equations for bending measures and shear. Since
this formulation is independent of the plate thickness and inherently shear-
locking-free, no high-order interpolation scheme and/or sophisticated
integration technique is required for the discretization and FE-
implementation, so the computational efficiency can be significantly improved.
However, our second goal is to ensure that the FE-implementation is
asymptotically accurate. According to [14, 17], the asymptotic accuracy
requires that both the transverse displacement and the rotation angles belong
to the $C^{1}$-function space. We will show by numerical simulations of the
circular and square plates and comparison with the analytical solutions of the
FSDT and the numerical solution of the three-dimensional elasticity theory
that the asymptotic accuracy is indeed achieved, provided that the
isogeometric elements guaranteeing the $C^{1}$-continuity for the primary
variables are used.

The paper is organized as follows. After this brief introduction, Section 2
gives the inherently shear-locking-free rescaled variational formulation of
the FSDT for plates. Section 3 is devoted to its FE-implementation. In Section
4, we consider two numerical examples where the developed FE-code is applied:
(i) clamped circular plates subjected to uniform loading, and (ii) rectangular
plates with one clamped edge and three free edges subjected to uniform
loading. Finally, Section 5 concludes the paper.

##  2 Rescaled variational formulation of FSDT for plates

Let $\Omega$ be a two-dimensional domain in the $(x_{1},x_{2})$-plane bounded
by a piecewise smooth closed curve $\partial\Omega$. We consider a linear
elastic homogeneous plate which in the undeformed stress-free state occupies
the 3-D region $\mathcal{V}=\Omega\times(-h/2,h/2)$. Its cross section in the
plane $(x_{1},x_{3})$ is shown in Fig. 1. We call $\Omega$ the plate mid-plane
and $h$ its thickness.

Figure 1: Cross section of a plate

Variational principle of the asymptotically exact two-dimensional first-order
shear deformation theory (FSDT) for linearly elastic homogeneous plate [14,
25] states that the true deflection and rotation angles of the plate minimize
the 2-D average energy functional

$J[u,\varphi_{\alpha}]=\int_{\Omega}\Bigl{\\{}\frac{\mu
h^{3}}{12}[\sigma(\rho_{\alpha\alpha})^{2}+\rho_{\alpha\beta}\rho_{\alpha\beta}]+\frac{5\mu
h}{12}\varphi_{\alpha}\varphi_{\alpha}\Bigr{\\}}\differential{a}-\int_{\Omega}\Bigl{(}fu+\frac{\sigma
h^{2}}{10}f\rho_{\alpha\alpha}\Bigr{)}\differential{a}$

among all deflections and rotation angles satisfying the kinematic boundary
conditions, with $\differential a=\differential x_{1}\differential x_{2}$ and
$\sigma=\frac{\lambda}{\lambda+2\mu}=\frac{\nu}{1-\nu}$. Note that
$\varphi_{\alpha}$ are the rotation angles due to pure shear deformation. Here
$\rho_{\alpha\beta}$ are the measures of bending and $f$ is the external
transverse load. They are given by

$\displaystyle\rho_{\alpha\beta}=u_{,\alpha\beta}-\varphi_{\alpha,\beta},$
$\displaystyle f=\tau_{3}|_{x_{3}=h/2}+\tau_{3}|_{x_{3}=-h/2},$

with $\tau_{3}$ being the normal traction. In the following, we use Greek
indices, running from 1 to 2, to refer to the plane coordinates $x_{1}$ and
$x_{2}$. The comma before an index denotes differentiation with respect to the
corresponding coordinate, the parentheses surrounding a pair of indices denote
symmetrization operation and summation over repeated indices is understood.

The FE-implementation of this variational problem is not straightforward due
to the different magnitudes of the bending and shear stiffnesses. Since the
shear stiffness $5\mu h/12$ is two orders of magnitude larger than the bending
stiffness $\mu h^{3}/12$ with respect to $h$, its multiplication by very small
rotation angles leads to a numerical instability called the shear-locking
effect when $h$ becomes small. To solve this problem, we first try to get rid
of the second derivatives in this variational problem. We introduce the new
unknown functions

$\psi_{\alpha}=-u_{,\alpha}+\varphi_{\alpha}$

which have the meaning of the total angles of rotation of the transverse
fibers due to both bending and shear deformation so that

$\displaystyle\rho_{\alpha\beta}=u_{,\alpha\beta}-\varphi_{\alpha,\beta}=-\psi_{(\alpha,\beta)},$
$\displaystyle\varphi_{\alpha}=u_{,\alpha}+\psi_{\alpha}.$

In terms of these new unknown functions, the functional becomes

$J[u,\psi_{\alpha}]=\int_{\Omega}\Bigl{\\{}\frac{\mu
h^{3}}{12}[\sigma(\psi_{\alpha,\alpha})^{2}+\psi_{(\alpha,\beta)}\psi_{(\alpha,\beta)}]+\frac{5\mu
h}{12}(u_{,\alpha}+\psi_{\alpha})(u_{,\alpha}+\psi_{\alpha})\Bigr{\\}}\differential{a}\\\
-\int_{\Omega}\Bigl{(}fu-\frac{\sigma
h^{2}}{10}f\psi_{\alpha,\alpha}\Bigr{)}\differential{a}.$

To avoid the shear-locking effect, we try to make the problem independent of
$h$ so that the orders of magnitude of the bending and shear stiffnesses
become the same. For this purpose, we introduce the rescaled coordinates and
rotation angles as follows

$\bar{x}_{\alpha}=\frac{x_{\alpha}}{h},\quad\bar{\psi}_{\bar{\alpha}}=h\psi_{\alpha}.$

Note that $\bar{x}_{\alpha}$ are dimensionless, while
$\bar{\psi}_{\bar{\alpha}}$ have the dimension of length and can be
interpreted as the longitudinal displacements of the positive face surface of
the plate [14, 17]. Since

$\displaystyle\rho_{\alpha\beta}=-\psi_{(\alpha,\beta)}=-\frac{1}{h^{2}}\bar{\psi}_{(\bar{\alpha},\bar{\beta})},$
$\displaystyle
u_{,\alpha}+\psi_{\alpha}=\frac{1}{h}(u_{,\bar{\alpha}}+\bar{\psi}_{\bar{\alpha}}),$
$\displaystyle\differential{a}=\differential{x_{1}}\differential{x_{2}}=h^{2}\differential{\bar{x}_{1}}\differential{\bar{x}_{2}}=h^{2}\differential{\bar{a}},$

we reduce the energy functional to

$J[u,\bar{\psi}_{\bar{\alpha}}]=h\int_{\bar{\Omega}}\Bigl{\\{}\frac{\mu}{12}[\sigma(\bar{\psi}_{\bar{\alpha},\bar{\alpha}})^{2}+\bar{\psi}_{(\bar{\alpha},\bar{\beta})}\bar{\psi}_{(\bar{\alpha},\bar{\beta})}]+\frac{5\mu}{12}(u_{,\bar{\alpha}}+\bar{\psi}_{\bar{\alpha}})(u_{,\bar{\alpha}}+\bar{\psi}_{\bar{\alpha}})\Bigr{\\}}\differential{\bar{a}}\\\
-\int_{\bar{\Omega}}h^{2}\Bigl{(}fu-\frac{\sigma}{10}f\bar{\psi}_{\bar{\alpha},\bar{\alpha}}\Bigr{)}\differential{\bar{a}}.$

Here $\bar{\Omega}=\\{(\bar{x}_{1},\bar{x}_{2})\,|\,(x_{1},x_{2})\in\Omega\\}$
denotes the rescaled 2-D domain, and components of vectors and their
derivatives having Greek indices with bar are related to the rescaled
coordinates $\bar{x}_{\alpha}$. We can further simplify this functional by
dividing it by $\mu h$. Then the minimization problem reduces to

$\bar{J}[u,\bar{\psi}_{\bar{\alpha}}]=\int_{\bar{\Omega}}\Bigl{\\{}\frac{1}{12}[\sigma(\bar{\psi}_{\bar{\alpha},\bar{\alpha}})^{2}+\bar{\psi}_{(\bar{\alpha},\bar{\beta})}\bar{\psi}_{(\bar{\alpha},\bar{\beta})}]+\frac{5}{12}(u_{,\bar{\alpha}}+\bar{\psi}_{\bar{\alpha}})(u_{,\bar{\alpha}}+\bar{\psi}_{\bar{\alpha}})\Bigr{\\}}\differential{\bar{a}}\\\
-\int_{\bar{\Omega}}\Bigl{(}\bar{f}u-\frac{\sigma}{10}\bar{f}\bar{\psi}_{\bar{\alpha},\bar{\alpha}}\Bigr{)}\differential{\bar{a}}\rightarrow\min_{u,\bar{\psi}_{\bar{\alpha}}},$
(1)

where

$\bar{f}=\frac{hf}{\mu}.$

Thus, $\bar{f}$ is equal to the plate thickness $h$ times the characteristic
strain $\varepsilon=f/\mu$, and since both the bending and shear stiffness in
the rescaled functional (1) have the order of unity, the minimizer must be of
the same order as $\bar{f}$. Returning to the original functions, we see that
the rotation angles $\psi_{\alpha}$ have the order of the characteristic
strain ($\varphi_{\alpha}$ are even much smaller [24, 25]), while the bending
measures $\rho_{\alpha\beta}$ have the order of the characteristic strain
divided by $h$, as discussed in the introduction. In this way, the rescaled
problem (1) provides an elegant and effective way to avoid the shear-locking
effect.

To solve problem (1) we need to pose the boundary conditions. If the part of
the plate’s edge, $\bar{\partial}_{k}$, is clamped, the admissible functions
must satisfy the following kinematical conditions

$u=0,\quad\bar{\psi}_{\bar{\alpha}}=0\quad\text{at $\bar{\partial}_{k}$}.$ (2)

If the remaining part of the plate’s edge is simply supported, then only the
kinematical condition $u=0$ should be fulfilled, while
$\bar{\psi}_{\bar{\alpha}}$ can have arbitrary variation at that part.
Finally, if the remaining part $\bar{\partial}_{s}$ of the plate’s edge is
free, then no constraints are imposed on $u$ and $\bar{\psi}_{\bar{\alpha}}$
there. After finding the solution, we want also to check the asymptotic
accuracy of the theory. For this purpose, we need to compute the true average
displacement of the plate by adding the correction term to the previously
determined deflection [14]

$\check{u}=u+\frac{h^{2}\sigma}{60}\rho_{\alpha\alpha}=u-\frac{\sigma}{60}\bar{\psi}_{\bar{\alpha},\bar{\alpha}}.$
(3)

The determination of $\varphi_{\alpha}$ also requires the first derivatives of
$u$. Therefore we would like to have the solution of (1) such that $u$ and
$\bar{\psi}_{\bar{\alpha}}$ are both $C^{1}$-functions. To check the
asymptotic accuracy we must compare functions $\check{u}$ and
$\check{\psi}_{\alpha}$ found by the 2-D FSDT with

$\langle
w_{3}(x_{\alpha},x_{3})\rangle\equiv\frac{1}{h}\int_{-h/2}^{h/2}w_{3}(x_{\alpha},x_{3})\differential{x_{3}}\quad\text{and}\quad\langle
w_{\alpha}(x_{\alpha},x_{3})x_{3}\rangle/(h^{3}/12),$ (4)

where $w_{i}(x_{\alpha},x_{3})$ are the displacements computed by the 3-D
exact theory of elasticity. If they agree with the accuracy up to
$h^{2}/l^{2}$, then the asymptotic accuracy of our FE-implementation is
guaranteed.

##  3 Finite element implementation

###  3.1 Weak and strong formulations

Since we will only deal with rescaled coordinates and rotation angles in this
Section and in the analytical parts of the next Section, we briefly omit all
bars in the functional (1). Calculating the first variation of (1) and
equating it to zero, we obtain the following necessary condition for the
minimum

$\variation{J}=\int_{\Omega}\Bigl{\\{}\frac{1}{6}[\sigma\psi_{\alpha,\alpha}\variation{\psi}_{\beta,\beta}+\psi_{(\alpha,\beta)}\variation{\psi}_{(\alpha,\beta)}]+\frac{5}{6}(u_{,\alpha}+\psi_{\alpha})(\variation{u}_{,\alpha}+\variation{\psi}_{\alpha})\Bigr{\\}}\differential{a}\\\
-\int_{\Omega}\Bigl{(}f\variation{u}-\frac{\sigma}{10}f\variation{\psi}_{\alpha,\alpha}\Bigr{)}\differential{a}=0.$
(5)

Introducing the notations

$\begin{split}&\variation{W}^{b}=\int_{\Omega}\frac{1}{6}[\sigma\psi_{\alpha,\alpha}\variation{\psi}_{\beta,\beta}+\psi_{(\alpha,\beta)}\variation{\psi}_{(\alpha,\beta)}]\differential{a},\\\
&\variation{W}^{s}=\int_{\Omega}\frac{5}{6}(u_{,\alpha}+\psi_{\alpha})(\variation{u}_{,\alpha}+\variation{\psi}_{\alpha})\differential{a},\\\
&\variation{W}^{e}=\int_{\Omega}(f\variation{u}-\frac{\sigma}{10}f\variation{\psi}_{\alpha,\alpha})\differential{a},\end{split}$
(6)

for the variations of the bending energy, shear energy, and the work of
external forces, respectively, we recast Eq. (5) in the form

$\variation{W}^{b}+\variation{W}^{s}=\variation{W}^{e}.$ (7)

Denote the space of admissible functions as
$\mathcal{K}=\\{(v,\chi_{\alpha})\,|\,(v,\chi_{\alpha})|_{\partial_{k}}=0\\}$.
The weak formulation of the problem is stated as follows: Given $f$, find
$(u,\psi_{\alpha})\in\mathcal{K}$ such that Eq. (7) is satisfied for all
$(\variation{u},\variation{\psi}_{\alpha})\in\mathcal{K}$. For the integrals
in (6) to be meaningful, $v$ and $\chi_{\alpha}$ must belong at least to the
Sobolev’s space of square integrable functions with the square integrable
first derivatives, $H^{1}(\Omega)$. However, since the desired asymptotic
accuracy of FSDT may require higher smoothness of $u$ and $\psi_{\alpha}$, the
continuity assumption in $\mathcal{K}$ still remains unspecified.

Although not relevant for the FE-implementation, the strong formulation of the
problem is also stated for the sake of completeness. Under the assumptions
that $u$ and $\psi_{\alpha}$ are doubly differentiable, we apply the partial
integration to Eq. (5) and put it in the form

$\variation{J}=\int_{\Omega}\Bigl{\\{}\Bigl{[}-\frac{\sigma}{6}\psi_{\beta,\beta\alpha}-\frac{1}{6}\psi_{(\alpha,\beta)\beta}\Bigr{]}\variation{\psi}_{\alpha}-\frac{5}{6}(u_{,\alpha}+\psi_{\alpha})_{,\alpha}\variation{u}+\frac{5}{6}(u_{,\alpha}+\psi_{\alpha})\variation{\psi}_{\alpha}\Bigr{\\}}\differential{a}\\\
-\int_{\Omega}\Bigl{(}f\variation{u}+\frac{\sigma}{10}f_{,\alpha}\variation{\psi}_{\alpha}\Bigr{)}\differential{a}+\int_{\partial_{s}}\Bigl{[}\Bigl{(}\frac{\sigma}{6}\psi_{\beta,\beta}n_{\alpha}+\frac{1}{6}\psi_{(\alpha,\beta)}n_{\beta}\Bigr{)}\variation{\psi}_{\alpha}\\\
+\frac{5}{6}(u_{,\alpha}+\psi_{\alpha})n_{\alpha}\variation{u}\Bigr{]}\differential{s}+\int_{\partial_{s}}\frac{\sigma}{10}fn_{\alpha}\variation{\psi}_{\alpha}\differential{s}=0,$
(8)

where $n_{\alpha}$ is the unit normal outward to the boundary and
$\differential s$ is the length element. It is assumed here that the remaining
part $\partial_{s}$ of the plate edge is free. Because of the arbitrariness of
the variations $\variation{\psi}_{\alpha}$ and $\variation{u}$ in $\Omega$ and
on $\partial_{s}$, from (8) follow the second order partial differential
equations

$\begin{split}&-\frac{\sigma+1}{6}\psi_{\beta,\beta\alpha}-\frac{1}{6}\psi_{\alpha,\beta\beta}+\frac{5}{6}(u_{,\alpha}+\psi_{\alpha})=\frac{\sigma}{10}f_{,\alpha},\\\
&-\frac{5}{6}(u_{,\alpha}+\psi_{\alpha})_{,\alpha}=f,\end{split}$ (9)

the kinematic boundary conditions (2) on $\partial_{k}$, and the natural
boundary conditions

$\begin{split}&\frac{\sigma}{6}\psi_{\beta,\beta}n_{\alpha}+\frac{1}{6}\psi_{(\alpha,\beta)}n_{\beta}+\frac{5}{6}(u_{,\alpha}+\psi_{\alpha})=-\frac{\sigma}{10}fn_{\alpha},\\\
&\frac{5}{6}(u_{,\alpha}+\psi_{\alpha})n_{\alpha}=0\end{split}$ (10)

on $\partial_{s}$. Eqs. (9), (2), and (10) constitute the strong formulation
of the problem.

###  3.2 Discretization

To discretize Eq. (7), interpolation spaces for displacement and rotation
angles are required. In the following, a general interpolation scheme is used
without regard to specific requirements on the order and continuity of the
function space as mentioned above. In this sense, $u$ and $\psi_{\alpha}$
($\alpha=1,2$) are considered as primary variables. Their interpolation is

$u=(\mathbf{N}^{u})^{T}\mathbf{u},\quad\psi_{\alpha}=(\mathbf{N}^{\psi})^{T}{\bf\it\psi}_{\alpha},$
(11)

where $\mathbf{N}^{u}$ and $\mathbf{N}^{\psi}$ are column vectors representing
the shape functions for $u$ and $\psi_{\alpha}$, respectively, while
$\mathbf{u}$ and ${\bf\it\psi}_{\alpha}$ are deflection and rotation angle
vectors, respectively. The vector notation is in column ordering unless
otherwise stated.

Based on Eq. (11) we find the variation of the primary variables

$\variation{u}=\variation{\mathbf{u}}^{T}\mathbf{N}^{u},\quad\variation{\psi}_{\alpha}=\variation{{\bf\it\psi}}^{T}_{\alpha}\mathbf{N}^{\psi},$
(12)

and their derivatives

$\begin{split}&\variation{u}_{,\alpha}=\variation{\mathbf{u}}^{T}\mathbf{N}^{u}_{,\alpha},\\\
&\variation{\psi}_{\alpha,\alpha}=\variation{{\bf\it\psi}}^{T}_{\alpha}\mathbf{N}^{\psi}_{,\alpha}=\variation{{\bf\it\psi}}^{T}_{1}\mathbf{N}^{\psi}_{,1}+\variation{{\bf\it\psi}}^{T}_{2}\mathbf{N}^{\psi}_{,2},\\\
&\variation{\psi}_{(\alpha,\beta)}=\frac{1}{2}(\variation{{\bf\it\psi}}^{T}_{\alpha}\mathbf{N}^{\psi}_{,\beta}+\variation{{\bf\it\psi}}^{T}_{\beta}\mathbf{N}^{\psi}_{,\alpha}).\end{split}$

Substituting Eqs. (12) into Eqs. (6), we obtain

$\begin{split}\variation{W}^{b}&=\frac{1}{6}\int_{\Omega}[\sigma(\variation{{\bf\it\psi}}^{T}_{1}\mathbf{N}^{\psi}_{,1}+\variation{{\bf\it\psi}}^{T}_{2}\mathbf{N}^{\psi}_{,2})(\psi_{1,1}+\psi_{2,2})+\variation{{\bf\it\psi}}^{T}_{1}\mathbf{N}^{\psi}_{,1}\psi_{1,1}\\\
&+\variation{{\bf\it\psi}}^{T}_{2}\mathbf{N}^{\psi}_{,2}\psi_{2,2}+(\variation{{\bf\it\psi}}^{T}_{1}\mathbf{N}^{\psi}_{,2}+\variation{{\bf\it\psi}}^{T}_{2}\mathbf{N}^{\psi}_{,1})\psi_{(1,2)}]\differential{a},\\\
\variation{W}^{s}&=\frac{5}{6}\int_{\Omega}[(\variation{\mathbf{u}}^{T}\mathbf{N}^{u}_{,1}+\variation{{\bf\it\psi}}^{T}_{1}\mathbf{N}^{\psi})(u_{,1}+\psi_{1})\\\
&+(\variation{\mathbf{u}}^{T}\mathbf{N}^{u}_{,2}+\variation{{\bf\it\psi}}^{T}_{2}\mathbf{N}^{\psi})(u_{,2}+\psi_{2})]\differential{a},\\\
\variation{W}^{e}&=\int_{\Omega}[\variation{\mathbf{u}}^{T}\mathbf{N}^{u}f-\frac{\sigma}{10}(\variation{{\bf\it\psi}}^{T}_{1}\mathbf{N}^{\psi}_{,1}+\variation{{\bf\it\psi}}^{T}_{2}\mathbf{N}^{\psi}_{,2})f]\differential{a}.\end{split}$

Consequently, the residual forces with respect to $u$ and $\psi_{\alpha}$ read

$\begin{split}\mathbf{R}^{u}&=\int_{\Omega}\mathbf{N}^{u}f\differential{a}-\frac{5}{6}\int_{\Omega}[\mathbf{N}^{u}_{,1}(u_{,1}+\psi_{1})+\mathbf{N}^{u}_{,2}(u_{,2}+\psi_{2})]\differential{a},\\\
\mathbf{R}^{\psi_{1}}&=-\int_{\Omega}\frac{\sigma}{10}\mathbf{N}^{\psi}_{,1}f\differential{a}-\frac{5}{6}\int_{\Omega}\mathbf{N}^{\psi}(u_{,1}+\psi_{1})\differential{a}\\\
&-\frac{1}{6}\int_{\Omega}[\sigma\mathbf{N}^{\psi}_{,1}(\psi_{1,1}+\psi_{2,2})+\mathbf{N}^{\psi}_{,1}\psi_{1,1}+\mathbf{N}^{\psi}_{,2}\psi_{(1,2)}]\differential{a},\\\
\mathbf{R}^{\psi_{2}}&=-\int_{\Omega}\frac{\sigma}{10}\mathbf{N}^{\psi}_{,2}f\differential{a}-\frac{5}{6}\int_{\Omega}\mathbf{N}^{\psi}(u_{,2}+\psi_{2})\differential{a}\\\
&-\frac{1}{6}\int_{\Omega}[\sigma\mathbf{N}^{\psi}_{,2}(\psi_{1,1}+\psi_{2,2})+\mathbf{N}^{\psi}_{,2}\psi_{2,2}+\mathbf{N}^{\psi}_{,1}\psi_{(1,2)}]\differential{a}.\end{split}$
(13)

The use of Eq. (11) in Eq. (13) allows us to derive the following blocks of
the stiffness matrix

$\begin{split}\mathbf{K}^{uu}&=\frac{5}{6}\int_{\Omega}\left[\mathbf{N}^{u}_{,1}(\mathbf{N}^{u}_{,1})^{T}+\mathbf{N}^{u}_{,2}(\mathbf{N}^{u}_{,2})^{T}\right]\differential{a},\\\
\mathbf{K}^{u\psi_{1}}&=\frac{5}{6}\int_{\Omega}\mathbf{N}^{\psi}_{,1}(\mathbf{N}^{\psi})^{T}\differential{a},\\\
\mathbf{K}^{u\psi_{2}}&=\frac{5}{6}\int_{\Omega}\mathbf{N}^{\psi}_{,2}(\mathbf{N}^{\psi})^{T}\differential{a},\end{split}$
$\begin{split}\mathbf{K}^{\psi_{1}u}&=\frac{5}{6}\int_{\Omega}\mathbf{N}^{\psi}(\mathbf{N}^{\psi}_{,1})^{T}\differential{a},\\\
\mathbf{K}^{\psi_{1}\psi_{1}}&=\int_{\Omega}\Bigl{[}\frac{5}{6}\mathbf{N}^{\psi}(\mathbf{N}^{\psi})^{T}+\frac{\sigma+1}{6}\mathbf{N}^{\psi}_{,1}(\mathbf{N}^{\psi}_{,1})^{T}+\frac{1}{12}\mathbf{N}^{\psi}_{,2}(\mathbf{N}^{\psi}_{,2})^{T}\Bigr{]}\differential{a},\\\
\mathbf{K}^{\psi_{1}\psi_{2}}&=\int_{\Omega}\Bigl{[}\frac{\sigma}{6}\mathbf{N}^{\psi}_{,1}(\mathbf{N}^{\psi}_{,2})^{T}+\frac{1}{12}\mathbf{N}^{\psi}_{,2}(\mathbf{N}^{\psi}_{,1})^{T}\Bigr{]}\differential{a},\end{split}$
$\begin{split}\mathbf{K}^{\psi_{2}u}&=\frac{5}{6}\int_{\Omega}\mathbf{N}^{\psi}(\mathbf{N}^{\psi}_{,2})^{T}\differential{a},\\\
\mathbf{K}^{\psi_{1}\psi_{2}}&=\int_{\Omega}\Bigl{[}\frac{\sigma}{6}\mathbf{N}^{\psi}_{,2}(\mathbf{N}^{\psi}_{,1})^{T}+\frac{1}{12}\mathbf{N}^{\psi}_{,1}(\mathbf{N}^{\psi}_{,2})^{T}\Bigr{]}\differential{a},\\\
\mathbf{K}^{\psi_{2}\psi_{2}}&=\int_{\Omega}\Bigl{[}\frac{5}{6}\mathbf{N}^{\psi}(\mathbf{N}^{\psi})^{T}+\frac{\sigma+1}{6}\mathbf{N}^{\psi}_{,2}(\mathbf{N}^{\psi}_{,2})^{T}+\frac{1}{12}\mathbf{N}^{\psi}_{,1}(\mathbf{N}^{\psi}_{,1})^{T}\Bigr{]}\differential{a}.\end{split}$

###  3.3 Isogeometric analysis

The weak form (8) contains the first order derivatives of deflection and
rotation angles and their variations and is therefore meaningful only if
$(u,\psi_{\alpha})$ and $(\variation{u},\variation{\psi}_{\alpha})$ belong to
the Sobolev space $H^{1}(\Omega)$. However, as mentioned in Section 2, for the
desired asymptotic accuracy, the discretization space must be at least $C^{1}$
to ensure continuity and smoothness of the solution. In this paper, we use
isogeometric analysis with the non-uniform rational B-splines (NURBS) shape
function [44]. This approach also allows to facilitate the refinement of the
mesh and the elevation of the discretization order if necessary.

In the context of geometric modeling, a plate can be represented by NURBS
surface patches, in which each patch can be described as a tensor product of
two univariate B-splines such as

$\mathbf{S}\left(\xi_{1},\xi_{2}\right)=\sum_{i=1}^{m}\sum_{j=1}^{n}N_{i}^{p}\left(\xi_{1}\right)N_{j}^{q}\left(\xi_{2}\right)\mathbf{P}_{ij}.$
(14)

In Eq. (14), the univariate B-spline basis function of order $p$ and $q$,
denoted as $N_{i}^{p}$ and $N_{j}^{q}$, is computed via recursive Cox-de-Boor
formula

$N_{i}^{p}(\xi)=\dfrac{\xi-\xi_{i}}{\xi_{i+p}-\xi_{i}}N^{p-1}_{i}(\xi)+\dfrac{\xi_{i+p+1}-\xi_{i}}{\xi_{i+p+1}-\xi_{i+1}}N^{p-1}_{i+1}(\xi),\;N_{i}^{0}(\xi)=\begin{cases}1&\xi_{i}\leq\xi\leq\xi_{i+1}\\\
0&\text{otherwise}\end{cases}.$

We use $\\{\mathbf{P}_{ij}\\}_{0\leq i\leq m,0\leq j\leq n}$ to denote the
control point grid. The definition of the control point in homogeneous
coordinates will facilitate the construction of NURBS patch. A surface can
comprise a single or multiple patches that are connected at the patch
interface. In the latter case, the concept is called multipatch analysis. In
the typical scenario, the patch information matches at the interface, i.e.
strong coupling, including the parametric (knot) information and the location
of the control points, the multipatch is naturally connected and the $C^{0}$
continuity is guaranteed on the patch interfaces. However, patch-wise
parametric matching is not the necessary condition. There are various methods
to maintain the patch continuity in a NURBS multipatch, such as penalty method
[45], Nitsche method [46] or bending strip method [47]. These methods impose
weak coupling condition and are more sophisticated to implement.

The design of the B-spline surface allows different interpolation orders in
each parametric direction. This feature is particularly useful for analyses of
materials with anisotropic behavior where accuracy needs to be improved in a
particular direction. In addition, rescaling of the mesh can be conveniently
performed by rescaling the control point coordinates. This eliminates the need
for remeshing.

The NURBS patches in a multipatch structure are macro elements where the shape
function needs the full knot vectors for values and derivatives evaluation.
This on one hand prevents the possibility of applying parallelization
algorithms to the numerical code, on the other hand it makes the element
matrices very dense, which reduces the performance of the sparse direct
solver. Bézier decomposition algorithm [48] is proposed to retain the local
characteristic of the finite element. In essence, the Bézier decomposition
strategy constructs a local operator for each knot span (in 1D) or knot cell
(in 2D and 3D) so that the shape function within the cell can be evaluated
independently with the knot vectors. The support domain, i.e. control points
influenced by the cell, of each cell is reduced to the minimum number along
with the number of shape functions required for each cell. This eliminates the
global access to patch-wise data and allows a typical finite element code to
work seamlessly with NURBS elements.

###  3.4 Structure of the finite element code

Figure 2: Structure of the computational code.

The computational code for simulating the proposed FSDT plate element is
implemented as an extension of the Kratos Multiphysics framework [49]. A brief
overview of the structure of the code is visualized in Fig. 2. The plate
element is implemented within the extension PlateAndShellApplication. By
combining with the existing extension IsogeometricApplication, which supports
NURBS patch modeling and Bézier extraction operator evaluation, a new
extension IsogeometricPlateAndShellApplication can be created. This module
enables isogeometric analysis with FSDT plate element. In addition, mesh
refinement and order elevation are supported via IsogeometricApplication.

##  4 Numerical examples

###  4.1 Clamped circular plate under uniform transverse load

As a first test problem, we apply the developed finite element code to the
numerical simulation of the axisymmetric bending of a circular plate clamped
at the edge and subjected to a constant transverse load $f$. We chose this 2-D
problem because it admits an analytical solution, so we can directly compare
the numerical solution with it to test the former’s divergence and investigate
the optimal choice of interpolation order. As the bending is axisymmetric, we
must have in the polar coordinates $(r,\theta)$ the vanishing rotation angle
$\psi_{\theta}=0$, while $u$ and $\psi_{r}$ are functions of $r$ only.
Consequently,

$\rho_{rr}=-\psi_{r,r},\quad\rho_{\theta\theta}=-\frac{1}{r}\psi_{r},$

while all other components of $\rho_{\alpha\beta}$ vanish. It is easy to show
that the problem (1) reduces to minimizing the following functional

$I[u,\psi_{r}]=\int_{0}^{R}\Bigl{\\{}\frac{1}{12}[\sigma(\psi_{r,r}+\frac{1}{r}\psi_{r})^{2}+(\psi_{r,r})^{2}+\frac{1}{r^{2}}(\psi_{r})^{2}]+\frac{5}{12}(u_{,r}+\psi_{r})^{2}\Bigr{\\}}r\differential{r}\\\
-\int_{0}^{R}\Bigl{[}fu-\frac{\sigma}{10}f\Bigl{(}\psi_{r,r}+\frac{1}{r}\psi_{r}\Bigr{)}\Bigr{]}r\differential{r}\rightarrow\min_{u(R)=\psi_{r}(R)=0},$

where $R$ is the dimensionless radius of the plate and all bars are still
dropped for short. The standard calculus of variations shows that the
minimizer must satisfy the equation

$\begin{split}&-\frac{\sigma}{6}\Bigl{[}(\psi_{r,r}+\frac{1}{r}\psi_{r})r\Bigr{]}_{,r}+\frac{\sigma}{6}(\psi_{r,r}+\frac{1}{r}\psi_{r})-\frac{1}{6}(\psi_{r,r}r)_{,r}\\\
&+\frac{1}{6r}\psi_{r}+\frac{5}{6}(u_{,r}+\psi_{r})r=0,\\\
&-\frac{5}{6}\left[(u_{,r}+\psi_{r})r\right]_{,r}=0,\end{split}$ (15)

subjected to the clamped boundary conditions $u(R)=\psi_{r}(R)=0$. Integrating
the second equation of (15) and using the fact that $\varphi=u_{,r}+\psi_{r}$
cannot be singular at $r=0$, we find that

$\varphi(r)=u_{,r}+\psi_{r}=-\frac{3}{5}fr.$ (16)

Substituting (16) into the first of (15), we reduce it to

$\psi_{r,rr}+\frac{1}{r}\psi_{r,r}-\frac{1}{r^{2}}\psi_{r}=-\frac{3}{\sigma+1}fr.$

Integrating this equation and using the boundary condition $\psi_{r}(R)=0$ as
well as the non-singularity of $\psi$ at $r=0$, we get

$\psi_{r}(r)=\frac{3f}{8(\sigma+1)}(R^{2}r-r^{3})=\frac{3(1-\nu)f}{8}(R^{2}r-r^{3}).$
(17)

Now, the deflection $u$ can be found from

$u_{,r}=-\frac{3}{5}fr-\psi_{r}=-\frac{3}{5}fr-\frac{3f}{8(\sigma+1)}(R^{2}r-r^{3}).$

Integrating this equation and using the boundary condition $u(R)=0$, we obtain
finally (cf. [50, 51])

$u(r)=\frac{3}{10}f(R^{2}-r^{2})+\frac{3f}{32(\sigma+1)}(R^{2}-r^{2})^{2}=\frac{3}{10}f(R^{2}-r^{2})+\frac{3(1-\nu)f}{32}(R^{2}-r^{2})^{2},$
(18)

Note that the deflection according to Kirchhoff’s plate theory reads [52]

$u_{K}(r)=\frac{3(1-\nu)f}{32}(R^{2}-r^{2})^{2}.$

while the angle of rotation, $-u_{K,r}$, coincides with (17).

Figure 3: Numerical example 4.1: Multipatch geometry (left) and contour plot
of scaled deflection $\left(\dfrac{u}{h\varepsilon\bar{R}^{4}}\right)$ at
$\bar{R}=10$ (right).

For the numerical simulations of the real plate let us return to the original
notation with bars. In the rescaled formulation, the problem for the circular
plate contains only three parameters: $\bar{R}=R/h$ (geometry), $\nu$
(material), and $\bar{f}=hf/\mu$ (load). The material parameter is chosen to
be $\nu=0.3$. Since the problem is linear, the solution depends linearly on
$\bar{f}$, so by setting $\bar{f}=1$ we calculate the normalized solution. To
get the real deflection, we just need to multiply our normalized $u$ by
$h\varepsilon$ (where $\varepsilon=f/\mu$). Similarly, to calculate the real
rotation angles, we need to multiply our normalized
$\bar{\psi}_{\bar{\alpha}}$ (or $\bar{\varphi}_{\bar{\alpha}}$) by
$\varepsilon$. To investigate the efficiency of the proposed plate element, we
perform the analysis with three dimensionless radii $\bar{R}=10$,
$\bar{R}=100$ and $\bar{R}=1000$. The circular plate is constructed with 5
NURBS fields, as shown in Fig. 3 (left).

Figure 4: Numerical example 4.1: Convergence rate of displacement ($L_{2}$
error) with $\bar{R}=10$ (left), $\bar{R}=100$ (middle) and $\bar{R}=1000$
(right).

The deflection agrees very well with the analytical solution (18) (with bars
being restored over the rescaled quantities), as shown in Table 1. The
numerical values of the normalized deflection (for $\bar{f}=1$) are large.
However, considering that the real deflection is obtained after multiplication
by $h\varepsilon$, where $\varepsilon$ is much smaller than $10^{-3}$, these
values become reasonable. The $L_{2}$ error shown in the right column of this
Table is computed according to

$\|e\|_{L_{2}}=\dfrac{\int_{\Omega}\left(u_{computed}-u_{ana}\right)^{2}\differential{a}}{\int_{\Omega}u^{2}_{ana}\differential{a}}.$

It can be seen that the cubic order discretization gives very good results
even at coarse mesh. Fig. 4 visualizes the convergence rate of analysis with
second order ($p=2,q=2$) and cubic order ($p=3,q=3$) discretization. We
observe that the error increases with increasing $\bar{R}$, yet the
convergence rate is close to optimal, showing that the proposed formulation
does not exhibit shear-locking behavior. Although rescaling inflates the
domain, the number of elements and degrees of freedom required to obtain a
small error is not significantly affected by this rescaling. A slight
reduction of convergence rate is observed for $\bar{R}=10$ of the cubic order
analysis. This can be attributed to the loss of accuracy when the $L_{2}$
error is very small ($\sim 10^{-9}$) and the machine precision starts to
affect the computation.

$\bar{R}=10$
---
| ndofs | $u$ |  $u_{e}$ % |  $L_{2}$ error
IGA-p2-q2 | 1752 | 6.8615046723e2 | 1.4503863949e-2 | 3.0660060196e-4
6492 | 6.8624375111e2 | 9.1058452571e-4 | 2.2378738174e-5
14232 | 6.8624876504e2 | 1.7995811726e-4 | 5.3214286086e-6
38712 | 6.8624983992e2 | 2.3326748357e-5 | 9.7389087961e-7
152412 | 6.8624999000e2 | 1.4579201519e-6 | 1.1120739899e-7
IGA-p3-q3 | 1986 | 6.8624841430e2 | 2.3106803895e-4 | 7.8773270723e-6
6936 | 6.8624990465e2 | 1.3894075265e-05 | 7.4005944177e-7
14886 | 6.8624998130e2 | 2.7245754934e-06 | 2.0315875429e-7
39786 | 6.8624999759e2 | 3.5171608593e-07 | 4.1964948985e-8
154536 | 6.8624999985e2 | 2.1741422667e-08 | 5.1396533350e-9
$\bar{R}=100$
| ndofs | $u$ | $u_{e}$ |  $L_{2}$ error
IGA-p2-q2 | 1752 | 6.4784590440e6 | 1.3257323276 | 2.4336817628e-2
6492 | 6.5593384992e6 | 9.3846634028e-02 | 1.7699590666e-3
14232 | 6.5642574481e6 | 1.8925472721e-02 | 3.5821555521e-4
38712 | 6.5653375422e6 | 2.4744165221e-03 | 4.6899163890e-5
152412 | 6.5654898199e6 | 1.5505456947e-04 | 2.9411995117e-6
IGA-p3-q3 | 1986 | 6.5644455940e6 | 1.6059797387e-2 | 2.9372761135e-4
6936 | 6.5654300922e6 | 1.0647752621e-03 | 1.8101151749e-5
14886 | 6.5654861376e6 | 2.1113940168e-04 | 3.5340802949e-6
39786 | 6.5654982029e6 | 2.7372178611e-05 | 4.5584241877e-7
154536 | 6.5654998877e6 | 1.7097563813e-06 | 2.8764590547e-8
$\bar{R}=1000$
| ndofs | $u$ | $u_{e}$ |  $L_{2}$ error
IGA-p2-q2 | 1752 | 4.6439563968e10 | 2.9235273640e1 | 4.4941755443e-1
6492 | 6.1501718808e10 | 6.2835235678e+00 | 1.0699473685e-1
14232 | 6.4595537398e10 | 1.5691548868e+00 | 2.8321455202e-2
38712 | 6.5473621571e10 | 2.3112797835e-01 | 4.3069485094e-3
152412 | 6.5615277512e10 | 1.5272292485e-02 | 2.8832105167e-4
IGA-p3-q3 | 1986 | 6.5493593646e10 | 2.0069447938e-1 | 3.4437565289e-3
6936 | 6.5582249763e10 | 6.5600061147e-02 | 1.1614533849e-3
14886 | 6.5613643904e10 | 1.7761589296e-02 | 3.1440539208e-4
39786 | 6.5623595741e10 | 2.5969542594e-03 | 4.4297675346e-5
154536 | 6.5625189222e10 | 1.6880430048e-04 | 2.8107143827e-6
Table 1: Numerical example 4.1: Deflection computed at the center of the plate
and error. $u_{e}$ is the error of deflection compared to the analytical
solution.

Due to the high accuracy of the solution, the distribution of the deflection
is almost radially symmetric as expected (see Fig. 3 (right) for the solution
of the cubic discretization with 39786 d.o.fs for $\bar{R}=10$). The
normalized deflection $u/(h\varepsilon\bar{R}^{4})$ together with the
normalized rotation angles, i.e. $\psi_{r}/(\varepsilon\bar{R}^{3})$ and
$\varphi_{r}/(\varepsilon\bar{R})$, for this particular discretization along
the radial direction are plotted in Fig. 5. It clearly shows perfect agreement
with the analytical formulas (16)-(18) (where the bars are recovered over the
scaled quantities) for both $\bar{R}=10$ and $\bar{R}=1000$. It is worth
noting that to evaluate $\varphi$, the derivatives of $u$ at the integration
point are calculated and then transferred to the control points using an
$L_{2}$-global projection algorithm [53].

Figure 5: Numerical example 4.1: Normalized deflection and rotation angle
along the radial direction for $\bar{R}=10$ (left) and $\bar{R}=1000$ (right).

###  4.2 Rectangular plate under uniform transverse load

For the second test problem, we apply the developed finite element code to the
numerical simulation of the bending of a rectangular plate subjected to a
constant load $f$. We assume that $\Omega$ occupies the domain
$(0,L)\times(-D/2,D/2)$ of the $(x_{1},x_{2})$-plane, the left edge at
$x_{1}=0$ is clamped while the three remaining edges of the plate are free
(see Fig. 6 (left)).

Figure 6: Numerical example 4.2: Geometry and boundary condition for (left)
the analysis using FSDT plate element and (right) the 2D analysis using solid
element.

When the depth of the rectangular plate, $D$, becomes large, the solution must
exhibit the plane strain state due to the almost translational invariance in
the $x_{2}$-direction. For FSDT in the rescaled formulation (1), this means
that the 2-D problem reduces to the following 1-D variational problem:
Minimize the energy functional of the beam-like model

$I[u,\psi]=\int_{0}^{L}\Bigl{[}\frac{1}{2}E_{b}(\psi_{,x})^{2}+\frac{1}{2}E_{s}(u_{,x}+\psi)^{2}\Bigr{]}\differential{x}\\\
-\int_{0}^{L}\Bigl{(}fu-\frac{\sigma}{10}f\psi_{,x}\Bigr{)}\differential{x}\rightarrow\min_{u(0)=\psi(0)=0},$
(19)

where the bars in this theoretical part are again dropped for short, $x\equiv
x_{1}$ and

$E_{b}=\frac{\sigma+1}{6}=\frac{1}{6(1-\nu)},\quad E_{s}=\frac{5}{6}.$

The vanishing first variation leads to the differential equations

$\begin{split}&-E_{b}\psi_{,xx}+E_{s}(u_{,x}+\psi)=0,\\\
&-E_{s}(u_{,x}+\psi)_{,x}=0,\end{split}$ (20)

which are subjected to the boundary conditions

$\begin{split}&u(0)=\psi(0)=0,\\\ &E_{b}\psi_{,x}(L)=-\frac{\sigma}{10}f,\quad
E_{s}(u_{,x}+\psi)|_{x=L}=0.\end{split}$

From the first equation of (20) follows

$u_{,x}+\psi=\varphi=\frac{E_{b}}{E_{s}}\psi_{,xx}.$

Substituting this into the second equation yield

$E_{b}\psi_{,xxx}=-f.$ (21)

Integrating Eq. (21) and using the boundary conditions, we obtain

$\psi(x)=\frac{f}{E_{b}}\Bigl{[}-\frac{1}{6}x^{3}+\frac{1}{2}Lx^{2}-(\frac{1}{2}L^{2}+\frac{\sigma}{10})x\Bigr{]}.$
(22)

The deflection can be found from the equation

$u_{,x}=\varphi-\psi=-\psi+\frac{E_{b}}{E_{s}}\psi_{,xx},$

which, together with the boundary condition $u(0)=0$, yields

$u(x)=\frac{f}{E_{b}}\Bigl{[}\frac{1}{24}x^{4}-\frac{1}{6}Lx^{3}+(\frac{1}{4}L^{2}+\frac{\sigma}{20})x^{2}\Bigr{]}+\frac{f}{E_{s}}(-\frac{x^{2}}{2}+Lx).$

However, this is still not the true average transverse displacement of the
plate. The latter should be computed in accordance with Eq. (3) giving

$\check{u}(x)=\frac{f}{E_{b}}\Bigl{[}\frac{1}{24}x^{4}-\frac{1}{6}Lx^{3}+(\frac{1}{4}L^{2}+\frac{\sigma}{20})x^{2}\Bigr{]}+\frac{f}{E_{s}}(-\frac{x^{2}}{2}+Lx)\\\
-\frac{\sigma}{60}\frac{f}{E_{b}}\Bigl{(}-\frac{1}{2}x^{2}+Lx-\frac{1}{2}L^{2}-\frac{\sigma}{10}\Bigr{)}.$
(23)

Note that the average transverse displacement according to the classical
Kirchhoff’s plate theory in this beam-like model is given by

$u_{K}(x)=\frac{f}{E_{b}}\Bigl{(}\frac{1}{24}x^{4}-\frac{1}{6}Lx^{3}+\frac{1}{4}L^{2}x^{2}\Bigr{)},$

while the angle of rotation reads

$\psi_{K}(x)=-u_{K,x}=-\frac{f}{E_{b}}\Bigl{(}\frac{1}{6}x^{3}-\frac{1}{2}Lx^{2}+\frac{1}{2}L^{2}x\Bigr{)}.$

The presence of the analytical solution given by Eqs. (22)-(23) for this
particular case allows a comparison with the solution of the 2-D FSDT and with
the solution of the exact elasticity theory. On one hand, the comparison with
the solution of the 2-D FSDT can show that our FE-code is free of shear-
locking and allows us to check its convergence and efficiency as $D$ becomes
large. On the other hand, the comparison with the solution of the exact
elasticity theory allows us to verify the asymptotic accuracy of our FE-code.
To solve this problem within the exact elasticity theory, we use the geometry
and boundary condition shown in Fig. 6 (right). The normal traction is applied
to the top and bottom of the plate, as shown in this Figure, to avoid its
elongation and at the same time maintain the bending combined with a pure
shear [14]. Under the plane strain condition the problem becomes 2-D, so we
can employ 2-D solid elements. Having the numerical solution within the plane
strain elasticity theory, the average transverse displacement and rotation
angles over the thickness are evaluated using Eq. (4). Within 2-D FSDT the
deflection and rotation angles are evaluated along the centerline $x_{2}=0$.
Then the true average transverse displacement $\check{u}$ is computed in
accordance with (3).

Both the plate analysis and the elastic solid analysis employ the cubic order
NURBS discretization with sufficiently fine mesh. Using again the bars for the
recalled quantities, we show the results of numerical simulation in Fig. 7, in
which two cases $\bar{L}=3$ (left) and $\bar{L}=10$ (right) are analyzed (with
$\bar{L}=L/h$).

Figure 7: Numerical example 4.2: Deflection and rotation angle for $\bar{L}=3$
(left) and $\bar{L}=10$ (right).

In Fig. 7, two trends are clearly visible. First, the numerical solution using
the 2-D FSDT plate element converges to the analytical solution of FSDT given
by Eqs. (22)-(23) (with the bars being recovered over the scaled quantities)
as $\bar{D}=D/h$ increases, as shown in all plots for three different cases
$\bar{D}=1$, $\bar{D}=10$, and $\bar{D}=100$. This is due to the fact that the
boundary effect becomes negligible for large $\bar{D}$ and agrees well with
the assumption of almost translational invariance discussed at the beginning
of this sub-Section. This also confirms that our rescaled formulation of the
FSDT is inherently shear-locking-free and that the FE-code is efficient.
Second, the solution of the 2-D FSDT approximates the solution of the exact
elasticity theory much better than that of the classical Kirchhoff plate
theory. For large $\bar{L}$ (thin plate), the difference between the solutions
of FSDT and exact elasticity theory is negligibly small. Even for moderate
$\bar{L}$ (moderately thick plate), the error is only noticeable near the end
point when the deflection reaches its maximum. This also confirms the
asymptotic accuracy of our FE-implementation of the FSDT. Since the difference
between the analytical solution of the FSDT (red solid line) and the solution
resulting from the classical Kirchhoff (beam-like) theory (blue dashed line)
becomes significant for moderate $\bar{L}$, the latter is not applicable for
moderately thick plates. A detail worth mentioning is that the true average
displacement calculated according to (3) does not exactly satisfy the
kinematic boundary condition: The last correction term in (3), although small,
does not vanish at $\bar{x}_{1}=0$ (cf. Eq. (23)). The reason for this is well
known: The FSDT correction might not work in a thin boundary layer near the
edge of the plate [14]. Far from the edge, this correction term of order
$h^{2}/l^{2}$ compared to unity is essential for ensuring the asymptotic
accuracy of the FSDT up to that order. Another detail worth mentioning is that
the elastic 2-D analysis does not use an approximation, but evaluates the
solution fields in a brute force method with a very fine mesh. The average
transverse displacements and rotation angles obtained by calculating integrals
(4) over the thickness in a post-processing represent only the integral
characteristics of the FSDT. The comparison between the detailed distributions
of displacements across the thickness is left for future investigations.

##  5 Conclusion

In this work, the inherently shear-locking-free formulation of the
asymptotically exact FSDT for homogeneous elastic plates in the rescaled
coordinates and rotation angles is found. This leads to an asymptotically
accurate and computationally efficient finite element code. The requirement
that the deflection and rotation angles belong to the $C^{1}$ function space
is only necessary to achieve asymptotic accuracy but is not relevant to the
shear-locking effect, which is absent in this formulation. Since our main
focus was on the inherently shear-locking-free formulation, the simplest
linear FSDT for homogeneous elastic plates was chosen intentionally. Once
understood, the results obtained in this paper can be extended in many
directions. Among others, we will develop FE-implementation based on the
inherently shear-locking-free formulation of the asymptotically exact: (i)
linear FSDT for homogeneous elastic shells, (ii) linear FSDT for FG-plates and
shells, (iii) dynamic FSDT for FG-plates and shells, (iv) nonlinear FSDT for
FG-plates and shells with application to buckling analysis. The asymptotically
accurate and shear-locking free finite element implementation of the linear
FSDT for homogeneous elastic shells will be addressed in our next publication.

## References

  * [1] E. Reissner, The effect of transverse shear deformation on the bending of elastic plates, J. Appl. Mech. 12 (2) (1945) 69–77.
  * [2] R. Mindlin, Influence of rotatory inertia and shear on flexural motions of isotropic, elastic plates, J. Appl. Mech. 18 (1) (1951) 31–38.
  * [3] K. J. Bathe, E. N. Dvorkin, A formulation of general shell elements – the use of mixed interpolation of tensorial components, J. Numer. Methods Eng. 22 (3) (1986) 697–722.
  * [4] D. N. Arnold, R. S. Falk, A uniformly accurate finite element method for the Reissner–Mindlin plate, SIAM J. Numer. Anal. 26 (6) (1989) 1276–1290.
  * [5] K. U. Bletzinger, M. Bischoff, E. Ramm, A unified approach for shear-locking-free triangular and rectangular shell finite elements, Comput. Struct. 75 (3) (2000) 321–334.
  * [6] H. Nguyen-Xuan, G. R. Liu, H. C. Thai, T. T. Nguyen, An edge-based smoothed finite element method (ES-FEM) with stabilized discrete shear gap technique for analysis of Reissner–Mindlin plates, Comput. Methods Appl. Mech. Eng. 199 (9-12) (2010) 471–489.
  * [7] C. H. Thai, H. Nguyen‐Xuan, N. Nguyen‐Thanh, T. H. Le, T. Nguyen‐Thoi, T. Rabczuk, Static, free vibration, and buckling analysis of laminated composite Reissner–Mindlin plates using NURBS‐based isogeometric approach, Int. J. Numer. Methods Eng. 91(6) (2012) 571-603.
  * [8] S. Yin, J. S. Hale, , T. Yu, T. Q. Bui, S. P. Bordas, Isogeometric locking-free plate element: a simple first order shear deformation theory for functionally graded plates, Compos. Struct. 118 (2014) 121-138.
  * [9] T. N. Nguyen, T. D. Ngo, H. Nguyen-Xuan, A novel three-variable shear deformation plate formulation: Theory and isogeometric implementation, Comput. Methods Appl. Mech. Eng. 326 (2017) 376–401.
  * [10] J. Videla, S. Natarajan, S. P. Bordas, A new locking-free polygonal plate element for thin and thick plates based on Reissner-Mindlin plate theory and assumed shear strain fields. Comput. Struct. 220 (2019) 32-42.
  * [11] C. M. Wang, J. N. Reddy, K. H. Lee (Eds.), Shear deformable beams and plates: Relationships with classical solutions, Elsevier, 2000.
  * [12] M. Batista, An elementary derivation of basic equations of the Reissner and Mindlin plate theories, Eng. Struct. 32 (3) (2010) 906–909.
  * [13] N. Challamel, I. Elishakoff, A brief history of first-order shear-deformable beam and plate models, Mech. Res. Commun. 102 (2019) 103389.
  * [14] V. L. Berdichevsky, Variational-asymptotic method of constructing a theory of shells, J. Appl. Math. Mech. 43 (4) (1979) 711–736.
  * [15] V. G. Sutyrin, Derivation of plate theory accounting asymptotically correct shear deformation, J. Appl. Mech. 64 (1997) 905–915.
  * [16] W. Yu, Mathematical construction of a Reissner–Mindlin plate theory for composite laminates, Int. J. Solids Struct. 42 (26) (2005) 6680–6699.
  * [17] K. C. Le, An asymptotically exact first-order shear deformation theory for functionally graded plates. Int. J. Eng. Sci. 190 (2023) 103875.
  * [18] K. C. Le, B. D. Nguyen, On bending of single crystal beam with continuously distributed dislocations, Int. J. Plasticity 48 (2013) 152–167.
  * [19] K. C. Le, L. T. K. Nguyen, Energy Methods in Dynamics, Springer, Berlin, 2014.
  * [20] K. C. Le, J.-H. Yi, An asymptotically exact theory of smart sandwich shells, Int. J. Eng. Sci. 106 (2016) 179–198.
  * [21] K. C. Le, An asymptotically exact theory of functionally graded piezoelectric shells, Int. J. Eng. Sci. 112 (2017) 42-62.
  * [22] K. C. Le, Introduction to Micromechanics (second edition). Nova Science Publishers, New York, 2020.
  * [23] K. C. Le, T. M. Tran, Asymptotically exact theory of fiber-reinforced composite beams, Compos. Struct. 244 (2020) 112279.
  * [24] K. C. Le, Vibrations of Shells and Rods, Springer Verlag, Berlin, 1999.
  * [25] V. L. Berdichevsky, Variational Principles of Continuum Mechanics, Vol. 2, Springer Verlag, Berlin, 2009.
  * [26] O. C. Zienkiewicz, R. L. Taylor, J. M. Too, Reduced and selective integration techniques in general analysis of plates and shells, Int. J. Numer. Methods Eng. 3 (1971) 275–290.
  * [27] T. J. R. Hughes, R. L. Taylor, W. Kano-Knukulchal, A simple and efficient finite element for plate bending, Int. J. Numer. Methods Eng. 11 (1977) 1529-1543.
  * [28] T. J. R. Hughes, M. Cohen, M. Haroun, Reduced and selective integration techniques in finite element method of plates, Nucl. Eng. Des. 46 (1978) 203-222.
  * [29] E. D. L. Pugh, E. Hinton, O. C. Zienkiewicz, A study of quadrilateral plate bending elements with reduced integration, Int. J. Numer. Methods Eng. 12 (1978) 1059-1079.
  * [30] D. S. Malkus, T. J. R. Hughes, Mixed finite element methods-reduced and selective integration techniques: a unification of concepts, Comput. Methods Appl. Mech. Eng. 15 (1978) 63-81.
  * [31] C. P. Jackson, Singular capacity matrices produced by low-order Gaussian interpolation in the finite element method, Int. J. Numer. Meth. Eng. 17 (1981) 871-877.
  * [32] L. J. Hayes, Practical stability test for finite elements with reduced integration, Int. J. Numer. Meth. Eng. 17 (1981) 1689-1695.
  * [33] L. J. Hayes, Ill-conditioning resulting from reduced integration, Int. J. Numer. Meth. Eng. 17 (1981) 1279-1280.
  * [34] T. J. R. Hughes, T. Tezduyar, Finite elements based upon Mindlin plate theory with particular reference to the four-node isoparametric element, J. Appl. Mech. 48 (1981) 587-596.
  * [35] M. A. Crisfield, A quadratic Mindlin element using shear constraints, Comput. Struct. 18 (1984) 833-852.
  * [36] S. W. Lee, T. H. H. Pian, Finite elements based upon Mindlin plate theory with particular reference to the four-node isoparametric element, AIAA J. (1978) 16, 29–34.
  * [37] S. W. Lee, C. Wong, Mixed formulation finite elements for Mindlin theory plate bending, Int. J. Numer. Methods Eng. 18 (1982) 1297-1311.
  * [38] O. C. Zienkiewicz, Z. Xu, L. F. Zeng, A. Samuelsson, N. E. Wiberg, Linked interpolation for Reissner–Mindlin plate element: Part I – a simple quadrilateral, Int. J. Numer. Methods Eng. 36 (1993) 3043-3056.
  * [39] J. Dolbow, T. Belytschko, Volumetric locking in the element free Galerkin method, Int. J. Numer. Methods Eng. 46(6) (1999) 925-942.
  * [40] J. C. Simo, M. S. Rifai, A class of mixed assumed strain methods and the method of incompatible modes, Int. J. Numer. Methods Eng. 29 (1990) 1595-1638.
  * [41] R. P. R. Cardoso, J. W. Yoon, M. Mahardika, S. Choudhry, R.J. Alves de Sousa, R. A. Fontes Valente, Enhanced assumed strain (EAS) and assumed natural strain (ANS) methods for one-point quadrature solid-shell elements, Int. J. Numer. Methods Eng. 75 (2008) 156-187.
  * [42] A. Tessler, T. J. R. Hughes, A three-node Mindlin plate element with improved transverse shear, Comput. Methods Appl. Mech. Eng. 50 (1985) 71–101.
  * [43] B. Oesterle, S. Bieber, R. Sachse, E. Ramm, M. Bischoff, Intrinsically locking‐free formulations for isogeometric beam, plate and shell analysis, PAMM 18(1) (2018) e201800399.
  * [44] T. J. R. Hughes, J. Cottrell, Y. Bazilevs, Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement, Comput. Methods Appl. Mech. Eng. 194 (2005) 4135–4195.
  * [45] B. Marussig, T. J. R. Hughes, A review of trimming in isogeometric analysis: challenges, data exchange and simulation aspect, Arch. Computat. Methods. Eng. 25 (2018) 1059-1127.
  * [46] M. Ruess, D. Schillinger, A. Özcan, E. Rank, Weak coupling for isogeometric analysis of non-matching and trimmed multi-patch geometries, Comput. Methods Appl. Mech. Eng. 269 (2014) 46-71.
  * [47] J. Kiendl, Y. Bazilevs, M. C. Hsu, R. Wüchner, K-U. Bletzinger, The bending strip method for isogeometric analysis of Kirchhoff–Love shell structures comprised of multiple patches, Comput. Methods Appl. Mech. Eng. 199 (2010) 2403-2416.
  * [48] M. J. Borden, M. A. Scott, J. A. Evans, T. J. R. Hughes, Isogeometric finite element data structures based on Bézier extraction of NURBS, Int. J. Numer. Methods Eng. 87 (2011) 15–47.
  * [49] P. Dadvand, R. Rossi, E. Onate, An object-oriented environment for developing finite element codes for multi-disciplinary applications, Arch. Computat. Methods. Eng. 17 (2010) 253-297.
  * [50] C. M. Wang, K. H. Lee, Deflection and stress-resultants of axisymmetric Mindlin plates in terms of corresponding Kirchhoff solutions, Int. J. Mech. Sci. 38(11) (1996) 1179-1185.
  * [51] J. N. Reddy, Theories and analyses of beams and axisymmetric circular plates, CRC Press, 2022.
  * [52] L. D. Landau, E. M. Lifshitz, Theory of elasticity, Elsevier, 1986.
  * [53] X. Jiao, M. T. Heath, Common-refinement-based data transfer between non-matching meshes in multiphysics simulations, Int. J. Numer. Methods Eng. 61 (2004) 2402-2427.