Two-scale optimization of graded lattice structures respecting buckling on micro- and macroscale

Interest in components with detailed structures increased with the progress in advanced manufacturing techniques in recent years. Parts with graded lattice elements can provide interesting mechanical, thermal, and acoustic properties compared to parts where only coarse features are included. One of these improvements is better global buckling resistance of the component. However, thin features are prone to local buckling. Normally, analyses with high-computational effort are conducted on high-resolution finite element meshes to optimize parts with good global and local stability. Until recently, works focused only on either global or local buckling behavior. We use two-scale optimization based on asymptotic homogenization of elastic properties and local buckling behavior to reduce the effort of full-scale analyses. For this, we present an approach for concurrent local and global buckling optimization of parameterized graded lattice structures. It is based on a worst-case model for the homogenized buckling load factor, which acts as a safeguard against pure local buckling. Cross-modes residing on both scales are not detected. We support our theory with numerical examples and validations on dehomogenized designs, which show the capabilities of our method, and discuss the advantages and limitations of the worst-case model.


Introduction
The ongoing progress in additive manufacturing allows structures with fine details to gain increasing focus.Lattice structures in particular, both homogeneous and graded, are utilized in many applications, e.g., thermal management, energy absorption, noise reduction, biomedical engineering, etc. (Rahman et al. 2022).Lattice infill is also recognized as potentially increasing global buckling resistance of a component (Clausen et al. 2016).However, fine features are prone to local buckling (Ferrari and Sigmund 2019).Two-scale optimization (Wu et al. 2021) can be used to design such structures without the need to resolve all the fine details of the full design in a single setting.The idea of this approach started with the work of Bendsøe and Kikuchi (1988), in which the design process is divided into two scales: the macroscopic scale, which describes the overall component, and the microscopic scale, which shows the fine details.Bendsøe and Kikuchi bridged the gap between these two scales by asymptotic homogenization (Allaire et al. 1997).This technique yields approximate material properties of microstructures on the macroscopic scale, which can then be used in macroscopic constitutive equations.Choosing a parameterized microstructure, Bendsøe and Kikuchi were able to effectively decouple the scales: Prior to any optimization procedure a discrete subset of the parameter space is chosen and homogenization is conducted for each of the microstructures gained from this parameter set.The obtained properties are then interpolated in the continuous parameter space and the interpolated material model can later be reused to solve various optimization problems.This technique introduces an interpolation error, but requires less computational effort when compared with on the fly homogenization, i.e., homogenization performed for each finite element in the discretized design domain and for each update step of the design during an iterative optimization procedure.Moreover the interpolation error can be controlled in a rather straight forward way during preprocessing.
Though there is exhaustive literature on optimal design considering the buckling behavior of structures using beam models (Ferrari and Sigmund (2019) and references therein), only a relatively small number of publications for continuum models exist.The initial problems evolved around finding optimal cross-sections for columns of fixed length and weight subject to uniaxial compression loads (Clausen 1851;Keller 1960;Tadjbakhsh and Keller 1962;Huang and Sheu 1968;Khot et al. 1976).Neves et al. (1995) were the first to conduct topology optimization with respect to buckling based on the method of Bendsøe and Kikuchi described above.However, they encountered various obstacles in the linear buckling analysis, which is stated as an eigenvalue problem (Bendsøe and Sigmund 2003).This includes localized modes in low density regions and clustering of eigenvalues when approaching the final design (Seyranian et al. 1994).Actions to alleviate issues with artificial, low-density modes have been proposed, e.g., different interpolation schemes for the stiffness and geometric stiffness matrices (Pedersen 2000;Bendsøe and Sigmund 2003) or applying an eigenvalue shift based on the last iteration and identifying artificial modes by their contribution to the total strain energy (Gao and Ma 2015).There are also methods for avoiding low-density regions like filtering and projection of the pseudo-density (Larsen et al. 2018), penalization of intermediate values in the objective (Allaire and Francfort 1993;Allaire and Kohn 1993), or element removal strategies (Behrou et al. 2021;Dalklint et al. 2020;Giele et al. 2021).Clustering of eigenvalues can be prevented by enforcing gaps between the eigenvalues (Bendsøe and Sigmund 2003).However, a large number of eigenvalues may still have to be computed to achieve good convergence (Bruyneel et al. 2008), which compels the use of efficient eigenproblem solvers (Dunning et al. 2016;Ferrari and Sigmund 2020).Nevertheless, topology optimization with respect to buckling still currently poses a challenging problem.
More recently, stability requirements have also been employed when tailoring microstructures (Neves et al. 2002b,a;Thomsen et al. 2018;Andersen et al. 2022), though models for the buckling of periodic microstructures have been investigated for decades.Homogenization theory for buckling load factors is well established (Neves et al. 2002b;Thomsen et al. 2018), but several challenges still arise in this context.Buckling modes can range from high-frequency modes with a wavelength shorter than the characteristic size of the microstructure to modes that span over multiple periods of the microstructure.Floquet-Bloch theory can be used to capture the latter in particular (Neves 2019).
For dehomogenized designs scale effects stemming from a finite cell size and effects due to grading of the microstructure might appear (Thomsen et al. 2018).Especially the latter might lead to an occurrence of undesired buckling modes at the boundary of the structure.
The aforementioned works only investigate the buckling behavior of a structure on a single scale, either macroscopic or microscopic.Even with the ongoing "competition for ultimately stiff and strong architected materials" (Andersen et al. 2021), these designed materials appear not to have been utilized so far in the optimization of components with stability requirements.
In this article, we present a two-scale topology optimization approach, where we incorporate buckling on both scales individually, but neglect cross-modes spanning over both scales.Assuming separation of length scales, we use homogenization theory to upscale the elastic behavior and buckling stability of a periodic lattice that resides on the microscopic scale.On the macroscopic scale we search for an optimal lattice grading with respect to both local and global buckling.We parameterize the lattice by its porosity, precompute homogenized properties for a selected set of porosities, and obtain a model on the continuous parameter space by piecewise cubic interpolation.This model is then used in macroscopic optimization problems with the local lattice porosity as design variable.Thus, our method can be classified as a multi-scale density optimization approach with a parameterized unit cell based on the microstructure porosity (category V/A in Wu et al. (2021)).
Our contribution is a new approach to integrate the upscaled microscopic buckling information into the optimization.From the macroscopic local stress, we obtain the buckling response on the microscale via a worst-case model, which acts as a safeguard against pure microscopic buckling.The worst case is realized by reducing the macroscopic local stress to its magnitude for reasons of computational effort.This comes at the cost of underestimation of predicted critical loads for the microstructure.
Recent work by Wang et al. (2021) is close to our approach.However, they deduce local stress constraints from the slenderness ratio of the lattice struts, while we directly use upscaled microscopic buckling information.Our approach has the advantage that essentially arbitrary microstructures can be treated, including those for which a slenderness ratio could not be uniquely defined.
The work of Christensen et al. (2022) is similar to ours, with some key differences: They assume isotropic buckling behavior of the microstructure, while our approach is applicable to arbitrary microstructures.They also fit a Willam-Warnke failure criterion to homogenized data, which introduces an approximation error of the local buckling factor.In contrast, we build a worst-case model from homogenized data, which also comes with an approximation error; however, applying a full interpolation model in lieu of the worst case, it is possible to render this error negligible (see Remark 1 in Section 4.1).Moreover, Christensen et al. (2022) employ a two term interpolation scheme to interpolate the buckling strength as a function of the relative porosity, whereas we use more accurate piecewise cubic interpolation.On the other hand, our approach is limited to sizing problems with a lower bound on the lattice density everywhere, while a major contribution of Christensen et al. (2022) is the possibility to allow for void regions in the design while still being able to maintain the lower bound on the density of the lattice using an auxiliary design field.
In this article, we restrict the analyses to linear elasticity and linearized buckling, though the general idea of the method can be extended to non-linear regimes.We also limit ourselves to a two-dimensional setting in this article to keep notations simple.However, it is straightforward to extend our approach to three dimensions.We emphasize that our method is applicable to arbitrary, parameterized microstructures.As an example, we choose a lattice consisting of equilateral triangles that is parameterized by its porosity.More sophisticated geometries and parameterizations are possible with higher computational effort without changing the general method.
The remainder of the article is structured as follows: In Section 2 we recap the state equations for linear elasticity and linearized buckling analysis.Section 3 briefly presents the main formulas for asymptotic homogenization of elasticity modulus and buckling load factors.We describe our exemplary microstructure, its parameterization, and important aspects for the upscaling process in Section 4. We present our method for designing a lattice unit cell by rounding corners of an equilateral triangle, showing finite element convergence on the microscopic scale, and obtaining homogenized properties of a lattice.Section 4.1 provides the novel approach to include the upscaled properties on the macroscopic scale via the suggested worst-case model for the homogenized buckling load factor.Section 4.2 demonstrates how to precompute homogenized properties on a discretized parameter space and apply an interpolation scheme to yield values on the continuous parameter space.Section 5 outlines three different two-scale optimization problems considering buckling on the macroscopic scale, on the microscopic scale and on both scales (without cross-modes).
Associated numerical examples are presented in Section 6 together with a pre-study, which helps to develop a better understanding of the optimized designs.A numerical validation on dehomogenized designs is given in Section 7.1 accompanied by a study on scale separation.We examine the impact of the worst-case error closely in Section 7.2 by means of two numerical examples.Finally, we complete with conclusions in Section 8.

Linear elasticity and buckling analysis
In this section we briefly recap the state equations for linear elasticity and buckling analysis on an algebraic level as given by Thomsen et al. (2018).
For continuous formulations in weak form, interested readers are referred to Neves (2019).We assume a linear elasticity setting and linear bifurcation buckling condition.This means that the prebuckling displacement, stresses, and strains vary linearly with the applied load and that the load factor, which indicates stability, appears linearly in the bifurcation eigenvalue problem.Linear buckling analysis consists of two steps: First, we solve the linear elasticity state equation for displacements of a structure under a given load.
Then, an eigenvalue problem is solved, where the eigenvalues correspond to bifurcation points, and eigenvectors are interpreted as buckling modes of deflection.
Throughout this article, we apply Voigt notation for tensors and assume plane stress conditions.
Let us consider a body Ω ⊂ R 2 and its discretization Ω = M e=1 Ω e by M finite elements.Then the state equation of linear elasticity for Ω reads as (Zienkiewicz et al. 2005): where u is the sought-for vector of displacements, ρ ∈ (0, 1] M is the pseudo-density field and f is the applied load vector (reference load).The stiffness matrix is given by with assembly operator denoted by and number of integration points N ip .E(ρ e ) is the elasticity tensor, which depends on the pseudo-density in element e. B e,k is the strain-displacement matrix of element e evaluated in the k-th integration point and contains derivatives of the finite element's shape functions.The factor c e,k = det(J e )w e,k contains the Jacobian determinant det(J e ) of element e and the integration weight w e,k associated with the k-th integration point of said element.
The buckling equation is given by the eigenvalue problem, and is solved for pairs of eigenvalues Λ and eigenvectors φ , ≤ M. Eigenvalues of (4) are also called load factors.The critical load, under which buckling occurs, can be computed by multiplying the reference load f with the eigenvalue with the smallest absolute value, which is also known as the critical load factor or buckling load factor (BLF).If the critical load factor is less than or equal to one, the structure buckles; otherwise, the structure is stable with respect to the applied reference load.The stress stiffness matrix (or geometric stiffness matrix) (Němec et al. 2016) in ( 4) is given by c e,k Be,i,k (σ e (u e (ρ))) i,j Be,j,k , (5) where summation over i and j is implied.Again, c e,k includes the Jacobian determinant of element e and the integration weight at the k-th integration point.The integrated stress in element e is numerically evaluated as and Be,i,k = ∂Ne ∂xi (x k ) (7) is the derivative of the matrix of the shape functions evaluated at integration point x k .Usually, solutions to (4) are considered ordered according to the absolute value of the load factors with Λ 1 as the smallest value and eigenvectors are G-normalized, i.e., φ G φ = 1. (8) In single-scale topology optimization, (6) can lead to artificial buckling modes in low density regions (Neves et al. 1995).However, in our approach the macroscopic pseudo-density represents the local relative lattice volume.Thus, low density regions represent thin lattice, and modes in such regions are not artificial, but rather comply with buckling of the lattice.Following Rodrigues et al. (1995), the sensitivity of an eigenvalue with algebraic multiplicity one with respect to density ρ e is given as (9) Here, w is the adjoint state obtained from the solution of the equation For eigenvalues with multiplicity greater one, the derivative of the load factor function does not exist in a strict sense.Nevertheless, using (9) still an element of the subdifferential can be computed (Rodrigues et al. 1995).Alternatively, derivations for multiple eigenvalues can be used (Seyranian et al. 1994).Note that G is in general not positive definite, so sometimes a reformulation of the eigenvalue problem (4) to G(ρ, u(ρ)) − 1 Λ K(ρ) φ = 0, (11) might be beneficial, as K is always positive definite for ρ > 0.

Asymptotic homogenization
To obtain the mechanical properties of a lattice structure on the macroscopic scale, we perform asymptotic homogenization.We briefly recall the homogenization formulas in the two-dimensional setting, which can be found, e.g., in Thomsen et al. (2018).Note that we do not perform topology optimization on the microscopic scale, but rather investigate given periodic lattice structures.Thus, there is only a finite element mesh for the solid material region.This prevents spurious modes, which typically arise in low density regions during topology optimization (Neves et al. 1995).However, we might get artificial modes owing to the finite element discretization.These modes are usually highly localized and can be filtered by setting a threshold on the minimal number of nodes that have to exhibit deflection, e.g., 5% of all nodes in the finite element mesh.
As we do not perform topology optimization on the microscopic scale and are only interested in obtaining the critical load factor of a given microstructure, sensitivities are not required here and thus clustering of eigenvalues as described in the work by Neves et al. (2002a) is not an issue.

Linear elasticity
We choose a representative volume element (RVE) Y , which is discretized with m finite elements Y e , i.e., Y = m e=1 Y e .We solve the equilibrium equations on the microscopic scale for the Y -periodic fields χ j using three test fields f j : with stiffness matrix K and f j given by and n ip is the number of integration points, E 0 is the elasticity tensor for the base material and ˜ j are unit strain fields with ˜ j i = δ ij .We note that the Y -periodicity is realized by assuming same values for the solutions on opposite sides of Y .The homogenized constitutive tensor E H of the microstructure is then given by (15)

Buckling
The equilibrium equation for buckling on the microscopic scale is an eigenvalue problem as in the macroscopic setting (cf.( 4)) and is solved for pairs of eigenvalues λ and associated eigenvectors ϕ , ≤ m: As in the macroscopic setting, the eigenvalue with the smallest absolute value denotes the buckling load factor.G is the microscopic initial stress stiffness matrix (cf.( 5)) with The microscopic initial stress tensor σ e describes how the macroscopic strain ¯ is distributed in unit cell Y .We note that (18) already takes stress amplification (see Ferrer et al. ( 2021)) into account.The macroscopic strain is given as solution of the macroscopic constitutive equation and X e = [χ 1 e χ 2 e χ 3 e ] contains the solutions of (12).The macroscopic stress σ is obtained via (1) and ( 6).Parameterization of the lattice geometry and linearity of the buckling analysis allows us to precompute homogenized properties.Thus, we can avoid solving homogenization problems during optimization (cf.Section 4.2).
We want to emphasize that the buckling analysis in (4) and ( 16) yields pure global and pure local modes, respectively.In particular, cross-modes spanning over both scales and defined as a solution to Eq. ( 40) in the work of Neves (2019), are not detected.

Microstructure parameterization and upscaling model
In the following section, we present our method to integrate the microscopic buckling load factor into a macroscopic scale optimization problem.The explanation is done for an exemplary lattice structure; however, we would like to stress again that the method can be applied to arbitrary base cell topologies.We first demonstrate an approach to obtain the homogenized load factors of the lattice.
The exemplary lattice we chose consists of equilateral triangles with edge length d, as the principal design of our microstructure.The primitive lattice vectors are thus given by R 1 = (1, 0) and R 2 = (1/2, (3)/2) (Fig. 1).This design is known in literature to provide good macroscopic buckling strength (Clausen et al. 2016).To get a periodic unit cell, we take two of these triangles, which together form a parallelogram, including its shorter diagonal (see Fig. 1).We parameterize the unit cell by one parameter ρ, which describes the relative volume or one minus porosity, respectively.With this parameterization the homogenization procedures from Section 3 can be written using two maps: (20) assigns each relative volume ρ a symmetric, positive definite homogenized elasticity tensor computed from (15), while ( 21) maps a macroscopic stress and a relative volume to the resulting homogenized buckling load factor λ(σ, ρ) obtained by solving (16).
As material parameters, we chose a Young's modulus of 10 Pa and a Poisson's ratio of 0.3.
Next, we present three aspects that are important for the upscaling process: the design of joints of lattice struts, the resolution of the finite element mesh on the microscopic level, and the number of unit cells in an RVE.We note that the latter is automatically covered if a Floquet-Bloch approach along with a sufficiently fine discretization of the Brillouin zone is used.

Rounded corners
A unit cell design with sharp corners where the lattice struts meet leads to local stress concentrations at those corners.To circumvent this issue, we round the corners with circular arcs with radius r, while keeping the overall volume of the structure constant (Fig. 2).A lattice with rounded unit cells can be seen in Fig. 8.It turns out that the choice of the radius parameterizing the arcs has only minor effect on the Poisson's ratio (Fig. 5).The Young's modulus is affected mainly for small relative lattice volume (Fig. 4); for a small volume, a large radius leads to very thin lattice struts due to volume preservation (cf.Fig. 2) and the lattice loses stiffness.The impact on the homogenized buckling load factor is more significant as can be seen in Fig. 3.For small radii, the buckling resistance increases with increasing radius, because stress concentrations are avoided, lattice struts get better supported, and their relative length to width ratio gets smaller.For larger radii, the struts get thinner, their length to width ratio increases, and the load factor decreases.For different volumes, the maximum of the smallest buckling load factor is achieved at different radii.The data in Fig. 3 suggests that rounding corners with a radius depending on the volume leads to good buckling resistance.We note that we only investigated an uniaxial loading with σ) = (−1, 0, 0), and the optimal radius might be different for other stress situations.To keep things simple, we thus decided to use a constant radius of r = 0.05d for our subsequent calculations.

Finite element convergence on the microscale
We discretize the microstructure by triangles with second order shape functions.A convergence study with respect to the number of finite elements for homogenization, i.e., on the microscopic scale, can be seen in Figs. 6 and 7. Observing that the convergence graphs are already very flat when approaching one million elements and taking into account that we want to keep the finite element error in the preprocessing small, we opt to choose for all subsequent homogenization procedures a discretization, which results for a density of ρ = 0.5 in approximately one million finite elements per RVE.Geometries corresponding to lower densities are resolved using slightly less elements.Given this number of finite elements, a single analysis is usually completed within a few minutes on a standard work station.

Choice of RVE size
Asymptotic homogenization can only detect highfrequency modes, i.e., modes with a wave length that is smaller than the size of the representative volume element.However, buckling modes  might span over more than one cell.To identify these modes also in a homogeneous lattice structure, we encompass more than one unit cell in the RVE.That is, the microstructure is considered Y κ -periodic, where κ 2 is the number of unit cells inside the RVE.RVEs containing various unit cells can be seen in Fig. 8.It is noted that all modes of an RVE with a given number κ of cell repetitions can also be found using an RVE comprised of κ × κ unit cells, where is an integer number, e.g., all modes of a 2 × 2 RVE appear also in the analysis of a 4 × 4 or 6 × 6 RVE.Using an analogous argument, modes with periodicity κ n1 along R 1 and periodicity κ n2 along R 2 for n 1 , n 2 ∈ N are captured using a κ × κ RVE.Note that in practice there is a natural lower bound; this value is given by the resolution of the underlying finite element discretization.
We observe that especially an RVE with only one unit cell results in quite high homogenized load factors.This is because the joint of the struts gets locked in rotation due to the periodicity of the RVE.In an RVE with more than one cell and in the analysis of dehomogenized designs in Section 7.1, joints of lattice structures are free to rotate (see Fig. 24).
As we are interested in the smallest, i.e., critical, load factor of the lattice, we pick the minimal load factor with respect to cell repetitions: In practice, one is not able to perform homogenization with κ ∈ N cell repetitions, but has to resort to κ ≤ K with an appropriately chosen upper bound K on the number of cell repetitions.We would like to stress that for a given cell layout, the optimal κ ∈ {1, 2, . . ., K} cannot, in general, be determined without testing essentially all choices and might even vary for fixed stress σ but different volume fractions ρ.
We want to remark that the presented method of repeating the unit cell within the RVE is analogous to applying Floquet-Bloch theory with a special discretization of the Brillouin zone (see Fig. 9).In Floquet-Bloch theory, periodicity conditions for modes φ are given as where i refers to the complex unity, k is the wave vector, and y is a location inside a lattice unit cell (Thomsen et al. 2018).Rather than solving cell problems with boundary conditions ( 23) for all k ∈ R 2 , in practice a number of k vectors on the boundary of the so called irreducible Brillouin zone (IBZ) is selected.This corresponds to a discretization of the IBZ.Now, a ( κ n1 , κ n2 )-periodic mode with n 1 , n 2 ≤ κ corresponds to a wave vector, which solves the system With that, we can plot all wave vectors corresponding to modes, covered by our κ × κ RVEs with κ from 1 to K, into a Floquet-Bloch diagram with an outline of the IBZ.Doing so for K = 7, we obtain the result depicted in Fig. 9.We emphasize that the worst-case model presented in the next Section 4.1 is independent of the method that is used to obtain the homogenized buckling load factor.

Worst-case model
This subsection describes our novel method to integrate the microscopic buckling load factor on the macroscopic scale.We stress that it is valid not only for our exemplary lattice but for any arbitrary, parameterized microstructure.
Due to linearized buckling analysis (( 16) to ( 19)) the homogenized load factor depends linearly on the macroscopic stress: That is, homogenization can be conducted with macroscopic unit stress (in some given norm) and later, the homogenized load factor has to be divided by σ .It is thus sufficient to examine stresses on the unit sphere surface S 2 instead of the whole three-dimensional stress space (σ xx , σ yy , σ xy ).We define which assigns a homogenized buckling load factor to each unit stress combined with a relative volume.
We use a spherical coordinate system to characterize unit stresses on S 2 : the zenith reference (z-axis) is the axis from the origin through the biaxial compression stress σ = (−1, −1, 0) and the azimuth reference (x-axis) is the axis from the origin through σ = (1, −1, 0) (Fig. 10).In this coordinate system, biaxial compression and tension stress conform to north and south pole, while all pure shear stresses rest on the equator.The inclination (or latitude if thinking of geographical coordinates) characterizes the type of stress: the rotation invariant biaxial compression and tension stresses conform to poles, and other special types, e.g., uniaxial and shear stresses form circles on the unit sphere surface (circle of latitude).The azimuthal angle (longitude) describes the rotation of the applied macroscopic stress σ relative to the RVE.For an additional explanation please refer to the video in Online Resource 1.
Common literature reduces the stress space even further.Under the assumption of isotropic buckling behavior, only biaxial loading without shearing component but varying the σxx σyy -ratio is investigated, or uniaxial loading with all possible load orientations relative to the investigated microstructure is applied (Bluhm et al. 2020).In contrast, we compute the buckling yield surface for the whole unit stress surface S 2 .Thus, our method is applicable to arbitrary, parameterized microstructures with isotropic or anisotropic buckling properties.
The homogenized load factors for uniaxial compression, biaxial compression and shear stress for an RVE with a relative volume of 30% and three cell repetitions are shown in Fig. 11.We can clearly see that the symmetry of the unit lattice cell is reflected in the load factors.The shape of the uniaxial loading case matches very well with results by Bluhm et al. (2020).
Next, we coalesce all data in a worst-case model.Having evaluated homogenized load factors for different cell repetitions κ ∈ N and all stresses on the unit stress sphere σ ∈ S 2 , i.e., for all stress types and directions, we select the smallest buckling load factor with respect to all unit stresses and number of cell repetitions: This worst-case model depends only on the local volume fraction and no longer on the local stress type or direction.We note that in practice the unit sphere is discretized and the number of cell repetitions bounded from above.The worst case is thus only a worst case with respect to the discretization resolution and maximal cell repetitions.We found, however, that the homogenized buckling load factors for our exemplary lattice cell have high regularity with respect to the stress variable.
If homogenized load factors were obtained via Floquet-Bloch theory (Neves et al. 2002a) instead of using different cell repetitions, (27) would contain a minimization over all possible wave vectors k in the Brillouin zone B in lieu of cell repetitions:

Decoupling of micro-and macroscale
The parameterization of the unit cell allows for effective decoupling of the micro-and macroscopic scales (Bendsøe and Kikuchi 1988).That is, we discretize the parameter space (0, 1] for the relative volume and precompute homogenized properties (elasticity tensor from (15) and microscopic buckling load factor from ( 16)) for the resulting parameter set (see also Section 6).For optimization of structures on the macroscopic scale, we then apply an interpolation model to these precomputed properties.In other words, we replace the maps E H (20) and λ (21) by where E I and λ I are approximations of E H and λ, respectively.
To construct λ I , we interpolate the worst-case load factors obtained from ( 27) with respect to the density variable: We obtain the worst-case microscopic buckling load factor associated with macroscopic stress σ by dividing by the norm of this stress (see (25)): For gradient based optimization, we need a differentiable interpolation model.Following Bendsøe and Kikuchi (1988) we apply a piecewise interpolation strategy.More precisely, we employ piecewise cubic Hermite interpolation (Birkhoff et al. 1968) for both, the homogenized elasticity tensor (29) as well as the buckling load factor (31).This scheme yields a continuously differentiable approximation that is composed of uniquely defined cubic polynomials between the provided data points and is exact at these points.However, other interpolation schemes using, e.g., tangent or RAMP (Rational Approximation of Material Properties, Stolpe and Svanberg ( 2001)) functionals are possible with the downside of lower accuracy, especially for high porosity values.
Under the piecewise Hermite approach, E I is constructed by individual interpolation of each coefficient of the elasticity tensor.The first order derivatives of the homogenized elasticity tensor and homogenized buckling load factor with respect to ρ are approximated by finite differences.For the homogenized load factors, we get numerical issues for high relative volumes.For these volumes, the unit cell is very resistant to buckling, but the iterative eigenvalue problem solver (ARPACK, version 3.7.0(Lehoucq et al. 1998)) yields artificial modes, which are caused by the buckling of individual finite elements.Hence, we only interpolate data points below 60% relative volume.At the boundary of this sample space, second-order central finite differences are not available, and thus only quadratic polynomials are used in the two outer subintervals, i.e., [0, 0.05] and [0.55, 0.6].For relative volumes above the chosen threshold of 60% we extrapolate using the quadratic function obtained for the subinterval [0.55, 0.6].It is noted, however, that in the optimized designs presented in Section 6 microbuckling occurs only for significantly lower densities of 60%.The interpolated worst-case model λ I wc for the homogenized microscopic buckling load factor is shown in Fig. 12.The presented interpolation scheme leads to sufficiently differentiable functions to perform gradient-based optimization.
The error of the approximated microscopic buckling load factor λ I comprises several individual errors: the error introduced by finite element discretization on the microscopic level to solve the homogenization problems, the error arising from restriction to the worst case, and the error resulting from the interpolation.The discretization and interpolation errors can easily be controlled by using finer finite element meshes on the microscopic scale and finer interpolation grids.Thus, the worst-case error has the highest significance.We will come back to this when we discuss numerical examples in Section 7.2.
Remark 1 The worst-case error can be avoided if the worst-case model ( 27) and univariate interpolation (31) are replaced by a discretization and trivariate interpolation of (26) in the combined stress and density space S 2 × (0, 1]. In gradient-based optimization context, a continuously differentiable interpolation scheme is essential.Possible techniques to realize that for this worst-case free approach include piecewise cubic Hermite interpolation (Birkhoff et al. 1968) or interpolation on sparse grids based on cubic B-splines (Valentin et al. 2020).When applying the worst-case free approach in a three-dimensional setting, due to the curse of dimensionality, the differentiable sparse grid approach is preferable, as stress has six independent entries there, which result in five parameters to represent unit stress, requiring interpolation to be carried out in a six dimensional space.

Sizing optimization
In this section, we formulate two-scale sizing optimization problems, which will be solved in Section 6.In contrast to topology optimization, where usually a solid/void design is of interest, we vary the local lattice volume fraction between 10% and 100%.This corresponds to a scaling of the width of the lattice struts.
We want to achieve structures that are resistant to given loadings both with respect to their stiffness and their buckling strength.To maximize buckling strength in our two-scale approach, we have to take into account buckling of the homogenized overall component as well as buckling of the microstructure.This leads to a multiobjective optimization problem.For a macroscopic domain, discretized by M finite elements, it can be stated that in terms of mechanical compliance c, macroscopic load factors Λ , = 1, . . ., L and microscopic buckling load factors λ I e , e = 1, . . ., M: The macroscopic load factors are given as solutions of ( 4) and a prediction for the microscopic buckling load factor λ I e for each element e is obtained from the interpolated worst-case model (32).The design variable ρ represents the local volume fraction of the homogenized lattice.Via the homogenization formulas, this can be mapped to mechanical properties of a (periodic) microstructure, e.g., elasticity tensor (15) and buckling load factor (32), which uses (16).g is an aggregating function and U ad is the admissible set (34) We include the first few L M macroscopic load factors in the objective to handle potential mode switching and multiple eigenvalues.Note that minimization of the compliance corresponds to maximization of its negative value.We treat the multi-objective problem (33) with the -constraint method (Mavrotas 2009).In this, a point on the Pareto front is obtained by solving the following problem: with an imposed compliance value c 0 .We want to maximize the smallest load factor to achieve a good buckling strength and thus choose g to be the min function.To get rid of the non-smooth character of the latter, we apply a bound formulation of this problem, as suggested by Bendsøe and Sigmund (2003), and end up with the following optimization problem: max s∈R,ρ∈U ad s, (36) λ I e (ρ) ≥ s, e = 1, . . ., M ( 38) With the bound formulation we experienced no problems with multiple eigenvalues or clustering of eigenvalues in our numerical experiments, if L was chosen sufficiently large.Alternatively, other formulations like a smooth minimum via Kreisselmeier and Steinhauser (1980) could be considered.
In Section 6.2, we look at three variations of this problem: A) We ignore (38), i.e., we maximize buckling on the macroscopic scale with a compliance constraint but do not take buckling on the microscopic scale into account during optimization.B) We drop (37), i.e., we only optimize microscopic buckling under a compliance constraint but disregard macroscopic buckling.C) We investigate the problem as given, i.e., we maximize buckling on both scales with respect to a given compliance.We compute a representation of the Pareto optimal set with 51 Pareto optimized solutions for each case.For this, we first compute a solution to the pure compliance minimization problem c * = arg min The representation of the Pareto optimal set is then given by solutions to A) -C) with compliance constraint c(ρ) = c j 0 in (39).

Numerical examples
In this section we present the results of numerical experiments for the two-scale sizing problems stated in Section 5. We first conduct a pre-study to develop a better understanding of later optimized designs.After that, solutions to the multi-objective optimization problems will be investigated.
To construct the worst-case model, we have to conduct buckling homogenization with applied macroscopic stress.For this, we discretize the unit stress sphere in a spherical coordinate system (see Section 4.1).We discretize the azimuthal angle, which describes the rotation of the applied stress around the RVE, with 5 • steps.Exploiting the symmetry of our unit cell, it suffices to investigate the interval from 0 • to 30 • degrees.For the inclination, which describes the type of stress, we choose 11.25 • steps in the interval [0 • , 180 • ).We discretize the relative volume with 0.05 steps and solve homogenization problems for four, five, six, and seven cell repetitions.Modes obtained with an RVE consisting of two and three unit cells can also be detected in an RVE with six unit cell repetitions.Likewise, the modes for an 1x1 RVE are covered by all other RVEs (see Section 4).This leads to 448 cell problems for each chosen relative volume and 9408 in total.Note that these simulations are independent of each other and can be run in parallel.We discretize the microscopic domain by approximately one million finite elements per RVE.As already mentioned in Section 4, the material properties are given by a Young's modulus of 10 Pa and a Poisson's ratio of 0.3.It is noted that the motivation for the choice of the Young's modulus was an improved numerical behaviour, in particular when solving the associated eigenvalue problem on the microscopic scale.Conversely, the results are invariant with respect to this choice.
The eigenvalue problems are solved by ARPACK, version 3.7.0(Lehoucq et al. 1998).To solve optimization problems, we apply SNOPT, version 7.2.8 (Gill et al. 2005), which employs a sequential quadratic programming method.Now, consider the setting shown in Fig. 13, left.A rectangular design domain is subject to a 1 relative local volume ρ pressure load at the top.At the bottom rolling boundary conditions are applied, i.e., the degree of freedom in horizontal direction is fixed for all nodes at the bottom edge.To prevent rigid body movement, the degree of freedom in vertical direction is also fixed for the central node at the bottom edge.This corresponds to Euler's case of a slender column with fixed/free boundary conditions, for which the buckling load is given by where L is the length of the column, E is the elastic modulus and I is the planar second moment of area.The design domain has a ratio of 1 : 5.2 and we discretize it with 25 × 130 bi-linear quadrilateral elements (Q 4 ), which is a sufficient resolution for this type of element (compare also (Ferrari and Sigmund 2019, Fig. 1)).Fig. 13 shows the optimized design of a mechanical compliance minimization with a SIMP (Solid Isotropic Material with Penalization) material model and a lower physical design bound of 1e−9.

Pre-study: Endoskeleton versus exoskeleton
Before we look into results of the optimization problems, we perform a study to develop a better understanding of later optimized designs.As a first step, consider ρ = 1 on Ω, i.e., solid material everywhere (Fig. 14, left).Assuming plane stress conditions, a width and virtual depth of 1 m each, length 5.2 m and E = 10 Pa, we get f crit = 0.076 N from ( 42).This fits well to our numerical result of f crit = 0.074 N. Next, we compare this solid column with two structures that have the same weight each as the solid column: weak material (ρ = 0.4) with either a solid reinforcement (ρ = 1) of width 0.5 m centered in the middle (endoskeleton, Fig. 14, mid) or a solid reinforcement that flanks the weak material at both sides with width 0.25 m (exoskeleton, Fig. 14, right).The force is applied at the top boundary via a solid plate with an elastic modulus that is 100 times higher than the solid material.As the weight is kept constant, these two structures are wider than the first one (their width is 1.75 m).The critical loads are 0.065 N and 0.22 N, respectively.Due to a larger second moment of inertia, the exoskeleton shows superior macroscopic buckling stiffness compared to the endoskeleton.Thus, we expect designs that are optimized with respect to macroscopic buckling to exhibit an exoskeleton-like structure.

Optimization
Next, we investigate the three different multiobjective optimization problems as given in Section 5. Let us briefly note that we did not run into the problem of switching eigenvalues in any of the following examples.The material models are given by the interpolation functions from Section 4.1 and a density filter (Borrvall and Petersson 2001) is applied with a radius of 1.6 times the edge length of a finite element for regularization purpose.We choose ρ min = 10% to obtain easily realizable lattice structures and apply a global volume constraint with V = 50% of the design domain's area.

A) Pure macroscopic load factor optimization
Ignoring buckling on the microscopic scale, we maximize the macroscopic load factor (37).The load factors for different values of the compliance constraint can be seen in Fig. 15.As reference, we include the value of a topology optimization with power law (ρ 3 ) and ρ min = 0.001.Selected designs are shown in Fig. 16.Design A 1 results from a compliance minimization without considering buckling.Relaxing the compliance constraint while maximizing the buckling load factor, we only see a rather small improvement of up to 15% in the macroscopic load factor.This is due to design A 1 having already a rather good macroscopic buckling resistance, as it forms an exoskeleton.Only in the upper part of the design domain, where the loaded edge has to be supported, is this alleviated in favor of a typical branching structure.With a larger compliance, bound diagonal bars appear (A 2 -A 4 ), which are known to increase macroscopic buckling stability (Bendsøe and Sigmund 2003).
The lower design bound of ρ min = 10% is active for all designs, see Fig. 15.We conclude from this, that the load factor could potentially be larger, i.e., better structures could be obtained, if we chose a lower ρ min .On the other hand, for very low ρ min , the problem would no longer be a sizing but rather a topology optimization problem, which would require special handling (see, e.g., Christensen et al. (2022)).

B) Pure microscopic load factor optimization
In the second optimization problem, we ignore macroscopic buckling and instead maximize the smallest microscopic load factor of all finite elements (38) subject to a compliance constraint.The microscopic buckling load factors are approximated by our worst-case model (32).In Fig. 17, a substantial improvement of the load factor can be seen when the compliance constraint is relaxed.
The optimal solution for the problem without compliance constraint (design B 5 ) is a fully homogeneous design (see Fig. 18).Design B 5 exhibits homogeneous stress under the given boundary conditions, which results in a homogeneous microscopic load factor.Thus, the further we relax the compliance constraint, the more homogeneous the C) Simultaneous optimization of pure macro-and microscopic load factors The obtained values for the simultaneous optimization of both macroscopic and microscopic buckling load factors ( 36)-( 39) can be seen in Fig. 19.As reference, the values of a topology optimization (TO) with ρ = 0.001 are given.To achieve a stiff (small compliance) design, thick solid structures are needed.Little material is left for the lattice part, which results in low local volume and a small minimal microscopic buckling load factor.Hence, only ( 38) and ( 39) are active and (37) remains inactive.For less restrictive (larger) compliance bounds, less solid material is necessary and material is redistributed to the lattice region (C 2 ).Thus, the microscopic buckling load factor can be improved.When it reaches the value of the macroscopic one, (37) becomes active (C 3 ).Therefore, raising only the microscopic buckling load factor further, i.e., steering towards homogeneous design as in B, yields no improvement in the objective, as the macroscopic load factor will define the value of the slack variable.For this reason, both micro-and macroscopic buckling load factors are raised simultaneously: the first by increasing the lattice density especially in the upper region of the design domain and the second by stiffening the exoskeleton (C 4 ,C 5 ).
The decreasing minimal local volume in Fig. 19 for more relaxed compliance appears nonintuitive.In Fig. 20, staircase structures on the solid parts can be seen due to discretization, see, e.g., C 5 .The resulting local stress field allows a slightly lower volume for individual finite elements, while preserving the microscopic buckling load factor.

Dehomogenization and validation
In this section, we want to give two distinct comparisons: First, we compare the predicted buckling behavior with high-resolution numerical analyses of dehomogenized designs.Second, we validate the worst-case model against a precise microscopic load factor evaluation.No diagonal bars appear as they do in the pure macroscopic buckling optimization, since the lattice already provides good macroscopic buckling resistance.Due to discretization, steps can be seen on the solid parts, as marked in C 5 by white ellipses.

Performance of dehomogenized designs
The number of lattice cells for dehomogenization of optimized results is chosen independently of the finite element resolution of the macroscopic model (cf.Fig. 13).To realize this, we proceed as  follows: First, the optimized density field ρ is interpolated and the number of lattice cells is chosen.
Then, an evaluation of this field at each lattice cell's center defines the width of the corresponding lattice struts.Note that our dehomogenization procedure is volume preserving, i.e., all dehomogenized designs have a volume of V = 50% of the design domain's area.Dehomogenized designs are shown in Figs.22 and 23.

Influence of cell size
We compare dehomogenizations with different cell sizes for design C 4 from Fig. 20.We investigate different buckling load factors: the one associated with a macroscopic buckling mode, i.e., deflection of the structure as a whole (low-frequency), and the one associated with a microscopic mode, i.e., deflection of the lattice (high-frequency).As expected, the predicted macroscopic buckling load factor is met better with a higher number of lattice cells (see Fig. 21).The microscopic buckling load factor stems from modes at the boundary (Fig. 23).These modes cannot be approximated well by homogenization theory, which ignores boundary effects.We therefore also search for the smallest load factor associated with an interior mode, i.e., a mode, that does not exhibit deflection at the structure's boundary (cf.Fig. 23).In order to avoid buckling of the microstructure occurring at the structural boundary, a coating strategy as where the pressure load is applied exhibit sag, which is reduced for a higher number of cells.
suggested by Christensen et al. ( 2022) may be pursued.
In the considered example, the cell size has almost no influence on this interior load factor.This might differ for other examples.Zooms for interior microscopic buckling modes are given in Fig. 24.
As can be seen in Fig. 21, the compliance tends to decrease with smaller cell size.The observed increase for 22 lattice cells compared to 20 is presumably caused by discretization effects.Fig. 22 shows that the displacement under the applied pressure load has a wave shape, where horizontal lattice struts at the top face sag between their diagonal supports.The more cells we have, the better these struts are supported, and thus the compliance is reduced.

Comparison between homogenized and dehomogenized designs
Next, we compare the homogenized design C 4 with its dehomogenized counterpart with 20 cells in a horizontal direction.For this design, we have a clear separation of length scales and only a small homogenization error.The macroscopic load factor is determined to be 0.0507 and differs by less  than 1% from the predicted macroscopic load factor of 0.0511 in the optimized homogenized design.A gap of 10% exists between the interior microscopic load factor of the dehomogenized design of 0.0563 and the predicted one of 0.0511.This is expected, as the predicted microscopic load factor is based on our worst-case model, which assumes the worst stress type and orientation (see Section 4.1).Note that the predicted microscopic buckling load factor is smaller than the microscopic buckling load factor, and hence our worst-case model acts as a safeguard against pure microscopic buckling in this example.Of course, such an observation is generally only valid up to remaining discretization and interpolation errors as well as the error introduced by homogenization.

Impact of worst-case model
Next, we investigate the gap between the load factor predicted by the worst-case model and the interior load factor obtained from the dehomogenized design.For all finite elements in the lattice region of design C 4 , we perform a posteriori homogenization by evaluation of (26).We extract the local relative volume and create RVEs with κ ∈ {4, 5, 6, 7} unit cell repetitions.Then, we conduct homogenization (Section 3) of the buckling load factor for each κ with the "real" macroscopic stress taken from the optimized design using ( 18) and ( 19) and use the minimal value with respect to κ.The resulting microscopic load factors are shown in Fig. 25.The mean relative difference between predicted and a posteriori load factor for all of the validated finite elements is 18%, with a standard deviation of 0.079.The marked element in Fig. 25 complies with the center of the interior microscopic mode in the dehomogenized design.For this element, the a posteriori homogenized buckling load factor is 0.0578, which differs less than 3% from the interior microscopic load factor of the dehomogenized design.This rather small difference demonstrates that homogenization is a valid tool for predicting microscopic buckling behavior.
The a posteriori homogenized buckling load factor for the mentioned element is 13% larger than predicted by the worst-case model.This can be explained by Fig. 11.The macroscopic stress for the chosen finite element, which is extracted from the macroscopic analysis, is uniaxial and acts in vertical direction.Thus, the a posteriori load factor corresponds to the point on the uniaxial curve at 90 • , which is 15% larger than the biaxial value.

Example with tension, compression, and shear stress
In the previously considered example, we mostly observe uniaxial compression stress.Next, we want to examine an example with compression, tension, and shear stress.For this, consider the setting shown in Fig. 26: A bow-tie shaped domain is fixed at the left edge in both degrees of freedom.At the right edge, movement in a horizontal direction is prevented and a force is pulling downwards.We simulate the behavior of a homogenized lattice with a constant relative local volume, i.e., In Fig. 27 we see, that the microscopic buckling load factor, which is obtained from our worstcase model, has its minimal value in regions A and B where the highest stress occurs.This highlights a drawback of the worst-case approach: the model cannot distinguish between tension and compression stress.This is because it is based on a minimization over all unit stresses and only depends on the local volume fraction and local stress magnitude but not on stress type or direction (cf.( 27)).As for our exemplary lattice, the computed load factors for compression stress are consistently smaller than for tension stress (see also (Christensen et al. 2022, Fig. 3) for a confirmation of this observation), and buckling load factors will always be underestimated by the model in regions of tension stress.This results in higher lattice volume fraction than actually necessary to prevent local buckling at these locations.Due to this expendable lattice over-sizing, microscopic buckling is not observed at the corresponding locations for the dehomogenized design, although predicted by the worst-case model.As in the previous example, we conduct a posteriori homogenization with the "real" macroscopic stress taken from the simulated design in Fig. 26.Fig. 28 shows that in the region of compression stress (B), the error of the prediction is quite low, but in the region of tension stress (A), the worst-case model clearly underestimates the actual load factor.For an exemplary finite element in the marked region C, the predicted buckling load factor is 0.100, while the homogenized value is 0.154.Thus, this element actually has 50% higher buckling strength than predicted by the worst-case model.This can be explained by Fig. 29.The load factor for shear stress at 45 • , which is the acting stress in the inspected element, is around 50% larger than the load factor for biaxial stress, on which our worst-case model is based.This explanation also leads to a simple a priori error estimate for the worst-case model: if we know the type of stress, the error is bounded by the difference between the maximal load factor on the curve for the specific stress type and the load factor for biaxial stress.For this last example, we showed a homogeneous design with non-homogeneous microscopic buckling load factor.In a simultaneous optimization of both macro-and microscopic buckling (cf.problem C in Section 5), only the microscopic constraint is active as the structure exhibits a comparatively high macroscopic buckling load factor.Optimization with respect to only microscopic buckling (problem B) leads to non-homogeneous lattice with homogeneous microscopic buckling load factor.For this optimized design, similar deviations in the predicted worst case versus the a posteriori buckling load factor have been found in regions that exhibit shear stress.(cf. Fig. 11).In region C of our example (Fig. 26), we have shear stress acting at an angle of 45 • .The buckling load factor for this stress situation (marked by a circle) is around 50% larger than the one for biaxial compression stress.
We recall that our worst-case model acts as a safeguard against pure microscopic buckling.The possibly excessive underestimation of the microscopic buckling load factor can lead to oversized lattice struts.This can be overcome by replacing the worst-case approach with other interpolation models, e.g., a C 1 -Interpolation in the threedimensional parameter space (ρ, σ)| σ∈S 2 .However, this comes with additional computational effort for the worst-case model and, as mentioned in Remark 1 in Section 4.1, requires sophisticated interpolation schemes if the proposed method is applied for three-dimensional structures, as the parameter space will become higher dimensional.

Conclusion
We presented a method to perform two-scale optimization of lattice structures while accounting for buckling on both scales using asymptotic homogenization.Based on a parameterization of our chosen lattice structure, we obtained homogenized elastic and buckling properties.We constructed a worst-case model for the homogenized buckling load factor (Section 4.1).Both elastic and buckling characteristics were upscaled by individual, continuously differentiable interpolation models.We provided numerical examples for optimization of only macroscopic or microscopic buckling and simultaneous optimization of both.
We demonstrated that our method to obtain homogenized properties is equivalent to a special discretization of the Brillouin zone in Floquet-Bloch theory.We showed that rounding corners can have a considerable influence on the buckling strength of lattices, as it prevents stress concentrations.
We compared the performance of optimized designs predicted by the worst-case model with their dehomogenized counterparts.Compliance and macroscopic buckling were predicted very well; the predicted microscopic buckling load factor was conservative.We explained the reason for a possible underestimation and how this can be avoided by replacing the worst-case model by interpolation in the whole parameter space.A posteriori homogenization matched dehomogenized results quite well.This shows that homogenization is a viable tool to upscale the buckling behavior of lattice structures.
Although we limit ourselves to a twodimensional setting in this article, an extension to three dimensions is straight forward.Further research could include the combination of lattice, solid, and void design in optimization problems while accounting for manufacturing constraints.

Fig. 1
Fig. 1 Principal design of an equilateral triangular lattice structure with highlighted unit cell.The lattice geometry is given by primitive lattice vectors R 1 and R 2 .For our enhanced version, see Fig. 8.

Fig. 2 Fig. 3
Fig.2To prevent stress concentrations, lattice corners are rounded by circular arcs, which are defined by a given radius.The figure shows the resulting lattice boundary for two configurations with different radii r 1 and r 2 .The arc midpoint M and lattice strut width w are determined by the radius via geometric calculations, forcing the relative volume fraction of the lattice constant.

Fig. 4 Fig. 5
Fig.4Normalized homogenized Young's modulus for different relative volumes as a function of the radius of rounded corners.The radius is given relative to a unit cells edge length d.

Fig. 7
Fig.7Homogenized buckling load factor for different finite element discretizations of an RVE with a relative volume of 30%.

Fig. 9
Fig. 9 Wave vectors in the Brillouin zone (see Thomsen et al. (2018)), which correspond to homogenization with 1 to 7 cell repetitions in an RVE.Only 7 cell problems have to be solved to capture all modes corresponding to the marked wave vectors.

Fig. 10
Fig. 10 Representation of unit stresses as a unit sphere surface (σ)xy out of plane axis).Special stress types form circles (uniaxial/shear stress) or single points (biaxial stress) on the surface.Unit stresses are parameterized by spherical coordinates with zenith reference z and azimuthal reference x.

Fig. 11
Fig. 11 Homogenized buckling load factors for an RVE with 3x3 unit cell repetitions and a relative volume of 30%.Note the logarithmic scaling of the radial axis.The displayed angle is the azimuthal angle from the spherical reference coordinate system (see Fig. 10) and corresponds to the rotation of the stress around the RVE.The symmetry of the unit cell is reflected in the load factors.

Fig. 12
Fig.12Homogenized worst-case buckling load factors and the interpolating curve.For high relative volumes, numerical issues arise (marked by an o).Thus, only values up to 0.6 are interpolated, and quadratic extrapolation is applied above.
A) -C) without the compliance constraint (39) and evaluate the compliance values c A , c B and c C of the resulting designs to get the extreme points on the Pareto set.Afterwards, we discretize the intervals [c * , c A ], [c * , c B ] and [c * , c C ] equidistantly, e.g., for A) c j 0 = c * + j (c A − c * )/50, j = 0, . . ., 50.(41)

Fig. 13
Fig. 13 Left: macroscopic design setting with loadings and bearings.Right: the solution of a standard compliance problem (40) in the design domain (dashed).

Fig. 14
Fig. 14 Left: buckling mode of a solid column.Center: column with weak (gray) material and a solid (black) reinforcement in endoskeleton configuration.Right: exoskeleton configuration.The force is applied via a stiff plate at the top (red).

Fig. 15 Fig. 16
Fig. 15 Optimized function values considering only macroscopic buckling subject to a compliance and a global volume constraint.As reference, we include a topology optimization result (TO) for pure compliance minimization with ρ min = 0.001 with compliance 1.10 Nm and load factor 0.014.Selected designs are shown in Fig. 16.

Fig. 17 Fig. 18
Fig. 17 Optimized function values considering only microscopic buckling subject to a compliance and a global volume constraint.The optimized designs show increasing minimal local volume (see also Fig. 18), which yields better local buckling resistance.

Fig. 20
Fig. 19 Optimized function values of concurrent load factor maximization of pure macro-and pure microscopic buckling with compliance and volume constraint (36)-(39).Up to compliance 1.33, only the microscopic buckling load factor is active.For higher compliance values, the macroscopic buckling load factor dominates.Selected designs are shown in Fig. 20.

Fig. 21
Fig.21Dehomogenization of design C 4 with increasing number of lattice cells.For the interior microscopic load factor, we searched for the first buckling mode that has no deflection at the boundary.The macroscopic buckling load factor matches our prediction well; the microscopic and the interior load factors are better than predicted by the worst-case model (cf.Fig.19).

Fig. 22
Fig.22Dehomogenization of design C 4 with 6 (left) and 12 (right) cells in horizontal direction with visualized prebuckling displacement.The horizontal struts at the top face where the pressure load is applied exhibit sag, which is reduced for a higher number of cells.

Fig. 25
Fig. 25 Left: microscopic buckling load factor predicted by the worst-case model.Right: visualization of a posteriori homogenized microscopic buckling load factor.Load factors are thresholded at 0.1.The marked finite element in the upper region corresponds roughly to the center of the mode from Fig.24.The a posteriori homogenized buckling load factor of this element is 13% larger than predicted and differs less than 3% from the dehomogenized design.

Fig. 26
Fig. 26 Setting for shear stress example, consisting of homogenized lattice with 50% porosity.The color shows the norm of the resulting mechanical stress.Region C exhibits shear, region A uniaxial tension, and region B uniaxial compression stress.

Fig. 27
Fig. 27 Microscopic buckling load factors obtained via our worst-case model.

Fig. 28
Fig. 28 Left: microscopic buckling load factors obtained by a posteriori homogenization with the "real" macroscopic stress.Right: relative difference of load factors predicted by the worst-case model compared with a posteriori load factors.The center region C is subject to pure shear stress and has a predicted load factor, which is 50% smaller than the one obtained by a posteriori homogenization.
Fig. 29 Homogenized buckling load factors for an RVE with 3x3 unit cell repetitions and a relative volume of 30%