\makeatletter
\@ifundefined{HCode}
{\documentclass[screen,CRMECA,Unicode,published]{cedram}
\newenvironment{noXML}{}{}
\def\tsup#1{$^{{#1}}$}
\def\tsub#1{$_{{#1}}$}
\def\sfrac#1#2{{#1}/{#2}}
\def\thead{\noalign{\relax}\hline}
\def\endthead{\noalign{\relax}\hline}
\def\tbody{\noalign{\relax}\hline}
\usepackage{etoolbox}
\usepackage{upgreek}
\usepackage[scr=boondoxo, scrscaled=1]{mathalfa}
\csdef{Seqnsplit}{\\}
\newcounter{runlevel}
\let\MakeYrStrItalic\relax
\def\xxboldsymbol#1{\boldsymbol{#1}}
\def\jobid{crmeca20240330}
%\graphicspath{{/tmp/\jobid_figs/web/}} 
\graphicspath{{./figures/}} 
\def\hl{}
\def\xsection#1{}
\def\0{\phantom{0}}
\def\mn{\phantom{$-$}}
\def\refinput#1{}
\def\back#1{}
\def\botline{\\\hline}
\def\tabnote#1{\vskip3pt\parbox{33pc}{#1}}
\newenvironment{inftab}{}{}
\usepackage{multirow} 
\def\nrow#1{\@tempcnta #1\relax%
\advance\@tempcnta by 1\relax%
\xdef\lenrow{\the\@tempcnta}}
\def\morerows#1#2{\nrow{#1}\multirow{\lenrow}{*}{#2}}
\def\xmorerows#1#2{#2}
\DOI{10.5802/crmeca.274}
\datereceived{2024-04-26}
\daterevised{2024-06-19}
\datererevised{2024-09-28}
\dateaccepted{2024-11-12}
\ItHasTeXPublished
\def\bPhi{\boldsymbol{\Phi}}
\allowdisplaybreaks
}
{\documentclass[crmeca]{article}
\usepackage{upgreek}
\usepackage{stmaryrd}
\def\hl{}
\newtheorem{rema}{Remark}
\def\CDRdoi{10.5802/crmeca.274}
\def\newline{\unskip\break}
\def\xmorerows#1#2{\morerows{#1}{#2}}
\def\xxboldsymbol#1{#1}
\def\hyperlink#1#2{#2}
\def\hypertarget#1#2{#2}
\def\sfrac#1#2{{#1}/{#2}}
}
\makeatother

\dateposted{2025-01-09}
\begin{document}

\begin{noXML}

\CDRsetmeta{articletype}{research-article}

\title{Large deflections of agglomerated carbon nanotubes  reinforced
sandwich cantilever beam partially embedded on foundation}

\alttitle{Grandes d\'{e}flexions d'une poutre sandwich renforc\'{e}e
par des nanotubes de carbone agglom\'{e}r\'{e}s, partiellement
encastr\'{e}e sur une fondation}

\author{\firstname{Thi Thu Hoai} \lastname{Bui}\IsCorresp\CDRorcid{0000-0002-0258-4121}}
\address{Graduate University of Science and Technology, VAST, 18 Hoang
Quoc Viet, Hanoi, Vietnam}
\address{Faculty of Vehicle and Energy Engineering, Phenikaa
University, Ha Dong, Hanoi 12116, Vietnam}
\email[T. T. H. Bui]{hoai.buithithu@phenikaa-uni.edu.vn}

\author{\firstname{Thi Thu Huong} \lastname{Tran}\CDRorcid{0000-0001-7035-4283}}
\address{School of Mechanical Engineering, Hanoi University of Science
and Technology, Hanoi, Vietnam}

\author{\firstname{Vu Nam} \lastname{Pham}}
\address{Thuyloi University, 175 Son Tay, Dong Da, Hanoi, Vietnam}

\author{\firstname{Dinh~Kien~} \lastname{Nguyen}\CDRorcid{0000-0001-9356-8401}}
\address{Institute of Mechanics, VAST, 18 Hoang Quoc Viet, Hanoi, Vietnam}

\thanks{Vietnam National Foundation for Science and Technology
Development (NAFOSTED) (grant number 107.02-2021.11).}

\begin{abstract}
The large deflections of a carbon nanotubes (CNTs) reinforced sandwich
cantilever beam partially embedded on Pasternak foundation are studied
considering the influence of CNTs agglomeration. The sandwich beam is
composed of a homogeneous core and composite face layers with effective
moduli being estimated by Eshelby--Mori--Tanaka model. The nonlinear
equilibrium equation is constructed using a first-order shear
deformable nonlinear beam element and solved by the Newton--Raphson
based iterative procedure. The result reveals that the effect of
slenderness ratio on the large deflections is dependent on the degree
of CNT agglomeration, and this effect is more significant when the
agglomeration degree is more severe. The effects of the CNT volume
fraction, the degree of CNTs agglomeration and the foundation
parameters on the large deflections are studied in detail. The
influence of the sandwich configuration and the porosities on the
behavior of the sandwich beam is also examined and discussed.
\end{abstract}

\begin{altabstract}
Les grandes d\'{e}flexions d'une poutre sandwich renforc\'{e}e par des
nanotubes de carbone (CNTs) partiellement encastr\'{e}e sur une
fondation de Pasternak sont \'{e}tudi\'{e}es en tenant compte de
l'influence de l'agglom\'{e}ration des CNTs. La poutre sandwich est
compos\'{e}e d'un noyau homog\`{e}ne et de couches externes composites,
dont les modules effectifs sont estim\'{e}s \`{a} l'aide du mod\`{e}le
d'Eshelby--Mori--Tanaka.

L'\'{e}quation d'\'{e}quilibre non lin\'{e}aire est construite en
utilisant un \'{e}l\'{e}ment de poutre non lin\'{e}aire d\'{e}formable
en cisaillement de premier ordre, puis r\'{e}solue \`{a} l'aide d'une
proc\'{e}dure it\'{e}rative bas\'{e}e sur l'algorithme de
Newton--Raphson.

Les r\'{e}sultats r\'{e}v\`{e}lent que l'effet de l'\'{e}lancement sur
les grandes d\'{e}flexions d\'{e}pend du degr\'{e} d'agglom\'{e}ration
des CNTs, et cet effet devient plus important lorsque le degr\'{e}
d'agglom\'{e}ration est \'{e}lev\'{e}. Les influences de la fraction
volumique des CNTs, du degr\'{e} d'agglom\'{e}ration des CNTs et des
param\`{e}tres de la fondation sur les grandes d\'{e}flexions sont
\'{e}tudi\'{e}es en d\'{e}tail. L'impact de la configuration sandwich
et des porosit\'{e}s sur le comportement de la poutre sandwich est
\'{e}galement examin\'{e} et discut\'{e}.
\end{altabstract}

\shortrunauthors

\keywords{\kwd{Sandwich beam}
          \kwd{Carbon nanotubes}
          \kwd{Agglomeration}
          \kwd{Partial foundation embedment}
          \kwd{Large deflection}}

\altkeywords{\kwd{Poutre sandwich}
          \kwd{Nanotubes de carbone}
          \kwd{Agglom\'{e}ration}
          \kwd{Encastrement partiel dans une fondation}
          \kwd{Grande d\'{e}flexion}}

\maketitle

\end{noXML}

\section{Introduction}\label{sec1}

Carbon nanotubes (CNTs), with low density but large surface area and
high stiffness~\cite{1,2}, have made them a promising material for
reinforcing polymer matrix. It was shown that the properties of
polymer-based composites are significantly improved by using small
amounts of {CNTs} as reinforcement~\cite{3,4}. CNT-reinforced
composite (CNTRC), therefore is an ideal material for fabricating
high-performance structural elements. A large number of studies on the
mechanical behavior of CNTRC structures are summarized in the review
papers~\cite{5,6}, the investigations that are close to the
present topic are discussed below.

Shen et~al.~\cite{7,8} employed the two-step perturbation technique to
study the thermal effect of bending and buckling of CNTRC plates, in
which the CNTs are graded in the thickness direction by different
patterns. Shen and Xiang~\cite{9} proposed the multi-scale approach for
studying nonlinear vibration, bending and post-buckling of functionally
graded (FG) CNTRC beams embedded on foundation. Wattanasakulpong and
Ungbhakorn~\cite{10} employed various beam theories to study vibration,
bending and buckling of CNTRC beams, considering different types of CNT
distribution. The effect of CNT gradation on behavior of beams made
from FG-CNTRC was also considered  by Mayandi et~al.~\cite{11} using the
finite element method. The Eshelby--Mori--Tanaka model was adopted by Lei
et~al.~\cite{12} to evaluate the properties of composite in their large
deflection study of composite plates reinforced by CNTs. The large
deformation and post-buckling of cylindrical panels made from CNT
composite was studied in~\cite{13,14} using the \textit{kp}-Ritz
method. The moving least square Ritz method was employed
in~\cite{15,16,17} to study buckling and large deformation of CNTRC
plates. Zghal et~al.~\cite{18} developed a $C^{0}$ isoparametric shell
element for studying large deformation of CNTRC plates and panels. The
effect of CNT distributions on the nonlinear response of the CNTRC
structures was numerically examined by the authors.

FG-CNTRC was also applied to sandwich structures to improve their
performance. Bending and free vibration of sandwich plates with
FG-CNTRC face layers were considered by Natarajan et~al.~\cite{19} using
an 8-node $C^{0}$ plate element. Their finding showed that
the deflections are decreased by increasing the CNT volume, while the
frequencies are increased. Buckling of Timoshenko beams and spherical
caps made from sandwich CNT reinforced composite was considered in
references~\cite{20,21}, considering the linear distribution of CNTs.
The influence of the CNT volume fraction on the instability of CNTRC
sandwich panels due to periodic loads was also examined in
reference~\cite{22} using a Q8 plate element. Mehar et~al.~\cite{23}
used a higher-order finite element to assess stresses, deflections, and
frequencies of sandwich plates with FG-CNTRC faces under
thermo-mechanical loading. Thermal buckling of sandwich plates with
CNTRC face sheets was studied by {Kamarian}~\cite{24} on the basis of a
piecewise low-order shear deformation theory. Their finding showed that
the critical temperature of sandwich plates is considerably increased
by CNT reinforcement.

{In the above} discussed references, the CNTs' orientations were
assumed to be straight and aligned, and the bulk mechanical properties
of the composite are estimated by the rule of mixture. However, due to
the low bending stiffness and high aspect ratio, CNTs tend to
agglomerate in polymer matrices~\cite{25,26}. Shi et~al.~\cite{27}
employed Eshelby's inclusion model to develop a micromechanics
approach for assessing CNTs agglomeration effects. They demonstrated
that CNT agglomeration significantly reduces the stiffness of the
composite. Using the method in {reference}~\cite{27}, Yas and
{Heshmati}~\cite{28} studied the vibration of an FG-CNTRC beam
traversed by a moving load, considering a random       
orientation of CNTs. Their result reveals that the random
distribution of CNTs leads to higher beam deflections. The effects of
CNT agglomeration and waviness on vibration of CNTRC beams were
considered by Heshmati et~al.~\cite{29,30}. The DQM was employed by
Mehrabadi and Aragh~\cite{31} to compute the stresses of agglomerated
FG-CNTRC open cylindrical shells due to mechanical loading. Also using
the DQM, Tornabene et~al.~\cite{32} studied the distribution of
stresses, strains, and displacements of composite plates and shells
reinforced by agglomerated CNTs. The dynamic buckling of viscoelastic
microplates reinforced by agglomerated CNTs is considered by Kolahchi
and Cheraghba  \cite{33} using the Bolotin method in conjunction
with the Navier method. Safaei et~al.~\cite{34} investigated the thermal
and mechanical buckling of porous sandwich plates formed from an
isotropic core and CNTRC face layers. Galerkin method was used in
{references}~\cite{35, 36} to investigate the vibration of agglomerated
CNTRC plates and stability of nanocomposite beams, respectively. Kassa
et~al.~\cite{37} used the finite element method to study bending of
CNTRC tapered panels, considering the CNT waviness and agglomeration
effects. Based on the nonlocal elasticity theory, Daghigh
et~al.~\cite{38} derived the governing equations for bending and
buckling analyses of CNTRC composite nanoplates. The two-parameter
model in {reference}~\cite{27} was adopted to estimate the plate moduli
and Navier's solution was employed to derive expressions for
deflections and buckling loads. 

In this paper, the large deflections of a CNTRC sandwich cantilever
beam partially embedded on Pasternak foundation under end forces is
studied considering the influence of CNTs agglomeration. The sandwich
beam is composed of a homogeneous core and two composite face layers.
The face layers are made from a polymer-based CNTs reinforced composite
with effective elastic moduli being predicted by the
Eshelby--Mori--Tanaka model. The discretized equilibrium is derived using
a total Lagrangian nonlinear beam element, and solved by Newton--Raphson
iterative method in combination with the arc-length technique. The
effects of the CNT reinforcement, the agglomeration degree, the
foundation support, and the porosities on the large deflection behavior
of the sandwich cantilever beam are studied in detail.

\section{Problem statement}\label{sec2}

Figure~\ref{fig1} shows a CNTRC sandwich cantilever with a rectangular
cross-section $(b\times h)$, partially embedded on a Pasternak
foundation. The beam is formed from a homogeneous core and two CNTRC
layers. The figure, $x$-axis is in the mid-plane; ${L}$ and $L_{F}$ are
the lengths of the beam and foundation supporting part;
$z_{0}$, $z_{3}$, $z_{1}$, $z_{2}$ are coordinates in the
${z}$-direction of the bottommost and topmost surfaces, interfaces
between the layers, respectively.

\begin{figure}
\includegraphics{fig01}
\caption{\label{fig1}CNTRC sandwich cantilever beam partially embedded
in Pasternak foundation.}
\end{figure}

CNTs with high elastic moduli are used to reinforce polymer to enhance
the stiffness of composite layers. The polymer core is fully
characterized by elastic moduli of the polymer matrix, a homogenization
scheme, namely the Eshelby--Mori--Tanaka model~\cite{39,40}, is used
herein to predict elastic moduli of the composite faces. As mentioned
above, CNTs tend to agglomerate in the polymer, creating a sort of
spherical shaped inclusions, as depicted in Figure~\ref{fig2} for a
representative volume element (RVE). CNTs in the RVE are found both in
bundle or clusters inside the inclusions and scattered in the matrix.
This particular distribution has been described by Shi et~al.~\cite{27}
by a two-parameter model, which is briefly summarized below.

\begin{figure}
\includegraphics{fig02}
\caption{\label{fig2}RVE with agglomerated CNT inclusions.}
\end{figure}

The total CNTs inside the RVE, $(V_{\mathrm{CNT}})$, is split into two
parts as follows
{\begin{equation}\label{eq1}
V_{\mathrm{CNT}}={V}_{\mathrm{CNT}}^{\mathrm{in}}+{V}_{\mathrm{CNT}}^{\mathrm{out}}
\end{equation}}\unskip
with ${V}_{\mathrm{CNT}}^{\mathrm{in}}$ and
${V}_{\mathrm{CNT}}^{\mathrm{out}}$ are, respectively, the volumes of
CNTs inside and outside the inclusions.

The following parameters represent the CNT agglomeration 
{\begin{equation}\label{eq2}
\xi =\frac{V_{\mathrm{in}}}{V},\quad \zeta =\frac{{V}_{\mathrm{CNT}}^{\mathrm{in}}}
{V_{\mathrm{CNT}}}\quad \text{with } (\xi ,\,\zeta)\in [0,1]
\end{equation}}\unskip
In Equation~(\ref{eq2}), the parameter $\xi $ defines the volume
proportion of inclusions $(V_{\mathrm{in}})$ to the element volume
$(V)$, while parameter $\zeta $ quantifies the volume of CNTs inside
inclusions $({V}_{\mathrm{CNT}}^{\mathrm{in}})$ in respect of the total CNT
volume. For $\xi < 1$, the agglomeration is partial with nanofillers
are both in the inclusions and polymer matrix. The case $\xi =1$ means
that all CNTs are uniformly dispersed in the polymer, while $\zeta =1$
corresponds to the case that all nanotubes are in the clusters.

The effective bulk and shear moduli inside the inclusions,
$(K_{\mathrm{in}})$ and $(G_{\mathrm{in}})$, respectively, and that
outside inclusions, $(K_{\mathrm{out}})$ and $(G_{\mathrm{out}})$, can
be calculated as {follow}s~\cite{27}
{\begin{equation}\label{eq3}
K_{\mathrm{in}}=K_{m}+\frac{f_{\mathrm{CNT}}\zeta 
(\delta _{r}-3K_{m}\alpha _{r})}{3(\xi -f_{\mathrm{CNT}}
\zeta +f_{\mathrm{CNT}}\zeta \alpha _{r})},\quad 
G_{\mathrm{in}}=G_{m}+\frac{f_{\mathrm{CNT}}
\zeta (\eta _{r}-2G_{m}\beta _{r})}{2\left(\xi -
f_{\mathrm{CNT}}\zeta +f_{\mathrm{CNT}}\zeta \beta _{r}\right)}
\end{equation}}\unskip
and
{\begin{equation}\label{eq4}
\begin{array}{rcl}
\displaystyle K_{\mathrm{out}}&=&\displaystyle K_{m}+\frac{f_{\mathrm{CNT}}(1-\zeta )
(\delta _{r}-3K_{m}\alpha _{r})}{3\left[1-\xi 
-f_{\mathrm{CNT}}(1-\zeta)+f_{\mathrm{CNT}}
(1-\zeta )\alpha _{r}\right]},\vspace*{5pt}\\
\displaystyle G_{\mathrm{out}}&=&\displaystyle G_{m}+\frac{f_{\mathrm{CNT}}
(1-\zeta )(\eta _{r}-2G_{m}\beta _{r})}
{2\left[1-\xi -f_{\mathrm{CNT}}(1-\zeta )+f_{\mathrm{CNT}}
(1-\zeta )\beta _{r}\right]}
\end{array}
\end{equation}}\unskip
with $f_{\mathrm{CNT}}=V_{\mathrm{CNT}}/V$, and
{\begin{equation}\label{eq5}
\begin{array}{l}
\displaystyle\alpha _{r}=\frac{3(K_{m}+G_{m})+k_{r}-l_{r}}
{3(G_{m}+k_{r})},\quad  \delta _{r}=\frac{1}{3}
\left[n_{r}+2l_{r}\frac{(2k_{r}+l_{r})(3K_{m}+2G_{m}-l_{r})}
{G_{m}+k_{r}}\right]\vspace*{5pt}\\
\displaystyle\beta _{r}=\frac{1}{5}\left\{\frac{4G_{m}+
2k_{r}+l_{r}}{3(G_{m}+k_{r})}+\frac{4G_{m}}{G_{m}+
p_{r}}+\frac{2G_{m}[(3K_{m}+G_{m})+(3K_{m}+7G_{m})]}
{G_{m}(3K_{m}+G_{m})+m_{r}(3K_{m}+7G_{m})}\right\}\vspace*{5pt}\\
\displaystyle\eta _{r}=\frac{1}{5}\left[\frac{2}{3}
(n_{r}-l_{r})+\frac{8G_{m}p_{r}}{G_{m}+p_{r}}+
\frac{8m_{r}G_{m}(3K_{m}+4G_{m})}{3K_{m}
(m_{r}+G_{m})+G_{m}(7m_{r}+G_{m})}+
\frac{2(k_{r}-l_{r})(2G_{m}+l_{r})}{3(G_{m}+k_{r})}\right]
\end{array}
\end{equation}}\unskip
In the above equation, the {subscript} ``$m$'' and ``$r$'' are used to
indicate the matrix and the reinforcement; $k_r$, $m_r$, $n_r$, $l_r$,
and $p_r$ are Hill's elastic moduli of CNT reinforcement phase. It is
noted that Equation~(\ref{eq3}) is not mathematically valid for 
${\xi}= {\zeta}=0$, while Equation~(\ref{eq4}) is not valid for ${\xi}=
{\zeta} =1$. 

The bulk modulus and shear modulus of the composite are obtained by
Mori--Tanaka homogenization model as
{\begin{equation}\label{eq6}
K=K_{\mathrm{out}}\left[1+\frac{\xi \left(\frac{K_{\mathrm{in}}}
{K_{\mathrm{out}}}-1\right)}{1+\alpha (1-\xi)
\left(\frac{K_{\mathrm{in}}}{K_{\mathrm{out}}}-1\right)}
\right],\quad G=G_{\mathrm{out}}\left[1+
\frac{\xi \left(\frac{G_{\mathrm{in}}}{G_{\mathrm{out}}}-1\right)}
{1+\beta (1-\xi)\left(\frac{G_{\mathrm{in}}}{G_{\mathrm{out}}}-1\right)}
\right]
\end{equation}}\unskip
where $\alpha $ and $\beta $ are defines as
{\begin{equation}\label{eq7}
\alpha =\frac{1+\nu _{\mathrm{out}}}{3(1-\nu _{\mathrm{out}})},
\quad \beta =\frac{2(4-5\nu _{\mathrm{out}})}{15(1-\nu _{\mathrm{out}})}
\end{equation}}\unskip
with $\nu _{\mathrm{out}}=(3K_{\mathrm{out}}-
2G_{\mathrm{out}})/2(3K_{\mathrm{out}}+G_{\mathrm{out}})$. The
effective Young's modulus $(E)$ and Poisson's ratio $(\nu )$ of the
composite are calculated as
{\begin{equation}\label{eq8}
E=\frac{9KG}{3K+G},\quad \nu =\frac{3K-2G}{6K+2G}
\end{equation}}\unskip
The effective moduli estimated by the above equations do not consider
the influence of porosities. The porosities have a negative effect on
the elastic properties of the composite. For instance, an effective
property $({\mathcal{P}}_{\mathrm{POR}})$ of the CNTRC
material {c}onsidering the porosity effect is given by~\cite{34} 
{\begin{equation}\label{eq9}
{\mathcal{P}}_{\mathrm{POR}}=(1-f_{\mathrm{p}})
{\mathcal{P}},\quad f_{\mathrm{p}}=1-\left(\frac{2}{\uppi}
\sqrt{1-e_{0}}-\frac{2}{\uppi }+1\right)^{2}
\end{equation}}\unskip
with $f_{\mathrm{p}}$ is the porosity parameter, $e_0$ is the porosity
c{oefficient, and} ${\mathcal{P}}$ is the effective property of the
composite without porosities.

\section{Solution method}\label{sec3}

A total Lagrangian beam element, initially proposed by Pacoste and
Eriksson~\cite{41} and employed by Almeida et~al.~\cite{42} as well
as Nguyen et~al.~\cite{43} to model geometric nonlinearity of
functionally graded beams, is extended herewith to the CNTRC sandwich
beam to construct the discretized equilibrium equation and to compute
the response of the sandwich beam.

\subsection{Total Lagrangian beam element}\label{sec31}

A two-node beam element, $(\mathrm{I},\,\mathrm{J})$, initially
straight as shown in Figure~\ref{fig3} is considered. The element is
embedded on a Pasternak foundation with stiffness $k_{\mathrm{W}}$ 
and $k_{G}$ of
the Winkler and shear layers, respectively. The vector of nodal
displacements and rotation has the form
{\begin{equation}\label{eq10}
\mathbf{d}={\{\begin{array}{@{}cccccc@{}}
u_{\mathrm{I}}&w_{\mathrm{I}}&\theta _{\mathrm{I}}&
u_{\mathrm{J}}&w_{\mathrm{J}}&\theta _{\mathrm{J}}\end{array}\}}^{T}
\end{equation}}\unskip
in which $u_{\mathrm{I}}$, $w_{\mathrm{I}}$, $\theta _{\mathrm{I}}$
respectively denote the displacements in the ${x}$ and ${z}$ directions
and rotation at the node $\mathrm{I}$, and the same for node
$\mathrm{J}$.

\begin{figure}
\includegraphics{fig03}
\caption{\label{fig3}Initial and deformed configurations of the beam
element on Pasternak foundation.}
\end{figure}

The current deformed configuration of the element, ($\mathrm{I}'$,
$\mathrm{J}'$), is defined through the position vector $\mathbf{r}$ of
point M (see Figure~\ref{fig3}), which is defined by
{\begin{equation}\label{eq11}
\mathbf{r}=[x+u(x)]\mathbf{i}+w(x)\mathbf{j}
\end{equation}}\unskip
In the above equation, \textbf{i} and \textbf{j} {denote} unit
vectors of the $x$- and $z$-{axes,} w{hile} $x$ is the abscissa of the
point M, measured from the node $\mathrm{I}$ in the undeformed
configuration; $u(x)$ and $w(x)$ are the displacements in ${x}$- and
${z}$-directions, respectively. The cross-section S associated with
point M in the undeformed configuration undergoes large deflection and
rotation to section S$'$ related to point M$'$ in the deformed
configuration as depicted in Figure~\ref{fig3}. The tangent vector
$\mathbf{r}'$ of the deformed line $\mathrm{I}'\mathrm{J}'$ at point
M$'$ is expressed in terms of the normal strain $\varepsilon$ and the
shear strain $\gamma$ as 
{\begin{equation}\label{eq12}
\mathbf{r}'=\frac{\partial \mathbf{r}}{\partial x}=
\left[1+\varepsilon \right]\mathbf{e}+\gamma \mathbf{k}
\end{equation}}\unskip
In the above equation, \textbf{e} and \textbf{k} are {the} unit
vectors, orthogonal and parallel to the current section S', and they
are related to the base vectors \textbf{i} and \textbf{j} as
{\begin{equation}\label{eq13}
\mathbf{e}=\cos \theta \mathbf{i}+\sin \theta \,\mathbf{j},\quad 
\mathbf{k}=-\sin \theta \,\mathbf{i}+\cos \theta \,\mathbf{j}
\end{equation}}\unskip
The curvature $\chi $ of the beam is defined as
{\begin{equation}\label{eq14}
\chi =\frac{\partial \theta }{\partial x}
\end{equation}}\unskip


From Equations~(\ref{eq11})--(\ref{eq13}), we can rewrite the normal and
shear strains as follows
{\begin{equation}\label{eq15}
\varepsilon =\left(1+\frac{\partial u}{\partial x}\right)
\cos \theta +\frac{\partial w}{\partial x}\sin 
\theta -1,\quad \gamma =\frac{\partial w}{\partial x}\cos 
\theta -\left(1+\frac{\partial u}{\partial x}\right)\sin \theta 
\end{equation}}\unskip
It should be emphasized that the strains $\varepsilon $, $\gamma$ and
the curvature $\chi $ although parameterized by the abscissa $x\in
[0,1]$, they take values on the current configuration.

The element strain energy $U$ is contributed to the beam bending
$(U_{\mathrm{B}})$ and foundation deformation $(U_{\mathrm{F}})$, 
$U=U_{\mathrm{B}}+U_{\mathrm{F}}$, in which 
{\begin{equation}\label{eq16}
U_{\mathrm{B}}=\frac{1}{2}{\int }_{0}^{l}\left(A_{11}\varepsilon ^{2}+
2A_{12}\varepsilon \chi +A_{22}\chi ^{2}+\psi A_{33}\gamma ^{2}
\right)\mathrm{d}x
\end{equation}}\unskip
with ${\psi}=5/6$ is the shear correction factor for the present beams
and the rigidities $A_{11}, A_{12}, A_{22}$ and $A_{33}$ are
defined as
{\begin{equation}\label{eq17}
\left(A_{11},A_{12},A_{22}\right)=\int _{A}
(1,z,z^{2})E\,\mathrm{d}A=b\sum_{k=1}^{3}
{\int }_{z_{k-1}}^{z_{k}}(1,z,z^{2})E\,\mathrm{d}z,\quad 
A_{33}=\int _{A}G\,\mathrm{d}A=b\sum_{k=1}^{3}{\int }_{z_{k-1}}^{z_{k}}
G\,\mathrm{d}z
\end{equation}}\unskip
The effective moduli ${E}$ and ${G}$ in the above equation are
dependent on the parameters $\xi$ and $\zeta$, as described in 
Section~\ref{sec2}. One can notice that the coupling rigidity $A_{12}$
in (\ref{eq17}) will vanish when the two face layers are symmetric with
respect to the mid-plane.

Since the beam undergoes large deflections, the foundation deformation
in both the axial and transverse directions should be
considered~\cite{44}. In this regard, the energy stored in the
Pasternak foundation for the large deflection analysis is 
{\begin{equation}\label{eq18}
U_{\mathrm{F}}=\frac{k_{\mathrm{W}}}{2}{\int }_{0}^{L_{F}}(u^{2}+w^{2})
\mathrm{d}x+\frac{k_{G}}{2}{\int }_{0}^{L_{F}}
\left(\theta -\gamma \right)^{2}\mathrm{d}x
\end{equation}}\unskip
where, as above mentioned, $k_{\mathrm{W}}$ and $k_{G}$ are the Winkler and
shear moduli of the Pasternak foundation, and ${L}_{\mathrm{F}}$ is the foundation
supporting length.

Interpolations are needed to be introduced for displacements and
rotations. Since the displacements and rotation are independent from
each another, the following linear interpolations can be employed 
{\begin{equation}\label{eq19}
\left\{\begin{array}{l}
u\\
w\\
\theta \end{array}\right\}=\left[\begin{array}{llllll}
\displaystyle\frac{l-x}{l} & 0 & 0 & \displaystyle\frac{x}{l} & 0 & 0\vspace*{5pt}\\
\displaystyle 0 & \displaystyle\frac{l-x}{l} & 0 & 0 & \displaystyle\frac{x}{l} & 0\vspace*{5pt}\\
\displaystyle 0 & 0 & \displaystyle\frac{l-x}{l} & 0 & 0 & \displaystyle\frac{x}{l}
\end{array}\right]\left\{\begin{array}{l}
u_{I}\\
w_{I}\\
\theta _{I}\\
u_{j}\\
w_{J}\\
\theta _{J}
\end{array}\right\}
\end{equation}}\unskip
The element derived from the linear interpolations (\ref{eq19}), however,
suffers from the problem of shear locking. To avoid this locking
problem, the reduced integration, namely one-point Gauss quadrature, is
used in this work to calculate the terms associated with shear
deformation in the strain energies in Equations~(\ref{eq16}) and 
(\ref{eq18}). In this
regard, the energy expression in Equations~(\ref{eq16}) 
and~(\ref{eq18}) can be written in
the following forms
{\begin{equation}\label{eq20}
\begin{array}{l}
\displaystyle U_{\mathrm{B}}=\frac{l}{2}\left(A_{11}\hat{\varepsilon }^{2}+2A_{12}\hat{\varepsilon }\hat{\chi }+A_{22}\hat{\chi }^{2}+\psi A_{33}\hat{\gamma }^{2}\right)\vspace*{5pt}\\
\displaystyle U_{\mathrm{F}}=\frac{lk_{\mathrm{W}}}{6}\left({u}_{I}^{2}+u_{I}u_{J}+{u}_{J}^{2}+{w}_{I}^{2}+w_{I}w_{J}+{w}_{J}^{2}\right)+\frac{lk_{G}}{2}\left(\hat{\theta }-\hat{\gamma }\right)^{2}
\end{array}
\end{equation}}\unskip
where
{\begin{equation}\label{eq21}
\begin{array}{l}
\displaystyle\hat{\varepsilon }=\left(1+\frac{u_{J}-u_{I}}{l}\right)\cos \hat{\theta }+\frac{w_{J}-w_{I}}{l}\sin \hat{\theta }-1\vspace*{5pt}\\
\displaystyle\hat{\gamma }=-\left(1+\frac{u_{J}-u_{I}}{l}\right)\sin \hat{\theta }+\frac{w_{J}-w_{I}}{l}\cos \hat{\theta }\vspace*{5pt}\\
\displaystyle\hat{\chi }=\frac{\theta _{J}-\theta _{I}}{l}\quad \text{with }\hat{\theta }=\frac{\theta _{I}+\theta _{J}}{2}
\end{array}
\end{equation}}\unskip
From Equation~(\ref{eq21}), one can compute
{\begin{equation}\label{eq22}
\begin{array}{l}
\displaystyle\frac{\partial \hat{\varepsilon }}{\partial u_{I}}=-\frac{1}{l}\cos \hat{\theta },\quad \frac{\partial \hat{\varepsilon }}{\partial u_{J}}=\frac{1}{l}\cos \hat{\theta },\quad \frac{\partial \hat{\varepsilon }}{\partial w_{I}}=-\frac{1}{l}\sin \hat{\theta },\quad \frac{\partial \hat{\varepsilon }}{\partial w_{J}}=\frac{1}{l}\sin \hat{\theta }\vspace*{5pt}\\
\displaystyle\frac{\partial \hat{\gamma }}{\partial u_{I}}=\frac{1}{l}\sin \hat{\theta },\quad \frac{\partial \hat{\gamma }}{\partial u_{J}}=-\frac{1}{l}\sin \hat{\theta },\quad \frac{\partial \hat{\gamma }}{\partial w_{I}}=-\frac{1}{l}\cos \hat{\theta },\quad \frac{\partial \hat{\gamma }}{\partial w_{J}}=\frac{1}{l}\cos \hat{\theta },\vspace*{5pt}\\
\displaystyle\frac{\partial \hat{\varepsilon }}{\partial \theta _{I}}=\frac{\partial \hat{\varepsilon }}{\partial \theta _{J}}=\frac{1}{2}\hat{\gamma },\quad \frac{\partial \hat{\gamma }}{\partial \theta _{I}}=\frac{\partial \hat{\gamma }}{\partial \theta _{J}}=-\frac{1}{2}\left(\hat{\varepsilon }+1\right),\quad \frac{\partial \hat{\chi }}{\partial \theta _{I}}=-\frac{\partial \hat{\chi }}{\partial \theta _{J}}=-\frac{1}{l}
\end{array}
\end{equation}}\unskip
The element vector of nodal internal forces $\mathbf{f}_{\mathrm{in}}$
can be split into two parts, ${\mathbf{f}}_{\mathrm{in}}^{\mathrm{B}}$ due to
beam bending and ${\mathbf{f}}_{\mathrm{in}}^{\mathrm{F}}$ due to foundation
deformation, as follows
{\begin{equation}\label{eq23}
\mathbf{f}_{\mathrm{in}}={\mathbf{f}}_{\mathrm{in}}^{\mathrm{B}}+
{\mathbf{f}}_{\mathrm{in}}^{\mathrm{F}}
\end{equation}}\unskip
These above internal force vectors are computed as derivatives of the
strain energies in~(\ref{eq19}) respecting the element vector of nodal
displacements as
{\begin{equation}\label{eq24}
{\mathbf{f}}_{\mathrm{in}}^{\mathrm{B}}=\frac{\partial U_{\mathrm{B}}}
{\partial \mathbf{d}}=\left\{{f}_{u_{I}}^{\mathrm{B}}\quad 
{f}_{w_{I}}^{\mathrm{B}}\quad {f}_{\theta _{I}}^{\mathrm{B}}\quad 
{f}_{u_{J}}^{\mathrm{B}}\quad {f}_{w_{J}}^{\mathrm{B}}\quad {f}_{\theta _{J}}^{\mathrm{B}}\right\}^{T},
\quad  {\mathbf{f}}_{\mathrm{in}}^{\mathrm{F}}=\frac{\partial U_{\mathrm{F}}}
{\partial \mathbf{d}}=\left\{{f}_{u_{I}}^{\mathrm{F}}\quad 
{f}_{w_{I}}^{\mathrm{F}}\quad {f}_{\theta _{I}}^{\mathrm{F}}\quad 
{f}_{u_{J}}^{\mathrm{F}}\quad {f}_{w_{J}}^{\mathrm{F}}\quad {f}_{\theta _{J}}^{\mathrm{F}}
\right\}^{T}
\end{equation}}\unskip
In substituting (\ref{eq20}) into (\ref{eq24}) and making use of
Equation~(\ref{eq22}), one can obtain explicit expressions for the
above element vectors of nodal forces. The detailed expressions for the
coefficients of these vectors are given by Equations~(\ref{eqA1})  and
(\ref{eqA2}) in the Appendix.

The tangent stiffness matrix $(\mathbf{k}_{t})$ of the element can also
be written in the form
{\begin{equation}\label{eq25}
\mathbf{k}_{t}={\mathbf{k}}_{t}^{B}+{\mathbf{k}}_{t}^{F}
\end{equation}}\unskip
where ${\mathbf{k}}_{t}^{B}$ and ${\mathbf{k}}_{t}^{F}$ are the tangent
stiffness matrices stemming from the deformation of the beam and the
deformation, respectively. These matrices are computed by twice
differentiation of the strain energy (\ref{eq20}) with respect to the vector
\textbf{d}, and they have the forms
{\begin{equation}\label{eq26}
{\mathbf{k}}_{t}^{B}=\frac{\partial ^{2}U_{\mathrm{B}}}{\partial 
\mathbf{d}^{2}}=\left[\begin{array}{llllll}
{k}_{u_{I}u_{I}}^{B} & {k}_{u_{I}w_{I}}^{B} & {k}_{u_{I}
\theta _{I}}^{B} & {k}_{u_{I}u_{J}}^{B} & 
{k}_{u_{I}w_{J}}^{B} & {k}_{u_{I}\theta _{J}}^{B}\vspace*{5pt}\\ 
& {k}_{w_{I}w_{I}}^{B} & {k}_{w_{I}\theta _{I}}^{B} & 
{k}_{w_{I}u_{J}}^{B} & {k}_{w_{I}w_{J}}^{B} & 
{k}_{w_{I}\theta _{J}}^{B}\vspace*{5pt}\\ 
& & {k}_{\theta _{I}\theta _{I}}^{B} & {k}_{\theta _{I}u_{J}}^{B} &
{k}_{\theta _{I}w_{J}}^{B} & {k}_{\theta _{I}\theta _{J}}^{B}\vspace*{5pt}\\ 
& & & {k}_{u_{J}u_{J}}^{B} & {k}_{u_{J}w_{J}}^{B} & {k}_{u_{J}
\theta _{J}}^{B}\vspace*{5pt}\\ 
& \mathrm{sym}. & & & {k}_{w_{J}w_{J}}^{B} & {k}_{w_{J}\theta _{J}}^{B}\\ 
& & & & & {k}_{\theta _{J}\theta _{J}}^{B}
\end{array}\right]
\end{equation}}\unskip
and 
{\begin{equation}\label{eq27}
{\mathbf{k}}_{t}^{F}=\frac{\partial ^{2}U_{\mathrm{F}}}{\partial 
\mathbf{d}^{2}}=\left[\begin{array}{llllll}
{k}_{u_{I}u_{I}}^{F} & {k}_{u_{I}w_{I}}^{F} & {k}_{u_{I}
\theta _{I}}^{F} & {k}_{u_{I}u_{J}}^{F} & {k}_{u_{I}
w_{J}}^{F} & {k}_{u_{I}\theta _{J}}^{F}\vspace*{5pt}\\ 
& {k}_{w_{I}w_{I}}^{F} & {k}_{w_{I}\theta _{I}}^{F} 
& {k}_{w_{I}u_{J}}^{F} & {k}_{w_{I}w_{J}}^{F} & 
{k}_{w_{I}\theta _{J}}^{F}\vspace*{5pt}\\ 
& & {k}_{\theta _{I}\theta _{I}}^{F} & {k}_{\theta _{I}
u_{J}}^{F} & {k}_{\theta _{I}w_{J}}^{F} & {k}_{\theta _{I}
\theta _{J}}^{F}\vspace*{5pt}\\ 
& & & {k}_{u_{J}u_{J}}^{F} & 
{k}_{u_{J}w_{J}}^{F} & {k}_{u_{J}\theta _{J}}^{F}\vspace*{5pt}\\ 
& \mathrm{sym}. & & & {k}_{w_{J}w_{J}}^{F} & {k}_{w_{J}
\theta _{J}}^{F}\\ & & & & & {k}_{\theta _{J}\theta _{J}}^{F}
\end{array}\right]
\end{equation}}\unskip
The detailed expressions for the coefficients of the matrices
${\mathbf{k}}_{t}^{B}$ and ${\mathbf{k}}_{t}^{F}$ are listed by 
Equations~(\ref{eqA3})--(\ref{eqA5}) in the Appendix.

\subsection{Equilibrium equation}\label{sec32}

After assembling the derived matrix $\mathbf{k}_{t}$ and vector
$\mathbf{f}_{\mathrm{in}}$ over the total number of elements, one can
establish the nonlinear equilibrium equation for the beam in the
following {for}m~\cite{45}
{\begin{equation}\label{eq28}
\mathbf{g}(\mathbf{p},\lambda )=\mathbf{q}_{\mathrm{in}}
(\mathbf{p})-\lambda \mathbf{f}_{\mathrm{ef}}
\end{equation}}\unskip
where $\mathbf{g}$ is the residual force vector, is dependent on the
vector of current nodal displacements $\mathbf{p}$ and parameter of the
load level $\lambda$; $\mathbf{q}_{\mathrm{in}}$ is the global vector
of nodal internal forces, obtained by merging vector
$\mathbf{f}_{\mathrm{in}}$ over the total elements, and
$\mathbf{f}_{\mathrm{ef}}$ is the vector of fixed external loads.

The nonlinear equations (\ref{eq28}) is solved herein by Newton--Raphson
iterative method, which contains predictor and corrector phases. In the
method, a new solution is firstly predicted from the previously
converged solution, and then corrected by iterative solutions. The
following Euclidean norm based criterion is used for the iterative
process
{\begin{equation}\label{eq29}
\| \mathbf{g}\| \leq \beta \| \lambda \mathbf{f}_{\mathrm{ef}}\| 
\end{equation}}\unskip
with $\beta $ is a tolerance factor, chosen by $10^{-4}$ herein. The
arc-length technique~\cite{45} is used with the above Newton--Raphson to
trace the equilibrium paths when the beam exhibits complex nonlinear
behavior. 

\section{Results and discussion}\label{sec4}

The sandwich cantilever beam is made from polymethyl pethacrylate
(PMMA) with $E^{m}=2.5$~GPa, $\nu ^{m}=0.34$ as a matrix, and $(10,10)$
single-walled CNTs with the Hill's elastic moduli in  Table~\ref{tab1}
as reinforcement is considered. The beam is loaded by a transverse load
$P$ or a moment $M$ at its free end. For the convenience of discussion,
the non-dimensional parameters are employed 
{\begin{equation}\label{eq30}
P^{*}=\frac{PL^{2}}{E^{m}I},\quad M^{*}=\frac{ML}{E^{m}I},
\quad u^{*}=\frac{u_{L}}{L},\quad w^{*}=
\frac{w_{L}}{L},\quad \theta ^{*}=\frac{\theta }{2\uppi }
\end{equation}}\unskip
where $u_{L}$ and $w_{L}$ are the axial and transverse displacements at
the free end. Also, the dimensionless parameters for foundation
stiffness supporting length are introduced as
{\begin{equation}\label{eq31}
k_{1}=\frac{L^{4}k_{\mathrm{W}}}{E^{m}I},\quad k_{2}=\frac{L^{2}
k_{G}}{\uppi ^{2}E^{m}I},\quad \alpha _{F}=\frac{L_{F}}{L}
\end{equation}}\unskip
with $k_{\mathrm{W}}$, $k_{G}$ and $L_{F}$, as above mentioned, are the
foundation stiffness and the length of the foundation supporting part.
Except for Section~\ref{sec46}, the results are reported for the perfect
beams (without porosities). 

\subsection{Verification}\label{sec41}

The derived formulation using the reduced integration has shown good
performance in analyzing homogeneous beams with different aspect ratios
in Ref.~\cite{41}, the influence of the integration type on the
solution is not presented herewith. The verification will focus on the
accuracy and convergence of the formulation in modeling the large
deflection of the CNT-reinforced composite beam. Firstly, the effective
Young's modulus of the CNT-reinforced composite obtained in the present
work is compared with the result of Daghigh et~al.~\cite{38} in
Figure~\ref{fig4} for a carbon nanotube volume fraction
$V_{\mathrm{CNT}}=0.1$. One can see from the figure that the present
result agrees well with that of the cited reference. It can also be
seen from the figure that Young's modulus considerably decreases by the
CNT agglomeration, and it attains the highest value when  $\xi =\zeta$,
which corresponds to the fully dispersed case.

\begin{figure}
\includegraphics{fig04}
\caption{\label{fig4}Comparison of effective Young's modulus of CNT
reinforced composite with $V_{\mathrm{CNT}}=0.1$.}
\end{figure}

\begin{table}
\caption{\label{tab1}Hill's elastic moduli for the carbon
nanotubes~\cite{27}}
\begin{tabular}{cccccc}
\thead
Radius of CNTs (\AA) & $k_{r}$ (GPa) & 
$l_{r}$ (GPa) & $m_{r}$ (GPa) & $n_{r}$ (GPa) & $p_{r}$ (GPa) \\
\endthead
10 & 30 & 10 & 1 & 450 & 1
\botline
\end{tabular}
\end{table}

In Table~\ref{tab2}, the dimensionless tip displacements of a (1-2-1)
CNTRC sandwich cantilever beam with $L/h=20$ under an end
moment $M^{*}=5$ obtained herein are compared
with the analytical solution of Refs.~\cite{46,47} for various
values of the agglomerated parameters. The tip displacements of
Refs.~\cite{46,47} are derived for an Euler--Bernoulli FG cantilever
beam, and they have the following forms
{\begin{equation}\label{eq32}
u_{L}=L\left[1-\frac{A_{22}}{ML}\sin \left(\frac{ML}{A_{22}}\right)
\right],\quad w_{L}=\frac{A_{22}}{M}\left[1-\cos \left(\frac{ML}
{A_{22}}\right)\right]
\end{equation}}\unskip
where $A_{22}$ is the bending stiffness, defined by
Equation~(\ref{eq17}). An excellent agreement between the present
result with the analytical solution of Refs.~\cite{46,47} is noted from
Table~\ref{tab2}, regardless of the agglomeration parameters. It is
noted that, for the symmetric (1-2-1) sandwich beam under
consideration, the coupling rigidity vanished, $A_{12}=0$, and the
analytical solution is derived by taking into account the neutral
surface position, which allows to elimination of the coupling term.

\begin{table}
\caption{\label{tab2}Comparison of tip displacements of symmetric
(1-2-1) sandwich cantilever beam under end moment $M^{*}=5$, $L/h=20$
and $V_{\mathrm{CNT}}=0.1$}
\tabcolsep=2pt
\begin{tabular}{cccccccccc}
\thead
\xmorerows{1}{$\zeta$} & & \multicolumn{2}{c}{$\xi =0.1$} & 
\multicolumn{2}{c}{$\xi =0.3$} &   \multicolumn{2}{c}{$\xi =0.5$} &
\multicolumn{2}{c}{$\xi =0.7$} \\\cline{3-10}
&& Present & Refs.~\cite{46,47} & Present & Refs.~\cite{46,47} &  
Present & Refs.~\cite{46,47} & Present & Refs.~\cite{46,47} \\
\endthead
\morerows{1}{0.1}  & $|u^*|$ & 0.2808 & 0.2810 & 0.3092 & 0.3093 &   0.3587 & 0.3589 & 0.4169 & 0.4171 \\
 & $w^{*}$ & 0.5815 & 0.5814 & 0.6025 & 0.6024 &   0.6342 & 0.6340 & 0.6644\vspace*{3pt} & 0.6642 \\
\morerows{1}{0.3}  & $|u^*|$ & 0.3238 & 0.3239 & 0.2808 & 0.2810 &   0.3004 & 0.3006 & 0.3464 & 0.3466 \\
 & $w^{*}$ & 0.6124 & 0.6123 & 0.5815 & 0.5814 &   0.5963 & 0.5962 & 0.6268\vspace*{3pt} & 0.6267 \\
\morerows{1}{0.5}  & $|u^*|$ & 0.4379 & 0.4381 & 0.3078 & 0.3079 &   0.2808 & 0.2810 & 0.3001 & 0.3002 \\
 & $w^{*}$ & 0.6737 & 0.6734 & 0.6015 & 0.6014 &   0.5815 & 0.5814 & 0.5960\vspace*{3pt} & 0.5959 \\
\morerows{1}{0.7}  & $|u^*|$ & 0.6551 & 0.6553 & 0.4096 & 0.4098 &   0.3084 & 0.3085 & 0.2808 & 0.2810 \\
 & $w^{*}$ & 0.7241 & 0.7237 & 0.6610 & 0.6608 &   0.6019 & 0.6018 & 0.5815\vspace*{3pt} & 0.5814\\
\morerows{1}{0.9}  & $|u^*|$ & 1.0323 & 1.0323 & 0.6942 & 0.6943 &   0.4609 & 0.4610 & 0.3289 & 0.3291 \\
 & $w^{*}$ & 0.6150 & 0.6143 & 0.7250 & 0.7246 &   0.6828 & 0.6826 & 0.6158 & 0.6157 
\botline
\end{tabular}
\end{table}

The convergence of the derived beam element in predicting the large
deflection response of the CNTRC beams is shown in Table~\ref{tab3},
where the dimensionless deflections at the free end of symmetric
(2-1-2) and nonsymmetric (1-1-2) sandwich cantilever beams with
different aspect ratios under the tip load obtained by different
{number of} elements {(nELE)} are given for $P^{*}=10$,
${V}_{\mathrm{CNT}}=0.3$, ${\alpha}_{\mathrm{F}}=0.3$, $k_{1}=50$ and $k_{2}=0.2$.
As seen from the table the convergence of both the symmetric (2-1-2)
and non-symmetric (1-1-2) beams is achieved by using 24 elements,
regardless of the agglomeration parameters and the aspect ratios.
Because of this convergence result, a mesh of 24 elements is used in
the computations reported below. 

\begin{table}
\caption{\label{tab3}Convergence of the formulation in predicting tip
deflection $w^{*}$ of the cantilever beam subjected to the tip load for
$P^{*}=10$, $V_{\mathrm{CNT}}=0.3$, $\alpha_{\mathrm{F}}=0.5$, $k_{1}=50$ and
$k_{2}=0.2$}
\tabcolsep=3pt
\begin{tabular}{cccccccccc}
\thead
\xmorerows{1}{$(\zeta,\xi)$}  & \xmorerows{1}{nELE}  & \multicolumn{4}{c}{(2-1-2) beam} & \multicolumn{4}{c}{(1-1-2) beam} \\\cline{3-10}
 &  & $L/h=50$ & $L/h=20$ & $L/h=10$ & $L/h=5$ & $L/h=50$ & $L/h=20$ & $L/h=10$ & $L/h=5$ \\
\endthead
\morerows{6}{(0.1, 0.4)}  & \08 & 0.5693 & 0.5703 & 0.5740 & 0.5889 & 0.7243 & 0.7276 & 0.7351 & 0.7565 \\
 & 12 & 0.5697 & 0.5707 & 0.5744 & 0.5892 & 0.7245 & 0.7277 & 0.7352 & 0.7578 \\
 & 16 & 0.5698 & 0.5708 & 0.5745 & 0.5894 & 0.7245 & 0.7278 & 0.7353 & 0.7581 \\
 & 20 & 0.5698 & 0.5709 & 0.5746 & 0.5894 & 0.7245 & 0.7278 & 0.7353 & 0.7581 \\
 & 22 & 0.5699 & 0.5709 & 0.5746 & 0.5895 & 0.7246 & 0.7278 & 0.7353 & 0.7581 \\
 & 24 & 0.5699 & 0.5709 & 0.5746 & 0.5895 & 0.7246 & 0.7278 & 0.7353 & 0.7581 \\
 & 26 & 0.5699 & 0.5709 & 0.5746 & 0.5895 & 0.7246 & 0.7278 & 0.7353\vspace*{3pt} & 0.7581 \\
\morerows{6}{(0.1,0.6)}  & \08 & 0.5945 & 0.5956 & 0.5965 & 0.6153 & 0.7397 & 0.7431 & 0.7508 & 0.7751 \\
 & 12 & 0.5948 & 0.5959 & 0.5986 & 0.6156 & 0.7398 & 0.7432 & 0.7511 & 0.7753 \\
 & 16 & 0.5955 & 0.5960 & 0.5998 & 0.6157 & 0.7415 & 0.7432 & 0.7512 & 0.7753 \\
 & 20 & 0.5958 & 0.5961 & 0.6000 & 0.6158 & 0.7419 & 0.7433 & 0.7512 & 0.7753 \\
 & 22 & 0.5978 & 0.5961 & 0.6001 & 0.6188 & 0.7424 & 0.7433 & 0.7512 & 0.7753 \\
 & 24 & 0.5978 & 0.5961 & 0.6001 & 0.6158 & 0.7429 & 0.7433 & 0.7512 & 0.7753 \\
 & 26 & 0.5978 & 0.5961 & 0.6001 & 0.6158 & 0.7429 & 0.7433 & 0.7512 & 0.7753
\botline
\end{tabular}
\end{table}

\subsection{Effect of CNT reinforcement}\label{sec42}

In Table~\ref{tab4}, the dimensionless tip displacements of (2-1-2) and
(1-1-2) sandwich beams corresponding to $P^{*}=10$ are given for
$P^{*}=10$, $L/h=10$, ${\alpha}_{\mathrm{F}}=0.4$,
$k_{1}=50$, and $k_{2}=0.5$. The influence of the CNT percentage
and the agglomeration parameters on the beam response can be observed
clearly from the table. As expected, the increase of
$V_{\mathrm{CNT}}$ leads to a noticeable decrease in the tip
displacements of both the symmetric and non-symmetric sandwich beams,
regardless of the agglomeration parameters. The agglomeration degree of
CNTs, indicated by the two parameters $\xi$ and $\zeta $, also has an
impact on the beam response. The tip displacements of both the
symmetric (2-1-2) and non-symmetric (2-2-1) beams attain the smallest
values when $\xi =\zeta $, which corresponds to the fully dispersed
case of CNTs. As mentioned above and shown in Figure~\ref{fig4}, the
effective Young's modulus of the composite {attains} the highest value
when the CNTs are fully dispersed, and thus the beam stiffness is
largest in this case. The tip response of both the sandwich beams is
significant when the degree of CNT agglomeration is severe, which is
the case of the two agglomerated parameters $\xi$ and $\zeta$ are far
from each other. A more detailed examination of the table shows that
the effect of the CNT agglomeration is, however dependent on the CNT
volume fraction, and this effect is more prominent for the beam
associated with a higher CNT volume fraction. For example, the
difference in the transverse displacements of the (2-1-2) beam
associated with $V_{\mathrm{CNT}}=0.05$ for the cases of severe
agglomeration (${\xi}=0.1$, ${\zeta}=0.9$) and the full dispersion of
CNTs (${\xi}={\zeta}$) is 6.56\%, while the corresponding value for
the beam associated with $V_{\mathrm{CNT}}=0.3$ is 33.41\%, which is
more than five times higher.

\begin{table}
\caption{\label{tab4}Dimensionless tip displacements of CNTRC
sandwich cantilever beam for $P^{*}=10$, $L/h=10$,
${\alpha}_{\mathrm{F}}=0.4$, $k_{1}=50$ and $k_{2}=0.5$}
\begin{tabular}{ccccccccccc}
\thead
\xmorerows{1}{$V_{\mathrm{CNT}}$}  & \xmorerows{1}{$\zeta $}  
&\multicolumn{1}{c}{}&  \multicolumn{4}{c}{(2-1-2)} &   \multicolumn{4}{c}{(1-1-2)} \\\cline{4-11}
 &  &  & $\xi =0.1$ & $\xi =0.3$ & $\xi =0.5$ & $\xi =0.7$ &   $\xi =0.1$ & $\xi =0.3$ & $\xi =0.5$ & $\xi =0.7$ \\
\endthead
\morerows{9}{0.05}  & \morerows{1}{0.1}  & $|u^*|$ & 0.5350 & 0.5421 & 0.5540 & 0.5679 &  0.6587 & 0.6636 & 0.6719 & 0.6815 \\
 &  & $w^{*}$ & 0.8088 & 0.8126 & 0.8190 & 0.8265 &   0.8774 & 0.8800 & 0.8843 & 0.8895 \\
 & \morerows{1}{0.3}  & $|u^*|$ & 0.5463 & 0.5350 & 0.5402 &  0.5522 &   0.6665 & 0.6587 & 0.6623 & 0.6706 \\
 &  & $w^{*}$ & 0.8149 & 0.8088 & 0.8116 & 0.8181 &   0.8815 & 0.8774 & 0.8793 &  0.8837 \\
 & \morerows{1}{0.5}  & $|u^*|$ & 0.5708 & 0.5418 & 0.5350 & 0.5402 &   0.6834 & 0.6634 & 0.6587 &  0.6623 \\
 &  & $w^{*}$ & 0.8280 & 0.8125 & 0.8088 & 0.8116 &   0.8905 & 0.8799 & 0.8774 & 0.8793 \\
 & \morerows{1}{0.7}  & $|u^*|$ & 0.6031 & 0.5625 & 0.5417 & 0.5350 &   0.7055 & 0.6777 & 0.6633 & 0.6587 \\
 &  & $w^{*}$ &  0.8449 & 0.8236 & 0.8124 & 0.8088 &   0.9027 & 0.8874 &  0.8798 & 0.8774 \\
 & \morerows{1}{0.9}  & $|u^*|$ & 0.6425 & 0.5987 & 0.5673 & 0.5452 &   0.7321 & 0.7025 & 0.6810 & 0.6658 \\
 &  & $w^{*}$ & 0.8656 & 0.8427 & 0.8261 & 0.8143 &  0.9182 & 0.9010 & 0.8892\vspace*{3pt} & 0.8811 \\
\morerows{9}{0.3}  & \morerows{1}{0.1}  & $|u^*|$ & 0.1953 & 0.2132 & 0.2391 & 0.2577 &   0.3729 & 0.3925 & 0.4194 & 0.4377 \\
 &  & $w^{*}$ & 0.5406 & 0.5618 & 0.5902 & 0.6093 &   0.7139 & 0.7279 & 0.7459 & 0.7577 \\
 & \morerows{1}{0.3}  & $|u^*|$ & 0.2214 & 0.1953 & 0.2073 & 0.2298 &   0.4013 & 0.3729 & 0.3862 & 0.4099 \\
 &  & $w^{*}$ & 0.5711 & 0.5406 & 0.5550 & 0.5803 &   0.7338 & 0.7139 & 0.7234 & 0.7397 \\
 & \morerows{1}{0.5}  & $|u^*|$ & 0.2830 & 0.2129 & 0.1953 & 0.2066 &   0.4616 & 0.3922 & 0.3729 & 0.3854 \\
 &  & $w^{*}$ & 0.6336 &  0.5614 & 0.5406 & 0.5541 &   0.7723 & 0.7276 & 0.7139 & 0.7229 \\
 & \morerows{1}{0.7}  & $|u^*|$ & 0.3823 & 0.2767 & 0.2142 & 0.1953 &   0.5460 & 0.4557 & 0.3936 & 0.3729 \\
 &  & $w^{*}$ & 0.7143 & 0.6277 & 0.5629 & 0.5406 &   0.8195 & 0.7688 &  0.7286 & 0.7139 \\
 & \morerows{1}{0.9}  & $|u^*|$ & 0.5405 & 0.4348 & 0.3240 & 0.2324 &   0.6625 & 0.5864 & 0.4980 & 0.4126 \\
 &  & $w^{*}$ & 0.8118 & 0.7499 & 0.6694 & 0.5831 &  0.8794 &  0.8404 & 0.7934 & 0.7414
\botline
\end{tabular}
\end{table}

To show the impact of CNT reinforcement on the beam response in more
further, the effects of CNT volume fraction on the large deflection
response of (2-1-2) symmetric and (1-1-2) non-symmetric sandwich
cantilever beam under the tip load ${P}$ and end moment $M$ are
respectively displayed in Figures~\ref{fig5} and~\ref{fig6} for 
$L/h=10$, ${\alpha}_{\mathrm{F}}=0.4$, $(\xi, \zeta)=(0.4,0.6)$ and
$(k_{1},k_{2})=(50,0.5)$. The effect of the degree of CNT agglomeration
on the large deflections is shown in Figures~\ref{fig7} and~\ref{fig8},
where the tip response of the symmetric (2-1-2) and non-symmetric
(1-1-2) beams due to the tip transverse load and the end moment are
depicted for $L/h=10$, ${V}_{\mathrm{CNT}}=0.1$, ${\alpha}_{\mathrm{F}}=0.4$,
${\zeta}=0.9$, $(k_{1}, k_{2})=(50,0.5)$ and different values of the
agglomeration parameter ${\xi}$. Figure~\ref{fig5} shows a decline in
the free end displacement of the sandwich cantilever beams by
increasing the volume fraction of CNTs, regardless of the load level
and the sandwich configuration. The influence of the CNT volume
fraction on the response of the sandwich beam loaded by the end moment
is seen in Figure~\ref{fig6}, and the beam is more conservative in the
case of a higher CNT volume fraction. The cantilever beam rolls up a
half circle when the applied moment reaches a certain value, and after
that the displacement $w^{*}$ starts to decrease, as seen in
Figure~\ref{fig6}b for the curve corresponding to
${V}_{\mathrm{CNT}}=0.05$. 

\begin{figure}
\includegraphics{fig05}
\caption{\label{fig5}Effect of $V_{\mathrm{CNT}}$ on tip response of
CNTRC sandwich beam under transverse load for $L/h=10$, ${\alpha}_{\mathrm{F}}=0.4$, 
$({\xi}, {\zeta})=(0.4,0.6)$ and $(k_{1},k_{2})=(50,0.5)$.}
\end{figure}

\begin{figure}
\includegraphics{fig06}
\caption{\label{fig6}Effect of $V_{\mathrm{CNT}}$ on tip response of
CNTRC sandwich beam under end moment for $L/h=10$, ${\alpha}_{\mathrm{F}}=0.4$,
$(\xi,\zeta)=(0.4,0.6)$ and $(k_{1},k_{2})=(50,0.5)$.}
\end{figure}

\begin{figure}
\includegraphics{fig07}
\caption{\label{fig7}The tip response of CNTRC sandwich beam to
transverse load for $L/h=10$, ${V}_{\mathrm{CNT}}=0.1$, ${\alpha}_{\mathrm{F}}=0.4$,
${\zeta}=0.9$, $(k_{1},k_{2})=(50,0.5)$ and different values of
${\xi}$.}
\vspace*{5pt}
\end{figure}

The tip response of the sandwich {cantilever} beam with the symmetric
and non-symmetric configurations, as seen from 7 and 8, is more
pronounced when the difference between the two agglomeration parameters
$\xi$ and $\zeta $ is larger. In other words, the more severe the
agglomeration degree is, the larger tip displacements are, regardless
of the loading type and the sandwich configuration. It is noted that
the cantilever beam associated with ${\xi}=0.3$ has rolled into a half
circle when ${M}^{*}=2$, and the displacement $w^{*}$ decreases by
increasing the moment some further, as seen in Figure~\ref{fig8}b.
Thus, though the displacement $w^{*}$ of the beam with ${\xi}=0.3$
loaded by ${M}^{*}>2$ is smaller, it still deforms more significantly
compared to the other cases (${\xi}=0.5$, 0.7, 0.9). The influence of
the {CNT} reinforcement can also be observed from Figures~\ref{fig9}
and~\ref{fig10}, where the configurations in the deformed state of
(2-1-2) and (1-1-2) beams due to the end moment are shown for $L/h=10$,
${\alpha}_{\mathrm{F}}=0.4$, and $(k_{1},k_{2})=(50,0.5)$. At the same applied
moment parameter $M^{*}=5$, the non-symmetric \mbox{(1-1-2)} beam with
$V_{\mathrm{CNT}}=0.05$ has already rolled into a circle, while the
deformed curve of the (1-1-2) beam with $V_{\mathrm{CNT}}=0.3$ is still
far from a haft circle (Figure~\ref{fig9}b). The influence of the CNT
volume fraction is seen more clearly from Table~\ref{tab5} where the
dimensionless tip displacements and rotation of the sandwich beam with
$L/h=10$, $(\xi,\zeta)=(0.4,0.7)$, $\alpha_{\mathrm{F}}=0.4$, 
$(k_{1},k_{2})=(50,0.5)$ are given for $M^*=3$ and $M^*=5$. One can see
that the beam associated with a lower ${V}_{\mathrm{CNT}}$ has larger
tip rotation, regardless of the beam type and the applied moment. In
other words, the beam with lower ${V}_{\mathrm{CNT}}$ deforms more
significantly. The deformed configuration of the sandwich cantilever
beams, as can be seen in Figure~\ref{fig10}, is also significantly
influenced by the agglomeration degree of CNTs. At $M^{*}=5$, the
non-symmetric \mbox{(1-1-2)} sandwich beam with fully dispersed CNTs rolls
into a haft circle, while the same beam with a severe degree of CNTs
agglomeration (${\xi}=0.3$, ${\zeta}=0.9$) has already rolled into a
circle (Figure~\ref{fig10}b).

\begin{figure}
\includegraphics{fig08}
\caption{\label{fig8}The tip response of CNTRC sandwich beam to the end
moment for $L/h=10$, ${V}_{\mathrm{CNT}}=0.1$, ${\alpha}_{\mathrm{F}}=0.4$,
${\zeta}=0.9$, $(k_{1},k_{2})=(50,0.5)$ and different values of ${\xi}$.}
\end{figure}

\begin{figure}
\includegraphics{fig09}
\vspace*{-5pt}
\caption{\label{fig9}Deformed configurations of cantilever
corresponding to $M^{*}=5$ for $L/h=10$,
$(\xi,\zeta)=(0.4,0.7)$, ${\alpha}_{\mathrm{F}}=0.4$,   $(k_{1},k_{2})=(50,0.5)$
and different values of ${V}_{\mathrm{CNT}}$.}
\end{figure}

\begin{figure}
\includegraphics{fig10}
\caption{\label{fig10}Deformed configurations of cantilever
corresponding to $M^{*}=5$ for $L/h=10$, $\zeta=0.9$, 
${\alpha}_{\mathrm{F}}=0.4$, $(k_{1},k_{2})=(50,0.5)$ and different values of
{${\xi}$}.}
\end{figure}

\begin{table}
\caption{\label{tab5}Dimensionless tip displacements and rotation of
CNTRC sandwich beam under end moment for $L/h=10$, $(\xi,\zeta)=(0.4,0.7)$, ${\alpha}_{\mathrm{F}}=0.4$,
$(k_{1},k_{2})=(50,0.5)$}
\begin{tabular}{ccccccc}
\thead
\xmorerows{1}{${M}^{*}$}  & \xmorerows{1}{Beam type}  &  &
\multicolumn{4}{c}{$V_{\mathrm{CNT}}$} \\\cline{3-7}
 &  &  & 0.05 & 0.07 & 0.1 & 0.3 \\
\endthead
\morerows{5}{3}  & \morerows{2}{(2-1-2)}  & $| u^{*}| $ & 0.4689 & 0.3414 & 0.2279 & 0.0441 \\
 &  & $w^{*}$ & 0.6721 & 0.6144 & 0.5301 & 0.2523 \\
 &  & ${\theta}^{*}$ & 0.2972 & 0.2461 & 0.1962 & 0.0830 \\
 & \morerows{2}{(1-1-2)}  & $| u^{*}| $ & 1.1029 & 0.9481 & 0.7280 & 0.1821 \\
 &  & $w^{*}$ & 0.4587 & 0.6101 & 0.6863 & 0.4730 \\
 &  & ${\theta}^{*}$ & 0.5997 & 0.4957 & 0.3959 & 0.1704\vspace*{5pt} \\
\morerows{5}{5}  & \morerows{2}{(2-1-2)}  & $| u^{*}| $ & 0.9734 & 0.7774 & 0.5580 & 0.1198 \\
 &  & $w^{*}$ & 0.6262 & 0.7019 & 0.7010 & 0.4040 \\
 &  & ${\theta}^{*}$ & 0.5008 & 0.4129 & 0.3282 & 0.1385 \\
 & \morerows{2}{(1-1-2)}  & $| u^{*}| $ & 0.9525 & 1.1242 & 1.1766 & 0.4513 \\
  &  & $w^*$ & 0.0022 & 0.0708 & 0.3322 & 0.6603 \\
 &  & ${\theta}^{*}$ & 1.0412 & 0.8550 & 0.6748 & 0.2848
\botline
\end{tabular}
\end{table}

\subsection{Effect of foundation support}\label{sec43}

Figure~\ref{fig11} shows the load-deflection curves of the (2-1-2)
symmetric and (1-1-2) non-symmetric sandwich cantilever beams under the
tip load $P$ for $L/h=10$, ${V}_{\mathrm{CNT}}=0.1$,
$(\xi,\zeta)=(0.4,0.7)$, ${\alpha}_{\mathrm{F}}=0.4$ and different foundation
stiffness parameters. The effect of the foundation supporting length is
shown in Figure~\ref{fig12} for $L/h=10$, ${V}_{\mathrm{CNT}}=0.1$,
$(\xi,\zeta)=(0.4,0.7)$, and $(k_{1},k_{2})=(50,0.5)$. As expected, the
tip displacements are declined by increasing the foundation stiffness,
irrespective of the load level and the sandwich configuration. The
influence of the foundation supporting length on the large deflections
of the sandwich beam is similar to that of the CNT volume fraction. One
can see from Figures~\ref{fig11} and~\ref{fig12} that the increase of
the supporting parameter $\alpha _{F}$ mitigates the tip displacements
of the sandwich cantilever beam significantly, regardless of the
loading type and the sandwich configuration.

\begin{figure}
{\vspace*{2pt}}
\includegraphics{fig11}
{\vspace*{2pt}}
\caption{\label{fig11}Effect of foundation stiffness on large
deflections of CNTRC sandwich beams under transverse load for $L/h=10$,
${V}_{\mathrm{CNT}}=0.1$, $(\xi,\zeta)=(0.4,0.7)$, ${\alpha}_{\mathrm{F}}=0.4$.}
\vspace*{5pt}
\end{figure}

\begin{figure}
\includegraphics{fig12}
{\vspace*{2pt}}
\caption{\label{fig12}Effect of foundation supporting parameter on
large deflections of sandwich beam under transverse load for $L/h=10,
V_{\mathrm{CNT}}=0.1$, $(\xi,\zeta)=(0.4,0.7)$,
$(k_{1},k_{2})=(50,0.5)$.}
{\vspace*{2pt}}
\end{figure}

\subsection{Effect of sandwich configuration}\label{sec44}

The sandwich configuration, as can be seen from Table~\ref{tab4}, has a
significant impact on the large deflection response of the CNTRC
sandwich cantilever beam. At the same volume fraction of the CNTs, the
tip displacements of a non-symmetric (1-1-2) beam are considerably
higher than that of the symmetric (2-1-2) beam, regardless of the CNT
volume fraction and the agglomerated parameters. The deformed
configurations of the sandwich beams under the end moment as depicted
in Figures~\ref{fig9} and~\ref{fig10} also show the more conservative
of the symmetric (2-1-2) beam compared to the non-symmetric (1-1-2)
one. As seen from Figure~\ref{fig9}, at the same end moment parameter
${M}^{*}=5$, the non-symmetric (1-1-2) beam with
${V}_{\mathrm{CNT}}=0.05$ has rolled into a full circle, while the
symmetric (2-1-2) beam just formed into a half circle only. The same
situation is seen for the configurations of the (2-1-2) and (1-1-2)
beams associated with ${\xi}=0.05$ under the end moment ${M}^{*}=5$ in
Figure~\ref{fig10}. Thus, at given values of the CNT volume fraction
and the agglomeration parameters, one can design a sandwich beam to
mitigate the large deformation by the appropriate choice of the layer
thickness ratio.

\begin{table}
\caption{\label{tab6}Tip displacements of CNTRC sandwich cantilever
beams with different length-to-height ratios for $P^{*}=10$,
${V}_{\mathrm{CNT}}=0.1$, ${\alpha}_{\mathrm{F}}=0.4$ and
$(k_{1},k_{2})=(50,0.5)$\vspace*{2pt}}
\begin{tabular}{ccccccccccc}
\thead
\xmorerows{1}{$L/h$}  & \xmorerows{1}{$\zeta $}  & \multicolumn{1}{c}{} &  \multicolumn{4}{c}{(2-1-2)} &   \multicolumn{4}{c}{(1-1-2)} \\\cline{4-11}
 &  &  & $\xi =0.1$ & $\xi =0.3$ & $\xi =0.5$ & $\xi =0.7$ &   $\xi =0.1$ & $\xi =0.3$ & $\xi =0.5$ & $\xi =0.7$ \\
\endthead
\morerows{7}{\05}  & \morerows{1}{0.1}  & $| u^{*}| $ & 0.4326 & 0.4459 & 0.4664 & 0.4871 &  0.5828 & 0.5927 & 0.6078 & 0.6229 \\
 &  & $w^{*}$ & 0.7707 & 0.7803 & 0.7947 & 0.8090 &   0.8749 &  0.8813 & 0.8913 & 0.9014 \\
 & \morerows{1}{0.3}  & $| u^{*}| $ & 0.4522 & 0.4326 & 0.4419 & 0.4616 &   0.5974 & 0.5828 & 0.5897 & 0.6043 \\
 &  & $w^{*}$ & 0.7848 & 0.7707 & 0.7774 & 0.7913 &   0.8844 & 0.8749 & 0.8794 &  0.8889 \\
 & \morerows{1}{0.5}  & $| u^{*}| $ & 0.4939 & 0.4452 & 0.4326 & 0.4418 &   0.6278 & 0.5922 & 0.5828 & 0.5896 \\
 &  & $w^{*}$ & 0.8137 & 0.7798 & 0.7707 & 0.7773 &   0.9047 & 0.8810 & 0.8749 & 0.8793 \\
 & \morerows{1}{0.7}  & $| u^{*}| $ & 0.5505 & 0.4847 & 0.4455 & 0.4326 &   0.6675 & 0.6211 & 0.5924 & 0.5828 \\
 &  & $w^{*}$ & 0.8519 & 0.8074 & 0.7800 & 0.7707 &   0.9337 & 0.9002 & 0.8812\vspace*{5pt} & 0.8749 \\
\morerows{7}{10}  & \morerows{1}{0.1}  & $| u^{*}| $ & 0.4280 & 0.4412 & 0.4617 & 0.4824 &  0.5813 & 0.5912 & 0.6063 & 0.6214 \\
 &  & $w^{*}$ & 0.7455 & 0.7540 & 0.7667 & 0.7791 &   0.8378 & 0.8429 & 0.8506 & 0.8582 \\
 & \morerows{1}{0.3}  & $| u^{*}| $ & 0.4475 & 0.4280 & 0.4373 & 0.4569 &   0.5959 & 0.5813 & 0.5883 & 0.6028 \\
 &  & $w^{*}$ & 0.7579 & 0.7455 & 0.7515 & 0.7637 &   0.8453 & 0.8378 & 0.8414 & 0.8488 \\
 & \morerows{1}{0.5}  & $| u^{*}| $ & 0.4892 & 0.4406 & 0.4280 & 0.4371 &   0.6263 & 0.5907 & 0.5813 & 0.5881 \\
 &  & $w^{*}$ & 0.7830 & 0.7536 & 0.7455 & 0.7513 &   0.8607 & 0.8426 & 0.8378 & 0.8413 \\
 & \morerows{1}{0.7}  & $| u^{*}| $ & 0.5457 & 0.4800 & 0.4408 & 0.4280 &   0.6661 & 0.6196 & 0.5909 & 0.5813 \\
 &  & $w^{*}$ & 0.8146 & 0.7776 & 0.7537 & 0.7455 &   0.8813 &  0.8574 & 0.8427\vspace*{5pt} & 0.8378 \\
\morerows{7}{20}  & \morerows{1}{0.1}  & $| u^{*}| $ & 0.4268 & 0.4400 & 0.4605 & 0.4812 &  0.5813 & 0.5912 & 0.6063 & 0.6214 \\
 &  & $w^{*}$ & 0.7391 & 0.7474 & 0.7596 & 0.7715 &   0.8265 & 0.8312 & 0.8384 & 0.8454 \\
 & \morerows{1}{0.3}  & $| u^{*}| $ & 0.4464 & 0.4268 & 0.4361 &  0.4557 &   0.5959 & 0.5813 & 0.5882 & 0.6028 \\
 &  & $w^{*}$ & 0.7512 & 0.7391 & 0.7449 & 0.7568 &   0.8335 & 0.8265 & 0.8298 &  0.8367 \\
 & \morerows{1}{0.5}  & $| u^{*}| $ & 0.4880 & 0.4394 & 0.4268 & 0.4359 &   0.6263 & 0.5907 & 0.5813 & 0.5881 \\
 &  & $w^{*}$ & 0.7753 & 0.7470 & 0.7391 & 0.7448 &   0.8477 & 0.8310 & 0.8265 & 0.8298 \\
 & \morerows{1}{0.7}  & $| u^{*}| $ & 0.5445 & 0.4788 & 0.4397 & 0.4268 &   0.6661 & 0.6196 & 0.5909 & 0.5813 \\
 &  & $w^{*}$ & 0.8052 & 0.7702 & 0.7472 & 0.7391 &   0.8660 & 0.8446 & 0.8311 & 0.8265
\botline
\end{tabular}
\vspace*{4pt}
\end{table}

\subsection{Effect of slenderness ratio}\label{sec45}

The length-to-height ratio, (${L}/{h}$), of {a rectangular
cross-section beam} represents the beam {slenderness ratio.} The
dimensionless displacements at the free end of the (2-1-2) and (1-1-2)
sandwich beams corresponding to a load parameter $P^{*}=10$ are
tabulated in Table~\ref{tab6} for $P^{*}=10$, ${V}_{\mathrm{CNT}}=0.1$,
${\alpha}_{\mathrm{F}}=0.4$, $(k_{1},k_{2})=(50,0.5)$ and different values of the
length-to-height ratio and the agglomeration parameters. The table
shows that the dimensionless tip displacements of both the sandwich
cantilever beams are higher for the beams associated with a smaller
length-to-height ratio, regardless of the agglomeration parameters and
the sandwich configuration. One can also see from Table~\ref{tab6} that
the effect of the slenderness ratio is dependent on the agglomeration
parameters, which is on the agglomeration degree of CNTs. For example,
for the fully dispersed case (${\xi}={\zeta}=0.1$), the difference in
the deflections of the symmetric (2-1-2) beam with $L/h=5$ and $L/h=20$
is 4.10\% but this number increases to 5.48\% for (${\xi}=0.1$,
${\zeta}=0.7$), the case that the agglomeration degree of CNTs is the
most severe. Since the length-to-height ratio indicates the shear
deformation effect, one can say that the influence of the shear
deformation is more remarkable when the agglomeration degree of CNTs is
more severe. Examining Table~\ref{tab5} in more detail one can see that
the effect of the slenderness ratio of the non-symmetric (1-1-2)
sandwich beam is more significant than that of the symmetric (2-1-2)
sandwich beam, regardless of the agglomeration degree. The difference
between the deflections of non-symmetric (1-1-2) beams with $L/h=5$ and
$L/h=20$ is 5.53\% for the case of fully dispersed CNTs
(${\xi}={\zeta}=0.1$), and that is 7.25\% for the case of severe degree
of CNTs agglomeration (${\xi}=0.1$, ${\zeta}=0.7$). These values are
slightly higher than that of the symmetric (2-1-2) beam as stated
above. The result in Table~\ref{tab5} reveals the ability of the
derived formulation to model the influence of shear deformation on the
large deflection response of the CNTRC sandwich beams.

\subsection{Effect of porosities}\label{sec46}

The influence of porosities on large deflections of the CNTRC sandwich
beam is examined in this Sub-section. Following the work in
Ref.~\cite{34}, the response of the beam is predicted for various
values of the porosity parameter $f_{\mathrm{p}}$. In Table~\ref{tab7}, the
dimensionless tip displacements and rotation of the porous beam with
various porosity parameters are given for $L/h=10$, 
${V}_{\mathrm{CNT}}=0.07$,  ${\alpha}_{\mathrm{F}}=0.4$, $(k_{1},k_{2})=(50,0.5)$,
and two values of the end moment, $M^*=4$ and $M^*=8$. The impact of
the porosities on the large deflections is seen from the table, and the
higher the porosity parameter is, the larger the tip rotation the beam
has. The effect of the porosities on the large deformation of the CNTRC
sandwich beam can also be seen in Figure~\ref{fig13}, where the tip
response and the deformed configurations of {the} (2-1-2) sandwich beam
under a tip moment $M^*=8$ are depicted for $L/h=10$,
${V}_{\mathrm{CNT}}=0.07$, ${\xi}=0.4, \zeta=0.7, \alpha_{\mathrm{F}}=0.4$,
$(k_{1},k_{2})=(50,0.5)$ and different values of the porosity
parameter. Figure~\ref{fig13}b shows that the beam associated with a
higher porosity parameter deforms more significantly, though the
deflection corresponding to $M^*=8$ in Figure~\ref{fig13}a decreases
by increasing $f_{\mathrm{p}}$. The increase of the tip rotation by increasing the
porosity parameter as seen in Table~\ref{tab7} also confirms the
influence of the porosities on the large deflections of the beam.   

\begin{table}
\caption{\label{tab7}Effect of porosities on the tip response of the
sandwich cantilever beam under end moment for $L/h=10$, 
${V}_{\mathrm{CNT}}=0.07$, ${\alpha}_{\mathrm{F}}=0.4$,
$(k_{1},k_{2})=(50,0.5)$}
\begin{tabular}{ccccccccc}
\thead
\xmorerows{1}{Beam type}  & \xmorerows{1}{$(\xi,\zeta)$}  &\multicolumn{1}{c}{} & \multicolumn{3}{c}{$M^*=4$} &  \multicolumn{3}{c}{$M^*=8$} \\\cline{4-9}
 &  & & $f_{\mathrm{p}}=0$ & $f_{\mathrm{p}}=0.1$ & $f_{\mathrm{p}}=0.2$ &   $f_{\mathrm{p}}=0$ & {$f_{\mathrm{p}}=0.1$} & $f_{\mathrm{p}}=0.2$ \\
\endthead
\morerows{5}{(2-1-2) }  & \morerows{2}{(0.1, 0.1)}  & $| u^{*}| $ & 0.4787 & 0.5638 & 0.6697 &   {1.1445} & 1.1928 & {1.1939} \\
 &  & $w^{*}$ & 0.6787 & 0.6994 & 0.7082 &   {0.4658} & 0.3426 & {0.1953} \\
 &  & ${\theta}^{*}$ & 0.2992 & 0.3317 & 0.3723 &   {0.6069} & 0.6754\vspace*{5pt} & {0.7618} \\
 & \morerows{2}{(0.4, 0.7)}  & $| u^{*}| $ & 0.5570 & 0.6506 & 0.7641 &   {1.1905} & 1.1982 & {1.1470} \\
 &  & $w^{*}$ & 0.6982 & 0.7079 & 0.7001 &   {0.3526} & 0.2205 & {0.0874} \\
 &  & ${\theta}^{*}$ & 0.3291 & 0.3649 & 0.4097 &   {0.6699} & 0.7459\vspace*{5pt} & {0.8419} \\
\morerows{5}{(1-1-2) }  & \morerows{2}{(0.1, 0.1)}  & $| u^{*}| $ & 1.1260 & 1.1707 & 1.1714 &   {0.8757} & 0.9295 & {1.0305} \\
 &  & $w^{*}$ & 0.4437 & 0.3256 & 0.1848 &   {0.1247} & 0.1896 & {0.1815} \\
 &  & ${\theta}^{*}$ & 0.6102 & 0.6782 & 0.7644 &   {1.2610} & 1.3954\vspace*{5pt} & {1.5612} \\
 & \morerows{2}{(0.4, 0.7)}  & $| u^{*}| $ & 1.1685 & 1.1755 & 1.1273 &   {0.9233} & 1.0137 & {1.0785} \\
 &  & $w^{*}$ & 0.3351 & 0.2089 & 0.0818 &   {0.1864} & 0.1905 & {0.1050} \\
 &  & ${\theta}^{*}$ & 0.6727 & 0.7486 & 0.8447 &   {1.3847} & 1.5310 & {1.7135} 
\botline
\end{tabular}
\end{table}

\begin{figure}
\includegraphics{fig13}
\caption{\label{fig13}The tip response and deformed configurations of
CNTRC (2-1-2) beam according to ${M}^{*}=8$ for $L/h=10$,
${V}_{\mathrm{CNT}}=0.07$, ${\xi}=0.4$, $\zeta=0.7$, ${\alpha}_{\mathrm{F}}=0.4$,
$(k_{1},k_{2})=(50,0.5)$ and different values of $f_{\mathrm{p}}$.}
\end{figure}

\section{Conclusions}

The large deflections of the CNTRC sandwich cantilever beam partially
embedded on a Pasternak foundation under the end forces have been
studied considering the influence of CNTs agglomeration. The sandwich
beam is composed of homogeneous and two face layers made from a
CNT-reinforced composite. Eshelby--Mori--Tanaka model has been applied to
predict the effective moduli of the composite. Based on the first-order
shear deformation theory, a total Lagrangian beam element was derived
and used to establish the equilibrium equation for the beam. The
Newton--Raphson method was used in conjunction with the arc-length
technique to compute the load-displacement curves. The effects of
various important factors, including the CNT volume fraction, CNT
agglomeration parameters, porosities, foundation stiffness and
supporting length on the large deflection response of the beam have
been studied in detail. The finding shows that the sandwich cantilever
beam is more conservative for a higher CNT volume fraction and a longer
supporting part. The agglomeration of CNTs has an important role in the
behavior of the sandwich beam, and it not only makes the beam more
flexible but also raises the effect of the shear deformation on the
large deflections. It was also shown that the sandwich configuration is
important in the behavior of the CNTRC sandwich beam and the symmetric
beam tends to be more conservative than the non-symmetric one.

\section*{Declaration of interests}

The authors do not work for, advise, own shares in, or receive funds
from any organization that could benefit from this article, and have
declared no affiliations other than their research organizations.

\section*{Data availability statement}

The paper is self-contained and uses no external data.

\section*{Funding}

This research is funded by Vietnam National Foundation for Science and
Technology Development (NAFOSTED) under grant number 107.02-2021.11.

\CDRGrant[NAFOSTED]{107.02-2021.11}

\back{}

\appendix{}

\def\appendixlabel{}

\setcounter{equation}{0}
\renewcommand\theequation{A\arabic{equation}}

\section{}

\def\thesection{Appendix A}

This Appendix lists the coefficients of the internal force vectors and
tangent stiffness matrices in Equations~(\ref{eq24}), (\ref{eq26}) and
(\ref{eq27}).
\begin{enumerate}[(1)]
\item~Components of the internal force vectors in Equation~(\ref{eq24})
{\begin{eqnarray}
\begin{array}{l}
{f}_{u_{I}}^{\mathrm{B}}=-\left(A_{11}\hat{\varepsilon }+A_{12}\hat{\chi
}\right)\cos \hat{\theta }+\psi A_{33}\hat{\gamma }\sin \hat{\theta
},\quad {f}_{w_{I}}^{\mathrm{B}}=-\left(A_{11}\hat{\varepsilon }+A_{12}\hat{\chi
}\right)\sin \hat{\theta }-\psi A_{33}\hat{\gamma }\cos \hat{\theta }\vspace*{5pt}\\
f_{\theta _{I}}^{\mathrm{B}}=\displaystyle\frac{l}{2}A_{11}\hat{\varepsilon }\hat{\gamma
}+\displaystyle\frac{l}{2}A_{12}\left(\hat{\gamma }\hat{\chi
}-\displaystyle\frac{2}{l}\hat{\varepsilon }\right)-A_{22}\hat{\chi
}-\displaystyle\frac{l}{2}\psi A_{33}\left(\hat{\varepsilon }+1\right)\hat{\gamma}\vspace*{5pt}\\
{f}_{u_{J}}^{\mathrm{B}}=\left(A_{11}\hat{\varepsilon }+A_{12}\hat{\chi
}\right)\cos \hat{\theta }-\psi A_{33}\hat{\gamma }\sin \hat{\theta
},\quad {f}_{w_{J}}^{\mathrm{B}}=\left(A_{11}\hat{\varepsilon }+A_{12}\hat{\chi
}\right)\sin \hat{\theta }+\psi A_{33}\hat{\gamma }\cos \hat{\theta }\vspace*{5pt}\\
{f}_{\theta _{J}}^{\mathrm{B}}=\displaystyle\frac{l}{2}A_{11}\hat{\varepsilon }\hat{\gamma
}+A_{12}\left(\displaystyle\frac{l}{2}\hat{\gamma }^{2}-\hat{\varepsilon
}\right)+A_{22}\hat{\chi }-\frac{l}{2}\psi A_{33}\left(\hat{\varepsilon
}+1\right)\hat{\gamma }
\end{array}
\label{eqA1}
\end{eqnarray}}\unskip
\noindent and
{\begin{eqnarray}
\begin{array}{l}
{f}_{u_{I}}^{\mathrm{F}}=\displaystyle\frac{l}{6}k_{\mathrm{W}}\left(2u_{I}+u_{J}\right)+k_{G}\left(\hat{\theta }-
\hat{\gamma }\right)\sin \hat{\theta },\quad  {f}_{w_{I}}^{\mathrm{F}}=
\frac{l}{6}k_{\mathrm{W}}\left(2w_{I}+w_{J}\right)+k_{G}\left(\hat{\theta }
-\hat{\gamma }\right)\cos \hat{\theta }\vspace*{5pt}\\
{f}_{\theta _{I}}^{\mathrm{F}}=lk_{G}\left(\hat{\theta }-\hat{\gamma
}\right)\left(1-\displaystyle\frac{l}{2}\hat{\varepsilon
}\right),\quad {f}_{u_{J}}^{\mathrm{F}}=\frac{l}{6}k_{\mathrm{W}}
\left(u_{I}+2u_{J}\right)+k_{G}\left(\hat{\theta }-\hat{\gamma }\right)
\sin \hat{\theta }\vspace*{5pt}\\
{f}_{w_{J}}^{\mathrm{F}}=\displaystyle\frac{l}{6}k_{\mathrm{W}}\left(w_{I}+2w_{J}\right)-k_{G}
\left(\hat{\theta }-\hat{\gamma }\right)\cos \hat{\theta },\quad 
\mbox{and}\quad {f}_{\theta _{J}}^{\mathrm{F}}={f}_{\theta _{I}}^{F}
\end{array} \label{eqA2}
\end{eqnarray}}\unskip
\item Coefficients of the tangent stiffness matrices in Equations~(\ref{eq26})
and (\ref{eq27})
\end{enumerate}
\begin{itemize}
\item The diagonal coefficients of the matrix
${\mathbf{k}}_{t}^{B}$ in (\ref{eq26})
{\begin{eqnarray}
\begin{array}{l}
{k}_{u_{I}u_{I}}^{B}=\displaystyle\frac{1}{l}\left(A_{11}\cos ^{2}\hat{\theta }+\psi
A_{33}\sin ^{2}\hat{\theta
}\right),\quad {k}_{w_{I}w_{I}}^{B}=\displaystyle\frac{1}{l}\left(A_{11}\sin
^{2}\hat{\theta }+\psi A_{33}\cos ^{2}\hat{\theta }\right)\vspace*{5pt}\\
{k}_{\theta _{I}\theta _{I}}^{B}=\displaystyle\frac{lA_{11}}{4}\left[\hat{\gamma
}^{2}-\hat{\varepsilon }\left(\hat{\varepsilon
}+1\right)\right]-\displaystyle\frac{lA_{12}}{4}\left[\left(\hat{\varepsilon
}+1\right)\hat{\chi }+\displaystyle\frac{4}{l}\hat{\gamma
}\right]+\displaystyle\frac{A_{22}}{l}+\displaystyle\frac{l\psi
A_{33}}{4}\left[\left(\hat{\varepsilon }+1\right)^{2}-\gamma
^{2}\right]\vspace*{5pt}\\
{k}_{\theta _{J}\theta _{J}}^{B}=\displaystyle\frac{lA_{11}}{4}\left[\hat{\gamma
}^{2}-\hat{\varepsilon }\left(\hat{\varepsilon
}+1\right)\right]-\displaystyle\frac{lA_{12}}{4}\left(\hat{\varepsilon
}+1\right)\hat{\chi }+\frac{A_{22}}{l}+\frac{l\psi
A_{33}}{4}\left[\left(\hat{\varepsilon }+1\right)^{2}-\gamma
^{2}\right]
\end{array}
\label{eqA3}
\end{eqnarray}}\unskip
\noindent and
${k}_{u_{J}u_{J}}^{B}={k}_{u_{I}u_{I}}^{B},\,\,{k}_{w_{J}w_{J}}^{B}=
{k}_{w_{I}w_{I}}^{B}$.\vspace*{3pt}

\item The non-diagonal coefficients of the stiffness matrix
${\mathbf{k}}_{t}^{B}$ in (\ref{eq26})
{\begin{eqnarray}
\begin{array}{l}
{k}_{u_{I}w_{I}}^{B}=\displaystyle\frac{A_{11}}{l}\sin \hat{\theta }\cos \hat{\theta
}-\frac{\psi A_{33}}{l}\sin \hat{\theta }\cos \hat{\theta }\vspace*{5pt}\\
{k}_{u_{I}\theta _{I}}^{B}=\displaystyle\frac{A_{11}}{2}\left(\hat{\varepsilon }\sin
\hat{\theta }-\hat{\gamma }\cos \hat{\theta
}\right)-\frac{A_{12}}{2l}\left(\hat{\chi }\sin \hat{\theta }+2\cos
\hat{\theta }\right)+\displaystyle\frac{\psi A_{33}}{2}\left[\hat{\gamma }\cos
\hat{\theta }-\left(\hat{\varepsilon }+1\right)\sin \hat{\theta
}\right]\vspace*{5pt}\\
{k}_{w_{I}\theta _{I}}^{B}=-\displaystyle\frac{A_{11}}{2}\left(\hat{\gamma }\sin
\hat{\theta }+\hat{\varepsilon }\cos \hat{\theta
}\right)+\displaystyle\frac{A_{12}}{2l}\left(l\hat{\chi }\cos \hat{\theta }-2\sin
\hat{\theta }\right)+\displaystyle\frac{\psi A_{33}}{2}\left[\hat{\gamma }\sin
\hat{\theta }+\left(\hat{\varepsilon }+1\right)\cos \hat{\theta
}\right]\vspace*{5pt}\\
{k}_{w_{I}\theta _{J}}^{B}=-\displaystyle\frac{A_{11}}{2}\left(\hat{\gamma }\sin
\hat{\theta }+\hat{\varepsilon }\cos \hat{\theta
}\right)+\frac{A_{12}}{2l}\left(l\hat{\chi }\cos \hat{\theta }+2\sin
\hat{\theta }\right)+\displaystyle\frac{\psi A_{33}}{2}\left[\hat{\gamma }\sin
\hat{\theta }+\left(\hat{\varepsilon }+1\right)\cos \hat{\theta
}\right]\vspace*{5pt}\\
{k}_{u_{I}\theta _{J}}^{B}=\displaystyle\frac{A_{11}}{2}\left(\hat{\varepsilon }\sin
\hat{\theta }-\hat{\gamma }\cos \hat{\theta
}\right)+\displaystyle\frac{A_{12}}{2l}\left(\hat{\chi }\sin \hat{\theta }-2\cos
\hat{\theta }\right)+\displaystyle\frac{\psi A_{33}}{2}\left[\hat{\gamma }\cos
\hat{\theta }-\left(\hat{\varepsilon }+1\right)\sin \hat{\theta
}\right]
\end{array}
\label{eqA4}
\end{eqnarray}}\unskip
\noindent and
{\begin{eqnarray*}
\begin{array}{l}
{k}_{u_{I}u_{J}}^{B}=-{k}_{u_{I}u_{I}}^{B},\quad 
{k}_{u_{I}w_{J}}^{B}=-{k}_{u_{I}w_{I}}^{B},\quad
{k}_{u_{I}\theta _{J}}^{B}={k}_{u_{I}\theta _{I}}^{B},
\quad {k}_{w_{I}w_{J}}^{B}=-{k}_{w_{I}w_{I}}^{B},\vspace*{5pt}\\
{k}_{w_{I}\theta _{J}}^{B}=-{k}_{u_{I}\theta
_{I}}^{B},\quad {k}_{u_{J}w_{J}}^{B}={k}_{u_{I}w_{I}}^{B},
\quad {k}_{u_{J}\theta _{J}}^{B}=-{k}_{u_{I}\theta _{I}}^{B},
\quad {k}_{w_{J}\theta _{J}}^{B}=-{k}_{w_{I}\theta _{J}}^{B}
\end{array}
\end{eqnarray*}}\unskip
\item~The coefficients of the stiffness matrix
${\mathbf{k}}_{t}^{F}$ in (\ref{eq27})
{\begin{eqnarray}
\begin{array}{l}
{k}_{u_{I}u_{I}}^{F}=\displaystyle\frac{l}{3}k_{\mathrm{W}}+\displaystyle\frac{1}{l}\sin ^{2}\hat{\theta
}k_{G},\quad {k}_{w_{I}w_{I}}^{F}=\displaystyle\frac{l}{3}k_{\mathrm{W}}+\displaystyle\frac{1}{l}\cos
^{2}\hat{\theta }k_{G},\quad {k}_{\theta _{I}\theta
_{I}}^{F}=\displaystyle\frac{l}{4}\left[\left(2-\hat{\varepsilon
}\right)^{2}-\left(\hat{\theta }-\hat{\gamma }\right)\hat{\gamma
}\right]k_{G},\vspace*{5pt}\\
{k}_{u_{I}w_{I}}^{F}=-\displaystyle\frac{1}{2l}\sin 2\hat{\theta }k_{G}\vspace*{5pt}\\
{k}_{u_{I}\theta _{I}}^{F}=-\displaystyle\tfrac{1}{2}\left[\left(2-\hat{\varepsilon
}\right)\sin \hat{\theta }-\left(\hat{\theta }-\hat{\gamma }\right)\cos
\hat{\theta }\right]k_{G},\quad {k}_{w_{I}\theta
_{I}}^{F}=\displaystyle\tfrac{1}{2}\left[\left(2-\hat{\varepsilon }\right)\cos
\hat{\theta }-\left(\hat{\theta }-\hat{\gamma }\right)\sin \hat{\theta
}\right]k_{G}
\end{array}\label{eqA5}
\end{eqnarray}}\unskip
\noindent and
{\begin{eqnarray*}             
\begin{array}{l}
{k}_{u_{J}u_{J}}^{F}={k}_{u_{I}u_{I}}^{F},\quad
{k}_{u_{I}u_{J}}^{F}=-{k}_{u_{I}u_{J}}^{F},\quad 
{k}_{w_{J}w_{J}}^{F}={k}_{w_{I}w_{I}}^{F},\quad 
{k}_{w_{I}w_{J}}^{F}=-{k}_{w_{J}w_{J}}^{F},\quad 
{k}_{\theta _{I}\theta _{J}}^{F}={k}_{\theta _{J}
\theta _{J}}^{F}={k}_{\theta _{I}\theta _{I}}^{F},\vspace*{5pt}\\
{k}_{u_{J}w_{J}}^{F}={k}_{u_{I}w_{I}}^{F},\quad 
{k}_{u_{I}w_{J}}^{F}=-{k}_{u_{I}w_{I}}^{F},\quad 
{k}_{u_{I}\theta _{J}}^{F}={k}_{u_{I}\theta _{I}}^{F},\quad 
{k}_{w_{I}\theta _{J}}^{F}={k}_{w_{I}\theta _{I}}^{F},\quad 
{k}_{w_{J}\theta _{I}}^{F}={k}_{w_{J}\theta _{J}}^{F}=-{k}_{w_{I}
\theta _{J}}^{F}
\end{array}
\end{eqnarray*}}\unskip
\end{itemize}

\bibliographystyle{crunsrt}
\bibliography{crmeca20240330}
\refinput{crmeca20240330-reference.tex}

\end{document}
