Topology optimization of differentiable microstructures

approach the theoretical limit.


Introduction
In recent years, alongside the advancements in additive manufacturing technologies, there has been a growing interest in the optimal design of multi-scale structures, i.e., structures that are composed of intricate sub-structures.Multi-scale structural optimization dates back to the seminal work of Bendsøe and Kikuchi [1].A comprehensive overview of the historical and latest developments in this field is available in a recent review article [2].Multi-scale structures are intrinsically lightweight and hold the promise of achieving superior functionalities, e.g., buckling resistance [3,4] and damage tolerance [5].Within a macroscale structural layout, the distribution of microstructures typically falls into one of two categories: periodic or spatially graded.In the former case, a single microstructure is repeated in the macroscale structural layout, resulting in a homogeneous distribution.In the latter, a set of varying microstructures is distributed within the macroscale structure, creating a spatial gradation of mechanical properties.The tailoring of local microstructure configurations to meet spatially varying functional requirements is of great importance in engineering design.For instance, in the design of orthopedic implants, it has been suggested to place highly porous microstructures at the bone-implant interface for promoting bony ingrowth, and allocate dense microstructures in the central region (where bony ingrowth is irrelevant) for maintaining structural integrity [6,7].
A systematic approach for the inverse design of microstructures is topology optimization [8,9].It involves optimizing the material distribution within a unit cell to tailor its effective properties, which can be determined through methods such as numerical homogenization.This inverse design approach, often referred to as inverse homogenization, was initially proposed to design periodic microstructures [8] and has since been extended to design a series of varying microstructures (e.g., [7,10]).For example, by applying inverse homogenization to maximize the bulk modulus of microstructures under a linearly varying volume fraction of solid material, it becomes possible to design a series of microstructures with varying porosity and bulk modulus.However, optimized microstructures that are close with respect to the volume fraction do not necessarily match when placed adjacently.This is due to the existence of numerous local optima in the tremendous solution space of microstructure topology optimization.The disconnectivity issue has received significant attention in the research community [7,[10][11][12][13][14][15][16][17][18].Zhou and Li [10] were among the first to investigate this problem.They proposed to optimize the summation of mechanical properties of individual cells, and introduced a non-linear diffusion term into the objective function to penalize disconnection.In a recent development, Garner et al. [7] proposed to optimize the mechanical properties measured on the extended domain of compound cells.This formulation stimulates mechanical compatibility across adjacent cells, since a poor connection is reflected by the suboptimal mechanical performance of the compound cells.
The set of graded microstructures obtained through inverse homogenization is typically parameterized by the volume fraction of solid material (or equivalently, porosity).Since inverse homogenization is applied to a discrete set of volume fractions, either individually or holistically, the obtained set of microstructures is inherently discrete.In other words, there is a finite difference in volume fractions between adjacent microstructures.This discrete nature implies abrupt changes in geometry between adjacent microstructures, even when the connectivity has been optimized.Consequently, stress concentration at the interface between adjacent microstructures may arise under loading conditions.In contrast to the discrete nature of the series of optimized microstructures, two-scale optimization with graded microstructures typically assumes a continuous parameter (e.g., the volume fraction of solid material) for the macroscale optimization problem [19,20].Thus, quantization of the continuous parameter is required for mapping the discrete microstructures to the optimized continuous macroscale layout.The error introduced by the quantization is often neglected.A strategy to reduce the quantization error is to increase the size of the discrete set of volume fractions.For instance, Garner et al. [7] demonstrated simultaneous optimization of up to 100 unit cells of varying volume fraction.This involves intensive computation, since numerical homogenization must be applied to each volume fraction.Alternatively, the quantization error can be mitigated by regarding the discrete microstructures as different materials in the macroscale optimization [21,22].This requires a multi-material optimization approach on the macroscale [23,24], which is more complex than optimizing with a single continuous parameter.
In this paper, we introduce the concept of differentiable microstructures, which entails a continuous parameterization of microstructures.In this framework, any parameter value within a valid interval corresponds to a unique microstructure instance.This innovative approach eliminates the need for quantization in two-scale optimization and removes the requirement for a multi-material optimization setup.Furthermore, the continuous parameterization ensures that an infinitesimal change in the parameter value results in an infinitesimal change in the geometry and mechanical properties of the microstructure.This theoretically eliminates the possibility of stress concentration occurring between adjacent microstructures.
Several known families of microstructures can be classified as differentiable microstructures.For instance, simple 2D X-shaped lattices parameterized by the in-plane thickness of bars, parameterized auxetic metamaterials [25], and 3D cellular structures based on parametric surfaces such as triply-periodic minimal surfaces [26].In fact, these and similar geometric parameterizations are not uncommon in two-scale structural optimization (e.g., [27][28][29][30], reviewed in [2]).These geometries offer a continuous parameterization, but at the price of a compromise in the mechanical properties of individual microstructures.The bulk modulus of these simply parameterized microstructures is lower than the theoretic limits [31].Differentiable microstructures may also be obtained by geometric interpolation.Cramer et al. [19] optimized isotropic microstructures for maximum bulk modulus at five volume fractions, and proposed the use of a signed distance field function for interpolation.This geometric approach works for microstructures of similar topology, which, however, is not always the case in inverse homogenization.Zobaer and Sutradhar [32] presented an optimization-based approach to interpolate orthotropic microstructures through supershape-based parametric modeling.Nevertheless, the complexity of microstructures that can be accurately represented by supershapes is limited.
In this paper, we present a novel formulation for designing differentiable microstructures through topology optimization.Our method takes a microstructure with a low volume fraction as the input and generates a series of microstructures by expanding it.To ensure the differentiability of microstructures, our approach starts with a continuous parameterization of microstructures.Specifically, we represent the continuous series of 2D microstructures by a height field defined over the microstructure design domain, where each height corresponds to a specific microstructure.We analyze the geometric feature of the height field and its implications on the mechanical properties of the resulting microstructures.In particular, it is found that the height field shall be smooth, and free of local maxima and saddle points to avoid sub-optimal microstructure configurations.To achieve this, we propose a novel regularization technique using a heat equation.We further analyze the optimality of the differentiable microstructures regarding different input microstructures, and provide heuristic principles to generate these input microstructures under specified conditions.Our method is demonstrated on the maximization of bulk modules of 2D microstructures.
The remainder of the paper is organized as follows.Section 2 introduces the definition of and the conceptualization towards differentiable microstructures.Section 3 presents the methodology for topology optimization of differentiable microstructures.Section 4 provides a comparative analysis of our results in relation to alternative approaches, and discusses the geometric features and differentiability of the optimized outcomes.The paper concludes with final remarks and suggestions for future research in Section 5.

Differentiable microstructures
Let  be a binary field, representing a microstructure that is defined on a domain  ⊂ R  ( = 2 or 3).At each point  ∈ ,  is 0 or 1, indicating that the material at this point is void or solid.We consider a series of varying microstructures that is parameterized by one or a few continuous parameters,  ∈ [, ].This series of varying microstructures, i.e., functionally graded microstructures, is compactly represented by ().
For multi-scale design optimization, it is essential to require that () are continuous functionally and geometrically.Here function refers to physical properties of (), e.g., its bulk modulus and average density.Functional continuity enables the use of gradient-based mathematical programming for solving multi-scale design optimization problems.Geometric continuity ensures that the functionally graded microstructures, when adjacently placed, can effectively deliver their physical properties, i.e., to transfer mechanical loads.Based on these requirements, we propose the concept of differentiable microstructures.
Differentiable microstructures.Differentiable microstructures are graded microstructures that are described by continuous parameters, and wherein both physical properties and the geometry of () are continuous and differentiable with respect to .Denoting a generic physical property by , the physical differentiability is related to the existence of lim →0   , where  is an infinitesimal change of the continuous parameter .The physical properties can be evaluated, for example, by numerical homogenization.To define geometric differentiability, a metric is needed for quantifying the geometric difference between () and ( + ).To this end, here we use the size of the absolute difference between the two fields, The geometric differentiability thus translates to the existence of lim →0   .We note that the metric of the geometric difference is independent of the topology of () or ( +).This is beneficial to allow for varying topologies in functionally graded microstructures.While often a small modification in geometry leads to a small change in physical properties, the change in physical properties is not guaranteed to be small.If the topology of ( + ) differs from that of (), the physical properties may differ unproportionally.Thus, both geometric and physical differentiabilities are needed to define differentiable microstructures.

Conceptualization towards differentiable microstructures
An approach to differentiable microstructures is to start with a simple geometry (e.g., an X-shaped microstructure in 2D), and uniformly change its thickness, e.g., by offsetting its boundary (shown in Fig. 1(a)).This achieves geometric and physical differentiability, but the physical property of the microstructures is sub-optimal.
To approach the optimal physical property, the main idea of our method is to non-uniformly morph the boundary of an initial microstructure.The uniform offsetting and non-uniform morphing are illustrated in Fig. 1(a) and (b) for a simple X-shaped input microstructure, respectively.The series of 2D microstructures can be seen as the cross-sections of a height field (cf. the bottom row in Fig. 1).The continuous morphing ensures geometric differentiability, while the non-uniform characteristic provides an opportunity to maximize the physical property for the entire sequence during morphing, (),  ≤  ≤ .The objective function is generically written as Here  is an operator for evaluating the physical property of interest, e.g., numerical homogenization.We note that, mathematically, maximizing the integral function does not guarantee the optimality of each individual microstructure.In our approach and in multi-scale optimization, the interest is to get an optimal series, and thus the integral form of the objective function.
In this work, we take the bulk modulus as an example.The input microstructure  0 corresponds to one with a low volume fraction.The optimization thus non-uniformly expands the boundary of  0 .During the expansion, the topology of the microstructure is allowed to change, offering flexibility to maximize the bulk modulus of the sequence.
We introduce a novel parameterization for encoding the non-uniform expansion.We conceptualize the expansion of the boundary as heat conduction in the design domain .The non-uniform expansion can thus be encoded by the temperature distribution resulting from the non-homogeneous thermal properties in the domain.Specifically, heat conduction is governed by where  is the thermal diffusivity. and  refer to temperature and time, respectively.The domain boundary () is insulated, i.e., (∇ ⋅ ⃖ ⃗ )|  = 0 where ⃖ ⃗  is the normal vector of boundary .The solid region of the input microstructure ( 0 = { |  0 () = 1}) serves as the heat source, with a constant temperature, i.e.,  |  0 ,= 0 = 1.The rest of the domain starts with a temperature of 0, i.e.,  | ⧵ 0 ,= 0 = 0.The temperature distribution at a specific time  =  1 is thus determined by the thermal diffusivity .We denote the temperature field at  1 by .The iso-contours of  are analogously interpreted as boundaries of varying microstructures.As a result of the heat conduction process, the obtained microstructures are geometrically differentiable.
This conceptualization is illustrated in Fig. 2. Given an input microstructure  0 and an initial thermal diffusivity field , a temperature distribution is computed by solving the heat equation.The iso-contours of the temperature field are indicated by black  lines in Fig. 2(c).These iso-contours are interpreted as varying microstructures, shown in (e-g).The temperature distribution can also be visualized by a height field (d), where the height corresponds to a temperature value.In the thermal diffusivity field for this illustration, the right-hand side of the domain has larger thermal diffusivities than the left-hand side.Consequently, the resulting microstructures have larger solid regions on their right-hand side (e-g).In order to generate optimal differentiable microstructures, the thermal diffusivity field is updated in an iterative numerical optimization process.The optimized microstructures are visualized in (h-l).
The height field  depends on the input microstructure  0 and the thermal diffusivity field .To obtain optimal microstructures, we take the heterogeneous thermal diffusivity field as design variables.The input microstructure is obtained by inverse homogenization to maximize its bulk modulus, at a small volume fraction.The inverse homogenization approach is able to find microstructures whose bulk moduli are close to the theoretical limit [7,8,33].The solution by the numerical approach however is not unique -there exist many microstructures with a similar bulk modulus but different topology and geometry.In the results section, we thoroughly examine the influence of geometric features of the input microstructure on the optimality of the optimized differentiable microstructures.

Topology optimization of differentiable microstructures
The design of differentiable microstructures is formulated as a topology optimization problem based on a finite element discretization of the domain .The integral objective (Eq.( 2)) is rewritten as the summation over  sample microstructures in the series, Here,  is a function to evaluate the physical property of an individual microstructure. [] represents the th sample microstructure, corresponding to   ∈ [, ].
As explained, the design variables are the thermal diffusivity field.The design variables  are first smoothed by a convolution filter, as commonly done in density-based topology optimization.The smoothed field α is then used to obtain a temperature field , by solving the heat equation.From the continuous temperature field, smoothed Heaviside projections with varying thresholds are applied to get the sample microstructures,  [] ,  = 1, … , .
This section starts with the computational steps to obtain the sample microstructures from the design variables.We then proceed to present the optimization problem and the constraints.The sensitivity analysis and pseudocode to solve the design optimization problem are also included.

Computational steps to obtain sample microstructures
Step 1: Smoothing of design variables  In numerical optimization, the integral objective (Eq.( 2)) is replaced by the summation over a few sample microstructures (Eq.( 4)).This means that the microstructures between sample microstructures are not explicitly part of the objective function.To ensure a continuous change between the sample microstructures, it is thus necessary to have a smooth thermal diffusivity field.The smoothness of the thermal diffusivity field affects the gradient of the resulting temperature field.To this end, a convolution filter commonly used in density-based topology optimization is adopted here to smooth the thermal diffusivity field.The smoothed diffusivity is where   is the volume of an element, which is constant in our implementations.The weighting factor decreases linearly with respect to the distance between   and   , where  is a predefined filter radius.  and   are respectively the centroid of element  and of an element that is in its close vicinity, Step 2: From the thermal diffusivity field α to a temperature field  The heat equation is used to generate a smoothly varying temperature field.In 2D, it is written as To solve this partial differential equation, the domain  is discretized into bi-linear quad elements.The temperature ( ) is defined on the nodes.Each element is associated with a thermal diffusivity value α .The left-hand side of Eq. ( 7) can be discretized using finite difference [34] or finite element schemes [35].A finite element scheme is used in this work for its applicability to irregular grids.The right-hand side of Eq. ( 7) is discretized using a one-step finite difference scheme.This is in accordance with the purpose of getting a smoothly varying temperature field.Note that we set the solid region of the input microstructure ( 0 = {| 0 () = 1}) as the heat source with a constant temperature, i.e.,  |  0 = 1, and the domain boundary is insulated.As  → ∞, the solution to Eq. ( 7) is a uniform temperature distribution.Using a finite element scheme for the space and a one-step finite difference scheme for the time dimension, the heat equation is written as is the global thermal diffusivity matrix, assembled from that of each element, i.e., ( α) = ∑    ( α ).The thermal diffusivity matrix for an element with α is linearly scaled from that of the element with a thermal diffusivity of 1, which is denoted by  =1 , We rearrange Eq. ( 8) as Here,  is an × identity matrix, where  is the number of nodes of the spatial discretization.Solving this linear system of equations to obtain the temperature distribution at nodes, the temperature of each element is then computed by bi-linear interpolation.This is written in the matrix form as a transformation by is a sparse matrix that transforms the temperature from nodes to elements.For structured quad meshes in 2D, as used in our examples,  , = 1 4 if the th node is one of the th element's nodes, and  , = 0, otherwise.
Step 3: From the temperature field  to sample microstructures  []  The temperature distribution is a continuous scalar field with values between 0 and 1.It can be visualized as a height field (see Fig. 2(d)).The cross-section of the height field by a horizontal plane is interpreted as a microstructure.In the microstructure corresponding to a cross-section at height   , a point at  is solid if   ≤   , or empty, otherwise.We use a smoothed Heaviside projection to approximate the thresholding, Here,  is a parameter to control the sharpness of the projection function [36].As  approaches infinity, the projection is equivalent to a step function.A continuation is adopted for stability.The initial value of  is set to 2. It is doubled every 50 iterations, until it reaches 128.

Optimization problem
With the encoding of the microstructures through a heat equation, the optimization problem is formulated as subject to The objective function and heat equation follow the previous section.Eq. ( 15) is a volume constraint that is introduced to regularize the varying volume fraction in the series of microstructures.  and | | represent the volume of an element and of the domain, respectively. is the number of elements in the domain. [] is the prescribed volume fraction of the th sample microstructure.
Recall the height field visualization of differentiable microstructures (Fig. 2(d)).A larger height  corresponds to a microstructure with a small volume fraction.We prescribe a linear relationship between the height parameter and the volume fraction, Eq. ( 16) sets design variables to be within [0, 1].The upper bound, together with  in Eq. ( 14), influences the frequency of occurrences of different temperature values in the resulting temperature field.For instance, a large uniform  and  will result in a temperature distribution that is dominated by temperature values close to 1. On the contrary, a small uniform  and  result in a temperature distribution that is dominated by temperature values close to 0. In both cases, the iso-contours of the temperature field are clustered, negatively affecting the optimization convergence.In our examples, setting the upper bound of   to 1 and  = 10 −3 is found to be able to obtain a temperature field that smoothly varies from 1 to 0, resulting in good convergence.
In this work, we optimize differentiable microstructures to maximize their bulk modulus.The bulk modulus is computed based on the homogenized elasticity tensor that is calculated using numerical homogenization [7,33].The optimization problem is solved using the method of moving asymptotes (MMA) [37].The optimized bulk modulus of the differentiable microstructures is compared to the Hashin-Shtrikman bounds [31].

Sensitivity analysis
The sensitivity of the objective function  () = ∑  =1 ( [] ) regarding the design variable is the number of elements.Take one of these elements for example, by using the chain law,    is defined as where   is the neighborhood of element .  []    is stated below: ( [] ) Here, ( [] ) = ∑ 3 ,=1    []   is a function of the homogenized elasticity tensors .  is constant, 1 or 0.
[]     []    is calculated based on the equilibrium equation and the modified solid isotropic material with penalization (SIMP) method (see e.g., [7]).
[]      is calculated from the Heaviside function.

X. Zhai et al.
is computed using the adjoint method.Let  0 =  − ( α).The Lagrangian function for the volume constraint is Then the derivative of  []   regarding α is calculated as: Let ]  , and then with  =   1 we have With  [] being solved from the following equation   0  [] +    [] = 0, []     α is rewritten as

Solution procedure
The pseudocode in Algorithm 1 details the optimization process of generating differentiable microstructures.This algorithm takes the filter radius , number of sample microstructures , and microstructure  0 as input.The output is the differentiable microstructures encoded by an optimized temperature field .5), (10), and (11).

Results and discussion
In this section, we thoroughly evaluate the proposed method by comparing optimized results with those obtained from alternative methods (Section 4.1), and by analyzing the geometric features and differentiability of optimized results (Section 4.2).This analysis is followed by validation of differentiable microstructures for two-scale topology optimization (Section 4.3).The proposed method is implemented using MATLAB, and tested on a standard desktop PC equipped with an Intel i7 CPU running at 2.90 GHz and with 64 GB RAM.The material model has Young's modulus of  = 1.0 and Poisson's ratio of  = 0.3.In all examples, the design domain is discretized into 200 × 200 bi-linear quad elements.

Results and comparison
Fig. 3(a) shows an input microstructure  0 , obtained by inverse homogenization for maximizing its bulk modulus at a volume fraction of 20%.The input microstructure  0 serves as the heat source in the optimization.Fig. 3(b) shows the optimized thermal diffusivity field .Fig. 3(c) and (d) visualize the corresponding temperature field based on iso-contours and a height field , respectively.In this optimization, the input microstructure, with a volume fraction of 20%, is morphed to reach a volume fraction of 90%.A total of 8 key microstructures is uniformly sampled from 20% to 90% with an interval of 10%.To evaluate the optimality of the obtained microstructures, we further sample the range to get 50 microstructures and plot their corresponding volume fraction and effective bulk modulus in Fig. 3(e).The theoretical limit from the Hashin-Shtrikman bound is plotted as a reference.
We further compare the optimized microstructures with those from two alternative methods.The results are visually compared in Fig. 4, using the height field as well as iso-contours.The slope of the height field is inversely proportional to the space between adjacent iso-contours.In Fig. 4(a), the input microstructure is uniformly dilated with an increasing dilation radius.This results in a constant slope in the height field and constant spacing between adjacent iso-contours.In Fig. 4(b), the input microstructure is morphed in a manner similar to the proposed approach, but with a uniform thermal diffusivity, i.e., no optimization of the thermal diffusivities is involved.Variations in the slope in the height field and the spacing between adjacent iso-contours can be observed.Fig. 4(c) corresponds to our approach.Here one can observe even larger variations in the slope in the height field and the spacing between adjacent iso-contours.This large variation is achieved with non-uniform thermal diffusivities.Three instances of the microstructures in each series are shown at the bottom.
To compare the mechanical properties of the three series of microstructures, we sampled 50 microstructures from each of them and plotted the bulk modulus of these microstructures under varying volume fractions (Fig. 5) .The theoretical limit from the Hashin-Shtrikman bound is plotted as a reference (the red curve in Fig. 5(a)).The bulk modulus of microstructures in the three series is shown as dots in blue, orange, and yellow, respectively.As can be seen from this plot, the results from our method (yellow dots) are closest to the H-S bound.This is followed by the orange dots corresponding to the uniform thermal diffusivity.The plot in Fig. 5(b) shows the relationship between the volume fraction  and the parameter .A linear relation between these two is prescribed in our method (Eq.( 17)).The results from our method (yellow dots in Fig. 5(b)) closely align with this linear relation.
To quantitatively measure the deviation from the theoretical limit, we define two metrics to evaluate the quality of the optimized microstructures. 1 measures the average of the difference between the bulk modulus of the optimized microstructure and the theoretical limit at different volume fractions. 2 measures the maximum of the differences.
Here Ṽ = { Ṽ1 , Ṽ2 , … , Ṽ } is the evenly sampled volume fractions in the predetermined range between   and   . is the number of samples. evaluates the bulk modulus of the optimized microstructure, while  is the Hashin-Shtrikman bound of the bulk modulus.In Fig. 6, we evaluate the quality of the optimized results by our method and those obtained by alternative methods.This evaluation is performed on ten different input microstructures (including the one in the above analysis).These ten input microstructures were obtained by inverse homogenization for maximizing their bulk modulus at a volume fraction of 20%.The filter radius and the power-law material interpolation in the inverse homogenization were modified to obtain diverse results.The three methods were then used to generate microstructures for the range of 20%∼90%.The evaluation of the optimized microstructures using  1 and  2 are summarized in Fig. 6.In each case, the bar graph, from left to right, corresponds to uniform dilation, uniform thermal diffusivity, and optimized thermal diffusivity.In all ten tested cases, the optimized thermal diffusivity results in  1 and  2 that are several times smaller compared to their counterparts from alternative approaches.

Differentiability and geometric features
The quality of the series of optimized microstructures depends not only on the optimization but also on the input microstructure.To understand the influence of the input microstructure, we look at results from the input shown in Fig. 7(a).The bulk moduli of 50 sample microstructures in this series are plotted in Fig. 7(d).It can be observed that there is a strong discontinuity in bulk modulus at a volume fraction close to 0.6.
The differentiability can be analyzed in relation to the geometric features shown in the height field visualization (Fig. 7(b)).The height field corresponds to the temperature field , representing an open surface.Three types of local surfaces may appear, as revealed by the cut-planes in Fig. 7(e-g), i.e., convex surfaces, concave surfaces, and saddle surfaces.These can be detected by analyzing the Hessian matrix ) at each point in the temperature field.Convex, concave, and saddle surfaces correspond to the determinant of the Hessian matrix being positive, negative, or zero, respectively.These three cases are analyzed in isolation in Fig. 8.For convex surfaces, the solid region of the parametric microstructure concentrates at the center.As the volume fraction increases, before the central solid region is in touch with the domain boundary, the bulk modulus of these microstructure instances is zero.In contrast, for concave surfaces, the solid region is initiated from the domain boundary, and the bulk modulus increases smoothly along with an increase in the volume fraction.At a low volume fraction, the solid regions corresponding to a saddle surface are disconnected, offering no stiffness in one of the two directions.As the volume fraction increases, the solid regions begin to merge, increasing the stiffness greatly.The topological change of the solid region at the saddle points results in discontinuity in the bulk modulus.
To ensure the differentiability of generated microstructures, we need to avoid local maximum points and saddle points in the temperature field.In the thermal diffusion approach, the maximum temperature only appears at the region defined by the input microstructure.It guarantees the temperature field is free of further maxima.The saddle points, however, may appear.Their appearance seems solely determined by the microstructure.To analyze the appearance of saddle points, we have performed a test on a simpler geometry.The input microstructure is shown in Fig. 9(a), while the microstructures by uniform dilation, uniform thermal diffusivity, and optimized thermal diffusivity are shown in Fig. 9(b)-(d), respectively.The saddle points are revealed by the cut planes in the subfigures at the bottom.The bulk modulus of the three series of microstructures is plotted in Fig. 10.The discontinuity can be observed in all three series, corresponding to the saddle points.However, they appear at different volume fractions.The discontinuity in our approach, with optimized thermal diffusivity, appears at the smallest volume fraction.This leaves flexibility for reaching optimal bulk modulus at higher volume fractions.The optimization of the thermal diffusivity achieves bulk modulus much closer to the theoretical limit.Our experiments have shown that saddle points are likely to appear if there are elongated holes in the input microstructure.This observation allows us to make a selection of input microstructures for design optimization.To confirm this observation, we have conducted further tests on input microstructures that have the same topology but different geometric features.The two input microstructures are shown in Fig. 11(a1)(b1).Case A has an elongated void region, in contrast to Case B. The correspondingly optimized height fields, section views, and iso-contours are displayed in Fig. 11 (a2)-(a4) and (b2)-(b4).Saddle points can be detected in the middle of the domain boundaries, highlighted in (a4).As a comparison, there is a local saddle patch in (a3), while at the same location, the patch in (b3) is concave.

Two-scale topology optimization
The concept of differentiable microstructures is introduced primarily for the design of two-scale structures with gradation.In this section, we validate their benefits using different two-scale topology optimization strategies.

Decoupled two-scale optimization
We follow the decoupled two-scale optimization approach as described in Ref. [7,19].In this approach, the series of microstructures is predefined, possibly by optimization.The relation between the elasticity tensor of microstructures, evaluated by numerical homogenization, and their volume fraction is established by polynomial curve fitting.This relation is then used to replace the power-law interpolation (i.e., SIMP) in topology optimization on the macroscale.Following optimization, each intermediate density is mapped to a microstructure in the predefined series which shares the same or the closest volume fraction.
We compare our method against ''uniform dilation'' and ''uniform thermal diffusivity'' (Fig. 12) as well as the compatibilityoptimized microstrucutres [7] (Fig. 13).The design domain and boundary conditions are illustrated in Fig. 12(a).The volume constraint is  * = 0.5.Under the distributed loads on the top and bottom of the domain, the density-based topology optimization with the power-law interpolation results in a large number of gray elements, as can also be seen in the histogram in Fig. 12(b).In this example, we take the microstructures from Fig. 4. Fig. 12(c) plots the relation between the elasticity tensor and the volume fraction for three families of microstructures, which are shown in sub-figures (d) -''uniform dilation'', (e) -''uniform thermal diffusivity'', and (f) -our approach.The resolution for the macroscale topology optimization is 60 × 30, with a relatively large filter size of 5 in order to promote a continuous gradation of density.The compliance of the optimized structure is evaluated using the homogenized elasticity tensor, as well as by discretizing the full scale using a resolution of 6000 × 3000.Here, each unit cell has 100 × 100 elements.In all three cases, the compliance by the full-scale analysis is larger than the compliance evaluated using the homogenized elasticity tensor, in alignment with previous results [7].The compliance from our approach is the smallest among the three, as a result of the optimized bulk modulus of microstructures by our approach.
The improvement is more pronounced when comparing our approach with the compatibility-optimized microstructures [7].Two series of microstructures are obtained by the compatibility optimization approach, using 5 and 10 key microstructures, shown respectively in Fig. 13(a) and (b).The microstructure with the smallest volume fraction in Fig. 13(b) is used to generate differentiable microstructures, as shown in Fig. 13(c).Again, macroscale topology optimization is performed based on the elasticity interpolation function from the three series of microstructures.In Fig. 13(a) and (b), the small set of microstructures results in volume deviation when mapping the intermediate densities to microstructures that have the closest volume fraction.From the compliance values, it is confirmed that the differentiable microstructures achieve the best results.

Concurrent two-scale topology optimization
In addition to the optimization of bulk modulus, the concept of differentiable microstructures can also be used in concurrent two-scale optimization where the microstructure and the macroscale layout are optimized simultaneously.The concurrent twoscale topology optimization employs two sets of design variables: the density distribution of the macroscale layout and the thermal diffusivity field of a single microstructure.Fig. 14 shows the test on a short cantilever beam, with a resolution of 30 × 30 and 100 × 100 on the macro and microscale respectively.The input is an X-shaped microstructure.While the input is rotationally symmetric, the obtained differentiable microstructures are anisotropic, due to the specific mechanical boundary conditions on the design domain.

Conclusions
In this paper, we have introduced the concept of differentiable microstructures -parameterized microstructures that are continuous both in terms of geometry and physical properties.To construct differentiable microstructures, we propose a novel  topology optimization formulation aimed at maximizing the bulk modulus of the entire series.The optimization morphs an input microstructure in a non-uniform manner to obtain a series of microstructures with an increasing volume fraction.When compared to uniform morphing approaches, our proposed optimization method achieves a series of microstructures with bulk moduli closely approaching the theoretical limit.
The quality of optimized differentiable microstructures depends on the input microstructure.We have analyzed critical geometric features of input microstructures and their impact on the differentiability.In our future work, we intend to explore the possibility of considering the input microstructure as a design variable.Furthermore, our current work is restricted to 2D differentiable microstructures optimized for maximum bulk modulus.While the concept of differentiable microstructures is potentially applicable to other physical properties such as Poisson's ratio, it may require different parameterizations of the optimization problem to effectively explore the solution space for these properties.

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

Data availability
Data will be made available on request.

Fig. 1 .
Fig. 1.Illustration of uniform offsetting (a) and non-uniform morphing (b) from an X-shaped 2D microstructure.They correspond to the cross-sections of different height fields (c and d).

Fig. 3 .
Fig. 3. (a) The input microstructure for topology optimization of differentiable microstructures.(b) The optimized thermal diffusivity .(c) The iso-contours of the optimized temperature field, interpreted as microstructure instances.(d) Visualization of the temperature field as a height field.(e) The bulk modulus of 50 uniformly sampled microstructures in the optimized differentiable microstructure.

Fig. 5 .
Fig. 5.The theoretical bulk modulus is plotted by the red line in (a) and real bulk modulus of 50 microstructures generated by #1, #2, #3 are represented by blue, red, and yellow dots.(b) shows the relationship between volume fraction and parameter .(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 6 .Fig. 7 .
Fig. 6.The bar graph of the evaluation factor  1 (blue bar) and  2 (pink bar) of three approaches under ten cases.The histograms from left to right show uniform dilation operation (1), uniform thermal diffusivity operation (2), and optimized thermal diffusivity operation (3), respectively.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

X
. Zhai et al.

Fig. 8 .
Fig. 8. Illustration of the influence of three types of surfaces on microstructure generation.

Fig. 9 .
Fig. 9.A simple input geometry and the microstructures generated by uniform dilation, uniform thermal diffusion, and optimized non-uniform thermal diffusion.

Fig. 10 .
Fig. 10.Comparison of the bulk moduli of microstructures generated with uniform dilation, uniform thermal diffusion, and optimized non-uniform thermal diffusion.The microstructures are shown in Fig. 9.

Fig. 12 .
Fig. 12. Decoupled two-scale topology optimization with differentiable microstructures.(a) Design domain and boundary conditions.The standard density-based topology optimization results in a large number of gray elements.(b) Histogram showing the number of intermediate densities between 0 and 1. (c) The elasticity interpolation function for three series of differentiable microstructures, obtained by ''uniform dilation'' (d), ''uniform thermal diffusivity'' (e), and our approach (f).The compliance is evaluated based on homogenization as well as full-scale analysis.

Fig. 14 .
Fig. 14.Differentiable microstructures applied to concurrent two-scale topology optimization.(a) Illustration of the design domain, boundary conditions, and the optimized two-scale structure.(b, c, d) Visualization of the optimized differentiable microstructures, showing the anisotropic feature due to the specific boundary conditions.