A level set-based topology optimization approach for thermally radiating structures

The need for efficient thermally radiating structures is apparent in many aerospace system designs including satellites, launch vehicles, and hypersonic aircraft. This paper presents a novel level set-based topology optimization approach for designing thermally efficient radiating structures. In this paper, we derive a shape sensitivity of the thermal heat power radiated objective function using the adjoint method. This sensitivity is a necessary ingredient for our gradient-based algorithm. We apply an augmented Lagrangian method to solve an example 2D problem where the goal is to maximize heat power rejected subject to a material volume constraint. The radiating surface is kept fixed during the optimization to maintain a design-independent boundary condition, while the conducting region is optimized. Several solutions are illustrated with varying initial conditions. We also present a case study indicating that maximizing the thermal compliance functional is not sufficient for solving this class of problems.


Introduction
Thermal radiation is a physical phenomenon affecting the design of many aerospace systems including hypersonic vehicles and spacecraft. The need to design lightweight, efficient structural elements within these systems has driven the development of many design optimization tools. Topology optimization has grown to become a popular method for light-weighting structural elements due to its ability to modify an existing design in a free-form manner.
Topology optimization methods identify an optimal material layout within a fixed design domain subject to a set of boundary conditions and other constraints to minimize a cost functional. While other structural optimization methods handle a fixed set of design variables such as the lengths, thicknesses, and other relevant dimensions of predetermined features, topology optimization allows for variation in the overall shape and connectivity of the structure (2014). This procedure is typically more advantageous in earlier design stages where a variety of geometries can be explored. Several papers provide overviews of the various topology optimization methods and current areas of research (Eschenauer and Olhoff 2001;Rozvany 2009;Sigmund and Maute 2013;Deaton and Grandhi 2014).
Much research has been done on optimizing thermally conductive structures using topology optimization methods (Dbouk 2017), including heat exchanger design applications with convective boundary conditions (Coffin and Maute 2016;Kambampati and Kim 2020), and transient thermal problems (Hyun and Kim 2021). However, less attention has been given to solving problems incorporating a radiative boundary condition. Two recent works include Castro et al. (2015) and Munk et al. (2017). Castro et al. (2015) investigated the design of radiative enclosures by redistributing reflective material using a solid isotropic material with penalization (SIMP)-based approach. The 2D element density variables were used to interpolate surface reflectivity within the enclosure and subsequently determine the temperature and heat flux fields. Munk et al. (2017) used an evolutionary structural topology optimization approach on a Responsible Editor: Seonho Cho three-dimensional panel model of a hypersonic aircraft wing incorporating aerothermoelastic effects. The evolutionary optimization method was a hard-kill bi-directional method; element densities were integer variables and void elements were excluded from the finite element and radiation analysis altogether. The objective was to minimize weight and stress.
In the approach used by Castro et al. (2015), the radiating surface does not change topologically, and in Munk et al. (2017), the radiative heat transfer was not included in either the objective or constraint functions. A level set method approach which explicitly incorporates the radiative heat transfer in the optimization problem would be advantageous for two reasons: (1) a design-dependent load problem, such as those involving structures subjected to hydrostatic pressure loads (Neofytou et al. 2020), would be easier to formulate since the structural boundaries are implicitly defined by the level set function and (2) the designer could manipulate the radiative capabilities of the final design more easily.
The remainder of this paper is organized as follows: Sect. 2 provides some background on level set methods and reviews how they implicitly represent structural boundaries. Section 3 presents a two-dimensional single-objective, single-constraint optimization problem to design a thermally efficient radiating structure. In Sect. 4, a level setbased topology optimization approach is defined using an augmented Lagrangian method. We then derive the shape sensitivity in Sect. 5 for maximizing thermal heat power radiated. In Sect. 6, we demonstrate the proposed methodology on a 2D radiating aluminum plate problem with varied initial material configurations to illustrate how the method is identifying locally optimal solutions. We also present several candidate designs to illustrate that using a thermal compliance objective function is insufficient for this class of problems. Finally in Sect. 7, we provide a summary of the contributions presented in this paper and discuss future work.

Background
Topology optimization offers the designer flexibility in the geometric shape of the structure under consideration. Typically, a finite element model is employed to evaluate the objective and constraint functions, in addition to any sensitivities needed for implementing a gradient-based optimizer. A density value is assigned to each element of the mesh which interpolates the stiffness, conductivity, or other relevant material property throughout the design domain. Because the number of design variables in topology optimization problems scales directly with the number of nodes or elements in the finite element model, they tend to be much larger in dimension than parametric design problems. However, topology optimization offers much greater design freedom since a nearly limitless number of shapes can be obtained.
Boundary variation methods are a subset of topology optimization routines which employ an implicit representation of the structure's boundaries rather than an explicit parameterization of the design domain (Deaton and Grandhi 2014). Several methods fall under this category including level set method (LSMs) and phase-field methods. Both utilize an implicit function to define the structural boundary as well as an evolution equation used to update between optimization iterations. Compared to the level set approach, phase-field methods tend to exhibit higher numerical instability and are computationally more expensive when tested on similar design problems Hua et al. (2014). Therefore, they have not been widely employed for topology optimization problems and the reader is referred to Takezawa et al. (2010), Choi et al. (2011), Jeong et al. (2014 for more information. In addition, Sobolev gradient flows such as the H1 and L2 gradient flows are other methods of representing and evolving shapes in a topology optimization framework (Burger 2003;Sundaramoorthi et al. 2007;Henning and Peterseim 2020).
In level set methods, is defined as the level set function and implicitly defines the boundary of the structure through the following conditions: where is the structure domain, is the structure boundary and D is the design domain. Note that D defines the physical limitations of how large the structure can become, and in all cases ⊂ D . Figure 1 illustrates an example level set function and the zero level set which identifies and defined in Eq. (1).
The level set method was developed by Osher and Sethian (1988) for modeling moving boundaries encountered in many applications including combustion, computer vision, and path planning (Sethian 1999). Sethian and Wiegmann (2000) introduced a topology optimization routine where the structural boundary is defined by a level set function and updated using a stress-based heuristic criteria. A complete review of level set methods can be found in (Osher and Fedkiw 2001;van Dijk et al. 2013). The traditional algorithm requires solving the Hamilton Jacobi equation (HJE) to update using a numerical procedure such as an upwind finite difference scheme. The HJE is written as follows: where t is a fictitious time and v n represents a velocity field normal to the structural boundary defined on x ∈ . Many (1) techniques have been implemented to extend the normal velocity field such that it is redefined on the entire design domain x ∈ D (Burger 2003). The velocity field is constructed using a shape sensitivity analysis to drive the design toward a local optimum Allaire et al. (2004). As the solution converges, the velocity field on the structural boundary will converge to zero, and the boundaries will stagnate.

Problem statement
While a substantial amount of research has been devoted to optimizing thermally conductive structures, relatively little attention has been given to problems where the boundary conditions are nonlinear. In Sect. 3.1, we introduce the nonlinear phenomenology of thermal radiation. Subsequently, we define a 2D, single-objective, single-constraint example optimization problem in Sect. 3.2.

Thermal radiation
Thermal radiation is a nonlinear mode of heat transfer where the heat power radiated is directly proportional to the fourth power of the surface temperature. The Stefan-Boltzmann equation quantifies the heat power radiated from a surface across all wavelengths in the electromagnetic spectrum.
In Eq. (3), Q rad is the heat power flux radiated in units of W m −2 , is the emissivity, is the Stefan-Boltzmann constant 5.670 × 10 −8 W m −2 K −4 , and surf is the temperature field of the radiating surface. The emissivity varies between 0 and 1 and accounts for the fact that real surfaces do not absorb or emit as much energy as a blackbody ( = 1 ). While the emissivity of real surfaces is also a function of wavelength, analysts oftentimes assume it to be constant as an approximation.

2D single-objective thermal radiation problems
We consider 2D problems where a surface is radiating solely to an ambient environment and a fixed temperature boundary condition is enforced. A Dirichlet boundary condition is applied to an edge segment on the structural boundary and the radiating surface is defined as a fixed subdomain. Recall the notation defined in Sect. 2 which is used to describe the geometry. Figure 2 illustrates the design domain D , structure domain , radiating subdomain rad , and structure boundary . The structural boundary consists of two disjoint parts: the Dirichlet Γ D, and insulating Γ 0, boundaries.
The structure domain always remains a subset of the design domain ( ⊂ D ) as does the radiating subdomain always remain a subset of the structure domain ( rad ⊂ ). These relations are illustrated in Fig. 2. While our approach would be capable of solving 3D problems, the added dimension significantly increases computational cost and complexity. We consider only 2D problems here in order to validate our gradient-based approach, rather than optimize the methodology for large and computationally expensive problems.  In this paper, we define a thermally efficient radiating structure as one that can exchange as much heat power between itself and its surroundings as possible provided some restriction on material or weight. The total heat power radiated (in units of watts W ) is calculated by integrating the Stefan-Boltzmann law from Eq. (3) over the radiating surface rad . The heat power radiated functional J( rad , ) is written as follows: where = (x) is the temperature field within the design domain x ∈ D . Since rad ⊂ , the objective function J( rad , ) can also be defined without loss of generality as follows: where is a domain flag defined in Eq. (7).
The reason why the objective is transformed in this way will be apparent when the notion of shape derivatives are introduced in Sect. 5.
To limit the amount of material that can be used to design the structure, a volume fraction equality constraint is enforced. The volume fraction is defined as the volume of the structure domain normalized by the volume of the design domain D . Thus, the equality constraint function is written as follows: where V req is the required material volume. The integrals in the numerator and denominator represent the material and total design domain volumes, respectively. Finally, the governing partial differential equation (PDE) defined by Fourier's law of conduction coupled with the Dirichlet boundary and radiation conditions is written in Eq. (9).
In Eq. (9), is the second-order thermal conductivity tensor, and 0 is the prescribed temperature on the Dirichlet boundary Γ D, . We assume that the 2D elements have a uniform thickness of 1 and thus the radiating term 4 on the first line of Eq. (9) is in units of W m −3 . For isotropic materials, the thermal conductivity tensor reduces to a single term defined in units of W m −1 K −1 and also assumed to be invariant with temperature. The above objective and constraint functions along with the governing PDE can now be formally written as the following optimization problem (SP): The volume fraction equality constraint is added to formulate a well-posed topology optimization problem. Without this constraint, the optimizer would select a fully filled material design domain to maximize the radiated heat power. While the volume fraction could also be posed as an inequality constraint, any solution with an inactive volume fraction inequality constraint could be improved by the addition of material to further increase the thermal conductivity of the structure.

Methodology
In order to solve optimization problem (SP), we employ an augmented Lagrangian method to handle the volume constraint. This method is detailed in Sect. 4.1. In addition, our level set method requires an explicit definition of the level set function over the design domain. This definition is detailed in Sect. 4.2. Section 4.3 provides an explanation of how the material properties are interpolated across the entire design domain using the Ersatz material approach. And finally Sect. 4.4 describes how the HJE is solved numerically at each step of the optimizer.

Augmented Lagrangian method
A common way to solve constrained optimization problems is to replace the original problem by a series of unconstrained subproblems which are iteratively solved until a converged solution is achieved. Each subproblem seeks to minimize a single equation constructed by adding terms derived from the constraint functions to the original objective function (Nocedal and Wright 2006). Example methods which uniquely define these subproblems include penalty methods, interior point methods, method of Lagrange multipliers, and augmented Lagrangian methods.
We select an augmented Lagrangian method to solve optimization problem (SP). The augmented Lagrangian combines both the Lagrangian function and a penalty function derived from the equality constraint. Equation (10) defines the augmented Lagrangian: where J( , ) is the original objective function found in Eq. (6), h( ) is the equality constraint function found in Eq. (8), is the Lagrange multiplier estimate, and is the penalty parameter. The unconstrained subproblem to be solved at every iteration k during the optimization routine is then: where the superscript k denotes the current iteration. Recall that the optimizer is changing the structural boundaries via updates to the level set function. However, for simplicity, we indicate as the design variable describing the structural shape. Also, we have defined the temperature field as a function of the structural domain = ( ) since at every iteration, the finite element model is updated and solved using the current shape . More details on this step are provided in Sect. 4.3. The Lagrange multiplier estimate and penalty parameter are also iteratively updated to drive the design toward a feasible solution where approaches the Lagrange multiplier * that satisfies the Karush-Kuhn-Tucker (KKT) conditions. Our approach for updating these parameters follows that of Nocedal and Wright (2006).

Defining the level set function
The level set function = (x) is a scalar valued function defined on the design domain x ∈ D illustrated in Fig. 2. Equation (1) defines the structure domain where an isotropic material exists and the void domain D ⧵ where no material exists. The boundary between the material and void domains can be calculated using a variety of interpolation techniques employing an approximate Heaviside function ( van Dijk et al. 2013). In our approach, however, the level set function is interpolated within each element e via the shape functions N e i from the finite element mesh used to analyze the structure: In Eq. (12) i = {1, 2, ..., } where is the total number of shape functions in element e, e i is the ith nodal value of the level set function in element e, and T e is the domain of For our approach, we use the finite element model to define the level set function. Thus, the structure boundary is defined within each element cut by = 0 by solving the following equation for x: For quadratic and higher order elements, a Heaviside function is a useful approximation for determining the structural boundary x ∈ within each element. However, this approximation will introduce numerical errors and care must be taken to devise a suitable function which balances numerical error and smoothness (since the Heaviside function itself is non-differentiable). For the example problem presented here, we use only linear finite elements and thus the structural boundary is piecewise linear and hence straightforward to calculate.

Ersatz material technique
Two strategies for handling elements cut by the zero level set = 0 are (1) the Eulerian approach and (2) the Lagrangian approach. In this paper, we adopt an Eulerian approach where the mesh for the design domain D is kept fixed for the entirety of the optimization procedure. The Ersatz material technique then dictates that elements cut by = 0 are assigned an interpolated value for the material thermal conductivity. In the Lagrangian approach, the structural domain is re-meshed at each iteration in the optimizer such that all elements are either solid material or completely void. Void elements are removed entirely from the finite element model. Both Dapogny et al. (2014) and Allaire et al. (2014) have shown that Lagrangian methods improve accuracy but suffer from the added computational effort needed to perform a meshing operation at every iteration. Using the Ersatz material approach, a density value e is assigned to each element e similar to the design variables used in a SIMP-based approach. These values are used to interpolate the thermal conductivity within cut elements and are calculated based on the volume fraction within element e lying inside the structure domain according to Eq. (14): where e is the material subdomain of element e defined as e ∶= x | x ∈ ∩ T e , | | e | | is the area of e , and | | T e | | is the total element area since the elements are two-dimensional. In the case of a 3D problem, the areas described above would be defined as volumes.
Equation (13) provides an analytic expression for defining the zero level set using the element shape functions. The densities e linearly interpolate the element thermal conductivity matrices e . The global thermal conductivity matrix is then assembled according to the following equation: where is a small positive number and m is the total number of elements in the mesh. If an element density lies in the range 0 ≤ e ≤ , then the element matrix is multiplied by to avoid ill-conditioning or generating a singular global matrix. The analyst must make a judicious choice for which balances numerical stability and model accuracy.

Solving the Hamilton Jacobi equation
The level set function is updated at each optimization iteration by solving the Hamilton-Jacobi equation found in Eq. (2). We initialize the level set function as a signed distance function, reproduced from Dapogny and Frey (2012) and presented in Eq. (16) below: where d(x, ) is the shortest distance between the point x and the boundary . From this distance function, several metrics can be computed such as the boundary normal n = ∇ |∇ | and curvature H = div ∇ |∇ | . The normal advection velocity v n can therefore be rewritten as a vector with magnitude v and direction n: Equation (17) is substituted into Eq. (2) to yield the following form of the HJE.
We periodically reinitialize the level set function to a signed distance function in order to accurately calculate the boundary normal n(x) and curvature H(x) for x ∈ D (Allaire et al. 2004;Wang et al. 2004; van Dijk et al. 2013). These fields are a necessary ingredient to constructing the velocity field discussed later in this paper. We accomplish these reinitializations by solving the Eikonal equation (Dapogny and Frey 2012), a boundary value problem stated below: A method for solving Eq. (19) was introduced by Sethian (1999) known as the fast marching method. We employ this method to reinitialize the level set function periodically because it is significantly faster than iterative algorithms. Since the level set function needs to be reinitialized often, this advantage will play an important role in speeding up our level set-based approach. Equation (18) can be solved using a variety of upwind finite difference schemes for structured meshes (Sethian 1999;Allaire et al. 2004;Carlini et al. 2005) as well as for unstructured meshes (Sethian 1999;Zhu and Qiu 2013) and more recently using the method of characteristics (Allaire et al. 2014). We choose a second-order upwind scheme to solve Eq. (18) on a structured Cartesian grid. Details regarding this scheme can be found in Sethian's book (1999). Neumann type boundary conditions are implemented and written as follows: where D is the boundary of the design domain. While there is no clear correct choice for the level set function boundary conditions, the Neumann type seems to be the most natural and most common (Allaire and Jouve 2005;van Dijk et al. 2013). This boundary condition alleviates any artificial effect that the design domain boundary would have on the structural boundary.
The velocity field v in Eq. (18) is constructed to minimize the augmented Lagrangian stated in Eq. (10). Recall in optimization problem (SP) that the design variable is the structural shape . In order to employ a gradient-based optimizer, we require sensitivity information of the augmented Lagrangian to perturbations in the structural boundary. This sensitivity information is commonly referred to as a shape derivative. In order to improve convergence and stability, the velocity field is also regularized. We accomplish this by projecting the velocity onto the Hilbert inner product space presented by Gournay (De Gournay 2006).

Shape sensitivity analysis
A shape sensitivity analysis is performed on the radiated heat power objective functional defined in Eq. (6) and restated below: We negate the function found in Eq. (6) since our goal is to maximize the heat power radiated. We also restate the governing PDE here for reference: where the domain flag is a Boolean operator used to specify the radiating ( = 1 ) and non-radiating ( = 0 ) domains explicitly written in Eq. (9). The derivation presented here closely follows that found in Allaire et al. (2004) which the authors credit to Murat and Simon (1975) and Simon (1980). The interested reader should refer to (Pironneau 2012) and (Sokołowski and Zolésio 1992) for a more detailed review of shape derivatives. We use the adjoint method to perform the necessary shape sensitivity analysis. The adjoint method allows for efficient calculation of sensitivities for large dimensional problems oftentimes encountered in PDE constrained optimization. This process is enumerated in the following steps: 1. Define the Lagrangian L( , , q) from the objective function in (21), governing PDE in (22) and corresponding Lagrange multiplier q. 2. Find the stationary point ( * , q * ) of the Lagrangian by setting the derivatives with respect to the Lagrange multiplier q and solution to zero. These equations define the optimality conditions for minimizing the Lagrangian. 3. Calculate the shape derivative of the Lagrangian at the stationary point L( , * , q * ) , where q * is the adjoint state vector.
Beginning with step (1), the Lagrangian is defined as follows: In Eq. (23), an insulating boundary condition on Γ 0, has been imposed. Recall that the boundary is composed of both the Dirichlet and insulating boundaries according to Fig. 2: After applying Green's first identity and grouping terms, the Lagrangian is transformed into Eq. (25).
In step (2), we find the stationary point ( * , q * ) where the partial derivatives of the Lagrangian with respect to q and are set equal to zero. The derivative of the Lagrangian with respect to q in the direction ∈ H 1 ( ) , where H 1 ( ) is the Sobolev space whose functions and first-order weak derivatives are both Lebesgue square-integrable, is given below: Using Green's first identity again to rewrite the first integral and collecting terms on the Dirichlet and insulating boundaries, Eq. (26) is transformed as follows: Setting Eq. (27) equal to zero yields the original state equation and boundary conditions found in Eq. (22). The first integral gives the state equation for any with compact support in . Varying the trace function on Γ 0, gives the insulating boundary condition and by varying ( ∇ ) ⋅ n on Γ D, gives the Dirichlet boundary condition. Next the derivative of the Lagrangian with respect to the state variable is used to define the adjoint equations.
Again by using Green's first identity and grouping terms on the Dirichlet and insulating boundary conditions, Eq. (28) is transformed as follows: Setting Eq. (29) equal to zero yields the adjoint state equations. The first integral gives the adjoint state equation for any with compact support in : Similarly, the insulating boundary condition is given by varying the trace function on the boundary Γ 0, : Finally, the Dirichlet boundary condition is given by varying ∇ on the boundary Γ D, : Note that q * is the adjoint state vector. In summary, the adjoint state is defined by the following set of equations: While the thermal radiation problem is shown to be well posed, it is not self-adjoint (i.e., q * ≠ ). In order to calculate the shape derivative, we need to first solve the original governing PDE in Eq. (22). Next, the solution is substituted into Eq. (33) as * and the equation is solved for q * . We use the finite element method to solve both PDEs numerically.
Finally, the shape derivative of the Lagrangian defined in Eq. (25) at the stationary point ( * , q * ) is calculated. Allaire et al. (2004) provide the general forms of shape derivatives for domain and boundary integral functionals which are used to construct the following: Splitting terms into the Dirichlet and insulating boundaries, and using Eqs. (22) and (33) to enforce * − 0 = 0 and q * = 0 on the boundary Γ D, , the above equation is simplified as follows: In practice, additional simplifications are made to shorten Eq. (35). First, the optimizer should not be removing material from the Dirichlet boundaries since this area is the source of the heat power being radiated by the structure. Isolating or removing material in this region would cause the temperature on the radiating subdomain to drop, causing the radiated heat power to decrease. Therefore, on Γ D, , the optimizer is forced to retain the structural boundary and the dot product ⋅ n goes to zero, eliminating the second integral in Eq. (35).
Using a similar reasoning, the radiation flag is only active ( = 1 ) on the subdomain rad . This subdomain should also remain fixed since removing radiating material will reduce the radiated heat power. As such, the optimizer is only affecting regions with (x) = 0 where there is no thermal radiation. Note that regions and boundaries that are not able to be modified by the optimizer are termed "non-design," meaning they are not subject to the removal or addition of material. By enforcing these practices, Eq. (35) is reduced to the following single term: As previously mentioned, solutions for both * and q * are required to calculate the shape derivative. These solutions are found numerically using the finite element method. The sensitivity field for the radiating heat power functional defined in Eq. (36) is also calculated using the finite element method according to Eq. (37): where the subscript e identifies the local element, e is the element adjoint vector, and e is the element nodal temperatures. Equation (37) is normalized by the element area |T e | due to the fact that the finite element equations are written in the weak form, and hence the term T e e e is an integral over the element volume.
where B is the strain matrix of shape function derivatives, D is the constitutive matrix, and Je is the element Jacobian matrix. Equation (38) must be normalized by |T e | to arrive at the shape sensitivity at the element centroid.

2D case study
We present a 2D problem to demonstrate that the optimizer is converging to locally optimal solutions and validate that the shape sensitivity is correct. Modeling, analysis, and optimization are all performed within the MATLAB ® environment. All numerical analyses are performed and timed on a Dell laptop equipped with an Intel ® Core TM i7-4700MQ processor running at 2.40 GHz , and 16.0 GB of DDR3 SDRAM. All codes were compiled and run using MATLAB ® version 8.6 (R2015b). Section 6.1 defines the problem geometry, boundary conditions, and material properties. Section 6.2 provides an overview of the finite element model used to evaluate the structural performance and sensitivities. Then in Sect. 6.3, we illustrate several candidate solutions, demonstrating that the thermal compliance objective function is insufficient for solving this class of problems.

Problem definition
The problem analyzes a thin square plate with dimensions 100 cm × 100 cm and thickness 1 cm . A Dirichlet boundary T e e e = T e ∫ T e T det e dV e condition is enforced on a 20 cm segment along the center of the left edge where the temperature is fixed at 300 K . The plate rejects heat via thermal radiation on a domain of size 4 cm × 100 cm located on the right side of the square plate. This subdomain rad is the only region where heat is being radiated. Thermal radiation is emitted from the surface of this subdomain rather than the edges, accounting for a total of 800 cm 2 when considering both front (out of the page) and back (into the page) surfaces. For simplicity, the radiating surface is modeled as a perfect blackbody, meaning the emissivity = 1 and the plate does not absorb any energy from the surrounding environment (assumes that the mean ambient temperature is 0 K ). Figure 3 illustrates the problem geometry. In Fig. 3, the gray region indicates the radiating subdomain rad as well as the non-design subdomain (i.e., the radiating material is fixed during the optimization and cannot change shape). This is a common tactic used in the topology optimization community to enforce a design-independent load condition. Some authors have also considered designdependent load conditions where the boundary conditions are affected by the structural shape (Bourdin and Chambolle 2003;Allaire et al. 2004;Iga et al. 2009;Yamada et al. 2011). However, design-dependent load conditions are computationally challenging to implement and are not necessary for the purposes of verifying the thermal heat-radiated sensitivity analysis.
The plate is made from an isotropic material (aluminum alloy) whose physical properties are listed in Table 1.  Note that the material properties are identical for both the design and non-design domains.

Finite element model
A finite element model is used to solve the governing PDE from Eq. (9) for the 2D radiating plate problem illustrated in Fig. 3. For an in-depth review of finite element methods, the interested reader is referred to Bathe (2014). We employ MATLAB ® 's Partial Differential Equation Toolbox to handle finite element modeling and analysis. A structured mesh is generated using 5000 linear triangle elements including 2601 nodes.
The finite element program solves the following global linear system of equations: where is the p × p global thermal conductivity matrix, is the p × 1 temperature vector, and is the p × 1 heat load vector. The dimension p is the number of nodes (and degrees of freedom) in the finite element model. MATLAB ® 's PDE Toolbox implements a Gauss-Newton solver since the heat load vector is a nonlinear function of the temperature according to the Stefan-Boltzmann law. The global conductivity matrix is assembled from the elemental conductivity matrices according to Eq. (15). The heat power radiated is calculated using Gauss quadrature to integrate the fourth power of the temperature in each radiating element.

Efficient thermal design using a compliance objective function
Various other researchers, while not tackling problems that incorporate thermally radiating boundary conditions in particular, have investigated the design of efficient thermally conductive structures (Gao et al. 2008;Iga et al. 2009;Burger et al. 2013;Liu and Tovar 2014;Zhu et al. 2016).
The common objective function used in these problems is a thermal compliance function written as follows: where ∇ is the temperature gradient vector field. Equation (40) is analogous to the structural compliance objective where ∇ and are replaced, respectively, by the secondorder strain tensor and fourth-order elasticity tensor (Allaire et al. 2002). The thermal compliance functional is a measure of how well the structure conducts heat. By minimizing the thermal compliance, the optimizer is minimizing the temperature gradients and hence keeping the structure as isothermal as possible. Thus, the goal of these problems is to design for thermal efficiency since a limited amount of material is allocated to minimize thermal resistance and hence transport heat as effectively as possible.
A classic example is the 2D design of a thermal heat sink in a square design domain. This problem has been solved using both SIMP (Bendsøe and Sigmund 2004;Maute 2014) and level set methods . In this problem, a uniformly distributed heat load is imparted on the design domain, with a fixed temperature on a segment of one edge. The solution is a symmetric branched tree structure which attempts to capture as much of the distributed heat load and transfer it to the Dirichlet boundary. This classic problem has also been solved in 3D by Burger et al. (2013).
We find that minimizing the thermal compliance of a thermally radiating structure does not necessarily maximize the heat power radiated. If these two objective functions were negatively correlated, we could pose the optimization problem using a thermal compliance objective function. However, we demonstrate a counterexample where the two objective functions are instead positively correlated. Consider the three test designs illustrated in Table 2, each with a volume fraction of 0.5, along with their corresponding finite element solutions. The boundary conditions and material properties are the same as those presented in Sect. 6.1. Each design is evaluated for thermal compliance and heat power radiated.
The material and void domains in the illustrated shapes are designated by the cyan and white colors, respectively. The radiating strip is shown in gray on the right side of the design domain as was depicted in Fig. 3. Note that although the third design (bottom row in Table 2) completely disconnects the material domain from the radiating domain, the structure is still able to radiate some heat. Using the Ersatz material approach described in Sect. 4.3, the void domain retains a small nonzero thermal conductivity. The thermal conductivity of the void region is 10 3 times smaller than the nominal material conductivity.
It is clear from Table 2 that there is a positive rather than negative correlation between the thermal compliance and heat power radiated functions of the three test designs. The designs in Table 2, which exhibit lower thermal compliance also have a reduced capacity to radiate heat. This finding debunks our hypothesis that these two objective functions are strictly negatively correlated, such that a lower thermal compliance design would always be radiating a larger amount of heat. While there are possible cases in which these two objectives do negatively correlate, it is clear from the above counterexamples that this is not always guaranteed. As such, it is crucial to develop a topology optimization framework which specifically defines a radiated heat power objective function and its corresponding shape derivative.

Numerical results
A finite element analysis is first performed on the plate with all of the design domain filled by material. The results of the finite element analysis are shown in Fig. 4. The temperature on the Dirichlet boundary is fixed at 300 K and heat flows from this boundary to the radiating domain, where it is then radiated out to the ambient surroundings. The average temperature of the radiating boundary is approximately 278 K and the total heat power radiated from this surface is 27.25 W.
We also evaluate two designs each with a volume fraction of 0.5 to serve as baselines to compare with those generated by the topology optimization method. They are shown in Fig. 5. Test Design 1 and Test Design 2 have total radiated heat powers of 22.58 W and 22.90 W , respectively. Because the boundary conditions are symmetric about the horizontal mid-line of the design domain, both design shapes also exhibit this symmetry. The intuition leading to these design shapes is that thermal conduction is a linear phenomenon, and the goal is to maximize the radiated heat power by designing a structure with minimal thermal gradients between the Dirichlet boundary and radiating subdomain. As such, keeping the thermal pathways as direct (i.e., as short) as possible to the radiating subdomain should yield a well-performing structure. Again these two designs will be revisited later to compare to the results from the topology optimization method.
The topology optimization problem (SP) is solved with a volume fraction equality constraint of 0.5. The first case is initiated with a fully filled design domain, and a Lagrange multiplier estimate and penalty parameter of 0 = 1 and 0 = 10 , respectively. The initial and final shapes are illustrated in Fig. 6. The optimal shape radiates 24.1 W of heat power. Thus, the structural weight is reduced by 50% while the heat rejected drops by 11.5% from the fully filled design.  The optimal design is analyzed to show the temperature and heat flow fields in Fig. 7.
The objective, constraint and augmented Lagrangian histories are given next in Fig. 8. From Fig. 8b, the shape converges after approximately 45 iterations.
The prior oscillations are due to several factors including that (1) the initial design is infeasible ( h( ) ≠ 0 ), (2) the Lagrange multiplier estimate is updated at each iteration using a monotonically increasing penalty parameter, and (3) the iterative updates to the structural boundary allow for a small increase in the augmented Lagrangian up until iteration 50. Allowing the augmented Lagrangian to increase in the first few iterations allows the optimizer to find a locally optimal solution that is sufficiently far, or topologically different, from the initial design. The number of iterations in which the augmented Lagrangian is allowed to increase is a heuristic parameter. Allaire et al. (2004) mention that different choices of this parameter can yield different locally optimal solutions. In the above analysis, we implement a backtracking line search method to speed convergence. Our backtracking routine iteratively reduces the number of time steps T and step size t to advance the level set function until the updated augmented Lagrangian is reduced. If the optimizer selects an update corresponding to the maximum number of time steps T = T max , the next optimization iteration will attempt to advance the level set function by a larger number of time steps. Conversely, if the optimizer selects a single time step, T = 1 , the next iteration will use a smaller time step. This technique is similar to the trust region methods discussed in Nocedal and Wright (2006). The time step size is constrained by the Courant-Friedrichs-Lewy (CFL) condition since we are using an explicit time integration scheme to solve the HJE.
For a converged solution, the gradient of the augmented Lagrangian should be approximately zero. This criterion is used in augmented Lagrangian methods to determine that a KKT point has been reached. Because the velocity field is constructed from the derivative of the augmented Lagrangian with respect to a perturbation in the normal direction of the boundary, we expect the velocity at the boundary of the shape to approach zero (Allaire et al. 2004). Figure 9 illustrates several metrics of the velocity field on the structural boundary over the optimization history.
In Fig. 9a, the metrics labeled as "Min(v )" and "Max(v )" refer to the minimum and maximum value of the velocity on the boundary, respectively. The normed metrics are defined using Eq. (41) for the L 1 norm and (42) for the L 2 norm: where | | is the perimeter of the structural boundary . Intuitively, as the velocity field approaches zero on the boundary , updates to the level set function from solving the Hamilton-Jacobi equations will cause smaller changes in the shape. Therefore, the structure domain will ultimately converge.
From the plots in Fig. 9b, it is clear that a converged solution has been achieved since the velocity field norms are asymptotically becoming smaller. Note that the velocity fields will never achieve zero due to numerical errors inherent in the level set method. One of the primary errors arises from interpolating the shape sensitivity field from the center of each mesh element to the node points (Dijk et al. (2012)). The interpolated sensitivity field is then used to construct the velocity field. Additionally, the reinitialization of the level set function using the fast marching method is not a volumepreserving procedure, meaning that the structural boundary will shift by some small amount. In the analyses considered here, the level set function is reinitialized every other time step taken by the HJE solver.
Since a gradient-based optimization method is employed, the solution shown in Fig. 6 is only guaranteed to be a local optimum. Two additional analyses are now presented, each with different initial conditions. Specifically, 2 × 2 and 4 × 4 arrays of circular holes are introduced in the interior of the material domain to allow for voids in the final design solution. Because the step sizes for our explicit HJE solver must be limited by the CFL condition coupled with the frequent reinitializations of the level set function, it is impossible to nucleate new holes in the interior of a material domain (Allaire et al. 2004). By initiating the level set function with voids in the interior, locally optimal solutions not achievable with a fully filled initial design domain can now be reached. We choose symmetric hole patterns for the initial material configurations since the boundary conditions and previous optimal design exhibit a line of symmetry at the horizontal mid-line of the design domain.
The optimal shape for the initial design containing a 2 × 2 array of holes is illustrated in Fig. 10. The corresponding finite element solution for the temperature field and heat flux of the optimal design are also presented in Fig. 11. The optimal solution retains several voids from the initial design. These voids align with the direction of heat flow and place material farther from the horizontal axis of symmetry. Although the complexity of this design has increased, the heat power radiated has also increased by 0.18 W from the optimal solution found using an initial fully filled design domain. The objective, constraint, and augmented Lagrangian histories are plotted in Fig. 12. These histories show that after approximately 15 iterations, the design converges. The same velocity field metrics previously introduced are also tracked to illustrate that convergence has been reached. Figure 13 plots the histories of these metrics Fig. 11 Finite element solution of optimal design for optimization problem (SP) with initial material configuration of 2 × 2 array of equally spaced circular holes Fig. 12 Objective, constraint, augmented Lagrangian, and Lagrange multiplier estimate histories for optimization problem (SP) with initial material configuration of 2 × 2 array of equally spaced circular holes and demonstrates that the velocity approaches zero on the structural boundaries. Next the results for the initial material configuration containing a 4 × 4 array of holes are displayed in Figs. 14, 15, 16, and 17. Our method converges to a final design which again retains some of the voids from the initial shape. Having increased the number of holes in the initial shape has resulted in more voids in the final design and more intricate structural features near the radiating surface.
A summary of the two test design shapes along with the three topologically optimized solutions, designated as "Design A," "Design B," and "Design C," is provided in Table 3.
All of the topologically optimized designs outperformed the two test designs presented earlier. Design A achieved a 3.5% and 2.2% improvement in radiated thermal heat power from Test Designs 1 and 2, respectively. Likewise, both Designs B and C achieved a 4.3% and 2.9% improvement in performance over Test Designs 1 and 2, respectively.
Another important finding here is that the added topological complexity in Design B and C offer only a 0.8% improvement in performance (in terms of radiated heat power) over Design A. Additionally, both Designs B and C are iso-performing, meaning they both radiate the same heat power. Therefore, the additional geometric complexity found in Design C compared to Design B in no way improves the performance. This finding suggests that design problems exist where too much geometric complexity could stunt or even degrade performance. Whether or not added geometric complexity affords better performing designs is problem-dependent.

Conclusions
In this paper, we presented a level set-based topology optimization method to maximize thermally radiated heat power. Using the adjoint method, we derived a shape derivative for this class of problems which was necessary for our gradientbased solver. The topology optimization method was then demonstrated on an example 2D problem of an aluminum plate radiating heat to the surrounding environment. The problem was well posed by including a single volume fraction constraint. We solved this problem using three different initial material configurations, each yielding a distinct locally optimal design. Our method was shown to converge smoothly by tracking the iteration histories of the objective, constraint, augmented Lagrangian, and Lagrange multiplier. The L 1 and L 2 norms of the boundary velocity also indicated that the solutions had converged. The resulting locally optimal geometries indicated that adding geometric complexity, in the form of voids and branching features, can in fact improve the performance. However, the performance benefit gained by increasing topological complexity can also stagnate. Because of these findings, it is important that we consider methodologies to explore the design tradespace in order to not only identify new unintuitive topological features, but also better understand their cost-benefit between geometric complexity and performance.
Acknowledgments Distribution Statement A. Approved for public release. Distribution is unlimited. This material is based upon work supported by the United States Air Force under Air Force Contract No. FA8702-15-D-0001. Any opinions, findings, conclusions, or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the United States Air Force. This work was also supported in part by the NSF Computational and