Data-driven topology optimization of spinodoid metamaterials with seamlessly tunable anisotropy

We present a two-scale topology optimization framework for the design of macroscopic bodies with an optimized elastic response, which is achieved by means of a spatially-variant cellular architecture on the microscale. The chosen spinodoid topology for the cellular network on the microscale (which is inspired by natural microstructures forming during spinodal decomposition) admits a seamless spatial grading as well as tunable elastic anisotropy, and it is parametrized by a small set of design parameters associated with the underlying Gaussian random field. The macroscale boundary value problem is discretized by finite elements, which in addition to the displacement field continuously interpolate the microscale design parameters. By assuming a separation of scales, the local constitutive behavior on the macroscale is identified as the homogenized elastic response of the microstructure based on the local design parameters. As a departure from classical FE 2 -type approaches, we replace the costly microscale homogenization by a data-driven surrogate model, using deep neural networks, which accurately and efficiently maps design parameters onto the effective elasticity tensor. The model is trained on homogenized stiffness data obtained from numerical homogenization by finite elements. As an added benefit, the machine learning setup admits automatic differentiation, so that sensitivities (required for the optimization problem) can be computed exactly and without the need for numerical derivatives – a strategy that holds promise far beyond the elastic stiffness. Therefore, this framework presents a new opportunity for multiscale topology optimization based on data-driven surrogate models


Introduction
Supported by developments in advanced manufacturing, mechanical metamaterials with tunable microstructure and controllable properties have made significant strides towards realizing materials by design [1][2][3][4][5][6].Despite all progress, key challenges have persisted, which can be exemplified by truss-, plate-, and shell-based cellular materials, which have dominated the design of metamaterials over the past decade.Most truss-based architectures exhibit poor scaling of stiffness and strength with relative density due to bending deformation of struts [7].Plate-based designs exhibit improved stiffness with relative density than trusses and have been shown to reach upper bounds on the achievable effective stiffness of cellular media [6,8].However, both types of architectures suffer from stress localization at the junctions of beams or plates, which results in early failure and poor recoverability [9][10][11].Smooth architectures like those based on triply-periodic minimal surfaces (TPMS) [12][13][14] address the aforementioned issue, but they do not overcome a common limitation in the fabrication of architectures based on periodic unit cells: classically, all truss, plate, and TPMS designs produce periodic structures with high sensitivity of mechanical properties to symmetry-breaking imperfections and defects and limited opportunity for introducing smoothly spatially variant structures.
In light of these challenges, spinodal metamaterials have emerged recently as a new class of non-periodic architected material [15][16][17][18][19][20].Their design is inspired by topologies observed during spinodal decomposition [21,22], which occurs, e.g., during rapid diffusion-driven phase separation in nanoporous alloys [23,24], polymer blends [25], and micro-emulsions [26].The computational design of spinodal microstructures relies on simulating the phase separation process by phase field methods whose kinetics are modeled by Cahn-Hilliard-type evolution equations [19,[27][28][29][30][31][32] or, as a shortcut, by Gaussian random fields [20].As a two-phase mixture spontaneously decomposes into two spatially-separated stable phases, the process is artificially arrested and solid-or shell-based topologies are extracted by, respectively, removing one of the two phases or by retaining the interfaces and eliminating both phases.The resulting topologies are composed of smooth, bi-continuous, and non-self-intersecting surfaces were shown to have intriguing properties.Portela et al. [16] demonstrated that spinodal topologies exhibit better (and close to optimal) stiffness scaling with respect to relative density, and they were shown to exhibit a significantly improved mechanical resilience over comparable truss and plate metamaterials, as seen in the recovery after repeated nonlinear compression.Unlike beam and plate architectures, spinodal metamaterials (like TPMS) avoid stress concentrations at beam or plate junctions, thus reducing stress concentrations, which considerably contributes to the observed high mechanical resilience.Compared to TPMS, the non-periodicity of spinodal architectures renders them insensitive to symmetry-breaking defects and fabrication-induced imperfections [17].Guell Izard et al. [18] further demonstrated ultra-high energy absorption characteristics of spinodal architectures.Spinodal metamaterials can also self-assemble across several length scales -from centimeters to nanometers [16], which is a promising avenue to overcome the scalability challenge of additive manufacturing.(Note that in this study we do not take into account the manufacturability of structures as, e.g., [33]).
We recently introduced a computational shortcut to generate spinodal-like topologies, referred to as spinodoid topologies [15].This approach replaces computationally expensive phase field simulations for topology generation and provides a simple parametrization based on anisotropic Gaussian random fields (GRFs) [34,35]; unlike the isotropic formulation of Soyarslan et al. [20], spinodoids allow for an efficient exploration of a wide design space of anisotropic mechanical properties.When designing functionally-graded metamaterials, the GRF-based approach admits seamless, spatially-variant topologies, which in contrast to periodic unit-cell-based designs bypasses discontinuities and tessellation-related limitations.
Spinodoid metamaterials bear potential for applications ranging from energy absorption and impact protection to heat exchange and to synthetic bone.An ongoing challenge in the design of patient-specific bone implants is to match the anisotropic topological and mechanical properties of bone -which can be highly heterogeneous across patients as well as within the same bone.Functionally-graded spinodoid metamaterials were shown to be promising candidates for inverse-designed synthetic bones [15] for improved biomechanical compatibility and reduced bone atrophy.Yet, aside from optimizing the properties of specific spinodoid topologies, they have not been used in any two-scale design challenge, such as, e.g., identifying the optimal macroscale shape of a bone implant while optimizing the local, spatially-varying spinodoid microstructure.To this end, we here address the systematic design of spatially-variant, functionally-graded bodies with a spinodoid microstructure through a data-driven topology optimization approach.
Topology optimization is a well-established technique (see [36] and [37] for detailed reviews).The classical Solid Isotropic Material Penalization (SIMP) method and its extensions [38][39][40][41][42] define a continuous volume fraction field ρ : Ω → [0, 1] over a body Ω ⊂ R d in d dimensions, such that the local linear elastic modulus tensor is approximated as where C 0 and C S are the stiffness tensors of void and solid regions, respectively, and p ≥ 3 is a penalization exponent to promote (approximately) purely void or solid states.Most SIMP-based methods optimize the material distribution within a macroscopic body, e.g., obtaining the fill-fraction field ρ(x) for all x ∈ Ω by minimization of the total compliance subject to given boundary conditions and loads.
Multiscale topology optimization additionally focuses on the microscale design (e.g., optimizing the architecture of a metamaterial, or the fiber orientation in composites) in a spatially-variant fashion along with the material distribution on the macroscale (unlike, e.g., multiresolution approaches which gain efficiency by introducing different mesh resolutions on the same (macro-)scale [43]).The two-scale optimization may be carried out sequentially [44,45] -e.g., identify optimal material stiffness tensor fields for the macroscale problem and then search the microstructural design space to meet the target properties using inverse homogenization methods [46,47].In practice, this approach faces challenges because the target material properties may be physically infeasible or unachievable by the chosen microstructural design space.As a remedy, simultaneous optimization at both scales (e.g., in an FE 2 setting) is more robust [48][49][50] but also computationally expensive.Here, we adopt the latter approach of simultaneous two-scale optimization in a new, computationally efficient fashion.
Topology optimization with anisotropic materials allows for manipulating the material orientation to generate structures with superior properties.In the context of compliance minimization problems, several works have shown that optimal orientations of anisotropic materials tend to align with the principal stress or strain directions [51][52][53][54][55][56][57][58][59].Within two-scale topology optimization, such alignment can be achieved by treating, e.g., the orientation angle(s) as additional design variables, which in composites is also known as continuous fiber angle optimization (CFAO) [60][61][62][63].In most such approaches, the inherent anisotropy of the base material is assumed constant.By contrast, recent works [33,[64][65][66] optimized the material anisotropy by tuning the microstructural architecture; yet, they did not consider the (spatially varying) orientation of the microstructure in a macroscopic body -primarily because such designs are based on periodic unit cells, which do not admit tessellations with arbitrary and/or spatiallyvariant orientations.For example, the optimized strut-based cuboidal unit-cells of Watts et al. [33] (while improving manufacturability) are always aligned with the Cartesian coordinate axes, which is sub-optimal compared to the case when the cells and the constituent struts are aligned locally with the principal stress or strain directions.Recently, Wu et al. [67] introduced conforming lattice optimization to address this issue; however, their strut-based design domain (rectangular or cuboidal cells with beams at the edges) strongly limits the microstructural and anisotropic tunability.Distinct from the above works with either fixed material anisotropy or fixed material orientation, our spinodoid topologies discussed here are simultaneously optimized for both the material anisotropy as well as the orientation distribution of the microstructure within a macroscopic body.
When designing metamaterials with spatially varying microstructure, multiscale topology optimization is computationally expensive, as on-the-fly homogenization methods using, e.g., the finite element method (FEM) must extract the effective properties at each material point on the macroscale from the underlying, local microstructure in each iteration of the optimization problem.Look-up tables are a convenient short-cut [68], yet they strongly limit the available design space and raise questions about interpolations between available data.More recently, machine learning (ML) techniques have attracted interest for accelerating topology optimization in two ways.First, by creating a dataset of solutions obtained using conventional multiscale topology optimization for several different boundary conditions, loads, microstructures, and material properties, a data-driven model can be learned for accelerated or even real-time prediction of the optimal solution as a function of these inputs without the need for an optimizer [69][70][71][72][73][74].However, since these approaches are essentially image-to-image learning (e.g., treating the boundary conditions or material distribution as multi-channel images), a large amount of training data is required to achieve reasonable accuracy.Additionally, generalization to unseen inputs (e.g., new loads or boundary conditions) is limited.The above methods are further challenging to extend to three-dimensional (3D) multiscale optimization with high-dimensional design parametrizations.
Second, topology optimization has been accelerated by employing homogenization surrogate models, which are typically developed through an offline training or interpolation of structure-property tables.For example, Watts et al. [33] developed polynomial surrogate models of homogenized elastic properties of open-truss micro-architectures for deployment in multiscale topology optimization.Others have used Gaussian processes [75], Numerical EXplicit Potentials (NEXP) [76,77], and single-layer neural networks (NNs) [66] in related contexts to approximate the effective material response to bypass expensive FEM simulations.However, complex design spaces like that of spinodoids (which have high-dimensional, highly-nonlinear, and multiply-connected parametrization domains) require surrogate models beyond simple interpolation methods.To this end, we employ deep NNs (DNNs), which have emerged as powerful tools to efficiently learn in high-dimensional spaces.In this contribution, a DNN pretrained on a structure-property dataset (obtained via FEM) is used to predict the anisotropic elastic stiffness of a given spinodoid topology.
In this study, we present a multiscale topology optimization framework for cellular structures based on spinodoid topologies, with simultaneous optimization of the macroscale material distribution and the microstructural (spinodoid) design and orientation, where the effective microscale response and associated sensitivities are provided by a data-driven surrogate model that replaces nested FE calculations on the microscale.We point out that the optimization problem requires computing the sensitivity of the compliance at a material point on the macroscale with respect to its design parameters.Lacking closed-form derivatives, this sensitivity analysis in our multiscale descripion requires numerical differentiation (ND), which involves perturbing the design parameters and recomputing the effective stiffness for each perturbation -requiring computationally expensive FE calculations.Moreover, its accuracy is strongly subject to round-off and truncation errors.While symbolic differentiation (SD) of the surrogate model can give exact derivatives, it leads to inefficient and redundant expressions when applied to complex nonlinear functions as is the case here.We address this issue by using automatic differentiation (AD), which is naturally supported by our NN-based surrogate model.Distinct from ND and SD, AD avoids the above limitations by creating a computational graph that tracks the series of mathematical operations between the inputs and outputs of an arbitrary function and applying chain rule recursively to compute its derivatives.Note that the derivatives are exact, even though AD does not provide the functional form of the derivatives like SD.An extensive introduction to AD can be found in [78].Examples of application of AD to design optimizations can be found in, e.g., [66,[79][80][81].
In the following, we present the topology generation and homogenization of spinodoids in Section 2. Section 3 defines the multiscale topology optimization problem, for which Section 4 introduces the DNN-based surrogate model and discusses its implementation in topology optimization.We use this model in Section 5 to present several benchmark examples of compliance minimization, before Section 6 concludes our study.

Spinodoid topology generation
Consider a homogeneous two-phase solution occupying a domain V ∈ R3 and undergoing spinodal decomposition.During the early stages, the separated phases are spatially-arranged in a stochastic bi-continuous topology, whose evolution is governed by the Cahn-Hilliard equation [21,27,35].Let ϕ : V → R be the phase field that describes the concentration fluctuation of one phase.Using Fourier analysis, Cahn showed that the phase field can be described by a GRF -a superposition of several plane waves with fixed wavelength but random wave vectors and phase shifts.That is, the concentration ϕ of either of the two phases at x ∈ V is given by where N is the number of waves, β > 0 is a constant wavenumber, and n i ∼ U(S2 ) and γ i ∼ U ([0, 2π )) denote the uniformly-distributed direction1 and phase angle of the ith wave vector, respectively.Without loss of generality, the amplitudes are assumed equal and constant to ensure that ϕ is normally distributed with zero mean and unitary standard deviation.Note that β prescribes a wavelength and hence directly controls the microstructural length scale.
The GRF in (2) describes the separated phases in the case of isotropic diffusion, whereas direction-dependent, anisotropic topologies are generally the result phase separation processes with directional preference in interfacial energy or diffusive mobility -which translates into a directional preference in the distribution of the wave vectors in (2).Here, we follow the anisotropic extension of Cahn's GRF solution presented in [15].Assuming the principle directions of mobility are aligned with the Cartesian basis {ê 1 , ê2 , ê3 }, the resulting anisotropic topologies can be approximated by a non-uniform orientation distribution function (ODF) parametrized by such that the wave vectors in (2) are restricted to lie within cone angles {θ 1 , θ 2 , θ 3 } from the orthogonal Cartesian basis vectors (see Fig. 1(a)).
For the purpose of topology generation, GRF models bypass the need for expensive simulations of the Cahn-Hilliard equation.Soyarslan et al. [20] showed that bi-continuous solid-void microstructures can be generated by applying a level set ϕ 0 to the phase field in (2), such that where ξ (x) = {1, 0} denotes the presence of solid and void, respectively, at x ∈ V. Since ϕ follows a standard normal distribution, the level set ϕ 0 is defined as the quantile evaluated at the average relative density ρ of the solid phase: . We note that, while we only consider solid networks for our multiscale topology optimization framework in the following, shell-type architectures can be generated in a similar fashion by choosing an isosurface of ϕ(x) = ϕ 0 .

Homogenization of the elastic stiffness
For a given set of design parameters, we use computational homogenization [82] via FEM to compute the effective mechanical stiffness of the corresponding spinodoid topology.We consider a cubic representative volume element (RVE), in which the GRF from (2) with level set ( 4) is used to generate a spinodoid architecture.The base material of the solid is assumed homogeneous, isotropic, linear elastic with Young's modulus E s and Poisson's ratio and ν s = 0.3 (E s is arbitrary as it scales linearly with the effective stiffness).For an RVE of size l × l × l, the wave number β = 10π/l was found to be sufficient for extracting converged homogenized stiffness values from the chosen RVE (larger RVE sizes at fixed β yielded only marginal variations in the computed effective stiffness values).Six numerical experiments are performed on the RVE, one uniaxial stretch and one simple-shear loading along each of the three principal axes.By applying average strains ⟨ε⟩ to the RVE for each of the six cases, the elastic stiffness tensor Ĉ is computed by solving the linear system of equations ⟨σ ⟩ = Ĉ ⟨ε⟩, where ⟨σ ⟩ = 1 |V| ∫ V σ dV is the volume-averaged stress obtained from FEM.Since the spinodoid topologies lack periodicity, we apply affine boundary conditions to the RVE and acknowledge that the computed response provides an upper bound to the actual effective stiffness (yet, spinodoids -like spinodals -retrieve their beneficial stiffness scaling properties from being stretching-dominated, so that affine boundary conditions can be expected to not affect the response as much as in, e.g., bending-dominated trusses).
The fourth-order stiffness tensor is visualized through the elastic surface, showing Young's modulus E(d) along all directions d ∈ S 2 , which is computed as Fig. 1 illustrates representative examples of the diverse anisotropic stiffnesses achievable by the spinodoid design space -from lamellar topologies (being highly soft in a dedicated direction) to columnar ones (stiff in two directions) to cubic and (approximately) isotropic topologies.In subsequent derivations, the strain, stress, and stiffness tensors will be used frequently in their respective Voigt notations, which we denote by subscript (•) v .
For further details on topology generation and homogenization of spinodoids, the reader is referred to [15].

Multiscale topology optimization
The overall multiscale topology optimization setup presented in the following is schematically summarized in Fig. 2.

Macroscale: boundary value problem and compliance minimization
The classical boundary value problem (BVP) of static equilibrium in a macroscale domain In our two-scale setup, we assume that the local elastic properties and mass density at a point x ∈ Ω on the macroscale stem from homogenization over a microscale RVE, whose architecture is defined by a local set of design parameters, χ (x).Assuming a separation of scales, the macroscale stresses σ and strains ε are related via the linear elastic constitutive relation where Ĉ represents the homogenized stiffness tensor obtained from an RVE with design parameters χ , so the spatial variability in microstructural topology is captured naturally by χ (x).Macro-and microscales are hence coupled by exchanging homogenized information such as χ and Ĉ.For the example of spinodoids, χ includes Θ = (ρ, θ 1 , θ 2 , θ 3 ) T among other design parameters which will be explained in Section 3.2.
We use FEM to solve the macroscale BVP (6) in a (Bubnov-)Galerkin setting, which discretizes the macroscale body Ω into a mesh containing n elements and n n nodes.In subsequent examples, we use constant-strain tetrahedra in 3D, having a single quadrature point per element.Of course, the following formulation can be analogously adopted for higher-order elements.Let X = {χ e , e = 1, . . ., n} be the set of all microstructural design parameters (one such set for each element e), such that the stiffness matrix of element Ω e is given by where B e denotes the strain-displacement matrix of element e [83].For given X , the global vector of nodal displacements U = {u i , i = 1, . . ., n n } is obtained from the following system of equations (assuming that essential boundary conditions have been imposed): The global displacement vector U, global stiffness matrix K , and the external force vector F are obtained via assembly of the element-wise displacement vectors u e , element stiffness matrices k e , and external forces (derived for each element from the applied tractions t), respectively.
The total compliance at the macroscale is given by where ( 8) and ( 9) have been substituted.The compliance minimization problem on the macroscale is formulated as where ρ prescribes a target on the overall relative density at the macroscale (effectively imposing a total weight constraint), and ( 9) ensures static equilibrium of the loaded structure.
The nonlinear optimization problem in ( 11) is solved using IPOPT [84], a primal-dual interior-point algorithm with a filter line-search method.IPOPT requires sensitivity information, i.e., the derivative of the compliance Φ with respect to the design parameters X , for determining the line search direction.Differentiating both sides of the equilibrium constraint (9) with respect to X and noting that K is symmetric, we obtain The compliance sensitivity hence follows as where (12) has been substituted to obtain the right-hand side.Using (10), the above may be expressed as an element-wise summation, viz.
To eliminate numerical instabilities such as checkerboard patterns, we impose a minimum length scale by a filtering technique.Commonly used filter techniques include density filters [85], sensitivity filters [86], and morphology-based filters [87].In our approach, the stiffness sensitivity in ( 14) is replaced by (no summation over i implied) with weight function where x ∈ Ω e , while x e denotes the position of the center of element e, and r filter > 0 is a cut-off radius that controls the length scale of the filter.The optimal solution of (11) yields the set of spatially-varying design parameters X , which can be used to generate the multiscale topology with fully resolved microstructures.Details on the computational generation of fully resolved topologies are presented in Appendix A.

Microscale: spinodoid microstructures
The spinodoid microstructure at every point x ∈ Ω on the macroscale depends on design parameters χ (x).As we are using simplicial elements (with a single quadrature point), we deal with element-wise constant χ e .Here, we consider the spinodoid topology within a single element and, for the sake of brevity, we drop the superscript (•) e in the following discussion, while tacitly implying that the microscale description applies for each element e ∈ {1, . . ., n}.For spinodoids, the microscale response is generally determined by design parameters χ = Θ = (ρ, θ 1 , θ 2 , θ 3 ) T .Unfortunately, this straightforward definition of χ presents technical challenges in the context of multiscale topology optimization, which are addressed here.
As discussed in Section 2, the spinodoid design space is restricted to ρ ≥ ρ min to avoid disjoint solid domains (at the microscale).With this constraint, however, the entire macroscale domain is filled with material (ρ > 0).As a remedy, we need to introduce macroscale holes by allowing ρ = 0 in the (macroscale) voids,2 while enforcing ρ ≥ ρ min in the non-void regions.A similar transformation was proposed by White et al. [66].Similar to the density ρ, the domain for angles θ 1 , θ 2 , θ 3 ∈ {0}∪[θ min , π/2] is also disjoint.From an optimization perspective, such disjoint parameter spaces are not favorable.Therefore, we transform the parameter space in a continuous fashion by defining the transformation (Fig. 3) Here, each transformation is a smooth approximation of a Heaviside step function (replacing the discontinuous design space by a C 0 -continuous function), which allows for continuous parameter spaces ρ ∈ [0, 1] and θ 1 , θ 2 , θ 3 ∈ [0, π/2] when replacing Θ by f (Θ) in the topology optimization problem.The smoothness, controlled by λ 1 and λ 2 , is necessary to ensure that gradient-based optimization methods such as IPOPT are applicable.The disallowed parameter values are approximately mapped to the closest allowed values.For example, using f (Θ) relaxes the bounds on the relative density to ρ ∈ [0, 1], as microstructures with ρ < ρ min are effectively mapped to void with zero stiffness.When weight or compliance minimization is the objective, this will force the density to be either zero (void) or greater than ρ min (spinodoid).
So far, the principal directions of the spinodoid topologies and their inherent mechanical anisotropy are aligned with the fixed Cartesian basis {ê 1 , ê2 , ê3 }.In this case, the effective stiffness Ĉ in (7) is the homogenized stiffness of the microstructure (Section 2.2).As the stiffness is direction dependent, we here further expand the design space by applying rigid-body rotations to the generated spinodoid topologies.For demonstration purposes (and since our subsequent benchmarks are in 2D), we restrict ourselves to rotations about a single axis, while noting that the framework can be extended to general 3D rotations.We introduce an expanded set of design parameters where Ĉv is the effective stiffness passed to the macroscale in (7), and C v ( f (Θ)) is obtained from microscale homogenization of the spinodoid topology with design parameters f (Θ).The components of rotation matrix T are given by The homogenized stiffness and its sensitivity, which involves computing ∂C v ( f (Θ))/∂Θ, are the computational bottlenecks (often computed by ND), especially since IPOPT is an iterative algorithm which requires computing the sensitivity of each element several times during the optimization.This motivates the use of data-driven surrogate models that can accelerate the optimization process.

Micro-to-macro transition
As in classical homogenization, we assume a separation of scales between the spinodoid architectures on the microscale and the characteristic lengths over which fields of interest vary in the macroscale boundary value problem.The spinodoid architectures introduced above are in general scale-independent (as long as no material-level size effects emerge), meaning that the size of the spinodoid features is irrelevant in the homogenization framework and can be chosen arbitrarily a posteriori when generating the spatially variant structures.Specifically, relative density and effective elastic properties are independent of the specific length scale of the spinodal features, which is governed by parameter β (the characteristic wavelength of the spinodal microstructures is 2π/β).Once an optimal two-scale topology has been identified in any of the subsequent benchmark problems, the spatially-variant structure can be realized with, in principle, arbitrary values of β, resulting in finer or coarser microscale features, as needed (and as limited by fabrication constraints).We point out that for spinodoid topologies whose features are not orders of magnitude smaller than the macroscale features, the chosen scheme, strictly speaking, presents an abuse of the homogenization assumption since a separation of scales may break down.Hence, for any given application, the choice of β should be a compromise between accuracy, efficiency, and manufacturability of the optimized structures.

Bypassing FEM-based homogenization
To bypass the computational expense of repetitive FEM simulations required for the microscale homogenization, we use a DNN-based surrogate model for the on-demand prediction of the homogenized stiffness tensor Ĉ. Due to the symmetries in the ODF defined in (3), the anisotropic stiffness tensor C is orthotropic with nine independent elastic moduli, which can be encoded into S = ( Ĉ1111 , Ĉ1122 , Ĉ1133 , Ĉ2222 , Ĉ2233 , Ĉ3333 , Ĉ2323 , Ĉ3131 , Ĉ1212 ) T .The DNN can be interpreted as a composition of several linear and nonlinear transformations that provide a map F ω : Θ → S from the spinodoid design parameters Θ to the compressed representation of the stiffness C.
We use the DNN architecture previously developed by Kumar et al. [15], which translates to the composite transformation Here, L i→ j ω k denotes the kth linear layer parametrized by the set of weights and biases ω k = { A k , b k } -and we write ω = {ω i , i = 1, . . ., 7} -such that any z ∈ R i is transformed according to R(•) = max(0, •) is the rectified linear activation unit (ReLU) that acts element-wise on the input and introduces nonlinearity to the series of linear transformations.The design parameters Θ are mapped into intermediate highdimensional spaces, where the surrogate model can efficiently and accurately learn the nonlinear relationship to the homogenized stiffness values S.This is particularly advantageous as the stiffness space is highly disjoint because of the discontinuities in the design space (recall that θ 1 , θ 2 , θ 3 ∈ {0} ∪ [θ min , π/2]); e.g., the effective stiffnesses of the lamellar and columnar topologies are quite different (cf.Fig. 1).
The above surrogate model is trained on a dataset { (Θ i , S i ), i = 1, . . ., n train } containing n train pairs of randomly sampled spinodoid design parameters Θ and their corresponding homogenized stiffnesses S (obtained via FEM).The optimal DNN parameters are obtained using back-propagation and the Adaptive Moment Estimation (Adam) [88] to minimize the mean squared prediction error over the training data Once trained, the DNN-based surrogate provides accurate stiffness predictions instantly compared to several minutes for FEM-based homogenization (for the meshes used in our study).As the focus here is on its integration into the multiscale topology optimization framework, we refer the interested reader to [15] for further details pertaining to the architecture, training, and dataset of the DNN surrogate model.

Bypassing ND for sensitivity computations
The DNN-based surrogate model beneficially provides the sensitivity of the stiffness with respect to the design parameters -an essential information for the multiscale topology optimization.Conventional approaches based on ND that perturb the design parameters are computationally expensive (requiring extensive FEM-based homogenization simulations) and prone to stability and precision issues.In contrast, our DNN-based surrogate model admits inexpensive computations of the sensitivity information, i.e., F ′ ω [Θ] = ∂ S/∂Θ, which can be rearranged tensorially into ∂ Ĉv /∂Θ.Since each transformation in the DNN (20) is differentiable almost everywhere, 3 the sensitivity is easily computed by propagating the gradients from S to Θ via the chain rule for differentiation.From an implementation perspective, this is efficiently achieved by AD, which tracks the computational graph from Θ and S. The back-propagation of gradients via AD is implemented in most contemporary software for NNs, as it forms the backbone for training NNs using gradient-based optimization methods such as Adam.Note that the gradients computed by AD are exact and hence do not experience stability or precision issues.

Performance validation of the surrogate model
Before applying the trained surrogate model for the homogenized stiffness to topology optimization problems, we first assess the accuracy of the surrogate model.For a quantitative assessment, the DNN is trained on n train = 19,170 pairs of design parameters and their corresponding effective elastic stiffnesses from FEM, i.e., (Θ, S), followed by testing on an independent test set of containing n test = 2,130 pairs.Ideally, the stiffness prediction F ω [Θ] from the DNN should agree with the true stiffness, i.e., F ω [Θ] ≈ S. Thus, in a plot of predicted vs. true stiffness (Fig. 5), the predictions are expected to lie on a line with zero-intercept and unit-slope.The trained DNN achieves an R 2 ≥ 0.996 accuracy for each stiffness component.We also perform a similar assessment (Fig. 6) for the stiffness sensitivities where we use the central finite difference scheme applied to FEM-computed stiffnesses to validate the derivatives obtained from the DNN via AD.The AD approach shows good agreement with R 2 ≥ 0.999 for each stiffness component.

Benchmarks
We investigate three design examples -a cantilever beam, an L-shaped structure, and two designs with simultaneous consideration of multiple load cases.As a baseline, for each example we also present the optimal topologies obtained from the classical SIMP method for a homogeneous and isotropic material with Young's modulus E s and ν s = 0.3 (which are the same values as those of the base material for spinodoids).The penalty factor for SIMP is chosen as p = 4.To reduce numerical artifacts due to checkerboard patterns, the filter radius r filter in ( 15) is set to 0.075 for all subsequent simulations (thresholding is applied at the end to achieve pure 0-1 solution).For simplicity, we avoid units in the following benchmarks by setting E s = 1 and normalizing all lengths.To avoid inadmissible design parameters, we choose λ 1 = 600 and λ 2 = 60 × (180/π) rad −1 in (17) (see Fig. 3).
In subsequent examples (with the exception of Benchmark IV), the small thickness of the chosen design domains (relative to all other dimensions) and the symmetry in boundary conditions results in uniform distributions of the optimal design parameters across the thickness of each design, thus effectively producing 2D topologies on the macroscale with constant thickness.

Benchmark I: Cantilever beam
We consider a cantilever beam (Fig. 7(a)) of size 1.5 × 1.0 × 0.1 discretized by a mesh of linear tetrahedral elements on a 64 × 48 × 2 uniform grid.The left vertical face of the beam is clamped (i.e., displacements are suppressed in all directions).A single point load of magnitude 2.0 × 10 −2 in the downward direction is applied at the center of the right vertical face.We seek to minimize the compliance subject to an average relative density constraint of ρ = 0.5.As an initial guess, all elements are initialized with spinodoid relative density ρ = 0.5, anisotropy parameters (θ 1 , θ 2 , θ 3 ) = (15 • , 0 • , 0 • ) and orientation α = 0 • .The initial design is intentionally chosen to be lamellar-type (with all lamellae normal to length of the beam).Such a design will strongly deform under axial loads (as evident from their low Young's moduli; see Fig. 1(b)) and therefore show a highly compliant behavior.
The optimal design with anisotropic spinodoid architectures is shown in Fig. 8.The material distribution (Fig. 8(a)) resembles the optimal topology obtained via SIMP (Fig. 7(b)), albeit the former is characterized by larger regions with intermediate density.The optimal compliances obtained by the single-scale SIMP design and our multiscale spinodoid design are shown in Fig. 7(c)).Contrary to the initial guess of a lamellar topology, the final design is dominated by cubic topologies with larger Young's moduli along the principal directions.Throughout the macroscale body, the spinodoid microstructures are rotated, so that their preferred orientations follow the material distribution (see Figs. 8(a) and 8(e)).Fig. 9 illustrates the seamlessly spatially-variant spinodoid topology (with fully resolved microstructure), which bypasses the challenge of incompatible microstructures in periodic metamaterials.To avoid the high computational expense of generating a large mesh with fully resolved microstructure, the microscale topology is illustrated for a representative subdomain only.

Benchmark II: L-shaped structure
In this example, an L-shape structure with the dimensions shown in Fig. 10(a) is optimized.The top face is fixed in all directions, while a uniformly-distributed vertical load is applied on the lower right edge, as shown.We seek to minimize the total compliance subject to an average relative density constraint of ρ = 0.5.As an initial guess, all elements are initialized with a spinodoid relative density ρ = 0.5, anisotropy parameters (θ 1 , θ 2 , θ 3 ) = (35 • , 15 • , 15 • ) (cubic topology), and orientation α = 0 • .The domain is discretized into a finite element mesh with 22,060 linear tetrahedral elements, yielding a total of 110,300 design variables.
The optimal topology with anisotropic spinodoid architectures is shown in Fig. 11 and resembles the SIMP result (Fig. 10(b)) in terms of the material distribution.The spinodoid-based design achieves 8.35% (Fig. 10(c)) improvement in the compliance relative to the SIMP-based design.This is possible due to the spatially-varying anisotropy; e.g., the zoomed-in microstructures in Fig. 11(a) show columnar features, which provide high stiffness in the principal stress directions -thus orienting material not only at the macroscale but also on the microscale for optimal compliance.Fig. 10(c) lists the optimal compliance for five different mesh resolutions.Increasing the number of finite elements does not affect the optimized compliance significantly.Counter-intuitively, higher mesh resolution may not necessarily yield a lower compliance in the context of a multiply-connected design parameter   space and non-convex property space -e.g. for spinodoids, ρ ∈ {0} ∪ [ρ min , 1] and θ 1 , θ 2 , θ 3 ∈ {0} ∪ [0, π/2].Similar observations were reported and explained previously by, e.g., [89] and [90].Benchmark III: (a) Schematic of the structure to be optimized for two symmetric load cases (shown by the two sets of applied forces).(b) Front view of the optimal topology obtained using the SIMP method.(c) Comparison of the optimal compliance with spinodoid microscale architectures (via the proposed method for different mesh resolutions) and solid material (via SIMP).

Benchmark III: Multiple load cases -symmetric loading
Here, we extend our multiscale topology optimization algorithm for compliance minimization to the simultaneous consideration of multiple load cases.The objective function in (10) is modified as the sum of the compliances for each individual load case, i.e., where the superscript (•) i denotes the ith load case, and M is the total number of load cases considered.
As an example, we consider the beam shown in Fig. 12(a) with M = 2 load cases: a uniformly-distributed force in the upward direction on the upper-right edge, and a uniformly-distributed force in the downward direction on the lower-right edge.The average relative density constraint remains ρ = 0.5.As an initial guess, all elements are initialized with a spinodoid relative density ρ = 0.5, anisotropy parameters (θ 1 , θ 2 , θ 3 ) = (35 • , 15 • , 15 • ) (cubic topology), and orientation α = 0 • .The domain is discretized into a uniform linear tetrahedral mesh with 25,350 elements and 126,750 design variables (a convergence study (Fig. 12(c)) revealed that increasing the mesh resolution leads to only marginal changes in accuracy; e.g., increasing the number of elements 23% (from 25,350 to 31,104) leads to a change in the optimal compliance of approximately 0.2%).
The optimized design obtained for spinodoid microstructures and SIMP are shown in Figs. 13 and 12(b), respectively.Similar to previous benchmarks, the former achieves an improvement of 4.88% in compliance (Fig. 12(c)) over the SIMP-based design due to the spatially-varying material distribution, anisotropy, and orientation on the microscale.Benchmark IV: (a) Schematic of the structure to be optimized for three non-symmetric load cases.(b) 3D view of the optimal topology obtained using the SIMP method.Gray regions here denote complete void.(c) Comparison of the optimal compliance with spinodoid microscale architecture (via the proposed method) and solid material (via SIMP).

Benchmark IV: Multiple load cases -non-symmetric loading
We modify Benchmark III (Section 5.3) by considering the three non-symmetric load cases shown in Fig. 14(a), which involve three point loads of the same magnitude but at three different corners (M = 3).An average volume constraint of ρ = 0.4 is imposed.The domain is discretized into a uniform linear tetrahedral mesh with 48,000 elements and 240,000 design variables.
Optimal SIMP-as well as spinodoid-based designs are illustrated in Figs.14(b) as well as in Figs. 15 and 16, respectively.Compared with the previous benchmarks, the spinodoid-based design shows only minor improvement (1.27%) in compliance over SIMP (Fig. 14(c)).This can, in fact, be expected for the simultaneous optimization for multiple load cases applied in different directions, since -under such constraints -the structure is likely to favor isotropic topologies as the best compromise between all design cases, thus leading to similar performance as SIMP.We note that the small improvement of 1.27% in compliance further needs to be considered with caution, as it can partially be due to numerical artifacts.

Conclusions
We have presented a two-scale topology optimization framework for macroscopic bodies made of a spatiallyvariant microscale architecture based on spinodoid topologies.Inspired by microstructures produced by spinodal decomposition, the spinodoid topologies are anisotropic and described by a set of four design parameters (the relative density and the orientation distribution of wave vectors in the underlying Gaussian random field).The topology optimization problem minimizes the linear elastic compliance of macroscopic bodies by solving both for the displacement field and for a continuous field of the design parameters.The effective material response at any point on the macroscale is identified with the homogenized, effective response of a representative volume element filled by the spinodoid topology defined by the local design parameters.To bypass costly computational homogenization simulations, we here introduce a new approach based on a deep neural network as a surrogate model that maps the design parameters onto the effective fourth-order stiffness tensor.Aside from significantly speeding up calculations, this approach also provides exact sensitivities (required for gradient-based optimization) at low numerical costs.Although we here focused on solid spinodoid topologies with relative densities larger than

Fig. 1 .
Fig. 1.Spinodoid topology generation and representative examples.(a) Schematic illustration of the three cone angles θ 1 , θ 2 , and θ 3 , which describe the anisotropic distribution of wave vectors n i .(b-e) Examples of spinodoid topologies with their respective anisotropic elasticity surfaces (and an illustration of the cone angles indicated in each caption).Each point on the elasticity surface denotes the Young's modulus in the corresponding direction.Light and dark gray spheres indicate, respectively, the elastic Voigt upper bound and the Hashin-Shtrikman upper bound (the latter being applicable only to the isotropic structure).Source: Adapted from [15].

Fig. 5 .
Fig. 5. DNN-predicted vs. true components of Ĉ in the test dataset.Dashed lines represent the ideal line with zero-intercept and unit-slope; the corresponding R 2 -deviations are indicated.

Fig. 6 .
Fig. 6.Ĉ1111 stiffness sensitivity computed via automatic differentiation vs. numerical differentiation in the test set.All dashed lines represent the ideal line with zero-intercept and unit-slope; the corresponding R 2 -deviations are indicated.

Fig. 7 .
Fig. 7. Benchmark I: (a) Schematic of the end-loaded cantilever beam whose compliance is to be optimized.(b) Front view of the optimal topology obtained using the SIMP method.(c) Comparison of the optimal compliance for the spinodoid architecture on the microscale (via proposed method) and a solid material (via SIMP).

Fig. 8 .
Fig. 8. Benchmark I: Front view of the optimal topology (material distribution, anisotropy, and orientation) of the cantilever beam with spinodoid microscale architectures.The spatially-varying design parameters are illustrated across the body, and two examples of distinct microscale topologies are shown in (a).

Fig. 9 .
Fig. 9. Benchmark I: Seamlessly spatially-variant spinodoid topology with fully resolved microstructure in a subdomain of the optimized cantilever beam design.The microscale GRF is generated with wavenumber β = 1600π/3 (arbitrarily chosen and sufficiently large under the assumption of an approximate separation of scales).Details on generating the microscale topology are presented in Appendix A.

Fig. 10 .
Fig.10.Benchmark II: (a) Schematic of the L-shaped structure whose compliance is to be optimized.(b) Front view of the optimal topology obtained using the SIMP method.(c) Comparison of the optimal compliance with spinodoid microscale architecture (via the proposed method for different mesh resolutions) and solid material (via SIMP).

Fig. 11 .
Fig. 11.Benchmark II: Front view of the optimal topology (material distribution, anisotropy, and orientation) of the L-shaped structure with spinodoid microscale architecture.Examples of the spatially-varying design parameters and the corresponding microscale topologies are shown in (a).

Fig. 12 .
Fig.12.Benchmark III: (a) Schematic of the structure to be optimized for two symmetric load cases (shown by the two sets of applied forces).(b) Front view of the optimal topology obtained using the SIMP method.(c) Comparison of the optimal compliance with spinodoid microscale architectures (via the proposed method for different mesh resolutions) and solid material (via SIMP).

Fig. 13 .
Fig. 13.Benchmark III: Front view of the optimal topology (material distribution, anisotropy, and orientation) of the spinodoid-based structure, optimized for two symmetric load cases.Examples of the spatially-varying design parameters and the corresponding microscale topologies are shown in (a).

Fig. 14 .
Fig.14.Benchmark IV: (a) Schematic of the structure to be optimized for three non-symmetric load cases.(b) 3D view of the optimal topology obtained using the SIMP method.Gray regions here denote complete void.(c) Comparison of the optimal compliance with spinodoid microscale architecture (via the proposed method) and solid material (via SIMP).

Fig. 15 .
Fig. 15.Benchmark IV: Optimal material distribution in the spinodoid-based structure with multiple non-symmetric load cases.Front-view of multiple cross-sectional planes is shown here.

Fig. 16 .
Fig. 16.Benchmark IV: Optimal topology (material distribution, anisotropy, and orientation) of the spinodoid-based structure optimized for three non-symmetric load cases.The front view of the cross-section at x 3 = 0.2 is shown here.Examples of spatially-varying design parameters and the corresponding microscale topologies are shown in (a).