In the literature, there are various simplifying assumptions adopted in the kinematic relations of the faces and the core when considering a geometrically nonlinear problem in sandwich structures. Most commonly, only one nonlinear term is included in the faces and the core nonlinearities are neglected. A critical assessment of these assumptions, as well as the effects of including the other nonlinear terms in the faces and the core, is the scope of this paper. The comprehensive investigation of all the nonlinear terms is accomplished by deriving and employing an advanced nonlinear high-order theory, namely, the recently developed “extended high-order sandwich panel theory” (EHSAPT). This theory, which was derived as a linear theory, is first formulated in this paper in its full nonlinear version for the simpler one-dimensional case of sandwich wide panels/beams. Large displacements and moderate rotations are taken into account in both faces and core. In addition, a nonlinear EHSAPT-based finite element (FE) is developed. A series of simplified models with various nonlinear terms included are derived accordingly to check the validity of each of these assumptions. Two sandwich panel configurations, one with a “soft” and one with a “hard” core, loaded in three-point bending, are analyzed. The geometric nonlinearity effects and the relative merits of the corresponding simplifications are analyzed with these two numerical examples. In addition to a relative comparison among all these different assumptions, the results are also compared to the corresponding ones from a commercial FE code.

## Introduction

Sandwich structures are widely used in the aerospace, marine, ground transportation, and civil industries due to their extraordinary mechanical properties combined with lightweight. A typical sandwich panel consists of two high stiffness thin faces, which are separated by a low density thick core. This configuration enables the sandwich structures to exhibit high bending stiffness with little resultant weight. Both low stiffness soft materials, such as nonmetallic honeycomb and polymeric foam or relatively higher stiffness metallic honeycomb or balsa wood, can be used as the core materials. In practice, the soft core sandwich structures are widely used due to their advantages of weight saving and efficient manufacturing processes compared to others with metallic core. Compared to other sandwich structures, the soft core makes sandwich structures to have several distinguish features, such as the transverse compressibility and the shear effect [1]. Due to these features, the sandwich structures are more sensitive to the localized effects, which always are a result of stress concentrations and large deformation gradients. Failures are most likely to happen in these areas. In general, these localized effects are associated with large deformations and moderate rotations. So, the geometric nonlinearities need to be considered and well understood when these problems are addressed.

Sandwich structures have been investigated by many researchers. Back to 1970, Pagano [2] gave the three-dimensional elasticity solution to the linear static behavior of sandwich plates, and in 2011, Kardomateas and Phan followed Pagano's work and extended the elastic solution to the negative discriminant case [3]. These elasticity analytical solutions serve as benchmarks but are only available for certain loading and boundary conditions and do not include nonlinearities. In addition, sandwich structural theories have also been developed. Plantema [4] and Allen [5] well summarized the work done in the 1960s. The earliest models are called the classical or first-order shear models. In these models, the core is assumed to be incompressible in its transverse direction and the in-plane rigidity is neglected, and only the shear resistance can be considered. These assumptions may be acceptable for sandwich structures with a high stiffness core, such as metallic core. But it has been shown to be inaccurate for sandwich structures with very soft core [3]. Later, models considering the flexible core were developed to overcome this limitation. In 1992, Frostig et al. [6] proposed the HSAPT, in which the transverse compressibility and the shear effects of the core sheet are included. However, the axial rigidity was neglected in this model. In 2012, Phan et al. [7] formulated the EHSAPT, in which the axial rigidity of core was included in addition to the transverse compressibility and shear effect. By comparison to the elasticity solution, the EHSAPT shows high accuracy in both displacement and stress distributions for a wide range of core materials [7].

The major part of the research works related to sandwich structures is performed under the framework of small displacement and small deformation. Indeed, the above models and theories were originally proposed and analyzed with the infinitesimal kinematic relation. This has enabled researchers to have a basic understanding of sandwich structures, but at the same time, it limits the further understanding involving phenomena that exhibit geometric nonlinearity. In 2003, Sokolinsky et al. [8] compared the linear HSAPT and nonlinear HSAPT with the experimental results from four-point bending tests of a sandwich beam. It was shown that the linear analysis is efficient in estimating the transverse displacements while the geometric nonlinearities should be considered when accurate predictions are required for the longitudinal displacements and peeling stresses.

The large deformation could be considered by imposing the geometric nonlinearities in the kinematic relations. The literature survey reveals that various simplifications are often made when including the geometric nonlinearities. These assumptions and simplifications are made mainly due to the complexity of the governing equations. Indeed, the complicated governing equations make it nearly impossible to consider the full nonlinear kinematic relations. Usually, in papers involving geometric nonlinearities, only some of the nonlinear terms appear in the kinematic relations. These simplifications are various in the literature, and no literature can be found that critically examines the effects and validations of these simplifications. This is actually one of the goals of this paper.

Most of the researches on this topic only consider large deformation in the face sheets, while small deformations and infinitesimal linear strain relations are used in the core. Frostig and Baruch's nonlinear HSAPT model [9] assumes the faces undergoing small displacements and large rotations while the core is assumed to deform with small displacements and rotations. Sokolinsky and Frostig [10,11] analyzed the nonlinear static behaviors of sandwich beams by following the nonlinear HSAPT. In 2008, Li and Kardomateas [12] presented a nonlinear higher-order theory for sandwich plates. The nonlinear kinematic relations were adopted in the faces while the core was considered to have large rotations with small displacements.

In some models, although the core uses large deformation and Green–Lagrange strain during the equation derivation, the simplifications are made to simplify the final equations and enable it to be solvable. In 2003, Hohe and Librescu [13] presented a nonlinear theory for double-curved sandwich shells with compressible core. In both faces and core, the large deformation is expressed by the Green–Lagrange strain tensor, and then, only the nonlinear strain terms with respect to the transverse displacements are kept while the other nonlinear terms are discarded. In 2005, Frostig et al. [14] formulated the nonlinear field equations which use the nonlinear kinematic relations of large displacement and moderate rotation in both faces and core, followed by two simplified models. The simplifications are all made in the core sheet. One model only includes the nonlinear kinematic relations of the shear angle while the other one assumes that the linear kinematic relations are still valid. These two models are the models used for the analysis eventually. The latter one coincides with the nonlinear HSAPT. In these two models, which can be thought of as an extension of the HSAPT, the axial rigidity of the core sheet is still neglected. In 2014, Dariushi and Sadighi [15] used the Green–Lagrange strain in both faces and core. Since the classical beam theory was adopted in the faces, the face sheets have nonzero nonlinear transverse normal strain component and shear strain component while the corresponding linear parts are equal to zero.

When geometric nonlinearities are considered, the analytical solutions are usually not available, and the numerical results are also not easy to obtain. Hence, a stable numerical approach is needed. Two methods have been traditionally used in this regard: one is the finite difference method [14], and the other is the FE approach, the latter being more commonly used in solid mechanics due to its advantages of easily handling the boundary conditions and presenting clear physical meaning and standardized formulation. Several special FEs for sandwich structures have been reported for linear or nonlinear analysis. In general, two kinds of such FEs have been proposed: One is the layer-by-layer element [16,17], which is based on the “zigzag” theory. The displacement has a piecewise variation through the thickness. Multiple elements through the thickness direction are need for laminates and sandwich composite plates. The other kind of FE is based on the sandwich theory [1821]. Only one element is required through the thickness direction, and it contains the information of all sheets. Hu et al. [18] used a one-dimensional element with the nonlinear kinematic equation in the skin and linear relation in the core to get the global and local buckling of sandwich beams. Several elements were formulated based on the high-order theory for the sandwich plate. The linear dynamic and linear static analyses were performed successfully [19,20]. Based on the EHSAPT, one novel element was recently proposed by Yuan et al. [21] in 2015. The linear static behavior of sandwich beam/wide panel was analyzed. It gives identical displacement and stress distribution as the elasticity solution.

This paper presents the nonlinear formulation of EHSAPT, in which large displacements and moderate rotations are considered in both faces and core. To assess the relative merits of the commonly used simplifications in the literature, various simplified models using different combinations of nonlinear terms in the kinematic relations will be considered. The EHSAPT-based element is further developed to include the nonlinear Green–Lagrange strain expressions, and it is then used to obtain the numerical results. A commercial finite element analysis (FEA) software is also used for comparison. Conclusions are finally drawn based on the results.

## Theory and Derivation

The nonlinear formulation is an extension of the linear formulation of the EHSAPT [7]. The faces are considered as Euler–Bernoulli beams with bending and axial rigidities while the core has a high-order displacement pattern.

Consider a sandwich panel of length a with a core of thickness 2c and top and bottom face sheets of thicknesses ft and fb, respectively. For convenience, a right-hand Cartesian coordinate system is set with the origin placed at the left edge of the beam. The x axis coincides with the middle line of the core, and the z axis is along the thickness direction, as shown in Fig. 1. We only consider loading and deformation in x–z plane, and we denote the displacement components in x and z coordinates as u and w, respectively.

Following the Euler–Bernoulli assumptions, the displacement field in the top face sheet $(c is expressed as
$ut(x,z)=u0t(x)−(z−c−ft2)w0,xt(x); wt(x,z)=w0t(x)$
(1)
and in the bottom face sheet $(−c−fb≤z<−c)$, it is
$ub(x,z)=u0b(x)−(z+c+fb2)w0,xt(x); wb(x,z)=w0b(x)$
(2)

where $u0t,b$ and $w0t,b$ are the axial and transverse displacements of the centroid of the top and bottom faces.

In the core sheet ($−c≤z≤c$), the axial displacement $uc(x,z)$ is a cubic polynomial in the coordinate z
$uc(x,z)=z22c2(1+zc)u0t(x)+ftz24c2(1+zc)w0,xt(x)+(1−z2c2)u0c(x)+z(1−z2c2)ϕ0c(x)+z22c2(1−zc)u0b(x)+fbz24c2(−1+zc)w0,xb(x)$
(3)
and the transverse displacement w(x, z) is assumed to be a quadratic function of z
$wc(x,z)=(z2c+z22c2)w0t(x)+(1−z2c2)w0c(x)+(−z2c+z22c2)w0b(x)$
(4)

where $u0c$ and $w0c$ are the axial displacement and the transverse displacement of centerline of the core sheet. $ϕ0c$ is the slop of the centroid of the core.

It is easy to prove that the displacement continuity conditions at the interfaces of the core and the two face sheets ($z=±c$) are fulfilled by the given displacement fields. Thus, the displacement field of the sandwich panel can be expressed by seven unknowns, $u0t(x), w0t(x), u0b(x), w0b(x), u0c(x), w0c(x)$, and $ϕ0c(x)$, which are functions of the x coordinate.

During the deformation, the materials are assumed within the linear elastic region and only geometric nonlinearities are considered. All sheets are considered undergoing large displacement and moderate rotation. If the Green–Lagrange strain is used directly, and recalling the displacement assumptions of the face sheet, listed by Eqs. (1) and (2), the strain components of the face sheets are
$ϵxxt,b(x,z)=∂ut,b(x,z)∂x+12[∂ut,b(x,z)∂x]2+12[∂wt,b(x,z)∂x]2$
(5a)

$ϵzzt,b(x,z)=12[∂ut,b(x,z)∂z]2$
(5b)

$γxzt,b(x,z)=∂ut,b(x,z)∂x∂ut,b(x,z)∂z$
(5c)
Compared to the general Green–Lagrange strain, some terms are neglected due to the zero value as a result of Eqs. (1) and (2). Notice that the Green–Lagrange strain results into nonzero transverse normal strain and shear strain in the face sheets. However, when an Euler–Bernoulli beam is considered, we usually consider only the axial normal strain ϵxx regardless of the analysis being linear or nonlinear. In the literature, some researchers include these two strain components ϵzz and γxz with zero linear part and nonzero nonlinear part in the face sheets, e.g., in Ref. [15]. The effects of these two additional strain components would be discussed later. For convenience, the strain components of the face sheets are written as
$ϵxxt,b(x,z)=∂ut,b(x,z)∂x+α12[∂ut,b(x,z)∂x]2+α22[∂wt,b(x,z)∂x]2$
(6a)

$ϵzzt,b(x,z)=α32[∂ut,b(x,z)∂z]2$
(6b)

$γxzt,b(x,z)=α4∂ut,b(x,z)∂x∂ut,b(x,z)∂z$
(6c)

where the values of αi can be 0 or 1. These coefficients are in the nonlinear terms to control whether the corresponding term is included or neglected.

Similarly, the strain components in the core sheet are written as follows:
$ϵxxc(x,z)=∂uc(x,z)∂x+β12[∂uc(x,z)∂x]2+β22[∂wc(x,z)∂x]2$
(7a)

$ϵzzc(x,z)=∂wc(x,z)∂z+β32[∂uc(x,z)∂z]2+β42[∂wc(x,z)∂z]2$
(7b)

$γxzc(x,z)=∂uc(x,z)∂z+∂wc(x,z)∂x+β5∂uc(x,z)∂x∂uc(x,z)∂z+β6∂wc(x,z)∂x∂wc(x,z)∂z$
(7c)

Again, coefficients $βj(j=1,...,6)$ are used, and their values are equal to either 0 or 1. When $βj=1$, the corresponding term is included in the kinematic relation and vice versa.

Hence, by controlling the different combinations of αi and βj, the kinematic relations can be reduced to the various simplified models. Some particular cases are listed as follows:

1. (1)

$αi=βj=0$: linear ESHAPT [7,21]

2. (2)

$α2=1, αi=βj=0, (i=2, j=1, ..., 6)$: only faces nonlinear

3. (3)

$α2=1, αi=βj=0, (i=2, j=1, ..., 6)$, and $E1c=0$: nonlinear HSAPT [8,10,11,14]

4. (4)

$β2=β3=1, βj=0, (j=1,4,5,6)$: neglecting the higher-order nonlinear terms of the core, e.g., the $[u,xc(x,z)]2$ in comparison to the $u,xc(x,z)$ in the $ϵxxc$, etc.

The kinematic relations are given by Eqs. (6) and (7). For convenience and consistency, the strain components of the faces and the core are expressed as a vector ${ϵt,b,c}=[ϵxxt,b,c ϵzzt,b,c γxzt,b,c]T$. The strain vector can be further divided into a linear part and a nonlinear part. Considering the EHSAPT displacement assumptions (Eqs. (1)(4)), the general strain expressions of the faces and the core are represented by the displacement vector, ${U¯}=[u0t w0t u0b w0b u0c w0c ϕ0c]T$. By collecting terms, it takes the following form:
${ϵk}=[LLk]{U¯}+12[LNLk({U¯})]{U¯} k=t,b,c$
(8)

where the first term of the right-hand side stands for the linear part of the strain vector, and $[LLt,b,c]$ is the same as the one used in the linear EHSAPT. The second term gives the nonlinear strain part, and $[LNLk({U¯})]$ is a function of the displacement vector ${U¯}$ and depends on the nonlinear coefficients αi and βj. It will be denoted as $[LNLk]$ for short.

Due to the square terms of the Lagrange strain, the variation form of strain vector is given as
$δ{ϵk}=([LLk]+[LNLk])δ{U¯} k=t,b,c$
(9)
The stress vector in the three sheets is ${σt,b,c}=[σxxt,b,c σzzt,b,c τxzt,b,c]T$. For sandwich panels made out of an orthotropic core, the stress in the core can be obtained from
${σxxcσzzcτxzc}=[C11cC13c0C13cC33c000C55c]{ϵxxcϵzzcγxzc}$
(10)
where $Cijc$ are the stiffness constants of the core. According to the two-dimensional elasticity, the material constants of the core are given by [1]
$C11c=E1cE3cE3c−E1cν31c2; C13c=ν31cE1cE3cE3c−E1cν31c2; C33c=E3c2E3c−E1cν31c2; C55c=G31c$

In the face sheets, when $α3=α4=0$, the stiffness constants of the face sheets are simply $C11t,b=E1t,b$, and the stresses are $σxxt,b=C11t,bϵxxt,b$. When $α3=1$, or $α4=1, Cijt,b$ are obtained in the same way as the core.

## FE Formulation

Earlier, a special FE was formulated based on the EHSAPT [21], and the linear static problem of sandwich panels was solved. It yielded identical results compared to the elasticity solutions. In this part, the EHSAPT-based FE will be further developed to include the nonlinear kinematic relations. This will be used to obtain numerical results that demonstrate the relative merits of the various nonlinear terms.

The EHSAPT-based element can have two nodes or more than two nodes [21]. The two-node element has shown very high accuracy in the linear static cases [21], and it will be used here for the nonlinear cases due to its simplicity. This element only requires a one-dimensional FE model based on the sandwich beam/panel, as shown in Fig. 2. The sandwich beam is discretized by m elements. Hence, there are $n=m+1$ nodes totally for a two-node element model. The modeling process is much easier when compared to the general FE analysis, where a two-dimensional FE model is required for a sandwich beam/wide panel. Each node of this EHSAPT element has ten degrees-of-freedom (DOFs), and these DOFs are placed through the thickness direction. One single node contains the displacement information of all sheets. The displacement vector is listed as

${Ui}=[uitwit(dwdxt)iuibwib(dwdxb)iuicwic(dwdxc)iϕic]T$
(11)

where $wt,b,c$ stand for the transverse displacements of the midline of the top face, bottom face, and core, and $dwt,b,c/dx$ are their first-order derivatives with respect to x. Also, $ut,b,c$ represent the axial displacements of the midline of three sheets, and $ϕc$ is the rotation angle of the core's centerline. The subscript i means that these unknowns are located at ith node.

The displacement field in each element is marked with tilde. Lagrange interpolation polynomials are used to interpolate the axial displacement field, $ũt,b,c$, and the rotation angle of the core, $ϕ̃c$, within each element. Hermite interpolating polynomials are used to interpolate the transverse displacement field $w̃t,b,c$. So, $ũt,b,c$ and $ϕ̃c$ have C0 continuity between adjacent elements, and $w̃t,b,c$ have C1 continuity. Details about this element can be found in Ref. [21]. By introducing the local coordinate $s=x−xi$, the displacement field within ith element is represented as
${Ũ(s)}=[N(s)]{Ue}$
(12)
where
${Ũ(s)}=[ũt(s)w̃t(s)ũb(s)w̃b(s)ũc(s)w̃c(s)ϕ̃c(s)]T$
and ${Ue}=[{Ui} {Ui+1}]T$ is a 20 × 1 column matrix representing the unknown displacement vector of the ith element. Also, $[N(s)]$ is the 7 × 20 displacement interpolation matrix (see the Appendix).
Using the element displacement vector ${Ũ(s)}$ of Eq. (12) instead of the midline displacement ${U¯}$, which appears in Eqs. (1)(4), and by substituting into Eq. (8), the strain is expressed in terms of the element node displacements, as follows:
${ϵk}=[BLk]{Ue}+12[BNLk]{Ue} k=t,b,c$
(13)

where $[BLk]=[LLk][N(s)]$ and $[BNLk]=[LNLk][N(s)]$. Notice that $[BLk]$ is the same as the one used in linear analysis and is independent of ${Ue}$. The second part stands for the nonlinear part ${ϵNL}$ and $[BNLk]$ depends on ${Ue}$.

Similarly, the variation is
$δ{ϵk}=([BLk]+[BNLk])δ{Ue} k=t,b,c$
(14)
The nonlinear element stiffness matrix comes from the minimization of the total element potential energy that consists of the strain energy of top face, core, and bottom face and the work of external load; thus, one has
$δΠe=δUt+δUc+δUb+δVe=0$
(15)
where $δUk$$(k=t,c,b)$ are the variations of strain energy, and $δVe$ is the variation of external potential energy associated with applied loads. In terms of strains and stresses, the variation of element total potential energy is
$δΠe=∫0hi∫cc+ftδ{ϵt}T{σt}dzds+∫0hi∫−ccδ{ϵc}T{σc}dzds+∫0hi∫−c−fb−cδ{ϵb}T{σb}dzds−δ{Ue}T{Re}$
(16)

where ${Re}$ is the nodal equivalent force.

Notice that, using Eqs. (13) and (14), and collecting the linear part, the variation of strain energy of each sheet is written as
$δUk=δ{Ue}T([KeLk]+[KeNLk]){Ue} k=t,b,c$
(17)
where $[KeLk]$ is the linear part of the element stiffness contributed by different sheets and independent of ${Ue}$, and $[KeNLk]$ is the nonlinear part stiffness matrix due to the geometric nonlinearities. These can be represented in the following form:
$[KeLk]=∬vk[BLk]T[D][BLk]dzds$

$[KeNLk]=∬vk(12[BLk]T[D][BNLk]+12[BNLk]T[D][BNLk]+[BNLk]T[D][BLk])dzds k=t,b,c$

where vk stand for the space domain of the top face, the bottom face, or the core within the element.

By substituting into Eq. (16), it leads to the element equilibrium equation
$([KeL]+[KeNL]){Ue}={Re}$
(18)
where
$[KeL]=[KeLt]+[KeLc]+[KeLb]$
and
$[KeNL]=[KeNLt]+[KeNLc]+[KeNLb]$

The secant stiffness contains a linear part and a nonlinear part. Both face sheets and the core contribute to these two parts. When all the nonlinear coefficients αi and βj are equal to zero, $[KeNL]$ is eliminated. In this case, only the linear part is left, and it is exactly the same as the one used in the linear analysis [21].

The stiffness matrix of the structure is obtained by assembling the element stiffness matrices. The arc-length method [22,23] is used to solve this nonlinear problem and to track the equilibrium path.

## Numerical Examples and Discussion of the Results

In this part, two sandwich panels loaded in three-point bending will be analyzed with the EHSAPT-based element including the geometric nonlinearities. Two different core materials will be used, in order to investigate the relative merits of various simplifications on the nonlinear kinematic relations. The results will also be compared to the numerical results obtained by the commercial FE code adina.

The geometry follows the sandwich panel used in Ref. [14]. The length of the panel is a = 300 mm. It contains two face sheets with same thickness, $ft=fb=0.50$ mm, and a core of thickness $2c=19.05$ mm. The top and bottom face sheets are made out of Kevlar with equivalent modulus of elasticity of 27.4 GPa. The core is made out of lightweight foam, Rohacell 50, with extensional modulus equal to 52.5 MPa and shear modulus $G31c=21.0$ MPa or Rohacell 200 WF with extensional modulus 350.0 MPa and shear modulus 150.0 MPa. A concentrated load is applied at the middle, and the beam/wide panel is loaded in three-point bending. The displacement constraint is applied at the lower surface, i.e., the edges of the bottom sheet. Since the geometry and loading condition exhibit symmetry about the midspan, only the left half-beam is analyzed. Together with the symmetry conditions, the boundary conditions applied to the EHSAPT-based FE model are
${w1b=0 (first node)unt=unb=unc=0 (last node)(dwtdx)n=(dwbdx)n=(dwcdx)n=ϕnc=0 (last node)$
(19)

To ensure accuracy, 340 EHSAPT-based two-node elements are used to build this half-beam model, that is, m = 340 and $n=m+1=341$. Furthermore, the elements are nonuniformly distributed and placed with a bias toward the edges and the center of the beam, where load concentrations occur.

The commercial fea software adina is used to compare the results. It is well known that numerous issues may happen when using commercial software, such as numerical difficulties and poor convergence. Thus, extra attention is needed to build the nonlinear FE model. A two-dimensional model with very fine mesh is built, as shown in Fig. 3. The nine-node plane stress element is used to model the face sheets and core. Through the thickness direction, four elements are used in the face sheets, and 15 elements are placed in the core. There are 19,787 nodes totally in the model. The large displacement and small strain analysis control option are chosen to consider the geometric nonlinearities.

### Example 1: Soft Core Sandwich Panel.

The first example to consider is a sandwich beam with a low density very soft core, Rohacell 50 foam. Various simplified models, which are reflected by the different combinations of coefficients αi and βj, are used. The corresponding applied load versus transverse displacement curves at the middle point of top face are given in Fig. 4. The middle point is where the load is applied and where the maximum transverse displacement is reached. These curves are denoted as curves 0–12 for convenience. Curve 0 is the result based on the linear version of EHSAPT, neglecting all nonlinear terms. The linear kinematic relation, infinitesimal strain, is used, and the displacement is proportional to the applied load. The result obtained by the commercial FE code adina is included in the figure as curve 12. Due to numerical issues, the adina cannot yield further results for displacements larger than 14 mm. The curve goes back because the elements are distorted. The element edges would penetrate each other under higher load. When changing the mesh size and the load step size, this issue cannot be overcome. This numerical issue always happens around that point, and it is interesting to notice that models with coarse elements can indeed go a little bit further. The other curves are results using EHSAPT with different coefficients αi and βj. Only the coefficients, whose value is equal to 1, are noted, and all others are assumed to be equal to zero. With the exception of the linear result, curve 0, all other curves exhibit nonlinear response when the load is higher than 500 N. When the load is under 500 N, and the maximum transverse displacement smaller than 3.5 mm, all curves agree well and are almost identical to the linear result; this is the region where the infinitesimal strain relation is valid. Under higher applied load, notable discrepancies among the different simplifications are observed.

Curves 1–4 consider the transverse normal strain $ϵzzt,b$ or/and transverse shear strain $γxzt,b$ in the faces. These two strain components are the consequence of using the Green–Lagrange strain directly under the assumed displacement field of the faces. These three curves show a distinctly different nonlinear behavior when compared to the other results, including adina. According to these curves, the sandwich panel becomes stiffer when the load is increased, which is unreasonable. This unusual result is because the nonzero nonlinear terms of transverse normal stain and transverse shear strain cause the sandwich panel to be “locked.” Recall that the Euler–Bernoulli assumptions used in the face sheets imply that the transverse rigidity and shear modulus are large enough that there is no transverse normal strain and shear strain in the face sheets (these assumptions are expressed in the displacement field). Any nonzero transverse normal strain and shear would cause a large stress component that would make the beam to be locked. Hence, only the axial normal strain should be considered in the face sheets. There is no problem in the linear analysis since $ϵzzt,b$ and $γxzt,b$ vanish automatically. However, when considering large displacements and using the Green–Lagrange strain in the faces, α3 and α4 should be zero in order to be in agreement with the Euler–Bernoulli assumptions. In addition, α1 can also be zero since $(u,xt,b)2$ usually can be neglected compared to $(w,xt,b)2$, which can be seen from curves 1 and 2. So, the strain that needs to be considered in the faces should the same as the ordinary Euler–Bernoulli beam, i.e., only $α2=1$.

Curves 5 and 6 are the results when the nonlinear kinematic relation is only used in the faces while the core is assumed to retain the linear kinematic relation. This is the most commonly used assumption in sandwich panel nonlinear theories that consider the geometric nonlinearities, e.g., Refs. [1012,14]. These two curves exhibit the same softening tendency as adina and the other models. Curve 6 neglects the axial rigidity of the core by setting $E1c=52.5×10−5$ MPa, as an approximation to zero. This degenerates the nonlinear EHSAPT to the nonlinear HSAPT. As expected, it yields similar results as using nonlinear HSAPT directly [14]. A limit point is observed in curve 6 around 1438 N while it is reported as 1423 N in Ref. [14]. Since the axial rigidity is neglected, curve 6 is the lowest curve among these results, even in the linear region.

Curves 7–11 list several representative simplified models that include full or part of the nonlinear kinematic relation in the core. Curve 7 only includes the same nonlinear term as the faces, $(w,xt,b)2$, and the other terms are neglected. Curve 8 neglects the higher-order nonlinear terms (i.e., it neglects $(u,xc)2$ in comparison to $u,xc$ in the $ϵxxc, (w,zc)2$ in comparison to $w,zc$ in the $ϵzzc, u,xcu,zc$ in comparison to $u,zc$ in $γxzc$, and $w,xcw,zc$ in comparison to $w,xc$ in $γxzc$). These two yield almost same results. Curve 9 takes the full nonlinear expressions of $ϵxxc$ and $ϵzzc$ and only the linear part in the shear strain $γxzc$. Curve 11 includes all the nonlinear terms in the core while curve 10 only includes the nonlinear terms related to the transverse displacement. These two are very close to each other. Among these results, curves 10 and 11 are the two results that are the closest to adina. To this point, one may draw the conclusion that the nonlinear strain terms of the core have a great effect and cannot be neglected. Also, it can be concluded that the nonlinear effects caused by the transverse deflection have a more pronounced effect. It can also be seen that the EHSAPT-based element can get the whole response curves without any numerical difficulties while the adina fails to converge at some point in the nonlinear region.

A comparison between the nonlinear EHSAPT and the linear EHSAPT appears in Fig. 5. It plots the axial and transverse displacements along the axial direction at load P = 1300 N. The nonlinear EHSAPT refers to the case of $α2=βj=1,(j=1,...,6)$ and is plotted with line. The linear EHSAPT is plotted with markers. The linear EHSAPT and nonlinear analysis give similar transverse displacement distribution while the axial displacements are completely different. The linear analysis predicts the axial displacement of the three sheets as a symmetric pattern; the top face and bottom face rotate with respect to the core. But when the geometric nonlinearities are included, the edges would move toward to the center. Although the transverse deflections along the x direction are similar in this analysis, it can be seen that the localized effect (local dimpling at the load application) can only be captured when considering the large displacement and nonlinear kinematic relation. Thus, it is necessary to include the geometric nonlinearities in order to obtain accurate displacement results that include the local dimpling.

Figures 6 and 7 give the through-thickness strain and stress distributions of the core at $P=1300 N$. Both EHSAPT ($α2=βj=1, j=1,...,6$) and adina results are given. Green–Lagrange strains are used by both EHSAPT and adina. The conjugate stress is the second Piola–Kirchhoff stress. Since the direct stress output of adina is the Cauchy stress, the stress results denoted as adina are the corresponding second Piola–Kirchhoff stress results obtained based on the constitutive law and the adina strain result. The linear EHSAPT results are also plotted. Figure 6 is the cross section located at $x=0.5a$ (load application), and Fig. 7 is at $x=0.4a$. Recall the displacement assumptions used in EHSAPT, namely, the axial displacement field of the core, uc, is assumed to be a cubic function of z, while the transverse displacement, wc, is a quadric polynomial, as shown in Eqs. (3) and (4). Hence, the linear part (major part) of the $ϵzzc$ can only be a linear function of z. The axial strain $ϵxxc$ and shear strain $γxzc$ along the z direction are cubic functions. So, the EHSAPT is not expected to yield the exactly same strain and stress distributions as adina in cases that the distribution exhibits a high nonlinearity but rather, an approximated result. For example, at $x=0.5a$, where the concentrated load is applied, the strains and stresses have extreme high values near load concentration. Hence, EHSAPT gives a reasonable approximation to the adina results, as shown in Fig. 6. It should be pointed out that only one element is used by EHSAPT through the thickness, while 15 elements are used in the adina model through the thickness. Figure 7 gives the strain and stress results at $x=0.4a$$(x=120$ mm). At this location, the axial normal strain and stress are exactly the same as adina. Also, a good approximated result is obtained for the shear strain and stress. The transverse normal stresses are not the same, due to the linear transverse strain implied by the EHSAPT displacement field assumption, plotted in Fig. 7(b). But the EHSAPT gives the same stress magnitude level. This is an unavoidable limitation of EHSAPT. Both Figs. 6 and 7 show that the linear EHSAPT is not adequate to predict the strains and stresses when $P=1300N$, and that the geometric nonlinearities have a very noticeable effect on the strain and stress distributions.

### Example 2: Moderate Core Sandwich Panel.

In this example, the geometry, loading conditions, and boundary conditions are the same as in the previous example. Kevlar is used to construct the faces and similar foam, but with higher stiffness, Rohacell 200 WF is used as the core material. Rohacell 200 WF has a higher density with higher elastic modulus, and its elastic modulus is 350.0 MPa and its shear modulus is 150.0 MPa. Thus, it is about seven times stiffer than Rohacell 50, which was used in the first example.

The transverse deflection at midspan of the top face versus the applied load curves are given in Fig. 9. Several selected simplified cases are included. Curve 0 is the linear EHSAPT result. Curve 5 is the result from adina. Curves 1–4 are the nonlinear EHSAPT results. In these curves, the face sheet only considers one nonlinear term, the $(w,xt,b)2$. Due to the stiffer material used in the core, the response predicted by adina is closer to the linear result (curves up to a deflection value of 20 mm). It should be noticed that the commonly used assumption that the face sheet undergoes large displacement only and the core employs the linear kinematic relation leads to a stiffer nonlinear response, which is unrealistic. Indeed, when the nonlinear kinematic relation is only included in the faces (curve 1), the slope becomes larger and the curve 1 is above the linear result. Curve 2 represents the case that only uses the nonlinear kinematic relation in the faces, but the axial rigidity of the core is neglected. It degenerates the EHSAPT to the nonlinear HSAPT. Curve 2 lies on the lower side of the linear results, but it is much lower than the adina results, even in the small deflection part. Curves 3 and 4 are the results when using the nonlinear kinematic relation in both faces and core. Again, they yield the closest result to the commercial software adina. Even when only one nonlinear term $(w,xc)2$ is considered in the core (in addition to the faces), the result is much improved, as shown as curve 3. This example shows that it is not sufficient to only consider the large displacement and geometric nonlinearities in the faces, especially when the core has high rigidity. In that case, the results may even be worse than the linear analysis result.

## Conclusions

The geometric nonlinear effects in sandwich panels are analyzed based on the EHSAPT. Large displacements with moderate rotations are considered in both faces and core, and the Green–Lagrange strain is used in all layers. In addition, the previously developed EHSAPT-based element (which was used to solve the linear static problem) was further developed to include the geometric nonlinearity effects. The nonlinear response was successfully obtained with the arc-length method. In the literature, various simplifications have been made when considering the large deformation of sandwich structures. These simplifications are critically assessed in this paper. Thus, a series of simplified models are derived, analyzed, and compared to verify the validity of these assumptions.

Two three-point bending examples are analyzed. These two examples use two similar sandwich panels with the same geometric properties, loading conditions, and boundary conditions. The only difference is the core material. The soft material and high density moderate materials are used in the core, respectively. For verification, the results are compared with the results obtained by the commercial fea software adina. The applied load versus midspan transverse displacement is used to analyze the nonlinear static response. The numerical examples reveal that these simplified kinematic relations yield similar results under small deformations, where the response is almost linear. When the sandwich panel exhibits nonlinear behavior, divergences are observed. It is concluded that in the faces, only the axial normal strain nonlinear term should be included. The transverse normal strain and shear strain nonlinear terms of the faces should not be included because they come from the nonlinear terms of the Green–Lagrange strain and violate the Euler–Bernoulli assumption adopted in the faces by EHSAPT. When the axial rigidity of the core is considered, as in the EHSAPT, the geometric nonlinearities of the core are important and cannot be ignored, especially when the core is of higher stiffness. In this case, considering that only the faces undergo large deformation can lead to a stiffer response, which is unrealistic. For other theories that neglect the axial rigidity of the core, such as the HSAPT, considering the large displacement only in the faces can still get the correct response tendency. But the response curve is lower, and the deformation is larger than it is expected, even at the small deformation region. In sum, the nonlinear EHSAPT, including the full Green–Lagrange strain in the core sheet and nonlinear axial normal strain in the faces, yields the closest results to adina, for both soft core material and harder core material. In the example with soft core, it successfully captures the limit point, whereas the top face buckles and cannot take more axial load.

The EHSAPT-based element is successful in performing the static analysis considering the geometric nonlinearities. The one-dimensional model used by this element shows high efficiency when compared to adina or similar commercial fea software. A two-dimensional model needs to be built when using adina, and special care is needed to avoid the numerical difficulties. When using adina, the numerical difficulties are observed in this paper and in other literatures [14]. The proposed EHSAPT-based element can yield the entire nonlinear static response curve.

## Acknowledgment

The financial support of the Office of Naval Research, Grant No. N00014-11-1-0597, and the interest and encouragement of the Grant Monitor, Dr. Y.D.S. Rajapakse, are both gratefully acknowledged. The authors would also like to acknowledge the valuable input of Professor Izhak Sheinman of Technion Israel Institute of Technology regarding adina.

### Appendix: Displacement Interpolation Matrix

The displacement interpolation matrix of the two-node EHSAPT-based element is given as follows:
$[N(s)]=[N1000000000N20000000000N3N400000000N5N60000000000N1000000000N20000000000N3N400000000N5N60000000000N1000000000N20000000000N3N400000000N5N60000000000N1000000000N2]$
where $Nk(k=1, ...,6)$ are functions of the local coordinate, s, and the length of ith element is $hi=xi+1−xi$
$N1=1−shi ; N2=shi$

$N3=1−3s2hi2+2s3hi3 ; N4=s−2s2hi+s3hi2$

$N5=3s2hi2−2s3hi3 ; N6=−s2hi+s3hi2$

## References

References
1.
Carlsson
,
L.
, and
Kardomateas
,
G. A.
,
2011
,
Structural and Failure Mechanics of Sandwich Composites
,
Springer
,
New York
.
2.
Pagano
,
N.
,
1970
, “
Exact Solutions for Rectangular Bidirectional Composites and Sandwich Plates
,”
J. Compos. Mater.
,
4
(
1
), pp.
20
34
.
3.
Kardomateas
,
G. A.
,
2011
, “
Three Dimensional Elasticity Solution for Sandwich Beams/Wide Plates With Orthotropic Phases: The Negative Discriminant Case
,”
J. Sandwich Struct. Mater.
,
13
(
6
), pp.
641
661
.
4.
Plantema
,
F. J.
,
1966
,
Sandwich Construction
,
Wiley
,
New York
.
5.
Allen
,
H. G.
,
1969
,
Analysis and Design of Structural Sandwich Panels
,
Pergamon
,
Oxford, UK
.
6.
Frostig
,
Y.
,
Baruch
,
M.
,
Vilnay
,
O.
, and
Sheinman
,
I.
,
1992
, “
High-Order Theory for Sandwich-Beam Behavior With Transversely Flexible Core
,”
J. Eng. Mech.
,
118
(
5
), pp.
1026
1043
.
7.
Phan
,
C. N.
,
Frostig
,
Y.
, and
Kardomateas
,
G. A.
,
2012
, “
Analysis of Sandwich Panels With a Compliant Core and With In-Plane Rigidity-Extended High-Order Sandwich Panel Theory Versus Elasticity
,”
ASME J. Appl. Mech.
,
79
(
4
), p.
041001
.
8.
Sokolinsky
,
V. S.
,
Shen
,
H.
,
Vaikhanski
,
L.
, and
Nutt
,
S. R.
,
2003
, “
Experimental and Analytical Study of Nonlinear Bending Response of Sandwich Beams
,”
Compos. Struct.
,
60
(
2
), pp.
219
229
.
9.
Frostig
,
Y.
, and
Baruch
,
M.
,
1993
, “
High-Order Buckling Analysis of Sandwich Beams With Transversely Flexible Core
,”
J. Eng. Mech.
,
119
(
3
), pp.
476
495
.
10.
Sokolinsky
,
V.
, and
Frostig
,
Y.
,
1999
, “
Nonlinear Behavior of Sandwich Panels With a Transversely Flexible Core
,”
AIAA J.
,
37
(
11
), pp.
1474
1482
.
11.
Sokolinsky
,
V.
, and
Frostig
,
Y.
,
2000
, “
Branching Behavior in the Nonlinear Response of Sandwich Panels With a Transversely Flexible Core
,”
Int. J. Solids Struct.
,
37
(
40
), pp.
5745
5772
.
12.
Li
,
R.
, and
Kardomateas
,
G. A.
,
2008
, “
Nonlinear High-Order Core Theory for Sandwich Plates With Orthotropic Phases
,”
AIAA J.
,
46
(
11
), pp.
2926
2934
.
13.
Hohe
,
J.
, and
Librescu
,
L.
,
2003
, “
A Nonlinear Theory for Doubly Curved Anisotropic Sandwich Shells With Transversely Compressible Core
,”
Int. J. Solids Struct.
,
40
(
5
), pp.
1059
1088
.
14.
Frostig
,
Y.
,
Thomsen
,
O. T.
, and
Sheinman
,
I.
,
2005
, “
On the Non-Linear High-Order Theory of Unidirectional Sandwich Panels With a Transversely Flexible Core
,”
Int. J. Solids Struct.
,
42
(
5–6
), pp.
1443
1463
.
15.
Dariushi
,
S.
, and
,
M.
,
2014
, “
A New Nonlinear High Order Theory for Sandwich Beams: An Analytical and Experimental Investigation
,”
Compos. Struct.
,
108
(
0
), pp.
779
788
.
16.
Ganapathi
,
M.
,
Patel
,
B. P.
, and
Makhecha
,
D. P.
,
2004
, “
Nonlinear Dynamic Analysis of Thick Composite/Sandwich Laminates Using an Accurate Higher-Order Theory
,”
Composites, Part B
,
35
(
4
), pp.
345
355
.
17.
,
S.
, and
Singha
,
M. K.
,
2013
, “
Geometrically Nonlinear Finite Element Analysis of Sandwich Plates Using Normal Deformation Theory
,”
Compos. Struct.
,
97
(
0
), pp.
84
90
.
18.
Hu
,
H.
,
Belouettar
,
S.
,
Potier-Ferry
,
M.
, and
,
A.
,
2009
, “
A Novel Finite Element for Global and Local Buckling Analysis of Sandwich Beams
,”
Compos. Struct.
,
90
(
3
), pp.
270
278
.
19.
Elmalich
,
D.
, and
Rabinovitch
,
O.
,
2012
, “
A High-Order Finite Element for Dynamic Analysis of Soft-Core Sandwich Plates
,”
J. Sandwich Struct. Mater.
,
14
(
5
), pp.
525
555
.
20.
Oskooei
,
S.
, and
Hansen
,
J. S.
,
2000
, “
Higher-Order Finite Element for Sandwich Plates
,”
AIAA J.
,
38
(
3
), pp.
525
533
.
21.
Yuan
,
Z.
,
Kardomateas
,
G. A.
, and
Frostig
,
Y.
,
2015
, “
Finite Element Formulation Based on the Extended High-Order Sandwich Panel Theory
,”
AIAA J.
,
52
(
10
), pp.
3006
3015
.
22.
Crisfield
,
M. A.
,
1981
, “
A Fast Incremental/Iterative Solution Procedure That Handles Snap-Through
,”
Comput. Struct.
,
13
(
1–3
), pp.
55
62
.
23.
Crisfield
,
M. A.
,
1983
, “
An Arc-Length Method Including Line Searches and Accelerations
,”
Int. J. Numer. Methods Eng.
,
19
(
9
), pp.
1269
1289
.