Topology Optimization of Multiscale Structures Considering Local and Global Buckling Response

Much work has been done in topology optimization of multiscale structures for maximum stiffness or minimum compliance design. Such approaches date back to the original homogenization-based work by Bends{\o}e and Kikuchi from 1988, which lately has been revived due to advances in manufacturing methods like additive manufacturing. Orthotropic microstructures locally oriented in principal stress directions provide for highly efficient stiffness optimal designs, whereas for the pure stiffness objective, porous isotropic microstructures are sub-optimal and hence not useful. It has, however, been postulated and exemplified that isotropic microstructures (infill) may enhance structural buckling stability but this has yet to be directly proven and optimized. In this work, we optimize buckling stability of multiscale structures with isotropic porous infill. To do this, we establish local density dependent Willam-Warnke yield surfaces based on local buckling estimates from Bloch-Floquet-based cell analysis to predict local instability of the homogenized materials. These local buckling-based stress constraints are combined with a global buckling criterion to obtain topology optimized designs that take both local and global buckling stability into account. De-homogenized structures with small and large cell sizes confirm validity of the approach and demonstrate huge structural gains as well as time savings compared to standard singlescale approaches.


Introduction
In recent years advances in additive manufacturing have facilitated the fabrication of multiscale or infill structures [1], which in turn has further promoted topology optimization of multiscale structures.The fact that structures with tailored microstructures can be manufactured means that multiscale structures are no longer restricted to theoretical research but can actually be utilized in real designs.A lot of research considering multiscale structural optimization following the original work of Bendsøe and Kikuchi [2] has been done in recent years as reviewed by Wu, Sigmund and Groen [3].Work in [4] showed that multiscale designs can be projected to singlescale using an implicit geometry description.This work was elaborated in [5,6] for 2D problems and [7,8] for 3D.
Previous multiscale structural optimization has mainly focused on compliance minimization without considering structural stability.Recently, singlescale studies of buckling optimization has become more popular.Work in [9] focused on the the use of linear buckling with the scope of solving large scale topology optimization.An extension of this is the work in [10], where buckling resistance and local ductile failure constraints are combined.Furthermore, the work in [11] considers singlescale buckling optimization, using the Topology Optimization of Binary Structures (TOBS) method, subject to design dependent pressure loads.This method is limited to singlescale optimization due to the binary nature of the TOBS method.Recent work in [12] focused on singlescale buckling optimization using the level set method.Work on multiscale stability optimization is only sparsely covered in literature.Recent work in [13] focuses on using filleted lattice structures considering a simplified local buckling formulation and yield stress constraints for compliance optimization.The method lacks the opportunity to generate true void and solid due to the local stress constraints and strut radii bounds used in the microstructure.As a result microstructures appear in the entire design domain even though global buckling is not considered and multiscale isotropic material is sub-optimal for compliance optimized designs [14].
Earlier work [15] has shown that isotropic infill may increase structural buckling resistance at the price of a small stiffness reduction even with a predefined constant infill density.To fully explore the potential buckling resistance, this study suggests to freely optimize the isotropic microstructure infill densities in the entire design domain via topology optimization.Previous studies have shown that microstructure buckling strength can be evaluated using Bloch-Floquet wave theory [16,17].Based on buckling stability analysis of a chosen microstructure evaluated at different densities it is possible to formulate interpolation functions for local buckling stress constraints.Combining this with a global buckling optimization provides a way of performing buckling optimization of multiscale structures.In this work, the global buckling load is estimated using the linearized buckling analysis framework presented in [18], as it increases computational efficiency, while producing sufficiently accurate buckling estimates as discussed in [9,19].Furthermore, the robust formulation by [20] and two-field formulation by [21] are used to allow control over the minimum density, thus allowing true void in the design.
We consider local and global buckling topology optimization of multiscale structures composed of a fixed isotropic triangular microstructure with spatially varying density.The stiffness of the isotropic triangular microstructure is very close to the theoretical maximum provided by the Hashin-Shtrikman (HS) upper bound [22] as found in e.g.[23].Thus, instead of performing a formal homogenization approach to compute the effective properties of the triangular microstructure, we can simply use the bulk and shear modulus provided by the HS upper bound to model its elastic properties [14,24,25].For pure single-load compliance minimization problems, the use of the isotropic (near-optimal) microstructure is suboptimal and results in 0-1, non-composite solutions, c.f. [14].For such problems, orthotropic rank-2 or rectangular hole microstructures are optimal, however, such microstructures have low shear stiffness and thus perform badly in terms of buckling stability [26].For this reason, we here opt for the stiffness sub-optimal but buckling-superior triangular microstructure.An added benefit of the stiffness isotropic triangular microstructure is that its buckling response is rather isotropic as well, thus making the establishing of a buckling "yield surface" considerably simpler.
A comparison of the stiffness interpolations using HS and standard SIMP (p=3) interpolations is visualized in Figure 1a.The difference between the two interpolation schemes is significant for lower densities and explains the advantage of using microstructures when considering buckling.Using the two schemes for pure compliance minimization confirms, contrary to optimally oriented rank-2 laminates or rectangular hole cells, that there is no advantage Figure 2: Flowchart illustrating the work flow of the method presented in this paper.The upper left box shows the main building blocks to set up the homogenized optimization problem.The lower left box shows an example of a design domain.The center column show a optimized homogenized design and its critical buckling mode.The second box from the right shows the de-homogenization method.Finally, the de-homogenized structure with its critical buckling mode is shown in the right most box. in using isotropic infill for compliance (C), see Figure 1c, c.f. [14].The small difference in compliances is simply a result of the different stiffness interpolations acting in the intermediate density regions that occur due to the density filter.However, for buckling load maximization large differences are visible in both design and performance as measured by the Buckling Load Factor (BLF).The design obtained using SIMP (seen in the upper half of Figure 1b) is black and white as intermediate densities are penalized.The multiscale design, based on the HS interpolation (seen in the lower part of Figure 1b), shows a large region of intermediate densities.The explanation lies in the higher stiffness for intermediate densities provided by the HS interpolation compared to SIMP as seen in Figure 1a.Comparing the performance of the two designs shows a 14.62% improvement in global buckling stability for the multiscale design.This indicates that isotropic microstructures are superior for buckling objectives.Of course this is only true if the structure does not reach local buckling of the infill before global buckling is reached.Therefore, local buckling stability must be taken into account when optimizing against buckling using topology optimization of multiscale structures.This paper will present the work flow for performing topology optimization of multiscale structures while preventing buckling on both local and global scale.This is done according to the work flow illustrated in Figure 2. The paper is organized as follows: In Section 2 the optimization problem for multiscale structures is presented.This includes a multifield method, interpolations and buckling stress constraint.Section 3 presents the results of two numerical examples demonstrating the advantage of the method.Finally, Section 4 concludes the work presented in this paper.

Buckling Optimization of Multiscale Structures
This study aims to optimize global structural stability while preventing local microstructure buckling via buckling stress constraints using a robust linearized buckling framework.The optimization problem using the robust formula-tion [20] and a two-field formulation that allows control over the minimum density [21] is written as The meaning of each term in (1) will be explained over the next several pages and subsections.First, g λ (ρ m ) represents the normalized aggregated inverse global stability objectives for m ∈ {e, b, d} being the eroded, blue print and dilated designs.The normalization factor J KS 0 is simply the value of J KS (γ i (ρ m )) for the initial design.g c and g V represent structural compliance and volume constraints, respectively.g s is the local buckling stress constraint, which is dependent on the global buckling load.The exact definition of g s is presented in the following subsections.v j and V Ω are the elemental and total structural volumes.The physical design field ρ m is constructed using a filtered density field x and a void indicator field s following the method by [21].This provides control over the minimum density while allowing void regions in the design.The detailed design procedure will be presented in the following subsection.
The compliance measure f T u depends on the external load f and the displacements u(ρ m ) for the design, m.The displacements are calculated by solving where K(ρ m ) is the linear symmetric global stiffness matrix for the considered design ρ m .The constant C * e is the compliance upper limit for the eroded design.The maximum allowed volume fraction is enforced through the dilated design using V * d .The value of V * d is continuously updated to ensure that the target volume V * b is enforced on the blue print design.This strategy is used to eliminate numerical artifacts and stabilize convergence of the optimization problem [20].
The critical buckling load is approximated using linearized buckling analysis (see [27]) based on the equilibrum in (2), defined as where the eigenpairs (γ i , ϕ i ) represent the eigenvalue and buckling mode vector, respectively.The actual Buckling Load Factors (BLF) λ i are recovered through the eigenvalues γ i as λ i = −1/γ i , ordered such that λ i ≤ λ i+1 .The global stress stiffness matrix G σ (ρ m , u) is built considering the design ρ m and the current stress state resulting from the displacements u(ρ m ).A Matlab implementation of (3) is presented in [18].There exist as many eigenpairs (γ i , ϕ i ), i ∈ B 0 as there are DOFs in the system.Only a subset of eigenpairs B ⊂ B 0 is included in the analysis.This subset should be large enough to capture all relevant buckling modes but not so large that the problem becomes computationally infeasible ( [28,29]).The critical buckling load is determined by the fundamental BLF λ 1 through f cr = λ 1 f .The Kreisselmeier-Steinhauser (KS) function [30] approximates the minimum structural buckling factor by aggregating the considered eigenvalues into one differentiable objective for each design m in the robust formulation.The version of the KS function used here differs from the one in [30] by normalizing the eigenvalues.This provides better scaling which improves stability of the optimization problem.The definition of the modified KS function is where P is the aggregation parameter and γ 0 is equal to the fundamental eigenvalue γ 1 from the initial design.Here (4) provides a smooth upper bound approximation of The result is a smooth lower bound for λ 1 .

Multifield Method and Material Interpolation
The multifield method by [21] is used in this work to generate the physical design field ρ m .It provides control over the minimum density of porous regions while also allowing fully void regions in the design.This is essential as both the buckling resistance and buckling stress limit goes to zero for ρ → 0. If nothing is done to deal with this issue the optimization will always require some amount of material in every elements to comply with the constraint.This essentially means that the method will not allow topology changes and is only useful for infill optimization of already known geometries.However, using the multifield method by [21] solves this issue as the stress limits are interpolated using the density field x while the stiffness interpolation is performed on the physical design field ρ m .
The construction procedure for the physical design field ρ m is illustrated in Figure 3.It shows how the indicator field, s, is filtered through F (s) then projected using a Heaviside function H m (s) at three different threshold levels.The density field is only filtered through F (x).The product of these two fields form the physical design field ρ m .
The purpose of the density filter is to eliminate numerical artifacts such as checker boarding and make the design mesh-independent [31].The definition of the density filter is [31,32] ỹ j = i∈N j,i w(n i )y i i∈N j,i w(n i ) where N j,i is the set of neighborhood elements within the filter radii R s and R x for the indicator field s and density field x respectively.In this work R x = 1.5∆x and R s = 4.5∆x is used throughout.The densities are weighted through w(n i ), where n i are the element center points, i.e.
The smoothed Heaviside function H m (s) (see [20,33]) is defined as where η m is the threshold value.In this work η b = 0.5 and the eroded and dilated values are determined by η b ± ∆η, where ∆η is specified explicitly for each of the numerical examples.The steepness of the projection is determined by β.When β increases the Heaviside function gets more and more non-linear.Therefore, a continuation method is used on β implying that it is initialized as β = 2 and slowly increased throughout the optimization process to β max = 256 to ensure a clear representation of void regions.The specific continuation scheme used in this work is stated in Section 3.
As discussed in the introduction, the stiffness of the considered triangular microstructure is very close to the theoretical optimum provided by the HS bounds.Therefore, instead of performing numerical homogenization to compute its properties, we may with just small errors use the stiffness propeties provided by the bounds.Hence, the HS upper bound is employed for elemental stiffness interpolation [22], i.e., where µ(ρ m j ) and κ(ρ m j ) are the effective shear and bulk moduli, respectively.The subscripts 0 and min refer to the base material and void.The shear and bulk modulus for the base material are defined as µ 0 = E 0 /(2(1 + ν 0 )) and κ 0 = E 0 /(2(1 − ν 0 )) and the stiffness of the void is determined similarly using E min = 10 −6 E 0 .The interpolation functions are used to build both the linear stiffness and stress stiffness matrices as Here K j 0,µ and K j 0,κ are the finite elemental shear and bulk stiffness matrices for element j with unit shear and bulk moduli.G j 0,µ (u j ) and G j 0,κ (u j ) are the finite elemental base stress stiffness matrices for unit shear and bulk modulus respectively.The definition of these matrices are based on the choice of finite element, hence they are not defined here for generality.For the 4-node quadrilateral elements used here we refer to the public Matlab code provided by [18].To avoid artificial buckling modes in the low stiffness regions, the stress stiffness matrix interpolation uses µ min = κ min = 0 (see [17,19,34]).

Buckling Stress Constraint
We employ a Willam-Warnke failure criterion to characterize the microstructure buckling failure under macroscopic stress situations in each structural element.The Willam-Warnke failure criterion, known from failure prediction of concrete, usually has a relatively small limit in tension but allows larger stresses in compression [35].In the case considered here, where buckling of microstructure should be prevented, the purpose of the Willam-Warnke failure surface is flipped to make it more suitable to characterize the buckling failure with limited compression stresses but allowing tension.The failure surface is determined by the stress limits in uniaxial tension σ t , uniaxial compression σ c and equibiaxial compression σ b .

Microstructure Buckling Stress Limits
In order to determine the three stress limits mentioned above, we characterize the microstructure buckling strength under arbitrary loading using the material buckling strength evaluations presented in [16,17].The material buckling strength under a given macroscopic stress situation is evaluated using linear buckling analysis together with Bloch-Floquet boundary conditions to capture all the possible buckling modes for an infinite periodic material.In short, Bloch-Floquet analysis numerically computes the unit cell buckling problem for all possible wave vectors, i.e. for buckling modes with all possible periodicities from cell periodic to infinite.The lowest computed buckling mode is the most critical one and thus used to create the buckling yield surface.The detailed formulations for the material buckling strength evaluation have been presented in [17] and are not repeated here for brevity.The general strength characteristics under arbitrary loading of a microstucture are evaluated by considering 21 evenly distributed macroscopic stress situations.This is done by interpolating between external stress states σ 0 = [σ xx , σ yy , σ xy ] ranging from . The considered isotropic stiffness-optimal microstructures are characterized (see Figure 4) for eight different relative densities (volume fractions), equally distributed in the interval [0.1, 0.8].Furthermore, the microstructure is evaluated at different angles to capture the worst case orientation.Due to symmetry in the microstructure 5 different orientation in the interval [0 • − 30 • ] are deemed sufficient to capture all orientations.The corresponding buckling strength surfaces are presented in Figure 4d by the square points.The data points for all orientation are only presented for the microstructure using x = 0.6.It is shown that the symmetry in the microstructure provides equal but flipped surfaces for the 0 • and 6  12) marked by the dashed red line.The solid blue line is the fit of (11) to the data points, marked by black points, obtained using Bloch-Floquet-based cell analysis.We remark here that the estimated cell buckling surface is scale independent (assuming infinite periodicity).This can a.o.be understood from considering the buckling yield stress of a pinned-pinned column, which is Thus, the estimated cell buckling values should be accurate for the cell size going to zero but be influenced by boundary conditions for cell sizes approaching the macroscale.
For each considered volume fraction, the buckling strengths under uniaxial σ c and biaxial compression σ b are directly obtained from the buckling surface.The tension strength, σ t , is estimated such that the failure surface has the best fit with the Willam-Warnke failure surface.Figure 4d shows the stress limits σ c , σ t and σ b for x = 0.4 by circles and the solid lines represent the fitted surfaces.As it is visible in the figure the fitted lines deviate more for the higher densities.This is a result of prioritizing accuracy for lower densities.However, the deviation results in a conservative approximation of the buckling stress limit and therefore still prevents local buckling for all densities.
Similar to [36], two term interpolation schemes are employed to interpolate the buckling strength as function of the relative density from the eight evaluated microstructures.Thus the function to be fitted is given as where x j is the elementwise relative density of the microstructures.The subscript k ∈ {c, t, b} indicates the strength for uniaxial compression, uniaxial tension and biaxial compression and the exponent n 0 = 3 is chosen to represent the buckling-to-relative-density exponent of a column for the microstructures at low density regimes.The coefficients b 0 and b 1 are obtained using curve fitting.Figure 4e shows the fitted curve (solid blue line) and data points (black points) for the buckling strengths under uniaxial compression.All coefficients, b 0 and b 1 , for all the buckling strengths are presented in Table 1.
Dense microstructures are more likely to fail due to plastic yielding than local buckling [37].In this work we do not consider plastic yielding.Therefore, the buckling stress constraint in the high volume fraction limit is relaxed using a Heaviside function equivalent to (7) Here ψ is the relaxation parameter, which in this work is defined as 10 times the stress limit for a solid unit cell using (11), i.e. ψ = 10 σk (1).The threshold η is chosen at the highest density of the known data points such that η = 0.8.Finally, a sharpness parameter of β = 50 is used to get a relatively smooth yet steep relaxation at the threshold value.
The relaxed buckling limit under uniaxial compression is shown by the dashed red line in Figure 4e.The effect is that the buckling stress constraint is made inactive for high densities by overestimating the buckling stress value.

Willam-Warnke Failure Surface
Using the interpolated stress limits defined by (11), the density dependent Willam-Warnke failure surface is defined using the unified approach from [35] in terms of the equivalent stress, σ eq (ρ m j , x j ), given as where J 2 and I 1 are the second invariant of the Cauchy tensor σ and invariant of the deviatoric stress tensor s, respectively.The term G(I 1 (ρ m j ), x j ) is a function determined from the first stress invariant of the deviatoric stress.All of these are defined in Appendix A. Finally, the term α(ρ m j , x j ) defines the shape of the yield surface.It is defined as with A( x j ), B( x j ), C( x j ), D( x j ) and E( x j ) determined from the stress limits and θ(ρ m j ) being the modified Lode angle.The definitions of all these are presented in Appendix A.

Buckling Stress Constraint Formulation
The p-norm stress measure proposed by [38] is used to approximate the maximum equivalent stress as max The p-norm approximation of the max-function depends on the value of p.As p → ∞ the p-norm approaches max ∀ j (σ eq (ρ m j , x j )) from above.This means that the p-norm will always overestimate the actual maximum constraint value and therefore act as a conservative constraint.
The actual buckling stress constraint is enforced on the p-norm stress measure under the critical global buckling load.This results in the buckling stress constraint where c is the correction presented by [39] for better approximation of max ∀ j (σ eq (ρ m j , x j )).The correction is defined as where n is the design iteration number and α ∈ (0, 1] controls the update speed between design iterations.In this work α = 0.1, ∀n, is used and the correction value is initialized as c 0 = 1.The final constraint is formulated by normalizing (16) with the initial value at n = 0, which yields where σ0 PN = −σ 0 PN /J KS ,0 .The optimization problem is solved using the method of moving asymptotes (MMA) [40] based on the sensitivity of the objective and constraints presented in Appendix B.

Numerical Results
The approach presented in this paper is verified and demonstrated on two examples.The results illustrate the superiority of multiscale structures using isotropic microstructures over single scale structures.Both examples use E 0 = 1 and ν 0 = ν min = 0.3 and the coefficients from Table 1 to define the buckling failure surface.All examples use a β-continuation where β is initialized as β = 2.After the first 125 iterations β is doubled every 75 iterations until β max = 256 is reached.

Uni-axial Compression
The first example compares singlescale with multiscale optimization.Therefore three optimizations are performed.The first is a singlescale optimization using SIMP as stiffness interpolation, performed using the code from [18].The second is a multiscale optimization using HS for the stiffness interpolation.Finally, a multiscale optimization using the local buckling stress constraint with HS stiffness interpolation is performed.This will be denoted HS-S where the added S indicates the use of the buckling stress constraint.All three optimizations are performed using the design domain in Figure 5 which is discretized using 200 × 200 bilinear quadrilateral elements.The robust formulation is set up using ∆η = 0.25 to generate the eroded and dilated designs.The center nodes on each side of the domain are restricted against vertical displacements.The node in the center of the domain is restricted against horizontal displacements.Two inward pointing distributed loads are applied around the horizontal symmetry axis.The widths of these loads are L load = 0.04 with a magnitude of 10 −3 at each load.Passive solid regions are located at the loads with a width of 0.025.The design is forced to be symmetric around the horizontal and vertical center axes, due to the symmetry in the design problem, but no symmetry is used on the physics to allow for both symmetric and unsymmetric buckling modes.The objective function is a weighted average between compliance and buckling using weight {ω | 0 ≤ ω ≤ 1}.This yields the optimization problem min where g w (ρ m ) is the weighted objective function.The objective function consists of the compliance C and aggregated eigenvalues J KS normalized using the values for the initial design, i.e.C 0 and J KS 0 .Furthermore, a global volume constraint g V (ρ d ) is used with a target volume on the blue print design of V * b = 0.15.Finally, the local buckling stress constraint is used only on the last of the three multiscale optimizations, hence the gray color.The optimization is performed by first optimizing the single-and multiscale problems for pure compliance i.e. ω = 1.Then ω is lowered step by step using the optimized design from the previous ω until ω = 0 is reached.Finally, the multiscale problem with the active buckling stress constraint is optimized using the singlescale designs as initial guesses.
For the case using ω = 0.5 the three designs in Figure 6 are obtained.The figure also presents the first two mode shapes with the strain energy density calculated on element level by φ j = ϕ T i j k j ϕ i j as indicated by the color scheme.It is seen how the three optimizations produce different designs performing at different levels for the BLF and almost similar level for compliance.All three designs have similar solid oval shells as the dominating part of the structures.These shells carry the majority of the load supported by an internal shear web that increases the buckling stability.The 10  difference in the designs lie in the infill that connects the shells.The singlescale material produces a design which uses thick bar connections between the shell walls.These help stabilize the structure by restricting the shell from being pushed outward during compression.The disadvantage of these bar connections is that the shell is only stabilized discretely.Therefore the shells have to be thicker to prevent buckling between these connections.This results in a BLF of λ 1 = 6.17.
Looking at the multiscale design in Figure 6b it is seen that a more homogeneous infill with intermediate density is used.This means that the shells are connected on their entire inner surfaces.Therefore the shells are restricted from being pushed outwards and the shell itself only has to carry the horizontal load.The first two modes shapes are very similar to the singlescale modes.The critical BLF has increased to λ 1 = 6.68 and the second mode is active at λ 2 = 6.82.Based on this, the multiscale design is superior for buckling stability.However, using infill with an intermediate density introduces a risk of having local buckling before global buckling.This is not prevented in the simple multiscale design in Figure 6b.The risk of having local buckling is tested using de-homogenization based on the approach described in [7] which is briefly reviewed in Section 3.2.1.The first buckling mode is presented in Figure 7a.Local buckling of the microstructure is clearly visible in the right side of the design.As a result the BLF is reduced to only λ 1 = 0.96, which is significantly lower than for the homogenized design.This BLF, together with the buckling mode, confirms that local buckling occurs if local stability is not taken into account during the optimization.
Finally, the multiscale design with the buckling stress constraint is presented in Figure 6c.The design is somewhat a mix of the multiscale design without the local buckling stress constraint and the singlescale design.The same bar structures as in the singlescale design are visible inside the shell.However, the bars are no longer solid, but instead they are wider and filled with material of intermediate density.This means that the distance between weak free hanging shells is reduced and therefore stability is increased.Looking at the first two buckling modes in Figure 6f and 6i it is seen that they are similar to the modes from the singlescale and unconstrained multiscale designs.The eigenvalues of the two modes are λ 1 = 6.54 and λ 2 = 6.63 respectively, which means that both are likely to be active at the same time.The critical BLF is λ 1 = 6.54 which is a slight improvement of 6% in buckling stability compared to the singlescale design.To verify that the buckling stress constraint works as anticipated, the design is de-homogenized and the two first buckling modes are presented in Figure 7b and 7c, respectively.The buckling modes are clearly global, which confirms that the buckling stress constraint prevents local buckling.The buckling modes are switched compared to the homogenized results.The reason is most likely found in the de-homogenization of the two thick bars with intermediate density.The buckling modes for the homogenized design shows that the strain energy density in the two bars is low for the first buckling mode, but large for the second buckling mode.This means that the second buckling mode is more sensitive to the size of the de-homogenized microstructure used in the two bars.Therefore, the BLF of the second buckling mode is likely to drop below the first buckling mode when de-homogenized.This is indeed the case with a BLF of λ 1 = 6.47 which is slightly lower than that of the homogenized design.The BLF of the second buckling mode is λ 2 = 6.65 which is actually a bit higher than the corresponding buckling mode for the homogenized design.Nevertheless, both buckling modes have BLF much closer to those of the homogenized design than it was the case when the buckling stress constraint was not included.The reason for the lower BLF value is examined thoroughly in the following section.Furthermore, a solution is presented to produce de-homogenized designs that perform at the same level or higher than the homogenized designs.
The results for all the optimization setups with all the values of ω are presented in Figure 8.A curve is fitted to the data points of each of the optimization setups to visualize the position of the Pareto front.Selected designs from all three test cases illustrating similarities and differences are also presented in the figure .The figure shows that for pure compliance optimization both single and multiscale optimizations reach similar designs with almost identical performance.Looking at the pure buckling case it is seen that the multiscale optimization without the local buckling constraint is the superior choice.However, the low density region in the center of that design is likely to buckle locally at significantly lower loads than the global buckling load at λ 1 = 10.81.The design which takes local buckling stability into account shows a lower global buckling resistance at λ 1 = 9.24.But it is still better than the single scale design which experiences global buckling at λ 1 = 6.86.To sum up, the Pareto fronts show a potential buckling stability increase of 57.6% when using multiscale material with isotropic microstructures instead of single scale.When local buckling stability is taken into account the increased performance is reduced to 34.7% which is still a significant improvement compared to the singlescale design.It does however, leave a potential room for improvements of 22.9 percentage points.This could possibly be exploited using multiscale microstructures with enhanced buckling stability [17,26].

Modified L-Beam
The second problem is a modified version of the frequently used benchmark L-beam problem.The domain is illustrated in Figure 9 with the domain rotated 90 • counter clockwise and the load in the vertical downward pointing direction.This change provides a more interesting buckling response as seen in the following.The width of the distributed load is 0.1 and its center is located 0.25 from the right side of the domain.The magnitude of the load is 10 −3 .The upper left quadrant of the domain is passive void.Furthermore, two passive solid regions are present in the domain.The first is located at the load with the same width as the load and a height of 0.05.The second is positioned at the lower left corner with the dimensions 0.05 × 0.05.The left boundary of the active domain is simply supported in the horizontal direction.The part at the passive solid domain is also supported in the vertical direction.The domain is discretized by 100 × 100 bilinear quadrilateral elements.The robust formulation is set up using ∆η = 0.1 to generate the eroded and dilated designs.Furthermore, all designs are obtained from uniform initial guesses.The problem is solved using both singlescale material (SIMP) and multiscale material with the buckling stress constraint (HS-S).The optimized designs are post-processed using de-homogenization where the multiscale material is projected onto a fine mesh and interpreted as a triangular microstructure with a varying bar width w using the method by [5,7].For this reason x min = 0.27, corresponding to a slenderness ratio of L/w = 10.3, is used to ensure that individual bars in the microstructure can always be resolved using 3 elements across the width of the bar when projected onto the fine mesh.This is done to prevent node connections and reduce artificial stiffness through shear locking.In reality the choice of x min should be based on the available manufacturing method and associated limitations.
All the optimized homogenized designs are presented in Figure 10 with the corresponding buckling modes presented in Figure 11.The homogenized design for the compliance minimization is presented in Figure 10a.The design has a characteristic vertical bar which obviously has low buckling resistance.This is visible in the first two buckling modes shown in Figure 11a and 11d where the critical buckling mode is clearly visible in the vertical bar.The critical BLF also confirms that the design is not buckling resistant with a value of λ 1 = 0.13.The compliance of the eroded, blue print and dilated designs are C m = [6.9× 10 −5 , 5.33 × 10 −5 , 4.68 × 10 −5 ], respectively.
Next a buckling optimization is performed using the singlescale material and a compliance constraint allowing 1% higher compliance on the eroded design.The optimized design is presented in Figure 10b.The design is slightly different from the compliance optimized design.The single vertical bar has thickened to increase buckling stability.The critical buckling mode in Figure 11b is distributed more evenly over the structure and as a result the critical BLF increases with 715% to λ 1 = 1.06.The higher order modes are active at significantly higher BLFs.Furthermore, they are mainly located at the thinner inclined bars in the lower left part of the design as seen in Figure 11e and 11h.The  compliance values are C m = [6.97× 10 −5 , 6.36 × 10 −5 , 5.77 × 10 −5 ], which for the eroded design is an increase of 1% but for the blue print and dilated designs correspond to larger increases.Finally, the optimization is performed using the multiscale material with the compliance constraint allowing 1% higher compliance on the eroded design.The optimized design is presented in Figure 10c.The design is visibly different from the singlescale designs.The single vertical bar has been split into two to increase buckling stability.The critical buckling mode in Figure 11c is now located in the two bars in the upper part of the domain.As a result the BLF increases with a factor of approximately 20 compared to the compliance optimized design, resulting in a critical BLF of λ 1 = 2.73.The second buckling mode in Figure 11f is distributed over the entire structure.The third buckling mode in Figure 11i shows buckling at the transverse bar originating from the lower left corner of the domain.This is very similar to the corresponding mode for the compliance optimized design.Common for the two higher order modes presented are that their BLFs are closer to the critical BLF.The difference between the first and third BLF is only 0.22.The compliance values are C m = [6.97× 10 −5 , 6.66 × 10 −5 , 6.40 × 10 −5 ] which again is a 1% increase for the eroded design.The blue print and dilated designs have larger increases, just as it was the case for the singlescale buckling optimized design.

De-Homogenized Buckling Analysis
Using multiscale material for any kind of optimization is only beneficial if the de-homogenized designs perform at the same level as the homogenization based designs.To verify that this is the case the multiscale design in Figure 10c is de-homogenized.The de-homogenization procedure is based on the approach described in [7].First the densities are mapped to a fine mesh.In the current work a refinement ratio of 30 from the coarse to fine mesh is used resulting in a fine mesh using a total of 9 × 10 6 elements.Then the densities are converted to widths w of the bars in the microstructure This is used to create an implicit geometry description ρi , i ∈ [1,3] for each of the three layers where H is the Heaviside function, S ∈ [−1, 1] is a triangular wave function and x 1 and x 2 are the spacial coordinates.The cell size is determined through the periodicity    where is the absolute length describing the periodicity between two parallel bars in the microstructure as illustrated in Figure 12.The orientation of the individual layers is determined through θ i .Given the triangular microstructure the layers are positioned in intervals of 60 • .Finally, the layers are combined to create an implicit geometry of the de-homogenized structure Examples of de-homogenization using four different cell sizes are presented in Figure 13.The example using = 0.01 has very fine microstructures which visually looks similar to the homogenized design in Figure 10c.For the higher values of the size of the microstructures increases to a level where only a few triangles are present.This obviously introduces a risk of loosing connectivity throughout the structure as cell sizes exceed the length scale used on the homogeneous design.Nevertheless, it is interesting to look at large cell sizes to gain insight into the validity of the homogenization assumptions.An example is presented in Figure 13c where a cell size of = 0.15 is used.The region in the red circle is disconnected as a result of the coarse microstructure.This will severely reduce stiffness and stability of the structure.Additionally, coarse microstructures will result in free hanging members which are not contributing to neither stiffness nor stability.This is especially pronounced in Figure 13c 12.This results in a total of 320 test cases that are used for the buckling post validation.The first buckling mode of four selected cell sizes are presented in Figure 14.For each cell size the minimum and maximum of the first buckling mode as well as the buckling mode with λ 1 closest to the average are presented.
Common for the structures using = 0.01 is that the BLFs are all within a very tight range.Furthermore, the buckling modes are very similar for all the cell positions.Comparing the buckling modes to those of the homogenized design it is seen that they do not correspond exactly to the first buckling mode.The upper part of the structure experiences similar buckling compared to the homogenized design.However, the de-homogenized design also experiences buckling of the two bars in the lower left part of the structure.The average BLF of the de-homogenized structures is λ 1 = 2.78.This is slightly higher than the BLF of the homogenized design which was λ 1 = 2.73.The reason for not seeing the same first buckling mode in the de-homogenized designs as the homogenized design is most likely found in the de-homogenization of the thinner bars.A small dilation or erosion to the width of these bars will change their stiffness significantly.Even though the buckling modes of the de-homogenized designs are not identical to those of the homogenized design the critical BLF is still higher.This means that for = 0.01 the de-homogenized design performs better than the homogenized design.In general the homogenized designs should act as an upper bound.However, added volume and too stiff finite elements in the coarsely discretized microstructure as well as possible local thickening of critical features can increase the buckling stability of the de-homogenized structures.
Looking at the coarser cell sizes reveals additional different buckling modes and an increasing range of BLF for each cell size.In general the critical BLF decreases in value resulting in a worse performance than the homogenized design.This is largely due to the size of the effective region of the microstructure being reduced.The center column of Figure 14 highlights how smaller cell sizes are able to effectively use a larger part of the structural domain.The problem is that larger cell sizes risk producing larger free hanging members at the boundary of the microstructure.These members will not carry any load and therefore not contribute to stability.This is seen by the blue color in these members indicating low strain energy density.The worst case scenario is that the effective region of the microstructure is reduced in all direction by the width of the cell .Therefore for → 0 the effective load carrying region of the de-homogenized design will be identical to the homogenized region.This effect highlights the assumption that the homogenized multiscale method only holds in the limit of infinitely periodic microstructure.However, based on the presented example it seems like a cell size of = 0.01, corresponding to 1/100 of the largest dimension of the design domain is small enough to provide consistent results with less than 4% deviation from the homogenized result.
In some cases the structure is subject to local buckling of the microstructure before global buckling occurs.An example of this is presented in Figure 14i.Here buckling occurs in a single bar in the microstructure at the BLF λ 1 = 2.22.The reason for this local buckling can be found in the open microstructure at the boundary next to the bar experiencing local buckling.The fact that the microstructure is incomplete weakens the surrounding microstructure.This leads to local buckling even though the homogenized local buckling stress constraint is satisfied.
As discussed earlier there is a risk of having disconnected structures when is increased to much.An example of this using = 0.15 was shown in Figure 13c and the critical buckling mode for that design is presented in Figure 14j.The structure has a critical BLF of λ 1 = 0.53 which is even lower than the buckling optimized single scale design.The color scheme clearly visualizes how the transverse bar is not carrying any load and therefore not contributing to the stability.
All the average values of the volume fractions, stiffness per volume and BLFs relative to the cell size are presented in Figure 15.The plots also indicate the minimum and maximum values obtained using the different cell positions at each cell size.The volume fractions of the de-homogenized designs are compared to the target for the homogenized design.The stiffness per volume and buckling values of the de-homogenized designs are compared to the singleand multiscale homogenized blue print designs.The comparison of the volume fraction in Figure 15a shows that the de-homogenized volumefraction V * b in general is higher than the target volume fraction on the homogenized design.As expected, the volume fraction varies increasingly when cell sizes are increased.
To filter the effect of volume change out of the discussion, we also compute the stiffness per volume ratio as presented in Figure 15b.Not unexpected, the general tendency is that structures get softer when cell sizes increase.For the finest microstructure using = 0.01 the variation is very small and the average value mean(1/(CV)) = 7.2452×10 4 is very close to that of the homogenized design 1/(C b V) = 7.5120 × 10 4 with a difference of only −3.55%.
The BLFs are presented in Figure 15c.Colored dots represent the individual critical BLFs visually identified as global buckling (green) or local buckling only present in the microstructure (red).The figure shows that for = 0.01 only global buckling modes are present.However for = [0.08,0.09] local buckling modes occur at some cell positions.For ≥ 0.1 only global buckling modes are present.This is partly because the length scale of the microstructure enters global length scale and partly because the microstructure is so poorly connected that only solid regions are left to buckle.Figure 15c also indicates that for increasing cell sizes the BLF varies more between individual cell positions.This demonstrates that the homogenized multiscale method holds, not only in the limit of infinite periodicity, but at least up to the length scale imposed via the filter radius R x ≈ = 0.015 on the homogenized 18  infill density as indicated by the red line in Figure 15.
Based on Figure 14 the cell size should be equal to the element size used for the homogenized optimization.Alternatively, it should be equivalent to the length scale imposed on the density field x used for the homogenized optimization, but this is yet to be investigated.

Closing the Boundary Shell
The previous section showed instabilities mainly related to boundary effects at void interfaces due to incomplete microstructures.One way to overcome this problem and make the de-homogenized structures more stable and independent of the cell size and position is to close the microstructure at void interfaces.This is now done using an erosion-based interface detection [41] with a local erosion threshold parameter ∆η j where w j is the width of the shell and R shell is the filter radius in the standard filter.In this work the PDE filter based on a Helmholtz-type partial differential equation [42] is used to filter the indicator field The relation between R shell and r shell is The shell width w j is determined locally according to the densities using (20).The shell emerges from the difference between the eroded and non-eroded indicator fields where η = 0.5 and H is a sharp heaviside step function.An illustration of the erosion-based interface identification is presented in Figure 16.
The method is applied to the de-homogenized structures presented in the previous section.Illustrations of the structures shown in Figure 13 with added shells are presented in Figure 17.Here the original structure is black and the added shells are shown in blue.The structure using the smallest cell size i.e. = 0.01, only has a very thin shell added.As a result the volume of the structure is only increased very little i.e. 0.73% at average.However, for increasing cell sizes the thickness of the shell increases and the added shell has a larger impact on the volume of the structure.
The case using = 0.15 presented in Figure 13c had a region where the structure was disconnected prior to the addition of the shell.This resulted in very poor stiffness and buckling performance as shown in Figure 14j and 15.The addition of the shell helps reconnect the otherwise disconnected region as shown in Figure 17c.This correction should result in increased stiffness and buckling stability.
The free hanging members, which were not contributing to neither stiffness nor stability before the shell was added, are now all connected to the shells.This helps activate them such that they can contribute to the performance of the structure.Especially for larger cell sizes, as illustrated in Figure 17c and 17d, this will reactivate a large region of the microstructure and thereby have a significant impact on the stiffness and buckling stability.
The same 320 test cases that were used in Section 3.2.1 are evaluated with the de-homogenized structures where shells are added.Figure 18 shows the first buckling mode of six selected cases using five different cell sizes.Common for all six example is that they experience global buckling.The example in Figure 18a uses = 0.01 and has a critical BLF of λ 1 = 2.96 which is an 6.47% increase compared to the same structure without the shell.The buckling mode of the structure with the shell also appears to be more global than the one for the structure without the shell.These results are mainly due to two factors.First, the shell adds material to the structure, which increases stiffness and stability.Second, the finite element discretization of the microstructure can increase stiffness and stability artificially.
The global buckling mode presented in 18c, using = 0.08, is a clear testimony to the improved stability obtained using the shells.The corresponding structure without the shell experienced local buckling of the microstructure due to incomplete microstructures resulting in weakening near void interfaces as seen in Figure 14i.The added shell closes off the microstruture which can then exploit the full potential of the multiscale material all the way to the void interface.For this very reason the de-homogenized structures experience better stability and are less sensitive to the position of the microstructure.
For very coarse microstructures it was seen that a risk of having disconnected regions in the structure was introduced, see Figure 13c and 14j.Adding the shell along void interfaces helped reconnect these regions as shown in Figure 17c.The first buckling mode of this structure is presented in Figure 18d.The mode is global and active at a BLF of λ 1 = 3.09.This is a significant improvement of 483% highlighting the stability gained by using the shell.
Even for the coarsest case using = 0.2, where the microstucture is almost entirely incomplete, stability is improved significantly using the shell.Figure 18e and 18f show the buckling modes with the minimum and maximum BLFs respectively.With a minimum BLF of λ 1 = 2.70 and a maximum BLF of λ 1 = 3.25 the difference is only 20.4%.This is not negligible but still a significant improvement compared to the structure without the shell.
The volume fraction, stiffness per volume and BLFs for all 320 cases are summarized in Figure 19.The volume constraint is increasingly violated with increasing cell sizes as seen in Figure 19a.For the finest case the average volume fraction is V = 0.2062 which is an increase of only 0.74% compared to the structure without the shell.However, for the coarsest cell size the average volume fraction is V = 0.2252 resulting in a significantly higher increase of 11.6%.Furthermore, the variation of the volume fraction increases making it difficult to predict the volume fraction after de-homogenization.With the added volume comes a benefit of added stiffness.However, this is only the case for cell sizes up to = 0.07 as shown in Figure 19b.For larger cell sizes the stiffness per volume measure decreases below that of the smallest cell size.Nevertheless, the average stiffness per volume never exceeds a variation above 2.82% even for the largest cell size.The variation in stiffness per volume relative to cell position increases in correlation with cell size.However, this is not as significant as it was the case when the shell was not added.
Finally, the critical BLF is presented in Figure 19c.Compared to the results without the shell the stability is significantly improved.The average BLF is almost constant for all cell sizes.Furthermore, the variation is significantly reduced compared to Figure 15c.Looking at the critical buckling modes shows that all 320 cases experience global buckling before local.The case using = 0.01 has an average BLF of λ 1 = 2.96 which is a 6.47% improvement compared to the structure without the shell.This is a significant improvement considering the small price of only adding 0.74% material.Furthermore, the compliance is reduced by −3.16% meaning that there is a stiffness gain associated with the shell.

Conclusion
This paper has presented an approach to obtain topology optimized multiscale designs that takes both local and global buckling stability into account.This was done using homogenization based topology optimization with isotropic multiscale material modeled by the Hashin-Shtrikman upper bound.Local buckling stability was ensured using a local buckling stress constraint.The constraint was based on the buckling stress limit obtained using Bloch-Floquet-based cell analysis on eight different relative cell densities.These limits where used to determine the shape of the Willam-Warnke yield surface to formulate an effective buckling stress surface.The combination of the local buckling stress surface and the global multiscale formulation provided a method capable of obtaining topology optimized designs that takes both local and global buckling stability into account.
The method was demonstrated on two examples.First a uni-axial compression setup was evaluated using both singlescale and multiscale material.This example confirmed, what has been postulated for years, that isotropic microstructures enhance structural stability and are superior to singlescale structures when buckling is considered.The example also showed that even though stability was improved using the isotropic microstructure, there is still room for further improvements.Here the use of hierarchical microstructures with enhanced buckling stability [17] is an obvious choice to examine.
The second example was the modified L-beam.A significant factor 2.5 improvement of the BLF was observed when using the multiscale instead of the singlescale design.The de-homogenization of the multiscale design without adding a shell showed very good performance in compliance and BLF with only a small increase in volume fraction for the finest microstructure.However, for the larger cell sizes growing variations in performance related to cell positions were observed for both volume fraction, compliance and BLF.These instabilities were reduced significantly by adding a shell at the void interfaces to the de-homogenized designs.No local buckling modes were observed in any of the 320 test cases and the variation related to cell position was reduced drastically for the BLF.In fact, all the de-homogenized structures perform better than the homogenized design in terms of BLF.Regarding the stiffness per volume measure all cell positions for cell sizes up to = 0.14 are within 5% of the homogenized design.Finally, the volume fraction of the de-homogenized designs increase with cell size.However, all cell positions for cell sizes up to = 0.04, i.e. more than twice the length scale of the infill in the homogenized design, are within 5% of the limit imposed on the homogenized design.Even for extremely large cells sizes, errors are moderate and confirm immense computational savings by formulating the complex buckling optimization problem as a multiscale problem that takes local buckling effects into account.
A( x j ), B( x j ), C( x j ), D( x j ) and E( x j ) which are defined as where r c ( x j ) and r t ( x j ) are defined as The modified Lode angle θ(ρ m j ) in ( 14) is defined as with θ(ρ m j ) being the Lode angle [43] defined as where the third invariant of the diviatoric stress is introduced.The defineition of the term G(I 1 (ρ m j ), x j ) in ( 13) depends on the yields surface.Here the Willam-Warnke yield surface is used but definitions for other yield surfaces can be found in [35].The definition of G(I 1 (ρ m j ), x j ) is where the parameter β( x j ) is defined as The three invariants are defined as with the Cauchy stress calculated using the HS upper bound for interpolation defined in (8).The diviatoric stress tensor is defined as

Appendix B. Sensitivity Analysis
The optimization problem (1) is solved by updating the design iteratively using the Method of Moving Asymptotes (MMA) [40] based on the gradients of the objective and constraints.The sensitivity of the objective is determined assuming a distinct eigenvalue, γ i and that the eigenvector ϕ i is K-normalized, the sensitivity of the eigenvalue γ i with respect to ρ m j is defined as where v is the adjoint vector which is determined solving The sensitivity of the KS function with respect to ρ m j is The sensitivities of the objective g λ with respect to x j and s j are determined using the chain rule The sensitivities of g c with respect to ρ m j is calculated by and equivalent to (B.4) the chain rule is applied to get the sensitivity with respect to x j and s j .The gradients of the buckling stress constraint in (18) with respect to x j and s j are also determined using adjoint sensitivity where σ sum is the sum in the p-norm such that σ PN = σ 1/p sum and v σ i is the adjoint vector.The sensitivity of J KS is defined by (B.3) and the sensitivities of σ eq , κ, µ, filters and projections are all presented in the following section.The element adjoint vector v σ i contains the degrees of freedom for each element i.The adjoint field is calculated by solving where K is the global stiffness matrix.The derivative of g s with respect to u on element level is defined in the following section.

Appendix B.1. Sensitivity of Stress Constraint Components
The sensitivities of σ eq with respect to x j and sm j used in (B.The gradient of β used in (B.9) is The gradients of the stress inveriants are defined as Finally, the gradients of I 1 , J 2 and J 3 are defined as

Figure 3 :
Figure 3: Illustration of how the physical design field ρ m is constructed from the density field x and the void indicator field s.

Figure 4 :
Figure 4: Illustration of triangular microstructures and stress limit fits: (a -c) Triangular microstructures at different relative densities with red line indicating one unit cell oriented at 0 • relative to the coordinate system in (a).(d) Fitted buckling stress surfaces compared to Block-Floquet-based cell analysis (BF).Squares are results from the cell analysis, black circles indicate the values of σ c , σ t and σ b for x = 0.4.Finally, solid lines are the fitted surfaces at different densities x.(e) Buckling stress limit in compression σc defined by (12) marked by the dashed red line.The solid blue line is the fit of(11) to the data points, marked by black points, obtained using Bloch-Floquet-based cell analysis.

30 •
as well as 10 • and 20 • angles.The worst case orientations are 0 • and 30 • , i.e. when one of the uniaxial load cases are aligned with one of the bars in the microstructure.Therefore, the Willam-Warnke surfaces have to be fitted accordingly.Only the data points for the 0 • angle are presented for the remaining densities, meaning that the fitted Willam-Warnke surfaces should match only for the upper left data points.

Figure 5 :
Figure 5: Illustration of the design domain used for the uni-axial compression optimization.

Figure 6 :
Figure 6: Optimized designs and the first two buckling modes using ω = 0.5 and color indication of the normalized strain energy density log(φ e /φ max ): (a) Singlescale design (SIMP), (b) Multiscale design (HS), (c) Design obtained using the multiscale formulation with local buckling stress constraint (HS-S), (d) First buckling mode for the SIMP design, (e) First buckling mode for the HS design, (f) First buckling mode for the HS-S design, (g) Second buckling mode for the SIMP design, (h) Second buckling mode for the HS design, (i) Second buckling mode for the HS-S design .

Figure 7 :
Figure 7: Comparison of de-homogenized results of the two multiscale designs from Figure 6.The color scale illustrates the normalized strain energy density log(φ j /φ max ).(a) First buckling mode without buckling stress constraint (HS), (b) First buckling mode with buckling stress constraint (HS-S), (c) Second buckling mode with buckling stress constraint (HS-S).

Figure 8 :
Figure 8: Pareto fronts of the three uni-axial optimization cases, i.e.SIMP is single scale, HS is multiscale using Hashin-Shtrikman and HS-S is the same but with the buckling stress constraint included.Results are similar for compliance dominated problems.However, both multiscale optimization's show superior performance in buckling.

Figure 9 :
Figure 9: Illustration of the design domain used for the L-beam optimization problem.

Figure 11 :
Figure 11: Buckling modes for designs using both single-and multiscale material.The normalized strain energy density log(φ j /φ max ) is indicated by the color scale.(a,d,g) Singlescale buckling modes of compliance optimized design, (b,e,h) Singlescale buckling modes of buckling optimized design, (c,f,i) Multiscale buckling modes.

Figure 12 :
Figure 12: Illustration of the cell size (right) and position (left).The positions is defined at the center of the cell.Here the red cell is located at (0, 0) and three remaining position are illustrated at (−0.25, 0.25), (0.25, −0.25) and (0.5, −0.5).

Figure 13 :
Figure 13: De-homogenization of the multiscale buckling optimized design at different unit cell sizes .

Figure 15 :
Figure 15: Convergence of the volume fraction, stiffness per volume and BLF for the de-homogenized structures relative to the cell size.SIMP and HS-S mark the values of the homogenized designs.The upper bound for the volume fraction on the blue print design is marked as Target: (a) Volume fraction, (b) Stiffness per volume, (c) BLF with colors indicating global (green) and local (red) buckling modes.

Figure 16 :
Figure 16: Illustration of the erosion-based interface detection method used to generate the shells with variable thickness.

2 Figure 17 :
Figure 17: De-homogenization for different unit cell sizes with shells used to close off the microstructure.Black indicates the original structure while blue color represents the added shells.

Figure 18 :Figure 19 :
Figure 18: Comparison of buckling modes obtained using different cell sizes and positions during the de-homogenization including added shell.The color scale illustrates the normalized strain energy density log(φ j /φ max ).

−
the gradient of σ being determined by the gradients of the HS interpolation in(8).The derivative of g s with respect to u j used on the right hand side of (B.8) is calculated on element level as ∂g s 2A cos( θ) sin( θ) C cos( θ) + D cos 2 ( θ) + E − A cos 2 ( θ) + B C cos( θ) + D cos 2 ( θ) + E