Journal of Computational Design and Engineering, 2025, 12, 107–126

DOI: 10.1093/jcde/qwaf067
Advance access publication date: 22 July 2025
Research Article

Enhancing buckling resistance of strut lattice structures

via three-dimensional topology optimization
Asha Viswanath1 , Wesley James Cantwell1,2 and Kamran Ahmed Khan 1,2 ,
*
1
Department of Aerospace Engineering, Khalifa University of Science and Technology, Abu Dhabi PO Box 127788, UAE
2
Advanced Digital & Additive Manufacturing Group, Khalifa University of Science and Technology, Abu Dhabi PO Box 127788, UAE
∗
Correspondence: kamran.khan@ku.ac.ae

Downloaded from https://academic.oup.com/jcde/article/12/8/107/8210400 by guest on 07 September 2025

Abstract
Buckling failure is the dominant failure mode in strut-based or thin-walled lattice structures, limiting their effectiveness in structural
applications requiring high stiffness and strength. This paper investigates a novel buckling-based 3D topology optimization method
to generate buckling-resistant strut-based lattice structure for the first time. A MATLAB code for 3D topology optimization, based
on buckling objective/constraint, is developed inspired by a 2D code available in literature. A single cell from a simple cubic lattice
with square cross-sectional struts, a body-centred cubic lattice with circular struts, and a multi-cell 2 × 2 simple cubic lattice are
investigated as reference structures. Through analysing their buckling modes, the buckling-prone members are actively optimized
via topology optimization exercise to maximize buckling load factors, whilst constraining the volume fraction and stiffness of the
structure. The results demonstrate that topology-optimized designs exhibit higher buckling load factors than the reference structure.
Numerical validation of linear buckling analysis using ABAQUS software is also performed. Few parametric studies with varying cell
sizes and relative densities are also discussed. This methodology can be easily extended to enhance the buckling resistance of any
lattice structure without altering the relative density of the corresponding reference structure. Future directions include numerical
post-buckling studies with experimental validation.

Keywords: simple cubic lattice, body-centred cubic lattice, lattice structure buckling, topology optimization, MATLAB code, finite-
element analysis

Nomenclature
ABBREVIATIONS
LS: Lattice structure
SC: Simple cubic
TO: Topology optimization
FE: Finite element
SIMP: Solid isotropic material with penalization
BLF: Buckling load factor
BC: Boundary condition
CAD: Computer-aided design
STL: Stereo lithography

NOTATIONS
Variables
: Structure design domain
m: number of FEs
e : Element domain
n: Total dofs
ρe : Element density
ρe : Filtered-element density
ρe : Pseudo-element density
η,β : Projection parameters
rmin : Density filter
pK : Penalization factor for stiffness
pG : Penalization factor for stress

Received: March 6, 2025. Revised: June 3, 2025. Accepted: June 7, 2025

© The Author(s) 2025. Published by Oxford University Press on behalf of the Society for Computational Design and Engineering. This is an Open Access article
distributed under the terms of the Creative Commons Attribution-NonCommercial License (https://creativecommons.org/licenses/by-nc/4.0/), which permits
non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact
reprints@oup.com for reprints and translation rights for reprints. All other permissions can be obtained through our RightsLink service via the Permissions
link on the article page on our site-for further information please contact journals.permissions@oup.com
108 | Enhancing buckling resistance via three-dimensional topology optimization

EK : SIMP interpolation of Young’s modulus for stiffness matrix

EG : SIMP interpolation of Young’s modulus for stress matrix
Emin : Void element Young’s modulus
Es : Solid element Young’s modulus
A: Active domain region
δ: Dirac delta function
Vf , c f : Maximum allowed volume fraction, compliance

Vectors
ρ: Element density vector
ρˆ : Pseudo-element density vector
ρˆ A : Active pseudo-density vector
ρˆ P1 : Passive pseudo-densities vector 1
ρˆ P0 : Passive pseudo-densities vector 0

Downloaded from https://academic.oup.com/jcde/article/12/8/107/8210400 by guest on 07 September 2025

F: Load vector
Fcr : Critical load vector
u: Displacement vector
v: Eigenvector
ϕ: Normalized eigenvector
wi : Adjoint vector

Sets
Rn : n tuple of real numbers
Rnxn : Size n real square matrices set

Matrices
K: Linear stiffness matrix
Kσ : Stress stiffness matrix

Subscripts
KS: Kreisselmeier and Steinhauser
e : Element
s : Solid

PSEUDO-CODE VARIABLES
nelx, nely: FE mesh size
penalK : Penalization factor for stiffness
penalG : Penalization factor for stress
rmin : Density filter radius
P1 : Passive domain region with 1 density
P0 : Passive domain region with 0 density
V: Total volume fraction
c: Compliance
λi : ith buckling mode
λ1 : Fundamental BLF
μ: Eigenvalue
JKS : Aggregation function
pAGG : Aggregation parameter
J0KS : Aggregation function for μ1
J1KS : Aggregation function for constraints
gv, gc : Constraints on volume fraction, compliance
B: Set of buckling modes
ft : Projection filtering scheme with beta, eta—projection parameters
ftBC : Filter BCs, Neumann or Dirichlect
ocPar : Optimality criteria parameters
maxit : Maximum number of iterations
Lx : Physical length of the domain
nEig : Number of BLFs in optimization
pAgg : Initial value of KS aggregation factor
prSel : Data structure specifying optimization problem
Cmax /C0 : Ratio of maximum allowable compliance to initial compliance
start, maxP, ds, dP : Continuation parameters
 Journal of Computational Design and Engineering, 2025, 12(8), 107–126 | 109

1. Introduction
Lattice structures (LSs) are increasingly used in various structural engineering disciplines over recent decades, due to their strength
and lightweight properties (Mantovani et al., 2021; Ramadani et al., 2021; Khan & Riccio, 2024; Khan et al., 2024; Nogueira et al., 2024;
Yang et al., 2024). Topology optimization (TO) is particularly useful for creating intelligent designs based on these lightweight lattice
components in the aerospace and automotive industries by providing a systematic mathematical design approach (Seharing et al., 2020;
Zhu et al., 2021; Prathyusha & Raghu Babu, 2022). However, TO applied to LSs having lower volume fractions can result in very thin or
slender members that are susceptible to buckling prior to material failure. This may be because topology-optimized designs, achieved
via a standard minimum compliance approach (Bendsoe & Sigmund, 2013; Gao et al., 2023; Lin et al., 2025), do not account for buckling
but result in tension/compression-dominated configurations, avoiding bending members. As the buckling load of a structure is directly
proportional to its bending stiffness (for the case of a simple Euler column), these structures fail from buckling rather than fracture.
Consequently, designing an LS for improved buckling strength is crucial to ensure that specific load-bearing requirements are satisfied
in lightweight structural design. The methodology of shape optimization (Zhang et al., 2018) and structural optimization of LSs (Wang
et al., 2021) members to improve buckling have been researched. However, TO studies with the goal of improving buckling strength are

Downloaded from https://academic.oup.com/jcde/article/12/8/107/8210400 by guest on 07 September 2025

limited in the literature.
Buckling-based TO remains a challenging topic, presenting a number of difficulties (Bruyneel et al., 2008), even in its linearized formu-
lation. Initial attempts at to include buckling into the TO of LSs involved including local stability constraints (Neves et al., 2002) in the
minimum compliance TO problem or maximizing the buckling strength of the buckling modes (Thomsen et al., 2018). However, due to
the large number of eigenvalues that must be computed, leading to a prohibitively high computational burden, the buckling analysis is
usually performed as a post-optimization step, rather than an integrated optimization constraint. If the elastic stability of the compo-
nent is found insufficient, post-processing is used to improve the minimum buckling load; however, this process may lead to suboptimal
components. More recently, progress has been made within the MATLAB environment by introducing multi-level solvers for eigenprob-
lem solutions and fully vectorized strategy for the stress stiffness matrix set-up, leading to MATLAB educational codes with improved
efficiency being developed for the 2D buckling-based TO for all types of structures (Ferrari et al., 2021). These codes can also be easily
extended to 3D cases and applied to both LSs and microstructures. However, since the solution of a 3D eigenproblem could become
prohibitive with built-in MATLAB tools, it was never published in literature. Subsequent studies investigating the stability of 3D LSs have
used 2D TO formulation to design buckling-resistant topologies, which are then extruded in the third dimension for manufacturing
and experimental testing (Bluhm et al., 2022; Viswanath et al., 2024). Christensen et al. (2023) investigated the global and local buckling
response of multiscale structures in a 2D framework. Xu et al. (2024) developed buckling-constrained 3D structures by considering a
linear material model for compliance minimization and reformulating the problem with a projection method with minimal parameter
tuning. This did not provide a detailed 3D implementation procedure for readers. Wang & Sigmund (2021) used 3D formulations to de-
sign architected materials with an enhanced buckling by minimizing the weighted stiffness and buckling response, however they did not
discuss any related MATLAB codes. This work stressed the relevance of a 3D buckling-based TO formulation of isotropic microstructures,
which becomes even more important for anisotropic LSs. Wang et al. (2023) created three plate unit cells at different hierarchical scales
and compared their performance with hollow truss-based structures, concluding that the non-hierarchical structure outperformed the
hollow truss structure by a factor of 2.4 in terms of buckling strength and by factors of 5 and 1.4 for first- and second-order structures,
respectively. Christensen et al. (2025) dealt with TO-based built-in pre-failure indicators within multiscale structures, which can provide
early warning of local buckling without compromising structural integrity.
TO with a buckling-focused objective, applied to LSs, can result in structures that are optimized in terms of material distribution, while
providing the maximum stiffness alongside improved buckling resistance of LS. Considering the computational burden of applying the
optimization to multiple cells of the LS, one should first ascertain if a single-cell LS can effectively be optimized for buckling and if this
microstructure design can be exploited to build the complete cellular structure. Recently, Viswanath et al. (2023) implemented an optimal
single-cell design for a square honeycomb cell structure subjected to in-plane buckling using 2D TO by modifying the cross-section of
the load-bearing walls while keeping the volume fraction of the structure constant. However, for strut-based lattices or plate lattices that
can buckle in all three dimensions, a 3D TO should be implemented, taking into account the buckling response in all three dimensions.
Hence, there is a need to develop a 3D extension of the general 2D buckling-based TO code, which can then be applied to microstructural
design. As the microstructure/single cell of a lattice has a low volume fraction with mostly hollow domains, the number of elements
involved in active optimization will be less, thus not presenting any computational hurdle.
The contribution of this paper is two-fold:

1. Development of 3D buckling-based TO code. Though this code was developed and applied by researchers previously for optimizing
continuum structures, the code was not made available in any reported works. The code, originally taken from its 2D counterpart,
topBuck250 (Ferrari et al., 2021), solves TO problems involving a combination of compliance, volume, and buckling load factors
(BLFs). This work not only provides the 3D matrices of shape function, stiffness, stress stiffness, sensitivities, density filters, ob-
jective functions, and constraints, but also the 3D extension of vectorized matrix products developed in the 2D code for compact
representation of stress matrices and hence efficient implementation of the code. The computational cuts due to these imple-
mentations developed by Ferrari et al. are more substantial for the large 3D matrices of this work. The MATLAB code developed is
provided for educational purposes, which will alleviate the need for re-inventing the wheel for future research in buckling-based
TO.
2. Next, we use this 3D MATLAB code to design a single-cell microstructure of simple cubic (SC) and body-centric cubic (BCC) LSs
along with multi-cell 2 × 2 and 3 × 3 SC LS. Though TO on 3D LS was attempted in the literature (Wang & Sigmund, 2021), none
have attempted to do so without altering the relative density of the structure. In continuation of our previous studies in this area of
improving the buckling resistance of LS cells with slender members (struts/walls) by various methods without altering the relative
110 | Enhancing buckling resistance via three-dimensional topology optimization

Downloaded from https://academic.oup.com/jcde/article/12/8/107/8210400 by guest on 07 September 2025

Figure 1: TO-driven approach to design buckling-resistant SC LS.

density of the lattice (Viswanath et al., 2023, 2024, 2025). The current work discusses the methodology of using 3D TO. Similar
to each of the previous works, this paper also introduces another independent methodology not tested in the literature for 3D
lattices to increase buckling resistance but maintain the relative density. This can be easily implemented on any structure for
future research by editing only the geometry section in the accompanying code.

In the following sections, the methodology section details the procedure adopted for the single-cell design, the 3D TO section discusses
the mathematical formulations involved in the 3D TO, and the following section elaborates the matrices to be implemented in the 3D
formulation. The design section discusses the part of the code specific to the LSs investigated in this work related to building their
finite-element (FE) model. The final section of results concludes the paper with key topologies and their mode shapes.

2. Methodology
The proposed strategy for designing lattice unit cells with improved BLFs is shown in Figure 1. As observed from the previous work by the
authors (Viswanath et al., 2023), a preliminary investigation with a linear buckling analysis on any given LS to identify buckling-prone
members is essential to identify the weakest members that are susceptible to buckling under compressive loading. The first step is to
select a strut-based LS unit cell with a specific volume fraction and a traditional strut cross-section and subject it to linear buckling
analysis on the commercial FE program ABAQUS. Using the buckling mode shapes obtained from the eigenvalue buckling analysis, the
prominent buckling regions are identified to select those struts at greatest risk of buckling, typically the vertical and inclined columns.
These lattice columns are then modified using a 3D TO algorithm based on BLF maximization objective, maintaining the weak members
as an active optimization domain and the remaining members passive, while maintaining a constant relative density of the single cell.
2D MATLAB codes of TO available in literature are developed in this study for this purpose.
Prior to discussing the design procedure, the required theory and related codes of the 3D TO based on optimizing the BLFs are first
established and explained since this is the first paper to publish the 3D TO code. The main advantage of applying 3D TO to porous
LSs lies in the significantly lower computational effort required, as only a relatively small volume of material is optimized compared
to solid macrostructures. While the 3D TO framework developed in this work is general and applicable to all types of structures, its
implementation on solid macrostructures would be computationally intensive, largely depending on the FE mesh density used to model
the structure. In the case of LSs, as described in our methodology, the TO domain consists largely of void regions that represent the
porous areas of the lattice, as well as passive regions. These passive regions correspond either to the constrained parts of the structure,
which do not undergo deformation, or to boundary topologies provided to form the external framework of the final design. Both void and
passive elements are maintained at a constant density throughout the optimization process, with only the active elements participating
in the material redistribution. When a lattice with a low relative density (e.g. 15%) is considered, less than 15% of the total elements (due
to the inclusion of passive elements) actually contribute to the optimization. This results in a significantly faster computational process
compared to a fully solid 3D TO model. This computational efficiency is a key advantage leveraged in this study to specifically explore
buckling-based 3D optimized topologies of LSs.
Next, two sections describe this theory of 3D TO and the MATLAB code modifications on the existing 2D codes to incorporate the
same.
 Journal of Computational Design and Engineering, 2025, 12(8), 107–126 | 111

3. Buckling-Based 3D Topology Optimization

In this study, we employ 3D TO in order to design the buckling-resistant strut LSs. We follow the methodology outlined by Ferrari et al.
(2021, extending it to 3D elements), who developed a 250-line MATLAB program for optimizing linearized buckling. In this section, the
theory underpinning the buckling-based TO process on a 2D design domain extended to a 3D domain is briefly described. The formulation
of the buckling problem with 3D elements will be considered in the next section.
Consider a 3D structure domain  discretized into m unit size FEs with total of n degrees of freedom. To perform a linearized buckling
analysis of this structure, we specify a reference load vector F Rn acting on the domain and solve the linear system.

Ku = F (1)
nxn n
where K R is the linear, symmetric, positive-definite stiffness matrix, and u R is displacement vector. The linear buckling analysis
assumes an initial perturbed shape in each element that induces a stress-dependent stiffening effect in addition to the original linear
static stiffness. Let Kσ (u ) Rnxn be that symmetric but indefinite stress/geometric stiffness matrix, then the interaction between Kσ and
K in the structure is analogous to the mass and stiffness matrices in a dynamic modal analysis and hence involves the solution of an

Downloaded from https://academic.oup.com/jcde/article/12/8/107/8210400 by guest on 07 September 2025

eigenvalue problem of the form

[K + λKσ (u)] ϕ = 0, ϕ = 0 (2)

This gives the eigen pairs (λi , ϕi ), i B, set of buckling modes, consisting of BLFs λi and associated modes ϕi normalized such that
ϕ Tj Kσ ϕi = −δ ji . These BLFs when applied to the nominal static load will scale the stress stiffening terms. We are particularly interested
in the fundamental BLF λ1 associated with the critical load Fcr =λ1 F and given by

vT Kv
λ1 (ρ, u ) = min − (3)
n
v R ,v=0 vT Kσ (u) v
among all eigenvectors v. Since the optimization problem is intended to maximize λ1 , to obtain a minimization problem, we rewrite
Equation 2 as

[Kσ (u) + μK] ϕ = 0, ϕ = 0 (4)

obtaining λi = 1/μi , where λ1 associates to μ1 , maximum algebraic value of μi s. The value of μ1 involves the ‘max’ operator just like
‘min’ in Equation 3 and can be approximated by the Kreisselmeier and Steinhauser (KS) smooth aggregation function (Kreisselmeier &
Steinhauser, 1980) as
 m 
1 
JKS [μi ] (ρˆ ) = μ1 (ρˆ ) + ln e pAGG (μi (ρ)
ˆ −μ1 (ρ)
ˆ )
(5)
pAGG
i=1

with parameter pAgg [1, ∞ ). This gives an upper bound for μ1 and a lower bound for λ1 .
Let us now relate the material distribution of the structure to μ. The density-based TO performed on a regular FE grid considers design
variables as the vector of elemental material densities, ρ = {ρe }m
e=1 , introduced to generate a material distribution in the structure with a
value of 0 showing no material and 1 showing solid material. These element densities are first filtered into ρe , and then using projection
parameters η [0, 1] and β [1, ∞ ), converted to the pseudo-densities ρe by
m
he,i ρi tanh βη + tanh (β (ρe − η) )
ρe = i=1
m , ρe = (6)
i=1 he,i tanh βη + tanh (β (1 − η) )
where the minimum filter radius rmin > 0 is used for calculating he,i = max(0, rmin − dist( i ,  j ) ), ∀ i, j [1, m] for element domain
e (Bourdin, 2001). ρˆ is divided into three parts: active densities (ρˆ A ) obtained from Equation 4 for region A, passive densities with ρe =1
(ρˆ P1 ) for region P1 and passive densities with ρe =0 (ρˆ P0 ) for region P0 . These three regions guide the optimization process wherein the
material is re-distributed for every iteration only in A, while P1 always remain solid and P0 always a void. The solid isotropic material with
penalization interpolation parametrization for the element Young’s modulus in TO with separate penalizations pK , pG >1 for stiffness EK
and stress EG , respectively is

p p
EK (ρe ) = Emin + (Es − Emin ) ρˆe K EG (ρe ) = (Es ) ρˆe G (7)
−9
where Es = 1 GPa, the solid element Young’s modulus and Emin = 1e GPa, Young’s modulus of a void, introduced to avoid singularity
in element linear stiffness matrix. From the element matrices calculated from the above Young’s modulus, the global matrices Kσ and
K are assembled. For implementation of Kσ , we refer to de Borst et al. (2012), steps discussed in detail in next section. Before defining
the objective function and constraints of the optimization problem, the physical quantities of volume fraction V(ρ) ˆ , compliance c(ρ)
ˆ are
defined respectively as follows:

1 
m
V (ρˆ ) = ρˆe c (ρˆ ) = FT u (ρˆ ) (8)
m
e=1

The optimization problem for BLF maximization with compliance and volume constraints can be now stated as minimization of JKS
in Equation 5 as
112 | Enhancing buckling resistance via three-dimensional topology optimization

Downloaded from https://academic.oup.com/jcde/article/12/8/107/8210400 by guest on 07 September 2025

Figure 2: Flowchart of the changes in topBuck250 code (Ferrari et al., 2021). All changes incorporated for 3D extension are shown in dashed boxes with
corresponding block names of the code as provided in the MATLAB code in the Supporting Information. All changes are marked with a #3D tag in the
code for differentiation.

min :

ρ ∈ [0, 1]m J0KS ([μi ] (ρˆ ) )

(9)
such that : J1KS ([gv, gc] (ρˆ ) ) ≤ 0

V(ρ)
ˆ c(ρ)
ˆ
where KS function used to aggregate multiple constraint functions gv = V f − 1, gc = c f − 1 is given by J1KS . V f , c f are the maximum
allowed volume fraction and compliance, respectively. The sensitivities of the objective function and constraints calculated during the
TO are given as

∂c
∂ ρˆe
= −uT ∂∂K
ρˆe
uδeA
∂V
∂ ρˆe
= m−1 δeA
∂μi (ρ)ˆ
m pAGG (μi (ρ) ˆ )
ˆ −μ1 (ρ) (10)
∂JKS [μi ] e ∂ ρˆe

∂ ρˆe
(ρˆ ) = − i=1
m p
AGG (μi (ρ) ˆ )
ˆ −μ1 (ρ)
 i=1 e
∂μi T ∂Kσ
∂ ρˆe
= − ϕi ∂ ρˆe ϕi + μi ϕiT ∂∂Kρˆe i
ϕ − wTi ∂K σ
∂ ρˆe
u , wi = K−1 ϕiT (∇u Kσ ) ϕi

The Dirac delta function δeA =1 if e is in A, else 0 and wi is the adjoint vector. The above sensitivities with respect to pseudo-densities
will be used in the chain rule to calculate the sensitivities with respect to the actual densities of the elements. Adopting this problem
set-up, the next section will offer an in-depth description of the 3D formulation of matrices and changes required in existing 2D MATLAB
code to incorporate these.

4. Set-up of Buckling Eigenvalue Problem

4.1 Formulation of the 3D stiffness matrices
As discussed in the previous section, the eigenvalue problem is solved to obtain the BLFs, which are optimized to obtain buckling-resistant
structures. Here, we modified the 250-line MATLAB program for optimizing linearized buckling to cater to 3D optimization by changing
the 2D Q4 bilinear FE to 3D H8 trilinear element. Consequently, the modification in matrices from 2D to 3D and the related equations are
presented in detail here and consequent changes made in the MATLAB code are shown using a flowchart in Figure 2 and highlighted in
the MATLAB code given in the Supporting Information S1. Readers are encouraged to refer to Ferrari et al.(2021) for a detailed formulation
of the TO process and general structure of the 2D code along with Zhao et al. (2023) for the 3D formulations involved in a geometric non-
linear problem, which will be utilized here. To avoid repetitions, this paper details the modification in the 250-line MATLAB code only. The
code is called by the function topBuck3D (nelx, nely, nelz, penalK, rmin, ft, ftBC, eta, beta, maxit, ocPar, Lx, Ly, penalG, nEig, pAgg, prSel,
and fil). First, the 3D extension of code involves the initial discretization of elements and degree of freedom matrices in z-dimension by
introducing a ‘nelz’ parameter and using it to calculate all FE-related variables. The lines of code changed for this in the MATLAB code
are shown in block Ia in the Supporting Information S1.
 Journal of Computational Design and Engineering, 2025, 12(8), 107–126 | 113

Next, the setting up of stiffness matrices of the eight-node element is performed as in the 3D MATLAB code of Ferrari & Sigmund
(2020) shown in block Ib in the Supporting Information S1. The shape function array N(ξ , ζ , η) with logical coordinates (ξ , ζ , η) from the
isoparametric FE formulation is shown in Figure A1 and Table A1 in Appendix A. Let Xe = [xei , yei , zei ]8i=1 be the co-ordinates of an element,
then the (x, y, z)-gradient of N is computed as ∇(x,y,z) N = J−1 ∇(ξ , ζ , η) N where (ξ , ζ , η) – gradient and Jacobian are

⎡  ⎤
⎡ ⎤ ∂ξ Nk xk ∂ξ Nk yk
∂ξ N1 , . . . , ∂ξ Nk ⎢k 
k ⎥
⎢ ⎥ ⎢ ⎥
∇(ξ , ζ , η) N = ⎣ ∂ζ N1 , . . . , ∂ζ Nk ⎦ , J (ξ , ζ , η) = ⎢ ∂ N x ∂ N y ⎥ e
⎢ k ζ k k k ζ k k ⎥ = ∇(ξ , ζ , η) NX (11)
∂η N1 , . . . , ∂η Nk ⎣  ⎦
∂η Nk xk ∂η Nk yk
k k

In buckling analysis, in addition to the linear stiffness matrix K, the stress stiffness matrix Kσ (denoted by G in MATLAB code) also
contributes to the eigenvalue problem as shown in Equation 2. The stress stiffness matrix of a 3D H8 trilinear element is derived from
the stress components matrix σ0e constituting the 3D stress components, σxe , σye , σze , τxy
e
, τxz
e
, τyz
e
, at the element level. This stress matrix
is obtained from the constitutive relation,σ0 = DB0 u , with D as an elasticity matrix and B0 the linearized strain displacement pseudo
e e

Downloaded from https://academic.oup.com/jcde/article/12/8/107/8210400 by guest on 07 September 2025

differential operator. For 3D H8 trilinear element, B0 is obtained from

 
B0 = L ∇(x,y,z) N ⊗ I3 (12)
⎡ ⎤
1 0 0 0 0 0 0 0 0
⎢ ⎥
⎢0 0 0 0 1 0 0 0 0⎥
⎢ ⎥
⎢0 0 0 0 0 0 0 0 1⎥
where placement matrix L = ⎢ ⎥
⎢0 1 0 1 0 0 0 0 0⎥ and I3 is the identity matrix of order 3. Consequently, stresses are evaluated
⎢ ⎥
⎢ ⎥
⎣0 0 1 0 0 0 1 0 0⎦
0 0 0 0 0 1 0 1 0
at super convergent point of element centroids for all elements. The part of the MATLAB code to be changed for above formulation is
shown by block Ic in the Supporting Information S1.
The stress stiffness matrix Kσ is then efficiently calculated using vectorized matrix products. For each element, the formulation
described in Ferrari et al. (2021) is

Kσ e [ue ] :− ∫ BT1 T [σ0e ] B1 d0 (13)

0

where B1 discretizes the deformation gradient and for an H8 trilinear element (de Borst René et al., 2012) given by
⎡ ⎤
∂x N1 0 0 ∂x N2 0 0 . . ∂x N8 0 0
⎢∂y N1 0 0 ∂y N2 0 0 . . ∂y N8 0 0 ⎥
⎢ ⎥
⎢∂ N 0 0 ∂z N2 0 0 . . ∂z N8 0 0 ⎥
⎢ z 1 ⎥
⎢ 0 ∂x N1 ∂x N2 . . ∂x N8 0 ⎥
⎢ 0 0 0 0 ⎥
⎢ ⎥
⎢ 0 ∂y N1 0 0 ∂y N2 0 . . 0 ∂y N8 0 ⎥
⎢ ⎥
⎢ 0 ∂z N1 0 0 ∂z N2 0 . . 0 ∂z N8 0 ⎥
⎢ ⎥
⎢ 0 0 ∂x N1 0 0 ∂x N2 . . 0 0 ∂x N8 ⎥
⎢ ⎥
⎣ 0 0 ∂y N1 0 0 ∂y N2 . . 0 0 ∂y N8 ⎦
0 0 ∂z N1 0 0 ∂z N2 . . 0 0 ∂z N8 9×24

written in operator notation similar to Equation 10 as

⎡ ⎤
∇(x,y,z) N ⊗ [1, 0, 0]T
⎢ ⎥
B1 = ⎣∇(x,y,z) N ⊗ [0, 1, 0]T ⎦ (14)
T
∇(x,y,z) N ⊗ [0, 0, 1]

The stress components matrix is

⎡ ⎤
σx τxy τxz 0 0 0 0 0 0
⎢ ⎥
⎢τyx σy τyz 0 0 0 0 0 0⎥
⎢ ⎥
⎢τzx τzy σz 0 0 0 0 0 0⎥
⎢ ⎥
⎢0 0 0 σx τxy τxz 0 0 0⎥
⎢ ⎥
⎢ ⎥
T [σ0e ] = I3 ⊗ σ0e = ⎢ 0 0 0 τyx σy τyz 0 0 0⎥ (15)
⎢ ⎥
⎢0 0 0 τzx τzy σz 0 0 0⎥
⎢ ⎥
⎢0 0 0 0 0 0 σx τxy τxz ⎥
⎢ ⎥
⎢ ⎥
⎣0 0 0 0 0 0 τyx σy τyz ⎦
0 0 0 0 0 0 τzx τzy σz 9×9
114 | Enhancing buckling resistance via three-dimensional topology optimization

Table 1: Mapping Mcol1 : {1 : 36} → (i, j ) for the 36 unique coefficients of Z and lower symmetric part of Kσ e . The consecutive columns of
Kσ e also shown as Mcol2 and Mcol3.

1 … 8 9 … 15 16 … 21 22 … 31 … 34 …
36

Mcol1 (1,1) . (22,1) (4,4). (22,4) (7,7).(22,7) (10,10).(16,16). (19,19).(22,22)

Mcol2 (2,2). (23,2) (5,5) . (23,5) (8,8).(23,8) (11,11).(17,17). (20,20).(23,23)
Mcol3 (3,3). (22,1) (4,4) . (22,4) (9,9).(24,9) (12,12).(18,18). (21,21).(24,24)

Expanding integrand of Equation 13, we obtain the following structure (de Borst René et al., 2012) in which we identify 36 independent
coefficients out of 576, highlighted in red.

⎡ ⎤
z11 0 0 z12 0 0 . . z18 0 0

Downloaded from https://academic.oup.com/jcde/article/12/8/107/8210400 by guest on 07 September 2025

⎢0 z11 0 0 z12 0 . . 0 z18 0⎥
⎢ ⎥
⎢ ⎥
⎢0 0 z11 0 0 z12 . . 0 0 z18 ⎥
⎢ ⎥
⎢z21 0 0 z22 0 0 . . z28 0 0⎥
⎢ ⎥
⎢0 . . 0⎥
⎢ z21 0 0 z22 0 0 z28 ⎥
⎢ ⎥
⎢0 0 z21 0 0 z22 . . 0 0 z28 ⎥
⎢ ⎥
Kσ e {H8} = ⎢
⎢z31 0 0 z32 0 0 . . z38 0 0⎥ ⎥ (16)
⎢0 . . 0⎥
⎢ z31 0 0 z32 0 0 z38 ⎥
⎢ ⎥
⎢0 0 z31 0 0 z32 . . 0 0 z38 ⎥
⎢ ⎥
⎢ . . . . . . . . . . . ⎥
⎢ ⎥
⎢ ⎥
⎢z81 0 0 z82 0 0 . . z88 0 0⎥
⎢ ⎥
⎣0 z81 0 0 z82 0 . . 0 z88 0⎦
0 0 z81 0 0 z82 . . 0 0 z88 24×24

where zik = σx ak ai + σy bk bi + σz ck ci + τxy (bk ai + ak bi ) + τxz (ak ci + ck ai ) + τyz (bk ci + ck bi ) is a linear combination of stress components with
ai = ∂x Ni , bi = ∂y Ni and ci = ∂z Ni . A compact representation of these independent coefficients is given by array Z Rm×36 , obtained
through the Hadamard matrix product Z = BT S.
S = [σxe , σye , σze , τxy
e
, τxz
e
, τyz
e
]e =1: m collects the stress components of all elements and array B[6 x 36] is

⎡ ⎤i=1:36
al(i,1) al(i,2)
⎢ ⎥
⎢ bl(i,1) bl(i,2) ⎥
⎢ ⎥
⎢ c c ⎥
B = ⎢
⎢b
l( i, 1 ) l( i, 2 ) ⎥
⎥ (17)
⎢ l(i,2) al(i,1) + al(i,2) bl(i,1) ⎥
⎢ ⎥
⎣ al(i,2) cl(i,1) + cl(i,2) al(i,1) ⎦
cl(i,2) bl(i,1) + bl(i,2) cl(i,1)

where l N36×2 is the set of indices mapping each Zei (i = 1, …, 36) to the corresponding zl (i,1)l (i,2) within the pattern of non-zero elements
of (), i.e.
 
1 2 . 7 8 2 3 . 8 3 4 . 8 . . 7 8 8 row index
l=
1 1 . 1 1 2 2 . 2 3 3 . 3 . . 7 7 8 col index

The columns of Z collect the coefficients highlighted in the reduced pattern shown below, sorted column-wise.

⎡ ⎤
z11 z12 z13 z14 z15 z16 z17 z18
⎢z z22 z23 z24 z25 z26 z27 z28 ⎥
⎢ 21 ⎥
⎢ ⎥
⎢z31 z32 z33 z34 z35 z36 z37 z38 ⎥
⎢ ⎥
⎢z41 z41 z43 z44 z45 z46 z47 z48 ⎥
=Z=⎢ ⎥
unique
Kσ e {H8} ⎢z (18)
⎢ 51 z52 z51 z54 z55 z56 z57 z58 ⎥
⎥
⎢ ⎥
⎢z61 z62 z63 z64 z65 z66 z67 z68 ⎥
⎢ ⎥
⎣z71 z72 z73 z74 z75 z76 z77 z78 ⎦
z81 z82 z83 z84 z85 z86 z87 z88 8×8

Hence, the global stress stiffness matrix Kσ can be assembled, as discussed in [19], by considering a mapping Mcol1 between the
columns of Z and the entries of Kσ e indicated in red in (16) as given in Table 1. The two consecutive columns of Kσ e are obtained by
adding 1 and 2, respectively, to indices of Mcol1 as shown in table. Considering the column-wise sorting of lower symmetric part of (16),
Mcol1 corresponds to 36 entries of {1, 4, 7, 10, 13, 16, 19, 22, 70, 73, 76, 79, 82, 85, 88, 130, 133, 136, 139, 142, 145, 181, 184, 187, 190, 193, 223,
226, 229, 232, 256, 259, 262, 280, 283, 295} of Kσ e . Finally, the lower symmetric part of the geometric matrix is assembled by using fsparse
(Engblom & Lukarski, 2016).
The block of code for the above formulation of geometric matrix is marked as if with indices vector of Mcol1 computed as ‘IkG’ in block
Ic of the Supporting Information S1.
 Journal of Computational Design and Engineering, 2025, 12(8), 107–126 | 115

Downloaded from https://academic.oup.com/jcde/article/12/8/107/8210400 by guest on 07 September 2025

Figure 3: Geometry setup of (A) SC and (B) BCC LS optimization showing cell size L, column size s, beam size s1 , active region A, passive region P1
densities with ρe = 1 (solid), passive region P0 with ρe = 0 (void).

The sensitivity expressions of

   2  
ϕT Kσ e0 ϕ = ϕ12 + ϕ22 + ϕ32 z11 + · · · + ϕ22 2
+ ϕ23 2
+ ϕ24 z88

+ [2 (ϕ1 ϕ4 + ϕ2 ϕ5 + ϕ3 ϕ6 ) z21 + · · · + 2 (ϕ1 ϕ22 + ϕ2 ϕ23 + ϕ3 ϕ24 ) z81 ]

+ [2 (ϕ4 ϕ7 + ϕ5 ϕ8 + ϕ6 ϕ9 ) z32 + · · · + 2 (ϕ4 ϕ22 + ϕ5 ϕ23 + ϕ6 ϕ24 ) z82 ]

+ [2 (ϕ7 ϕ10 + ϕ8 ϕ11 + ϕ9 ϕ12 ) z43 + · · · + 2 (ϕ7 ϕ22 + ϕ8 ϕ23 + ϕ9 ϕ24 ) z83 ]

+ · · · + [2 (ϕ16 ϕ19 + ϕ17 ϕ20 + ϕ18 ϕ21 ) z76 + 2 (ϕ16 ϕ22 + ϕ17 ϕ23 + ϕ18 ϕ24 ) z86 ]

+ [2 (ϕ19 ϕ22 + ϕ20 ϕ23 + ϕ21 ϕ24 ) z87 ] (19)

]can be obtained with the Hadamard product Z˜ p where p R36×m collects all coefficients of φ as shown in Table A2 given in Appendix
A. p differs from 2D formulation that it has three terms corresponding to Mcol1 , Mcol2 and Mcol3 . To have consistency with (17), the
elements of columns {2,3,4,5,6,7,8, 10,11,12,13,14,15, 17,18,19,20,21, 23,24,25,26, 28,29,30, 32,33, 35} of array Z must be doubled, obtaining
Z˜ . This part of code is given in block Ic for computation of ‘dZdu’ to be used for sensitivity computation in block Ig and Ih of the
Supporting Information S1.
Apart from the changes in matrices and equations, the MATLAB code is changed for different structures based on geometry, loading
conditions, and its boundary conditions marked by code block II in the Supporting Information S1 and Figure 2. This part of the code will
change for each LS example as presented in th Supporting Information S2 and will be discussed separately in subsequent subsections
in the next section.

5. Examples of Lattice Design

The next part of this study involved designing LS using the TO software tool developed in the sections above. As strut-based lattices are
being investigated, both vertical strut and inclined strut lattices are studied to investigate the variations in buckling of the two. Once the
methodology with changes to the MATLAB code for each geometry investigated is established, to show the scalability of the methodology
to multi-cell lattices, two examples of vertical strut multi-cell lattices of 2 × 2 and 3 × 3 are also presented.

5.1 SC structure
The first primary structure chosen for this study was a 3D SC single cell having a prescribed desired relative density (Figure 3A).
A cell size L = 20 mm, with beam (s1 ) of 2 mm and column thickness (s) of 4 mm to achieve a specific relative density was
adopted. The physical length Lx and Ly in the code are kept as unity, as the aspect ratio of a square. In the previous work by
the authors (Viswanath et al., 2023), a linear buckling analysis of the SC structure, performed as a preliminary investigation of
its buckling modes, indicates that the columns are most susceptible to buckling under compressive loading. Consequently, the
columns of the SC were modified using TO to achieve an optimal material distribution while preserving the relative density of
the SC, as elaborated in the next section. The BCs on FE model were not periodic, as we consider the SC to be an independent
structure and instead applied as a compressive load of 1 N on top surface with fixity on bottom surface schematically shown in
Figure 3A.

5.1.1 Optimization setup

As described in the associated MATLAB code description in Section 4.1, various TO parameters must first be defined before modelling
the structure. The main parameters affecting the buckling optimization process are the filter radius and optimality criteria parameters,
along with penalization exponents, projection parameters, aggregation parameter, and all their continuation schemes. These parameters
are substituted while calling the code as described in Section 4.1 as
116 | Enhancing buckling resistance via three-dimensional topology optimization

(20,20,20,3,1.5,2,‘N’,0.5,2 500, [0.1,0.7,1.2],1,1,3,12 200,{[‘B’,‘C’,‘V’], [2.5,0.15]},‘ig.mat’).

and the continuation schemes at the start of the code as

penalCntK = { 25, 1, 25, 0.25}; % continuation scheme on penalK

penalCntG = {25,1,25,0.25}; % “ ” on G-penal
betaCnt = { 300, 18, 50, 2}; % “ ” on beta
pAggCnt = { 2e5,1,25,2}; % “ ” on the KS aggregation factor

The final argument ig.mat refers to the initial design supplied to the TO code in the first iteration of buckling-based TO. Here, a
minimum volume design is obtained by first running the code with [‘V,’C’] for prsel argument without any buckling computations. The

Downloaded from https://academic.oup.com/jcde/article/12/8/107/8210400 by guest on 07 September 2025

densities obtained in this run are stored in ig.mat in the above buckling-based TO. The function call of of the code for the volume
minimization TO run performed apriori is
(20,20,20,3,1.5,2,‘N’,0.5,2 500, [0.1,0.7,1.2],1,1, [], [], [],{[‘V’,‘C’],2.5})
and the densities obtained are saved with the command ‘save ig.mat xPhys’.
Next comes the domain partitioning by setting active and passive regions. The domain of active optimization (A) will be concentrated
on columns based on information obtained from the preliminary investigation of the reference LS, and the rest of the structure is treated
under passive regions (P1 ). So as to provide an outer boundary for the topology that will evolve through TO, three edge lines are also
provided for the columns in P1, as shown in Figure 3. The inner edges of the columns remain active and hence change during TO. This
method of providing box boundaries is inspired by previous work of authors in which the edge-bounded 2D TO yielded the best buckling-
resistant design (Viswanath et al., 2024). The void space in the centre of the SC LS is also assigned (P0 ).
The code changes in the Supporting Information S2 show this domain set-up and the application of loading and BCs on the LS as
shown in Figure 3A.
The model is then analysed, and the optimized designs are obtained. The BLFs and buckling modes obtained from the MATLAB code
are investigated.

5.2 BCC structure

Another LS investigated in this study was the BCC strut lattice with circular columns. The lattice was first voxelized using an algorithm
developed by Dong et al. (2018) for ease of application of the TO algorithm to lattices with curved surfaces. When voxelized with cubes,
the inclined circular columns look like square columns, and by increasing the FE mesh (voxels), a smoother topology can be obtained
but at the cost of computational time consumed by the TO. The TO domain will show a rough surface for the circular strut with coarse
mesh which gets smoother as we refine the mesh. This paper considers only 20 and 40 mm unit cell size. Hence, slight roughness is
evident on walls of circular struts in figures. This lattice has inclined members as high buckling regions, modelled as active regions for
optimization. Passive elements include elements upon which loads and BCs are applied. No passive edges are provided for the inclined
columns as the cross-section is circular.

5.2.1 Optimization setup

Similar to the SC LS setup in Section 5.1.1, the BCC lattice with the active-passive-void regions and the BCs were established. First, the
voxelization of the lattice was carried out with the line of code obtained from Dong et al. (2018), that is,
[voxel, Density] = GenerateVoxel(40,‘bcc.txt’,0.1);
and the voxel matrix is stored in a mat file named ‘initial.mat’. The function call with all the TO parameters would then be
(40,40,40,3,1.5,2,‘N’,0.8,4 500, [0.01,0.1,2.2],1,1,3,12 200,{[‘B’,‘C’,‘V’], [2.5,0.10]},‘initial.mat’).
The model setup follows code shown in the Supporting Information S2.

5.3 Multi-cell lattice structures

Inspired by the single-cell studies of SC and BCC lattices, the work was extended to multiple-cell LS. The SC LS of 20 mm × 20 mm ×
20 mm was chosen to create the 2 × 2 × 1 matrix of total cell size 40 mm × 40 mm × 20 mm consisting of four SC cells and 3 × 3 × 1
matrix of total cell size 60 mm × 60 mm × 20 mm consisting of nine SC cells. The column and beam sizes and the unit cell size remain
the same as those considered in Section 5.1, including the loading, boundary conditions, the active-passive-void regions of each cell and
the TO parameters.

5.3.1 Optimization setup

The minimum volume TO design run was performed and the obtained topology used as the initial design with the following function
calls:
(40,40,20,3,1.5,2,‘N’,0.5,2 500, [0.1,0.7,1.2],1,1, [], [], [],{[‘V’,‘C’],2.5})
(40,40,20,3,1.5,2,‘N’,0.5,2 500, [0.1,0.1,1.2],1,1,3,12 200,{[‘B’,‘C’,‘V’], [2.5,0.15]},‘ig.mat’)
for 2 × 2 × 1 lattice and
(60,60,20,3,1.5,2,‘N’,0.5,2 500, [0.1,0.7,1.2],1,1, [], [], [],{[‘V’,‘C’],2.5})
(60,60,20,3,1.5,2,‘N’,0.5,2 500, [0.1,0.1,1.2],1,1,3,12 200,{[‘B’,‘C’,‘V’], [2.5,0.15]},‘ig.mat’)
for 3 × 3 × 1 lattice. The model setup code is given in the Supporting Information S2.
 Journal of Computational Design and Engineering, 2025, 12(8), 107–126 | 117

6. Numerical Analysis
Further to obtaining the optimized topology, a numerical linear buckling analysis on one of the investigated LSs was performed for
evaluating the first eigenvalue of the optimized topology relative to a reference structure based on the same relative density to determine
the improvement in buckling resistance. To achieve this, the optimized topology was then exported as an STL file from MATLAB and
processed in ABAQUS for a further linear buckling analysis validation to compare its first buckling mode with the reference structure.
Post-processing of the CAD files in SOLIDWORKS software was performed for the ease of creating refined meshes with quality checks
since most optimized topologies have surfaces that may prove difficult to mesh using analysis software. The meshed model is exported
as a part file to ABAQUS software for linear buckling analysis validation.
The linear buckling modes obtained from buckling analysis in ABAQUS are compared with those of reference LS previously analysed
in ABAQUS, maintaining the same volume fraction for the columns in both designs employing 10-node tetrahedral elements (C3D10)
for elements. The material properties were taken as Young’s modulus of 1 GPa and a Poisson’s ratio of 0.3. The same size of FE mesh
and BCs application that the authors have developed in their previous works used for validating experimental tested lattices is used for
numerical modelling for comparison for a future experimental validation. The post-buckling numerical analysis, required to calculate

Downloaded from https://academic.oup.com/jcde/article/12/8/107/8210400 by guest on 07 September 2025

the critical buckling strength of the structure, and experimental validation are works in progress and will contribute to the subject of
the next research and are out of the scope of this work.

7. Results and Discussions

A unit cell size of 20 mm × 20 mm × 20 mm was chosen for the SC structure. The beam width (s1 ) was reduced to 2 mm for a longer
column length. The column thickness (s) was maintained at 4 mm. The three outer edges of the columns were maintained as passive
regions, which would have a solid density throughout optimization and maintain the outer boundary of the topology after optimization.
In contrast, the rest of the element densities in the column would change. The structure (as shown in Figure 4A) is first topology optimized
for minimum volume design subject to a compliance restraint without considering any buckling-based TO parameters as explained in
Section 5.1.1. This is performed to obtain an initial design optimized in terms of weight and compliance, stored in the ig.mat file and
used a base design for buckling-based optimization. Starting the buckling-based TO with a minimum volume design also expedites this
time-consuming TO, since it involves an added solution of the eigenvalue problem. Figure 4A shows the minimum volume design and
buckling optimized topology design obtained for the columns along with the BLF plot showing progression of the first four BLFs through
the optimization. The minimum volume design in the left shows buckling prone struts and unconnected regions obtaining volume
fraction of 18%. Hence, for the buckling TO, a volume fraction constraint of 18% was chosen so as to obtain a well-connected topology.
The BLF plot in centre highlights an increasing trend during the optimization process, showing an increased buckling capacity of the
final design from the initial minimum volume design. The arrows indicate the BLFs corresponding to the two designs. Next, the first
three buckling modes obtained from the TO code for the final design are plotted in Figure 4B. The computational time taken for each
iteration of buckling-based TO is approximately 15 s and hence a TO with maximum of 300 iterations would be completed in 1.25 h and
did not pose any computational burden.
To validate the linear buckling analysis, the CAD is now exported from the STL format obtained in MATLAB to Solidworks software
for conversion to a part file, which can be imported and analysed using ABAQUS. The first mode obtained from the numerical analysis is
shown in Figure 4C compared to a base structure with the same volume fraction. The eigenvalue for the optimized SC was much higher
than that of base structure, validating that TO design performs better in buckling resistance. The eigenvalues obtained from MATLAB
and ABAQUS cannot be compared as both use different FE types and discretization. Hence, we obtain the topology using MATLAB while
its verification is realized through ABAQUS.
Next, different relative densities for the 20 mm cell size were tested, and the obtained BLFs and final topologies are shown in Figure 5.
For lower (12%) and higher (30%) volume fractions, the passive regions of the columns were eliminated as they rendered very thin
members, which lowered the BLF. This may be the reason why the 12% lattice showed a higher BLF than the 18% lattice which had
passive parts in the column. More parametric studies include further testing various cell sizes for the SC structure with a relative
density of 18%. The obtained topologies are shown in Figure 6 for cell sizes 40 and 60 mm, which both showed BLFs with increasing
trend from their respective initial minimum volume designs. The volume fraction constraint for the buckling TO was considered as the
minimum volume fraction obtained in the minimum volume design. Therefore, it is noteworthy that both the designs, shown in left and
right of the BLF plots, are of the same volume fraction, but the buckling design is superior due to its higher buckling resistance. For these
higher cell sizes, the computational time taken an iteration for 40 mm was 225 s and 60 mm was 1200 s, limiting the study of higher cell
sizes.
Similarly, for the BCC structure of 40 mm × 40 mm × 40 mm cell size, the minimum volume design, the optimized topology and the
trend in the BLF are shown in Figure 7A for a volume fraction of 10%. The high stress concentration is normally observed in the centre of
the BCC and the optimized topology where material has been removed from the centre to decrease stress concentrations and improve
buckling resistance. The computational time taken for each iteration was only 10 s. The buckling modes are shown in Figure 7B. The
optimized topology indicated an increasing BLF from the initial topology of the optimization. These results were numerically validated
by performing a linear buckling analysis in ABAQUS and comparing the first eigenvalues against the traditional lattice of circular BCC
cell similar to SC LS analysis. As shown in Figure 7C, the first eigenvalue is slightly higher for the optimized topology than the reference
BCC with uniform circular cross-section.
However, all the topologies obtained in this study have to be additively manufactured for experimental validation. These tasks are
out of the scope of this work but are definitely the future scope of this research.
Similar to SC LS, parametric studies were conducted for BCC lattice by testing two different relative densities for two cell sizes of the
BCC lattice. The topologies are shown in Table 2 and the BLFs obtained are shown in Figure 8. A strut diameter of 5 and 8 mm was used
118 | Enhancing buckling resistance via three-dimensional topology optimization

Downloaded from https://academic.oup.com/jcde/article/12/8/107/8210400 by guest on 07 September 2025

Figure 4: (A) Minimum volume design and buckling reinforced design for SC LS of relative density 18% with the plot of first four BLFs. The
buckling-based TO design has a higher BLF compared to minimum volume TO design, (B) first three buckling modes for TO design obtained from
MATLAB TO algorithm, and (C) first buckling mode for base square column SC (left) and TO SC LS (right) from numerical simulations on ABAQUS.
Eigenvalue of modified SC structure is much greater than base structure of same relative density.

for 40 mm cell size and 3 and 6 mm for 20 mm cell size to obtain 10% and 25%, respectively. The TO parameter mainly affected by cell
size variations is the filter radius, which was reduced for lower cell sizes of 20 mm to rmin = 2, while it was 3 for 40 mm size and increased
to 4 for 60 mm cell size of SC lattice. The BLFs show vast variation among relative densities, unlike the SC lattice. This may be due to
the inclined member buckling, which is different from the straight column buckling observed in SC LS. The lower cell size of 20 mm also
developed completely hollow centres in TO LS (top views in Table 2), whereas the 40 mm cell size only had partial hollow regions in the
centre in both relative densities. However, the BLFs are approximately similar for any specific relative density of both cell sizes. The BLFs
obtained from the TO do not indicate comparison of buckling resistance among cell sizes but indicates how TO generates topologies
comparable to their base design.
 Journal of Computational Design and Engineering, 2025, 12(8), 107–126 | 119

Downloaded from https://academic.oup.com/jcde/article/12/8/107/8210400 by guest on 07 September 2025

Figure 5: Variation of BLF with relative densities in SC lattice.

Figure 6: TO for different cell sizes for relative density of 15% with minimum volume design (left), BLF (centre), and optimized design (right) of (A)
40 × 40 × 40 mm and (B) 60 × 60 × 60 mm.

The geometrical set-up of the multi-cell lattice optimized for buckling is shown in Figure 9A. The top beams carry the compressive
loads, and the bottom beams are kept fixed, similar to the single cell. The minimum volume design is obtained with angle-shaped beams,
which are improved for BLF with the buckling-based TO (Figure 9C). The volume fraction constraint is 15%. The BLF shows an increasing
trend, and the first three modes of the buckling-resistant structure are shown in Figure 9B. The topology of each cell is similar to the
single cell obtained with the same TO parameters. The topology is symmetrical about the centre due to the symmetrical loading and
symmetry of geometry. In multi-LSs, the buckling is global as well as local. The compressive loads are transmitted from top to bottom
cells, along with the weight of the top cells. Hence, they buckle differently and possess different buckling-based topologies, with the top
cells having wider holes than the lower cells.
Numerical simulation on ABAQUS (Figure 9D) indicates the first mode of TO SC lattice to be much greater than the base structure
with square columns of 3 mm width, both possessing the same relative density. This result once again reinforces the findings of this
work that 3D TO-based buckling helps in obtaining buckling resistant single and multi-cell LS. A 3 × 3 × 1 structure is also studied to
reinforce the effectiveness of the algorithm in multi-cell structures. The active-passive domain of optimization, the loads and boundary
conditions are similar to the 2 × 2 × 1 multi-cell structure. Figure 10A shows the optimized topology obtained, the increasing trend in
BLF plot shown in Figure 10B and the first buckling mode shown in Figure 10C.
120 | Enhancing buckling resistance via three-dimensional topology optimization

Downloaded from https://academic.oup.com/jcde/article/12/8/107/8210400 by guest on 07 September 2025

Figure 7: (A) Initial design and Buckling reinforced TO design of BCC lattice for 10% relative density and BLF plot from initial to final design, (B) the first
three buckling modes obtained from TO algorithm, and (C) first buckling mode for base circular column BCC (left) and optimized BCC LS (right) from
numerical simulations on ABAQUS. Eigenvalue of modified BCC structure is greater than base structure of same relative density.

The next lattice of application would be a periodic LS used in real-world structures, considering periodic boundary conditions (PBC)
to a unit cell of the lattice (microstructure). This extension is straightforward, which only requires incorporating PBC in the MATLAB
code instead of the compression boundary conditions in the single-cell study of this work. However, periodic lattice studies should
ideally include multi-cell periodic structures with the study of cell interactions, local buckling of unit cells, and the impact of cell
size and number of cell arrays. A different active-passive-domain combination of optimization is also to be ascertained by studying
 Journal of Computational Design and Engineering, 2025, 12(8), 107–126 | 121

Table 2: Study of buckling reinforced TO designs varying cell sizes and relative densities.

Cell
Relative density Front-side view Top view
Size
20 x 20 x 20

10

25

Downloaded from https://academic.oup.com/jcde/article/12/8/107/8210400 by guest on 07 September 2025

40 x 40 x 40

25

Figure 8: BLF for buckling reinforced TO designs of varying cell sizes and relative densities.

the extent of buckling in different members in multiple cells of the modes shapes of the structure under investigation and selectively
assigning members with lower buckling (closer to supports) as active domains and those with higher buckling (centre) as passive-active
combinations, to arrive at the best buckling resistant design. Hence, a periodic cell structure study is out of the scope of this research
and will be pursued by the authors in future research.
The scope of research in this work is first to develop the 3D MATLAB tool (not available in literature) and next to use it on
few representative LSs, such as simple single-cell lattices with vertical and inclined columns and a multi-cell lattice for under-
standing multi-cell buckling behaviour, to improve their buckling response without changing the relative density. The optimized
structures obtained with the mathematically formulated code are then validated with numerical studies by modelling the opti-
mized geometry with the same size of FEM element mesh and BCs application that the authors have developed in their previ-
ous works used for validating experimental tested lattices. While manufacturing and experimental validation of these optimized
designs are essential next steps, they fall beyond the current scope of this paper. Such efforts will require careful consideration
of design for additive manufacturing constraints and will be addressed in a dedicated future publication. Accordingly, the present
work focuses solely on numerical validation of the improved buckling-optimized designs generated by the developed MATLAB
code.

8. Conclusion
In this study, we have introduced a 3D TO-based methodology for designing buckling-resistant LSs. The MATLAB code for buck-
ling objective/constraint based on a 3D TO is provided for the first time, inspired by the 2D code available in literature. It is then
tested on a single cell of an SC lattice with square cross-section struts and a BCC lattice with circular struts and on a multi-
cell 2 × 2 and 3 × 3 SC LS. By analysing their initial buckling modes, the buckling-prone members are identified and actively
optimized using the newly developed 3D TO algorithm to maximize the BLFs while constraining the volume fraction and stiff-
122 | Enhancing buckling resistance via three-dimensional topology optimization

Downloaded from https://academic.oup.com/jcde/article/12/8/107/8210400 by guest on 07 September 2025

Figure 9: (A) Geometry setup multi-cell SC LS optimization, (B) first three modes obtained from TO, (C) minimum volume design, buckling reinforced
design for 15% volume fraction with BLF plot, and (D) first buckling mode for base square SC LS (left) and optimized SC LS (right) from numerical
simulations on ABAQUS. Eigenvalue of TO-modified SC structure is greater than base structure of same relative density 15%.

ness of the structure. The results demonstrate that the topology-optimized designs have higher BLFs than the reference struc-
ture, which must be numerically validated with post-buckling studies and experimental validation by 3D printing. Numerical val-
idations via linear buckling analysis is also performed in this study and the optimized designs are found to have a higher first
eigenvalue than their reference structure. Parametric studies with varying cell sizes and relative densities were also conducted
to study the effect of TO parameters on these variations. The major contributions of this paper are: (1) the developed MATLAB
code, which is of educational value to any research requiring a 3D approach to buckling-based TO, and (2) the methodology, which
can be easily extended to enhance the buckling resistance of any LS while maintaining the relative density of the corresponding
reference LS.
 Journal of Computational Design and Engineering, 2025, 12(8), 107–126 | 123

Downloaded from https://academic.oup.com/jcde/article/12/8/107/8210400 by guest on 07 September 2025

Figure 10: (A) Buckling reinforced design for 15% volume fraction, (B) BLF plot, and (C) first buckling mod.

Supplementary Data
Supplementary data are available at JCDENGN Journal online.

Conflicts of Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to
influence the work reported in this paper.

Author Contributions
Asha Viswanath: Data curation, Conceptualization, Methodology, Investigation, Visualization, Software, Validation, Writing—original
draft. Wesley Cantwell: Writing—review & editing, Funding Acquisition. Kamran Khan: Supervision, Conceptualization, Methodology,
Investigation, Writing—review & editing.

Funding
This publication is based on work supported by the Khalifa University of Science and Technology under award no. CIRA-2021–109 /
8474000423.

Data Availability
The TO code used in this work is available in the Supporting Information. The authors reserve all rights to the code. The code may be
distributed and used for academic and educational purposes. The authors do not guarantee that the code is free from errors, and they
shall not be liable in any event caused by the use of the code.

Acknowledgments
The authors thank Prof. Ole Sigmund and Dr Federico Ferrari of the Technical University of Denmark for the helpful discussions that
helped in the preparation of the manuscript.

References
Bendsoe, M. P., & Sigmund, O. (2013). Topology Optimization: Theory, Methods, and Applications. (5th ed.). Springer Science & Business Media.
Bluhm, G. L., Christensen, K., Poulios, K., Sigmund, O., & Wang, F. (2022). Experimental verification of a novel hierarchical lattice material with
superior buckling strength. APL Materials, 10, 090701. https://doi.org/10.1063/5.0101390.
Bourdin, B. (2001). Filters in topology optimization. International Journal for Numerical Methods in Engineering, 50, 2143–2158. https://doi.org/10.1
002/nme.116.
Bruyneel, M., Colson, B., & Remouchamps, A. (2008). Discussion on some convergence problems in buckling optimisation. Structural and Multi-
disciplinary Optimization, 35, 181–186. https://doi.org/10.1007/s00158- 007- 0129- z.
Christensen, C. F., Engqvist, J., Wang, F., Sigmund, O., & Wallin, M. (2025). Extremal structures with embedded pre-failure indicators. Proceedings
of the National Academy of Sciences, 122, e2412285122. https://doi.org/10.1073/pnas.2412285122.
Christensen, C. F., Wang, F., & Sigmund, O. (2023). Topology optimization of multiscale structures considering local and global buckling response.
Computer Methods in Applied Mechanics and Engineering, 408, 115969. https://doi.org/10.1016/j.cma.2023.115969.
124 | Enhancing buckling resistance via three-dimensional topology optimization

de Borst René, C. M., Remmers, J., & Verhoosel, C. (2012). Nonlinear Finite Element Analysis of Solids and Structures(2nd edn), John Wiley and Sons Ltd.
https://www.wiley.com/en-us/Nonlinear+Finite+Element+Analysis+of+Solids+and+Structures%2C+2nd+Edition-p-9781118376010.
Dong, G., Tang, Y., & Zhao, Y. F. (2018). A 149 line homogenization code for three-dimensional cellular materials written in matlab. Journal of
Engineering Materials and Technology, 141, 011005. https://doi.org/10.1115/1.4040555.
Engblom, S., & Lukarski, D. (2016). Fast Matlab compatible sparse assembly on multicore computers. Parallel Computing, 56, 1–17. https://doi.or
g/10.1016/j.parco.2016.04.001.
Ferrari, F., & Sigmund, O. (2020). A new generation 99 line Matlab code for compliance topology optimization and its extension to 3D. Structural
and Multidisciplinary Optimization, 62, 2211–2228. https://doi.org/10.1007/s00158- 020- 02629- w.
Ferrari, F., Sigmund, O., & Guest, J. K. (2021). Topology optimization with linearized buckling criteria in 250 lines of Matlab. Structural and
Multidisciplinary Optimization, 63, 3045–3066. https://doi.org/10.1007/s00158- 021- 02854- x.
Gao, J., Wu, X., Xiao, M., Nguyen, V. P., Gao, L., & Rabczuk, T. (2023). Multi-patch isogeometric topology optimization for cellular structures with
flexible designs using Nitsche’s method. Computer Methods in Applied Mechanics and Engineering, 410, 36. https://doi.org/10.1016/j.cma.2023
.116036.
Khan, N., & Riccio, A. (2024). A systematic review of design for additive manufacturing of aerospace lattice structures: Current trends and

Downloaded from https://academic.oup.com/jcde/article/12/8/107/8210400 by guest on 07 September 2025

future directions. Progress in Aerospace Sciences, 149, 101021. https://doi.org/10.1016/j.paerosci.2024.101021.
Khan, S. A., Rahman, M. A., Khraisheh, M., & Hassan, I. G. (2024). Advances in 3D printed periodic lattice structures for energy research: energy
storage, transport and conversion applications. Materials & Design, 239, 112773. https://doi.org/10.1016/j.matdes.2024.112773.
Kreisselmeier, G., & Steinhauser, R. (1980). Systematic control design by optimizing a vector performance index. In Cuenod M. A. (ed.), Computer
Aided Design of Control Systems, Proceedings of the IFAC Symposium, (pp. 113–117). Pergamon. https://doi.org/10.1016/B978- 0- 08- 024488- 4.500
22-X.
Lin, D., Gao, L., & Gao, J. (2025). The Lagrangian-Eulerian described Particle Flow Topology Optimization (PFTO) approach with isogeometric
material point method. Computer Methods in Applied Mechanics and Engineering, 440, 117892. https://doi.org/10.1016/j.cma.2025.117892.
Mantovani, S., Barbieri, S., Giacopini, M., Croce, A., Sola, A., & Bassoli, E. (2021). Synergy between topology optimization and additive manu-
facturing in the automotive field. Proceedings of the Institution of Mechanical Engineers, Part B: Journal of Engineering Manufacture, 235, 555–567.
https://doi.org/10.1177/0954405420949209.
Neves, M. M., Sigmund, O., & Bendsøe, M. P. (2002). Topology optimization of periodic microstructures with a penalization of highly localized
buckling modes. International Journal for Numerical Methods in Engineering, 54, 809–834. https://doi.org/10.1002/nme.449.
Nogueira, P., Lopes, P., Oliveira, L., Alves, J. L., Magrinho, J. P. G., Deus, A. M. d., Vaz, M. F., & Silva, M. B. (2024). Evaluation of lattice structures for
medical implants: a study on the mechanical properties of various unit cell types. Metals, 14, Article 7:780. https://doi.org/10.3390/met140
70780.
Prathyusha, A. L. R., & Raghu Babu, G. (2022). A review on additive manufacturing and topology optimization process for weight reduction
studies in various industrial applications. Materials Today: Proceedings, 62, 109–117. https://doi.org/10.1016/j.matpr.2022.02.604.
Ramadani, R., Pal, S., Kegl, M., Predan, J., Drstvenšek, I., Pehan, S., & Belšak, A. (2021). Topology optimization and additive manufacturing in
producing lightweight and low vibration gear body. The International Journal of Advanced Manufacturing Technology, 113, 3389–3399. https:
//doi.org/10.1007/s00170- 021- 06841- w.
Seharing, A., Azman, A. H., & Abdullah, S. (2014). A review on integration of lightweight gradient lattice structures in additive manufacturing
parts. Advances in Mechanical Engineering, 12. https://doi.org/10.1177/1687814020916951.
Thomsen, C. R., Wang, F., & Sigmund, O. (2018). Buckling strength topology optimization of 2D periodic materials based on linearized bifurcation
analysis. Computer Methods in Applied Mechanics and Engineering, 339, 115–136. https://doi.org/10.1016/j.cma.2018.04.031.
Viswanath, A., Khalil, M., Al Maskari, F., Cantwell, W. J., & Khan, K. A. (2023). Harnessing buckling response to design lattice structures with
improved buckling strength. Materials & Design, 232, 112113. https://doi.org/10.1016/j.matdes.2023.112113.
Viswanath, A., Khalil, M., Khan, M. K. A., Al Maskari, F., Cantwell, W. J., & Khan, K. A. (2024). A novel design strategy to enhance buckling
resistance of thin-walled single-cell lattice structures via topology optimisation. Virtual and Physical Prototyping, 19, e2345390. https://doi.
org/10.1080/17452759.2024.2345390.
Viswanath, A., Khalil, M., Khan, M. K. A., Cantwell, W. J., & Khan, K. A. (2025). Hierarchical cubic lattice structures with bending- and stretching-
dominated cellular designs for enhanced buckling resistance. International Journal of Lightweight Materials and Manufacture, 8, 310–320. https:
//doi.org/10.1016/j.ijlmm.2025.02.002.
Wang, F., Brøns, M., & Sigmund, O. (2023). Non-hierarchical architected materials with extreme stiffness and strength. Advanced Functional
Materials, 33, 2211561. https://doi.org/10.1002/adfm.202211561.
Wang, F., & Sigmund, O. (2021). 3D architected isotropic materials with tunable stiffness and buckling strength. Journal of the Mechanics and
Physics of Solids, 152, 104415. https://doi.org/10.1016/j.jmps.2021.104415.
Wang, X., Zhu, L., Sun, L., & Li, N. (2021). Optimization of graded filleted lattice structures subject to yield and buckling constraints. Materials
& Design, 206, 109746. https://doi.org/10.1016/j.matdes.2021.109746.
Xu, T., Huang, X., Lin, X., & Xie, Y. M. (2024). Topology optimization of continuum structures for buckling resistance using a floating projection
method. Computer Methods in Applied Mechanics and Engineering, 429, 117204. https://doi.org/10.1016/j.cma.2024.117204.
Yang, F., Li, P., Guo, Z., Li, X., Zhao, J., Wang, L., & Zhong, Z. (2024). Lattice metamaterials with controllable mechanical properties inspired by
projection of four-dimensional hypercubes. International Journal of Solids and Structures, 305, 113091. https://doi.org/10.1016/j.ijsolstr.2024.
113091.
Zhang, L., Feih, S., Daynes, S., Wang, Y., Wang, M. Y., Wei, J., & Lu, W. F. (2018). Buckling optimization of Kagome lattice cores with free-form
trusses. Materials & Design, 145, 144–155. https://doi.org/10.1016/j.matdes.2018.02.026.
Zhao, Y., Guo, G., & Zuo, W. (2023). MATLAB implementations for 3D geometrically nonlinear topology optimization: 230-line code for SIMP
method and 280-line code for MMB method. Structural and Multidisciplinary Optimization, 66, 146. https://doi.org/10.1007/s00158- 023- 03590- 0.
 Journal of Computational Design and Engineering, 2025, 12(8), 107–126 | 125

Zhu, J., Zhou, H., Wang, C., Zhou, L., Yuan, S., & Zhang, W. (2021). A review of topology optimization for additive manufacturing: Status and
challenges. Chinese Journal of Aeronautics, 34, 91–110. https://doi.org/10.1016/j.cja.2020.09.020.

Appendix A: Formulations

Downloaded from https://academic.oup.com/jcde/article/12/8/107/8210400 by guest on 07 September 2025

Figure A1: Isoparamteric formulation of 3D H8 trilinear element.

Table A1: Shape functions of 3D H8 trilinear element.

N1 = 1
8
(1 − ξ)(1 − η)(1 − ζ)

N2 = 1
8
(1 + ξ )(1 − η)(1 − ζ )
1
N3 = 8
(1 + ξ )(1 + η)(1 − ζ )
1
N4 = 8
(1 − ξ )(1 + η)(1 − ζ )
1
N5 = 8
(1 − ξ )(1 − η)(1 + ζ )
1
N6 = 8
(1 + ξ )(1 − η)(1 + ζ )
1
N7 = 8
(1 + ξ )(1 + η)(1 + ζ )
N8 = 1
8
(1 − ξ )(1 + η)(1 + ζ )
126 | Enhancing buckling resistance via three-dimensional topology optimization

Table A2: Mapping eigenvector components into array p for ϕT Keσ0 ϕ ≡ Z˜ p.

p1 (ϕ12 + ϕ22 + ϕ32 )

p2 (ϕ1 ϕ4 + ϕ2 ϕ5 + ϕ3 ϕ6 )
p3 (ϕ1 ϕ7 + ϕ2 ϕ8 + ϕ3 ϕ9 )
p4 (ϕ1 ϕ10 + ϕ2 ϕ11 + ϕ3 ϕ12 )
p5 (ϕ1 ϕ13 + ϕ2 ϕ14 + ϕ3 ϕ15 )
p6 (ϕ1 ϕ16 + ϕ2 ϕ17 + ϕ3 ϕ18 )
p7 (ϕ1 ϕ19 + ϕ2 ϕ20 + ϕ3 ϕ21 )
p8 (ϕ1 ϕ22 + ϕ2 ϕ23 + ϕ3 ϕ24 )
p9 (ϕ42 + ϕ52 + ϕ62 )
p10 (ϕ4 ϕ7 + ϕ5 ϕ8 + ϕ6 ϕ9 )
p11 (ϕ4 ϕ10 + ϕ5 ϕ11 + ϕ6 ϕ12 )
p12 (ϕ4 ϕ13 + ϕ5 ϕ14 + ϕ6 ϕ15 )
p13 (ϕ4 ϕ16 + ϕ5 ϕ17 + ϕ6 ϕ18 )

Downloaded from https://academic.oup.com/jcde/article/12/8/107/8210400 by guest on 07 September 2025

p14 (ϕ4 ϕ19 + ϕ5 ϕ20 + ϕ6 ϕ21 )
p15 (ϕ4 ϕ22 + ϕ5 ϕ23 + ϕ6 ϕ24 )
p16 (ϕ72 + ϕ82 + ϕ92 )
p17 (ϕ7 ϕ10 + ϕ8 ϕ11 + ϕ9 ϕ12 )
p18 (ϕ7 ϕ13 + ϕ8 ϕ14 + ϕ9 ϕ15 )
p19 (ϕ7 ϕ16 + ϕ8 ϕ17 + ϕ9 ϕ18 )
p20 (ϕ7 ϕ19 + ϕ8 ϕ20 + ϕ9 ϕ21 )
p21 (ϕ7 ϕ22 + ϕ8 ϕ23 + ϕ9 ϕ24 )
p22 2
(ϕ10 + ϕ11
2
+ ϕ12
2
)
p23 (ϕ10 ϕ13 + ϕ11 ϕ14 + ϕ12 ϕ15 )
p24 (ϕ10 ϕ16 + ϕ11 ϕ17 + ϕ12 ϕ18 )
p25 (ϕ10 ϕ19 + ϕ11 ϕ20 + ϕ12 ϕ21 )
p26 (ϕ10 ϕ22 + ϕ11 ϕ23 + ϕ12 ϕ24 )
p27 2
(ϕ13 + ϕ14
2
+ ϕ15
2
)
p28 (ϕ13 ϕ16 + ϕ14 ϕ17 + ϕ15 ϕ18 )
p29 (ϕ13 ϕ19 + ϕ14 ϕ20 + ϕ15 ϕ21 )
p30 (ϕ13 ϕ22 + ϕ14 ϕ23 + ϕ15 ϕ24 )
p31 2
(ϕ16 + ϕ17
2
+ ϕ18
2
)
p32 (ϕ16 ϕ19 + ϕ17 ϕ20 + ϕ18 ϕ21 )
p33 (ϕ16 ϕ22 + ϕ17 ϕ23 + ϕ18 ϕ24 )
p34 2
(ϕ19 + ϕ20
2
+ ϕ21
2
)
p35 (ϕ19 ϕ22 + ϕ20 ϕ23 + ϕ21 ϕ24 )
p36 2
(ϕ22 + ϕ23
2
+ ϕ24
2
)

Received: March 6, 2025. Revised: June 3, 2025. Accepted: June 7, 2025

© The Author(s) 2025. Published by Oxford University Press on behalf of the Society for Computational Design and Engineering. This is an Open Access article distributed
under the terms of the Creative Commons Attribution-NonCommercial License (https://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use,
distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact reprints@oup.com for reprints and
translation rights for reprints. All other permissions can be obtained through our RightsLink service via the Permissions link on the article page on our site-for further
information please contact journals.permissions@oup.com