# An asymptotically exact first-order shear deformation theory for
functionally graded plates

K. C. Lea,b111Phone: +84 93 4152458, email<EMAIL_ADDRESS>aDivision
of Computational Mathematics and Engineering, Institute for Computational
Science, Ton Duc Thang University, Ho Chi Minh City, Vietnam
bFaculty of Civil Engineering, Ton Duc Thang University, Ho Chi Minh City,
Vietnam

###### Abstract

An asymptotically exact first-order shear deformation theory for functionally
graded elastic plates is derived using the variational-asymptotic method. As
an application, an analytical solution to the problem of wave propagation in a
sandwich plate is found in accordance with this refined theory. Comparison
between the dispersion curves obtained by 2-D plate theory and 3-D elasticity
theory reveals that the former is accurate up to the order of $h^{2}/l^{2}$,
where $h$ is the plate thickness and $l$ the wavelength.

###### keywords:

functionally graded, plate, variational-asymptotic method, wave propagation,
dispersion curves.

††journal: International Journal of Engineering Science

##  1 Introduction

Reissner (32) was the first to establish a refined plate theory accounting for
transverse shear that bears his name and is now called the first-order shear
deformation plate theory (FSDT).222One should also acknowledge an important
contribution to this topic by Mindlin (27), although his motivation for
including transverse shear and rotatory inertia was somewhat different: he
wanted to capture the high-frequency thickness vibration of the plates. He
derived his theory from Castigliano’s variational principle by making several
hypotheses about the distribution of stresses on the transverse coordinate and
by determining the energy in terms of stress resultants, bending moments, and
shear forces. His pioneering work led to a large number of subsequent
publications, mainly because of the applicability of his theory to moderately
thick plates and the development of numerical methods (3, 1, 10, 29, 30), but
also because of the logic behind his derivation, which can be applied to other
problems (39, 4, 11). There have been a number of attempts to justify or re-
derive Reissner’s plate theory based on Kantorovic-Krylov’s semi-discrete
method for the approximate solution of variational problems (see, for example,
the review by Wang et al. (39)). However, the real advance in understanding
the asymptotic nature of the Reissner derivation was made by Berdichevsky (5,
6, 7). Using the variational-asymptotic method developed by himself, he gave
the asymptotically exact derivation of the refined plate and shell theories
accounting for transverse shear, which reduces to Reissner’s theory for
plates. He has shown that the integral characteristics of Reissner’s theory
are asymptotically accurate, but the distribution of displacements over the
thickness are not.

Considering the wide application of functionally graded (FG) materials and
structures, there is a large body of research proposing FSDT and higher order
shear deformation theories (HSDT) for FG-plates (31, 28, 37, 34, 12, 13),
among others333The literature on this topic is huge due to the variety of the
2-D theories for FG-plates: single-layer, multi-layer, refined theories
including rotatory inertias and transverse shears et cetera. It is therefore
impossible to cite all references. For the overview the reader may consult
(38, 34) and the references therein.. All the references cited above use the
Kantorovic-Krylov’s semi-discrete method, based on some assumptions about the
kinematics of the deformation or a stress state through the thickness of the
plates, to proceed with dimension reduction, so none of the obtained 2-D
theories guarantees the asymptotic accuracy. The dynamic problem of wave
propagation in FG-waveguides has been analyzed within the 3-D elasticity
theory by Lefebvre et al. (26), Baron and Naili (2), Kuznetsov (16, 17). Since
this analysis is not limited to low frequencies or long waves, its advantage
is that all dispersion curves in the entire domain can be studied, yielding
limits for low frequencies, points with zero group velocity, and zones with
negative group velocity. The disadvantage of this 3-D approach is that the
boundary conditions at the side edge of the finite plate lead to an infinite
number of coupled equations for the amplitudes of modes of vibrations, which
cannot be satisfied exactly. As for the asymptotic approaches, unfortunately,
the author of this article is not currently aware of any asymptotically
accurate FSDT for FG-plates. Berdichevsky (7) has tried to construct the
refined theory for anisotropic and inhomogeneous shells, which takes into
account the transverse shear under a constraint for the elastic moduli, but
the constraint of the constant Poisson ratio is too restrictive, so almost all
real FG-plates do not satisfy it. For laminated plates, which are a special
case of FG-plates, an asymptotic analysis of the energy functional was
provided by Sutyrin (36) and Yu (40). However, since in the general case of
the laminated plates the dimension reduction does not lead to an FSDT, they
have tried to optimize the parameters so that a derived theory is as close as
possible to asymptotic correctness while being an FSDT. It is also worth
mentioning some recent applications of the variational-asymptotic method to
plates and shells and to homogenization in (22, 24, 20, 21, 23, 35, 33).

The objective of this paper is to construct the asymptotically exact FSDT for
FG-plates by means of the variational-asymptotic method. Our goal is to
construct the refined plate theory with the asymptotic accuracy up to the
order of $h^{2}/l^{2}$, where $h$ is the plate thickness and $l$ is the
wavelength. Since we consider only the inhomogeneity in the transverse
direction, we assume that the mass density and elastic moduli of the isotropic
elastic FG-plate vary in this direction such that their distributions are even
functions of the transverse coordinate. The dimension reduction is based on
the asymptotic analysis of the action functional containing a small parameter
$h/l$, which allows to obtain the distribution of displacements from the
solution of the thickness problem. The obtained average energy must also be
extrapolated to ensure the correct behavior at short waves. As a result of
this asymptotic analysis, the FSDT is obtained, where all the coefficients of
the theory are given in closed analytical form and can be calculated for any
FG-plate satisfying the above assumption. We apply this theory to the problem
of wave propagation in a sandwich plate and compare the dispersion curves
resulting from the proposed theory and the elasticity theory. We show that the
former is exact up to the order of $h^{2}/l^{2}$, which justifies our
asymptotic analysis.

The paper is organized as follows. After this brief introduction, the
variational formulation of the problem is given in Section 2. Section 3 is
devoted to the asymptotic analysis of the action functional. Section 4
presents the two-dimensional FSDT for FG-plates. Section 5 analyzes the wave
propagation in a sandwich plate and shows the comparison of the dispersion
curves. Finally, Section 6 concludes the paper.

##  2 Variational formulation for FG-plates

Let $\Omega$ be a two-dimensional domain in the $(x_{1},x_{2})$-plane bounded
by a smooth closed curve $\partial\Omega$. We consider a plate which in the
undeformed 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

We analyze the forced vibration of the plate made of a functionally graded
isotropic elastic material whose mass density and elastic moduli depend on
$x_{3}$. The action functional reads

$I[\mathbf{w}(\mathbf{x},t)]=\int_{t_{0}}^{t_{1}}\int_{\Omega}\int_{-h/2}^{h/2}[T(x_{3},\dot{\mathbf{w}})-W(x_{3},{\bf\it\varepsilon})]\differential{x_{3}}\differential{a}\differential{t}\\\
+\int_{t_{0}}^{t_{1}}\int_{\Omega}({\bf\it\tau}\dotproduct\mathbf{w}|_{x_{3}=h/2}+{\bf\it\tau}\dotproduct\mathbf{w}|_{x_{3}=-h/2})\differential{a}\differential{t},$
(1)

where $\differential a=\differential x_{1}\differential x_{2}$ denotes the
area element of the mid-plane, $\mathbf{w}(\mathbf{x},t)$ is the displacement
field, while ${\bf\it\tau}$ is the traction acting on the upper and lower
faces of the plate. The explicit dependence of the Lagrangian on the
transverse coordinate $x_{3}$ of this functionally graded material is
precisely indicated. The kinetic energy density equals

$T(x_{3},\dot{\mathbf{w}})=\frac{1}{2}\rho(x_{3})\dot{\mathbf{w}}\dotproduct\dot{\mathbf{w}}=\frac{1}{2}\rho(x_{3})(\dot{w}_{\alpha}\dot{w}_{\alpha}+\dot{w}^{2}),$

with dot denoting the time derivative. The stored energy density $W$ reads

$W(x_{3},{\bf\it\varepsilon})=\frac{1}{2}\lambda(x_{3})(\varepsilon_{ii})^{2}+\mu(x_{3})\varepsilon_{ij}\varepsilon_{ij},$

with the components of the strain tensor ${\bf\it\varepsilon}$ being given by

$\varepsilon_{ij}=\frac{1}{2}(w_{i,j}+w_{j,i})\equiv w_{(i,j)}.$

We additionally assume that $\rho(x_{3})$, $\lambda(x_{3})$, and $\mu(x_{3})$
are even functions of $x_{3}$. Here and in the following, we use Latin
indices, running from 1 to 3, to refer to the spatial co-ordinates and the
Greek indices, running from 1 to 2, to refer to the plane co-ordinates $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 the symmetrization operation, and summation over repeated indices is
understood.

Hamilton’s variational principle states that the true displacement
$\check{\mathbf{w}}(\mathbf{x},t)$ of a functionally graded plate change in
space and time in such a way that the action functional (1) becomes stationary
at $\check{\mathbf{w}}(\mathbf{x},t)$ among all continuously differentiable
functions $\mathbf{w}(\mathbf{x},t)$ satisfying the initial and end conditions
as well as the kinematic boundary conditions (see, e.g., (7)). The problem is
to replace the three-dimensional action functional by an approximate two-
dimensional action functional for a thin FG-plate, whose functions depend only
on the longitudinal co-ordinates $x_{1},x_{2}$ and time $t$. The possibility
of dimension reduction is related to the smallness of the ratio between the
thickness $h$ and the characteristic scale of change of the deformation state
in the longitudinal directions, $l$, (see (19) and Section 3). We assume that

$\frac{h}{l}\ll 1.$

Additionally, we assume that

$\frac{h}{c\tau}\ll 1,$ (2)

where $\tau$ is the characteristic scale of change of the function
$\check{\mathbf{w}}(\mathbf{x},t)$ in time (see (19)) and $c$ the minimal
velocity of plane waves in the elastic material under consideration. This
means that we consider only statics or low-frequency vibrations of the
functionally graded elastic plate. By using the variational-asymptotic method,
the two-dimensional action functional will be constructed below in which terms
up to the order $h^{2}/l^{2}$ are exact compared to unity (the refined plate
theory accounting for transverse shear, now called FSDT).

In order to fix the domain of the transverse co-ordinate in the passage to the
limit $h\to 0$, we introduce the dimensionless co-ordinate

$\zeta=\frac{x_{3}}{h},\quad\zeta\in[-1/2,1/2],$

and transform the action functional to

$I[\mathbf{w}(\mathbf{x},t)]=\int_{t_{0}}^{t_{1}}\int_{\Omega}\int_{-1/2}^{1/2}h[T(\zeta,\dot{\mathbf{w}})-W(\zeta,{\bf\it\varepsilon})]\differential{\zeta}\differential{a}\differential{t}\\\
+\int_{t_{0}}^{t_{1}}\int_{\Omega}(\tau_{i}w_{i}|_{\zeta=1/2}+\tau_{i}w_{i}|_{\zeta=-1/2})\differential{a}\differential{t}.$
(3)

Now $h$ enters the action functional explicitly through the components of the
strain tensor $\varepsilon_{ij}$

$\varepsilon_{\alpha\beta}=w_{(\alpha,\beta)},\quad 2\varepsilon_{\alpha
3}=\frac{1}{h}w_{\alpha|\zeta}+w_{,\alpha},\quad\varepsilon_{33}=\frac{1}{h}w_{|\zeta}.$

The vertical bar followed by $\zeta$ indicates the partial derivative with
respect to $\zeta$ and not with respect to $x_{\zeta}$.

##  3 Asymptotic analysis of the action functional

We restrict ourselves to the low frequency vibrations of the FG-plate for
which assumption (2) is valid. Based on this assumption we may neglect the
kinetic energy density in the variational-asymptotic procedure of finding the
asymptotic expansion of the displacement field.444For the high-frequency
vibrations of elastic shells and rods where the kinetic energy density should
be kept in the variational-asymptotic analysis see (8, 9, 18, 19, 14).
Likewise, the work of external traction ${\bf\it\tau}$ can also be neglected
during this procedure. The average kinetic energy density and the work of
external traction can then be calculated once the asymptotic formulas for the
displacements are found. Before applying the variational-asymptotic procedure
to functional (3) let us transform the stored energy density to another form
more convenient for the asymptotic analysis. We note that, among terms of
$W(\varepsilon_{ij})$, the derivatives $w_{\alpha|\zeta}/h$ and $w_{|\zeta}/h$
in $\varepsilon_{\alpha 3}$ and $\varepsilon_{33}$ are the main ones in the
asymptotic sense. Therefore it is convenient to single out the components
$\varepsilon_{\alpha 3}$ and $\varepsilon_{33}$ in the stored energy density.
To this end we represent the latter as the sum of two quadratic forms
$W_{\parallel}$ and $W_{\perp}$ corresponding to longitudinal and transverse
stored energy densities, respectively. These are defined by

$W_{\parallel}=\min_{\varepsilon_{\alpha 3},\varepsilon_{33}}W,\quad
W_{\perp}=W-W_{\parallel}.$ (4)

Simple calculations give

$\displaystyle W_{\parallel}$
$\displaystyle=\mu[\sigma(\varepsilon_{\alpha\alpha})^{2}+\varepsilon_{\alpha\beta}\varepsilon_{\alpha\beta}],$
$\displaystyle W_{\perp}$
$\displaystyle=\frac{1}{2}(\lambda+2\mu)(\epsilon_{33}+\sigma\varepsilon_{\alpha\alpha})^{2}+2\mu\varepsilon_{\alpha
3}\varepsilon_{\alpha 3},$

where $\sigma=\lambda/(\lambda+2\mu)=\nu/(1-\nu)$. For FG-plate $\lambda$,
$\mu$, and $\sigma$ are functions of $\zeta$. The $\zeta$-argument is omitted
briefly in these functions when a precise specification is not required.

We could start the variational-asymptotic procedure with the determination of
the set $\mathcal{N}$ according to its general scheme (19, 7). As a result, it
would turn out that, at the first step, the function ${\bf w}$ does not depend
on the transverse co-ordinate $\zeta$: ${\bf w}={\bf u}(x_{\alpha},t)$; at the
second step additional degrees of freedom $\varphi_{\alpha}(x_{\alpha},t)$
associated with transverse shear occur in ${\bf w}^{*}$; and at the next step
${\bf w}^{**}$ is completely determined through ${\bf u}$ and
$\varphi_{\alpha}$. Thus, the set $\mathcal{N}$ according to the variational-
asymptotic scheme consists of functions ${\bf u}(x_{\alpha},t)$ and
$\varphi_{\alpha}(x_{\alpha},t)$. We will pass over these long, but otherwise
standard, derivations and make the following change of unknown functions

$\begin{split}&w_{\alpha}(x_{\alpha},\zeta,t)=u_{\alpha}(x_{\alpha},t)+h\varphi_{\alpha}(x_{\alpha},t)\zeta-
hu_{,\alpha}\zeta+hy_{\alpha}(x_{\alpha},\zeta,t),\\\
&w(x_{\alpha},\zeta,t)=u(x_{\alpha},t)-h\mathcal{I}[\sigma]u_{\beta,\beta}+hy(x_{\alpha},\zeta,t).\end{split}$
(5)

In these formulas $u_{\alpha},u$ correspond to the mean displacements of the
plate, while $\psi_{\alpha}=\varphi_{\alpha}-u_{,\alpha}$ describe the mean
rotation angles555It is easy to see that $\varphi_{\alpha}$ are the rotation
angles due to the transverse shear. If $\varphi_{\alpha}=0$, the shear strains
$\varepsilon_{\alpha 3}$ vanish in the first approximation, and the transverse
fibers remain “on average” perpendicular to the deformed mid-plane as in the
Kirchhoff’s plate theory. in the following sense

$\begin{split}&u_{\alpha}(x_{\alpha},t)=\langle
w_{\alpha}(x_{\alpha},\zeta,t)\rangle,\quad u(x_{\alpha},t)=\langle
w(x_{\alpha},\zeta,t)\rangle,\\\ &\psi_{\alpha}(x_{\alpha},t)=\langle
w_{\alpha}(x_{\alpha},\zeta,t)\zeta\rangle/(h/12),\end{split}$ (6)

where $\langle.\rangle$ denotes the integration over $\zeta$ within the limits
$[-1/2,1/2]$. By $\mathcal{I}[f]$ we mean the following transformation applied
to an arbitrary function $f(\zeta)$

$\mathcal{I}[f](\zeta)=\int_{0}^{\zeta}f(\xi)\differential{\xi}-\langle\int_{0}^{\zeta}f(\xi)\differential{\xi}\rangle.$

By this definition, $\mathcal{I}[f](\zeta)$ is an odd function if $f(\zeta)$
is an even function and vice versa. Moreover,
$\langle\mathcal{I}[f]\rangle=0$.

Because of definitions (6) functions $y_{\alpha}$ and $y$ should satisfy the
following constraints

$\langle y_{\alpha}\rangle=0,\quad\langle
y_{\alpha}\zeta\rangle=0,\quad\langle y\rangle=0.$ (7)

Equations (5)-(7) set up a one-to-one correspondence between $w_{\alpha},w$
and the set of functions $u_{\alpha},u,\varphi_{\alpha},y_{\alpha},y$ and
determine the change in the unknown functions
$\\{w_{\alpha},w\\}\to\\{u_{\alpha},u,\varphi_{\alpha},y_{\alpha},y\\}$.

Asymptotic analysis enables one to determine the order of smallness of
$y_{\alpha},y$. If these terms are neglected, then the deformation state of a
plate is characterized solely by the measures of extension
$A_{\alpha\beta}=u_{(\alpha,\beta)}$ and the measures of bending
$B_{\alpha\beta}=u_{,\alpha\beta}-\varphi_{(\alpha,\beta)}$. We introduce the
following notation

$\displaystyle\varepsilon_{A}=\max_{\Omega}\sqrt{A_{\alpha\beta}A^{\alpha\beta}},\quad\varepsilon_{B}=h\max_{\Omega}\sqrt{B_{\alpha\beta}B^{\alpha\beta}},$
$\displaystyle\Delta_{\alpha}=\max_{\mathcal{V}}|y_{\alpha|\zeta}|,\quad\Delta=\max_{\mathcal{V}}|y_{|\zeta}|.$

Consider a certain point of the mid-plane $\Omega$. The best constant $l$ in
the inequalities

$\begin{split}\left|A_{\alpha\beta,\gamma}\right|\leq\frac{\varepsilon_{A}}{l},\quad
h\left|B_{\alpha\beta,\gamma}\right|\leq\frac{\varepsilon_{B}}{l},\quad\left|\varphi_{\alpha}\right|\leq\frac{h}{l}(\varepsilon_{A}+\varepsilon_{B})\\\
\left|\varphi_{\alpha,\beta}\right|\leq\frac{h}{l^{2}}(\varepsilon_{A}+\varepsilon_{B}),\quad\max_{\zeta}\left|y_{\alpha,\beta}\right|\leq\frac{\Delta_{\alpha}}{l},\quad\max_{\zeta}\left|y_{,\alpha}\right|\leq\frac{\Delta}{l}\end{split}$

is called the characteristic scale of change of the deformation state in the
longitudinal directions. We define the inner domain $\Omega_{i}$ as a
subdomain of $\Omega$ in which the following inequality holds:

$h/l\ll 1.$

We assume the domain $\Omega$ to consist of the inner domain $\Omega_{i}$ and
a domain $\Omega_{b}$ abutting on the contour $\partial\Omega$ with width of
the order $h$ (boundary layer). Then functional (3) can be decomposed into the
sum of two functionals, an inner one for which an iteration process will be
applied, and a boundary layer functional. We concentrate first on the inner
functional.

We now fix $u_{\alpha},u,\varphi_{\alpha}$ and seek $y_{\alpha},y$.
Substituting (5) into the action functional (3), we will keep in it the
asymptotically principal terms containing $y_{\alpha},y$ and neglect all other
terms. The estimations based on the above inequalities lead to the asymptotic
formulas

$\begin{split}&\varepsilon_{33}=-\sigma A_{\beta\beta}+y_{|\zeta},\quad
2\varepsilon_{\alpha
3}=\varphi_{\alpha}+y_{\alpha|\zeta}-h\mathcal{I}[\sigma]A_{\beta\beta,\alpha}+hy_{,\alpha},\\\
&\varepsilon_{\alpha\beta}=A_{\alpha\beta}-hB_{\alpha\beta}\zeta+hy_{(\alpha,\beta)}.\end{split}$
(8)

According to formulas (8) the longitudinal stored energy density does not
contain asymptotically principal terms containing $y_{\alpha},y$ and can be
neglected. Since the asymptotically principal terms in the transverse stored
energy density contains only the derivatives of $y_{\alpha},y$ with respect to
$\zeta$, we drop the integration over $\Omega$ and $t$ and reduce the
thickness problem to finding minimum of the functional

$I_{\perp}=\frac{h}{2}\langle[(\lambda+2\mu)(y_{|\zeta}-\sigma
hB_{\beta\beta}\zeta)^{2}\\\
+\mu(\varphi_{\alpha}+y_{\alpha|\zeta}-h\mathcal{I}[\sigma]A_{\beta\beta,\alpha}+hy_{,\alpha})(\varphi_{\alpha}+y_{\alpha|\zeta}-h\mathcal{I}[\sigma]A_{\gamma\gamma,\alpha}+hy_{,\alpha})]\rangle.$

This problem can be solved in two steps. We first minimize the functional

$\frac{h}{2}\langle(\lambda+2\mu)(y_{|\zeta}-\sigma
hB_{\beta\beta}\zeta)^{2}\rangle$

under constraint (7)3. Obviously, the minimum is equal to zero and is attained
at

$y=h\mathcal{I}[\sigma\zeta]B_{\beta\beta}.$ (9)

Now, substituting $y$ from here into the remaining term of $I_{\perp}$, we
arrive at the second step: Minimize

$\frac{h}{2}\langle\mu(\varphi_{\alpha}+y_{\alpha|\zeta}-h\mathcal{I}[\sigma]A_{\beta\beta,\alpha}+h^{2}\mathcal{I}[\sigma\zeta]B_{\beta\beta,\alpha})(\varphi_{\alpha}+y_{\alpha|\zeta}-h\mathcal{I}[\sigma]A_{\gamma\gamma,\alpha}+h^{2}\mathcal{I}[\sigma\zeta]B_{\gamma\gamma,\alpha})\rangle$

among $y_{\alpha}$ satisfying the first two constraints in (7). To simplify
this minimization problem we change the unknown functions

$y_{\alpha}=f(\zeta)hA_{\beta\beta,\alpha}+g(\zeta)h^{2}B_{\beta\beta,\alpha}+z_{\alpha}$
(10)

such that constraints (7)1,2 are fulfilled and

$y_{\alpha|\zeta}=h\mathcal{I}[\sigma]A_{\beta\beta,\alpha}-h^{2}\mathcal{I}[\sigma\zeta]B_{\beta\beta,\alpha}+\alpha
h^{2}B_{\beta\beta,\alpha}+z_{\alpha|\zeta},$

where $\alpha$ is still an unknown coefficient. Then

$f(\zeta)=\mathcal{I}[\mathcal{I}[\sigma]],\quad
g(\zeta)=-\mathcal{I}[\mathcal{I}[\sigma\zeta]]+\alpha\zeta.$

Since $\sigma(\zeta)$ is even function, $f(\zeta)$ is also even function and
$\langle f\rangle=0$. Therefore the first term on the right-hand side of (10)
satisfies the constraints (7)1,2. Function $g(\zeta)$ is odd function, so the
constraint (7)1 is automatically satisfied. To fulfill the second constraint
$\langle g(\zeta)\zeta\rangle=0$ the coefficient $\alpha$ must be

$\alpha=12\langle\zeta\mathcal{I}[\mathcal{I}[\sigma\zeta]]\rangle.$ (11)

With (10) we reduce the problem to minimizing the following functional

$J_{\perp}=\frac{h}{2}\langle\mu(z_{\alpha|\zeta}+\bar{\varphi}_{\alpha})(z_{\alpha|\zeta}+\bar{\varphi}_{\alpha})\rangle$

among $z_{\alpha}$ satisfying the constraints $\langle
z_{\alpha}\rangle=\langle z_{\alpha}\zeta\rangle=0$, where

$\bar{\varphi}_{\alpha}=\varphi_{\alpha}+\alpha h^{2}B_{\beta\beta,\alpha}.$
(12)

To solve this minimization problem we consider the dual variational problem
(7)

$\min J_{\perp}=\max J_{\perp}^{*},\quad J_{\perp}^{*}=h\langle
p_{\alpha}\bar{\varphi}_{\alpha}-\frac{1}{2\mu}p_{\alpha}p_{\alpha}\rangle.$

The maximum is sought among $p_{\alpha}$ and the numbers $a_{\alpha}$ and
$b_{\alpha}$ such that

$p_{\alpha|\zeta}=a_{\alpha}+b_{\alpha}\zeta,\quad p_{\alpha}(\pm 1/2)=0.$

The numbers $a_{\alpha}$ and $b_{\alpha}$ are the Lagrange multipliers for the
constraints $\langle z_{\alpha}\rangle=0$ and $\langle
z_{\alpha}\zeta\rangle=0$. Integrating the equations for $p_{\alpha}$ and
using the boundary conditions for them, we find that $a_{\alpha}=0$ and
$p_{\alpha}=b_{\alpha}k(\zeta)$, where $k(\zeta)=\frac{1}{2}(\zeta^{2}-1/4)$.
Substituting $p_{\alpha}$ into $J_{\perp}^{*}$ and maximizing with respect to
$b_{\alpha}$, we get

$\Phi_{\perp}=\min J_{\perp}=\max
J_{\perp}^{*}=\frac{1}{2}h\mu^{*}\bar{\varphi}_{\alpha}\bar{\varphi}_{\alpha},\quad\mu^{*}=\frac{\langle
k\rangle^{2}}{\langle k^{2}/\mu\rangle}.$

Since the $b_{\alpha}$ that give maximum to $J_{\perp}^{*}$ are

$b_{\alpha}=\frac{\langle k\rangle}{\langle
k^{2}/\mu\rangle}\bar{\varphi}_{\alpha},$

we have

$p_{\alpha}=b_{\alpha}k(\zeta)=\frac{\langle k\rangle k}{\langle
k^{2}/\mu\rangle}\bar{\varphi}_{\alpha}.$

Using the relation between minimizer and maximizer

$z_{\alpha|\zeta}+\bar{\varphi}_{\alpha}=\frac{1}{\mu}p_{\alpha},$

we find also $z_{\alpha}$ from here

$z_{\alpha}=\mathcal{I}\Bigl{[}\frac{\langle k\rangle k}{\langle
k^{2}/\mu\rangle\mu}-1\Bigr{]}\bar{\varphi}_{\alpha}.$ (13)

We want now to compute the average longitudinal energy density. Taking into
account (10), it is sufficient to approximate $\varepsilon_{\alpha\beta}$ in
this refined theory by

$\varepsilon_{\alpha\beta}=A_{\alpha\beta}-hB_{\alpha\beta}\zeta+\mathcal{I}[\mathcal{I}[\sigma]]h^{2}A_{\gamma\gamma,\alpha\beta}$

because other terms do not bring correction of order $h^{2}/l^{2}$.
Substituting $\varepsilon_{\alpha\beta}$ into $W_{\parallel}$ and integrating
over the thickness, we obtain

$\Phi_{\parallel}=h[\langle\mu\sigma\rangle(A_{\alpha\alpha})^{2}+\langle\mu\rangle
A_{\alpha\beta}A_{\alpha\beta}+\langle\mu\sigma\zeta^{2}\rangle
h^{2}(B_{\alpha\alpha})^{2}+\langle\mu\zeta^{2}\rangle
h^{2}B_{\alpha\beta}B_{\alpha\beta}\\\
+2\langle\mu\sigma\mathcal{I}[\mathcal{I}[\sigma]]\rangle
h^{2}A_{\delta\delta}A_{\gamma\gamma,\alpha\alpha}+2\langle\mu\mathcal{I}[\mathcal{I}[\sigma]]\rangle
h^{2}A_{\alpha\beta}A_{\gamma\gamma,\alpha\beta}].$

Note that, due to the evenness of elastic moduli, there is no cross term
between $A_{\alpha\beta}$ and $B_{\alpha\beta}$. However, in contrast to the
refined theory for homogeneous plates, there are two additional cross terms
between $A_{\alpha\beta}$ and $A_{\gamma\gamma,\alpha\beta}$ due to the
dependence of $\mu$ and $\sigma$ on $\zeta$. Adding the transverse and
longitudinal stored energy density together, we obtain the total average
stored energy density in the form

$\Phi=h[\langle\mu\sigma\rangle(A_{\alpha\alpha})^{2}+\langle\mu\rangle
A_{\alpha\beta}A_{\alpha\beta}+\langle\mu\sigma\zeta^{2}\rangle
h^{2}(B_{\alpha\alpha})^{2}+\langle\mu\zeta^{2}\rangle
h^{2}B_{\alpha\beta}B_{\alpha\beta}\\\
+2\langle\mu\sigma\mathcal{I}[\mathcal{I}[\sigma]]\rangle
h^{2}A_{\delta\delta}A_{\gamma\gamma,\alpha\alpha}+2\langle\mu\mathcal{I}[\mathcal{I}[\sigma]]\rangle
h^{2}A_{\alpha\beta}A_{\gamma\gamma,\alpha\beta}+\frac{1}{2}\mu^{*}\bar{\varphi}_{\alpha}\bar{\varphi}_{\alpha}].$
(14)

We turn next to the average kinetic energy density. To find the latter we must
first compute the velocity based on the asymptotic expansion (5). Within the
desired accuracy of the refined theory (7) we can neglect $\dot{y}_{\alpha}$
and $\dot{y}$ as compared to other time derivatives, so

$\dot{w}_{\alpha}=\dot{u}_{\alpha}+\dot{\psi}_{\alpha}\zeta,\quad\dot{w}=\dot{u}-h\mathcal{I}[\sigma]\dot{u}_{\alpha,\alpha}.$

Substituting these formulas into the kinetic energy density and integrating
over the thickness, we find that

$\Theta=\frac{h}{2}[\langle\rho\rangle\dot{u}_{\alpha}\dot{u}_{\alpha}+\langle\rho\rangle\dot{u}^{2}+\langle\rho\zeta^{2}\rangle
h^{2}\dot{\psi}_{\alpha}\dot{\psi}_{\alpha}+\langle\rho(\mathcal{I}[\sigma])^{2}\rangle
h^{2}(\dot{u}_{\alpha,\alpha})^{2}].$

Again, the evenness of $\rho$ and the oddness of $\zeta$ and
$\mathcal{I}[\sigma]$ ensure that no cross term occurs between
$\dot{u}_{,\alpha}$ and $\dot{\psi}_{\alpha}$ as well as between $\dot{u}$ and
$\dot{u}_{\alpha,\alpha}$. The average kinetic energy density turns out to be
simple enough.

Combining the average stored and kinetic energy densities, we obtain the
average 2-D action functional in the form

$J[u_{\alpha},u,\psi_{\alpha}]=\int_{t_{0}}^{t_{1}}\int_{\Omega}[\Theta(\dot{u}_{\alpha},\dot{u},\dot{\psi}_{\alpha})-\Phi(A_{\alpha\beta},B_{\alpha\beta},A_{\gamma\gamma,\alpha\beta},\bar{\varphi}_{\alpha})]\differential{a}\differential{t},$

with $\Theta$ and $\Phi$ given above. This action functional can still be
simplified. First, if we introduce a new unknown function

$\bar{u}=u+\alpha h^{2}B_{\alpha\alpha}.$ (15)

then the rotation angles $\psi_{\alpha}=\varphi_{\alpha}-u_{,\alpha}$ do not
alter their form after these changes of unknown functions (12) and (15),

$\psi_{\alpha}=\varphi_{\alpha}-u_{,\alpha}=\bar{\varphi}_{\alpha}-\bar{u}_{,\alpha},$

while $B_{\alpha\beta}$, expressed in terms of $\psi_{\alpha}$, do not contain
the second derivatives

$B_{\alpha\beta}=-\psi_{(\alpha,\beta)}.$

Therefore, it is convenient to choose $\psi_{\alpha}$ together with
$u_{\alpha}$ and $\bar{u}$ as the primary degrees of freedom and express the
total average stored energy density in their terms. However, by making the
change of unknown functions $u\to\bar{u}$ according to (15) $\Theta$ becomes
after neglecting terms of higher order of smallness

$\Theta=\frac{h}{2}[\langle\rho\rangle\dot{u}_{\alpha}\dot{u}_{\alpha}+\langle\rho\rangle\dot{\bar{u}}^{2}-2\langle\rho\rangle\alpha
h^{2}\dot{\bar{u}}\dot{\bar{u}}_{,\alpha\alpha}+\langle\rho\zeta^{2}\rangle
h^{2}\dot{\psi}_{\alpha}\dot{\psi}_{\alpha}+\langle\rho(\mathcal{I}[\sigma])^{2}\rangle
h^{2}(\dot{u}_{\alpha,\alpha})^{2}].$

This average kinetic energy density is asymptotically accurate (up to the
order $h^{2}/l^{2}$) for long waves but not suitable for short waves because
it is not positive definite. If we integrate the third term by parts,
neglecting the divergence term going to the boundary (as a null Lagrangian),
and adding small terms involving $\dot{\bar{\varphi}}_{\alpha}$, we get the
meaningful short wave extrapolation (5, 19) of the previous formula

$\Theta=\frac{h}{2}[\langle\rho\rangle\dot{u}_{\alpha}\dot{u}_{\alpha}+\langle\rho\rangle\dot{\bar{u}}^{2}+\langle\rho\zeta^{2}+2\rho\alpha\rangle
h^{2}\dot{\psi}_{\alpha}\dot{\psi}_{\alpha}+\langle\rho(\mathcal{I}[\sigma])^{2}\rangle
h^{2}(\dot{u}_{\alpha,\alpha})^{2}].$

Second, there is another short-wave extrapolation where the fifth and sixth
term in $\Phi_{\parallel}$, which are not positive definite and behave “badly”
at short waves, can be transformed and combined with the last term in
$\Theta$. Indeed, by integration by parts and neglecting the null-Lagrangian,
we reduce these two terms to

$2(\langle\mu\sigma\mathcal{I}[\mathcal{I}[\sigma]]\rangle+\langle\mu\mathcal{I}[\mathcal{I}[\sigma]]\rangle)h^{2}A_{\beta\beta}A_{\gamma\gamma,\alpha\alpha}.$
(16)

By observing that, in the first approximation (see next Section)

$\langle\rho\rangle\ddot{u}_{\alpha}=2\langle\mu\sigma\rangle
A_{\gamma\gamma,\alpha}+2\langle\mu\rangle A_{\alpha\beta,\beta},$

so that

$\langle\rho\rangle\ddot{u}_{\alpha,\alpha}=2(\langle\mu\sigma\rangle+\langle\mu\rangle)A_{\gamma\gamma,\alpha\alpha},$

we can replace the last factor $A_{\gamma\gamma,\alpha\alpha}$ of the
correction term in (16) by

$\frac{\langle\rho\rangle}{2(\langle\mu\sigma\rangle+\langle\mu\rangle)}\ddot{u}_{\alpha,\alpha}.$

Then, doing the integration by parts in time and neglecting the null-
Lagrangian, we can combine the last term in $\Theta$ with this term to
$(h/2)\rho^{*}h^{2}(\dot{u}_{\alpha,\alpha})^{2}$, where

$\rho^{*}=\langle\rho(\mathcal{I}[\sigma])^{2}\rangle-\frac{2\langle\rho\rangle(\langle\mu\sigma\mathcal{I}[\mathcal{I}[\sigma]]\rangle+\langle\mu\mathcal{I}[\mathcal{I}[\sigma]]\rangle)}{\langle\mu\sigma\rangle+\langle\mu\rangle}.$

Thus, the final formulas for the average stored and kinetic energy density
after this short-wave extrapolation are as follows

$\Phi=h[\langle\mu\sigma\rangle(A_{\alpha\alpha})^{2}+\langle\mu\rangle
A_{\alpha\beta}A_{\alpha\beta}+\langle\mu\sigma\zeta^{2}\rangle
h^{2}(B_{\alpha\alpha})^{2}+\langle\mu\zeta^{2}\rangle
h^{2}B_{\alpha\beta}B_{\alpha\beta}\\\
+\frac{1}{2}\mu^{*}(\psi_{\alpha}+\bar{u}_{,\alpha})(\psi_{\alpha}+\bar{u}_{,\alpha})],$
(17)

and

$\Theta=\frac{h}{2}[\langle\rho\rangle\dot{u}_{\alpha}\dot{u}_{\alpha}+\langle\rho\rangle\dot{\bar{u}}^{2}+\langle\rho\zeta^{2}+2\rho\alpha\rangle
h^{2}\dot{\psi}_{\alpha}\dot{\psi}_{\alpha}+\rho^{*}h^{2}(\dot{u}_{\alpha,\alpha})^{2}].$
(18)

We turn finally to the work of external tractions $\tau_{\alpha}$ and $\tau$.
Substituting $w_{\alpha}$ and $w$ from (5) into the last integral of (3) and
taking into account Eqs. (9), (10), (11), and (13), we find that

$\mathcal{A}=\int_{t_{0}}^{t_{1}}\int_{\Omega}\Bigl{\\{}P_{\alpha}u_{\alpha}+Pu+\frac{1}{2}hQ_{\alpha}\psi_{\alpha}-h\mathcal{I}[\sigma](1/2)QA_{\beta\beta}+h\mathcal{I}[\mathcal{I}[\sigma]](1/2)P_{\alpha}A_{\beta\beta,\alpha}\\\
+h\mathcal{I}[\sigma\zeta](1/2)PB_{\alpha\alpha}+h^{2}(\frac{1}{2}\alpha-\mathcal{I}[\mathcal{I}[\sigma\zeta]](1/2))Q_{\alpha}B_{\beta\beta,\alpha}\\\
+h\mathcal{I}\Bigl{[}\frac{\langle k\rangle k}{\langle
k^{2}/\mu\rangle\mu}-1\Bigr{]}(1/2)Q_{\alpha}\bar{\varphi}_{\alpha}\Bigr{\\}}\differential{a}\differential{t},$

where

$\begin{split}P_{\alpha}=\tau_{\alpha}|_{\zeta=1/2}+\tau_{\alpha}|_{\zeta=-1/2},\quad
P=\tau|_{\zeta=1/2}+\tau|_{\zeta=-1/2},\\\
Q_{\alpha}=\tau_{\alpha}|_{\zeta=1/2}-\tau_{\alpha}|_{\zeta=-1/2},\quad
Q=\tau|_{\zeta=1/2}-\tau|_{\zeta=-1/2}.\end{split}$

The fifth and seventh terms can be transformed by integration by parts. In
terms of the new unknown functions $\bar{u}$ and $\psi_{\alpha}$ the work
becomes finally

$\mathcal{A}=\int_{t_{0}}^{t_{1}}\int_{\Omega}\Bigl{\\{}P_{\alpha}u_{\alpha}+P\bar{u}+\frac{1}{2}hQ_{\alpha}\psi_{\alpha}-h\mathcal{I}[\sigma](1/2)QA_{\beta\beta}-h\mathcal{I}[\mathcal{I}[\sigma]](1/2)P_{\alpha,\alpha}A_{\beta\beta}\\\
+h\mathcal{I}[\sigma\zeta](1/2)PB_{\beta\beta}-h^{2}\alpha
PB_{\beta\beta}-h^{2}(\frac{1}{2}\alpha-\mathcal{I}[\mathcal{I}[\sigma\zeta]](1/2))Q_{\alpha,\alpha}B_{\beta\beta}\\\
+h\mathcal{I}\Bigl{[}\frac{\langle k\rangle k}{\langle
k^{2}/\mu\rangle\mu}-1\Bigr{]}(1/2)Q_{\alpha}(\psi_{\alpha}+\bar{u}_{,\alpha})\Bigr{\\}}\differential{a}\differential{t}.$
(19)

It is now appropriate to make a remark about the energy of the boundary layer.
Unfortunately, this energy is currently unknown, and the neglect of the null-
Lagrangian in the shortwave extrapolation makes its determination even more
difficult. As a result, the asymptotically exact boundary conditions for the
FSDT cannot be established. However, numerous numerical simulations for
homogeneous plates show that the natural boundary conditions give sufficiently
accurate results (19), so it is reasonable to expect the same for the FG-
plates.

##  4 Two-dimensional FSDT

We formulate now the variational principle of the 2-D FSDT of free vibration
of the FG-plate: the true displacements and rotation angles of the freely
vibrating FG-plate change in space and time in such a way that the 2-D average
action functional

$J[u_{\alpha},\bar{u},\psi_{\alpha}]=\int_{t_{0}}^{t_{1}}\int_{\Omega}[\Theta(\dot{u}_{\alpha},\dot{\bar{u}},\dot{\psi}_{\alpha})-\Phi(A_{\alpha\beta},B_{\alpha\beta},\psi_{\alpha}+\bar{u}_{,\alpha})]\differential{a}\differential{t},$

becomes stationary among all continuously differentiable functions
$u_{\alpha}(x_{\alpha},t)$, $\bar{u}(x_{\alpha},t)$, and
$\psi_{\alpha}(x_{\alpha},t)$. We assume that these functions are fixed at the
initial and end time. If the edge of the plate is free, then no constraints
are imposed on these functions at the boundary. The stored and kinetic energy
densities are taken from (17) and (18), respectively, with
$A_{\alpha\beta}=u_{(\alpha,\beta)}$ and
$B_{\alpha\beta}=-\psi_{(\alpha,\beta)}$. The standard calculus of variation
shows that the stationarity condition $\delta J=0$ implies the following two-
dimensional equations

$\begin{split}&\langle\rho\rangle
h\ddot{u}_{\alpha}-\rho^{*}h^{3}\ddot{u}_{\beta,\beta\alpha}=n_{\alpha\beta,\beta},\\\
&\langle\rho\rangle h\ddot{\bar{u}}=q_{\alpha,\alpha},\\\
&\langle\rho\zeta^{2}+2\rho\alpha\rangle
h^{3}\ddot{\psi}_{\alpha}=-m_{\alpha\beta,\beta}-q_{\alpha}.\end{split}$ (20)

The 2-D tensors $n_{\alpha\beta}$, $q_{\alpha}$, and $m_{\alpha\beta}$ are
given by the constitutive equations

$\begin{split}&n_{\alpha\beta}=\partialderivative{\Phi}{A_{\alpha\beta}}=2\langle\mu\sigma\rangle
hA_{\gamma\gamma}\delta_{\alpha\beta}+2\langle\mu\rangle hA_{\alpha\beta},\\\
&q_{\alpha}=\partialderivative{\Phi}{\bar{\varphi}_{\alpha}}=\mu^{*}h(\psi_{\alpha}+\bar{u}_{,\alpha}),\\\
&m_{\alpha\beta}=\partialderivative{\Phi}{B_{\alpha\beta}}=2\langle\mu\sigma\zeta^{2}\rangle
h^{3}B_{\gamma\gamma}\delta_{\alpha\beta}+2\langle\mu\zeta^{2}\rangle
h^{3}B_{\alpha\beta},\end{split}$ (21)

with $\delta_{\alpha\beta}$ being the Kronecker delta. For the plate with the
free edge, these equations are subjected to the natural boundary conditions

$n_{\alpha\beta}\nu_{\beta}=0,\quad q_{\alpha}\nu_{\alpha}=0,\quad
m_{\alpha\beta}\nu_{\beta}=0,$ (22)

with $\nu_{\alpha}$ being the components of the unit vector in the
$(x_{1},x_{2})$-plane normal to the curve $\partial\Omega$. For the plate with
the clamped or simply supported edge, the natural boundary conditions (22)
must be replaced by the corresponding kinematical boundary conditions. If
there are tractions acting on the faces of the plate so that we have forced
vibration, then the work $\mathcal{A}$ should be added to the action
functional. This leads to additional terms on the right-hand sides of Eqs.
(20) that can easily be obtained by varying the work (19) (see, e.g.,
Berdichevsky (7)). Note that the system of equations (20), (21) and boundary
conditions (22) reduces to those obtained by Reissner (32) and Berdichevsky
(5) for homogeneous plates with constant $\rho$, $\lambda$, and $\mu$.

To complete the 2-D FSDT of FG-plates we should also indicate the method of
restoring the 3-D stress and strain states by means of the 2-D fields. To do
this, the strain tensor field $\boldsymbol{\varepsilon}(x_{\alpha},\zeta,t)$
should be found from (8). Using the asymptotic formulas (9), (10), and (12),
we obtain

$\begin{split}&\varepsilon_{\alpha\beta}=A_{\alpha\beta}-h\zeta
B_{\alpha\beta},\quad 2\varepsilon_{\alpha 3}=\frac{\langle k\rangle
k}{\langle k^{2}/\mu\rangle\mu}\bar{\varphi}_{\alpha},\\\
&\varepsilon_{33}=-\sigma A_{\alpha\alpha}+\sigma\zeta
hB_{\alpha\alpha}.\end{split}$ (23)

In contrast to $\varepsilon_{\alpha\beta}$, the shear strains
$\varepsilon_{\alpha 3}$ and the normal strain $\varepsilon_{33}$ turn out to
be nonlinear functions of $\zeta$ due to the dependence of $\mu$ and $\sigma$
on $\zeta$. The stress tensor field $\boldsymbol{\sigma}(x_{\alpha},\zeta,t)$
is then determined by the 3-D constitutive equation (Hooke’s law). While doing
so, it is convenient to use the decomposition (4) for the stored energy
density. Within this approximation we find

$\begin{split}&\sigma_{\alpha\beta}=2\mu\sigma
A_{\gamma\gamma}\delta_{\alpha\beta}+2\mu A_{\alpha\beta}-2\mu\sigma\zeta
hB_{\gamma\gamma}\delta_{\alpha\beta}-2\mu\zeta hB_{\alpha\beta},\\\
&\sigma_{\alpha 3}=\frac{\langle k\rangle k}{\langle
k^{2}/\mu\rangle}\bar{\varphi}_{\alpha},\quad\sigma_{33}=0,\end{split}$ (24)

We see that, due to the dependence of $\mu$ and $\sigma$ on $\zeta$,
$\sigma_{\alpha\beta}$ and $\sigma_{\alpha 3}$ are not polynomials of $\zeta$,
as in the case of homogeneous plates. However, integrating these stresses over
the thickness of the plate, we get

$\langle\sigma_{\alpha\beta}\rangle=\frac{n_{\alpha\beta}}{h},\quad\langle\sigma_{\alpha\beta}\zeta\rangle=-\frac{m_{\alpha\beta}}{h^{2}},\quad\langle\sigma_{\alpha
3}\rangle=\frac{q_{\alpha}}{h}.$ (25)

Eqs. (25) give a clear physical interpretation of these integral
characteristics: $n_{\alpha\beta}$ are the stress resultants,
$m_{\alpha\beta}$ are the bending moments, and $q_{\alpha}$ are the shear
forces. For homogeneous plates the distributions (23) and (24) reduce to those
obtained by Berdichevsky (5), but differ from those by Reissner (32).
Concerning the integral characteristics (25) one can see the full agreement
with both cited works for homogeneous plates.

With the purpose of applying the theory to special FG-plates let us introduce
the short-hand notations

$\begin{split}&\rho_{1}=\langle\rho\rangle,\quad\rho_{2}=\rho^{*},\quad\rho_{3}=\langle\rho\zeta^{2}+2\rho\alpha\rangle,\quad\mu_{1}=\langle\mu\rangle,\\\
&\mu_{2}=\langle\mu\sigma\rangle,\quad\mu_{3}=\langle\mu\zeta^{2}\rangle,\quad\mu_{4}=\langle\mu\sigma\zeta^{2}\rangle,\quad\mu_{5}=\mu^{*}=\frac{\langle
k\rangle^{2}}{\langle k^{2}/\mu\rangle}.\end{split}$ (26)

Thus, together with $\alpha$ the total number of coefficients of the 2-D FSDT
is nine. Using these short-hand notations we present the system (20) and (21)
in the form

$\begin{split}&\rho_{1}h\ddot{u}_{\alpha}-\rho_{2}h^{3}\ddot{u}_{\beta,\beta\alpha}=n_{\alpha\beta,\beta}\\\
&\rho_{1}h\ddot{\bar{u}}=q_{\alpha,\alpha},\\\
&\rho_{3}h^{3}\ddot{\psi}_{\alpha}=-m_{\alpha\beta,\beta}-q_{\alpha},\end{split}$
(27)

and

$\begin{split}&n_{\alpha\beta}=2\mu_{2}hA_{\gamma\gamma}\delta_{\alpha\beta}+2\mu_{1}hA_{\alpha\beta}\\\
&q_{\alpha}=\mu_{5}h(\psi_{\alpha}+\bar{u}_{,\alpha}),\\\
&m_{\alpha\beta}=2\mu_{4}h^{3}B_{\gamma\gamma}\delta_{\alpha\beta}+2\mu_{3}h^{3}B_{\alpha\beta}.\end{split}$
(28)

##  5 Wave propagation in sandwich plates

Figure 2: Cross section of a sandwich plate

Since the sandwich plate is a special case of FG-plates that admits an
analytical solution to the 3-D problem of wave propagation (25) (see also
(15)), we can verify the asymptotic accuracy of our theory by applying it to
this test problem. Consider an infinite sandwich plate, shown in Fig. 2, whose
elastic moduli are given by

$\lambda(\zeta)=\begin{cases}\lambda_{c},&\text{$|\zeta|<b/2$},\\\
\lambda_{s},&\text{$b/2<|\zeta|<1/2$}\end{cases},\quad\mu(\zeta)=\begin{cases}\mu_{c},&\text{$|\zeta|<b/2$},\\\
\mu_{s},&\text{$b/2<|\zeta|<1/2$}.\end{cases}$

The mass density and the coefficient $\sigma=\lambda/(\lambda+2\mu)$ have
similar distributions

$\rho(\zeta)=\begin{cases}\rho_{c},&\text{$|\zeta|<b/2$},\\\
\rho_{s},&\text{$b/2<|\zeta|<1/2$}\end{cases},\quad\sigma(\zeta)=\begin{cases}\sigma_{c},&\text{$|\zeta|<b/2$},\\\
\sigma_{s},&\text{$b/2<|\zeta|<1/2$}.\end{cases}$

Here $b=h_{c}/h$, with $h_{c}$ being the height of the core layer. Based on
these distributions we can now compute all coefficients of the FSDT. Long, but
otherwise elementary calculations according to (26) give

$\begin{split}&\rho_{1}=b(\rho_{c}-\rho_{s})+\rho_{s},\\\
&\rho_{2}=\frac{1}{12\sigma_{s}}[b^{3}\sigma_{c}^{2}(\rho_{c}\sigma_{s}-\rho_{s}\sigma_{s})+\rho_{s}(b(\sigma_{c}-\sigma_{s})+\sigma_{s})^{3}]-\frac{\rho_{1}}{12(\mu_{1}+\mu_{2})}(1-b)b\\\
&\times[b^{2}(\sigma_{c}-\sigma_{s})-\sigma_{s}+b(2\sigma_{s}-3\sigma_{c})](\mu_{c}\sigma_{c}-\mu_{s}\sigma_{s}+\mu_{c}-\mu_{s}),\\\
&\rho_{3}=\frac{1}{12}[b^{3}(\rho_{c}-\rho_{s})+\rho_{s}]+2\alpha[b(\rho_{c}-\rho_{s})+\rho_{s}],\\\
&\mu_{1}=b(\mu_{c}-\mu_{s})+\mu_{s},\\\
&\mu_{2}=b(\mu_{c}\sigma_{c}-\mu_{s}\sigma_{s})+\mu_{s}\sigma_{s},\\\
&\mu_{3}=\frac{1}{12}[b^{3}(\mu_{c}-\mu_{s})+\mu_{s}],\\\
&\mu_{4}=\frac{1}{12}[b^{3}(\mu_{c}\sigma_{c}-\mu_{s}\sigma_{s})+\mu_{s}\sigma_{s}],\\\
&\mu_{5}=\frac{20\mu_{c}\mu_{s}}{3(1-b)^{3}(8+9b+3b^{2})\mu_{c}+3b(15-10b^{2}+3b^{4})\mu_{s}}\\\
&\alpha=\frac{1}{120}[-5b^{3}(\sigma_{c}-\sigma_{s})+3b^{5}(\sigma_{c}-\sigma_{s})-2\sigma_{c}].\end{split}$

If $b=1$ and the $\rho$, $\mu$, and $\sigma$ values of the core and skin
layers are the same (homogeneous plate), it is easy to verify that these
coefficients reduce to those in (5).

For sandwich plates (and generally for FG-plates with the even distributions
of $\rho$, $\mu$, and $\sigma$), the system of equations (25) and (26)
decomposes into the equations of the longitudinal and flexural waves. Let us
first consider the equations (25)1 and (26)1, which describe the longitudinal
waves. In terms of $u_{\alpha}$ the equation of motion becomes

$\rho_{1}h\ddot{u}_{\alpha}-\rho_{2}h^{3}\ddot{u}_{\beta,\beta\alpha}=h[(2\mu_{2}+\mu_{1})u_{\beta,\beta\alpha}+\mu_{1}\laplacian
u_{\alpha}],$ (29)

where $\laplacian$ is 2-D Laplace operator. We introduce the dimensionless
time and coordinates

$\tau=\frac{t}{h}\sqrt{\frac{\mu_{c}}{\rho_{c}}},\quad\zeta_{\alpha}=\frac{x_{\alpha}}{h}$
(30)

and rewrite Eq. (29) in the form

$r_{1}u_{\alpha|\tau\tau}-r_{2}u_{\beta|\beta\alpha\tau\tau}=(2m_{2}+m_{1})u_{\beta|\beta\alpha}+m_{1}\laplacian
u_{\alpha}.$ (31)

Here and in the sequel the vertical bar preceding indices denotes the
derivatives with respect to the corresponding dimensionless variables and, for
simplicity, we use the same $\nabla^{2}$ to denote the “dimensionless” Laplace
operator. The dimensionless coefficients are

$r_{i}=\frac{\rho_{i}}{\rho_{c}},\quad i=1,2,3;\quad
m_{i}=\frac{\mu_{i}}{\mu_{c}},\quad i=1,2,3,4,5.$

Let us seek the solutions of (31) in form of the harmonic plane waves
propagating in the direction $\zeta_{1}$,
$u_{\alpha}=a_{\alpha}e^{i(\kappa\zeta_{1}-\vartheta\tau)}$. It is easy to
show that there are two possible types of waves corresponding to

$\displaystyle u_{2}=a_{2}e^{i(\kappa\zeta_{1}-\vartheta\tau)},\quad
u_{1}=0\quad\text{SS-wave}$ $\displaystyle
u_{1}=a_{1}e^{i(\kappa\zeta_{1}-\vartheta\tau)},\quad
u_{2}=0\quad\text{L-wave}.$

For the SS-waves (symmetric shear waves) the dispersion relation reads

$r_{1}\vartheta^{2}=m_{1}\kappa^{2}$

that coincides with the exact dispersion relation of the theory of elasticity
(25). For L-waves we have from (31)

$r_{1}\vartheta^{2}+r_{2}\vartheta^{2}\kappa^{2}=2(m_{1}+m_{2})\kappa^{2}.$
(32)

For homogenous plates this equation becomes

$\vartheta^{2}+\frac{\sigma^{2}}{12}\vartheta^{2}\kappa^{2}=2(\sigma+1)\kappa^{2}.$

Bringing the second correction term on the left-hand side to the right-hand
side, and replacing $\vartheta^{2}$ in it by $2(\sigma+1)\kappa^{2}$, which is
allowed in this approximation, we obtain the asymptotic formula for the
dispersion relation of the refined plate theory

$\vartheta^{2}=2(\sigma+1)\kappa^{2}-\frac{\sigma^{2}(\sigma+1)}{6}\kappa^{4}=\frac{2}{1-\nu}\kappa^{2}-\frac{\nu^{2}}{6(1-\nu)^{3}}\kappa^{4}.$

This formula can also be derived from the Rayleigh-Lamb dispersion relation of
elasticity theory (19), which means that the FSDT in this particular case is
asymptotically accurate up to terms of order $h^{2}/l^{2}$. As for the
dispersion relation (32) for sandwich plates, we find that the first and last
terms in this equation agree with the asymptotic formula for the dispersion
relation from elasticity theory, which is the vanishing determinant of the
$6\times 6$ matrix ($V$-matrix) (25). It is also possible to check the
validity of the second correction term in equation (32) by expanding the
determinant after the Taylor series expansion of the sine and cosine functions
of the elements of this matrix. Since this is quite tedious and difficult to
find the main asymptotic terms in this huge expression, we decide to check
this by numerical simulation of the dispersion curves. Looking at Lee-Chang
dispersion relation and Eq. (32), we see that there are five parameters
determining each dispersion curve:

$b,\quad\sigma_{c}=\frac{\nu_{c}}{1-\nu_{c}},\quad\sigma_{s}=\frac{\nu_{s}}{1-\nu_{s}},\quad
r_{s}=\frac{\rho_{s}}{\rho_{c}},\quad m_{s}=\frac{\mu_{s}}{\mu_{c}}.$

Fig. 3 shows a representative result of the numerical simulations of the
dispersion curves, where the bold line corresponds to the dispersion curve
obtained from the elasticity theory and the dashed and dotted lines represent
the dispersion curves of the FSDT and the classical lamination theory for
sandwich plates, respectively. The parameters chosen for these numerical
simulations are: $b=0.8$, $\nu_{c}=0.3$, $\nu_{s}=0.35$, $r_{s}=0.8$,
$m_{s}=0.8$. As we can see, the FSDT agrees much better with the dispersion
curve of the elasticity theory than the classical lamination theory for
sandwich plates. The excellent match even extends to the shortwave range,
where $\kappa>1$. The variation of the parameters does not change this
feature.

Figure 3: Dispersion curve of longitudinal waves in sandwich plate: (i) Bold
line: theory of elasticity, (ii) Dashed line: FSDT for sandwich plate, (iii)
Dotted line: Classical lamination theory for sandwich plate. The chosen
parameters: $b=0.8$, $m_{s}=0.8$, $r_{s}=0.8$, $\nu_{c}=0.3$, $\nu_{s}=0.35$.

Let us now consider equations (27)2,3 and (28)2,3 describing the flexural
waves. We look for the harmonic plane waves propagating in $x_{1}$-direction
such that $\psi_{1}$ and $\bar{u}$ are functions of $x_{1}$ and $t$, while
$\psi_{2}=0$. Then the system of equations (27)2,3 and (28)2,3 reduces to

$\begin{split}&\rho_{1}h\ddot{\bar{u}}=\mu_{5}h(\psi_{1,1}+\bar{u}_{,11}),\\\
&\rho_{3}h^{3}\ddot{\psi}_{1}=2(\mu_{3}+\mu_{4})\psi_{1,11}-\mu_{5}h(\psi_{1}+\bar{u}_{,1}).\end{split}$
(33)

We introduce the dimensionless variables $\tau$ and $\zeta_{1}$ in accordance
with (30) and $\bar{\psi}=h\psi$ and rewrite Eqs. (33) in the form

$\begin{split}&r_{1}\bar{u}_{|\tau\tau}=m_{5}(\bar{\psi}_{1|1}+\bar{u}_{|11}),\\\
&r_{3}\bar{\psi}_{1|\tau\tau}=2(m_{3}+m_{4})\bar{\psi}_{1|11}-m_{5}(\bar{\psi}_{1}+\bar{u}_{|1}).\end{split}$
(34)

We seek solutions of these equations in form of the harmonic waves

$(\bar{u},\bar{\psi}_{1})=(a,b)e^{i(\kappa\zeta_{1}-\vartheta\tau)}$

and substitute it into (34). The condition of non-triviality of the solutions
yields the dispersion relation

$(m_{5}\kappa^{2}-r_{1}\vartheta^{2})[2(m_{3}+m_{4})\kappa^{2}-r_{3}\vartheta^{2}+m_{5}]-m_{5}^{2}\kappa^{2}=0.$
(35)

Expanding the left-hand side, we get

$2m_{5}(m_{3}+m_{4})\kappa^{4}-[m_{5}r_{3}+2(m_{3}+m_{4})r_{1}]\kappa^{2}\vartheta^{2}-m_{5}r_{1}\vartheta^{2}+r_{1}r_{3}\vartheta^{4}=0.$

Within the desired approximation the last term on the left-hand side of this
equation can be neglected, while in the second correction term $\vartheta^{2}$
can be replaced by $2(m_{3}+m_{4})\kappa^{4}/r_{1}$. This leads to the
following asymptotic formula

$\vartheta^{2}=\frac{2(m_{3}+m_{4})}{r_{1}}\kappa^{4}-2(m_{3}+m_{4})\Bigl{(}\frac{r_{3}}{r_{1}}+2\frac{m_{3}+m_{4}}{m_{5}}\Bigr{)}\kappa^{6}.$
(36)

For homogenous plates this asymptotic formula becomes

$\vartheta^{2}=\frac{1}{6(1-\nu)}\kappa^{4}-\frac{17-7\nu}{360(1-\nu)^{2}}\kappa^{6}.$

Note that this same formula can also be derived from the Rayleigh-Lamb
dispersion relation of elasticity theory (19), which means that the FSDT for
bending and flexural vibration of homogeneous plates is asymptotically exact
up to terms of order $h^{2}/l^{2}$.

Figure 4: Dispersion curve of flexural waves in sandwich plate: (i) Bold line:
theory of elasticity, (ii) Dashed line: FSDT for sandwich plate, (iii) Dotted
line: Classical lamination theory for sandwich plate. The chosen parameters:
$b=0.8$, $m_{s}=0.8$, $r_{s}=0.8$, $\nu_{c}=0.3$, $\nu_{s}=0.35$.

As for the dispersion relation (36) for sandwich plates, we note that the
first term on the right-hand side agree with the asymptotic formula for the
dispersion relation from elasticity theory (25). It is also possible to check
the validity of the second correction term in Eq. (36) by expanding the
determinant of the $6\times 6$ matrix ($W$-matrix) after the Taylor series
expansion of the trigonometric functions of the elements of this
matrix.666There is a misprint in (25) for $W$-matrix: in the element $w_{56}$
of this matrix $b_{1}$ should be replaced by $b_{2}$. Similar considerations
as in the previous case convinced us that the better way to check this is to
do numerical simulations of the dispersion curves. Fig. 4 shows a
representative result of these numerical simulations, where the bold line
corresponds to the dispersion curve obtained from the elasticity theory and
the dashed and dotted lines represent the dispersion curves of the FSDT in
accordance with (35) and the classical lamination theory for sandwich plates.
The parameters chosen for these numerical simulations remain the same as in
the previous case. As we can see, the FSDT agrees much better with the
dispersion curve of the elasticity theory than the classical lamination theory
for sandwich plates. The excellent agreement between the exact curve and that
of the FSDT is observed even in the short-wave range, up to $\kappa=3$. Again,
the variation of the parameters does not change this feature. This
verification confirms the asymptotic exactness of the FSDT up to order
$h^{2}/l^{2}$ for sandwich plates, and to our conviction also for FG-plates.
It is also interesting to note that, as observed by Kuznetsov (17), the
dispersion curves of guided waves in plates with even (asymmetric layout) or
odd (symmetric layout) number of alternating layers are significantly
different and remain different even with a large number of alternating layers.
Since such a plate can be considered as an FG-plate, the developed
variational-asymptotic method leads to asymptotic accuracy of FSDT only for
the symmetric layout, but not for the asymmetric one.

##  6 Conclusion

In this work, it has been demonstrated that the asymptotically exact FSDT of
FG-plates whose mass density and elastic moduli are even functions of the
transverse coordinate can be derived from the 3-D elasticity theory by the
variational-asymptotic method. Comparison of the dispersion curves obtained
from the elasticity theory and the FSDT for sandwich plates shows that the
latter is accurate up to an order of magnitude $h^{2}/l^{2}$ compared to
unity. The generalization to anisotropic elastic FG-shells, which brings
further terms describing the cross effect between extension and bending, is
straightforward. It is also not difficult to include the geometrically
nonlinear effect when the displacements become large. The open question is how
to derive the FSDT for FG-plates and shells whose mass density and elastic
moduli are not even functions of the transverse coordinate. As discussed at
the end of Section 4, the difference between symmetric and asymmetric layouts
limitates our asymptotic analysis, with the likely consequence that there is
no FSDT for FG-plates/shells with asymmetric layouts, and one must construct
the FSDT that most closely approximates the asymptotically exact theory of the
plate/shell (36, 40). Another extremely important question is the numerical
implementation of the refined plate theory that is free of shear-locking
effect. This issue will be addressed in our next paper.

## References

  * Arnold and Falk (1989) Arnold, D. N., Falk, R. S. (1989). A uniformly accurate finite element method for the Reissner–Mindlin plate. SIAM J. Numer. Anal. 26 (6), 1276–1290.
  * Baron and Naili (2010) Baron, C., Naili, S. (2010). Propagation of elastic waves in a fluid-loaded anisotropic functionally graded waveguide: Application to ultrasound characterization. J. Acoust. Soc. Am. 127 (3), 1307–1317.
  * Bathe and Dvorkin (1986) Bathe, K. J., Dvorkin, E. N. (1986). A formulation of general shell elements – the use of mixed interpolation of tensorial components. J. Numer. Methods Eng. 22 (3), 697–722.
  * Batista (2010) Batista, M. (2010). An elementary derivation of basic equations of the Reissner and Mindlin plate theories. Eng. Struct. 32 (3), 906–909.
  * Berdichevsky (1979a) Berdichevsky, V. L. (1979a). Variational-asymptotic method of constructing a theory of shells. J. Appl. Math. Mech. 43 (4), 711–736.
  * Berdichevsky (1979b) Berdichevsky, V. L. (1979b). Variational-asymptotic method of constructing the nonlinear shell theory. In G. K Mikhailov and W. T. Koiter (Eds.) Proceedings of IUTAM Symposium on Shell Theory (pp. 137–161). Amsterdam, North Holland.
  * Berdichevsky (2009) Berdichevsky, V. L. (2009). Variational Principles of Continuum Mechanics, Vol. 2. Springer Verlag, Berlin.
  * Berdichevsky and Le (1980) Berdichevsky, V. L., Le, K. C. (1980). High-frequency long-wave shell vibration. J. Appl. Math. Mech. 44 (4), 520–525.
  * Berdichevsky and Le (1982) Berdichevsky, V. L., Le, K. C. (1982) High-frequency vibrations of shells. Soviet Physics Doklady 27, 988-990.
  * Bletzinger et al. (2000) Bletzinger, K. U., Bischoff, M., Ramm, E. (2000). A unified approach for shear-locking-free triangular and rectangular shell finite elements. Comput. Struct. 75 (3), 321–334.
  * Challamel and Elishakoff (2019) Challamel, N., Elishakoff, I. (2019). A brief history of first-order shear-deformable beam and plate models. Mech. Res. Commun. 102, 103389.
  * Daikh et al. (2020) Daikh, A. A., Drai, A., Houari, M. S. A., Eltaher, M. A. (2020). Static analysis of multilayer nonlocal strain gradient nanobeam reinforced by carbon nanotubes. Steel Comp. Struct. 36 (6), 643-656.
  * Hirane et al. (2021) Hirane, H., Belarbi, M. O., Houari, M. S. A., Tounsi, A. (2021). On the layerwise finite element formulation for static and free vibration analysis of functionally graded sandwich plates. Eng. Comput. 38, 3871-3899.
  * Kaplunov et al. (1998) Kaplunov, J. D., Kossovitch, L. Y., Nolde, E. (1998). Dynamics of thin walled elastic bodies. Academic Press.
  * Kaplunov et al. (2017) Kaplunov, J., Prikazchikov, D. A., Prikazchikova, L. A. (2017). Dispersion of elastic waves in a strongly inhomogeneous three-layered plate. Int. J. Solids Struct. 113, 169–179.
  * Kuznetsov (2019) Kuznetsov, S. V. (2019). Abnormal dispersion of flexural Lamb waves in functionally graded plates. ZAMP 70 (3), 70–89.
  * Kuznetsov (2021) Kuznetsov, S. V. (2021). Guided waves in periodically stratified plates: Deficiency of homogenizations. Compos. Struct. 276, 114579.
  * Le (1997) Le, K. C. (1997). High frequency vibrations and wave propagation in elastic shells: variational-asymptotic approach. Int. J. Solids Struct. 34 (30), 3923–3939.
  * Le (1999) Le, K. C. (1999). Vibrations of shells and rods. Berlin: Springer Verlag.
  * Le (2017) Le, K. C. (2017). An asymptotically exact theory of functionally graded piezoelectric shells. Int. J. Eng. Sci. 112, 42-62.
  * Le (2020) Le, K. C. (2020). Introduction to Micromechanics (second edition). New York: Nova Science Publishers.
  * Le and Nguyen (2013) Le, K. C., Nguyen, B. D. (2013). On bending of single crystal beam with continuously distributed dislocations. Int. J. Plasticity 48, 152–167.
  * Le and Tran (2020) Le, K. C., Tran, T. M. (2020). Asymptotically exact theory of fiber-reinforced composite beams. Compos. Struct. 244, 112279.
  * Le and Yi (2016) Le, K. C., Yi, J.-H. (2016). An asymptotically exact theory of smart sandwich shells. Int. J. Eng. Sci. 106, 179–198.
  * Lee and Chang (1979) Lee, P. C. Y., Chang, N. (1979). Harmonic waves in elastic sandwich plates. J. Elast. 9 (1), 51–69.
  * Lefebvre et al. (2001) Lefebvre, J. E., Zhang, V., Gazalet, J., Gryba, T., Sadaune, V. (2001). Acoustic wave propagation in continuous functionally graded plates: an extension of the Legendre polynomial approach. IEEE Trans. Ultrason. Ferroelectr. Freq. Control 48 (5), 1332-1340.
  * Mindlin (1951) Mindlin, R. (1951). Influence of rotatory inertia and shear on flexural motions of isotropic, elastic plates. J. Appl. Mech. 18 (1), 31–38.
  * Nguyen et al. (2008) Nguyen, T. K., Sab, K., Bonnet, G. (2008). First-order shear deformation plate models for functionally graded materials. Compos. Struct. 83 (1), 25-36.
  * Nguyen et al. (2010) Nguyen-Xuan, H., Liu, G. R., Thai, H. C., Nguyen, T. T. (2010). 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), 471–489.
  * Nguyen et al. (2017) Nguyen, T. N., Ngo, T. D., Nguyen-Xuan, H. (2017). A novel three-variable shear deformation plate formulation: Theory and isogeometric implementation. Comput. Methods Appl. Mech. Eng. 326, 376–401.
  * Reddy (2000) Reddy, J. (2000). Analysis of functionally graded plates. J. Numer. Methods Eng. 47 (1‐3), 663–684.
  * Reissner (1945) Reissner, E. (1945). The effect of transverse shear deformation on the bending of elastic plates. J. Appl. Mech. 12 (2), 69–77.
  * Phanendra Kumar et al. (2022) Phanendra Kumar, A., Khajamoinuddin, S. M., Burela, R.G., Mahesh, V., Harursampath, D. (2022). Snap-through analysis of multistable laminate using the variational asymptotic method. Mech. Based Des. Struct. Mach. 1-26.
  * Shen (2016) Shen, H.-S. (2016). Functionally graded materials: nonlinear analysis of plates and shells. CRC press.
  * Shi et al. (2021) Shi, Z., Zhong, Y., Yi, Q., Peng, X. (2021). High efficiency analysis model for composite honeycomb sandwich plate by using variational asymptotic method. Thin-Walled Struct. 163, 107709.
  * Sutyrin (1997) Sutyrin, V. G. (1997). Derivation of plate theory accounting asymptotically correct shear deformation. J. Appl. Mech. 64, 905–915.
  * Thai and Choi (2013) Thai, H. T., Choi, D. H. (2013). A simple first-order shear deformation theory for the bending and free vibration analysis of functionally graded plates. Compos. Struct. 101, 332-340.
  * Thai and Kim (2015) Thai, H. T., Kim, S. E. (2015). A review of theories for the modeling and analysis of functionally graded plates and shells. Compos. Struct. 128, 70-86.
  * Wang et al. (2000) Wang, C. M., Reddy, J. N., Lee, K. H. (Eds.) (2000). Shear deformable beams and plates: Relationships with classical solutions. Elsevier.
  * Yu (2005) Yu, W. (2005). Mathematical construction of a Reissner–Mindlin plate theory for composite laminates. Int. J. Solids Struct. 42 (26), 6680–6699.