<figure><img src="image_1.png" /></figure>

HAL
open science

Analytical solutions for the buckling/wrinkling of
anisotropic sandwich structures: Application to
honeycomb cores

Robin Turlin, Philippe Le Grognec

► To cite this version:

Robin Turlin, Philippe Le Grognec. Analytical solutions for the buckling/wrinkling of anisotropic sandwich structures: Application to honeycomb cores. Composite Structures, 2025, 351, pp.118618. 10.1016/j.compstruct.2024.118618. hal-05344954

HAL Id: hal-05344954

https://hal.univ-brest.fr/hal-05344954v1

Submitted on 3 Nov 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Analytical solutions for the buckling/wrinkling of anisotropic sandwich structures: Application to honeycomb cores

Robin Turlinª, Philippe Le Grognecª,\*

ªENSTA Bretagne, UMR CNRS 6027, IRDL, F-29200 Brest, France

## Abstract

Sandwich structures are commonly used in a wide range of industrial sectors, due to their very good compromise between strength and lightness. However, owing to the presence of thin/slender parts, they are prone to buckling under global or local compression, which becomes primarily responsible for the collapse of such sandwiches. This paper thus focuses on the elastic buckling/wrinkling behavior of sandwich columns under overall axial compression. This work follows on from a previous analysis restricted to homogeneous and isotropic skin and core materials. In the present study, the more general case of laminated composite faces is investigated, and anisotropic core materials are also considered. It allows one to deal with heterogeneous core structures in an homogenized way. The modeling approach is based on the use of a laminated beam model for the faces whereas the core layer is represented as a 2D equivalent continuous solid. With such a formulation, analytical solutions can be obtained for the critical loadings and the associated buckling modes, by means of a bifurcation analysis. These solutions are finally confronted to the numerical results derived from 2D linearized buckling finite element computations, for validation purposes, and the particular case of a honeycomb sandwich is addressed.

**Keywords:** Buckling/Wrinkling, Sandwich structures, Anisotropic behavior, Analytical modeling, Finite element validation, Honeycomb cores

<footnote>\*Corresponding author.</footnote>

<footnote>Email address: philippe.le_grognec@ensta-bretagne.fr (Philippe Le Grognec)</footnote>
# 1. Introduction

Sandwich composites are plate-like structures composed of two thin and stiff faces, separated by a thicker and softer core layer. This core layer is often made of a foam-type material, which provides the extreme lightweight property of the sandwich. Sometimes, this core material is strengthened using transverse fibrous reinforcements or replaced by a thin-walled core layer (with a honeycomb or corrugated structure, for instance) so as to reinforce the through-thickness behavior of the whole sandwich. Conversely, the skins and their distance to the middle surface of the composite contribute to the tensile properties and particularly to the flexural rigidity. Owing to this beneficial compromise, sandwich materials are commonly used in applications of aerospace, marine or transportation industries, where the question of specific properties is a key issue. However, the geometric configuration of such composite materials (the thinness of the skins and possibly of the thin-walled core layer, or the slenderness of the core reinforcements) and also the stiffness ratios between the different components contribute to the emergence of instabilities under overall or even localized compressive loads. The buckling phenomenon is therefore one of the major causes for the failure of such composite materials, and it has been greatly investigated in the last decades.

This paper deals with the buckling of sandwich columns under longitudinal compression. Whatever the core configuration, it is assumed to be such that no particular local buckling of reinforcements or thin-walled members arises in the first place, in such a way that the core structure can be replaced by an equivalent homogeneous core material. As a result, only two types of geometric instabilities may occur, namely the global buckling of the sandwich structure and the so-called wrinkling (or local buckling) of the faces, which may appear insofar as they undergo compressive stresses (see Figure 1 for the different mode types observed in sandwich columns under compression). On one hand, the global buckling of a sandwich material can easily be viewed as the classical buckling of an homogeneous structure as soon
as the equivalent properties have been properly derived. On the other hand, the symmetric or antisymmetric wrinkling analysis of sandwich beams requires the use of advanced models, since classical buckling solutions admissible for homogeneous beams are no more valid.

<figure><img src="image_3.png" /><figcaption>Figure 1: Different buckling modes of a compressed sandwich column</figcaption></figure>

The pioneering analytical works dealing with the buckling of sandwich structures date back to the early 1940s. The first contributions on the buckling of sandwich columns under compression (see [1] or [2], for instance) are restricted to a single type of buckling mode (global or local, symmetric or antisymmetric) or deal with several ones, but generally in an uncoupled manner. The critical load predictions are often accurate as far as global buckling is concerned, whereas the numerous expressions derived for the critical wrinkling loads are very scattered. Many analyses are conducted on simplified models made up of a single beam resting on an elastic foundation (which is supposed to stand for the core material), limited thus to symmetric wrinkling, or two beams separated by spring distributions, allowing for the obtention of antisymmetric modes. The foundation is usually described by a one- or two-parameter model (accounting for the transverse extension and possible shear effects) which most often does not convey the actual behavior of the core layer.

As an alternative, many authors have tried to achieve unified higher-order models capable of describing both global and local modes (both symmetric and antisymmetric) while improving the accuracy of the critical load predictions. Benson and Mayers [3] have been the
first to suggest a unified approach to solve the overall buckling and wrinkling problems simul-
taneously, in a linearized framework. More recently, Léotoing et al. [4] also investigated the
overall buckling and local wrinkling of sandwich columns with a unified formulation where
the core material is represented by a higher-order beam model. Jasion et al. [5] performed
an analytical, numerical and experimental study of the global buckling/local wrinkling of
sandwich rectangular beams and circular plates under compression or bending. Simplifying
assumptions in the core kinematics allow them to define simple expressions for the critical
forces or bending moments. In every instance, relative errors of about 10-20% are frequently
observed between the analytical solutions and numerical reference critical values obtained
from finite element computations. Moreover, these errors may change the mode sequence,
and particularly the type (global/local, symmetric/antisymmetric) and the wave number
of the first dominant mode. Let us mention finally Niu and Talreja [6] who have repre-
sented the core material as an elastic continuous medium (instead of spring distributions
or higher-order beam models) characterized by an Airy's function, whose determination has
been proved however to be rather tricky.

As far as anisotropy is concerned, two situations may arise. First, the faces may be
composed of laminated plies with a transverse isotropic (or even orthotropic) constitutive
behavior. This implies defining both traction and bending equivalent longitudinal moduli,
together with the possible coupling terms in the case of asymmetric stacking sequences.

Secondly, the core layer may also behave in an anisotropic way. Either the homogeneous
foam core exhibits a non-isotropic constitutive behavior, or the core layer is made up of a
thin-walled structure (or reinforced) in such a way that the macroscopic equivalent material
is very likely to be anisotropic. The pioneering formulae based on elastic foundation models
should be used in certain particular circumstances, provided that the one or two parameters
related to the elastic foundation are conveniently identified among the constitutive param-
eters of the core material. Conversely, the more advanced solutions based on higher-order
beam or 2D models for the core representation are no more valid and new developments should be performed so as to extend to these new configurations.

At this time, one can already cite a number of authors who have taken an interest in this issue. Dealing specifically with honeycomb structures, one can mention first Norris et al. [7] who conducted many wrinkling experiments on sandwich panels and analyzed the influence of the honeycomb cell size on the critical loads. Later on, Pearce and Webber [8] and Webber et al. [9] also performed some experimental buckling experiments on honeycomb panels under axial compression. As stated by many authors, there is not a lot of experimental data in the literature concerning the buckling/wrinkling of sandwich beams, especially with honeycomb cores, since it is very difficult to achieve such buckling tests and produce repeatable results in terms of critical loadings. The reason for these difficulties arises from the complex geometric and material nature of a sandwich structure, displaying classically composite laminated faces and an heterogeneous (foam or thin-walled) core layer. As highlighted in [10], the core material response is often non-linear long before local buckling occurs in the sandwich beam. Moreover, the standard experimental set-ups often induce eccentric compressive loads on the faces which give rise to complex stress/strain fields in the core layer and cause premature failure of the structure. More generally, sandwich structures are very imperfection-sensitive, due to the possible interactive buckling phenomena between global and local (symmetric and antisymmetric) modes. For all these reasons, despite maximum precautions, the rare experimental results available in the literature are thus very scattered, as can be seen in [8, 9, 11, 12].

From a theoretical point of view, as long as the honeycomb thin-walled structure is not supposed to buckle locally, the case of honeycomb cores can be viewed as a particular case of orthotropic core material. Many studies deal thus with general anisotropic (most often orthotropic) core materials which are likely to reflect both cases of isotropic or transversely isotropic foam and honeycomb cores. All the aspects discussed above (non-linearities, im-
perfections) are still very difficult to consider using analytical approaches, so that they are
generally disregarded, as will be the case in this paper. It is therefore hazardous to rely
on experimental results in order to validate new analytical solutions, so that in most recent
studies, the new expressions derived for the critical load of a sandwich beam under longitu-
dinal compression are essentially compared to reference numerical results. One can mention
Kim and Hong [13], who investigated the general case of unbalanced anisotropic sandwich
plates, using the Rayleigh-Ritz approximate method. Following Benson-Mayers theory, Hadi
and Matthews [14] implemented a general approach for the buckling and wrinkling analysis
of anisotropic sandwich panels, where the core is assumed to be anti-plane (the longitudinal
stiffness is neglected). Their numerical results were confronted to previously mentioned ex-
periments [8, 9]. Vonach and Rammerstorfer [15] proposed a general analytical procedure so
as to deal with the wrinkling behavior of orthotropic sandwich plates under various loading
conditions, but considering only infinitely thick cores. Fagerberg and Zenkert [16], as for
them, used a strip model to analyze the wrinkling of highly anisotropic sandwich panels
under multiaxial loading. More recently, D'Ottavio and Polit [17] analyzed the global and
local buckling problem of sandwich struts, by means of a general approximate 2D model.
They particularly emphasized the influence of the orthotropy ratio of the core material on
the occurrence of symmetric or antisymmetric wrinkling modes. A few years later, Vescovini
et al. [18] carried out a comprehensive study of the buckling/wrinkling behavior of sand-
wich plates with both laminated faces and an orthotropic core material, using Ritz method.
One can also mention Fagerberg [19], who discussed the relative influence of the traction
and bending stiffnesses of the faces on the buckling design of sandwich panels, thanks to
appropriate experimental analyses, and suggested the use of minisandwich faces in order to
increase the wrinkling loads in a cost-effective manner.

Besides, Wadee and Hunt [20] focused on the unstable post-buckling response of sand-
wich structures, after localization of the deformed shape due to the interactive buckling
phenomenon, in the general case of orthotropic core materials. Stiftinger and Rammerstofer [21] also investigated the non-linear post-buckling response of anisotropic sandwich shells, from a numerical and experimental point of view, introducing material non-linearities either in the faces or in the crushable foam core. The interested reader may refer to the review paper of Ley et al. [12] for an overview of theoretical and experimental wrinkling analyses of sandwich structures with isotropic or laminated facesheets and solid or cellular cores.

This paper follows on from a previous article from one of the authors dealing with the buckling analysis of isotropic sandwich beam-columns [22]. In this earlier work, original closed-form solutions have been obtained for the global/local antisymmetric/symmetric critical loads. A very good accuracy was achieved, by comparison to other simplified models based on elastic foundations or higher-order beam models, thanks to the full 2D description of the core layer using no specific hypotheses for the through-thickness kinematics. In the present paper, the idea is to investigate the case of anisotropic sandwich columns under axial compression, using the same kinematics as before. While the faces are classically considered as slender beams, the core is represented by a 2D continuous solid, without any simplification regarding the deformation field. The main improvement comes from the fact that one now considers the more general case of laminated (possibly asymmetric) composite faces (made up of isotropic or transverse isotropic plies) together with an orthotropic equivalent core material. A classical elastic bifurcation analysis is then carried out, based on the so-called bifurcation equation of the problem. This equation is particularized, using the above kinematic assumptions and considering a longitudinal pre-critical compressive stress state. This general method was already applied in [22] but also in many other geometric, material and loading configurations, including the case of composite beams with partial interaction [23], reinforced sandwich materials under through-thickness compression [24] and sandwich panels under biaxial loading [25]. The final set of partial differential equations is analytically solved for given classical boundary conditions. Original solutions are obtained in the general case.
In particular conditions, the previous solutions can be simplified by cancelling negligible moduli, leading to very simple expressions which sometimes coincide with earlier formulae from the literature. The new solutions are validated against numerical finite element results (from linearized buckling analyses) in different configurations, including the particular case of honeycomb cores. In most cases, 2D finite element models are developed in line with the analytical model, but without any simplification, in order to obtain reference results that are as general as possible. In the case of honeycomb cores, a 3D representation of the sandwich beam (using shell elements) is particularly retained. Parametric analyses are finally achieved to show the relative influence of various material parameters on the predictions of critical loadings and modes.

## 2. Analytical study

### 2.1. Problem definition

The present study deals with sandwich columns (namely narrow sandwich panels) whose response is assumed to be invariant in the width direction (along the y-axis), in such a way that a 2D representation of the sandwich structure in the xz-plane is sufficient, as displayed in Figure 2. The length and width of the sandwich structure are denoted by $L$ and $b$, whereas the total thicknesses of the core layer and the two skins are respectively $2h_c$ and $2h_s$. The two skins are supposed to be identical (in which case the sandwich material is symmetric) like in most practical cases.

The core material is considered to be an homogeneous 2D continuous medium satisfying the plane stress hypothesis (due to the small lateral dimensions). It is assumed to be linearly elastic and orthotropic in the most general case. It allows one to deal with many homogeneous cores like foams, but also with thin-walled core structures like honeycombs, in an homogenized way. In the present 2D context, only four material parameters will
<figure><img src="image_4.png" /><figcaption>Figure 2: Two-dimensional representation of the sandwich column under longitudinal compression</figcaption></figure>

be required for the following calculations, namely the longitudinal and transverse Young's moduli $E_x$ and $E_z$, the Poisson's ratio $\nu_{xz}$ and the “in-plane” shear modulus $G_{xz}$, which actually corresponds to an out-of-plane shear modulus in relation to the plane of the sandwich faces.

As for the two skins, they consist in laminated composites whose stacked layers are all orthotropic. As they are further represented by Euler-Bernoulli beams in the adopted modeling approach (transverse shear effects are neglected due to their low thickness-to-length ratios), only the longitudinal stiffnesses (along the x-axis) will be useful in the sequel. These are the extensional stiffness $A$, the bending stiffness $D$ and possibly the coupling stiffness $B$ (for asymmetric stacking sequences), which can all be obtained by the appropriate through-thickness integration of the longitudinal reduced moduli of the successive layers. The classical laminate theory is briefly recalled in Appendix A, from which it is possible to deduce easily the desired stiffnesses.

Lastly, the left and right end sections of the sandwich column are respectively subjected to zero and non-zero (compressive) displacement boundary conditions in the longitudinal direction, what leads to buckling (this particular choice of displacement boundary conditions has been retained as it gives rise to piecewise uniform stress/strain states in the sandwich column during the pre-critical stage).

In the sequel, a total Lagrangian formulation of the buckling problem is first implemented
(from a general 3D framework) in a weak form, using the so-called bifurcation equation [26]. Introducing the appropriate kinematics and stress states results in the achievement of a set of partial differential governing equations as well as associated natural boundary conditions. Finally, by solving this eigenvalue problem, one can derive the sought critical displacements and corresponding bifurcation modes. Such a procedure has already been successfully used in many other contexts, as in the case of isotropic sandwich beam-columns and plates or composite beams with partial interaction, for example [22, 23, 25].

## 2.2. Governing equations

The critical loading $\lambda_{cr}$ and the bifurcation mode $\mathbf{X}$ of a 3D body can generally be obtained by solving the following bifurcation equation:

$$ \forall \delta \mathbf{u}, \quad \int_{\Omega} \nabla^T \delta \mathbf{u} : \mathbf{K}(\lambda_{cr}) : \nabla \mathbf{X} d\Omega = 0 \quad (1) $$

where $\delta \mathbf{u}$ can be considered as a kinematically admissible test function or, more physically, as the virtual variation of the unknown displacement field $\mathbf{u}$.

In Equation (1), $\mathbf{K}$ represents the fourth-order nominal tangent elastic tensor, which can be written as follows:

$$ \mathbf{K} = \mathbf{F} \cdot \mathbf{D} \cdot \mathbf{F}^T + (\mathbb{I} \cdot \mathbf{\Sigma})^T \quad (2) $$

In the above equation, $\mathbb{I}$ stands for the fourth-order unit tensor and the superscript $^T$ the transposition of a second-order tensor and the major transposition of a fourth-order tensor, depending on the case. $\mathbf{F}$ denotes the deformation gradient and $\mathbf{D}$ represents the fourth-order material tangent elastic tensor which may depend on the current point if ever the structure is heterogeneous.

When dealing with the skins, a 1D model is intended, so that the transverse normal material stresses $\Sigma_{yy}$ and $\Sigma_{zz}$ must be taken equal to 0 so as to enforce the appropriate
stress state in the beams. This amounts to replacing tensor **D**<sup>s</sup> with the associated reduced tensor **C**<sup>s</sup> defined by the following relations:

$$C_{ijkl}^{s} = D_{ijkl}^{s} + \frac{D_{ijyy}^{s}(D_{yzz}^{s}D_{zzkl}^{s}-D_{zzzz}^{s}D_{yykl}^{s})+D_{ijzz}^{s}(D_{zzyy}^{s}D_{yykl}^{s}-D_{yyyy}^{s}D_{zzkl}^{s})}{D_{yyyy}^{s}D_{zzzz}^{s}-D_{yzzz}^{s}D_{zyyy}^{s}} \quad (3)$$
$$(i,j) \neq (y,y), (z,z) \quad (k,l) \neq (y,y), (z,z)$$

As regards the core layer, the plane stress assumption leads one to replace tensor **D**<sup>c</sup> by the appropriate reduced tensor **C**<sup>c</sup>, by enforcing the zero out-of-plane stress condition $\Sigma_{yy} = 0$. The following plane components will be used subsequently:

$$C_{xxxx}^{c} = \frac{E_x^2}{E_x - E_z \nu_{xz}^2} \qquad C_{zzzz}^{c} = \frac{E_x E_z}{E_x - E_z \nu_{xz}^2} \qquad (4)$$
$$C_{xxzz}^{c} = \frac{E_x E_z \nu_{xz}}{E_x - E_z \nu_{xz}^2} \qquad C_{xzxz}^{c} = G_{xz}$$

the others being zero or deduced by major or minor symmetries.

In the present study, the (symmetric) second Kirchhoff stress tensor $\Sigma$ involved in Equation (2) reflects the uniaxial pre-critical stress state due to the uniform compressive displacement applied on the right-hand side of the sandwich column, in both the facings and the core layer. However, owing to the respective moduli of the two components, these stresses are far larger in the skins than in the core material, so that they will be neglected in the latter.

The skins alone are therefore subject to a compressive force. In the general case of laminated faces, each ply is submitted to a nominal axial compressive stress $\Pi_{xx} = -P < 0$ in its longitudinal direction, so that the first Kirchhoff stress tensor **Π** can be expressed in the orthonormal basis ($\mathbf{e}_x, \mathbf{e}_y, \mathbf{e}_z$) as:

$$\mathbf{\Pi} = -P \mathbf{e}_x \otimes \mathbf{e}_x = \begin{bmatrix} -P & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} \quad (P > 0) \qquad (5)$$
Assuming small pre-critical deformations (||∇u|| ≪ 1), the corresponding stress tensor Σ writes then:

$$ \Sigma = \mathbf{F}^{-1} \cdot \Pi \approx \Pi \qquad (6) $$

and the nominal tangent elastic tensor in the ply becomes:

$$ \mathbf{K}^s \approx \mathbf{C}^s - P\mathbf{e}_i \otimes \mathbf{e}_x \otimes \mathbf{e}_x \otimes \mathbf{e}_i \qquad (7) $$

using the implicit summation convention on repeated indices with i = x, y, z.

Similarly, the nominal tangent elastic tensor in the core layer simply writes:

$$ \mathbf{K}^c \approx \mathbf{C}^c \qquad (8) $$

since initial stresses have been excluded in this case.

Eventually, the bifurcation equation (1) of the whole sandwich structure writes:

$$ \begin{aligned} & \forall \delta \mathbf{u}_a, \delta \mathbf{u}_b, \delta \mathbf{u}_c, \quad \int_{\Omega_a} \nabla^T \delta \mathbf{u}_a : (\mathbf{C}^s - P_{cr} \mathbf{e}_i \otimes \mathbf{e}_x \otimes \mathbf{e}_x \otimes \mathbf{e}_i) : \nabla \mathbf{X}_a d\Omega_a \\ & + \int_{\Omega_b} \nabla^T \delta \mathbf{u}_b : (\mathbf{C}^s - P_{cr} \mathbf{e}_i \otimes \mathbf{e}_x \otimes \mathbf{e}_x \otimes \mathbf{e}_i) : \nabla \mathbf{X}_b d\Omega_b + \int_{\Omega_c} \nabla^T \delta \mathbf{u}_c : \mathbf{C}^c : \nabla \mathbf{X}_c d\Omega_c = 0 \end{aligned} \qquad (9) $$

In Equation (9), $\mathbf{u}_a, \mathbf{u}_b, \mathbf{u}_c, \mathbf{X}_a, \mathbf{X}_b$ and $\mathbf{X}_c$ represent the displacement field and bifurcation mode components, respectively, relative to the two skins (with indices $\bullet_a$ and $\bullet_b$) and the core (with index $\bullet_c$). Besides, $P_{cr}$ corresponds to the critical value of the compressive stress in the faces, which may depend on the considered ply.

Let us now consider the bending problem of the skins in the $xz$-plane, using the Euler-Bernoulli beam theory (since transverse shear effects are negligible in practice). The kinematics of an Euler-Bernoulli beam is defined by two scalar displacement fields $u(x)$ and $w(x)$, respectively the axial and transverse displacements of the centroid axis. The following expressions are thus retained for the bifurcation modes $\mathbf{X}_i$ as well as the displacement
variations $\delta \mathbf{u}_i$ (with $i = a, b$ depending on the skin considered):

$$ \mathbf{X}_i = \begin{vmatrix} U_i - zW_{i,x} & & \delta u_i - z\delta w_{i,x} \\ 0 & & 0 \\ W_i & & \delta w_i \end{vmatrix} \qquad (10) $$

Besides, the bifurcation mode $\mathbf{X}_c$ and the displacement variation $\delta \mathbf{u}_c$ in the core layer are classically expressed in the orthonormal basis ($\mathbf{e}_x, \mathbf{e}_y, \mathbf{e}_z$) as follows:

$$ \mathbf{X}_c = \begin{vmatrix} U_c & & \delta u_c \\ 0 & & 0 \\ W_c & & \delta w_c \end{vmatrix} \qquad (11) $$

The global bifurcation equation then writes:

$$ \begin{aligned} & \forall \delta u_a, \delta w_a, \delta u_b, \delta w_b, \delta u_c, \delta w_c, \quad \int_{\Omega_a} [C_{xxxx}^s(U_{a,x} - zW_{a,xx})(\delta u_{a,x} - z\delta w_{a,xx}) \\ & \qquad -P_{cr}(U_{a,x} - zW_{a,xx})(\delta u_{a,x} - z\delta w_{a,xx}) - P_{cr}W_{a,x}\delta w_{a,x}] d\Omega_a \\ & \qquad + \int_{\Omega_b} [C_{xxxx}^s(U_{b,x} - zW_{b,xx})(\delta u_{b,x} - z\delta w_{b,xx}) \\ & \qquad -P_{cr}(U_{b,x} - zW_{b,xx})(\delta u_{b,x} - z\delta w_{b,xx}) - P_{cr}W_{b,x}\delta w_{b,x}] d\Omega_b \\ & \qquad + \int_{\Omega_c} (C_{xxxx}^c U_{c,x}\delta u_{c,x} + C_{xxzx}^c W_{c,z}\delta u_{c,z} + C_{zzzx}^c U_{c,x}\delta w_{c,z} + C_{zzzz}^c W_{c,z}\delta w_{c,z} \\ & \qquad + C_{xxzx}^c W_{c,x}\delta u_{c,z} + C_{xxzx}^c U_{c,z}\delta u_{c,z} + C_{zzzx}^c W_{c,x}\delta w_{c,z} + C_{zxzx}^c U_{c,z}\delta w_{c,z}) d\Omega_c = 0 \end{aligned} \qquad (12) $$

where $z$ stands for the $z$-coordinate of a current point relative to the centroid axis of the considered skin.

Some quantities above must be first integrated through the thickness of the faces. The only component of the reduced elastic tensor involved in the present calculations is easy to
integrate across all the plies of the upper laminated face:

$$ 
\begin{aligned}
\int_{-h_s}^{h_s} C_{xxxx}^s dz &= A \\
\int_{-h_s}^{h_s} C_{xxxx}^s z dz &= B \\
\int_{-h_s}^{h_s} C_{xxxx}^s z^2 dz &= D
\end{aligned} 
\qquad (13) $$

where A, B and D are the reduced extensional, coupling and bending stiffnesses of the upper face, which can be obtained through the classical laminate theory (see Equation (A.9) in Appendix A). Owing to the sandwich global symmetry, the same integrals carried out on the lower face give rise to the same values, except B which is replaced by $-B$. Next, the local compressive stress ($P > 0$) can be related to the enforced displacement $\lambda > 0$ (which will act further as the bifurcation parameter) as follows:

$$ P = \frac{C_{xxxx}^s \lambda}{L} \qquad (14) $$

and integrating the critical stresses throughout the thickness of the faces leads thus to the following expression:

$$ \int_{-h_s}^{h_s} P_{cr} dz = \frac{A\lambda_{cr}}{L} \qquad (15) $$

Furthermore, the integration over the width of the sandwich column is straightforward, since all the fields are invariant along this direction. One can thus choose arbitrarily a unit width ($b = 1$) for the subsequent calculations, regardless of the actual width, as long as we are talking about critical displacements.

Afterwards, integrating by parts with respect to $x$ and $z$, and eliminating negligible higher-order terms (presupposing that $\lambda_{cr} \ll L$) yields six local partial differential equations
for the components $U_a, W_a, U_b, W_b, U_c$ and $W_c$ of the eigenmode:

$$
\begin{align}
AU_a,xx + BW_a,xxx + C_{xxxz}^c (U_c,z + W_c,x)|_{z=-h_c} &= 0 \tag{a} \\
DW_a,xxxx + BU_a,xxx + \frac{A\lambda_{cr}}{L} W_a,xx - C_{xxzz}^c U_c,x|_{z=-h_c} \nonumber \\
&\quad - C_{zzzz}^c W_c,z|_{z=-h_c} - h_s C_{xxxz}^c (U_c,xz + W_c,xx)|_{z=-h_c} = 0 \tag{b} \\
AU_b,xx - BW_b,xxx - C_{xxxz}^c (U_c,z + W_c,x)|_{z=h_c} &= 0 \tag{c} \\
DW_b,xxxx - BU_b,xxx + \frac{A\lambda_{cr}}{L} W_b,xx + C_{xxzz}^c U_c,x|_{z=h_c} \nonumber \\
&\quad + C_{zzzz}^c W_c,z|_{z=h_c} - h_s C_{xxxz}^c (U_c,xz + W_c,xx)|_{z=h_c} = 0 \tag{d} \\
C_{xxxx}^c U_c,xx + C_{xxxz}^c U_c,zz + (C_{xxzz}^c + C_{xxxz}^c) W_c,xz &= 0 \tag{e} \\
C_{zzzz}^c W_c,zz + C_{xxxz}^c W_c,xx + (C_{xxzz}^c + C_{xxxz}^c) U_c,xz &= 0 \tag{f}
\end{align}
$$

together with the natural boundary conditions, which will be further detailed.

The last two equations (16(e)) and (16(f)) identify with the classical local equilibrium equations of the core region in a 2D framework. Similarly, the first four equations (16(a)), (16(b)), (16(c)) and (16(d)) look like buckling differential equations of laminated beams, but include new quantities characterizing the influence of stresses at the interface between the skin considered and the core structure during the buckling phenomenon. Let us mention that these additional terms are here naturally obtained through integrations by parts performed in the core layer.

Let us now specify the boundary conditions in order to solve the previous differential system. First, connecting conditions between the different components of the bifurcation mode must be satisfied at the interfaces between the core and the two facings, namely:

$$
\begin{align}
& U_a - h_s W_{a,x} - U_c|_{z=-h_c} = 0 \nonumber \\
& U_b + h_s W_{b,x} - U_c|_{z=h_c} = 0 \tag{17} \\
& W_a - W_c|_{z=-h_c} = 0 \nonumber \\
& W_b - W_c|_{z=h_c} = 0 \nonumber
\end{align}
$$
Owing to the enforced displacement boundary conditions, the two ends of each face act as if they were guided, what leads to the following kinematical constraints: $U_a(0) = U_a(L) = U_b(0) = U_b(L) = 0$ and $W_{a,x}(0) = W_{a,x}(L) = W_{b,x}(0) = W_{b,x}(L) = 0$. Taking into account $\delta u_a(0) = \delta u_a(L) = \delta u_b(0) = \delta u_b(L) = 0$ and $\delta w_{a,x}(0) = \delta w_{a,x}(L) = \delta w_{b,x}(0) = \delta w_{b,x}(L) = 0$ in the bifurcation equation (12) leads one, after integration by parts, to the remaining stress boundary conditions at the ends $x = 0$ and $x = L$:

$$ 
\begin{aligned}
& DW_{a,xxx} + BU_{a,xx} + \frac{A\lambda_{cr}}{L} W_{a,x} - h_s C_{xzxz}^c (U_{c,z} + W_{c,x}) \Big|_{z=-h_c} = 0 \\
& DW_{b,xxx} - BU_{b,xx} + \frac{A\lambda_{cr}}{L} W_{b,x} - h_s C_{xzxz}^c (U_{c,z} + W_{c,x}) \Big|_{z=h_c} = 0
\end{aligned}
\qquad (18) 
$$

The last boundary conditions refer to the two remaining edges of the core zone. The enforced displacement boundary conditions induce kinematical constraints for the longitudinal component of the bifurcation mode at $x = 0$ and $x = L$, namely $U_c(= U_a = U_b) = 0$, $\forall z \in ] -h_c, h_c[$. Since these two edges are free in the $z$-direction, they satisfy the following stress boundary condition:

$$ C_{xzxz}^c (U_{c,z} + W_{c,x}) = 0 \qquad (19) $$

### 2.3. General solution procedure

According to the classical sinusoidal shapes of the buckling modes of a single beam under axial compression with the boundary conditions defined above, and in line with numerical observations, the following expressions are introduced for the beam components of the bifurcation modes of the sandwich column, assuming then that the core layer does not alter the transverse modal displacement shape:

$$ 
\left\{
\begin{aligned}
U_a &= \alpha \sin \frac{n\pi x}{L} \\
U_b &= \beta \sin \frac{n\pi x}{L} \\
W_a &= \pm \cos \frac{n\pi x}{L} \\
W_b &= \cos \frac{n\pi x}{L}
\end{aligned}
\right.
\qquad (20) 
$$
where $n$ is the associated half-wave number.

Owing to the symmetry of the sandwich column, the same unit amplitude is retained for the buckling mode components $W_a$ and $W_b$ of both faces. However, two cases are considered, depending on the relative sign of the two fields $W_a$ and $W_b$. The bifurcation mode of the sandwich column may be antisymmetric (with $W_a = W_b$) or symmetric (with $W_a = -W_b$). Sinusoidal shapes are also retained for the longitudinal components $U_a$ and $U_b$, consistently with the associated boundary conditions. Unknown amplitudes $\alpha$ and $\beta$ are introduced here, which will turn out to be equal or opposite to each other in the symmetric and antisymmetric cases, respectively.

Concerning the core modal displacement field, a separation of variables is performed and the following forms are presupposed, according to Equation (20):

$$ \left\{ \begin{array}{l} U_c = \zeta(z) \sin \frac{n\pi x}{L} \\ W_c = \xi(z) \cos \frac{n\pi x}{L} \end{array} \right. \qquad (21) $$

Solving first Equations (16(e)) and (16(f)), one obtains the following general expressions for functions $\zeta(z)$ and $\xi(z)$:

$$ \left\{ \begin{array}{l} \zeta(z) = \zeta_1 e^{-\frac{n\pi\eta-z}{L}} + \zeta_2 e^{\frac{n\pi\eta-z}{L}} + \zeta_3 e^{-\frac{n\pi\eta+z}{L}} + \zeta_4 e^{\frac{n\pi\eta+z}{L}} \\ \xi(z) = \xi_1 e^{-\frac{n\pi\eta-z}{L}} + \xi_2 e^{\frac{n\pi\eta-z}{L}} + \xi_3 e^{-\frac{n\pi\eta+z}{L}} + \xi_4 e^{\frac{n\pi\eta+z}{L}} \end{array} \right. \qquad (22) $$

with:

$$ \eta_{\pm} = \sqrt{\frac{E_z (E_x - 2\nu_{xz} G_{xz}) \pm \sqrt{E_z (E_z (E_x - 2\nu_{xz} G_{xz})^2 - 4E_x G_{xz}^2)}}{2E_z G_{xz}}} \qquad (23) $$

In Equation (22), $\zeta_1$, $\zeta_2$, $\zeta_3$, $\zeta_4$, $\xi_1$, $\xi_2$, $\xi_3$ and $\xi_4$ depend only on four coefficients which can be determined by satisfying the four connecting conditions (17). The remaining boundary conditions are automatically verified. Then, Equations (16(a)) and (16(c)) are solved in order to determine $\alpha$ and $\beta$. Lastly, Equations (16(b)) and (16(d)) give rise to the same
expression of the critical displacement $\lambda_{cr}$. This solution depends on whether the modes are symmetric or antisymmetric. Explicit formulae are obtained, but they are too cumbersome in the general case to be presented as closed-form expressions.

## 2.4. Particular case of an isotropic core material

In many cases, the core layer can be described as an isotropic equivalent homogeneous material, whereas the faces are still laminated composites. If so, one can use the same equations as before, setting $E_x = E_z = E_c$, $\nu_{xz} = \nu_c$ and $G_{xz} = \frac{E_c}{2(1+\nu_c)}$. In this common case, simple closed-form solutions can be derived for the critical displacements. In the sequel, the two antisymmetric and symmetric cases are treated separately.

### 2.4.1. Antisymmetric case

The solutions of Equations (16(e)), (16(f)), (16(a)) and (16(c)), combined with the four connecting conditions (17), when assuming $W_a = W_b = \cos \frac{n\pi x}{L}$, take the following degenerated form:

$$ \left\{ \begin{aligned} \zeta(z) &= \frac{[-k_1 n\pi(1+\nu_c)+k_2 L(\nu_c-3)]\sinh \frac{n\pi z}{L} - k_2 n\pi(1+\nu_c)z \cosh \frac{n\pi z}{L}}{n\pi(1+\nu_c)} \\ \xi(z) &= k_1 \cosh \frac{n\pi z}{L} + k_2 z \sinh \frac{n\pi z}{L} \end{aligned} \right. \qquad (24) $$

with:

$$ \left\{ \begin{aligned} k_1 &= \frac{-n\pi L h_c(1+\nu_c) \cosh \frac{n\pi h_c}{L} + [L^2(\nu_c-3)+n^2\pi^2 h_s h_c(1+\nu_c)+\alpha n\pi L h_c(1+\nu_c)] \sinh \frac{n\pi h_c}{L}}{L[L(\nu_c-3) \cosh \frac{n\pi h_c}{L} \sinh \frac{n\pi h_c}{L} - n\pi h_c(1+\nu_c)]} \\ k_2 &= \frac{n\pi(1+\nu_c) [-n\pi h_s + \alpha L] \cosh \frac{n\pi h_c}{L} + L \sinh \frac{n\pi h_c}{L}}{L[L(\nu_c-3) \cosh \frac{n\pi h_c}{L} \sinh \frac{n\pi h_c}{L} - n\pi h_c(1+\nu_c)]} \end{aligned} \right. \qquad (25) $$

given that:

$$ \left\{ \begin{aligned} \alpha &= \frac{2E_c L^2 n \pi h_s \cosh^2 \frac{n \pi h_c}{L} + [B n^2 \pi^2 (1 + \nu_c)(\nu_c - 3) - E_c L^2 (\nu_c - 1)] L \cosh \frac{n \pi h_c}{L} \sinh \frac{n \pi h_c}{L} - n \pi h_c (1 + \nu_c) [B n^2 \pi^2 (1 + \nu_c) - E_c L^2]}{L [-2 E_c L^2 \cosh^2 \frac{n \pi h_c}{L} + A L n \pi (1 + \nu_c)(\nu_c - 3) \cosh \frac{n \pi h_c}{L} \sinh \frac{n \pi h_c}{L} - A n^2 \pi^2 h_c (1 + \nu_c)^2]} \\ \beta &= -\alpha \end{aligned} \right. \qquad (26) $$
Then, Equations (16(b)) and (16(d)) lead to the same critical displacement:

$$
\lambda_{cr}^{A} = \frac{(-2E_c L^2 n\pi [n^2\pi^2(Ah_s^2 + 2Bh_s + D) + AL^2] \cosh^2 \frac{n\pi h_c}{L} + L [n^4\pi^4(1+\nu_c)(\nu_c-3)(AD-B^2) + 2E_c L^2 n^2\pi^2(\nu_c-1)(Ah_s + B) - E_c^2 L^4] \cosh \frac{n\pi h_c}{L} \sinh \frac{n\pi h_c}{L} - n\pi [n^4\pi^4 h_c(1+\nu_c)^2(AD-B^2) + 2E_c L^2 n^2\pi^2 h_c(1+\nu_c)(Ah_s + B) - E_c L^4(2A + E_c h_c)]}{(ALn\pi [-2E_c L^2 \cosh^2 \frac{n\pi h_c}{L} + ALn\pi(1+\nu_c)(\nu_c-3) \cosh \frac{n\pi h_c}{L} \sinh \frac{n\pi h_c}{L} - An^2\pi^2 h_c(1+\nu_c)^2])} \quad (27)
$$

2.4.2. Symmetric case

Assuming now $W_a = -\cos \frac{n\pi X}{L} = -W_b$ leads to somewhat different expressions relative to the so-called symmetric buckling solution. The same procedure is employed and gives rise to the following expressions:

$$
\left\{
\begin{aligned}
\alpha &= \frac{-2E_c L^2 n \pi h_s \cosh^2 \frac{n \pi h_c}{L} - [B n^2 \pi^2 (1+\nu_c)(\nu_c-3) - E_c L^2 (\nu_c-1)] L \cosh \frac{n \pi h_c}{L} \sinh \frac{n \pi h_c}{L} - n \pi [B n^2 \pi^2 h_c (1+\nu_c)^2 - E_c L^2 (h_c (1+\nu_c) + 2h_s)]}{L [-2E_c L^2 \cosh^2 \frac{n \pi h_c}{L} + A L n \pi (1+\nu_c)(\nu_c-3) \cosh \frac{n \pi h_c}{L} \sinh \frac{n \pi h_c}{L} + A n^2 \pi^2 h_c (1+\nu_c)^2 + 2 E_c L^2]} \\
\beta &= \alpha \\
k_1 &= \frac{[L^2 (\nu_c-3) + n^2 \pi^2 h_s h_c (1+\nu_c) - \alpha n \pi L h_c (1+\nu_c)] \cosh \frac{n \pi h_c}{L} - n \pi L h_c (1+\nu_c) \sinh \frac{n \pi h_c}{L}}{L [L (\nu_c-3) \cosh \frac{n \pi h_c}{L} \sinh \frac{n \pi h_c}{L} + n \pi h_c (1+\nu_c)]} \\
k_2 &= \frac{n \pi (1+\nu_c) [L \cosh \frac{n \pi h_c}{L} - (n \pi h_s - \alpha L) \sinh \frac{n \pi h_c}{L}]}{L [L (\nu_c-3) \cosh \frac{n \pi h_c}{L} \sinh \frac{n \pi h_c}{L} + n \pi h_c (1+\nu_c)]} \\
\zeta(z) &= \frac{[-k_1 n \pi (1+\nu_c) + k_2 L (\nu_c-3)] \cosh \frac{n \pi Y}{L} - k_2 n \pi (1+\nu_c) z \sinh \frac{n \pi z}{L}}{n \pi (1+\nu_c)} \\
\xi(z) &= k_1 \sinh \frac{n \pi z}{L} + k_2 z \cosh \frac{n \pi z}{L}
\end{aligned}
\right.
\qquad (28)
$$

In this case, the critical displacement writes:

$$
\lambda_{cr}^{S} = \frac{(-2E_c L^2 n\pi [n^2\pi^2(Ah_s^2 + 2Bh_s + D) + AL^2] \cosh^2 \frac{n\pi h_c}{L} + L [n^4\pi^4(1+\nu_c)(\nu_c-3)(AD-B^2) + 2E_c L^2 n^2\pi^2(\nu_c-1)(Ah_s + B) - E_c^2 L^4] \cosh \frac{n\pi h_c}{L} \sinh \frac{n\pi h_c}{L} + n\pi [n^4\pi^4 h_c(1+\nu_c)^2(AD-B^2) + 2E_c L^2 n^2\pi^2(Ah_s^2 + [A(1+\nu_c)h_c + 2B]h_s + B(1+\nu_c)h_c + D) - E_c^2 L^4 h_c])}{(ALn\pi [-2E_c L^2 \cosh^2 \frac{n\pi h_c}{L} + ALn\pi(1+\nu_c)(\nu_c-3) \cosh \frac{n\pi h_c}{L} \sinh \frac{n\pi h_c}{L} + An^2\pi^2 h_c(1+\nu_c)^2 + 2E_c L^2])} \quad (29)
$$
## 2.5. Particular case of a thin-walled core structure

In the case of thin-walled core structures such as honeycomb cores, the longitudinal modulus $E_x$ is often negligible when compared to the transverse modulus $E_z$. By setting to zero $E_x$ as well as the Poisson's ratio $\nu_{xz}$, the general governing equations are simplifying considerably and lead to very simple closed-form solutions. The two antisymmetric and symmetric cases are again treated separately in the sequel.

### 2.5.1. Antisymmetric case

In this new configuration, the solutions of Equations (16(e)), (16(f)), (16(a)) and (16(c)), combined with the four connecting conditions (17), when assuming $W_a = W_b = \cos \frac{n\pi x}{L}$, take the following form:

$$ 
\left\{ 
\begin{aligned}
\alpha &= \frac{n\pi [n^4\pi^4 BG_{xz}h_c^3 + 3n^2\pi^2 L^2 BE_z h_c - 3L^4 E_z G_{xz}(h_c + h_s)]}{L[n^4\pi^4 AG_{xz}h_c^3 + 3n^2\pi^2 L^2 AE_z h_c + 3L^4 E_z G_{xz}]} \\
\beta &= -\alpha \\
k_1 &= \frac{n^2\pi^2 G_{xz} [n\pi(h_c+h_s)+\alpha L]}{2Lh_c(n^2\pi^2 G_{xz}h_c^2 + 3L^2 E_z)} \\
k_2 &= -\frac{n^3\pi^3 G_{xz}h_c^2(h_c+3h_s)+3n^2\pi^2\alpha L G_{xz}h_c^2 + 6n\pi L^2 E_z h_s + 6\alpha L^3 E_z}{2Lh_c(n^2\pi^2 G_{xz}h_c^2 + 3L^2 E_z)} \\
\zeta(z) &= k_1 z^3 + k_2 z \\
\xi(z) &= \frac{L}{n^3\pi^3 G_{xz}} [n^2\pi^2 G_{xz}(3k_1 z^2 + k_2) + 6L^2 k_1 E_z]
\end{aligned}
\right. 
\quad (30) 
$$

Then, Equations (16(b)) and (16(d)) lead to the new critical displacement:

$$ 
\lambda_{cr}^A = (n^2\pi^2 [n^4\pi^4 G_{xz}h_c^3(AD - B^2) + 3n^2\pi^2 L^2 E_z h_c(AD - B^2) \\
+ 3L^4 E_z G_{xz} (Ah_c^2 + 2h_c(Ah_s + B) + Ah_s^2 + 2Bh_s + D)]) / (AL(n^4\pi^4 AG_{xz}h_c^3 + 3n^2\pi^2 L^2 AE_z h_c + 3L^4 E_z G_{xz})) 
\quad (31) 
$$
### 2.5.2. Symmetric case

With $W_a = -\cos \frac{n\pi X}{L} = -W_b$, one obtains the following expressions:

$$ 
\left\{ 
\begin{aligned}
\alpha &= -\frac{n\pi B}{AL} \\
\beta &= \alpha \\
k_1 &= \frac{n\pi}{2Lh_c} \\
k_2 &= \frac{2\alpha L - n\pi(h_c + 2h_s)}{2L} \\
\zeta(z) &= k_1 z^2 + k_2 \\
\xi(z) &= \frac{2k_1 n^2 \pi^2 L G_{xz} z}{n^3 \pi^3 G_{xz}}
\end{aligned}
\right. 
\qquad (32) 
$$

and the following critical displacement:

$$ \lambda_{cr}^{S} = \frac{n^{4}\pi^{4}h_{c}(AD - B^{2}) + AL^{4}E_{z}}{n^{2}\pi^{2}A^{2}Lh_{c}} \qquad (33) $$

In this particular case, where the general expression of the critical displacement is quite simple, it is possible, by derivation of the critical displacement with respect to $n$, to get a closed-form expression of the real value of $n$ which cancels the derivative and thus minimizes $\lambda_{cr}^S$:

$$ n_{min} = \frac{L}{\pi} \sqrt[4]{\frac{AE_z}{h_c (AD - B^2)}} \qquad (34) $$

When dealing with symmetric modes, the minimum critical values of interest correspond to somewhat high wave numbers, so that a good approximation of the minimum critical displacement can be obtained by inserting this value of $n$ in Equation (33), which results in:

$$ \lambda_{cr-min}^{S} = 2L \sqrt{\frac{E_z (AD - B^2)}{A^3 h_c}} \qquad (35) $$

For comparison with existing formulae, the specific case of symmetric laminated faces is further addressed. In such a configuration, $B = 0$ and one can note $A = 2E_t h_s$ and
$D = 2E_b h_s^3/3$, where $E_t$ and $E_b$ represent the longitudinal traction and bending equivalent moduli of a face, respectively. If so, the critical displacement reduces to:

$$ \lambda_{cr-min}^{S} = \frac{L}{E_t} \sqrt{\frac{2E_z E_b h_s}{3h_c}} \qquad (36) $$

and the critical longitudinal stress can also be evaluated as follows:

$$ \sigma_{cr-min}^{S} = \frac{\lambda_{cr-min}^{S} E_t}{L} = \sqrt{\frac{2E_z E_b h_s}{3h_c}} \qquad (37) $$

This last result turns out to coincide perfectly with the critical stress of symmetric wrinkling obtained with a Winkler's model, defined by a simple beam on an elastic foundation (made of a distribution of transverse springs). More precisely, in the case of symmetric modes, the critical stress only depends on the thickness ratio, on the bending stiffness of the faces and on the transverse modulus of the core (and not of the shear modulus).

## 3. Numerical results, validation and analysis

### 3.1. Numerical computation of the critical displacements and bifurcation modes

Numerical finite element computations (linearized buckling analyses) have been performed, using Abaqus software, in order to validate the previous analytical solutions. In most cases, the sandwich column is described according to a 2D modeling approach. The core layer is represented by an (equivalent) homogeneous material with an orthotropic elastic constitutive law in the most general case, in accordance with the theoretical analysis above. The faces are also represented by a 2D continuous solid, separated in several layers (one for each ply) in the case of laminates, each with its own set of constitutive parameters depending on the isotropic or orthotropic nature of the plies. In both core layer and faces, the plane stress assumption is retained. The resulting finite element mesh is made up of 8-noded quadrangular elements with reduced integration. In each configuration, a short mesh convergence
analysis is achieved so as to retain the optimal mesh size, realizing a compromise between accuracy and computational cost. Similarly to the theoretical study, the loading is applied by means of a uniform compressive displacement at the right end of the column, whereas the left end is fixed in the longitudinal direction. A supplementary transverse displacement is fixed at an arbitrary point in order to prevent the structure from rigid modes.

## 3.2. Influence of the core orthotropy

A first general configuration is retained for validation purposes, with asymmetric laminated faces and an orthotropic core material. The geometric dimensions of the sandwich column are indicated in Table 1. The core material is orthotropic with the useful parameters listed in Table 2. The faces consist in two layers: one made up of an isotropic material, such as aluminium, and another one of adhesive material (also isotropic), at the interface with the core material. The thicknesses and constitutive parameters of the two materials are specified in Table 3, and the resulting stiffnesses of the upper face in Table 4 (those of the lower face are similar, except for *B* which changes sign, see Appendix A).

<table><caption>Table 1: Geometric dimensions</caption><thead><tr><th>L (mm)</th><th>b (mm)</th><th>h<sub>c</sub> (mm)</th></tr></thead><tbody><tr><td>500</td><td>100</td><td>25</td></tr></tbody></table>

<table><caption>Table 2: Core material properties</caption><thead><tr><th>E<sub>x</sub> (MPa)</th><th>E<sub>z</sub> (MPa)</th><th>ν<sub>xz</sub></th><th>G<sub>xz</sub> (MPa)</th></tr></thead><tbody><tr><td>10</td><td>300</td><td>0.01</td><td>115</td></tr></tbody></table>

The critical displacements of a hundred or so modes are depicted in Figure 3(b) against the half-wave number, resulting from both analytical and numerical calculations. In each
<table><caption>Table 3: Material and geometric parameters of the faces</caption><thead><tr><th>E<sub>alu</sub> (MPa)</th><th>ν<sub>alu</sub></th><th>t<sub>alu</sub> (mm)</th><th>E<sub>adh</sub> (MPa)</th><th>ν<sub>adh</sub></th><th>t<sub>adh</sub> (mm)</th></tr></thead><tbody><tr><td>70000</td><td>0.33</td><td>1</td><td>3000</td><td>0.3</td><td>0.25</td></tr></tbody></table>

<table><caption>Table 4: Stiffnesses of the upper laminated face (aluminium with adhesive)</caption><thead><tr><th>A (N/mm)</th><th>B (N)</th><th>D (N.mm)</th></tr></thead><tbody><tr><td>70751</td><td>8375</td><td>7119</td></tr></tbody></table>

case, two curves may be distinguished, corresponding respectively to antisymmetric and symmetric modes. In any case, three values are of particular interest, namely the minimum antisymmetric and symmetric wrinkling displacements (which correspond to a large half-wave number), and the global buckling value (associated to the antisymmetric mode with $n = 1$) which is likely to be lower than the previous values, depending on the sandwich configuration. In the present case, the minimum critical displacement corresponds to a wrinkling mode, what is more a symmetric mode. In the case of isotropic core materials, it was shown in [22] that the mode with the minimum critical value was systematically antisymmetric. This is clearly no longer the case when considering orthotropic core materials. More generally, the analytical and numerical results are in very good accordance, whatever the half-wave number. The new analytical solution related to the general case of an orthotropic core layer proves to be satisfactory in both cases of global buckling and wrinkling.

Next, the longitudinal core modulus $E_x$ is varied for further validation, which makes it possible to analyze the influence of this modulus on the buckling response of the sandwich structure. For convenience purposes, the Poisson's ratio $\nu_{xz}$ is also varied proportionally. Four values of $E_x$ are considered, including the previous one, but also the two particular cases $E_x = 0$ and $E_x = E_z$. The former case is typical of a thin-walled core structure (displaying a negligible longitudinal modulus) and the latter case corresponds to the limit case of an isotropic core material, due to the appropriate initial choice of the shear modulus
<figure><img src="image_5.png" /><figcaption>Figure 3: Analytical and numerical critical displacements (antisymmetric and symmetric modes) for various core longitudinal moduli</figcaption></figure>

$G_{xz}$. In the last two cases, expressions (31)-(33) and (27)-(29) can be used for the analytical derivation of the antisymmetric-symmetric critical displacements, respectively.

Figure 3 describes the same results as before for all the configurations considered. Generally speaking, the analytical and numerical results are still in very good accordance. There are only a few discrepancies in the isotropic case and even then, only for small half-wave numbers. It is certainly due to the greater value of the longitudinal modulus which makes the pre-critical longitudinal stresses no more negligible in the core layer, as it is assumed in the theoretical formulation. Besides, it can be shown that the critical displacements increase with the longitudinal modulus, especially in the range of low values. Conversely, the critical displacements are barely affected by the variation of this modulus in the range of the highest
values considered, when the constitutive behavior approaches that of an isotropic material. The critical values associated with symmetric modes show greater variations than the other ones, in such a way that a first symmetric wrinkling mode is observed with low values of $E_x$, whereas an antisymmetric mode occurs in the first instance with greater values. In the present case, the first mode even turns out to be global in the isotropic configuration. Actually, antisymmetric and symmetric wrinkling values always tend to coincide with increasing half-wave numbers, and in some cases (like this one), the minimum antisymmetric and symmetric critical values are similar, as long as the longitudinal modulus is sufficiently high. Figure 4 focuses on the three values of interest (the minimum antisymmetric and symmetric wrinkling displacements and the global buckling one), displaying their evolution with the longitudinal modulus, for clarity purposes. These results clearly show the importance of considering the longitudinal modulus in the derivation of the critical response of the sandwich structure.

<figure><img src="image_6.png" /><figcaption>Figure 4: Evolution of the minimum critical displacements according to the core longitudinal modulus</figcaption></figure>

In a second stage, the core shear modulus $G_{xz}$ is varied, starting from the previous reference case. Figure 5 displays again all the antisymmetric and symmetric critical displacements.
ments. Contrary to the longitudinal modulus, the shear modulus has a greater influence on the antisymmetric modes. Consequently, for large values of $G_{xz}$, a symmetric wrinkling mode is likely to occur in the first instance, whereas an antisymmetric mode (and even the global buckling mode) is observed first when considering lower values of the shear modulus. Again, the three values of interest are depicted in Figure 6 according to the value of $G_{xz}$, for the sake of clarity. In any case, the analytical and numerical results are still in very good accordance.

<figure><img src="image_7.png" /><figcaption>Figure 5: Analytical and numerical critical displacements (antisymmetric and symmetric modes) for various core shear moduli</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 6: Evolution of the minimum critical displacements according to the core shear modulus</figcaption></figure>

### 3.3. Influence of the adhesive layer

Finally, the reference configuration is considered again, and the adhesive layer is removed, so as to measure its influence (and more generally the influence of the stiffnesses of the laminated faces) on the buckling response of the sandwich column. In the absence of adhesive, the new stiffnesses take the values listed in Table 5, to be compared to the reference values of Table 4. In this new configuration, the symmetric faces give rise to a null value of the coupling stiffness *B*. Besides, the extensional stiffness *A* is quite similar to the previous one (since the adhesive material has a very low modulus, relative to aluminium). Conversely, the bending stiffness is comparatively lower.

<table><caption>Table 5: Stiffnesses of the isotropic faces (aluminium without adhesive)</caption><thead><tr><td>A (N/mm)</td><td>B (N)</td><td>D (N.mm)</td></tr></thead><tbody><tr><td>70000</td><td>0</td><td>5833</td></tr></tbody></table>

For further validation, the same problem is analyzed (with and without the adhesive layer), considering a thickness of the aluminium layer equal to half of the previous one, namely $t_{alu} = 0.5$ mm. All the critical displacements are plotted again versus the half-wave
number in Figure 7 and the relative influence of the adhesive layer on the possible minimum critical values is quantified in Figure 8.

<figure><img src="image_9.png" /><figcaption>Figure 7: Analytical and numerical critical displacements (antisymmetric and symmetric modes) for various face thicknesses, with or without adhesive</figcaption></figure>

It appears that, with this new half thickness, the critical half-wave numbers have risen sharply. However, the critical displacements are about the same value (but the critical forces should be significantly reduced due to the reduction of the face stiffnesses by approximately half). Furthermore, the global critical displacement is shown to be rather not influenced by the presence of the adhesive layer, whatever the face thickness. Conversely, despite the low modulus of the adhesive material, the adhesive layer has a non negligible influence on the wrinkling minimum values, especially when the thickness decreases.
<figure><img src="image_10.png" /><figcaption>Figure 8: Influence of the adhesive layer on the critical displacements</figcaption></figure>

### 3.4. Case study: honeycomb core sandwich

Finally, the particular case of a honeycomb core sandwich is investigated. The objective is to confront the present analytical solutions to reference numerical results and still discuss the merits of the various formulae. In this instance, the thin-walled core material is described explicitly using shell finite elements, and no longer in an homogenized way. A 3D model of
the sandwich column is thus developed where the faces are also represented by shell elements.

The column is 335 *mm* long and 69.6 *mm* wide. The honeycomb core layer, with a total thickness of 80 *mm*, is made up of regular hexagons with sides of 3.87 *mm* and a uniform thickness of cell walls equal to 0.25 *mm* (see Figure 9 for a general representation of the sandwich column). The honeycomb material is aramid, which is supposed to be isotropic with Young's modulus 5000 *MPa* and Poisson's ratio 0.25. From all these geometric and material parameters, one can deduce easily the required equivalent core moduli for use in the analytical expressions. The usual expressions of these moduli, which can be found in [27], are recalled in Appendix B. The resulting values are listed in Table 6. For verification purposes, the four equivalent moduli were also computed using the same finite element model as for the buckling analysis (limited here to the core layer). Macroscopic forces were applied to the appropriate faces of the part, under suitable boundary conditions, so as to induce a longitudinal traction, transverse traction and transverse shear state, respectively, as illustrated in Figure 10. The longitudinal traction case (a) allows one to derive the longitudinal Young's modulus $E_x$ from the longitudinal maximum displacement. The transverse traction case (b) gives rise similarly to the transverse Young's modulus $E_z$, together with the Poisson's ratio $\nu_{xz}$ which can be deduced from the longitudinal shortening. Finally, the last loading case (c) provides the equivalent transverse shear modulus $G_{xz}$. The corresponding numerical values are also listed in Table 6, for comparison purposes, but they are not used in the sequel. One can mention, however, that the agreement between the analytical and numerical equivalent moduli is very good, despite unavoidable edge effects in the numerical computations.

The faces are made up of 4 plies of carbon/epoxy unidirectional composites, whose geometric and material parameters are given in Table 7, with a symmetric stacking sequence
<figure><img src="image_11.png" /><figcaption>Figure 9: Exploded view of the honeycomb core sandwich</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 10: Numerical derivation of the equivalent moduli of the core layer</figcaption></figure>

<table><caption>Table 6: Honeycomb core homogenized properties</caption><thead><tr><td>Equivalent core moduli</td><td>E<sub>x</sub> (MPa)</td><td>E<sub>z</sub> (MPa)</td><td>ν<sub>xz</sub></td><td>G<sub>xz</sub> (MPa)</td></tr></thead><tbody><tr><td>Analytical</td><td>3.12</td><td>373</td><td>0.00209</td><td>74.6</td></tr><tr><td>Numerical</td><td>3.24</td><td>374</td><td>0.00212</td><td>70.5</td></tr></tbody></table>

(0/90/90/0). Using the classical laminate theory (see Appendix A), one can deduce the relevant stiffnesses for further analytical use, which are indicated in Table 8.
<table><caption>Table 7: Material and geometric parameters of the plies</caption><thead><tr><th>E<sub>1</sub> (MPa)</th><th>E<sub>2</sub> (MPa)</th><th>ν<sub>12</sub></th><th>G<sub>12</sub> (MPa)</th><th>t<sub>ply</sub> (mm)</th></tr></thead><tbody><tr><td>230000</td><td>13000</td><td>0.32</td><td>5200</td><td>0.5</td></tr></tbody></table>

<table><caption>Table 8: Stiffnesses of the laminated faces (carbon/epoxy)</caption><thead><tr><th>A (N/mm)</th><th>B (N)</th><th>D (N.mm)</th></tr></thead><tbody><tr><td>244128</td><td>0</td><td>135748</td></tr></tbody></table>

Analytical and numerical calculations are performed and a set of first modes are examined. From the analytical point of view, the general solutions are used, involving the longitudinal core modulus and the Poisson's ratio, and the simplified closed-form solutions are also used, considering these two parameters to be zero. The antisymmetric and symmetric critical values are plotted against the half-wave number in Figure 11. For illustration purposes, typical antisymmetric and symmetric modes are also represented in Figure 12.

First of all, the influence of the longitudinal modulus $E_x$ on the critical values depends on the nature of the mode. As observed earlier, the antisymmetric critical displacements are virtually unaffected by this parameter, whereas the symmetric ones are more likely to be impacted. Nevertheless, the maximum relative difference between the two approaches (with and without considering this modulus in the formulation) hardly reaches 10%. Besides, the analytical solutions turn out to be in good agreement with the numerical results. The minimum critical value is associated to a symmetric mode, as could be expected with such an orthotropic ratio of the core material. For this particular mode, the analytical value matches perfectly with the numerical one. Furthermore, the analytical and numerical curves tend to move away from each other as the half-wave number increases. This could be due to edge effects, due to the small number of cells considered in the model, or to the mesh size, whose convergence was solely based on the first mode. It is important to mention that other
<figure><img src="image_13.png" /><figcaption>(a) Antisymmetric modes</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(b) Symmetric modes</figcaption></figure>

Figure 11: Comparison between analytical and numerical critical values of a honeycomb sandwich
<figure><img src="image_15.png" /><figcaption>Figure 12: Representative modes of a honeycomb structure</figcaption></figure>

types of mode, mainly involving the thin-walled core structure, are likely to occur between
the wrinkling modes of the faces. What is more, the antisymmetric and symmetric wrinkling
modes, such as the ones depicted in Figure 12, can probably also be seen as more general
buckling modes involving both faces and cell walls. In such circumstances, the non-linear
homogenized response of the core layer cannot be theoretically represented using only the
estimated equivalent moduli in a linear framework. However, such a simplified approach
may be sufficient in many configurations, as illustrated in the present case study.

**4. Conclusions**

In this paper, the buckling behavior of sandwich columns under axial compression was
investigated in an analytical way, from a rigorous bifurcation analysis. A 2D model was
defined, where the skins were assumed to behave like Euler-Bernoulli beams while the core
material was modeled as a 2D continuous solid, without considering any simplified defor-
mation field. The general case of heterogeneous (laminated) faces and an orthotropic core
material was considered, giving rise to original analytical solutions. Then, these solutions
were particularized in the simple case of an isotropic core material and the current case of
a thin-walled core structure, considering after homogenization that the associated longitu-
dinal modulus is negligible. Closed-form solutions were obtained in these specific cases for the antisymmetric and symmetric wrinkling modes, also encompassing the global buckling mode.

The overall objective of the paper is to provide some original analytical solutions that can replace cumbersome and time-consuming finite element calculations. All the previous solutions were thus validated against 2D numerical finite element results. In every case, the analytical and numerical critical values were shown to be very similar. More specifically, the same order of appearance of the wrinkling modes was observed, with the same half-wave number (and the same critical displacement) for the first and only mode of interest (with the minimum critical value).

Concerning the core material, a short parametric analysis was performed, by varying successively the longitudinal and shear moduli. It was shown that, contrary to the isotropic case, the first mode is likely to be symmetric, as soon as the longitudinal modulus is sufficiently low or the shear modulus quite significant. The evolution of the possible minimum critical values (corresponding to the global buckling mode, or to the minimum antisymmetric and symmetric wrinkling ones) with respect to these particular moduli was emphasized, showing transitions between the prediction of global buckling and antisymmetric or symmetric wrinkling. In most circumstances, the consideration of the full 2D anisotropic constitutive law of the (equivalent) core material is required to achieve reliable solutions.

Dealing now with the skins, the influence of an adhesive layer was illustrated, with various face thicknesses. This influence is naturally greater the thinner the skin, especially for local wrinkling modes.

Lastly, a case study was addressed about the buckling response of a honeycomb sandwich column. In this particular configuration, a numerical shell model was defined so as to describe explicitly the thin-walled structure of the core layer. Again, the analytical and numerical results were compared to each other. Despite many approximations and the fail-
ure to take account of non-linear local buckling phenomena involving the honeycomb cell walls, the present general and simplified analytical solutions allow one to predict successfully the minimum critical value and the associated mode of such a strongly heterogeneous and orthotropic structure.

## Data availability

The raw/processed data required to reproduce these findings cannot be shared at this time as the data also forms part of an ongoing study.

## References

[1] Hoff, N.J. and Mautner, S.E., Buckling of sandwich type panels, *Journal of the Aeronautical Sciences* **12** Issue 3 (1945) 285-297.

[2] Allen, H.G., *Analysis and Design of Structural Sandwich Panels* (Pergamon, 1969).

[3] Benson, A.S. and Mayers, J., General instability and face wrinkling of sandwich plates: Unified theory and applications, *AIAA Journal* **5** Issue 4 (1967) 729-739.

[4] Léotoing, L., Drapier S. and Vautrin, A., First applications of a novel unified model for global and local buckling of sandwich columns, *European Journal of Mechanics - A/Solids* **21** Issue 4 (2002) 683-701.

[5] Jasion, P., Magnucka-Blandzi, E., Szyc, W. and Magnucki, K., Global and local buckling of sandwich circular and beam-rectangular plates with metal foam core, *Thin-Walled Structures* **61** (2012) 154-161.

[6] Niu, K. and Talreja, R., Modeling of wrinkling in sandwich panels under compression, *Journal of Engineering Mechanics* **125** Issue 8 (1999) 875-883.

[7] Norris, C.B., Boller, K.H. and Voss, A.W., Wrinkling of the facings of sandwich construction subjected to edgewise compression, *FPL Report No. 1810-A* (1953).

[8] Pearce, T.R.A. and Webber, J.P.H., Experimental buckling loads of sandwich panels with carbon fibre faceplates, *Aeronautical Quarterly* **24** (1973) 295-312.

[9] Webber, J.P.H., Kyriakides, S. and Lee, C.T., On the wrinkling of honeycomb sandwich columns with laminated cross-ply faces, *The Aeronautical Journal* **80** Issue 786 (1976) 264-272.
[10] Koissin, V., Shipsha, A. and Skvortsov, V., Effect of physical nonlinearity on local buckling in sandwich beams, *Journal of Sandwich Structures and Materials* **12** Issue 4 (2010) 477–494.

[11] Fagerberg, L., Wrinkling and compression failure transition in sandwich panels, *Journal of Sandwich Structures and Materials* **6** Issue 2 (2004) 129–144.

[12] Ley, R.P., Lin, W. and Mbanefo, U., Facesheet wrinkling in sandwich structures, *Technical Report NASA/CR-1999-208994* (1999).

[13] Kim, C.G. and Hong, C.S., Buckling of unbalanced anisotropic sandwich plates with finite bonding stiffness, *AIAA Journal* **26** Issue 8 (1988) 982–988.

[14] Hadi, B.K. and Matthews, F.L., Development of Benson-Mayers theory on the wrinkling of anisotropic sandwich panels, *Composite Structures* **49** Issue 4 (2000) 425–434.

[15] Vonach, W.K. and Rammerstorfer, F.G., Wrinkling of thick orthotropic sandwich plates under general loading conditions, *Archive of Applied Mechanics* **70** (2000) 338–348.

[16] Fagerberg, L. and Zenkert, D., Effects of anisotropy and multiaxial loading on the wrinkling of sandwich panels, *Journal of Sandwich Structures and Materials* **7** (2005) 177–194.

[17] D'Ottavio, M. and Polit, O., Linearized global and local buckling analysis of sandwich struts with a refined quasi-3D model, *Acta Mechanica* **226** (2015) 81–101.

[18] Vescovini, R., D'Ottavio, M., Dozio, L. and Polit, O., Buckling and wrinkling of anisotropic sandwich plates, *International Journal of Engineering Science* **130** (2018) 136–156.

[19] Fagerberg, L., The effect of local bending stiffness on the wrinkling of sandwich panels, *Proceedings of the Institution of Mechanical Engineers, Part M: Journal of Engineering for the Maritime Environment* **217** (2003) 111–119.

[20] Wadee, M.A. and Hunt, G.W., Interactively induced localized buckling in sandwich structures with core orthotropy, *Journal of Applied Mechanics* **65** (1998) 523–528.

[21] Stiftinger, M.A. and Rammerstorfer, F.G., Face layer wrinkling in sandwich shells – Theoretical and experimental investigations, *Thin-Walled Structures* **29** Issues 1-4 (1997) 113–127.

[22] Douville, M.A. and Le Grognec, P., Exact analytical solutions for the local and global buckling of sandwich beam-columns under various loadings, *International Journal of Solids and Structures* **50** Issues 16-17 (2013) 2597–2609.

[23] Le Grognec, P., Nguyen, Q.H. and Hjiaj, M., Exact buckling solution for two-layer Timoshenko beams with interlayer slip, *International Journal of Solids and Structures* **49** Issue 1 (2012) 143–150.

[24] Lainé, C., Le Grognec, P., Comas Cardona, S. and Binetruy, C., Analytical, numerical and experimental
study of the bifurcation and collapse behavior of a 3D reinforced sandwich structure under through-thickness compression, *International Journal of Mechanical Sciences* **67** (2013) 42–52.

[25] Sad Saoud, K. and Le Grognec, P., A unified formulation for the biaxial local and global buckling analysis of sandwich panels, *Thin-Walled Structures* **82** (2014) 13–23.

[26] Nguyen, Q.S., *Stability and Non-Linear Solid Mechanics* (Wiley, 2000).

[27] Gibson, L.J. and Ashby, M.F., *Cellular Solids: Structure and Properties* (Cambridge University Press, 1997).

## Appendix A. Derivation of the extensional, bending and coupling stiffnesses of the faces

Let us consider the upper face of the sandwich column, which is supposed to be composed of *n* layers made up of orthotropic materials (or, equivalently, of transversely isotropic materials, as long as one assumes plane stresses). The particular case of isotropic materials can naturally be investigated by setting properly the respective moduli. The thickness of each layer *k* is defined as $h_{k+1} - h_k$, where $h_1 = -h_s$, $h_{n+1} = h_s$ and $h_k$ for $k = 2, n$ refers to the through-thickness position of the interface between layers $k-1$ and $k$ relative to the mid-surface (see Figure A.13).

<figure><img src="image_16.png" /><figcaption>Figure A.13: Stacking sequence of a laminated composite face</figcaption></figure>

Each layer *k* (from 1 to *n*) can be viewed as an unidirectional composite plate in plane stresses, with an angle $\theta_k$ between global axis *x* and local axis 1 in the fiber direction (see
<figure><img src="image_17.png" /><figcaption>Figure A.14: Oriented unidirectional ply</figcaption></figure>

Figure A.14). In the in-plane material basis (1, 2), the reduced 2D constitutive law of the UD material can be written as follows, using the classical notations:

$$ \left\{ \begin{matrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{12} \end{matrix} \right\}_k = \begin{bmatrix} Q_{11} & Q_{12} & 0 \\ Q_{12} & Q_{22} & 0 \\ 0 & 0 & Q_{66} \end{bmatrix}_k \left\{ \begin{matrix} \varepsilon_{11} \\ \varepsilon_{22} \\ \gamma_{12} \end{matrix} \right\} \quad (A.1) $$

with:

$$ Q_{11} = \frac{E_1}{1 - \frac{E_2}{E_1} \nu_{12}^2} \qquad Q_{22} = \frac{E_2}{E_1} Q_{11} \quad (A.2) $$

$$ Q_{12} = \nu_{12} Q_{11} \qquad Q_{66} = G_{12} $$

where $E_1$ and $E_2$ are respectively the Young's moduli in the fiber and transverse directions, and $\nu_{12}$ and $G_{12}$ are the associated Poisson's ratio and shear modulus.

After rotation in the global basis $(x, y)$, the same constitutive law can be rewritten in the following form:

$$ \left\{ \begin{matrix} \sigma_{xx} \\ \sigma_{yy} \\ \sigma_{xy} \end{matrix} \right\}_k = \begin{bmatrix} Q'_{11} & Q'_{12} & Q'_{16} \\ Q'_{12} & Q'_{22} & Q'_{26} \\ Q'_{16} & Q'_{26} & Q'_{66} \end{bmatrix}_k \left\{ \begin{matrix} \varepsilon_{xx} \\ \varepsilon_{yy} \\ \gamma_{xy} \end{matrix} \right\} = [\mathbf{Q}']_k \left\{ \begin{matrix} \varepsilon_{xx}^0 \\ \varepsilon_{yy}^0 \\ \varepsilon_{xy}^0 \end{matrix} \right\} + z[\mathbf{Q}']_k \left\{ \begin{matrix} \kappa_{xx} \\ \kappa_{yy} \\ \kappa_{xy} \end{matrix} \right\} \quad (A.3) $$
The components of matrix $[Q']_k$ involved in Equation (A.3) can be expressed in terms of their counterparts in the local basis:

$$
\begin{align*}
Q'_{11} &= Q_{11} \cos^4 \theta_k + Q_{22} \sin^4 \theta_k + 2(Q_{12} + 2Q_{66}) \cos^2 \theta_k \sin^2 \theta_k \\
Q'_{22} &= Q_{11} \sin^4 \theta_k + Q_{22} \cos^4 \theta_k + 2(Q_{12} + 2Q_{66}) \cos^2 \theta_k \sin^2 \theta_k \\
Q'_{12} &= Q_{12} (\cos^4 \theta_k + \sin^4 \theta_k) + (Q_{11} + Q_{22} - 4Q_{66}) \cos^2 \theta_k \sin^2 \theta_k \tag{A.4} \\
Q'_{66} &= Q_{66} (\cos^4 \theta_k + \sin^4 \theta_k) + (Q_{11} + Q_{22} - 2(Q_{12} + Q_{66})) \cos^2 \theta_k \sin^2 \theta_k \\
Q'_{16} &= (Q_{11} - Q_{12} - 2Q_{66}) \cos^3 \theta_k \sin \theta_k + (Q_{12} - Q_{22} + 2Q_{66}) \cos \theta_k \sin^3 \theta_k \\
Q'_{26} &= (Q_{11} - Q_{12} - 2Q_{66}) \cos \theta_k \sin^3 \theta_k + (Q_{12} - Q_{22} + 2Q_{66}) \cos^3 \theta_k \sin \theta_k
\end{align*}
$$

Besides, vectors {$\boldsymbol{\epsilon}^0$} and {$\boldsymbol{\kappa}$} in Equation (A.3) contain the in-plane membrane and
curvature strain components, respectively, which are the same for all layers. According to
the classical laminate theory, one can define the dual generalized stress vectors, namely the
normal force and bending moment vectors:

$$
\{\mathbf{N}\} = \sum_{k=1}^{n} \int_{h_k}^{h_{k+1}} \left\{ \begin{matrix} \sigma_{xx} \\ \sigma_{yy} \\ \sigma_{xy} \end{matrix} \right\}_k dz \qquad \{\mathbf{M}\} = \sum_{k=1}^{n} \int_{h_k}^{h_{k+1}} \left\{ \begin{matrix} \sigma_{xx} \\ \sigma_{yy} \\ \sigma_{xy} \end{matrix} \right\}_k z \, dz \quad (\text{A.5})
$$

At the end, the generalized constitutive law of the composite plate writes:

$$
\left\{ \begin{matrix} \mathbf{N} \\ \mathbf{M} \end{matrix} \right\} = \left[ \begin{matrix} \mathbf{A} & \mathbf{B} \\ \mathbf{B} & \mathbf{D} \end{matrix} \right] \left\{ \begin{matrix} \boldsymbol{\epsilon}^0 \\ \boldsymbol{\kappa} \end{matrix} \right\} \qquad (\text{A.6})
$$
namely:

$$ \left\{ \begin{array}{c} N_{xx} \\ N_{yy} \\ N_{xy} \\ M_{xx} \\ M_{yy} \\ M_{xy} \end{array} \right\} = \left[ \begin{array}{cccccc} A_{11} & A_{12} & A_{16} & B_{11} & B_{12} & B_{16} \\ A_{12} & A_{22} & A_{26} & B_{12} & B_{22} & B_{26} \\ A_{16} & A_{26} & A_{66} & B_{16} & B_{26} & B_{66} \\ B_{11} & B_{12} & B_{16} & D_{11} & D_{12} & D_{16} \\ B_{12} & B_{22} & B_{26} & D_{12} & D_{22} & D_{26} \\ B_{16} & B_{26} & B_{66} & D_{16} & D_{26} & D_{66} \end{array} \right] \left\{ \begin{array}{c} \varepsilon_{xx}^0 \\ \varepsilon_{yy}^0 \\ \varepsilon_{xy}^0 \\ \kappa_{xx} \\ \kappa_{yy} \\ \kappa_{xy} \end{array} \right\} \quad (A.7) $$

$\{A\}$, $\{B\}$ and $\{D\}$ are respectively the membrane, coupling, and bending stiffness matrices, whose components are defined as follows:

$$ \begin{aligned} A_{ij} &= \sum_{k=1}^{n} (Q'_{ij})_k (h_{k+1} - h_k) \\ B_{ij} &= \frac{1}{2} \sum_{k=1}^{n} (Q'_{ij})_k (h_{k+1}^2 - h_k^2) \\ D_{ij} &= \frac{1}{3} \sum_{k=1}^{n} (Q'_{ij})_k (h_{k+1}^3 - h_k^3) \end{aligned} \quad (A.8) $$

The coupling traction-bending terms $B_{11}$, $B_{22}$, $B_{12}$, $B_{66}$, traction-shear terms $A_{16}$, $A_{26}$, traction-torsion terms $B_{16}$, $B_{26}$ and bending-torsion terms $D_{16}$, $D_{26}$ should be zero in practical cases of symmetric and/or balanced stacking sequences. Nevertheless, in the present case of a composite beam along the x-axis, it is necessary to reduce again the constitutive law, due to the zero out-of-plane generalized stresses ($N_{yy} = N_{xy} = 0$ and $M_{yy} = M_{xy} = 0$). A simple way of doing this is to invert the total stiffness matrix in Equation (A.7), to truncate the resulting flexibility matrix only to the longitudinal terms $\bullet_{xx}$ (or $\bullet_{11}$) and to invert it again. The resulting reduced constitutive law writes finally:

$$ \left\{ \begin{array}{c} N_{xx} \\ M_{xx} \end{array} \right\} = \left[ \begin{array}{cc} A & B \\ B & D \end{array} \right] \left\{ \begin{array}{c} \varepsilon_{xx}^0 \\ \kappa_{xx} \end{array} \right\} \quad (A.9) $$

The reduced stiffnesses A, B and D will serve as the only material parameters of the
upper face in the buckling analyses. As the sandwich column is symmetric, the lower face is the same, but with a stacking sequence in the opposite direction. In practice, the same traction and bending stiffnesses A and D as for the upper face can be used, whereas an opposite coupling term −B must be considered.

## Appendix B. Equivalent moduli of a simple honeycomb structure

In the case study of this paper, a simple honeycomb structure is considered with a periodic microstructure and a uniform wall thickness over the entire core layer. The geometry of the considered unit cell is represented in Figure B.15. It is defined by the wall thickness t, the dimensions h and l, and the characteristic angle θ. In the present study, use is made of the equivalent longitudinal modulus E<sub>x</sub>, transverse (out-of-plane) modulus E<sub>z</sub>, Poisson's ratio ν<sub>xz</sub> and transverse shear modulus G<sub>xz</sub> of the honeycomb structure, which classically take the following expressions [27]:

$$
\begin{align*}
E_x &= E_s \left(\frac{t}{l}\right)^3 \frac{\cos \theta}{\left(\frac{h}{l} + \sin \theta\right) \sin^2 \theta} \\
E_z &= E_s \left(\frac{t}{l}\right) \frac{\frac{h}{l} + 2}{2\left(\frac{h}{l} + \sin \theta\right) \cos \theta} \tag{B.1} \\
\nu_{xz} &= \nu_s \left(\frac{t}{l}\right)^2 \frac{2 \cos^2 \theta}{\left(\frac{h}{l} + 2\right) \sin^2 \theta} \\
G_{xz} &= G_s \left(\frac{t}{l}\right) \frac{\cos \theta}{\left(\frac{h}{l} + \sin \theta\right)}
\end{align*}
$$

where E<sub>s</sub> and ν<sub>s</sub> are the Young's modulus and the Poisson's ratio of the cell walls, respectively.

When the unit cell is a regular hexagon (h = l and θ = π/6), as is the case in this study, the structure is isotropic in the xy-plane and the formulae above are getting simpler.
<figure><img src="image_18.png" /><figcaption>Figure B.15: Unit cell of the honeycomb structure</figcaption></figure>