Structural stability and artificial buckling modes in topology optimization

Abstract

This paper demonstrates how a strain energy transition approach can be used to remove artificial buckling modes that often occur in stability constrained topology optimization problems. To simulate the structural response, a nonlinear large deformation hyperelastic simulation is performed, wherein the fundamental load path is traversed using Newton’s method and the critical buckling load levels are estimated by an eigenvalue analysis. The goal of the optimization is to minimize displacement, subject to constraints on the lowest critical buckling loads and maximum volume. The topology optimization problem is regularized via the Helmholtz PDE-filter and the method of moving asymptotes is used to update the design. The stability and sensitivity analyses are outlined in detail. The effectiveness of the energy transition scheme is demonstrated in numerical examples.

1 Introduction

Designs optimized for weight or compliance often consist of slender members which are inherently prone to buckle. This failure mode has been considered in the optimization of truss and beam structures, dating back to the work of Khot et al. (1976) where the weight of space truss structures was minimized subject to stability constraints and Olhoff and Rasmussen (1977) who presented closed form expressions for the sufficient optimality condition of statically indeterminate columns. More recently, Madah and Amir (2017, 2019) improved the accuracy of the buckling predictions by using geometrically nonlinear beam theory to optimize stable truss structures.

Although stability is crucial for a robust design, it is often neglected in topology optimization of continuum structures due to its complicated analysis. Indeed, an unstable equilibrium configuration is characterized by a singular tangent stiffness matrix, and therefore the most accurate buckling load estimate is obtained by performing an incremental loading until at least one eigenvalue of the tangent stiffness matrix becomes zero (see, e.g. Suleman and Sedaghati (2005) and Pedersen and Pedersen (2018)). In topology optimization, this is a costly computation as it must be performed for each design iteration, and its accuracy is adversely affected by the artificial buckling modes that arise due to the ersatz material model.

To lessen the computational burden associated with the buckling analysis, it is common to employ a so-called “linear prebuckling analysis”, in which the buckling load estimate is obtained from a generalized eigenvalue problem about the structure’s undeformed configuration. Research on topology optimization that incorporates linear prebuckling analysis include the seminal work of Neves et al. (1995) where the lowest critical load is maximized. Zhou (2004) and Bruyneel et al. (2008) discuss convergence issues encountered in density based topology optimization. The necessity of including a large number of buckling constraints to obtain monotonic convergence is scrutinized by Dunning et al. (2016), as well as the need to use an efficient eigenvalue solver. Recently, Ferrari and Sigmund (2019, 2020) investigated the use of nonconforming finite elements and an aggregated stability constraint to solve large scale buckling constrained optimization problems.

Unfortunately, designs that are based on linear prebuckling analysis may not behave as expected due to the crude nature of the approximation. For more accurate buckling load estimates, a fully nonlinear analysis is required. In the nonlinear setting, there are two ways to identify the buckling load levels. The first, known as the extended system approach, directly computes the critical configurations by finding the unstable equilibrium solutions, cf. Wriggers (2008). The second approach, known as the one point approach, utilizes the information at a stable equilibrium configuration and approximates the buckling load level by solving a generalized eigenvalue problem based on a decomposition of the tangent stiffness matrix, cf. Brendel and Ramm (1980). The extended system approach provides an exact measure of the critical loads up to a numerical tolerance; however, it is only able to identify a single critical load. The one point approach does not suffer from this restriction and is therefore adopted in this work, despite the fact that it is an approximation. Reitinger and Ramm (1995) and Kemmler et al. (2005) employ the extended system approach to obtain critical buckling loads in their topology optimization work. The one point approach is used in continuum topology optimization by, e.g. Lindgaard and Dahl (2013), where different buckling load and compliance objectives are investigated, and by Lindgaard and Lund (2010, 2011) to optimize stable composite structures. For reasons discussed below, a large number of buckling constraints must be considered to obtain a smooth convergence of the optimization problem and hence we employ the one point approach.

Unfortunately, the solution to eigenvalue problems in topology optimization with an ersatz material model is complicated due to the existence of the low volume fraction regions with minimal stiffness. Nonphysical, artificial instabilities appearing in the low volume fraction regions must be distinguished from the actual physical modes. Several methods for accommodating the artificial buckling modes exist. A simple remedy used in the linear buckling analysis is presented by Neves et al. (1995), wherein the geometrical stiffness of low stiffness elements is ignored. This method might, however, cause convergence issues in the optimization due to discontinuities in the gradient information, and therefore Bendsøe and Sigmund (2003) propose an interpolation scheme to gradually remove the low stiffness elements. Both Zhou (2004) and Gao and Ma (2015) discuss the difficulty of choosing the interpolation parameters such that only the physical buckling loads are obtained.

The fictitious domain–based topology optimization considering natural frequencies in Du and Olhoff (2007) and Pedersen (2000) use separate interpolations of the stiffness and mass matrices. This idea could potentially be applied in the one point buckling analysis. However, Lindgaard and Dahl (2013) found that different interpolations for different stiffness matrix contributions result in loss of convergence in the equilibrium equations. Instead, Lindgaard and Dahl (2013) used the same interpolation for all components of the stiffness matrix, and eliminated the artificial modes from the analysis by increasing the stiffness of the void regions. Of course, this method results in loss of precision in the sense that void regions will be artificially stiff. As an alternative, Gao and Ma (2015) distinguished between artificial and physical buckling modes by using a strain energy density criterion. We note that this approach is computationally intensive.

The low volume fraction regions pose additional problems when subjected to finite deformations, as excessively distorted finite elements might result in loss of convergence in the Newton equilibrium iterations. Therefore, in a fully nonlinear analysis, it is attractive to find an approach which deals with both the artificial modes and the equilibrium divergence simultaneously. One such approach is the element removal strategy, as done by Dalklint et al. (2020) in their topology optimization of preloaded elastic continuum structures subject to natural frequency constraints. Herein, we use the energy transition scheme presented by Wang et al. (2014), whereby the physical strain energy is replaced by a linearized strain energy in void regions. As such, we eliminate the artificial modes and avoid equilibrium convergence issues.

Summarizing, in our work, we obtain accurate estimates of the buckling loads of compressible neo-Hookean hyperelastic structures by using the one point approach wherein a generalized eigenvalue problem about the deformed equilibrium configuration is solved to approximate the buckling loads. The deformed configuration is found by incrementally loading the structure to approach, but not reach, an unstable configuration. Newton’s method facilitates this analysis. The optimization minimizes the end-compliance at the operating load, subject to lower bound constraints on the n-lowest critical loads and an upper bound constraint on the mass. The method of moving asymptotes (MMA) nonlinear programming algorithm is used to update the design, and the sensitivity analysis is performed using the adjoint approach. To address the nondifferentiability of degenerate eigenvalues, we follow the work of Chin and Kennedy (2016) and Ferrari and Sigmund (2019) and replace the, possibly nondifferentiable, n eigenvalue constraints with a differentiable aggregated constraint.

2 Governing equations

To consider finite deformations, we need to distinguish between the undeformed, reference configuration, Ω, and the deformed, current configuration, Ω_c. Assuming Lagrangian elasto-statics, each material point identified by the position vector X ∈Ω in our n_dim-dimensional space, is displaced to the position x ∈Ω_c, by the deformation mapping $\boldsymbol {\varphi }: {{\varOmega }} \rightarrow {{\varOmega }}_{c}$, i.e.

$$ \boldsymbol{x}=\boldsymbol{\varphi}(\boldsymbol{X})=\boldsymbol{u}(\boldsymbol{X})+\boldsymbol{X}, $$

(1)

where $\boldsymbol {u}:{{\varOmega }} \rightarrow \mathbb {R}^{n_{dim}}$ represents the displacement field. To describe the strain in this nonlinear theory, we define the deformation gradient

$$ \boldsymbol{F}= \mathbf{1} + \boldsymbol{\nabla}{\boldsymbol{u}}, $$

(2)

where ∇ is the material gradient operator on Ω and 1 is the second order identity tensor. The local deformation is described by the right Green-Cauchy deformation tensor C = F^TF. For later purposes we also introduce the linearized strain $\boldsymbol {\epsilon }=\frac {1}{2}\left (\left (\boldsymbol {\nabla }\boldsymbol {u}\right )^{T} + \boldsymbol {\nabla }\boldsymbol {u}\right )$.

Under the assumption of vanishing body forces, the potential energy of our hyperelastic body is

$$ {{\varPi}}(\boldsymbol{u})= {\int}_{{{\varOmega}}} W \ dV -{\int}_{\partial {{\varOmega}}_{t}}\bar{\boldsymbol{t}} \cdot \boldsymbol{u}\ dS , $$

(3)

where W = W(∇u) is the specific strain energy and the boundary ∂Ω with unit normal n is partitioned into complementary subsets ∂Ω_t and ∂Ω_u, over which the traction $\boldsymbol {t}=\bar {\boldsymbol {t}}$ and the displacement ${\boldsymbol {u}}=\boldsymbol {\varphi }-\boldsymbol {X}=\bar {\boldsymbol {u}}$, are prescribed, respectively. The equilibrium configurations of Ω are described by the admissible displacement fields $\boldsymbol {u} \in \{\boldsymbol {u}\in H^{1}: \boldsymbol {u}(\boldsymbol {X}) =\bar {\boldsymbol {u}} \ \text {for} \ \boldsymbol {X} \in \partial {{\varOmega }}_{u}\}$ that minimize the energy, i.e. satisfy the stationary condition

$$ \begin{array}{@{}rcl@{}} \delta{{\varPi}}(\boldsymbol{u};\delta \boldsymbol{u} )&=&{\int}_{{{\varOmega}}} \frac{\partial W}{\partial \boldsymbol{\nabla} \boldsymbol{u}}: \boldsymbol{\nabla}\left( \delta \boldsymbol{u}\right) dV \\ &&-{\int}_{\partial {{\varOmega}}_{t}}\bar{\boldsymbol{t}} \cdot \delta \boldsymbol{u}\ dS = 0, \end{array} $$

(4)

for all smooth virtual displacements $\delta \boldsymbol {u} \in \{\delta \boldsymbol {u}\in H^{1}: \delta \boldsymbol {u}(\boldsymbol {X}) = \textsf {\textbf {0}} \ \text {for} \ \boldsymbol {X} \in \partial {{\varOmega }}_{u}\}$ and where H¹ is a Hilbert space.

3 Material interpolation and design variables

The goal of topology optimization is to optimally distribute material in a design domain Ω. To quantify the material distribution, we use the material indicator function ${{\varGamma }}: {{\varOmega }} \rightarrow \{0,1\}$. However, as is standard in topology optimization, we replace the binary valued material indicator function with a continuous nondimensional volume fraction field $z : {{\varOmega }} \rightarrow [0,1]$ and to obtain a well-posed optimization problem, we restrict the design space. Here, the restriction is performed by penalizing fine-scale oscillations of z via the Helmholtz PDE-filter proposed in Lazarov and Sigmund (2011). As such, z is replaced by the smoothed field $\nu : {{\varOmega }} \rightarrow [0,1]$, obtained from the solution of the boundary-value problem

$$ \begin{array}{@{}rcl@{}} -{l}^{2}{{\varDelta}} {\nu} + {\nu}&=&z \ \ \text{in} \ \ {{\varOmega}}, \\ \text{s.t.} \ \ \boldsymbol{\nabla} \nu \cdot \boldsymbol{n}&=&0 \ \ \text{on} \ \ \partial{{\varOmega}}, \end{array} $$

(5)

where Δ is the Laplacian with respect to Ω and the oscillations of z are further smoothed by increasing the value of $l\in \mathbb {R}^{+}$. The homogeneous Neumann boundary condition in (5) is well-known to cause an artificial “attraction” of material to the domain boundaries. This anomaly can be mitigated by applying a Robin condition, cf. Wallin et al. (2020). In our discretization, we parameterize z via the piece-wise uniform design field $\textsf {z}: {{\varOmega }} \rightarrow [0,1]$ over the finite elements Ω^e ∈Ω, such that z(X) = z_e for X ∈Ω^e. We use standard Lagrange finite elements to solve (5) for the continuous smoothed field ν, cf. Lazarov and Sigmund (2011) for details.

Unfortunately, the filtering produces gray regions wherein ν(X) ∈ (0,1). To limit the extent of these regions, we implement the thresholding technique proposed by Guest et al. (2004) and Wang et al. (2014), wherein a smoothed Heaviside function $H_{\beta ,\eta } : \mathbb {R} \rightarrow [0,1]$ is employed to define the thresholded filtered volume fraction

$$ \bar{\nu}=H_{\beta,\eta}(\nu)=\frac {\tanh \left( \beta\eta \right) + \tanh \left( \beta \left( \nu-\eta \right) \right) }{\tanh \left( \beta\eta \right) +\tanh \left( \beta \left( 1-\eta \right) \right) }, $$

(6)

where β and η are numerical parameters defined such that $\lim _{\beta \rightarrow \infty } H_{\beta ,\eta }(\nu ) = u_{s}(\nu -\eta )$ where u_s is the unit step function. The filter (5) and the Heaviside approximation (6) are used to obtain the differentiable field $\bar {\nu }$. Indeed, although z is a piecewise uniform function, ν and $\bar {\nu }$ are continuous functions that vary spatially throughout the domain Ω.

To further reduce the size of the gray regions for which $\bar {\nu }(\boldsymbol {X})\in (0,1)$ we employ the SIMP penalization scheme, cf. Bendsøe (1989), where a nonlinear scaling of the elastic material parameters, K and G, which describe the isotropic hyperelastic material model is introduced as

$$ G=\chi(\bar{\nu})G_{o}, \quad K=\chi(\bar{\nu})K_{o}. $$

(7)

In (7), G_o and K_o denote the shear and bulk moduli of the optimized structure’s isotropic material, in the limit of infinitesimal strain. The SIMP scaling function $\chi : [0,1] \rightarrow [\delta _{o},1]$, is defined as

$$ \chi({{\bar{\nu}}})= \bar{\nu}^{q}(1-\delta_{o})+ \delta_{o}, $$

(8)

whereby increasing values of q > 1 enforce increasing levels of penalization and the lower threshold δ_o = 10^− 9 in this ersatz material model is used to avoid numerical problems otherwise caused by a vanishing material stiffness.

The aforementioned material interpolation allows for nearly void regions to be part of the physical domain. As such, even minimal stress leads to excessive distortion in these regions resulting in loss of convergence in the finite element equilibrium simulations. To mitigate this issue, we use the energy transition strategy proposed by Wang et al. (2014). In this scheme, the deformation gradient is redefined as F = 1 + γ∇u, where $\gamma : [0,1]\rightarrow [0,1]$ is another penalty function defined as

$$ \gamma(\bar{\nu})= H_{\tilde{\beta},\tilde{\eta}}(\bar{\nu}^{q}), $$

(9)

where $\tilde {\beta }=500$ and $\tilde {\eta }=0.01$ replaces β and η in (6).

To formulate stress and stiffness relations consistent with the strain energy transition, we weight the physical strain energy, W^PH, against the linearized energy, W^L, such that the body’s internal energy is

$$ \begin{array}{@{}rcl@{}} W(\boldsymbol{\nabla}\boldsymbol{u};\gamma)& =&{W}^{PH}(\gamma\nabla\boldsymbol{u})+{W}^{L}(\nabla\boldsymbol{u})-{W}^{L}(\gamma\nabla\boldsymbol{u} ) \\ &=& {W}^{PH}(\gamma\nabla\boldsymbol{u})+\left( 1-\gamma^{2}\right){W}^{L}(\nabla\boldsymbol{u} ). \end{array} $$

(10)

By assigning points devoid of material γ = 0, we therefore obtain W = W^L. The physical strain energy is modelled using an isotropic compressible neo-Hookean hyperelastic constitutive model such that

$$ \begin{array}{@{}rcl@{}} {W}^{PH}&=&\frac{1}{4}K_{o} \left( \left( J^{2} -1\right) - \ln\left( J^{2}\right) \right) \\ &&+ \frac{1}{2}G_{o}\left( J^{-2/3} \text{tr}(\boldsymbol{C}) -3\right). \end{array} $$

(11)

Its linearization yields the linearized energy

$$ \begin{array}{@{}rcl@{}} W^{L} &=& \frac{1}{2} \boldsymbol{\epsilon}:\left( \left.4\frac{\partial^{2} W^{PH}}{\partial \boldsymbol{C}\partial \boldsymbol{C}}\right|_{\boldsymbol{C}=\mathbf{1}}\right):\boldsymbol{\epsilon} \\ &= &\frac{1}{2} \boldsymbol{\epsilon}:\boldsymbol{D}^{L}:\boldsymbol{\epsilon}, \end{array} $$

(12)

where D^L is the linear stiffness tensor corresponding to an isotropic material, with bulk and shear modulii K_o and G_o.

4 Total Lagrangian FE-formulation and solution procedure

Our equilibrium equations are discretized using standard finite elements, see, e.g. Crisfield (1993). We use to represent the finite element assembly operator, lower case indices $e\in \mathbb {N}_{e}$ to indicate quantities associated with finite element Ω^e, and lower case $n \in \mathbb {N}_{n}$ to denote the finite element degrees-of-freedom. The sets of finite element degrees-of-freedoms and finite elements, are denoted $\mathbb {N}_{n}$ and $\mathbb {N}_{e}$, respectively.

We seek the discretized displacement field $\textsf {\textbf {u}} \in \mathbb {R}^{n_{dim}}$ that approximates the displacement u via standard Lagrange shape functions N within the element Ω^e such that u(X) = N(X)a_e, where a_e is the element nodal degrees-of-freedom vector. The finite element formulation is based on the Galerkin approach, implying δu(X) = N(X)δa_e, where δa_e is the virtual element nodal degrees-of-freedom vector.

To discretize (4), we first note that from (10) the virtual stress power is

$$ \frac{\partial W}{\partial \boldsymbol{\nabla}\boldsymbol{u}}: \boldsymbol{\nabla}\delta\boldsymbol{u} = \boldsymbol{S}^{PH}:\delta {\boldsymbol{E}}+\left( 1-\gamma^{2}\right) \boldsymbol{\sigma}^{L}:\delta\boldsymbol{\epsilon}, $$

(13)

where

$$ \boldsymbol{S}^{PH}(\gamma\boldsymbol{\nabla} \boldsymbol{u}) = 2\frac{\partial W^{PH}}{\partial \boldsymbol{C}}, \quad \boldsymbol{\sigma}^{L}(\boldsymbol{\nabla} \boldsymbol{u}) =\boldsymbol{D}^{L}\boldsymbol{\epsilon}, $$

(14)

are the physical and linearized stress tensors. In (13), we define the virtual Green-Lagrange strain $\delta \boldsymbol {E}=\frac {\gamma }{2}\left (\left (\boldsymbol {\nabla } \delta \boldsymbol {u}\right )^{T} \boldsymbol {F}+\boldsymbol {F}^{T} \boldsymbol {\nabla } \delta \boldsymbol {u}\right )$ and the virtual linearized strain $\delta \boldsymbol {\epsilon } = \frac {1}{2}\left (\left (\boldsymbol {\nabla } \delta \boldsymbol {u}\right )^{T} +\boldsymbol {\nabla } \delta \boldsymbol {u}\right )$. The discretization of δu is used to express these tensors as $\delta \textsf {\textbf {E}} = \textsf {\textbf {B}}(\gamma \textsf {\textbf {a}}_{e}) \gamma \delta \textsf {\textbf {a}}_{e} = \left (\textsf {\textbf {B}}^{c} + \textsf {\textbf {B}}^{l}(\gamma {\textsf {\textbf {a}}}_{e})\right )\gamma \delta {\textsf {\textbf {a}}}_{e}$ and $\delta \boldsymbol {\epsilon }=\textsf {\textbf {B}}^{c}\delta {\textsf {\textbf {a}}}_{e}$, where B is expressed as the sum of a constant B^c and a linear function of the displacement $\textsf {\textbf {B}}^{l}(\gamma {\textsf {\textbf {a}}}_{e})$, cf. Crisfield (1993).

Under the assumption of dead loading, we express the external load as λP where $\textsf {\textbf {P}}\in \mathbb {R}^{n}$ is a unit magnitude external load vector and $\lambda \in \mathbb {R}^{+}$ is the scalar load intensity factor. Using the arbitrariness of the virtual nodal displacements, δa, and (13), the discretization of (4) is expressed as

where $\textsf {\textbf {r}}(\textsf {\textbf {a}})\in \mathbb {R}^{n}$ is the residual vector.

We find the degree-of-freedom vector a in (??) for the given λ. To this end, we employ Newton’s method to solve (??). As such, we iterate solving K(a)Δa = −r(a) for Δa and updating ${\textsf {\textbf {a}}} \rightarrow {\textsf {\textbf {a}}} + {{\varDelta }} {\textsf {\textbf {a}}}$ until convergence, where the tangent stiffness matrix $\textsf {\textbf {K}}({\textsf {\textbf {a}}})\in \mathbb {R}^{n \times n}$ is

and where $\boldsymbol {D}= 4\frac {\partial ^{2} W^{PH}}{\partial \boldsymbol {C}\partial \boldsymbol {C}}$, G contains the gradient of the element shape functions, Y is a symmetric block matrix of S^PH(γa) and $\left (\textsf {\textbf {D}}, \textsf {\textbf {D}}^{L}\right )$ are the Voigt representations of $\left ({{{\boldsymbol {D}}}}, {{{\boldsymbol {D}}}}^{L}\right )$ respectively. After convergence, we increase λ and again evaluate a. We repeat this process until we reach the operating load, i.e. until λ = λ^op.

5 Structural stability

A conservative system is stable if the potential energy has a (local) isolated minimum. This condition is equivalent to the requirement that the second variation of the potential energy is positive definite. For our system, this implies that all eigenvalues of the tangent stiffness matrix (15) are strictly positive. On the other hand, the tangent stiffness matrix has negative eigenvalues if the system is unstable. A system that transitions from being stable to unstable must pass a critical, i.e. singular point, for which at least one eigenvalue of K vanishes, i.e. when

$$ \textsf{\textbf{K}} \pmb{\upphi}_{j} = \textsf{\textbf{0}}, \ \ j \in \mathbb{N}_{N}, $$

(17)

has a nontrivial solution ϕ_j. Unfortunately, the zero eigenvalue condition is difficult to enforce in topology optimization due to the appearance of artificial eigenmodes attributed to the ersatz material. As such, the search for the zero eigenvalue would require many eigenvalue solves and a strategy to distinguish between artificial and physical eigenmodes.

In our alternate strategy, we find the critical physical eigenvalues by approximating the eigenvalue problem (17), as a generalized eigenvalue problem that is based on a decomposition of the stiffness matrix. Motivated by Wriggers (2008), we assume that K is the sum of the stress dependent contribution K_g that varies linearly with respect to the load level in the vicinity of the critical point, and the constant contribution K_o. To leverage this assumption, a scalar $\uppi \in \mathbb {R}^{+}$ is introduced such that λ^c =π λ^op where $\lambda ^{\text {op}}\in \mathbb {R}$ denotes the operating load level at which we evaluate

$$ \textsf{\textbf{K}}\approx \textsf{\textbf{K}}_{o}+{\uppi}_{j} \textsf{\textbf{K}}_{g}, $$

(18)

where

and

In this work, we limit the analysis to buckling load factors π > 1, i.e. we assume that the decomposition (18) is performed in a stable equilibrium position, prior to reaching any instabilities.

We use the assumption (18), to approximate the eigenvalue problem (17) via the generalized eigenvalue problem¹

$$ \textsf{\textbf{K}}_{o}\pmb{\upphi}_{j} = -{\uppi}_{j} \textsf{\textbf{K}}_{g}\pmb{\upphi}_{j} \ \ j \in \mathbb{N}_{n}, $$

(21)

where $({\uppi }_{j},\pmb {\upphi }_{j}) \in (\mathbb {R}\times \mathbb {R}^{n})$ are the eigenpairs. The eigenvalues π_j are ordered in ascending order such that π₁ ≤π₂ ≤… ≤π_j ≤π_n and the eigenvectors are normalized such that

$$ \pmb{\upphi}_{i}^{T} \textsf{\textbf{K}}_{g} \pmb{\upphi}_{j} = -\delta_{ij},\ \ i,j \in \mathbb{N}_{n}. $$

(22)

It is emphasized that in the case of N-fold degenerate eigenvalues, the eigenvectors are not uniquely defined as they span an N-dimensional hyperplane, i.e. any linear combination of ϕ_j, $j\in \mathbb {N}_{N}$, satisfy (21) where the set $\mathbb {N}_{N}\subseteq \mathbb {N}_{n}$ contains the N indices j for which π_j = π_j+ 1 = … = π_j+N− 1.

6 Optimization formulation

The optimization problem is stated as finding the design z that solves

$$ \mathcal{P} \left\{ \begin{array}{l} \underset{\textsf{\textbf{z}}}{\min} \ g_{0}(\textsf{\textbf{z}}) \\ \text{s.t.} \left\{\begin{array}{l} g_{i}(\textsf{\textbf{z}})\leq 0, \ \ \ \ i \in \mathbb{N}_{c} \\ 0 \leq \textsf{z}_{e} \leq 1, \ \ e \in \mathbb{N}_{e} \end{array} \right. \end{array} \right. $$

(23)

where g₀ is the cost function to be minimized and g_i(z) ≤ 0 are the $i\in \mathbb {N}_{c}=\{1,2,\hdots \}$ constraint inequalities.

A common cost function to minimize in topology optimization, and the one used herein, is the end-compliance defined by

$$ g_{0}=\textsf{\textbf{P}}^{T}\textsf{\textbf{a}}. $$

(24)

Since P is independent of the load magnitude, the objective minimizes the displacement over the regions where the load is applied. Alternatively, one could maximize the tangent stiffness as done in Wallin et al. (2018).

For all examples considered, the amount of available material is limited by imposing the constraint

$$ g_{1}={\int}_{{{\varOmega}}} \bar{\nu} dv- {V}_{max}\leq 0, $$

(25)

where ${V}_{max} \leq V = \displaystyle {\int \limits }_{{{\varOmega }}} dv$ is the maximum volume available for the design.

Ideally, the stability constraints are defined as

$$ g_{i}=s\lambda^{\text{op}} - {\lambda_{i}^{c}}\leq 0, \ \ \Leftrightarrow \ \ {\uppi}_{i} \geq s, \ \ i \in \mathbb{N}_{\lambda}, $$

(26)

where ${\lambda ^{c}_{i}} = {\uppi }_{i}\lambda ^{\text {op}}$ represent the critical load level, and $\mathbb {N}_{\lambda } \subset \mathbb {N}_{n}$ is a subset of eigenvalues, i.e. buckling load levels, considered. In (26), we introduce constraints on the $\mathbb {N}_{\lambda }$ first critical loads to avoid convergence issues associated with mode-crossing, cf. Dunning et al. (2016). The operating load level λ^op and the safety factor s ≥ 1 are problem specific parameters. We follow Ferrari and Sigmund (2019) and Maharaj and James (2020) and employ an alternative representation of (26), whereby the $\mathbb {N}_{\lambda }$-constraints (26) are replaced by the single constraint

$$ g_{2} = s\lambda^{\text{op}} - \tilde{\lambda}^{c} \leq 0, $$

(27)

where $\tilde {\lambda }^{c} = \frac {1}{\left \lVert \boldsymbol {{\varLambda }}\right \rVert _{p}}$, $\boldsymbol {{\varLambda }} = \left [\frac {1}{{\lambda ^{c}_{1}}}, \frac {1}{{\lambda ^{c}_{2}}}, \hdots , \frac {1}{{\lambda ^{c}_{n}}} \right ]$ and ∥Λ∥_p is the p-norm,² i.e.

$$ \left\lVert\boldsymbol{{\varLambda}}\right\rVert_{p} = \left( \sum\limits_{i\in \mathbb{N}_{\lambda}} \left( \frac{1}{{\lambda^{c}_{i}}} \right)^{p} \right)^{1/p}, \ \ p \in \mathbb{N}. $$

(28)

In (28), we take advantage of the strictly positive ${\lambda ^{c}_{i}}$, i.e. $\left \vert \frac {1}{{\lambda ^{c}_{i}}}\right \vert = \frac {1}{{\lambda ^{c}_{i}}}$. By denoting the smallest buckling load level $\lambda ^{c}_{\min \limits } =\underset {i \in \mathbb {N}_{\lambda }}{\min \limits } {\lambda ^{c}_{i}}$ we note that $\lambda ^{c}_{\min \limits }= \underset {p \rightarrow \infty }{\lim } \tilde {\lambda }^{c}$ and $\tilde {\lambda }^{c} \leq \lambda ^{c}_{\min \limits }$, i.e. $\tilde {\lambda }^{c}$ approaches $\lambda ^{c}_{\min \limits }$ from below. Thus, the (27) constraint is conservative. The use of (27) instead of (26) will be justified when the sensitivity analysis is discussed.

7 Sensitivity analysis

To evaluate the gradients of the objective and constraint functions, we use the adjoint procedure. In the followin, we obtain the sensitivities $\frac {\partial g_{i}}{\partial \pmb {\upnu }}$. But remember that the optimization requires the sensitivities with respect to the design variables z, i.e. $\frac {\partial g_{i}}{\partial \textsf {\textbf {z}}}$. These are obtained from the chain rule

$$ \frac{\partial g_{i}}{\partial \textsf{\textbf{z}}} = \frac{\partial g_{i}}{\partial {\pmb{\upnu}}}\frac{\partial {\pmb{\upnu}}}{\partial \textsf{\textbf{z}}}, \quad i \in[0,\mathbb{N}_{c}], $$

(29)

where the product $\frac {\partial g_{i}}{\partial {\pmb {\upnu }}}\frac {\partial {\pmb {\upnu }}}{\partial \textsf {\textbf {z}}}$ is obtained from an adjoint sensitivity analysis using (5), cf. Lazarov and Sigmund (2011).

7.1 Sensitivity of objective

The sensitivity of the objective function, i.e $\frac {\partial g_{o}}{\partial {\pmb {\upnu }}}$ is obtained by augmenting g_o with the zero valued residual vector by means of the adjoint vector $\boldsymbol {\upvarphi }\in \mathbb {R}^{n}$, such that

$$ g_{0}=\textsf{\textbf{P}}^{T}\textsf{\textbf{a}}-{\boldsymbol{\upvarphi}}^{T}\textsf{\textbf{r}}. $$

(30)

Differentiating (30) with respect to ν results in

$$ \frac{\partial g_{0}}{\partial {\pmb{\upnu}}}= \left( \textsf{\textbf{P}}^{T}-{\boldsymbol{\upvarphi}}^{T}\textsf{\textbf{K}}\right) \frac{\partial \textsf{\textbf{a}}}{\partial {\pmb{\upnu}}}-{\boldsymbol{\upvarphi}}^{T}\frac{\partial \textsf{\textbf{F}}^{int}}{\partial {\pmb{\upnu}}}. $$

(31)

To eliminate the implicitly defined derivative $\frac {\partial \textsf {\textbf {a}}}{\partial {\pmb {\upnu }}}$, we require φ to solve the adjoint problem

$$ \textsf{\textbf{K}}{\boldsymbol{\upvarphi}}=\textsf{\textbf{P}}, $$

(32)

where we utilize the symmetry of K. The final expression for the sensitivity of g₀ is thus reduced to

$$ \frac{\partial g_{0}}{\partial {\pmb{\upnu}}}=-\frac{\partial \left( \hat{{\boldsymbol{\upvarphi}}}^{T} \textsf{\textbf{F}}^{int}\right)}{\partial \pmb{\upnu}}, $$

(33)

where the notation $\hat {(\cdot )}$ here and henceforth denotes quantities which are treated as constants in the differentiation.

7.2 Sensitivity of critical loads

To obtain the sensitivity of the stability constraint, g₂, we start by differentiating (27)

$$ \frac{\partial g_{2}}{\partial {\pmb{\upnu}}}= - \frac{\partial \tilde{\lambda}^{c}}{\partial {\pmb{\upnu}}} = \sum\limits_{j\in \mathbb{N}_{\lambda}} \left( \frac{{\varLambda}_{j}}{\|{\boldsymbol{{\varLambda}}}\|_{p}} \right)^{p+1} \frac{{\partial\lambda_{j}^{c}}}{\partial \pmb{\upnu}} , $$

(34)

which follows from (28), and where ${{\varLambda }}_{j} = \frac {1}{{\lambda _{j}^{c}}}$. As the operating load is constant, the sensitivity $\frac {\partial {\lambda _{j}^{c}}}{\partial \pmb {\upnu }}$ is expressed as

$$ \frac{\partial {\lambda_{j}^{c}}}{\partial {\pmb{\upnu}}} = \lambda^{\text{op}}\frac{\partial {\uppi}_{j}}{\partial {\pmb{\upnu}}} , \ \ j\in\mathbb{N}_{\lambda}, $$

(35)

cf. the discussion preceding (17). The eigenvalue derivative $\frac {\partial {\uppi }_{j}}{\partial {\pmb {\upnu }}}$ is obtained by manipulating (21) and using the equilibrium residual equation to obtain the equality

$$ \pmb{\upphi}_{j}^{T} \left( \textsf{\textbf{K}}_{o}+ {\uppi}_{j}\textsf{\textbf{K}}_{g}\right) \pmb{\upphi}_{j}-{\boldsymbol{\upmu}}_{j}^{T}\left( \textsf{\textbf{F}}^{int}-\lambda^{\text{op}} \textsf{\textbf{P}}\right)=0, \quad j\in\mathbb{N}, $$

(36)

where ${\boldsymbol {\upmu }} \in \mathbb {R}^{n}$ is another adjoint vector. Differentiation of (36) with respect to ν results in

$$ \begin{array}{@{}rcl@{}} &&\pmb{\upphi}^{T}_{j} \left[\frac{\partial \textsf{\textbf{K}}_{o}}{\partial {\pmb{\upnu}}} +{\uppi}_{j}\frac{\partial \textsf{\textbf{K}}_{g}}{\partial {\pmb{\upnu}}} +\textsf{\textbf{K}}_{g}\frac{\partial {\uppi}_{j}}{\partial {\pmb{\upnu}}} +\left( \frac{\partial \textsf{\textbf{K}}_{o}}{\partial \textsf{\textbf{a}}}+{\uppi}_{j}\frac{\partial \textsf{\textbf{K}}_{g}}{\partial \textsf{\textbf{a}}}\right)\frac{\partial \textsf{\textbf{a}}}{\partial \pmb{\upnu}}\right] \pmb{\upphi}_{j} \\ &&-{\boldsymbol{\upmu}}^{T}_{i}\left( \frac{\partial \textsf{\textbf{F}}^{int}}{\partial {\pmb{\upnu}}} +\textsf{\textbf{K}} \frac{\partial \textsf{\textbf{a}}}{\partial {\pmb{\upnu}}} \right)=\textsf{\textbf{0}}. \end{array} $$

(37)

Similar to (32), we require the adjoint μ_j vector to solve

$$ \textsf{\textbf{K}} {\boldsymbol{\upmu}}_{j}= \frac{\partial \left( \hat{\pmb{\upphi}}_{j}^{T}\left( \textsf{\textbf{K}}_{o}+\hat{\uppi}_{j}\textsf{\textbf{K}}_{g}\right)\hat{\pmb{\upphi}}_{j}\right)}{\partial {\textsf{\textbf{a}}}}, $$

(38)

which annihilates the implicit sensitivity, $\frac {\partial \textsf {\textbf {a}}}{\partial {\pmb {\upnu }}}$, in (36). Insertion of (38) into (37) and making use of the normalization (22) results in

$$ \frac{\partial{\uppi}_{j}}{\partial {\pmb{\upnu}}}= \frac{\partial\left( \hat{\pmb{\upphi}}_{j}^{T}\left( \textsf{\textbf{K}}_{o} + {\uppi}_{j}\textsf{\textbf{K}}_{g}\right)\hat{\pmb{\upphi}}_{j}\right)}{\partial \pmb{\upnu}} -\frac{\partial \left( \hat{{\boldsymbol{\upmu}}}_{j}^{T}\textsf{\textbf{F}}^{int}\right)}{\partial {\pmb{\upnu}}}. $$

(39)

The explicit computations in the right hand sides of (38) and (39) appear in the 1.

It is well known that (39) does not hold for degenerate eigenvalues as they are not differentiable. Nonetheless, the simple eigenvalue result is adequate because both K_o and K_g are real, symmetric and continuously differentiable with respect to ν and g₂ is a symmetric polynomial of the eigenvalues π_j for $j \in \mathbb {N}_{\lambda }$. As such, the gradient $\frac {\partial \tilde {\lambda }_{j}^{c}}{\partial \pmb {\upnu }}$ exists even for the case of degenerate eigenvalues, by evaluating the sensitivity $\frac {\partial g_{2}}{\partial \pmb {\upnu }}$ via (34), (35) and (39), cf. Torii and De Faria (2017) and Gravesen et al. (2011).

8 Artificial buckling modes

If the material properties in void regions are allowed to strictly assume the zero lower bound, i.e δ_o = 0 in (8), the equilibrium equation cannot be solved. To avoid this, it is common to rely on an ersatz material model such that material points outside the physical domain are assigned a “small” stiffness, i.e. 0 < δ_o ≪ 1. For eigenvalue problems in topology optimization, the small stiffness in void regions renders nonphysical, i.e. artificial eigenmodes.

This artificial mode anomaly is readily shown in the following example, wherein a portal frame structure, similar to that studied in Gao and Ma (2015), is analyzed, cf. Fig. 1. The design domain is discretized using 40 × 20 × 1 3D 8-node fully integrated displacement based finite elements. Following Gao and Ma (2015), plane stress conditions are simulated by assigning the out-of-plane displacement on the front x-y plane to zero and enforcing traction free boundary conditions on the opposing face. The operating load level is λ^op = 60 N and the elastic parameters are E = 3 GPa and ν = 0.4 which provide G_o = E/(2(1 + ν)) and K_o = E/3(1 − 2ν). In this example, we disable the filter (5), and assign the outer frame, i.e. the dark grey physical domain of thickness 5 mm, a fully dense, i.e. ρ = 1, material, whereas ρ = 0.001 for the inner “void” region (light grey). The SIMP penalization parameter is q = 3 and the Heaviside projection parameters are β = 1 and η = 0.5, cf. (8) and (6). We search for eigenvalues π_j in the interval I = [1,100]. The physical buckling modes and load factors, are evaluated in three different ways. Firstly, the computations are performed on a conforming mesh over the physical domain, i.e. excluding the inner grey region. Secondly, the grey void domain is modelled via the SIMP ersatz model, cf. (8), without the energy transition scheme, i.e. by setting γ = 1. Lastly, the void domain is modelled via the combined SIMP and strain energy transition schemes of (8) and (10). Various matrix penalizations in the generalized eigenvalue problem (21) exist, cf., e.g. Lindgaard and Dahl (2013). However, as noted by Gao and Ma (2015) and Lindgaard and Dahl (2013), these methods deteriorate the convergence of the equilibrium analysis and we therefore exclude them.

Fig. 1

Portal frame structure. The thickness is 10 mm

The smallest buckling load factor, π₁, together with the number of modes, j_max, found in the interval I, are presented in Table 1. As expected, we find a large number of artificial modes in the inner “void” region when using the second method. This behaviour is readily predicted by inspecting the Rayleigh quotient of (21)

Table 1 The magnitude of the smallest load factor π₁ and the number of eigenvalues j_max found in the interval I for the three modeling schemes

When employing the strain energy transition scheme, the denominator in the Rayleigh quotient (40) tends to zero as $\gamma \rightarrow 0$, thereby removing the artificial modes from the interval I. The idea of identifying artificial modes using the Rayleigh quotient also appears in Gao and Ma (2015). The first buckling modes for the three simulations are depicted in Fig. 2. We emphasize that Fig. 2b displays an out-of-plane artificial buckling mode.

Fig. 2

The first buckling mode when using a conforming mesh (a), the SIMP ersatz material (b) and combined SIMP and strain energy transition models (c)

To conclude the investigation, we compute the ratio between the greatest positive physical eigenvalue and smallest artificial eigenvalue obtained using the proposed method. To this end, we extend I = [1,10²⁰] and compute all positive physical eigenvalues using the first method, whereby we find 391 positive physical eigenvalues. Subsequently, we compute the eigenvalues in I using the third method. As expected, the first 391 eigenvalues are physical, whereas the rest are artificial. The first artificial eigenvalue is approximately 11 magnitudes larger than the greatest physical eigenvalue, which further demonstrates the efficacy of the proposed method.

9 Numerical examples

We solve the two optimization problems using the strain energy transition model and the MMA-scheme with default parameters, cf. Svanberg (1987). The problems are discretized with 2D quadrilateral plane strain fully integrated elements. The constitutive parameters are the same as in the previous example. The 10 lowest eigenvalues, i.e. buckling loads are considered in the (26) constraint. The filter (5) is enabled and a continuation scheme for the penalty factor, q, is utilized wherein the initial value of q = 2 is increased to the final value q = 6 in increments of 0.05 every 5th design iteration. We chose a conservative continuation scheme of q as this can be used for all numerical examples. After the terminal value of q = 6 is reached, we wait 20 iteration before the Heaviside parameter β, cf. (6), is increased from β = 1 to β ≈ 12 following the exponential continuation scheme $\beta \rightarrow 1.3\cdot \beta$ every β_incr design iteration. The threshold parameter η = 0.5 in (6) remains fixed throughout the iterations.

9.1 Stubby cantilever example

The first optimization problem is the stubby cantilever beam design illustrated in Fig. 3 subject to a transverse load. The filter length scale is $l=40/2\sqrt {3}$ mm, cf. (5), and the maximum allowed volume is V_max = 0.2V. The operating load is initially λ^op = 20 kN, and increased by 20 kN each design iteration to the terminal load of λ^op = 140 kN. After the terminal load is reached, we wait three additional iterations before introducing the stability constraint (27). The reason for not initially activating the stability constraint is to limit the occurrence of artificial buckling modes associated with stress concentrations in the support and load regions, see Ferrari and Sigmund (2019) for further discussion. The β update increment parameter is β_incr = 10. The geometry is discretized using 90 × 200 = 18,000 elements. The load vector, P, is uniformly applied to the 11 vertical degrees-of-freedom symmetrically located about the right center point.

Fig. 3

Stubby cantilever beam design domain and boundary conditions. The angle $\alpha = \frac {\pi }{4}$ rad

To investigate the influence of the safety factor, we provide 10 examples with s in the range s ∈ [1.1,2.0], cf. (27). The optimized designs appear in Fig. 4. The results show that as the safety factor increases, i.e. minimum buckling load increases, the lower beam maintains its volume while increasing its second moment of inertia with respect to its neutral axis by developing a web-like structure.

Fig. 4

Optimized stubby beam designs for varying safety factor s

Figure 5 depicts the ten first buckling modes and Fig. 6 plots the histories of the ten load scaling factors π_j of the Fig. 5s = 2.0 design. It is clear that the strain energy transition scheme successfully removes artificial buckling modes from the optimization. The small buckling load oscillations that occur before iteration 400 are due to the q-continuation scheme, whereas the more prominent oscillations thereafter are associated with the β-continuation scheme.

Fig. 5

Scaled first ten buckling modes ϕ_j, $j \in \mathbb {N}_{\lambda }$ for the Fig. 4, s = 2.0, design

Fig. 6

The design histories of the ten load scaling factors π_j for the Fig. 4, s = 2.0, design

Fig. 7

Arch design domain and boundary conditions. The angle $\alpha = \frac {\pi }{6}$ and the radii R_in = 10000 mm, R_out = 10400 mm

9.2 Arch example

In this example, we design the arch depicted in Fig. 7. The filter length scale is $l=32/2\sqrt {3}$ mm, cf. (5), the maximum allowed volume is V_max = 0.4V and β_incr = 30. The operating load is initially λ^op = 1.5 kN, and increased by 0.2 kN each design iteration until λ^op = 2.7 kN. The artificial modes associated with stress localizations in load and support regions are not problematic in the arch compared to the stubby cantilever problem, and therefore the stability constraint is enabled after the first design iteration. The geometry is discretized using 666 × 50 elements. The structure is subject to the load defined by the load vector, P, which is uniformly applied to the 19 vertical degrees-of-freedom symmetrically located about the top center point.

We again investigate the influence of the safety factor s on the optimized design, cf. Fig. 8, where it is again apparent that as s increases so too does the second moment of inertia with respect to the neutral axes of the beams.

Fig. 8

Optimized arch designs for two safety factor s values

Figure 9 depicts the first 10 buckling modes of the Fig. 8, s = 2.0, design. Once again we emphasize the effectiveness of the strain energy transition scheme in removing artificial modes.

Fig. 9

Scaled first ten buckling modes ϕ_j, $j \in \mathbb {N}_{\lambda }$ for the Fig. 8, s = 2.0, design

Figure 10 displays the histories of the ten first load scaling factors π_j, $j \in \mathbb {N}_{\lambda }$ for the Fig. 8, s = 2.0, design. We again observe oscillatory behaviour due to the q and β continuation strategies, but otherwise a stable convergence.

Fig. 10

The design histories of the ten first load scaling factors π_j for the Fig. 8, s = 2.0, design

10 Conclusions

In this work, we demonstrate how the energy transition scheme presented by Wang et al. (2014) can be applied in buckling constrained topology optimization to eliminate artificial eigenmodes. When considering finite strain loadings, this scheme also improves convergence in the finite element simulations. The nondifferentiability of degenerate critical points is handled by p-norm aggregation, as done by Torii and De Faria (2017).

The numerical examples validate the efficacy of the energy transition scheme to remove artifical eigenmodes and produce stiff designs that are stable.

Appendix : A: Sensitivity identities

The first part of the right hand side of (38) is evaluated using (5) and the identity $\frac {\partial \left (\textsf {\textbf {B}} \pmb {\upphi }_{e,j}\right )}{\partial {\textsf {\textbf {a}}}_{e}} = \textsf {\textbf {B}}^{l}(\pmb {\upphi }_{e,j})$, such that

where the explicit evaluation of the last term in (41) appears in Wallin et al. (2018). The previously presented identity, and (5), are used to evaluate the second part of (38)

The combined result of (41) and (41) is denoted h for later use, i.e.

$$ \begin{array}{@{}rcl@{}} \gamma\textsf{\textbf{h}}(\gamma\textsf{\textbf{a}}_{e}, \gamma\pmb{\upphi}_{e,j}) &=& \gamma\textsf{\textbf{B}}(\gamma{\textsf{\textbf{a}}}_{e})^{T}\Big({\uppi}_{j}\textsf{\textbf{D}}\textsf{\textbf{B}}^{l}(\gamma\pmb{\upphi}_{e,j}) \gamma \pmb{\upphi}_{e,j} \\&&+ 2 \frac{\partial}{\partial \textsf{\textbf{C}}} \left[\gamma\hat{\pmb{\upphi}}_{e,j}^{T} \textsf{\textbf{B}}(\gamma \hat{{\textsf{\textbf{a}}}}_{e})^{T} \textsf{\textbf{D}} \textsf{\textbf{B}}(\gamma \hat{{\textsf{\textbf{a}}}}_{e})\gamma\hat{\pmb{\upphi}}_{e,j} \right] \Big) \\&& + 2\gamma\textsf{\textbf{B}}^{l}\left( \gamma\pmb{\upphi}_{e,j}\right)^{T}\textsf{\textbf{D}}\textsf{\textbf{B}}(\gamma{\textsf{\textbf{a}}}_{e})\gamma\pmb{\upphi}_{e,j}. \end{array} $$

(43)

The first part on the right hand side of (39) is evaluated as

using (7), (8), (9) and (??). The chain rule results $\frac {\partial \bar {\nu }}{\partial \nu }$ and $\frac {\partial \nu }{\partial \pmb {\upnu }}$ appear in, e.g. Wallin et al. (2018). Similarly, the second part on the right hand side of (39) is determined as