Static and Free-Boundary Vibration Analysis of Egg-Crate Honeycomb Core Sandwich Panels Using the VAM-Based Equivalent Model

1. Introduction

Sandwich panels have become indispensable in modern engineering due to their exceptional combination of high stiffness, strength, and low weight [1,2,3]. A typical sandwich structure consists of two thin, stiff face sheets bonded to a lightweight core, which can occupy more than half the volume and weight of the panel [4,5]. The core thus plays a critical role in mechanical performance, and extensive research has been devoted to optimizing core configurations for improved strength-to-weight ratios [6,7]. Early studies focused on honeycomb cores – named for their biomimetic hexagonal cells – which endow panels with high out-of-plane compression and shear stiffness at minimal weight [8,9].

Recently, numerous analytical models have been developed to characterize sandwich panels, ranging from classical first-order shear deformation theory [10,11] to higher-order and layer-wise (zig-zag) theories [12,13]. These models improve accuracy in predicting bending and vibration by accounting for transverse shear and layer discontinuities, but they also introduce complexity (e.g. requiring additional degrees of freedom for each layer) and may not easily accommodate complex core architectures [14,15]. As a result, the analysis and design of novel core geometries often rely on numerical methods. Finite element analysis (FEA) is frequently used in both industry and academia to evaluate sandwich panel behavior, yet full three-dimensional (3D) modeling of a sandwich structure with intricate cellular cores is computationally expensive and time-consuming [16,17].

To alleviate this issue, homogenization and equivalent single-layer approaches have been explored, wherein a periodic unit cell of the core is analyzed to derive effective properties for an equivalent continuum panel [18,19]. Such micromechanical approaches greatly reduce computational cost by replacing the detailed core with a homogeneous plate model [20,21]. However, standard homogenization can lose microstructural detail and may sacrifice accuracy in predicting local stress concentrations or higher-order effects, especially for advanced core designs [22]. This limitation motivates the pursuit of more refined yet efficient modeling techniques.

In parallel with advances in modeling, a rich variety of innovative core designs has emerged to overcome the limitations of conventional honeycombs and lattice grids [23,24]. Researchers have proposed architected cellular materials and metamaterial cores that achieve superior mechanical performance through geometry [25,26]. For example, auxetic re-entrant honeycombs (with negative Poisson’s ratio) have demonstrated enhanced energy absorption and stiffness compared to regular hexagonal cells [27,28]. Novel re-entrant chiral and hybrid honeycomb patterns have been shown to improve impact resistance and specific energy absorption under compression [29,30]. Elliptical and other non-periodic cell geometries have likewise been tailored to enhance in-plane stiffness or to achieve desirable anisotropic properties [31]. Beyond honeycombs, lattice and grid-like cores manufactured via advanced methods (e.g. additive manufacturing) are gaining attention [32,33].

Periodic truss lattices and orthogrid structures can be optimized to align with load paths, and they offer the promise of multi-functionality – for instance, open-cell lattice cores can double as heat exchange surfaces for thermal management [34,35]. Despite these advantages, conventional grid and lattice sandwich structures face several limitations. First, purely open lattice cores are prone to local instabilities: under concentrated or transverse loads, their slender struts can lead to highly localized stress and early failure in the core [36,37]. This makes lattice cores less effective in carrying impact loads without additional reinforcement.

Second, while grid-stiffened panels (e.g. orthogrid or iso-grid structures) provide high in-plane stiffness, they often exhibit complex buckling modes and stress concentrations at nodal intersections or cut-outs. Orthogonal grid patterns increase stiffness but may require curved or hierarchical modifications to delay buckling in critical regions [38]. In general, fully exploiting these novel core architectures is challenging because their mechanical behavior is not as well understood as that of traditional honeycombs. Designers lack simple formulas or design guidelines, and instead must resort to extensive FEA for each new topology [39,40]. This trial-and-error approach is costly and time-consuming, underscoring the need for efficient analytical or equivalent models to guide the design of advanced sandwich cores.

The present work introduces a novel “egg-crate” honeycomb core sandwich panel (SP-EHC), along with a high-fidelity yet computationally efficient modeling framework for its static and dynamic analysis. The SP-EHC core concept is derived from conventional lattice and grid structures, aiming to integrate their strengths while mitigating weaknesses. The core consists of an orthogonal grid network with intersecting walls, but unlike a simple rectangular grid, the intersection regions are chamfered to form a semi-enclosed cell geometry as shown in Figure 1. This egg-crate-like configuration creates a pattern of contiguous cells analogous to a honeycomb, but with open corners and pathways. The chamfered grid walls serve to enhance structural strength and buckling resistance by eliminating sharp corner stress concentrations and providing smoother load transfer between intersecting ribs. Meanwhile, the semi-enclosed cells preserve the panel’s lightweight characteristics and introduce practical benefits such as ventilation and drainage through the core.

To rigorously analyze the SP-EHC design, the two-dimensional equivalent plate model (2D-EPM) using the variational asymptotic method (VAM) is developed. The VAM, originally formulated by Cesnik and Hodges for composite structures [41,42], is a multiscale asymptotic technique that reduces a 3D elasticity problem to a lower-dimensional model without ad-hoc assumptions. The core idea of VAM is to exploit the small geometric parameter (the core thickness-to-panel length ratio, in this case) and systematically minimize the 3D strain energy functional to obtain a simplified model that still captures the essential physics [43]. In practice, this involves decomposing the displacement field into a 2D (mid-plane) part plus a warping field, and deriving effective plate stiffness by averaging the 3D behavior over a representative unit cell. Unlike classical plate theories, the VAM-based approach inherently accounts for the actual microstructure of the core, enabling accurate prediction of both global responses and local fields.

This method has been successfully applied to periodic lattice and honeycomb structures in recent studies, showing that a VAM-derived plate model can reproduce the full 3D behavior with dramatically reduced computational cost [44,45]. For example, Wang et al. [46] used a VAM-based equivalent model for triangular honeycomb panels and achieved excellent agreement with experimental and 3D FEA results (for bending deflections, buckling loads, and vibration frequencies) while improving computation efficiency by orders of magnitude. In this study, the Variational Asymptotic Method (VAM) is employed to derive the three-dimensional strain energy of the SP-EHC cell, followed by an asymptotic expansion to extract the governing two-dimensional plate equations and construct the equivalent stiffness matrices. Based on the 2D solution, recovery relations are established to reconstruct local stress and strain fields within the core. A key advantage of the VAM framework lies in its multiscale capability: it enables efficient computation of global deflections and stresses, while also allowing detailed reconstruction of micro-scale stress distributions—such as those in cell walls—for failure and fatigue analysis. This level of resolution is not readily achievable through conventional homogenization techniques.

2. VAM-Based Equivalent Model of SP-EHC

2.1. Three-Dimensional Energy Expression of SP-EHC

From the perspective of the equivalent single-layer theory, the 3D displacement field $u_i$ of SP-EHC can be expressed using the 2D displacement field ${\overline{u}}_i$ and the warping functions ${\chi }_i$, as

$$\begin{array}{cc}\hfill u_1(x_{\alpha },y_i)& =\underline{{\overline{u}}_1\left(x_{\alpha }\right)-\langle y_3\rangle {\overline{u}}_{3,1}\left(x_{\alpha }\right)}+{\chi }_1(x_{\alpha };y_i)\hfill \\ \hfill u_2(x_{\alpha },y_i)& =\underline{{\overline{u}}_2\left(x_{\alpha }\right)-\langle y_3\rangle {\overline{u}}_{3,2}\left(x_{\alpha }\right)}+{\chi }_2(x_{\alpha };y_i)\hfill \\ \hfill u_3(x_{\alpha },y_i)& =\underline{{\overline{u}}_3\left(x_{\alpha }\right)}+{\chi }_3(x_{\alpha };y_i)\hfill \end{array}$$

(1)

where the underlined terms represent the displacements of classical plate theory.

The explicit form of $u_i$ can be derived from Eq. (1) by applying a volume-averaging procedure over the representative unit cell:

$${\overline{u}}_1=\langle u_1\rangle +\langle y_3\rangle {\overline{u}}_{3,1},{\overline{u}}_2=\langle u_2\rangle +\langle y_3\rangle {\overline{u}}_{3,2},{\overline{u}}_3=\langle u_3\rangle$$

(2)

Considering that the local coordinate system originates from an arbitrary point within the unit cell, the 3D displacement field can be projected to obtain an equivalent 2D displacement field using Eq. (2). To ensure consistency with the macroscopic kinematics, the warping functions are subject to the following constraint:

$$\langle {\chi }_i\rangle =0$$

(3)

By introducing the concept of strain decomposition with local rotation, the 3D strain field accounting for small rotations can be expressed as:

$${\epsilon }_{ij}={\displaystyle \frac{1}{2}}\left({\displaystyle \frac{\partial u_i}{\partial x_j}}+{\displaystyle \frac{\partial u_j}{\partial x_i}}\right)-{\delta }_{ij},$$

(4)

where ${\delta }_{ij}$ is the Kronecker delta symbol.

To derive the explicit expression of the 3D strain field, Eq. (1) is substituted into Eq. (4). Neglecting higher-order terms that do not contribute to the total strain energy yields:

$$\begin{array}{cc}\hfill {\epsilon }_{11}& ={\overline{\epsilon }}_{11}+\xi y_3{k}_{11}+x_1,1+{\displaystyle \frac{1}{\xi }}{\chi }_{11},\hfill \\ \hfill 2{\epsilon }_{12}& =2{\overline{\epsilon }}_{12}+\xi y_3(k_{12}+k_{21})+x_1,2+x_2,1+{\displaystyle \frac{1}{\xi }}(x_{211}+x_{12}),\hfill \\ \hfill {\epsilon }_{22}& ={\overline{\epsilon }}_{22}+\xi y_3{k}_{22}+{\chi }_{2,2}+{\displaystyle \frac{1}{\xi }}{\chi }_{212},\hfill \\ \hfill 2{\epsilon }_{13}& =x_{113}+{\chi }_{3,1}+{\displaystyle \frac{1}{\xi }}{\chi }_{313},\hfill \\ \hfill 2{\epsilon }_{23}& =x_{213}+{\chi }_{3,2}+{\displaystyle \frac{1}{\xi }}{\chi }_{312},\hfill \\ \hfill {\epsilon }_{33}& ={\chi }_{313}\hfill \end{array}$$

(5)

where

$${\epsilon }_{\alpha \beta }(x_1,x_2)={\displaystyle \frac{1}{2}}\left({\overline{u}}_{\alpha ,\beta }+{\overline{u}}_{\beta ,\alpha }\right),k_{\alpha \beta }(x_1,x_2)=-{\overline{u}}_{3,\alpha \beta }$$

(6)

The strain energy of the SP-EHC can be expressed as:

$$U={\int }_{x_1}{\int }_{x_2}{\displaystyle \frac{U_{\Omega }}{\Omega }}d{x}_{1}d{x}_2,$$

(7)

where ${\displaystyle \frac{U_{\Omega }}{\Omega }}$ denotes the strain energy density, and $\Omega$ is the representative unit cell domain in the $x_1$-$x_2$ plane.

The explicit form of $U_{\Omega }$ is given by:

$$\begin{array}{cc}\hfill U_{\Omega }=& {\int }_{-h_2/2}^{t_f}{\int }_{L_2}^{L_2}{\int }_{L_2}^{L_2}{\sigma }_t{\epsilon }_{t}d{y}_{1}d{y}_{2}d{y}_3+{\int }_{-h_2/2}^{h_2/2}2{\int }_{L_1/2}^{L_1/2}{\sigma }_b{\epsilon }_{b}d{y}_{1}d{y}_2\hfill \\ \hfill & +4\times ({\int }_{h_1/2}^{h_2+b_1}{\int }_{-t/2}^0{\int }_{b_1}^{b_2}{\sigma }_{\mathrm{I}}{\epsilon }_{\mathrm{I}}d{y}_{1}d{y}_{2}d{y}_3-{\int }_{h_1/2}^{h_2+b_2}{\int }_{-t/2}^{t/2}{\sigma }_{\mathrm{II}}{\epsilon }_{\mathrm{II}}d{y}_{1}d{y}_{2}d{y}_3\hfill \\ \hfill & -{\displaystyle \frac{1}{2}}{\int }_{h_1/2+b_2}^{h_2+b_2}{\int }_{-t/2}^{t/2}{\sigma }_{\mathrm{III}}{\epsilon }_{\mathrm{III}}d{y}_{1}d{y}_{2}d{y}_3)\hfill \end{array},$$

(8)

where the subscript b and t refer to the bottom and top face sheets, respectively; subscripts I, II, and III denote distinct regions within the core cell, as illustrated in Figure 2.

The virtual work $W_{3D}$ done by external loads on the SP-EHC can be decomposed into two parts: the virtual work $W_{2D}$ related to the 2D-EPM, and the remaining part $W^{*}$, such that

$$\begin{array}{cc}\hfill W_{3D}& =W_{2D}+W^{*}\hfill \\ \hfill & ={\int }_s(p_i{\overline{u}}_i+q_{\alpha }\delta u_3,\alpha )ds+{\int }_s\left(\langle f_i{\chi }_i\rangle +{\tau }_i{\chi }_i^{+}-{\beta }_i{\chi }_i^{-}\right)ds\hfill \end{array}$$

(9)

where ${\chi }_i^{-}$ and ${\chi }_i^{+}$denote the warping displacements on the bottom and top surfaces, respectively; $p_i$ and ${\tau }_i$ represent the surface tractions on the bottom and top surfaces; $f_i$ is the body force. The distributed force contributions to the virtual work are expressed as:

$$p_i=\left(f_i\right)+{\tau }_i+{\beta }_i,q_{\alpha }={\displaystyle \frac{h}{2}}\left({\beta }_{\alpha }-{\tau }_{\alpha }\right)-〈y_3{f}_{\alpha }〉$$

(10)

The original 3D kinetic energy of SP-EHC can be expressed as:

$$K_{3D}=K_{2D}+K^{*}$$

(11)

where

$$K_{2D}={\displaystyle \frac{1}{2}}{\int }_s\left(\rho {\overline{V}}^T\overline{V}+2{\alpha }^T\rho \xi \overline{V}+{\alpha }^T\Phi \alpha \right)ds$$

(12)

$$K^{*}={\displaystyle \frac{1}{2}}{\int }_{\nu }\rho \left[{(\dot{\chi }+\dot{\chi })}^T(\alpha \xi +\dot{\chi })+2{(V+\alpha \xi )}^T(\alpha \chi +\dot{\chi })\right]dv$$

(13)

where V and $\alpha$ denote the absolute velocity and angular velocity, respectively, and v is the volume occupied by the sandwich panel. The corresponding quantities are defined as:

$$\dot{x}=\partial x/\partial t,\xi ={\left[00x_3\right]}^T,\rho \xi ={\left[00\langle x_3\rho \rangle \right]}^T,\Phi =\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& \langle x_3^2\rho \rangle \end{array}\right]$$

(14)

The parameters $\rho$ and $\overline{\rho }$ represent the mass densities of 3D-FEM and 2D-EPM, respectively. The equivalent mass density $\overline{\rho }$ can be determined by assuming ideal face sheet-core bonding, while neglecting the mass contribution of internal walls within the cellular structure, such as

$$\overline{\rho }={\displaystyle \frac{m_c+m_f}{h·l_1^2}}={\displaystyle \frac{4{\rho }_c{V}_c+2{\rho }_f{V}_f}{(h_1+2t_f)·l_1^2}}$$

(15)

where ${\rho }_c$ and ${\rho }_f$ are the mass densities of the core and face sheets, respectively; $V_f=l_1^2·t_f$ represents the volume of face sheet. Based on geometric modeling, the core volume $V_c$ can be expressed as:

$$V_c=\left[{\displaystyle \frac{l_1}{2}}\times h_1-b_2(h_1-h_2)-{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{l_1}{2}}-b_2-{\displaystyle \frac{b_1}{2}}\right)(h_1-h_2)\right]·t$$

(16)

The elastic dynamic behavior of SP-EHC satisfies Hamilton’s principle, which can be stated as:

$${\int }_{t_1}^{t_2}\left[\delta (K_{2D}+K^{*}-U)+\delta W_{2D}+\delta W^{*}\right]dt=0$$

(17)

The governing equation for the variation of ${\chi }_i$ can be expressed as:

$$\underset{{\chi }_i\in Eq.\left(3\right)}{min}\langle {\epsilon }^{T}D\epsilon \rangle$$

(18)

where D is the elasticity matrix and $\epsilon$ is the strain vector.

Directly solving the minimization problem in Eq. (18) poses significant challenges, particularly in the context of fully 3D elastic analyses. However, by conducting an asymptotic analysis of the transformed expressions within the equivalent formulation, a potential solution for ${\chi }_i$ can be derived without imposing restrictive assumptions on the displacement field or material behavior. The detailed procedure for this asymptotic approach will be presented in the following section.

2.2. Equivalent Analysis Based on the VAM

To accurately reduce the original 3D problem to a 2D equivalent model, it is essential to recover the strain energy stored in the 3D structure within the 2D formulation. This is most effectively achieved by leveraging small geometric or material parameters inherent to the structure and applying a stepwise transformation via asymptotic analysis. Within this framework, the orders of magnitude of all quantities involved in the expansion can be expressed as:

$${\epsilon }_{\alpha \beta }\sim h{\kappa }_{\alpha \beta }\sim \eta ,f_3\sim \mu {\left({\displaystyle \frac{h}{L_1}}\right)}^2\eta ,f_{\alpha }\sim \mu \left({\displaystyle \frac{h}{L_1}}\right)\eta ,m_{\alpha }\sim \mu h\left({\displaystyle \frac{h}{L_1}}\right)\eta$$

(19)

where $\mu$ and h represent the characteristic material property scale and geometric scale, respectively; $\eta$ is a small perturbation parameter.

By analyzing the scaling orders in Eq. (19), it becomes evident that the applied loads contribute at higher orders. As a result, terms corresponding to zero-order energy equivalence can be reasonably neglected. Referring to Eq. (5), the strain field under the zero-order approximation can be expressed in matrix form as:

$$\epsilon ={\Gamma }_{\chi }\chi +{\Gamma }_{\epsilon }\epsilon$$

(20)

where $\chi ={\left[\begin{array}{ccc}{\chi }_1& {\chi }_2& {\chi }_3\end{array}\right]}^T,\epsilon ={\left[\begin{array}{cccccccc}{\epsilon }_{11}& 2{\epsilon }_{12}& {\epsilon }_{22}& {\kappa }_{11}& 2{\kappa }_{12}& {\kappa }_{22}& 2{\gamma }_{13}& 2{\gamma }_{23}\end{array}\right]}^T$, and

$${\Gamma }_{\chi }=\left[\begin{array}{ccc}{\displaystyle \frac{\partial }{\partial y_1}}& 0& 0\\ {\displaystyle \frac{\partial }{\partial y_2}}& {\displaystyle \frac{\partial }{\partial y_1}}& 0\\ 0& {\displaystyle \frac{\partial }{\partial y_2}}& 0\\ {\displaystyle \frac{\partial }{\partial x_3}}& 0& {\displaystyle \frac{\partial }{\partial y_1}}\\ 0& {\displaystyle \frac{\partial }{\partial x_3}}& {\displaystyle \frac{\partial }{\partial y_2}}\\ 0& 0& {\displaystyle \frac{\partial }{\partial x_3}}\end{array}\right],{\Gamma }_{\epsilon }=\left[\begin{array}{ccccccccc}1& 0& 0& \zeta y_3& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& \zeta y_3& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& \zeta y_3& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 1\end{array}\right]$$

(21)

To obtain a more general approximation, the warping functions ${\chi }_i\left(x_1,x_2;y_1,y_2,y_3\right)$ are discretized using a numerical method (such as finite element interpolation) as:

$$\chi (x_1,x_2;y_1,y_2,y_3)=S(y_1,y_2,y_3)N(x_1,x_2)$$

(22)

where S is the shape function matrix depending on microscopic coordinates; N is the shape function matrix formed by the nodal values of ${\chi }_i$.

By substituting Eqs. (22) into (9) and (20), the weak form of the zero-order approximate energy can be expressed as:

$${\Pi }_{\Omega }={\displaystyle \frac{1}{2\Omega }}\left(N^T{D}_{h\chi }N+2N^T{D}_{h\epsilon }\epsilon +{\epsilon }^T{D}_{\epsilon }\epsilon \right)$$

(23)

where

$$\begin{array}{cc}\hfill D_{h\chi }& ={\int }_{\Omega }{\left({\Gamma }_{h}S\right)}^{T}D\left({\Gamma }_{h}S\right)d\Omega ,\hfill \\ \hfill D_{h\epsilon }& ={\int }_{\Omega }{\left({\Gamma }_{h}S\right)}^{T}D{\Gamma }_{\epsilon }d\Omega ,\hfill \\ \hfill D_{\epsilon }& ={\int }_{\Omega }{\Gamma }_{\epsilon }^{T}D{\Gamma }_{\epsilon }d\Omega .\hfill \end{array}$$

(24)

By applying the averaging constraint condition in Eq. (3) and minimizing Eq. (23), a linear system is obtained as:

$$D_{h\chi }N=-D_{h\epsilon }\epsilon$$

(25)

Since N depends linearly on the 2D plate strain $\epsilon$, and this relationship is independent of the strain magnitude, the solution for the warping parameters can be symbolically expressed as:

$$N=N_0\epsilon$$

(26)

Substituting Eq. (26) into Eq. (23) yields the final form of the strain energy density stored within the unit cell:

$${\Pi }_{\Omega }={\displaystyle \frac{1}{2\Omega }}{\epsilon }^T\left(N_0^T{D}_{h\epsilon }+D_{\epsilon }\right)\epsilon ={\displaystyle \frac{1}{2}}{\epsilon }^T{D}_E\epsilon$$

(27)

where $D_E$ denotes the effective stiffness matrix, characterizing the equivalent elastic properties of the unit cell that incorporates complex microstructural materials and geometry

2.3. Transformation into the Reissner–Mindlin Model

Although the zero-order asymptotic strain energy correction in Eq. (27) is mathematically valid, an additional step is required to map the current approximate energy expression to the classical Reissner–Mindlin plate model. The generalized form of the Reissner–Mindlin strain energy is given by:

$$2{\Pi }_R=R^{T}XR+{\gamma }^{T}G\gamma +2R^T{F}_R$$

(28)

In Eq. (28), the Reissner–Mindlin model is obtained by eliminating the coupling terms between the membrane–bending strain vector $\epsilon ={\left[{\epsilon }_{11},2{\epsilon }_{12},{\epsilon }_{22},{\kappa }_{11},{\kappa }_{12},{\kappa }_{21},{\kappa }_{22}\right]}^T$ and the transverse shear strain vector $\gamma ={\left[2{\gamma }_{13},2{\gamma }_{23}\right]}^T$. The full strain energy expression can be rewritten as:

$$\begin{array}{c}2{\Pi }_a={\epsilon }^T\overline{D}\epsilon ={\epsilon }^T\overline{A}\epsilon +2{\overline{\epsilon }}^T\overline{B}\gamma +{\gamma }^T\overline{C}\gamma \\ ={\overline{\epsilon }}^T\left(\overline{A}-{\overline{BC}}^T{\overline{B}}^T\right)\overline{\epsilon }+{\left(\gamma +{\overline{C}}^{-1}{\overline{B}}^T\overline{\epsilon }\right)}^T\dot{C}\left(\gamma +{\overline{C}}^{-1}{\overline{B}}^2\right.\end{array}$$

(29)

In addition, the strain energy expression must include the work conjugate to the external load term F, which is expressed as:

$$F=-N_0^T\tau -N_0^T\beta -\langle N_0^{T}f\rangle$$

(30)

The total energy of the 2D-EPM can be expressed using the variable transformations $R=\overline{\epsilon }$ and $\overline{\gamma }=\gamma +C^{-1}{B}^T\overline{\epsilon }$, leading to:

$$2{\Pi }_R=R^{T}XR+{\overline{\gamma }}^{T}G\overline{\gamma }+2R^T{F}_R$$

(31)

where

$$X=A-B{C}^{-1}{B}^T,G=C,F_R=F$$

(32)

Due to the axial and planar symmetry of the SP-EHC, certain stiffness components are eliminated, leading to the final constitutive relation of the 2D-EPM:

$$\left[\begin{array}{c}{N}_{11}\\ N_{22}\\ N_{12}\\ M_{11}\\ M_{22}\\ M_{12}\\ Q_1\\ Q_2\end{array}\right]=\left[\begin{array}{cccccccc}{A}_{11}\hfill & A_{12}\hfill & 0\hfill & B_{11}\hfill & B_{12}\hfill & 0\hfill & 0\hfill & 0\hfill \\ \hfill & A_{22}\hfill & 0\hfill & B_{12}\hfill & B_{22}\hfill & 0\hfill & 0\hfill & 0\hfill \\ \hfill & \hfill & A_{66}\hfill & 0\hfill & 0\hfill & B_{66}\hfill & 0\hfill & 0\hfill \\ \hfill & \hfill & & D_{11}\hfill & D_{12}\hfill & 0\hfill & 0\hfill & 0\hfill \\ \hfill & \mathrm{sym}\hfill & \hfill & & D_{22}\hfill & 0\hfill & 0\hfill & 0\hfill \\ \hfill & \hfill & & \hfill & & D_{66}\hfill & 0\hfill & 0\hfill \\ \hfill & \hfill & & \hfill & & \hfill & C_{11}\hfill & 0\hfill \\ \hfill & \hfill & & \hfill & & \hfill & & C_{22}\hfill \end{array}\right]\left[\begin{array}{c}{\epsilon }_{11}\\ {\epsilon }_{22}\\ 2{\epsilon }_{12}\\ {\kappa }_{11}\\ {\kappa }_{22}\\ 2{\kappa }_{12}\\ 2{\gamma }_{13}\\ 2{\gamma }_{23}\end{array}\right]$$

(33)

where $N_{\alpha \beta }$, $M_{\alpha \beta }$, and $Q_{\alpha }$ ($\alpha ,\beta =1,2$) denote the in-plane forces, bending moments, and shear force, respectively.

The principal advantage of this method lies in its ability to construct the Reissner–Mindlin plate model directly from the zero-order asymptotic approximation, thereby eliminating the need for higher-order expansions or external transformation procedures. This approach enables the macroscopic behavior to be explicitly described through variable expressions derived from the asymptotic system, which inherently capture all 2D field variables in the in-plane coordinates $x_1$ and $x_2$. Consequently, the overall problem formulation is greatly simplified while preserving the essential mechanical characteristics of the original 3D structure.

2.4. Local Field Recovery

In addition to capturing macroscopic responses, evaluating the local field behavior within the original microstructure is often essential. Based on the results of the macroscopic structural analysis, the local displacement field can be reconstructed as

$$u_i={\overline{u}}_i+\left[\begin{array}{ccc}{\overline{u}}_{1,1}& {\overline{u}}_{1,2}& {\overline{u}}_{1,3}\\ {\overline{u}}_{2,1}& {\overline{u}}_{2,2}& {\overline{u}}_{2,3}\\ {\overline{u}}_{3,1}& {\overline{u}}_{3,2}& {\overline{u}}_{3,3}\end{array}\right]\left[\begin{array}{c}{y}_1\\ y_2\\ y_3\end{array}\right]+{\zeta }_i$$

(34)

The local strain field can be reconstructed as:

$$\begin{array}{cc}\hfill & {\epsilon }_e={\left[{\epsilon }_{11}{\epsilon }_{22}2{\epsilon }_{12}\right]}^T=\epsilon +\zeta y_{3}k\hfill \\ \hfill & 2{\xi }_s={\left[2{\epsilon }_{13}2{\epsilon }_{23}\right]}^T=-\left(\epsilon +\xi y_{3}k\right){\overline{C}}_{es}^{*}{C}_s^{-1}\hfill \\ \hfill & {\epsilon }_t={\epsilon }_{33}=-\left(\epsilon +\zeta y_{3}k\right){\overline{C}}_{et}^{*}{\overline{C}}_t^{-1}\hfill \end{array}$$

(35)

According to the 3D constitutive relationship (Hooke’s law), the local stress field within the unit cell can be obtained as:

$$\sigma =K\epsilon$$

(36)

In summary, the local stress and strain fields within the unit cell can be accurately reconstructed using the generalized macroscopic strain $\epsilon$, curvature $\kappa$, and displacement ${\overline{u}}_i$ within the equivalent model framework. This multiscale reconstruction approach ensures that detailed microscale mechanical responses are captured while leveraging the computational efficiency of the reduced two-dimensional formulation.

3. Dynamic Verification

In practical engineering applications, sandwich panels are often subjected to stochastic dynamic loads, including wind pressure, seismic excitation, and vehicular vibrations, which can critically affect their dynamic behavior, buckling resistance, and fatigue life. Consequently, an accurate assessment of their vibrational response across different modes is essential to ensure structural safety and reliability. However, full three-dimensional (3D) modeling of the SP-EHC entails significant computational effort due to the complexity of its geometry and material composition, thereby limiting analysis efficiency and model scalability. To address these challenges, a VAM-based equivalent model was employed to predict the free-boundary vibration behavior of the SP-EHC.

3.1. Material Property Testing

The mechanical properties of the 3D-printed PLA material were determined through uniaxial tensile tests conducted in accordance with the GB/T 1040.2-2006 standard using a universal testing machine. Dog-bone-shaped specimens, with dimensions shown in Figure 3(a), were carefully aligned and mounted to ensure vertical positioning and prevent slippage. A small preload was applied to eliminate potential seating effects. The tests were performed at a constant displacement rate of 5 mm/min until failure. Load and displacement were recorded in real time, while strain was measured using an extensometer to ensure high accuracy. Nominal stress–strain curves were obtained from the load–displacement data, as illustrated in Figure 3(c). The linear elastic region was identified by excluding initial seating and post-yield nonlinearity, and the elastic modulus was calculated as the slope of a least-squares linear fit over this range. The average elastic modulus, determined from three specimens, was 2,413 MPa. The Poisson’s ratio and density were measured separately, yielding values of 0.39 and 1.30 g/cm³, respectively.

3.2. Free-Boundary Vibration Analysis

In modal testing, boundary constraints are typically applied using mechanical fixtures; however, the elastic deformation of these fixtures can introduce additional compliance, adversely affecting the accuracy of measured modal frequencies and mode shapes. This deviation from ideal boundary conditions complicates high-precision modal analysis. To mitigate such boundary interference and ensure reliable validation of the VAM-based 2D-EPM, all-PLA specimens were fabricated using 3D printing technology. Modal vibration tests were conducted under free boundary conditions to eliminate fixture-induced artifacts. The specimens consisted of 8 and 4 unit cells in the $x_1$ and $x_2$ directions, respectively, with unit cell dimensions listed in Table 1. The overall dimensions of the specimen were L=360 mm, W=180 mm, and H=24 mm, as illustrated in Figure 4(a). This dense measurement configuration significantly improves the spatial resolution of modal data, thereby enhancing the accuracy of numerical model validation.

The free-boundary vibration test setup for the SP-EHC specimen is shown in Figure 4(c). To simulate free boundary conditions, the specimen was suspended from a gantry frame using two nylon cords. An accelerometer was mounted at a fixed location on the specimen surface to capture real-time response signals. The structure was excited using an impact hammer, and both force and acceleration time histories were recorded via a dynamic signal analyzer and transmitted to a computer for processing. A Fast Fourier Transform (FFT) was then applied to convert the time-domain data into the frequency domain, allowing extraction of the natural frequencies. In the numerical simulation, a frequency analysis step within the linear perturbation module was used to compute the first seven eigenfrequencies. The model was assigned completely free boundary conditions, and the interaction between the honeycomb core and face sheets was represented using “Tie” constraints. The three-dimensional model comprised 84,125 C3D10 elements, while the two-dimensional model employed 1,525 S4R elements.

Figure 5 shows seven distinct energy peaks, corresponding to the first through seventh modal frequencies of the SP-EHC under free boundary conditions, within the 310–430 Hz range. Each modal frequency was accurately identified using the half-power bandwidth method. The prominent peaks reflect resonance responses at specific frequencies, highlighting the structure’s characteristic modal behavior. Figure 6 shows the first seven elastic mode frequencies, which increase progressively with mode order. Notably, the seventh mode appears at 420.27 Hz, approximately 1.12 times the frequency of the first mode at 375.73 Hz.

A comparison of the modal frequencies obtained from the three methods indicates that the experimental results are the highest, followed by the 2D-EPM, while the 3D-FEM yields the lowest frequencies. The elevated experimental frequencies are primarily due to the additional stiffness introduced by the nylon cord suspension. In contrast, the lower frequencies predicted by the 3D-FEM are attributed to localized face sheet vibrations, which are more pronounced due to their relatively low stiffness at lower modes. The 2D-EPM results lie between the two, capturing a balance between global structural stiffness and localized effects. Error analysis shows that the maximum deviation between the 2D-EPM and the other two methods occurs at the fourth mode, with an error of 10.23%, while all other modes exhibit errors below 10%. These results confirm the accuracy and reliability of the VAM-based equivalent model in predicting the free-boundary modal characteristics of the SP-EHC.

4. Static Verification

Static analysis plays a critical role in evaluating the structural integrity and load-bearing capacity of sandwich panels under various boundary and loading conditions. In this chapter, the static performance of the SP-EHC is examined using both 3D-FEM and the VAM-based 2D-EPM. The investigation covers buckling behavior under in-plane compression, out-of-plane bending with and without initial geometric imperfections, and localized stress and displacement field reconstruction.

4.1. Structural Modeling

The structural parameters used for finite element validation are provided in Table 1. Figure 7 illustrates that the 3D-FEM consists of a metallic core layer sandwiched between upper and lower carbon fiber reinforced polymer (CFRP) face sheets. The unit cell is periodically replicated 21 times along the $x_1$ direction and 16 times along the $x_2$ direction, forming a full-scale SP-EHC. Each face sheet comprises eight CFRP plies with a thickness of 0.25 mm per ply and a symmetric layup sequence of ${[0/90/90/0]}_s$, designed to enhance stiffness and load-bearing capacity in both principal directions. The matrix material is epoxy resin with a fiber volume fraction of 60%. The core layer is fabricated from aluminum, and its material properties are summarized in Table 2. The equivalent stiffness matrix, derived using the VAM, serves as the primary input for constructing the 2D-EPM, as illustrated in Figure 8.

4.2. Buckling Analysis Under In-Plane Loading

Buckling performance analysis of composite sandwich panels is of considerable engineering significance. Owing to their high strength-to-weight ratio, these structures are extensively employed in aerospace, civil engineering, and transportation, offering superior load-bearing capacity and structural stability. Buckling analysis facilitates the evaluation of instability risks under diverse loading and boundary conditions, providing a theoretical basis for structural optimization and safety evaluation.

To evaluate the applicability of the 2D-EPM for buckling prediction, analyses were performed using both 3D-FEM and 2D-EPM under four representative boundary conditions, as shown in Figure 9. In Cases 1 (CSFF) and 2 (CCFF), compressive loads were applied along the short edge, whereas in Cases 3 (SSSS) and 4 (CCSS), loading was applied along the long edge. Symmetric compression along the long edge is commonly employed to investigate the buckling behavior of sandwich panels under longitudinal loading. In practical applications, these panels often experience biaxial compression in the longitudinal direction due to self-weight, wind loads, and other external forces, potentially leading to in-plane buckling and a significant reduction in load-carrying capacity. Simulations under this loading condition enable effective evaluation of critical buckling loads, deformation characteristics, and failure modes, offering valuable insights for structural design.

Conversely, symmetric compression along the short edge examines the panel’s response to transverse loading. Scenarios such as lateral wind pressure or accidental impacts can induce local crushing or buckling near the short edges. This analysis facilitates a detailed assessment of local stiffness, compressive strength, and stability, thereby supporting safety evaluations under extreme loading conditions. A total compressive load of 1 N was applied at the boundaries in all cases; hence, the buckling eigenvalues obtained from the finite element analysis directly correspond to the critical buckling loads.

Table 3 compares the predicted first-order critical buckling loads and corresponding buckling modes of the SP-EHC under four typical boundary conditions, obtained using 2D-EPM and 3D-FEM. The results indicate that the critical buckling loads predicted by the two models show minimal discrepancies, with a maximum deviation of no more than 6%. Moreover, the buckling modes exhibit high morphological agreement, demonstrating that the 2D-EPM offers reliable accuracy and applicability for predicting the first-order buckling behavior of sandwich panels.

Detailed comparisons reveal that the critical buckling load in Case 2 is approximately twice that of Case 1, indicating that enhanced boundary constraints significantly improve structural stability. In addition, the buckling mode shapes reveal that stress concentration regions tend to shift toward the less constrained edges, indicating a higher likelihood of local instability initiating in areas with weaker restraint. A similar trend is observed between Cases 3 and 4, further underscoring the critical role of boundary conditions in influencing buckling behavior.

Notably, in the comparison between Cases 2 and 3, Case 3 exhibits a higher critical buckling load despite having weaker boundary constraints. This seemingly counterintuitive result is primarily attributed to differences in loading direction. When compression is applied along the long edge, the panel benefits from greater bending stiffness and improved buckling resistance. Additionally, the increased structural flexibility in this direction promotes stress redistribution, thereby reducing the likelihood of local buckling. In contrast, although Case 2 features stronger constraints under short-edge compression, it is more prone to stress concentrations and localized deformations, increasing its susceptibility to instability. These findings indicate that critical buckling performance is governed not only by boundary conditions but also by loading direction and structural geometry. Consequently, enhancing boundary support in conjunction with optimizing the loading path offers an effective strategy for improving the load-bearing capacity and buckling stability of sandwich panels.

Considering only the first-order buckling load may not fully capture the structural stability under realistic loading scenarios. Higher-order buckling modes often emerge at increased load levels and can significantly affect the overall instability mechanism, particularly under specific boundary or geometric conditions. Figure 10 shows that the first four buckling mode shapes predicted by the 2D-EPM and 3D-FEM demonstrates strong consistency across various boundary conditions. The relative errors in the corresponding buckling loads remain within 10%, confirming the accuracy and effectiveness of the 2D-EPM in capturing higher-order buckling behavior.

The buckling response of the SP-EHC is highly sensitive to both boundary conditions and loading direction. The number and arrangement of clamped edges significantly influence the critical buckling load and mode shapes: increasing the number of fixed boundaries enhances the overall structural stiffness, elevates the critical load, and tends to localize failure near the constrained regions. In contrast, a higher proportion of free edges leads to more complex buckling deformations and a more dispersed failure pattern. The loading direction also plays a critical role. Under long-edge compression with free boundaries, buckling modes are more evenly distributed across the panel. In comparison, short-edge compression induces instability localized near the free edges, resulting in more intricate and concentrated deformation patterns.

4.3. Out-of-Plane Bending Analysis Without Initial Displacement

Figure 11 illustrates four typical boundary conditions used to evaluate the out-of-plane bending deformation of the SP-EHC under a uniformly distributed load. In the figure, C, S, and F represent clamped, simply supported, and free edges, respectively. Cases 5 to 8 correspond to CCCC, CCCF, CCFF, and SSFF boundary conditions, all subjected to a uniform pressure of 100 kPa with zero initial displacement. A general static analysis step was performed using the ABAQUS software to valuate the maximum displacement and overall deformation behavior of the sandwich panels. Table 4 presents the resulting displacement contours and the relative errors in maximum displacement for each case.

The analysis results reveal that maximum displacement increases progressively as boundary constraints are relaxed—from CCCC to SSFF BCs—indicating that the CCCC BCs provides the greatest out-of-plane bending stiffness, while SSFF BCs offers the least. Specifically, the maximum displacement under SSFF BCs is nearly five times greater than that under CCCC BCs. This is attributed to the lack of moment restraint in the SSFF BCs, where simply supported edges allow free rotation and bending, resulting in significantly larger deformations. Comparison of displacement contours shows strong agreement between the 3D-FEM and 2D-EPM in both deformation patterns and numerical values. The relative error in maximum displacement across all boundary conditions ranges from 2.05% to 4.15%, remaining within acceptable engineering limits. These results confirm that the 2D-EPM provides accurate and reliable predictions of the out-of-plane bending behavior of the SP-EHC under various boundary conditions.

To systematically evaluate the accuracy of the 2D-EPM in predicting the out-of-plane mechanical performance of the SP-EHC, representative displacement extraction paths were designed in addition to the maximum displacement analysis. As illustrated in Figure 12, two characteristic paths were selected for detailed evaluation. Path 1 follows the midline along the longitudinal axis of the sandwich panel, capturing deformation along the primary load-bearing direction. Path 2 extends diagonally across the panel, providing a comprehensive view of global deformation patterns under complex loading conditions.

These two paths enable systematic extraction of displacement–path distribution curves under four representative boundary conditions. This path-based analysis offers several advantages. First, comparing displacement profiles along different paths allows rigorous validation of the 2D-EPM’s ability to capture global deformation patterns, thereby reinforcing its engineering reliability. Second, the curves intuitively reveal spatial variations in deformation gradients, aiding in the detection of localized stress concentrations—critical for structural safety assessment and design refinement. Finally, the method facilitates quantitative comparison of structural responses along symmetric and asymmetric directions, enabling evaluation of the model’s capability to capture anisotropic behavior and enhancing its adaptability across diverse design requirements.

Figure 13 shows the displacement distribution curves along Path 1 and Path 2 under four representative boundary conditions, comparing results from the 3D-FEM and 2D-EPM. Analysis of the curves reveals that for Cases 5 (CCCC), 7 (CCFF), and 8 (SSFF), the displacement profiles along Path 1 exhibit pronounced symmetry in the $x_2$ direction. This behavior is attributed to the centrally symmetric boundary and loading configurations, which promote a uniform deformation response. These results demonstrate that under symmetric boundary conditions, the displacement path curves effectively capture the global deformation characteristics, reinforcing the predictive accuracy of the 2D-EPM for structural behavior prediction.

In contrast, Case 6 (CCCF) exhibits clear asymmetry in the displacement distribution due to the uneven combination of clamped (C) and free (F) boundaries. The displacement curve along Path 2 shows a pronounced gradient variation, with a significantly steeper slope near the clamped edge compared to the free edge. This nonuniform deformation arises because the clamped edge fully restricts rotation, thereby limiting displacement in that region, while the free edge permits greater deformation. As a result, the displacement profile becomes asymmetric, reflecting the influence of boundary condition imbalance on structural response.

Despite the complexity of deformation patterns under asymmetric boundary conditions, the 2D-EPM accurately captures the overall trends of the displacement curves, demonstrating strong predictive capability. These results underscore the model’s robustness and adaptability in addressing complex boundary scenarios. Across all four cases, the displacement–path curves predicted by the 2D-EPM show excellent agreement with those from the 3D-FEM, with relative errors in maximum displacement remaining within acceptable engineering margins. In summary, the 2D-EPM not only reliably predicts the global out-of-plane deformation behavior of the SP-EHC but also effectively accounts for the influence of varying boundary conditions. Its consistent performance—even under complex and asymmetric configurations—establishes a solid theoretical foundation and provides valuable support for structural analysis and engineering design.

4.4. Out-of-Plane Bending Analysis with Initial Displacement

In practical engineering applications, structural components rarely conform to the idealized assumptions of perfect initial geometry. Sandwich panels, in particular, are susceptible to initial displacements or geometric imperfections introduced during manufacturing, transportation, and installation. These imperfections may arise from fabrication errors, material inhomogeneities, environmental factors such as temperature and humidity fluctuations, or external disturbances during handling and assembly. Such geometric deviations, present before the panel enters service, can significantly influence its mechanical response under out-of-plane loading. This effect is especially pronounced in thin-walled or lightweight structures, where even minor imperfections can reduce load-carrying capacity, trigger premature buckling, or induce localized instability.

To more accurately simulate real-world service conditions, this section incorporates initial geometric imperfections into the out-of-plane loading analysis of the SP-EHC under CCFF boundary conditions. These imperfections account for practical factors such as manufacturing deviations and residual stresses, which are commonly encountered in engineering applications. More importantly, the analysis captures the structure’s nonlinear response and potential stability issues under out-of-plane loading, providing valuable insights for the optimization and safety evaluation of sandwich panel designs. A systematic numerical investigation was conducted by introducing prescribed initial displacements of varying amplitudes into the model to evaluate their influence on the transverse load-bearing behavior and ultimate strength of the structure. The simulation procedure is illustrated in Figure 14, where $\Delta x_1$ and $\Delta x_2$ denote the imposed initial displacement amplitudes along the $x_1$ and $x_2$ directions, respectively. The range of initial imperfections considered includes: $-9\%,-6\%,-3\%,0\%,3\%,6\%$, and $9\%$.

Figure 15 illustrates the influence of initial deformation along the $x_1$ direction on the displacement response of the SP-EHC under CCFF boundary conditions. As the initial deformation amplitude increases from -9% to 9%, the total displacement of the SP-EHC exhibits a progressive upward trend. Among the three displacement components, $U_1$ exhibits a nonlinear response: it first decreases and then increases with varying initial deformation. Notably, at zero imperfection (0%), $U_1$ is nearly zero, indicating that displacement in this direction is primarily load-induced, with minimal contribution from geometric imperfections. Furthermore, introducing moderate compressive deformation in the -$x_1$ direction reduces the out-of-plane displacement component $U_3$, effectively enhancing the panel’s out-of-plane stiffness.

Figure 16 illustrates the influence of initial deformation along the $x_2$ direction on the displacement components in all three coordinate directions. As the initial deformation varies from -9% to 9%, the total displacement initially decreases and then increases, though the overall variation remains relatively small. For individual components, $U_1$ reaches a minimum at zero imperfection and increases with greater initial deformation, indicating a clear positive correlation. A similar trend is observed for $U_2$, indicating that in-plane displacements in both directions are sensitive to initial geometric imperfections. In contrast, $U_3$ remains relatively stable throughout the deformation range, implying that imperfections in the $x_2$ direction have minimal impact on the out-of-plane displacement response.

In summary, the SP-EHC under CCFF boundary conditions exhibits notable sensitivity to initial geometric imperfections. Disturbances along the $x_1$ and $x_2$ directions significantly affect the in-plane displacement components ($U_1$ and $U_2$), while their effects on the out-of-plane component ($U_3$) are direction-dependent. Specifically, compressive deformation along the $x_1$ axis contributes to a reduction in $U_3$, whereas deformation along the $x_2$ axis has negligible impact. These findings offer important guidance for controlling initial imperfections and optimizing the stability design of sandwich panels.

4.5. Local Field Recovery

Local field reconstruction is an effective approach for accurately evaluating localized structural performance and informing targeted design optimizations. By enabling precise adjustments to structural configurations or material distributions, this method enhances the strength and stiffness of critical regions, thereby improving overall system reliability and stability. Compared to global response analysis, local field reconstruction offers higher spatial resolution and predictive accuracy, making it particularly well-suited for capturing localized stress, strain, and displacement fields in complex structures such as SP-EHC.

In this section, local field reconstruction is conducted for two representative cases: (1) Case 3 without initial deformation and clamped (CCCC) boundary conditions; and (2) Case 3 with 9% initial deformation along the $x_1$-axis. By comparing the localized stress and displacement fields predicted by the 2D-EPM with those obtained from 3D-FEM, the accuracy and applicability of the dimensional reduction method in reconstructing local fields are systematically evaluated.

Table 5 compares the local stress field distributions within the central cell of the panel in Case 3, as predicted by 3D-FEM and 2D-EPM. The results indicate that the normal bending stresses (${\sigma }_{11}$ and ${\sigma }_{22}$) are primarily concentrated in the bottom face sheet and at the face sheet–core interface, with the maximum tensile stress reaching 71.676 MPa. This stress concentration is likely caused by pronounced bending deformation in the overhanging region, driven by geometric features and external loading. The coupling between face sheet bending and core shear deformation further amplifies this localized stress.

The shear stresses (${\sigma }_{13}$, ${\sigma }_{23}$, and ${\sigma }_{12}$) are primarily distributed along the cell walls of the core in the $x_1$-direction, with a maximum value of 18.372 MPa. Due to the high stiffness and small thickness of the face sheets, they predominantly resist bending, while the core layer bears the majority of the shear load. The observed stress heterogeneity highlights the strong coupling between load transfer and internal structure within the sandwich panel. Notably, shear stresses accumulate near the center of the face–core interface. In combination with elevated shear stresses, this significantly increases the risk of interfacial delamination. Comparison of the local stress contours reveals strong agreement between the 2D-EPM and 3D-FEM predictions, with a maximum relative error of only 1.641%. This confirms the accuracy and effectiveness of the equivalent model in capturing localized stress field behavior.

Table 6 compares the local displacement field distributions obtained from 3D-FEM and 2D-EPM. The results reveal a pronounced displacement concentration at the face–core interface, attributed to the modulus gradient between the face sheets and the core, with the core exhibiting significantly lower stiffness. Under out-of-plane loading, the displacement distribution displays clear biaxial symmetry, consistent with the geometric symmetry of the unit cell and the applied boundary conditions. Notably, the displacement along the $x_1$-direction is concentrated within the face sheet, while the displacement along the $x_2$-direction is concentrated near the face sheet-core interface. This behavior is attributed to the stiffness contrast between the constituent materials and the applied loading conditions. The high-stiffness face sheets dominate the deformation response along the $x_1$-axis, whereas the softer core material experiences greater displacement along the $x_2$-axis—particularly near the face–core interface—due to local stress concentrations induced by the mismatch in material stiffness.

The out-of-plane displacement component $U_3$ shows a red-shifted region in the contour map skewed toward the cell wall of unit cell. This asymmetry results from stress redistribution caused by the free boundary condition in Case 3: the lack of constraint in the $x_1$-direction allows the cell wall to bear a larger portion of the lateral load. Quantitative comparison of the reconstructed displacement fields ($U_1$, $U_2$, $U_3$) shows that the average relative errors of the 2D-EPM with respect to the 3D-FEM reference are 2.7%, 1.7%, and 1.5%, respectively. These results confirm the accuracy and applicability of the proposed dimensional reduction approach in capturing localized bending behavior of sandwich panels.

Table 7 compares the local stress fields within the central cell of the panel under the CCFF boundary condition (Case 3), with an applied initial displacement of 9%. Compared with the results in Table 6, the normal stress ${\sigma }_{11}$ increases by nearly an order of magnitude, while ${\sigma }_{22}$ remains unchanged. This indicates that the imposed initial displacement in the $x_1$-direction significantly amplifies ${\sigma }_{11}$, primarily due to the lack of constraint along this direction. The resulting tensile force induced by the initial imperfection leads to elevated stress concentrations.

Moreover, high tensile stress is observed along the upper outer arm in the $x_1$-direction, especially at the interface with adjacent unit cells, indicating a critical region for stress transmission. The combined effects of out-of-plane loading and initial geometric imperfection leads to stress concentration in this area. Shear stresses ${\sigma }_{13}$ and ${\sigma }_{12}$ also increase significantly—by approximately an order of magnitude—while ${\sigma }_{23}$ remains unchanged. This behavior may result from structural distortion near the boundary induced by the initial displacement, leading to a redistribution of internal forces. Under out-of-plane loading, the amplified local deformation intensifies both stress concentration and shear effects in this region.

Furthermore, the altered geometry introduces discontinuities in local boundary conditions, exacerbating shear stress gradients near the edges. Consequently, structural deformation and stress redistribution propagate inward from the boundaries, further intensifying the shear response. These findings highlight the importance of accounting for initial imperfections in structural core design. For cores susceptible to manufacturing or assembly defects, particular attention should be given to the strength and stiffness at key constraint regions—especially near interface boundaries and bonding lines—to mitigate stress concentrations and deformation amplification induced by initial displacements.

Table 8 compares the local displacement fields within the unit cell located at the same position under CCFF boundary condition (Case 3) with imposed initial displacement of 0.9%. By comparing these results with those in Table 6 and the corresponding displacement contours, several key observations emerge: $U_1$ exhibits a notable increase, $U_2$ shows a slight rise, while $U_3$ exhibits a slight reduction. This suggests that the influence of initial geometric imperfection exceeds that of out-of-plane loading, resulting in a decrease in $U_3$ compared to the perfect, unperturbed structure.

Therefore, the impact of initial geometric imperfections must be carefully considered in the design of sandwich panels. Reinforcing or modifying regions susceptible to such defects is recommended to mitigate their adverse effects. In particular, local stiffening along directions with pronounced imperfections can help reduce displacement concentrations and enhance overall structural robustness.

5. Conclusions

In this study, a novel egg-crate honeycomb core sandwich panel (SP-EHC) was proposed to address the limitations of conventional sandwich panels. A two-dimensional equivalent plate model (2D-EPM) based on the Variational Asymptotic Method (VAM) was developed and systematically validated through experimental and finite element analyses. Based on the findings, the following conclusions can be drawn:

(1) The proposed SP-EHC exhibits superior mechanical performance, including enhanced buckling resistance and improved out-of-plane stiffness. The 2D-EPM accurately predicts both static and dynamic responses of the structure. The relative error in modal frequency prediction is within 10.23%, and the critical buckling load prediction error remains below 5.23% across various boundary conditions.

(2) The equivalent model demonstrates high fidelity in simulating bending deformations under various boundary conditions, both with and without initial geometric imperfections. Maximum displacement errors remain below 4.15%, highlighting the model’s strong applicability and reliability for practical engineering applications.

(3) The 2D-EPM effectively reconstructs local stress and displacement fields, with maximum stress prediction errors below 2.5%, demonstrating its capability in local field analysis. This offers a powerful tool for identifying stress concentrations and guiding structural optimization.

Despite the effectiveness of the VAM-based equivalent model and the novel SP-EHC structure, several limitations persist. Future work will extend the model to incorporate nonlinear material behavior, thermo-mechanical coupling, and dynamic impact response. Moreover, experimental validation under complex loading conditions, including fatigue and multi-axial stress states, is essential to evaluate the SP-EHC’s performance in practical applications.