Adaptive thermodynamic topology optimization

The benefit of adaptive meshing strategies for a recently introduced thermodynamic topology optimization is presented. Employing an elementwise gradient penalization, stability is obtained and checkerboarding prevented while very fine structures can be resolved sharply using adaptive meshing at material-void interfaces. The usage of coarse elements and thereby smaller design space does not restrict the obtainable structures if a proper adaptive remeshing is considered during the optimization. Qualitatively equal structures and quantitatively the same stiffness as for uniform meshing are obtained with less degrees of freedom, memory requirement and overall optimization runtime. In addition, the adaptivity can be used to zoom into coarse global structures to better resolve details of interesting spots such as truss nodes.


Introduction
Topology optimization has been introduced several decades ago and it has been established as a powerful tool during engineering design processes. Review papers are provided by Rozvany (2009), van Dijk et al. (2013, and Huang and Xie (2010). Different target functions can be defined, among which the optimization of the mechanical stiffness of a system is probably the most prominent one (Sigmund and Maute 2013). The goal of this optimization problem is to find the spatial description of the topology, i.e., a distinct void/full material distribution. This can be expressed in terms of the so-called material density χ = χ(x) = {χ min , 1}, where the spatial coordinate is given by x.
Intermediate configurations, i.e., χ ∈]χ min , 1[, are difficult to be interpreted (foam) and even more challenging to be manifactured. Consequently, these "gray" solutions are to be avoided. A simple yet powerful strategy is the Continuum Mechanics, Ruhr University Bochum, Bochum, Germany multiplication of the compliance energy with a non-linear function in χ, e.g., = χ 3 0 is a very famous approach (Sigmund and Maute 2013). The effect of this non-linear interpolation between void and full material configurations turns the energy being non-convex, which, of course, renders the problem inherently ill-posed. A numerical artifact of the ill-posedness is the phenomenon of patterns of repeated black/white distributions that represent in average the gray solution. Due to its appearance, it is referred to as the checkerboard phenomenon (Diaz and Sigmund 1995). A prominent approach to prevent checkerboarding are filter schemes of which a huge variety can be found in literature (e.g., Sigmund and Petersson 1998;Bourdin 2001;Zhou et al. 2001;Lazarov and Sigmund 2011;Wadbro and Hägg 2015).
Along with the "pure" optimization of the mechanical stiffness, the inclusion of the precise material behavior is an emerging branch in the field of topology optimization (cf. Zhou and Wang 2007;Klarbring et al. 2017;Petrovic et al. 2018). It is obvious that accounting for the correct (non-linear) material behavior allows for using the potentials of design space and material at its best. In this regard, the so-called thermodynamic topology optimization has been introduced (Junker and Hackl 2015;Jantos et al. 2018). It makes use of mathematical schemes known from the field of material modeling to derive a suitable equation for the evolution of the material density. The inclusion of a specific material model is straightforward since the very same procedures can be employed to derive the corresponding evolution equations for the material description. Examples have been presented for anisotropic materials in Jantos et al. (2017) and systems of tension-and compression-affine materials, e.g., steel/concrete structures, in Gaganelis et al. (2019).
One advantage of the thermodynamic topology optimization is that the regularization is directly included in the governing equation, meaning that additional filter techniques are obsolete. Consequently, the thermodynamic topology optimization is immediately eligible for advanced numerical treatment. To be more precise, adaptive strategies for the spatial resolution can be employed which are of particular importance when large-scale engineering structures are supposed to be computed within one optimization run: details of the resultant topology can be resolved directly which allows to optimize, e.g., a bridge in the dimensions of meters but simultaneously accounting for details at the supports or connecting adapters. Since the adapters are usually problematic for designing, an iteration of large-scale and small-scale simulations whose respective results are used for the other computation, demanding a precise inclusion of boundary conditions, can be avoided: it is replaced by one single holistic optimization run that makes use of adaptive meshing.
Adaptive techniques have been used for a variety of topology optimization strategies to improve the quality and sharpness of the solution as well as to address efficiency issues for finely resolved meshes. A popular approach for optimization based on finite element (FEM) approaches consists in adaptive mesh refinement (AMR) that is employed in an h-refinement fashion to increase the accuracy where needed while saving computational cost in other regions. One of the earliest works in this direction has been presented by Maute and Ramm (1995). Starting with an uniform coarse mesh, they perform a density-based optimization which is then used to create a finer mesh with refinement at indicated spots. The new mesh is employed in a subsequent optimization and the process iterated until the desired degree of accuracy is reached. Derose (1996) employs octree data structures to implement adaptive meshing strategies for optimization. Lin and Chou (1999) use a two-stage approach to reduce computational costs. Arantes and Alves (2003) use the h-refinement topology optimization approach employing mesh quality and an analysis error estimator as criteria for refinement. Stainko (2006) employs an adaptive multilevel approach refining towards the material-void interface and a multigrid method to efficiently solve the elasticity problem. De Sturler et al. (2008) (see also Wang et al. 2010) use a dynamic adaptive mesh strategy with mesh refinements and coarsening in every optimization step to obtain the same designs as for a uniformly refined mesh. Guest and Genut (2010) separate the design variable field from the analysis mesh and thereby reduced the computational cost. Nguyen et al. (2010Nguyen et al. ( , 2012 employ multiresolution topology optimization. Bruggi and Verani (2011) present a fully adaptive algorithm with adaptive coarsening and refinement based on two goal-oriented error estimators. Wallin et al. (2012) employ the phase-field method and adaptive finite element formulations. Wang et al. (2013Wang et al. ( , 2014 separate the density field from the FEMcomputed displacements. Using a point-based background representation of the density and Shepard interpolants, they can obtain a high resolution of the solid-void interface. Nana et al. (2016) employ h-adaptivity for tetrahedral meshes. Nguyen-Xuan (2017) presents topology optimization on 2d polygonal adaptive meshes, which is generalized to multi-material by Chau-Nguyen et al. (2017). Panesar et al. (2017) use hierarchical remeshing strategies and study different mesh adaption techniques. Lambe and Czekanski (2018) employ a continuous density field and adaptive mesh refinement to sharply resolve the solidvoid interface. Salazar de Troya and Tortorelli (2018) use AMR in stress-constrained topology optimization. For the use of filter techniques during AMR, we exemplary refer to Nguyen-Xuan (2017) for polytree approaches and to Lambe and Czekanski (2018) which use a Helmholtz equation. In addition, the employment of the multigrid method (Hackbusch 1985) can considerably reduce the computational cost to solve the finite element analysis. This is in particular of interest since adaptive mesh refinement can be directly employed to create a nested mesh hierarchy that is used for the multigrid algorithm (Bramble et al. 1991;Bastian and Wittum 1994). In this direction, Amir et al. (2014) employ a multigrid-CG for efficient topology optimization. Aage et al. (2015) present an open-source framework based on PETSc employing multigrid. Kennedy (2015) uses geometric multigrid to solve large-scale multimaterial topology optimization. Chin and Kennedy (2018) use a parallel geometric multigrid for octree-based AMR optimization.
The paper is structured as follows: we begin with a short recall of the thermodynamic topology optimization for convenience. Afterwards, we present the employed numerical treatment, including a discretization of the Laplace term that is particularly useful when dealing with mesh adaptivity where hanging nodes might appear. As main contribution, we analyze in detail the numerical behavior of the thermodynamic topology optimization approach when used in conjunction with adaptive meshing techniques. We highlight that a locally chosen regularization parameter is necessary to avoid checkerboarding and show that qualitatively the same structures are obtained for adaptive and full refinement strategies. We show the difference between mesh adaption in every timestep from the beginning of the optimization and finer mesh adaption only at later stages. Four different test cases are considered and results in terms of obtained stiffness, stability and runtime presented.

Thermodynamic topology optimization
The fundamentals of thermodynamic topology optimization are recalled here for convenience. For more details, we refer to Junker and Hackl (2015) for the general idea of thermodynamic topology optimization, to Jantos et al. (2018) for an advanced numerical treatment, and to Jantos et al. (2019) for a comparison of the method with SIMP.
For the derivation of the thermodynamic topology optimization, we make use of Hamilton's principle. It is formulated for quasi-static, elastic continua as where G denotes the Gibbs energy and δu the virtual displacements. Consequently, the expression δG refers to a Gâteaux derivative. The Gibbs energy is given by where the Helmholtz free energy is given by , the body forces by b, and the tractions by t. The body's volume and surface are denoted by and ∂ = u ∪ σ , respectively, with the Dirichlet boundary u and the Neumann boundary σ . For inelastic continua, so-called internal variables υ are introduced that describe the current microstructure. Examples are plastic strains and damage variables. Then, Hamilton's principle is expanded to read with the dissipation functional δD := p δυ dV (4) and the non-conservative forceŝ The choice of the dissipation potential diss = diss (υ) determines the character of the Euler and Helmholtz equation, respectively, which results from (3). A great benefit of Hamilton's principle is its character of a mathematical potential, allowing to account for constraints by adding further potentials, e.g., δR and δC, viz For the case of topology optimization, the internal variable is the local material density χ = χ(x) for which the dissipation potential is chosen as with the viscosity η. As seen later, this approach yields a transient term in the final equation for χ.
The Helmholtz free energy is formulated in the mechanical stresses σ and reads where E denotes the elasticity tensor of order four. The exponent 3 is chosen in accordance to standard approaches from literature. This exponent, of course, yields on the one hand distinct void/full material distributions; on the other hand, it renders the problem ill-posed due to the loss of convexity. It is thus mandatory to introduce a functional for regularization, defined by Finally, the constraints are accounted for by the constraint functional with the Lagrange parameter λ ensuring that the mass of the structure equals the given mass , and the Kuhn-Tucker parameter γ for the interval constraint χ ∈ [χ min , 1]: whereγ ≡γ (x) is chosen such that an overhead of the driving force which would cause χ ∈ [χ min , 1] is corrected, i.e., the parameter is locally chosen such that χ remains in the admissible interval. In principal, the setvalued character of γ complicates the evaluation of the equations since a coupled system of inequalities has to be solved. However, instead of computingγ numerically, this constraint can equivalently be implemented by a minmax condition applied to a tentative density. The mass conservation is then ensured by using a bisection approach (cf. Sigmund 2001), as indicated in Algorithm 1 to find an appropriate offset for the driving force. More details can be found in Jantos et al. (2018).
The first stationarity condition in (12), δG [ε, χ], gives the balance of linear momentum in its weak form; the second condition provides the Helmholtz equation for the material density: Here, the constitutive relation for the stress with the effective stiffness χ 3 E has been used and the driving force has been introduced. The Laplace operator is defined by

Numerical treatment
For the numerical treatment of (13), we discretize the displacements u by linear finite element functions and the density field χ employing an elementwise constant approximation. To this end, the domain is meshed with quadrilaterals or hexahedrons and we denote the collection of elements by h . For two-dimensional approximations, the plain stress assumption is used. Starting from an initial mesh, we allow the refinement of individual elements and thereby also hanging nodes in the mesh. A 2:1 balance is, however, enforced such that only one hanging node per side can appear. Displacement degrees of freedom located on a hanging node are handled by continuity constraints in order to guarantee an overall continuous solution. The maximal number of refinements required to create a mesh from an initially coarse mesh is denoted as the maximal level of refinement L. For the computation of the Laplacian for the piecewise constant density field, we use a fitted stencil method based on the set of the nearest neighbors of an element e and the element itself, denoted by N e , illustrated in Fig. 1. Using the barycenters xẽ of elements connected via a mesh node to the element e, the stencil coefficients are chosen based on a second-order Taylor expansion, i.e., in the way that = ẽ∈Ne a eẽ χ(x e ) + ∇χ | xe · (x e − xẽ) Adaptive mesh hierarchy with hanging nodes: Shown are the neighbors (orange) of an element e and their barycenters (red cross). Neumann zero boundary condition are handled via a mirror element at boundary sides (green) is fulfilled with a first order error term, where H is the Hessian matrix. We used the same treatment for a gradientenhanced damage model in Vogel and Junker (2020) where we present a more detailed description. An analogous treatment has been used in Jantos et al. (2018) employing a fixed number of neighbors for topology optimization and similar discretization methods are discussed, e.g., in Coirier (1994) and Sadat and Prax (1996). Please note that the matrix A solely depends on the meshing and has to be recomputed only after a mesh adaption. The Neumann boundary conditions are incorporated by mirror elements at the boundary that are forced to posses the same value as the element inside the domain. For every optimization step, at first the displacement u for the current density field χ is computed. We use a geometric multigrid solver (Hackbusch 1985;Reiter et al. 2013) on the mesh hierarchy produced by successive refinement starting from a coarse uniform initial mesh. In case of adaptivity, an adaptive multigrid approach (Bramble et al. 1991;Bastian 1996) is employed without smoothing of the patch rims. All simulations have been carried out with a symmetric Gauss-Seidel smoother, V-cycle, 3 pre-and postsmoothing steps, canonical prolongation and restriction, RAP product for coarse matrices, and LU base solver on the coarsest mesh level. The multigrid is used in a CG solver as preconditioner. Residua have been reduced to an absolute value of 5 · 10 −8 .
Subsequently, given the displacements, the density is updated as presented in Algorithm 1. As suggested in Jantos et al. (2018), we solve a number of inner timesteps for stability reasons and also use the problem-independent control parameters η * and β * . In most cases, we employ the suggested value of β * = 2h 2 min (Jantos et al. 2018) which corresponds to an internal length scale (Peerlings et al. 1996;De Borst and Mühlhaus 1992) and thereby prescribes the minimal size of representable structures. This choice and the typically chosen value of η * = 15 then imply a single inner timestep. In order to compute the Lagrange parameter λ and fulfill the average density requirementρ for χ, a bisection algorithm (Sigmund 2001) is used. Density constraints are enforced explicitly after computing a trial density. Most notably, for the adaptive handling, we employ an elementwise choice for the stabilization parameter β → β e in the update of the density. This parameter is chosen at least large enough to ensure a local stability criterion (β e = p avg 2h 2 e ) or to satisfy a global requirement (β e = p avg β * ). On coarse elements, the local stability criterion will dominate and ensure stability. For regions with fine elements, the global criterion will control the minimal length scale.
The mesh adaption is controlled based on the updated density. The employed heuristic treatment is shown in Algorithm 2. We only allow a refinement up to a prescribed level of detail L. Elements are selected for refinement, if the density value χ e is still in the gray-zone between some bounds away from 0 or 1. In addition, elements are selected for refinement if the density difference of neighbored elements is large, i.e., if the gradient in the density field is large. On the contrary, elements are selected for coarsening, if the difference in density to the neighbors is small and the density values are sufficiently close to the black-white solution, i.e., close to 0 or 1. This approach is motivated by the observation that high resolutions are required to capture in detail the material-void interfaces where steep density gradients are present. For parts with only material or only void, the density is constant and no high resolution is required there.

Numerical results
We provide a detailed study for the behavior of the thermodynamic topology optimization method in conjunction with adaptive meshing employing the symmetrically reduced Messerschmidt-Bölkow-Blohm (MBB) beam as a running example. Subsequently, we present two additional 2d settings and one 3d example to show that the obtained results carry over to these boundary value problems. If not otherwise stated, the simulation parameters presented in Table 1 are used. The adaptive topology optimization has been implemented in the UG4 simulation software , meshing is performed with ProMesh (2020) and visualization employs Paraview (2020).
The main questions to be addressed are as follows: In Section 4.1.1, we compare different methods to adapt the mesh for the nonuniform case. A straightforward approach is to account for the slightly changed mass density after every optimization step and adjust the mesh accordingly. This approach will be compared with respect to stiffness and runtime to the full refinement approach. In particular, it will be of interest if the adaptive meshing restricts the design space too much in order to allow for finding a suitable structure compared with the uniform meshing. As a second approach, we let the density evolve for a certain number of optimization steps on a coarse mesh refinement, and then start with adapting to finer resolutions. This way, a coarse structure is first optimized which can subsequently be refined by the finer resolution. Since the optimization path has already evolved towards a certain optimum before allowing finer meshing, the final result may differ due to the different initial guess at fine level. In Section 4.1.2, we study in detail the influence of the stabilization parameter β. In particular, for the adaptive strategy, the difference between a locally or globally chosen parameter will be highlighted. Finally, in Section 4.1.3 we turn our attention to the mesh independence of the method by fixing the internal length parameter β and refine fully or adaptively far smaller than this size.

MBB beam
The domain specification and boundary conditions for the symmetrically reduced Messerschmidt-Bölkow-Blohm (MBB) beam are shown in Fig. 2a. The coarsest meshing for the domain, the level L = 0, is shown in Fig. 2b (top). The mesh is either refined uniformly (full refinement) or employing an adaptive strategy. An exemplary adaptive mesh occurring during the optimization is presented in Fig. 2b (bottom).

Adaption strategies and runtimes
We compare different meshing strategies with respect to their obtained optimal structure, the final stiffness of the structure, and the required runtime for the result. Motivated by the empirical finding for fully refined meshes in Jantos et al. (2018), we will use the choice β = 2h 2 for the uniform meshing and a locally chosen value of β e = 2h 2 e on every mesh element e in order to avoid numerical instabilities like checkerboarding. In particular, this implies a shrinking regularization parameter with mesh refinement such that the regularization is chosen just large enough to ensure stability. By this, the imposed structural length scale decreases with finer meshes and finer structures can be resolved. Fig. 3, an example for a typical optimization run is presented. The full and adaptive meshing strategy is compared for a finest mesh level L = 4, i.e., we allow the adaptive elements to be as fine as the full refinement, however, also to use coarser elements if this seems appropriate due to the adaption strategy. Several optimization steps are shown and in all cases the same density distribution is observed for full and adaptive meshing. This already provides us with the confidence that the meshing can be coarsened appropriately at certain locations while maintaining the design space broad enough to reach the same optimized structure. Taking a look at the element sizes for the different optimization steps, we find a quite large number of small elements at early stages where the structure cannot be clearly anticipated and thus fine resolution is required to allow for many possible further optimization directions. However, at later stages, only the emerged beams have to be resolved and this is achieved by quite a few fine elements at boundaries between void and material. In both meshing strategies, several smaller structures appear in early stages, but these structures are later removed resulting in the same computed optimum with only thick beam structures at the end. This is due to the relatively large regularization parameter β that imposes a bound on the minimally resolvable structure as an internal length. In conclusion, mesh adaption in every optimization step produces qualitatively the same topology as uniform meshing, i.e., the coarsening only restricts the discrete design space where nothing interesting takes place. In Fig. 4, we compare the computed structures between adaptive and uniform meshing for different levels of maximal refinement. The density distribution after 500 optimization steps is presented together with the employed adaptive mesh to resolve the structure. On all refinement levels, a very similar final structure is computed in both cases. Even for the very detailed resolutions, there is hardly seen a difference in the resulting structures. This again shows that mesh adaption in every optimization step can be used to obtain an optimization result very close to the fully refined meshing.

Mesh adaption in every optimization step In
We investigate these findings systematically with respect to the obtained stiffness in Fig. 5a on different refinement levels. To the left, the evolution of the structural stiffness during the optimization is shown for adaptive and full meshing with a maximal level for refinement up to level 7. To better visualize the difference at the end of the optimization, a zoom has been added. In all cases, a quite steep increase in stiffness is obtained until approximately 100 steps. Subsequently, the stiffness for further mesh for the optimized structure at optimization step 500. On every mesh level, the regularization parameter is chosen according to the mesh size allowing more detailed structures on finer meshes refinement is higher for finer meshes since the structures can be resolved in more detail and the design space is larger. Comparing the meshing strategies, the results for adaptive and full refinement are very similar with a slightly better stiffness for the adaptive case. This can also be seen in Fig. 5a (right) which shows the finally obtained stiffness for the different mesh levels. Employing a finer resolution provides a better stiffness in all cases. Correspondingly, the amount of the black-white parts in the obtained solutions increase. To see this, the measure of non-discreteness (MOD; Sigmund 2007) defined as is presented in Fig. 5b and decreases with every mesh refinement. For adaptive refinement and full meshing, a very similar measure of non-discreteness is obtained especially for the final structure. This is due to the mostly constant density being full material or void such that only the transition zones have to be resolved with finer resolution. In Fig. 5c, we take a look at the evolution of degrees of freedom and the overall runtime for 500 optimization steps. All measurements have been performed on the computing cluster at the Department of Civil and Environmental Engineering of the Ruhr University Bochum employing one core of a 2.4 GHz Intel Xeon Skylake Gold 6148 processor. The required degrees of freedom are less for the adaptive case for all steps. At the beginning of the optimization, they are slightly higher to account for the not yet evolved structure and then decrease towards the end of the simulation. For the runtime, we see that the density update requires substantially less time than the solution of the elasticity problem. Comparing the meshing strategies, the adaptive case is faster with respect to the overall simulation time. This gap increases with allowed mesh resolution and is close to an order of magnitude at level 7. The nonneglectable overhead for mesh adaption is overcompensated here by the faster elasticity solver due to the reduced number of degrees of freedom. Since we employ a multigrid solver with linear complexity w.r.t. the unknowns, a factor of ten in savings for the degrees of freedom directly translates into a saving of the same factor in runtime with a small reduction due to the mesh handling. For other solvers with suboptimal complexity like conjugate gradient or Gaussian elimination, this gap might even be larger. For the uniform case, the simulation runtime is mainly given by the elasticity solver and the plotted curves therefore coincide. In Fig. 6a, the iteration counts for the multigrid method to solve the elasticity problem are shown. While they are pretty constant over a large range of steps, the required number  of iterations increase significantly at certain optimization steps. A detailed investigation revealed that this is due to particular density patterns occurring when thin beams are removed from the structure during the optimization. In these cases, diagonal density chains of elements can occur that are only connected via vertices. For such a difficult situation, the solution is not easily computed which is reflected in the iteration count. However, such intermediate structures  The average iteration count increases slightly under mesh refinement. We attribute this to the density field with different values on each element which poses a challenge to the multigrid solver. In Fig. 6b, the required number of bisection steps is shown for the different mesh levels and meshing strategies. In all cases, roughly 40 steps are required which still holds true even under mesh refinement.

Mesh adaption using initially coarse structure
In the previous section, we employed mesh adaption in every optimization step right from the beginning of the optimization. As an important result, it has been shown that the employment of mesh adaption in every optimization step keeps the design space broad enough to produce very similar structures as a full meshing. In contrast, we now ask the question if differ-ent structures can also be obtained on purpose by a different remeshing technique. This might be useful, e.g., if a coarse overall structure is already known or desired, but the details of this structure are still to be optimized. For example, a truss structure can be desired as overall structure but the optimal design of the truss nodes might still be unknown. Our strategy to address such a setting is to first employ a couple of optimization steps with a coarse mesh to find a coarse initial structure. This coarse structure is then used as initial guess for the subsequent highly resolved adaptive meshing strategy as before. This way, the optimization already starts closer to a certain coarse optimum and we want to investigate how this strategy changes the obtained results for fine meshing.
In Fig. 7 steps on a coarse mesh resolution to compute a coarse inital structure and only use fine, adaptively refined meshes for the remaining optimization steps. Comparing these results with Fig. 4, where mesh adaption in every timestep has been employed from the beginning, we see that the overall coarse structure remains the same for all mesh levels and the optimization algorithm is indeed not capable to leave the already approached local optimum even if fine meshing is added by adaptivity at the later stages. However, the structure is altered in the details. A zoom to such a finer substructuring of the material in a truss node has been added together with the adaptive mesh. Such a strategy can thus be used if the global overall structure should be close to a coarse structural resolution but the details in truss nodes have to be further known in detail.

Influence of locally varying regularization parameter
For all previous results, the regularization paramater has been chosen as β = 2h 2 for the uniform meshing and a local value of β e = 2h 2 e on every mesh element e for the adaptive case. For the uniform meshing case, it has been reported in Jantos et al. (2018) that smaller values lead to checkerboarding solutions. Here, we investigate the choice of the regularization parameter for the adaptive case. Intuitively, the regularization parameter penalizes steep gradients in the density function and larger values of β smear out the interface between void and material. For stability reasons, however, the regularization should at least enforce a transition of the density within the range of one finite element or more. On the other hand, it should not be much larger to allow for a sharp representation of the obtained structure. In particular, adaptive meshing is employed to better resolve the solution at the material boundaries and a large value of β, although clearly stable, would be counterproductive since it would limit the obtainable structures to a large internal length and render adaptivity useless. Therefore, we have used the elementwise chosen regularization parameter that enforces a gradient of the solution only for the locally seen mesh size. As an alternative, one could opt for a globally chosen regularization that is as small as suggested by the smallest element of the mesh. This allows to maintain very thin structures at the material-void boundaries, however, is too small for larger elements encountered on other parts of the domain.
In Fig. 8, we present the obtained density field at step 20 of the optimization for all three cases: full refinement as well as adaptive refinement with local stabilization and with a globally chosen value of β. While there is no checkerboarding observed for the full refinement and the adaptive meshing with local β e = 2h 2 e , a significant checkerboarding is observed for the globally χ (full mesh) χ (adaptive, localβ ) (adaptive, global β ) On every mesh level, the regularization parameter is chosen according to the mesh size allowing more detailed structures on finer meshes chosen small regularization β = 2h 2 min in the adaptive case. These patterns are mainly observed for regions where the density is still in the gray regime between material and void.
In contrast, the finally obtained structures after 500 optimization steps are pretty similar as shown in Fig. 9. While it is interesting that the same structures are computed although a checkerboarding instability is encountered during the optimization path, the absence of checkerboarding for the final structure is intuitively reasonable: the regularization only has to enforce its gradient penalization purposes where gradients are actually present. In the given density fields, the solution is constant (0 or 1) in a coarsely resolved region. These are the spots where the regularization-although too small in principle-is not required at all. For regions with a non-zero gradient in the density function, the mesh is adapted to a very fine resolution and the small regularization matches with the mesh resolution. In conclusion, for the final structure without almost no intermediate values for the solution, the regularization is either not required for coarsely resolved regions with a constant density field or the mesh is adequately resolved for the global regularization at spots with steep density gradients. The unphysical density with a checkerboarding pattern can also be observed in the iteration counts for the solvers in Fig. 10. While the bisection count is not affected, the multigrid count increases dramatically at early optimization steps since the highly jumping coefficients pose a severe challenge to the geometric multigrid. Towards the end of the optimization this effect is mitigated as the checkerboarding disappears.

Mesh independence
Finally, we take a look at the mesh independence of the obtainable structures for the adaptive meshing. We therefore computed an optimized density field for a fixed regularization parameter β = 9.74 · 10 −4 with different levels of mesh refinement for the full and adaptive case. The results at the end of the optimization are shown in Fig. 11 and all show qualitatively the same structure regardless of the mesh resolution. This highlights that the stabilization parameter enforces the minimally resolvable structure sizes and that the obtained solution is independent of the mesh size as long as the mesh size is equal or smaller than the stable mesh resolution β e = 2h 2 e . In Fig. 12, the corresponding stiffness and measure of non-discreteness is shown for the mesh refinement levels. They are quite constant as expected from the qualitatively same structure.

L-shaped cantilever
The next example is the 2d L-shaped cantilever (see Fig. 13). Here, all displacements at the upper part of the cantilever are fixed and a force at the top of the L-tip is pointing downwards. The resultant optimized structure is presented in Fig. 14 after 200 time steps. Again, it can be observed that the adaptive mesh yields almost identical results as the full mesh resolution. However, at refinement level L = 5, different structures are obtained for the full and the adaptive mesh. This is due to the inherent small numerical differences for adaptive and full treatment which can lead to different local minima in principle. For most examples, only very small deviations are observed, but in this case the evolutionary thermodynamic topology optimization finds a different local minimum for the adaptive case compared with full mesh refinement. As expected, the stiffness increases with increasing mesh refinement as seen in Fig. 15a. Even for mesh refinement with L = 5, the adaptive mesh yields a result with higher stiffness than the one obtained for full mesh refinement. The measure of non-discreteness shrinks for higher levels of mesh refinement (see Fig. 15b). The benefit of adaptive finite element strategies is presented in Fig. 15c: the number of unknowns is smaller for adaptive simulations and this gap enlarges with mesh refinement. The advantage of a reduced number of unknowns is shaded for coarse levels due to the additional numerical costs caused by adaptivity. However, for finer meshes, this drawback is overcompensated and the adaptive simulation is faster.

Bridge-like setting
The last 2d example is the bridge-like setting as shown in Fig. 16. Here, we make use of symmetry conditions at the left-hand side. Furthermore, the displacements are fixed at the support corner at the lower right. On a centered horizontal line, a constant force is acting in negative y direction which might be caused by a roadway. For this problem, the parameters in Table 1 are employed except for the density which is chosen asρ = 0.2.
The optimal structure is presented for various mesh refinements in Fig. 17 after 100 time steps: an arch connects the support points and makes use of the entire design height available. The roadway, thus the horizontal line with line loads, is joint to the arch by pillars where the roadway lies above the arch and by hanging structures where the roadway lies under the arch. The hanging structures are both vertical lines as well as rope-like structures close to the center. While the arch remains more or less unmodified during mesh refinement, the joints between arch and roadway can be resolved with much higher accuracy.
The impact of mesh refinement in terms of stiffness, measure of non-discreteness, degrees of freedom and runtime is collected in Fig. 18. The already observed characteristics are present here also: finer levels of mesh refinement yield stiffer structures by simultaneously reducing the measure of non-discreteness. Furthermore, the numerical advantage of adaptive strategies is more pronounced for larger levels of refinement, resulting in a runtime saving of more than one order of magnitude for level 7.

Three-dimensional cantilever
We present a three-dimensional setting for a wall-fixed cantilever. The domain specifications and boundary conditions as illustrated in Fig. 19 are used. Half of the domain is computed and extended by symmetry boundary condition. The domain is fixed at the rear while a downward-pulling force is applied at the bottom front center. The simulation parameters are chosen as indicated in Table 1. In Fig. 20, the results after 200 optimization steps are shown for full and adaptive refinement strategies. Very similar final structures can be observed while the computational mesh consists of less elements in the adaptive case. Employing a regularization parameter adapted to the mesh size allows to resolve finer structures on finer meshes such as a thin material plane. In Fig. 21, the corresponding statistics for the computation on the different mesh levels are shown. Better stiffness and measure of non-discreteness are obtained for finer meshes similar to the two-dimensional cases. The runtimes on fine meshes are better for adaptive refinement than for full refinement. This is in agreement with the savings in degrees of freedom which become χ (full mesh) χ (adaptive mesh) adaptive mesh The results for a fixed regularization parameter β = 0.08 are presented in Fig. 22. Again, the adaptive meshing is capable to produce qualitatively the same results as a fully refined mesh. The fixed regularization parameter imposes a minimal structural size as can be observed by the computed topologies on finer meshes. The refinement only sharpens the resolution in this case but does not allow to resolve more detailed structures.

Conclusions
We presented a sophisticated study for the adaptive treatment of the thermodynamic topology optimization. As main finding, the stabilization parameter has to be chosen in an elementwise fashion in order to allow small gradient penalizations locally while maintaining stability globally. This choice leads to qualitatively the same results as a fully refined mesh with the same minimal mesh size but the required degrees of freedom, memory consumption and runtime are significantly reduced. Thus, by adaptive mesh refinement and a locally chosen regularization, the design space can be limited during the thermodynamic topology optimization for computational savings while still maintaining the search space broad enough to allow for the same optimized structure as fully refined meshing. The savings due to adaptivity become more dominant with finer mesh resolution and we expect even larger savings for very fine resolutions. The employment of fine meshes, as a general rule of thumb, improves the stiffness of the obtained structures.
In contrast, if the regularization parameter is fixed globally to the smallest element size, checkerboarding is observed and the iteration count of the employed geometric multigrid solver is disadvantageously impacted. If the regularization is chosen globally as a large value, this limits the obtainable structure sizes to an internal length scale although it provides mesh independent results. However, the employment of fine meshes is not meaningful with a large regularization since fine structures can not be computed due to the imposed large internal length scale. From a computational point of view, it is thus advisable to choose for a given mesh the regularization large enough to maintain stability but as small as admissible to compute as many structural features as possible. This strategy requires to shrink the regularization parameter with mesh refinement and has been employed in this work. From a manufacturing point of view, it is advisable to choose the regularization parameter as fine as the minimal structural size which can still be manufactured. The mesh size can then be chosen accordingly and further refinement will not add more structural details but only sharpen the structural boundary representation. Adaptivity is used for computational savings in this case.  As an interesting mesh adaption variation, the interplay between coarse and fine resolutions can be used to first obtain a coarse initial global structure and then zoom into relevant spots, e.g., at truss nodes, to resolve details with higher mesh accuracy.
is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creativecommonshorg/licenses/by/4.0/.