Composite Structures 284 (2022) 115161

Contents lists available at ScienceDirect

Composite Structures
journal homepage: www.elsevier.com/locate/compstruct

Modeling and simulation of static and dynamic behavior in composite

sandwich plates with hourglass lattice cores based on reduced-order model
Shi Jie, Zhong Yifeng ∗, Liu Rong, Shi Zheng
School of Civil Engineering, Chongqing University, Chongqing, PR China
Key Laboratory of New Technology of Construction of Cities in Mountain Area (Chongqing University), Ministry of Education, Chongqing 400045, PR China

ARTICLE INFO ABSTRACT

Keywords: The composite sandwich plate with hourglass lattice cores (CSP-HLC) is a novel cellular structure that can
Variational asymptotic method increase the width-to-length ratio and reduce the inter-node spacing. However, its static and dynamic behaviors
Hourglass lattice core are difficult to analyze due to the complex microstructures. In this work, a computational homogenization
Composite sandwich plate
method is proposed for calculating the constitutive parameters. A reduced-order plate model is then derived
Multiscale analysis
through dimensional reduction from the three-dimensional orthotropic thermoelasticity framework. The
Free vibration
effectiveness and accuracy of the proposed model were verified by comparing with the static-displacement
and free-vibration results of 3D direct numerical simulations. The parameter analysis showed that the different
material and structural parameters of the hourglass lattice core had different effects on the equivalent
stiffnesses and natural frequencies of the CSP-HLC. Compared with the composite sandwich plate with
pyramidal lattice cores, the displacement of the CSP-HLC was smaller under the same load and boundary
conditions, and the natural frequencies were also smaller, which may have been because the additional
equivalent density of the CSP-HLC influenced the frequency more predominantly than the extra equivalent
stiffness.

1. Introduction promising new generation of advanced, lightweight, and ultra-strength

materials [4].
A lattice sandwich structure is a kind of biomimetic structure with In the past few years, various truss lattice sandwich plates have
a statically determinate/indeterminate porous ordered microstructure
been produced to develop high-specific-strength, lightweight plates for
designed by simulating the molecular lattice configuration [1–3]. For
multifunctional applications, such as pyramidal, tetrahedral, and three-
the same structure weight, the in-plane Young’s modulus of lattice
materials can be more than two orders of magnitude higher than those dimensional (3D) Kagome truss lattices [5–7]. Pyramid and tetrahedral
of honeycombs and other lightweight materials, and the out-of-plane truss lattices are the most popular, showing superior structural charac-
strength can be more than one order of magnitude higher than those teristics. However, the results have shown that the strengths of these
of honeycomb materials. At the same time, the unique microscopic lattice sandwich structures depend on the compressive strength of the
periodic three-dimensional grid structure allows for the application core and the element size (controlling the inter-node spacing), and their
of finite element technology to determine the optimal configuration improvement seems to be feasible. To overcome the weakness of the
design.
traditional truss lattice, a novel hourglass truss lattice was developed.
In addition to the carrying capacity, a lattice material itself is an
Compared with the pyramid truss lattice with a similar relative density,
excellent compact heat exchanger. A large hole rate and hole design
scheme increases the heat dissipation area and is conducive to convec- the hourglass truss lattice not only has a higher width-to-length ratio
tion. The combination of mechanics and heat exchange performance of its core struts, but it also has a smaller inter-node spacing, which is
will be helpful to determine the best heat dissipation structures of approximately half that of the pyramid truss lattice with the equivalent
aircraft. At the same time, if the transmission transistor network is em- height [8].
bedded in the face sheets of the lattice sandwich structure, the bending Due to the complex geometry of lattice sandwich structures, in
and torsional deformation of the macro-structure can be realized by most cases, analytical solutions are difficult to obtain, and the calcu-
locally driving some lattice micro-rods with minimal energy. There-
lation methods become expensive when all of the geometric details
fore, materials with lattice sandwich structures have become the most

∗ Corresponding author.
E-mail address: [email protected] (Y. Zhong).

https://doi.org/10.1016/j.compstruct.2021.115161
Received 24 October 2021; Received in revised form 10 December 2021; Accepted 20 December 2021
Available online 7 January 2022
0263-8223/© 2022 Elsevier Ltd. All rights reserved.
J. Shi et al. Composite Structures 284 (2022) 115161

continuum theory. Later, Khakalo et al. [23] developed a pair of two-

Nomenclature scale plate models that relied on the anisotropic form of Mindlin’s
𝝐̄ global strains strain-gradient thermoelasticity theory. A framework was created by
𝑫 6 × 6 material matrix Goncalves et al. [24] for linear buckling and free-vibration analyses of
∗ sandwich beams using a microstructure-dependent Timoshenko beam
𝛿 2𝐷 , 𝛿 virtual work independent of and related to
model based on a modified coupled-stress theory. The stiffness pa-
the fluctuating function, respectively
rameters of a structural web-core sandwich panel were determined
𝛾𝛼𝛽 , 𝜅𝛼𝛽 in-plane strains and bending curvatures
by unit cell analysis. Karttunen et al. [25] derived a micropolar Tim-
𝜆𝑖 Lagrange multiplier
oshenko beam formulation and used it to model web-core sandwich
⟨⋅⟩ volume integral over the unit cell beams, in which the bending and shear stiffness coefficients were de-
𝑒 , 𝑠 , 𝑡 in-plane strains, transverse strains, and termined through unit cell analysis. A geometrically exact formulation
transverse normal strains, respectively was developed by Chowdhury and Reddy [26] to model the nonlinear
 kinetic energy elastic responses of sandwich beams that were made of architected
𝛺 volume domain occupied by unit cell lattice cores. Experiments and 3D modeling conducted by Lai [27]
𝜔 inertial angular velocity suggested that the lattice core could be modeled as an equivalent ho-
𝜌∗ equivalent density mogeneous orthotropic material in global response analyses. To apply
𝜌𝑐 , 𝜌𝑓 densities of the lattice core and the face the continuum model for analyzing the composite lattice sandwich
sheet, respectively plate, however, it is important for reliable values of the effective plate
𝜏𝑖 , 𝛽𝑖 traction forces on the top and bottom properties to be predicted.
surface of the plate, respectively As mentioned above, the development and design of lattice sand-
wich plates with enhanced mechanical properties has become an in-
𝜃 inclination angle of lattice rod
teresting and hot topic for researchers. In this regard, more efforts are
𝜁 small parameter denoting scale ratio
needed to investigate the static and dynamic behavior to obtain better
𝑎, 𝑏 length and width of unit cell mechanical properties than traditional sandwich plates. In addition, it
𝑓𝑖 body force is necessary to reveal the quantitative relationship between mechanical
ℎ height of the plate properties and various design variables to help the optimal design of
ℎ𝑐 height of lattice core such structure.
ℎ𝑓 thickness of face sheet Recently, Zhong and Yu [28–31] proposed a novel multiscale mod-
𝑙1 , 𝑙2 length and width of the plate, respectively eling technique for composite structures and materials based on the
𝑙𝑐 length of the lattice rod variational asymptotic method (VAM) developed by Hodges and Ces-
𝑟𝑐 lattice rod radius nik [32]. The proposed model of sandwich plate has the advantage of
using small structural parameters (such as the width–thickness ratio)
𝑢𝑖 , 𝑢̄ 𝑖 displacements of 3D-DNS and 2D-RPM
to asymptotically expand the energy functional, and the high-order
𝑉 absolute velocity
terms are removed to obtain the approximate energy of different orders
𝑤 natural frequency and the corresponding dimension reduction model, achieving a good
𝑤𝑖 fluctuation functions tradeoff between accuracy and effectiveness. The disadvantage is that
𝑥𝑖 , 𝑦𝑖 global and local coordinates the width–thickness ratio of the plate should be high enough to meet
𝑓𝑖 , 𝑚𝑖 generalized forces and moments the requirements of the variational asymptotic principle, and only the
large deformation and small rotations are considered.
In the present study, the original geometric nonlinear issue of the
composite sandwich plate with hourglass lattice cores (CSP-HLC) was
of the lattice core are taken into account. To reduce the computa- decomposed into a linear constitutive model over the unit cell and
tional cost and obtain a precision equivalent to that provided by the geometric nonlinear analysis over a 2D reduced-order plate model (2D-
complete three-dimensional computational analysis, various techniques RPM). In the linear constitutive model, the effective plate properties
have been developed to analyze lattice sandwich structures. A critical (such as the matrices 𝐴, 𝐵, and 𝐷) were calculated by homogenizing
review of the literature was presented by Ghugal and Sayyad [9] on the over the unit cell and input into the 2D-RPM to implement global
bending, buckling, and free-vibration analyses of sandwich beams that analysis. If the homogenization result is effective, the macro response
were based on equivalent single layer (ESL) theories, zig-zag theories, of the two-dimensional plate should be the same as that of the three-
layerwise theories, and exact elastic solutions. Vigliotti et al. [10] dimensional model. The accuracy of the homogenization method can be
applied a computational homogenization method to derive a nonlinear verified by comparing the macro response of two models (such as static
constitutive model for lattice cores, and a representative volume ele- displacement and free vibration). Fig. 1 illustrates the analysis process
of the CSP-HLC using a VAM-based reduced-order model. To the best
ment (RVE) of the lattice core was modeled by using discrete structural
of the authors’ knowledge, this technique has never been applied to
elements. Moderate-rotation theory was employed by Cohen et al. [11]
analyze the static displacements and free vibrations of CSP-HLCs.
to derive a fundamental multi-scale model that captured the non-linear
This paper is arranged as follows. The derivation of the 2D reduced-
response of infinitely periodic elastic lattice materials. Abrate and Di
order model of a CSP-HLC based on the VAM is presented in Section 2.
Sciuva [12] provided a comparable single-layer theory characterized by
Section 3 presents the verification of the results and a discussion. In
a single approximation of displacements along the thickness direction.
Section 4, the influences of the layup configurations and structural
TQ Quan and ND Duc [13–15] investigated the nonlinear vibration and parameters on the equivalent stiffnesses and eigenfrequencies of the
dynamic response of functionally graded double curved shallow shells CSP-HLC are investigated. In Section 5, the static and dynamic be-
using Reddy’s third-order shear deformation shell theory, and later ND haviors of a CSP-HLS and a composite sandwich plate with pyramidal
Duc et al. [16–21] provided the nonlinear dynamic analysis of auxetic lattice cores are compared. Finally, the conclusions are given in Sec-
composite and structures by using Runge–Kutta and Galerkin methods. tion 6. Because the procedure is fairly similar that presented in the
Khakalo et al. [22] developed a model of the size-dependent bend- authors’ previous publications, some text and formulas are repeated
ing, buckling, and vibration phenomena of two-dimensional (2D) tri- from these previous publications to make the present article more
angular lattices with the aid of a simplified strain-gradient elastic self-contained.

2
J. Shi et al. Composite Structures 284 (2022) 115161

Fig. 1. Analysis process of composite sandwich plate with hourglass lattice cores (CSP-HLC) using VAM-based equivalent model: (a) three-dimensional (3D) direct numerical
simulation model (3D-DNS) of CSP-HLC; (b) 3D unit cell; (c) two-dimensional (2D) reduced-order plate model (2D-RPM).

2. Theoretical formulation The 3D strain field can be obtained by substituting Eq. (2) into
Eq. (5) and ignoring higher-order terms according to the VAM:
2.1. Reduced-order modeling of CSP-HLC 𝜀11 = 𝛾11 + 𝜁 𝑦3 𝜅11 + 𝑤1,1 ,
2𝜀12 = 2𝛾12 + 2𝜁 𝑦3 𝜅12 + 𝑤1,2 + 𝑤2,1 ,
Two sets of coordinates are presented for derivation convenience: 𝜀22 = 𝛾22 + 𝜁 𝑦3 𝜅22 + 𝑤2,2 ,
(6)
global coordinates 𝑥𝑖 and local coordinates 𝑦𝑖 . For the reduced-order 2𝜀13 = 𝑤1,3 + 𝑤3,1 ,
model, the field variable is a function that is defined on the reference 2𝜀23 = 𝑤2,3 + 𝑤3,2 ,
plane, and 𝑥3 vanishes (as illustrated in Fig. 1(c)). The relationship 𝜀33 = 𝑤3,3 ,
between the local and global coordinates can be written as 𝑦𝑖 = 𝑥𝑖 ∕𝜁 where
(𝜁 is a small parameter, 𝑖, 𝑗 = 1, 2, 3, and 𝛼 = 1, 2). In the reduced-order ( ) 1( ) ( )
model, the field variables of the original CSP-HLC can be represented 𝛾𝛼𝛽 𝑥1 , 𝑥2 = 𝑢̄ + 𝑢̄ 𝛽,𝛼 , 𝜅𝛼𝛽 𝑥1 , 𝑥2 = −𝑢̄ 3,𝛼𝛽 . (7)
2 𝛼,𝛽
by the remaining global coordinates 𝑥𝛼 and local coordinates 𝑦𝑖 , and
The three-dimensional strain field  can be expressed in matrix
its partial derivative is
( ) ( ) ( ) forms as
𝜕𝑓 𝑥𝛼 , 𝑦𝑖 , 𝑡 𝜕𝑓 𝑥𝛼 , 𝑦𝑖 , 𝑡 || 1 𝜕𝑓 𝑥𝛼 , 𝑦𝑖 , 𝑡 || [ ]T
= | + | 𝑒 = 𝜀11 𝜀22 2𝜀12 = 𝜸 + 𝑥3 𝜿 + 𝑰 𝛼 𝒘∥,𝛼 ,
𝜕𝑥𝛼 𝜕𝑥𝛼 | | [ ]T
|𝑦𝑖 =𝑐𝑜𝑛𝑠𝑡 𝜁 𝜕𝑦𝑖 |𝑥𝛼 =𝑐𝑜𝑛𝑠𝑡 (1) 2𝑠 = 2𝜀13 2𝜀23 = 𝒘∥,3 + 𝒆𝛼 𝑤3,𝛼 , (8)
1 𝑡 = 𝜀33 = 𝑤3,3 ,
≡ 𝑓,𝛼 + 𝑓;𝑖 .
𝜁
where 𝑒 , 𝑠 , and 𝑡 are the in-plane strains, transverse shear strains,
To deduce the reduced-order model of the CSP-HLC using the
and transverse normal strains of the 3D-DNS, respectively, ()∥ =
VAM, it is necessary to represent the 3D displacement field by using [ ]T [ ]T [ ]T
()1 ()2 , 𝜸 = 𝛾11 2𝛾12 𝛾22 , 𝜿 = 𝜅11 𝜅12 + 𝜅21 𝜅22 , and
two-dimensional plate variables, as follows:
( ) ( ) ( ) ( )
𝑢1 𝑥𝛼 , 𝑦𝑖 , 𝑡 = 𝑢̄ 1 𝑥𝛼 , 𝑡 − 𝜁 𝑦3 𝑢̄ 3,1 𝑥𝛼 , 𝑡 + 𝜁 𝑤1 𝑥𝛼 , 𝑦𝑖 , 𝑡 , { } { }
( ) ( ) ( ) ( ) ⎡ 1 0 ⎤ ⎡ 0 0 ⎤
1 0
𝑢2 𝑥𝛼 , 𝑦𝑖 , 𝑡 = 𝑢̄ 2 𝑥𝛼 , 𝑡 − 𝜁 𝑦3 𝑢̄ 3,2 𝑥𝛼 , 𝑡 + 𝜁 𝑤2 𝑥𝛼 , 𝑦𝑖 , 𝑡 , (2) 𝒆1 = , 𝒆2 = , 𝑰1 = ⎢ 0 1 ⎥,𝑰 = ⎢ 1
2 0 ⎥. (9)
( ) ( ) ( ) 0 1 ⎢ ⎥ ⎢ ⎥
𝑢3 𝑥𝛼 , 𝑦𝑖 , 𝑡 = 𝑢̄ 3 𝑥𝛼 , 𝑡 + 𝜁 𝑤3 𝑥𝛼 , 𝑦𝑖 , 𝑡 , ⎣ 0 0 ⎦ ⎣ 0 1 ⎦

where 𝑢̄ 𝑖 and 𝑢𝑖 are the displacements of the 2D-RPM and the original The strain energy of the CSP-HLC can be written as
3D-DNS, respectively, and 𝑤𝑖 are fluctuating functions to be solved. The
1 𝑙 ∕22 𝑙 ∕2
underlined terms are the deformations of the reference surface of the 1 1
 =  𝑑𝑥 𝑑𝑥 , (10)
CSP-HLC, which are constrained by the following conditions: 2 ∫−𝑙1 ∕2 ∫−𝑙2 ∕2 𝛺 𝛺 2 1
( ) ( ) where
ℎ𝑢̄ 𝑎 𝑥𝑎 = ⟨𝑢𝑎 ⟩ + ⟨𝜂𝑦3 ⟩ 𝑢̄ 3,2 , ℎ𝑢̄ 3 𝑥𝛼 = ⟨𝑢3 ⟩ , (3)
ℎ𝑐 𝑏 𝑎
2
+𝑡𝑓 2 2
where ⟨⋅⟩ denotes the integral over the volume of the unit cell. 𝛺 = 𝑡𝑇 𝐷𝑡 𝑡 𝑑𝑦1 𝑑𝑦2 𝑑𝑦3
∫ ℎ𝑐 ∫− 𝑏 ∫− 𝑎
From Eq. (3), we can conclude that the fluctuating functions were 2 2 2
ℎ 𝑏 𝑎
constrained by the following conditions: − 2𝑐 2 2
+  𝑇 𝐷  𝑑𝑦 𝑑𝑦 𝑑𝑦 (11)
∫ ℎ
− 2𝑐 −𝑡𝑓 ∫− 𝑏 ∫− 𝑎 𝑏 𝑏 𝑏 1 2 3
⟨𝜁 𝑤𝑖 ⟩ = 0. (4) 2 2
√
𝑎2 +𝑏2
√ ′ 𝑟
2 cos 𝜃
𝑟2 −𝑦1 2
If we address the issue of large deformations and small local rota- 𝑟𝑇 𝐷𝑟 𝑟 𝑑𝑦′1 𝑑𝑦′2 𝑑𝑦′3 .
∫−√𝑟2 −𝑦′ ∫− 𝑟
+4 × √
∫− 𝑎2 +𝑏2
tions, the 3D strain can be approximated in a linear form based on the 2 cos 𝜃 1 2

decomposition of the rotation tensor, as follows: with the subscripts 𝑡, 𝑏 and 𝑟 denoting the top face sheet, the bot-
( ) tom face sheet and the lattice rod, respectively; 𝑦′𝑖 denoting the local
1 𝜕𝑢𝑖 𝜕𝑢𝑗
𝜀𝑖𝑗 = + . (5) coordinate system along each lattice rod, as illustrated in Fig. 2.
2 𝜕𝑥𝑗 𝜕𝑥𝑖

3
J. Shi et al. Composite Structures 284 (2022) 115161

The common method to solve Eq. (19) is to assume the unknown

fluctuating function 𝑤𝑖 to directly reduce the original 3D continuum
model into a 2D plate model. However, for a CSP-HLC made of gen-
eral anisotropic and heterogeneous materials, the imposition of this
assumption may introduce significant errors. Fortunately, the solution
of 𝑤𝑖 can be obtained through the asymptotic analysis of the variational
statement in Eq. (19), which will be discussed in the next section.

2.2. Reduced-order analysis of CSP-HLC

Following the VAM, the zeroth-order approximation of the varia-

tional statement in Eq. (19) can be obtained as follows:
𝑡2 [ ( ) ]
𝛿 2𝐷 − 0 𝑑𝛺 + 𝛿 2𝐷 𝑑𝑡 = 0, (20)
∫ 𝑡1 ∫𝛺

where 0 can be obtained from Eq. (12) by omitting the derivatives

with respect to 𝑥𝛼 in Eq. (8), as follows:
Fig. 2. Diagram of lattice unit cell for strain energy integral. ⟨ ( )T ( ) ⟩
𝜸 + 𝑥3 𝜿 𝑫 𝑒 𝜸 + 𝑥3 𝜿 + 𝑤T∥,3 𝑫 𝑠 𝒘∥,3 + 𝑤T3,3 𝑫 𝑡 𝑤3,3
20 = ( )T ( ) . (21)
+ 2 𝜸 + 𝑥3 𝜿 𝑫 𝑒𝑠 𝒘∥,3 + 𝑫 𝑒𝑡 𝑤3,3 + 2𝒘T∥,3 𝑫 𝑠𝑡 𝑤3,3
Eq. (11) can be written compactly as
Following the principle of minimum total potential energy, the
T fluctuating function can be solved by minimizing the zeroth-order
⟨⎧  ⎫ ⎡ 𝑫𝑒 𝑫 𝑒𝑠 𝑫 𝑒𝑡 ⎤ ⎧  ⎫⟩
1⟨ T ⟩ 1 ⎪ 𝑒 ⎪ ⎢ T ⎥⎪ 𝑒 ⎪ approximate strain energy under the constraint in Eq. (4), as follows:
 =  𝑫 =
2 ⎨ ⎬ ⎢ 𝑫 𝑒𝑠
2𝑠 𝑫𝑠 𝑫 𝑠𝑡 ⎥ ⎨ 2𝑠 ⎬ ,
2 ⎪ 𝑡 ⎪ ⎢ 𝑫T
⎩ ⎭ ⎣ 𝑒𝑡 𝑫 T𝑠𝑡 𝑫 𝑡 ⎥⎦ ⎪
⎩ 𝑡
⎪
⎭ 𝛿0 = 0. (22)
(12) The corresponding Euler–Lagrange equation can be obtained by
where 𝑫 𝑒 , 𝑫 𝑒𝑠 , 𝑫 𝑒𝑡 , 𝑫 𝑠 , 𝑫 𝑠𝑡 , and 𝑫 𝑡 are the corresponding sub-matrices partially integrating Eq. (21), as follows:
of the 3D 6 × 6 material matrix. [( )T ( )T ]
𝜸 + 𝑥3 𝜿 𝑫 𝑒𝑠 + 𝒘∥,3 𝑫 𝑠 + 𝑤3,3 𝑫 𝑠𝑡 = 𝝀∥ ,
The virtual work that is done by the external load can be written as [( ] ,3
)T ( )T (23)
𝜸 + 𝑥3 𝜿 𝑫 𝑒𝑡 + 𝒘∥,3 𝑫 𝑠𝑡 + 𝑤3,3 𝑫 𝑡 = 𝜆3 ,
∗ ,3
𝛿 3𝐷 = 𝛿 2𝐷 + 𝛿 , (13) [ ]T
∗ where 𝝀∥ = 𝜆1 𝜆2 and 𝜆3 are introduced Lagrange multipliers.
where 𝛿 2𝐷 and 𝛿 represent the virtual work independent of and For a free surface condition, the expressions in square brackets in
related to the fluctuating function, respectively. These are defined as
Eq. (23) should be zero at the bottom and top of the plate, and the
follows:
⟨ ⟩ ∗ ⟨ ⟩ boundary conditions of the bottom and top of the plate can be defined
𝛿 2𝐷 = 𝑝𝑖 𝛿𝑣𝑖 + 𝑞𝑎 𝛿𝑣3,𝑎 , 𝛿 = ⟨𝑓𝑖 𝛿𝑤𝑖 ⟩ + 𝜏𝑖 𝛿𝑤+ −
𝑖 + 𝛽𝑖 𝛿𝑤𝑖 , (14) as
[( )T ]+∕−
where (⋅)+ and (⋅)− denote the items that act on the top and bottom 𝜸 + 𝑥3 𝜿 𝑫 𝑒𝑠 + 𝑤T∥,3 𝑫 𝑠 + 𝑤3,3 𝑫 T𝑠𝑡 = 0,
surfaces of the plate, respectively; 𝜏𝑖 , 𝛽𝑖 , and 𝑓𝑖 are the traction forces [( ]+∕− (24)
)T T
on the top and bottom surfaces of the plate ( and) the body force, 𝜸 + 𝑥3 𝜿 𝑫 𝑒𝑡 + 𝑤∥,3 𝑫 𝑠𝑡 + 𝑤3,3 𝑫 𝑡 = 0,
respectively; 𝑝𝑖 = ⟨𝑓𝑖 ⟩ + 𝑎𝑖 + 𝛽𝑖 ; and 𝑞𝑎 = ℎ∕2 𝛽𝑎 − 𝑎𝑎 − ⟨𝑥3 𝑓𝑎 ⟩.
To calculate the kinetic energy of the CSP-HLC, the absolute velocity where the superscript ‘‘+/−’’ denotes the item at the bottom/top surface
of a point in the CSP-HLC can be evaluated as of the plate.
The solutions of 𝒘∥ and 𝑤3 can be obtained by substituting Eq. (24)
𝑣 = 𝑉 + 𝜔(𝜉 ̇
̃ + 𝑤) + 𝑤, (15)
into Eq. (23), as follows:
where 𝑤̇ = 𝜕𝑤∕𝜕𝑡; 𝜔 is the inertial angular velocity, 𝑉 is the absolute ⟨ ( ) ( )−1 ⟩T ⟨ ( ) ( )−1 ⟩
velocity of a point in the deformed reference surface, ()̃ 𝑖𝑗 = −𝑒𝑖𝑗𝑘 ()𝑘 𝒘∥ = − 𝜸 + 𝑥3 𝜿 𝑫 𝑒𝑠 𝑫 𝑠 , 𝑤3 = − 𝜸 + 𝑥3 𝜿 𝑫 𝑒𝑡 𝑫 𝑡 , (25)
with 𝑒𝑖𝑗𝑘 as the permutation symbol, and 𝜉 = [ 0 0 𝑥3 ]𝑇 .
The kinetic energy of the plate can be written as where
( )T ( )−1 ( )−1
1 𝑫 𝑒𝑠 = 𝑫 𝑒𝑠 − 𝑫 𝑒𝑡 𝑫 𝑠𝑡 𝑫𝑡 , 𝑫 𝑒𝑡 = 𝑫 𝑒𝑡 − 𝑫 𝑒𝑠 𝑫 𝑠 𝑫 𝑠𝑡 ,
= 𝜌𝑣𝑇 𝑣𝑑,  = 2𝐷 + ∗ , (16) (26)
2 ∫ ( )T ( )−1
𝑫 𝑡 = 𝑫 𝑡 − 𝑫 𝑠𝑡 𝑫𝑠 𝑫 𝑠𝑡 .
where 𝜌 is the mass of the plate and
( ) Eq. (25) is inserted into Eq. (21), yielding the following:
1
2𝐷 = ̄ 𝑇 𝑉 + 2𝜔𝑇 𝜇𝜉𝑉
𝜇𝑉 ̃ + 𝜔𝑇 𝛷𝜔 𝑑𝛺, (17) ⟨ )⟩
2 ∫𝛺 1 ( )T (
𝑈2𝐷 = 𝜸 + 𝑥3 𝜿 𝑫 𝑒 𝜸 + 𝑥3 𝜿
1 [ ] 2
∗ = 𝜌 (𝜔𝑤 ̇ 𝑇 (𝜔𝑤
̃ + 𝑤) ̃ + 𝑤) ̃ 𝑇 (𝜔𝑤
̇ + 2(𝑉 + 𝜔𝜉) ̇ 𝑑,
̃ + 𝑤) (18) { }T [ ]{ }, (27)
2 ∫ 1 𝜸 𝑨 𝑩 𝜸
⟨ 2 ⟩ =
2 𝜿 𝑩𝑇 𝑫 𝜿
⌊ ]𝑇 ⎡ 𝑥3 𝜌 0
⟨ 2 ⟩
0 ⎤
𝜇̄ = ⟨𝜌⟩, 𝜇𝜉 = 0 0 ⟨𝑥3 𝜌⟩ , and 𝛷 = ⎢ 0 𝑥3 𝜌 0 ⎥. where 𝑨, 𝑩, and 𝑫 are the 3 × 3 tensile stiffness matrix, bending–
⎢ ⎥
⎣ 0 0 0 ⎦ torsional coupled stiffness matrix, and bending stiffness matrix, respec-
The elastodynamic behavior of the CSP-HLC is determined by tively, which can be calculated as follows:
Hamilton’s principle: ⟨⟨ ⟩⟩ ⟨⟨ ⟩⟩ ⟨⟨ ⟩⟩
𝑡2 [ ( ) ∗
] 𝑨= 𝑫𝑒 , 𝑩 = 𝑥3 𝑫 𝑒 , 𝑫 = 𝑥23 𝑫 𝑒 ,
𝛿 2𝐷 + ∗ −  + 𝛿 2𝐷 + 𝛿 𝑑𝑡 = 0. (19) (28)
∫𝑡1 𝑫 𝑒 = 𝑫 𝑒 − 𝑫 𝑒𝑠 𝑫 −1 T T
𝑠 𝑫 𝑒𝑠 − 𝑫 𝑒𝑡 𝑫 𝑒𝑡 ∕𝑫 𝑡 .

4
J. Shi et al. Composite Structures 284 (2022) 115161

The constitutive relation of the reduced-order model for the CSP- ∑

∞ ∑
∞
HLC can be expressed as 𝑢̄ 3,2 = 𝜙2 𝑚𝑛 (𝑡) sin(𝛼𝑥1 ) cos(𝛽𝑥2 ), (41)
𝑚=1 𝑛=1
⎧ 𝑁11 ⎫ ⎡ 𝐴11 𝐴12 𝐴16 𝐵11 𝐵12 𝐵16 ⎤⎧ 𝛾11 ⎫
⎪ ⎪ ⎢ ⎥⎪ ⎪ where 𝛼 = 𝑚𝜋∕𝑙1 , 𝛽 = 𝑛𝜋∕𝑙2 , 𝑈𝑚𝑛 (𝑡), 𝑉𝑚𝑛 (𝑡), 𝑊𝑚𝑛 (𝑡), 𝜙1 𝑚𝑛 (𝑡), and 𝜙2 𝑚𝑛 (𝑡)
⎪ 𝑁22 ⎪ ⎢ 𝐴12 𝐴22 𝐴26 𝐵12 𝐵22 𝐵26 𝛾22
⎪ ⎪ ⎢ ⎥⎪ ⎪ are functions of time 𝑡.
⎨
𝑁12
⎬=⎢
𝐴16 𝐴26 𝐴66 𝐵16 𝐵26 𝐵66 ⎥⎪ 2𝛾12 ⎪
(29)
⎥⎨ ⎬. For free vibrations, the periodic solution is assumed to be
⎪ 𝑀11 ⎪ ⎢ 𝐵11 𝐵12 𝐵16 𝐷11 𝐷12 𝐷16 𝜅11
⎪ ⎪ ⎢ ⎥⎪ ⎪
⎧ 𝑈𝑚𝑛 (𝑡) ⎫ ⎧ 𝑈𝑚𝑛 0
⎫ 0
⎧ 𝑈𝑚𝑛 ⎫
⎪
𝑀22
⎪ ⎣
𝐵12 𝐵22 𝐵26 𝐷12 𝐷22 𝐷26 ⎥⎪ 𝜅22 ⎪
⎩ 𝑀12 ⎭ 𝐵16 𝐵26 𝐵66 𝐷16 𝐷26 𝐷66 ⎦⎪
⎩ 2𝜅12 ⎪
⎭
⎪ 𝑉𝑚𝑛 (𝑡) ⎪ ⎪ 0
𝑉𝑚𝑛 ⎪ ⎪ 𝑉𝑚𝑛 ⎪
0
⎪ ⎪ ⎪ 0 ⎪ 𝑖𝑤𝑡 ⎪ 0 ⎪
⎨ 𝑊𝑚𝑛 (𝑡) ⎬ = ⎨ 𝑊𝑚𝑛 ⎬ e ⎨ 𝑊𝑚𝑛 ⎬ ≠ 0, (42)
The three-dimensional strain field in the zeroth-order approxima- ⎪ 𝜙1 𝑚𝑛 (𝑡) ⎪ ⎪ 𝜙1 𝑚𝑛 ⎪
0 ⎪ 𝜙0 ⎪
tion can be recovered as ⎪ ⎪ ⎪ ⎪ ⎪ 10 𝑚𝑛 ⎪
⎩ 𝜙2 𝑚𝑛 (𝑡) ⎭ ⎩ 𝜙02 𝑚𝑛 ⎭ ⎩ 𝜙2 𝑚𝑛 ⎭
𝑒0 = 𝜸 + 𝑥3 𝜿, 2𝑠0 = −𝒘∥,3 , 𝑡0 = 𝑤3,3 . (30) √
where 𝑖 = −1, and 𝑤 is the natural frequency of the structure.
The local stress field can be recovered according to the generalized The elastic curved surface differential equation of the 2D-RPM
Hooke’s law: under free flexural vibrations is
𝜕 4 (𝑢)
̄ 𝜕 4 (𝑢)
̄ 𝜕 4 (𝑢)
̄
𝝈 = 𝑫 𝑒 . (31) 𝐷11 + 𝐷22 + 𝐷66 = 𝑞𝑖 , (43)
𝜕𝑥41 𝜕𝑥42 𝜕𝑥21 𝜕𝑥22
The macroscopic plate model is governed by the variational state- 2
ment in Eq. (20), as it only involves the 2D field variables in terms where 𝑞𝑖 = −𝜌∗ ℎ 𝜕𝜕𝑡2𝑢̄ is the inertial force.
of the remaining macro-coordinates, 𝑥𝛼 . Therefore, the 2D-RPM can Eqs. (34) and (37)–(42) are inserted into Eq. (43), yielding the
be used to represent the original CSP-HLC in the free-vibration and following:
static-displacement analysis and can be implemented in a finite element √
software package, such as ABAQUS/Standard. 𝐷11 𝛼 4 + 𝐷22 𝛽 4 + 2𝐷66 𝛼 2 𝛽 2
𝜔= . (44)
𝜌∗ ℎ
2.3. Equivalent density of CSP-HLC The natural frequencies of the plates with other boundary con-
ditions can be obtained by solving the characteristic equation after
It is assumed that the CSP-HLC consists of top and bottom face substituting the corresponding displacement into the partial differential
sheets, a cemented layer, and a core. Because the cemented layer is equation in Eq. (42).
very thin and light, the mass and volume of the cemented layer can be
ignored. The total mass of the unit cell within the CSP-HLC is 3. Validation example

𝑚 = 𝑚𝑐 + 𝑚𝑓 = 𝜌𝑐 ⋅ 4𝜋𝑟2𝑐 𝑙𝑐 + 2𝜌𝑓 ⋅ 𝑡𝑓 ⋅ 2𝑎𝑏, (32) In this section, the proposed model was applied to the static and dy-
where 𝜌𝑐 and 𝜌𝑓 are the density of the lattice core and the face sheet, namic analysis of the CSP-HLC, and the results were compared with the
respectively, and 𝑡𝑓 is the thickness of the face sheet. 3D-DNS results, as illustrated in Fig. 3. The relative error between the
The total volume occupied by the unit cell is 2D-RPM and 3D-DNS was defined as error = |2D-RPM 3D-DNS
results−3D-DNS results|
results
× 100%. As shown in Fig. 4(a), the 3D-DNS had 14 unit cells in the
𝑉 = 𝑎𝑏(ℎ𝑐 + 2𝑡𝑓 ). (33) 𝑥1 -direction and 10 unit cells in the 𝑥2 -direction, and the overall
dimensions of the plate were 560 mm × 400 mm.
The equivalent density of the CSP-HLC can be obtained as
The dimensions of the unit cell are shown in Fig. 4(b). The thickness
𝑚 𝜌𝑐 ⋅ 4𝜋𝑟2𝑐 𝑙𝑐 + 2𝜌𝑓 ⋅ 𝑡𝑓 ⋅ 2𝑎𝑏 of the face sheet was 𝑡𝑓 = 1 mm, and the layup configuration was
𝜌∗ = = . (34)
𝑉 𝑎𝑏(ℎ𝑐 + 2𝑡𝑓 ) [±45◦ ∕0◦ ∕90◦ ]s . The structural parameters of the unit cells were 𝑎 =
𝑏 = 40 mm, 𝑟𝑐 = 1.5 mm, ℎ𝑐 = 10 mm, and 𝜃 = 10◦ . The material
2.4. Calculation of natural frequencies using 2D-RPM of the face sheet was T300/epoxy composite, with 𝐸1 = 208.8 GPa,
𝐸2 = 𝐸3 = 43 GPa, 𝐺12 = 𝐺13 = 7.42 GPa, 𝐺23 = 3.98 GPa, 𝑣12 = 𝑣13 =
The boundary conditions of a rectangular plate with four edges 0.3, and 𝑣23 = 0.42 [33]. The material of the core was aluminum, with
simply supported are assumed to be 𝐸𝑐 = 70 GPa and 𝑣𝑐 = 0.3. The equivalent plate properties obtained by
homogenizing the unit cell are shown in Table 1 for reference.
̄ 11 = 0
𝑢̄ 2 = 𝑢̄ 3 = 𝑢̄ 3,2 = 𝑁11 = 𝑀 at 𝑥1 = 0, 𝑙1 , (35) The typical boundary conditions for static-displacement and free-
vibration analysis are shown in Fig. 5. The naming convention begins
̄ 22 = 0
𝑢̄ 1 = 𝑢̄ 3 = 𝑢̄ 3,1 = 𝑁22 = 𝑀 at 𝑥2 = 0, 𝑙2 . (36) with a pair of boundaries on the vertical side of the CSP-HLC, followed
by a pair of boundaries on the horizontal side. For instance, CCSS
According to the Navier method, to satisfy the boundary conditions, indicates that the plate was clamped on the left and right and simply
it is assumed that the displacement components 𝑢̄ 1 , 𝑢̄ 2 , 𝑢̄ 3 , 𝑢̄ 3,1 , and 𝑢̄ 3,2 supported at the bottom and top.
can be respectively expressed by the following double trigonometric
series: 3.1. Static-displacement analysis
∑
∞ ∑
∞
𝑢̄ 1 = 𝑈𝑚𝑛 (𝑡) cos(𝛼𝑥1 ) sin(𝛽𝑥2 ), (37)
To verify the accuracy of the model, a 0.5 MPa uniform load
𝑚=1 𝑛=1
was applied on the top surface of the plate, and the displacement
∑∞ ∑ ∞
𝑢̄ 2 = 𝑉𝑚𝑛 (𝑡) sin(𝛼𝑥1 ) cos(𝛽𝑥2 ), (38) distributions along path 1 predicted by the 2D-RPM and 3D-DNS are
𝑚=1 𝑛=1 compared in Fig. 6. The displacement distribution predicted by the 2D-
∑∞ ∑ ∞ RPM was basically consistent with that predicted by the 3D-DNS. The
𝑢̄ 3 = 𝑊𝑚𝑛 (𝑡) sin(𝛼𝑥1 ) cos(𝛽𝑥2 ), (39) difference was due to the different meshes used in the two models.
𝑚=1 𝑛=1 The vertical displacement error of 4.3% under the CCCC boundary
∑∞ ∑ ∞
condition was the largest, and the vertical displacement error of 1.3%
𝑢̄ 3,1 = 𝜙1 𝑚𝑛 (𝑡) cos(𝛼𝑥1 ) sin(𝛽𝑥2 ), (40)
under the SSSS boundary condition was the smallest.
𝑚=1 𝑛=1

5
J. Shi et al. Composite Structures 284 (2022) 115161

Fig. 3. Comparative analysis process of 2D-RPM and 3D-DNS.

Fig. 4. Geometry and dimensions of CSP-HLC and its unit cell.

Table 1
Effective plate properties of CSP-HLC after homogenization (unit: SI).
2.676 E+05 8.026 E+04 −4.208 E+00 1.096 E+06 4.013 E+05 −2.645 E+01
2.676 E+05 −1.263 E+00 4.013 E+05 1.338 E+06 −7.937 E+00
7.999 E+04 −1.824 E+01 −5.443 E+00 3.999 E+05
1.344 E+07 4.357 E+06 −1.586 E+02
sym. 1.344 E+07 −4.735 E+01
4.083 E+06

Fig. 5. Boundary conditions used in the static-displacement and free-vibration analysis.

6
J. Shi et al. Composite Structures 284 (2022) 115161

Fig. 6. Comparison of vertical displacements along path 1 under uniform load of 0.5 MPa and different boundary conditions.

The vertical displacement error decreased gradually from Fig. 6(a) 3.2. Local field recovery
to (f), which may have been because the clamp support had the greatest
The local field distribution within the unit cell closet to the center
influence on the constraint of the model in the applied boundary
point (as shown by the blue dot in Fig. 7) of the 2D-RPM was recovered
conditions. The displacement comparison showed that the vertical by Eqs. (30)–(31) under the SSSS boundary condition.
displacement error under each boundary condition was within 5%, Fig. 8 shows that the maximum values of 𝑈 and 𝑈3 were located
indicating that the 2D-RPM could predict the static displacements of at the intersection of the inclined lattice rod and the face sheet closest
the CSP-HLC with high precision and efficiency. It also verified that to the geometric center of the plate. The local displacement gradually
decreased along the 𝑦1 - and 𝑦2 -directions from the maximum value,
the equivalent plate properties calculated by the VAM had sufficient
which agreed with the variation trend of the global displacement. The
accuracy, and thus, it could be used to evaluate the macro-performance maximum vertical displacement 𝑈3 within the unit cell was 9.625 mm,
of the CSP-HLC. and the relative error was 0.997% compared with the 3D-DNS result

7
J. Shi et al. Composite Structures 284 (2022) 115161

Table 2
Comparison of first four natural frequencies (Hz) and modal shapes of CSP-HLC under CCCC boundary condition predicted by 2D-RPM and
3D-DNS. Mesh was removed for clear observation (values in red brackets denote relative errors).
1 2
Mode
3D-DNS 2D-RPM 3D-DNS 2D-RPM

3 4
Mode
3D-DNS 2D-RPM 3D-DNS 2D-RPM

consistent between the 3D-DNS and 2D-RPM. All of the relative errors,
however, were less than 5%, indicating that the 2D-RPM had high
precision in predicting the natural frequencies of the CSP-HLC, and the
equivalent stiffnesses that were obtained by the VAM were reasonable.

To further verify the efficiency and accuracy of the 2D-RPM in

vibration mode analysis, the free-vibration characteristics (fundamental
frequencies and related modes) of the CSP-HLC under different bound-
ary conditions are shown in Table 3. Different boundary conditions had
a significant impact on the mode shape. For instance, the mode shapes
of the CCCC and CCSS boundary conditions were one half-wave in the
horizontal and vertical directions, while the mode shapes of the CCFF
Fig. 7. Schematic diagram of the recovered point in the 2D-RPM. and CSFF were one half-wave in the horizontal direction. However, the
mode shapes of the 3D-DNS and 2D-RPM were basically the same under
the same boundary conditions.
(the maximum value of 𝑈3 was 9.53 mm), indicating that the recovered It is worth noting that there were certain differences in the rel-
displacement distribution had high accuracy and could be used to ative errors of the fundamental frequency under different boundary
evaluate the location of the maximum local displacement. conditions. For example, the frequency error under the CCCC boundary
The recovered local stress distributions within the unit cell are condition was 2.76%, while the error under the SSSS boundary condi-
shown in Fig. 9. There was a large stress concentration at the connec- tion was only 0.98%, which may have been because the degree to which
tion between the inclined lattice rods and the face sheets, and the von the constraints were not completely consistent between the 3D-DNS
Mises stress distribution in the other inclined lattice rods was relatively and 2D-RPM under the CCCC boundary condition was greater than that
uniform. It is worth noting that the maximum values of 𝜎12 and 𝜎13 were under the SSSS boundary condition. Although there were differences,
located in the inclined lattice rods, indicating that the inclined lattice the maximum error of the fundamental frequency was less than 3%.
rods played an important role in the load transfer process. With the enhancement of the boundary constraints from SSSS to CCCC,
the fundamental frequency increased from 260.64 to 916.52 Hz when
3.3. Free-vibration analysis predicted by the 2D-RPM and from 255.46 Hz to 891.90 Hz when
predicted by the 3D-DNS. These results are consistent with the common
To further verify the accuracy of the proposed model, the first four sense notion that the stronger the boundary constraint is, the greater
natural frequencies and modal shapes of CSP-HLC under CCCC bound- the fundamental frequency becomes.
ary conditions predicted by the 2D-RPM and 3D-DNS are compared in Table 4 lists the first four natural frequencies of the CSP-HLC
Table 2. The reason for choosing the CCCC boundary condition is that predicted by the 2D-RPM and 3D-DNS under different boundary con-
the relative error was the largest in the static-displacement analysis, as ditions. The data inside and outside of the brackets were the results
shown in Fig. 6(a). of the 2D-RPM and 3D-DNS, respectively. The relative errors of the
Table 2 shows that the vibrational modes predicted by 2D-RPM natural frequencies were less than 5%, indicating that the 2D-RPM
were highly consistent with those predicted by the 3D-DNS, and the can be used for high-order natural frequency analysis under different
half-wave numbers of each mode along the horizontal direction were boundary conditions.
the same. The relative errors of first-order to fourth-order natural Regarding the computational efficiency, the proposed model was
frequencies were 2.76%, 2.99%, 3.96%, and 4.55% respectively, which more cost-efficient than the 3D-DNS during modeling. The proposed
may have been due to the boundary constraints being not completely method required 1160 shell elements (S4R) in total compared to the

8
J. Shi et al. Composite Structures 284 (2022) 115161

Fig. 8. Local displacement distribution within the unit cell closet to the center point of the 2D-RPM, and the top face sheet was removed for clear observation (unit: mm).

Fig. 9. Local stress distribution within the unit cell closet to the center point of the 2D-RPM, and the top face sheet was removed for clear observation (unit: MPa).

9
J. Shi et al. Composite Structures 284 (2022) 115161

Table 3
Comparison of fundamental frequencies (Hz) and modal shapes of CSP-HLC under different boundary conditions predicted by 2D-RPM and
3D-DNS. Mesh was removed for clear observation (values in red brackets denote relative errors).
CCCC CCFF
B.C.a
3D-DNS 2D-RPM 3D-DNS 2D-RPM

CCSS CSFF
B.C.
3D-DNS 2D-RPM 3D-DNS 2D-RPM

CSSS SSSS
B.C.
3D-DNS 2D-RPM 3D-DNS 2D-RPM

a
Boundary conditions.

Fig. 10. Unit cell of CSP-HLC with different 𝑡𝑓 ∕ℎ𝑐 (𝑟𝑐 = 1.5 mm, 𝑎 = 𝑏 = 40 mm, 𝜃 = 10◦ ).

Table 4 Table 5
Comparison of first four natural frequencies (Hz) of CSP-HLC under different boundary Geometric parameters of unit cell used for parametric study.
conditions predicted by 2D-RPM and 3D-DNS (the values in brackets were predicted Geometric parameters
by 2D-RPM). Layup Annotations
𝑡1 ∕ℎ𝑐 𝜃 𝑟𝑐 𝑎∕𝑏
Mode CCCC CCFF CCSS CSFF CSSS SSSS
0.1–0.18 10◦ 1.5 mm 1.0 [±45◦ ∕90◦ ∕0◦ ]𝑠 Fig. 10
891.90 355.60 660.83 255.46 618.50 592.46
1 0.1 10◦ –14◦ 1.5 mm 1.0 [±45◦ ∕90◦ ∕0◦ ]𝑠 Fig. 11
(916.52) (362.84) (668.08) (260.64) (627.51) (598.29)
0.1 10◦ 1.0–2.0 mm 1.0 [±45◦ ∕90◦ ∕0◦ ]𝑠 Fig. 12
1416.99 461.75 1269.57 393.04 1138.03 1026.45 0.1 10◦ 1.5 mm 1.0–1.4 [±45◦ ∕90◦ ∕0◦ ]𝑠 Fig. 13
2
(1459.40) (479.73) (1297.00) (405.74) (1167.00) (1052.70)
2096.00 967.05 1488.00 778.53 1467.70 1453.70
3
(2179.10) (995.15) (1547.20) (816.23) (1520.50) (1497.60)
2255.00 1014.10 2094.80 962.07 1948.20 1783.00
of RAM. In short, compared to the 3D-DNS, the 2D-RPM reduced the
4 computational effort greatly but still achieved very high precision.
(2357.70) (1052.40) (2178.60) (1013.60) (2037.30) (1852.20)

4. Parametric study

182,118 solid elements (C3D10) of the 3D-DNS. The 2D-RPM was

4.1. Structural parameters
also more time-efficient than the 3D-DNS in performing free-vibration
analysis, requiring 28 min with one CPU as opposed to 2 h with four To further study the static and dynamic behaviors of the CSP-HLC,
CPUs. The workstation was a Dell Precision Tower 7910 powered by the unit cell with different structural parameters was homogenized
Intel Xeon CPU E5-2660 with a clock rate of 2.60 GHz and 256 GB and input into the 2D-RPM for static-displacement and free-vibration

10
J. Shi et al. Composite Structures 284 (2022) 115161

Fig. 11. Unit cell of CSP-HLC with different inclination angles of lattice rod (𝑟𝑐 = 1.5 mm, 𝑎 = 𝑏 = 40 mm, ℎ𝑐 = 10 mm, 𝑡𝑓 = 1.0 mm).

Fig. 12. Unit cell of CSP-HLC with different lattice rod radii (𝜃 = 10◦ , 𝑎 = 𝑏 = 40 mm, ℎ𝑐 = 10 mm, 𝑡𝑓 = 1.0 mm).

Fig. 13. Unit cell of CSP-HLC with different aspect ratios (𝜃 = 10◦ , 𝑟𝑐 = 1.5 mm, ℎ𝑐 = 10 mm, 𝑡𝑓 = 1.0 mm).

Fig. 14. Effect of structural parameters on equivalent stiffnesses of CSP-HLC.

11
J. Shi et al. Composite Structures 284 (2022) 115161

Fig. 15. Effect of structural parameters on equivalent stiffnesses of CSP-HLC.

analysis under SSSS boundary condition. Table 5 list the geometric and the changes of the equivalent tensile stiffness was more significant.
parameters of unit cell used for parametric study, and Fig. 10 ∼ Fig. 13 This was because the equivalent tensile stiffness was directly propor-
show the geometric configuration of unit cell with different structural tional to the cross-sectional area, which was positively correlated with
parameters. the lattice rod radius, as shown in Fig. 12. The equivalent bending
stiffness was directly proportional to the moment of inertia, which was
4.1.1. Effect of structural parameters on equivalent stiffnesses also positively correlated with the lattice rod radius.
Fig. 14(a) shows the effect of 𝑡𝑓 ∕ℎ𝑐 on the equivalent stiffnesses Fig. 14(d) shows the effect of the aspect ratio on the equivalent
of the CSP-HLC when the other parameters remained unchanged. The stiffnesses when the other parameters remained unchanged. The larger
equivalent stiffnesses increased linearly with the increase in this ratio, the aspect ratio was, the smaller the equivalent stiffnesses became. The
in which the equivalent tensile stiffnesses 𝐴11 and 𝐴22 and equivalent reason was that the suspended area of the face sheet not supported
bending stiffnesses 𝐷11 and 𝐷22 changed significantly. The reason was by the lattice core increased with increasing aspect ratio, as shown
that the equivalent tensile stiffness was directly proportional to the in Fig. 13, resulting in a downward trend in the equivalent tensile
cross-sectional area, which was linearly proportional to the ratio of stiffness. The equivalent bending stiffness was directly proportional to
𝑡𝑓 ∕ℎ𝑐 , as shown in Fig. 10. The equivalent bending stiffness was directly the moment of inertia, which decreased with the increase in the aspect
proportional to the moment of inertia, which also changed linearly with ratio.
𝑡𝑓 ∕ℎ𝑐 .
Fig. 14(b) shows the effect of the inclination angle of the lattice 4.1.2. Effect of structural parameters on static displacement
rod on the equivalent stiffnesses when the other parameters remained To further study the influence of the structural parameters on the
unchanged. The equivalent stiffnesses decreased with the increase in effective performance of the CSP-HLC, the effects of the structural pa-
the inclination angle. The reason was that the cross-sectional area and rameters on the static displacements of the CSP-HLC were determined
moment of inertia became smaller with the increase in inclination and are shown in Fig. 15. The vertical displacement of the CSP-HLC
angle, as shown in Fig. 11, which further resulted in the decrease in decreased with the increase in 𝑡𝑓 ∕ℎ𝑐 and the lattice rod radius, while
the equivalent stiffnesses. However, the equivalent stiffnesses basically it increased with the increase in the inclination angle and aspect ratio,
remained unchanged after the inclination angle changed from 10◦ to which was consistent with the variation trend of the equivalent bending
11◦ , indicating that the inclination angle had the greatest impact on the stiffness. It should be noted that 𝑡𝑓 ∕ℎ𝑐 and the aspect ratio had the
equivalent stiffnesses when the endpoint of the inclined rod coincided greatest impact on the static displacements, and the inclination angle
with the corner point of the face sheet. of the lattice rod had the least impact. In addition to the significant
Fig. 14(c) shows the effect of the lattice rod radius on the equivalent variation trend of the vertical displacement, the displacements 𝑈1 and
stiffnesses when the other parameters remained unchanged. The equiv- 𝑈2 were almost unchanged, indicating that the equivalent bending
alent stiffnesses increased with the increase in the lattice rod radius, stiffness had a significant impact on the vertical displacement.

12
J. Shi et al. Composite Structures 284 (2022) 115161

Fig. 16. Effect of structural parameters on equivalent stiffnesses of CSP-HLC.

4.1.3. Effect of structural parameters on natural frequencies of CSP-HLC and 𝐷11 showed nonlinear downward trends and 𝐴22 and 𝐷22 showed
The effects of the structural parameters on the first three natural nonlinear upward trends with the gradual increase in the ply angle,
frequencies of the CSP-HLC are shown in Fig. 16. The first three natural while the other stiffnesses showed a cyclic trend of increasing and
frequencies of the CSP-HLC were positively correlated with 𝑡𝑓 ∕ℎ𝑐 and decreasing. Thus, the layup configuration had a significant impact on
the inclination angle of the lattice rod and negatively correlated with the equivalent stiffnesses 𝐴11 , 𝐴22 , 𝐷11 , and 𝐷22 .
the lattice rod radius and the aspect ratio. The equivalent density Fig. 18(b) shows the first three natural frequencies and static dis-
increased with 𝑡𝑓 ∕ℎ𝑐 and the lattice rod radius and decreased with placements of the CSP-HLC with different layup configurations of the
the increase in the inclination angle and aspect ratio. 𝑡𝑓 ∕ℎ𝑐 and the face sheet. The static displacements decreased with the increase in the
lattice rod radius had the greatest influence on the natural frequencies. ply angle, while the influence of the layup configurations on the nat-
It should be noted that the influence of the equivalent density on the ural frequencies was irregular. The first and third natural frequencies
natural frequencies was not consistent in the analysis of the four struc- were positively correlated with the ply angle, and the second natural
tural parameters, because the natural frequencies depended on both the frequency increased first and then decreased with increasing ply angle.
equivalent density and the equivalent stiffness in a reciprocal manner.
For example, the extra equivalent stiffness influenced the frequencies 5. Comparative analysis of CSP-HLC and composite sandwich plate
more predominantly than the additional equivalent density with the with pyramidal lattice cores
increase in 𝑡𝑓 ∕ℎ𝑐 , so the natural frequencies of the plate increased.
In this section, the effects of different lattice forms (pyramidal and
4.2. Material parameters hourglass lattices) on the static and dynamic behaviors of composite
lattice sandwich plates were investigated. The schematic diagram of the
Due to the anisotropy and inhomogeneity of the composites, the two lattice forms are shown in Fig. 19. The materials and dimensions
layup configurations will affect the effective properties of the face sheet used in the composite sandwich plate with pyramidal lattice cores (CSP-
and further affect the effective plate properties of the CSP-HLC. In PLC) were the same as those of the CSP-HLC. The lattice rod radius
this section, the influence of the five incremental layup configurations was 1.5 mm, the thickness of the face sheet was 1 mm, the core height
(ID1: [±15◦ ∕0◦ ∕90◦ ]𝑠 , ID2: [±30◦ ∕0◦ ∕90◦ ]𝑠 , ID3: [±45◦ ∕0◦ ∕90◦ ]𝑠 , ID4: was 10 mm, and the dimensions of the unit cell were length × width
[±60◦ ∕0◦ ∕90◦ ]𝑠 , and ID5: [±75◦ ∕0◦ ∕90◦ ]𝑠 ) shown in Fig. 17 on the natu- = 40 mm × 40 mm. The dimensions of the CLSP-PLC obtained using a
ral frequencies and static displacements of the plate were investigated, 14 × 10 unit cell array were the same as those of the CSP-HLC (560 mm
in which the ply angles of the first two layers gradually increased from × 400 mm).
ID1 to ID5. The vertical displacements of the CSP-HLC and CSP-PLC along the
Fig. 18(a) shows the effect of the layup configurations on the centerline under SSSS boundary conditions and a uniform load of
equivalent stiffnesses of the CSP-HLC. The equivalent stiffnesses 𝐴11 0.5 MPa are compared in Fig. 20. The vertical displacement of the

13
J. Shi et al. Composite Structures 284 (2022) 115161

Fig. 17. Five incremental layup configurations of composite face sheets.

Fig. 18. Effect of layup configurations on equivalent stiffnesses, natural frequencies, and static displacements of CSP-HLC.

CSP-HLC (9.407 mm) was slightly less than that of CSP-PLC (9.59 mm) estimate the static displacements and natural frequencies of the CSP-
under the same uniformly distributed load, indicating that the CSP-HLC HLC quickly and accurately. By studying the influence of the material
had a greater equivalent bending stiffness, as expected. and structural parameters on the effective performance of the CSP-HLC,
Fig. 21 shows the load–displacement curves of these two composite the following conclusions can be drawn.
sandwich plates with increasing uniform distribution load in the elastic
(1) Increasing the ratio of the face sheet thickness to the core height
stage. The vertical displacements of the two composite sandwich plates
were directly proportional to the uniformly distributed load, and the and the lattice rod radius would increase the equivalent stiffnesses of
variation trends were basically the same, indicating that the static the CSP-HLC, while increasing the inclination angle of the lattice rod
mechanical properties of the CSP-HLC and CSP-PLC were similar. and the aspect ratio would reduce the equivalent stiffnesses.
Table 6 compares the first four natural frequencies of the CSP-HLC (2) The ratio of the face sheet thickness to the core height had
and CSP-PLC under SSSS boundary conditions. The mode shape of each the greatest influence on the static displacements of the CSP-HLC.
order predicted by the 2D-RPM and 3D-DNS were basically the same, The natural frequencies of the CSP-HLC were closely related to the
but the natural frequencies of the CSP-PLC were greater than those of
equivalent stiffnesses and equivalent density, and the variation of the
the CSP-HLC from the first order to the fourth order, which may have
lattice rod radius had the greatest influence on its natural frequencies.
been because the additional equivalent density of the CSP-HLC shown
in Fig. 19 influenced the frequency more predominantly than the extra (3) Different layup configurations had significant effects on the
equivalent stiffness of the CSP-PLC. natural frequencies, static displacements, and equivalent stiffnesses of
the CSP-HLC. The variation trends did were not regular but were
6. Conclusions related to the equivalent stiffnesses. The static displacements of the
CSP-HLC were smaller than those of the CSP-PLC, while the natural fre-
Based on the VAM, a reduced-order plate model was established quencies were similar to those of the CSP-PLC. This indicated that the
to analyze the static and dynamic behaviors of the CSP-HLC. The additional equivalent density of the CSP-HLC influenced the frequency
calculation results showed that the proposed model could be used to more predominantly than the extra equivalent stiffnesses.

14
J. Shi et al. Composite Structures 284 (2022) 115161

Fig. 19. Schematic diagram of unit cell with hourglass lattice core and pyramidal lattice core.

Table 6
Comparison of first four natural frequencies (Hz) and mode shapes of CSP-HLC and CSP-PLC under SSSS boundary conditions (values in red
brackets denote relative errors).
CSP-HLC CSP-PLC
Mode 3D-DNS 2D-RPM 3D-DNS 2D-RPM

15
J. Shi et al. Composite Structures 284 (2022) 115161

[2] Chen Y, Zhang L, Zhao Y, et al. Mechanical behaviors of C/SiC pyrami-

dal lattice core sandwich panel under in-plane compression. Compos Struct
2019;214(4):103–13.
[3] Gdx A, Tao Z, Su CA, et al. Free vibration of composite sandwich beam with
graded corrugated lattice core. Compos Struct 2019;229:111466.
[4] Sebaey TA, Mahdi E. Crushing behavior of a unit cell of CFRP lattice core for
sandwich structures’ application. Thin-Walled Struct 2017;116:91–5.
[5] Wang J, Evans AG, Dharmasena K, et al. On the performance of truss panels
with kagome cores. Int J Solids Struct 2003;40(25):6981–8.
[6] Hyun S, Karlsson AM, Torquato S, et al. Simulated properties of kagome and
tetragonal truss core panels. Int J Solids Struct 2003;40(25):6989–98.
[7] Hutchinson RG, Wicks N, Evans AG, et al. Kagome plate structures for actuation.
Int J Solids Struct 2003;40(25).
[8] Feng L, Wu L, Yu GC. An hourglass truss lattice structure and its mechanical
performances. Mater Des 2016;99:581–91.
[9] Sayyad AS, Ghugal YM. Bending, buckling and free vibration of laminated
composite and sandwich beams: A critical review of literature. Compos Struct
2017;171(7):486–504.
[10] Vigliotti A, Deshpande VS, Pasini D. Non linear constitutive models for lattice
materials. J Mech Phys Solids 2014;64(mar.):44–60.
Fig. 20. Comparison of vertical displacements of CSP-HLC and CSP-PLC under SSSS [11] Cohen N, Mcmeeking RM, Begley MR. Modeling the non-linear elastic response
boundary conditions and 0.5 MPa uniform load. of periodic lattice materials. Mech Mater 2018;129(1):159–68.
[12] Abrate S, Sciuva M Di. Equivalent single layer theories for composite and
sandwich structures: A review. Compos Struct 2017;179(11):482–94.
[13] Quan TQ, Duc ND. Nonlinear vibration and dynamic response of shear de-
formable imperfect functionally graded double curved shallow shells resting on
elastic foundations in thermal environment. J Therm Stresses 2016;39(4):437–59.
[14] Duc ND, Quan TQ. Transient responses of functionally graded double curved
shallow shells with temperature-dependent material properties in thermal
environment. Eur J Mech A Solids 2014;47:101–23.
[15] Duc ND, Tung HV. Nonlinear analysis of stability for functionally graded
cylindrical shells under axial compression. J Comput Mater Sci 2010;49(4):313–
6.
[16] Cong PH, Duc ND. Nonlinear dynamic analysis of porous eccentrically stiffened
double curved shallow auxetic shells in thermal environments. Thin-Walled
Struct 2021;163:107748.
[17] Quyen NV, Thanh NV, Quan TQ, Duc ND. Nonlinear forced vibration of sandwich
cylindrical panel with negative Poisson’s ratio auxetic honeycombs core and
CNTRC face sheets. Thin Walled Struct 2021;162:107571.
[18] Quan TQ, Quyen NV, Duc ND. An analytical approach for nonlinear thermo-
electro-elastic forced vibration of piezoelectric penta - graphene plates. Eur J
Mech A 2021;85:104095.
[19] Quan TQ, Anh VM, Vinyas M, Duc ND. Vibration and nonlinear dynamic response
Fig. 21. Comparison of load–displacement curves of CSP-HLC and CSP-PLC with of imperfect sandwich piezoelectric auxetic plate. J Mech Adv Mater Struct 2020.
increasing uniform distribution load in elastic stage. [20] Cong PH, Quyet PK, Duc ND. Effects of lattice stiffeners and blast load on
nonlinear dynamic response and vibration of auxetic honeycomb plates. Proc
Inst Mech Eng 2021.
[21] Duong NT, Duc ND. Evaluation of elastic properties and thermal expansion
CRediT authorship contribution statement coefficient of composites reinforced by randomly distributed spherical particles
with negative Poisson’s ratios. J Compos Struct 2016;153:569–77.
Shi Jie: Software, Methodology, Writing – review & editing. Zhong [22] Khakalo S, Balobanov V, Niiranen J. Modelling size-dependent bending, buck-
Yifeng: Writing – original draft, Supervision, Investigation. Liu Rong: ling and vibrations of 2D triangular lattices by strain gradient elasticity
models: Applications to sandwich beams and auxetics. Internat J Engrg Sci
Formal analysis, Investigation. Shi Zheng: Conceptualization, Supervi-
2018;127(6):33–52.
sion. [23] Khakalo S, Niiranen J. Anisotropic strain gradient thermoelasticity for cellular
structures: Plate models, homogenization and isogeometric analysis. J Mech Phys
Declaration of competing interest Solids 2019;134:103728.
[24] Goncalves BR, Karttunen A, Romanoff J, et al. Buckling and free vibration of
The authors declare that they have no known competing finan- shear-flexible sandwich beams using a couple-stress-based finite element. Compos
cial interests or personal relationships that could have appeared to Struct 2017;165:233–41.
[25] Karttunen AT, Reddy JN, Romanoff J. Two-scale micropolar plate model for
influence the work reported in this paper.
web-core sandwich panels. Int J Solids Struct 2019.
[26] Chowdhury SR, Reddy JN. Geometrically exact micropolar Timoshenko beam
Acknowledgments and its application in modelling sandwich beams made of architected lattice
core. Compos Struct 2019;226:111228.
The work is supported by the National Natural Science Founda- [27] Lai L. Study of free vibration of aluminium honeycomb panels. Toronto:
tion of China (Grant Numbers: 52073036, 51778088), and Natural University of Toronto; 2002, p. 102–30.
[28] Yu W, Qian M, Li H. Elastic and plastic properties of epoxy resin syntactic
Science Foundation of Chongqing, China (Grant Number: cstc2021jcyj-
foams filled with hollow glass microspheres and glass fibers. J Appl Polym Sci
msxmX0035).
2016;133(46):385–97.
[29] Zhong Y, Qin W, Yu W, et al. Variational asymptotic homogeniza-
References tion of magneto-electro-elastic materials with coated fibers. Compos Struct
2015;133:300–11.
[1] Mcshane GJ, Radford DD, Deshpande VS, et al. The response of clamped [30] Yifeng Z, Lei C, Yu W, et al. Variational asymptotic micromechanics mod-
sandwich plates with lattice cores subjected to shock loading. Eur J Mech A eling of heterogeneous magnetostrictive composite materials. Compos Struct
Solids 2006;25(2):215–29. 2013;106(12):502–9.

16
J. Shi et al. Composite Structures 284 (2022) 115161

[31] Yifeng Z, Lei C, Yu W. Variational asymptotic modeling of the thermomechanical [33] Peng X, Zhong YF, Shi Z, et al. Estimation of effective properties of composite
behavior of composite cylindrical shells. Compos Struct 2012;94(3):1023–31. sandwich panels with negative Poisson’s ratio by using variational asymptotic
[32] Cesnik CE, Hodges DH. VABS: A new concept for composite rotor blade multiscale method. Mater Today Commun 2020;23:101072.
cross-sectional modeling. J Am Helicopter Soc 1997;42(1):27–38.

17