Temperature-constrained topology optimization of nonlinear heat conduction problems

This paper presents topology optimization of nonlinear heat conduction problems with multiple domains and multiple constraints, including regional temperature and material volume for reducing temperature. Maximum approximation temperatures in the constraint regions are accurately and dynamically calculated, though temperature and temperature-dependent thermal conductivity change with the update of material distribution. A temperature measure with constant error to approximate regional maximum temperature is adaptive to different temperature ranges. A strategy of hole nucleation generation combined with the regional temperature constraints is presented for the level set-based topology optimization. The solid isotropic material with penalization (SIMP) and parametrized level set methods are compared for the temperature-constrained topology optimization. Finally, several numerical examples are solved by the SIMP and parametrized level set methods. The results demonstrate that the proposed approach can obtain intricate topological details and reduce regional temperatures for the nonlinear heat conduction problems.


Introduction
Temperature-constrained topology optimization of heat conduction problems has many promising industrial applications such as printed circuit board substrate design and electric motor covers in the power electronics field (Dbouk, 2017;Lohan & Allison, 2017). With the development of electronic equipment, the temperature of electronic components below an acceptable value was crucial to ensure the performance of equipment (Ye et al., 2002;Hinton, 2007;. Attaching the conductivity materials was efficient to control the temperatures surrounding the material region (Bejan, 1997). Zeidan et al. (2016) discussed the optimal distribution of a limited amount of phasechange material to keep the constant temperature of a CPU. Li et al. (2004) discussed the evolutionary topology optimization (ESO) for temperature reduction. The optimal results can yield substantial improvements in the designs of thermal systems such as a wing box, a thermal diffuser, and cooling fins. The design temperature within several constraint regions must be reduced to a fixed value, which belongs to the topology optimization problems with multiple domains and multiple constraints. In this paper, the nonlinear problems with the temperature constraints are investigated for topology optimization of heat conduction problems. The main contributions of this work are as follows: (1) reduce temperatures in the constraint regions using topology optimization for the nonlinear heat conduction problems; (2) regional temperature measure for approximating maximum temperature as a constraint condition; (3) a strategy of hole nucleation generation considering the regional temperature constraints; and (4) adaptive bounds of design variables for strict control of volume changes between successive hole nucleation generations to avoid local minimum.
Topology optimization related to the linear heat conduction was widely reported in recent years. The nonlinear problems have not received as much attention since the sensitivity analysis and finite element iteration procedure are more complex than the linear problems. The numerical examples related to the linear steady-state heat conduction problems were given with the solid isotropic material with penalization (SIMP) model (Sigmund, 2001;Bendsøe & Sigmund, 2003). A framework for topology optimization of nonlinear steady-state heat conduction was provided based on the element densities (Brunes, 2007). Kim et al. (2009) combined the topological derivative into the level set model for the nonlinear heat conduction problems. Tang et al. (2019) considered the nonlinear problems involving large temperature gradients. In the most literature reviews, the objective function of nonlinear heat conduction problems was defined as thermal compliance, and the constraint was the material volume. The temperature constraints were not considered in the design domain. The optimization design with multiple conditions, including the regional temperature and volume constraints, is not a trivial matter with the level set-based shape optimization.
A regional temperature measure for the regional temperature constraints was used for the constrained topology optimization of transient heat conduction problems (Zhuang & Xiong, 2015). Lohan and Allison (2017) discussed the temperature constraint formulations. The softmax function enforced a global temperature constraint. It also evaluated regional restrictions in two power electronics motivated design problems, which illustrated topology optimization with the regional temperature constraints. This paper emphasizes the nonlinear problems, and thermal compliance is taken as the objective function subject to the volume and regional temperature constraints. The temperature-dependent thermal conductivity is given with exponential law.
The homogenization or its simplified version called SIMP (Suzuki & Kikuchi, 1991;Sigmund, 2001;Bendsøe & Sigmund, 2003), ESO (Xie & Steven, 1993;Huang & Xie, 2010), and level setbased approaches (Wang et al., 2003) are the main topology optimization methods. These methods are applied to the structural design of elasticity and heat conduction problems. The SIMP and ESO use the element information as the design variables. The level set-based method uses the shape of the material domain as the design variable. Therefore, the density-based SIMP method as the discrete framework is adopted to topology optimization for the nonlinear problems since it is popular and physically permissible (Sigmund, 2001). The level set-based method as the continuum framework is selected for topology optimization, and the optimal results are compared with the results of the SIMP method. A range of topology optimization problems are solved by the following approaches such as the optimality criteria (OC; Kikuchi et al., 1998;Nishiwaki et al., 1998), the method of moving asymptotes (MMAs; Svanberg, 2002), and the sequential linear programming method (Gomes & Senne, 2014). Most works on topology optimization of heat conduction problems treated the optimization problems with the linear heat conduction and volume constraints (Joo et al., 2017;Yan et al., 2018), so the SIMP model with the OC algorithm was popular to solve these optimization problems. The level set-based topology optimization was also applied to the heat conduction problems. Ha and Cho (2005) derived the shape sensitivity using the adjoint sensitivity method for the steady-state heat conduction problems. Zhuang et al. (2007) discussed the topology derivative and the multiple load cases for the level set-based topology optimization of the steady-state heat conduction problems. Since the MMA algorithm has proven itself to be versatile and is an effective largescale constrained nonlinear optimization solver (Le et al., 2010;Kuo & Cheng, 2017;, the MMA algorithm is used to solve the temperature-constrained topology optimization for the nonlinear problems. Therefore, the optimization models and the numerical solver selected in this work are typical and efficient, ensuring the effectiveness of comparison between the proposed topology optimization methods. Unlike the gradient-based algorithms, the metaheuristic algorithms based on successful evolutionary behavior of natural systems are nature inspired, which has been applied to structural optimization. All metaheuristic algorithms use some trade-off of local search and global exploration (Gandomi et al., 2013). The metaheuristic algorithms include genetic algorithms (GAs), simulated annealing, and particle swarm optimization. Their performances are better than simple heuristics for topology optimization. The metaheuristic algorithms have been applied to topology optimization of elasticity and heat conduction problems. These algorithms have the merits of being reliable and robust in terms of global convergence. Savsani et al. (2017) proposed truss topology optimization with the static and dynamic constraints by integrating the random mutation-based search technique. The GA-based topology optimization was performed to minimize the weight and the compliance (Wang & Tai, 2005). It achieved better performance compared with the density-based approach due to the strategy of global search. The binary bat algorithm for structural topology optimization achieved the stiffness design of the structure (Jaafer et al., 2020). A global optimization metaheuristic based on the probabilistic learning, namely estimation of distribution algorithm, was proposed for topology optimization. The algorithm did not require initialization and a priori information (Valdez et al., 2018). Boichot and Fan (2016) used the GA to solve the classical volumeto-point heat conduction problems. The material distributions of the conductive paths were designed by GA with maximum and mean temperature minimizations as the objective functions. A new GA-based topology optimization combined with computational fluid dynamics was applied to create fin shapes for the heat exchangers in the aerospace applications (Mekki et al., 2019). Firefly algorithm (FA), initially developed by Yang (2009Yang ( , 2010, was a stochastic, metaheuristic algorithm. This method was applied to solve the optimization problems. An FA-based hybrid method for the stress-based topology optimization was proposed to generate a global convergence solution (Gebremedhen et al., 2020). For more details related to the metaheuristic algorithms with their applications, the reader may refer to Gandomi et al. (2013).
As discussed above, the topology optimization related to the steady-state heat conduction problems with volume constraint has been widely investigated in recent literature. In contrast, the nonlinear heat conduction problems with temperature constraints have not received as much attention, even though the regional temperature control has a variety of potential applications in industrial fields. Topology optimization has been applied to aeronautics and aerospace engineering (Zhu et al., 2016). Thus, large temperature variations should be considered. For example, a single part can endure a temperature gradient ranging from 100 to 500 • C/cm (Tang et al., 2019). The thermal conductivity measured at 600 • C is nearly three times higher than that measured at room temperature for material TC4 (Tang et al., 2019). The linear and nonlinear dynamic characteristics of plates with thermal effect have been widely investigated in the spacecraft and nuclear industrial fields, because of their excellent performances of withstanding severe high-temperature gradient while maintaining structural integrity (Yang & Huang, 2007). Therefore, the nonlinearity of material properties should be considered to achieve good designs with the desired performance (Ha & Cho, 2008;Tang et al., 2019). This paper focuses on topology optimization of the nonlinear heat conduction problems with regional temperature constraints. The multiple regional maximum temperatures and material volume are taken as the constraint conditions. The conventional level setbased topology optimization cannot easily be extended to the numerous constraint optimization problems, when the level set equation is solved by an upwind scheme (Wang et al., 2003). Thus, the parametrized level set method combined with the MMA algorithm is applied to solve the topology optimization for the nonlinear heat conduction problems with multiple constraints. Wei et al. (2018) used the first-order forward Euler's method to update design variables. The derivatives of the coefficients of radial basis functions (RBFs) concerning time were employed for the optimization algorithm. Different from the above algorithms, this paper employs the derivatives of the objective function and regional temperature concerning the coefficients of RBFs. Then, the global MMA algorithm is used to obtain the optimal results for updating material boundaries.

Nonlinear heat conduction problems
The model of the temperature-dependent thermal conductivity is selected according to the numerical examples in references (Lienemann et al., 2003;Filipov & Farago, 2018;Sun & Li, 2020). Such a temperature-dependent nonlinear model conforms to real physical systems, e.g. for silicon (Filipov & Farago, 2018). The thermal conductivity is approximated by an exponential law (Lienemann et al., 2003;Filipov & Farago, 2018;Sun & Li, 2020): where κ 0 , β, and u β are three constants and u is the temperature in the design domain. The nonlinear heat conduction is transformed into linear heat conduction when β is set to 0, as shown in Fig. 1. The thermal conductivity κ(u) is a nonlinear temperature-dependent function. The material distribution changes with the optimization iterations, which generates the temperature field variation over the entire design domain. At each iteration, the thermal conductivity should be recomputed according to the temperature field until a stopping criterion is satisfied. As shown in Fig. 2, the formulations of heat conduction are given as where u is the solution of nonlinear heat conduction problems, κ(u) is the thermal conductivity coefficient of the unknown temperature u, and q v is the inner heat generation rate per unit volume. Three kinds of boundary conditions, namely a temperature condition on 0 , a heat flux conduction on 1 , and a convective condition on 2 , are imposed; q is a heat flux imposed on 1 , is the heat transfer coefficient, u 0 and u f are constants, and n is an outward unit normal vector on the boundary. For topology optimization, new holes may be created, and the interior boundaries may merge. Thus, the newly created convection boundary should be specially considered, which usually belongs to the design-dependent problems (Chen & Kikuchi, 2001;Joo et al., 2017;Yoon & Koo, 2019). Furthermore, when the convection boundary conditions on the designable structural boundaries are involved, the normal directions and the curvatures of boundaries must be calculated on the boundaries. The conventional element-density-based framework cannot handle this design-dependent problem easily. For such issues, the independent point-wise density interpolation method (Kang & Wang, 2011;Wang et al., 2013), or the level set-based method (Zhuang et al., 2007), can be employed to provide finite elementsindependent and boundary description, which enables the extraction of smoother structural boundaries and normal directions/curvatures. The convection boundaries are not considered in this work. For completeness, the optimization formulations and derivatives include the convection boundary condition in the following sections.
Using the virtual temperature field ω, the weak form of governing equation is denoted as where a(u, ω) and (ω) are the thermal energy form and the linear load form, respectively: The relation of equation (3) is written as the finite element equation over the element domain e : where N is the shape function, B is the temperature differentiation matrix with B = ∇ N, and u e denotes values at the selected element nodes. The global finite element equation can be assembled from the finite element equations. Using equations (6) and (7), the nonlinear global finite element equation of equation (3) is expressed as where K(U) is the nonlinear function concerning the unknown nodal temperatures U, U is the solution to be determined, R E is the vector that contains all information concerning the heat loads acting on the design domain, and R ε (U) is the residual. The approximation to the solution can be obtained by solving the assembled equations with the Newton-Raphson method.

Topology optimization problems
When the element densities are taken as the design variables, the topology optimization problem is formulated as where V and V 0 are the material volume and design domain volume, respectively, γ is the volume ratio, ρ is the vector of design variables, ρ min is a vector of minimum relative densities, ζ max is the maximum temperature of the prescribed constraint regions, and u ζ describes an upper limit of the permissible temperature. The second inequality in equation (9) illustrates the re-gional temperature constraints. The maximum approximation measure originated from Wen (1989) is defined as where u i is the selected nodal temperature with i = 1, 2, . . . , n c in the constraint regions, is a constant that makes ζ max approximate the maximum temperature as → +∞, and ζ 0 is an error compensation term. More details related to equation (10) would be discussed at the end of this section.
In the continuum framework, the objective function is defined on a physical domain: where represents the physical domain; the maximum approximation temperature in the continuum framework is given as where ϕ is the signed function to represent the temperature constraint regions and H(ϕ) is the Heaviside function (Wang et al., 2003). Equation (12) is the maximum approximation temperature described as the domain integration of equation (10). The temperature constraints in equations (9) and (11) usually belong to the min-max problems, restricting the temperature below a prescribed value in the constraint regions. Thus, the maximum approximation of temperature is crucial and affects the optimal results.

Maximum temperature approximation
The global maximum approximation function is developed initially to tackle the min-max problems. Here, three different global measures to approximate the maximum temperature are investigated and compared. The V a approximation (Wen, 1989), smooth approximation (Tsoukalas et al., 2009), and p-norm (Le et al., 2010) are compared in Fig. 3. Tsoukalas et al. (2009) generalized the smoothing technique with a maximum function: A global solution to approximate the maximum temperature can be obtained as → +∞. The global p-norm measure for maximum approximation function is differentiable and smooth (Le et al., 2010): The p-norm measure has been used for the stress-based topology optimization. The three different global measures are tested with varying ranges of temperature and different .
To show the performances of the maximum approximation functions, the temperature ranges from 0 to 20 and is evenly divided into 20 intervals. Three maximum approximation functions are tested in each interval. Each interval has 300 random temperatures that are generated using a uniform distribution function. From Fig. 3, the error decreases with increment of the parameter from 15 to 90. However, the error of maximum temperature approximation of the Le's method with equation (14) increases with the increment of temperature marked as the  (13) corresponding to the Wen's and Tsoukalas' methods generate stable errors for maximum temperature approximation. However, Tsoukalas' function returns the scalar representation of positive infinity when specifying u > 10 and ≥ 45, as shown in Fig. 3c-f. Using random temperatures with the Gaussian distribution, the performances of the maximum approximation functions are similar to the performances using the random temperatures with uniform distribution, as shown in Fig. 3. According to the above analysis, this paper uses the Wen's model to approximate the maximum temperature. For approximation, a constant ζ 0 is introduced in equation (10) for theoretical completeness. However, the compensation constant ζ 0 does not arise during the densitybased and continuum shape-based sensitivity analyses.

Density-based sensitivity of nonlinear heat conduction problems
Using the Lagrange multiplier method, a completely equivalent problem to the original optimization problem is written as where λ is the Lagrange multiplier. Using equations (1) and (6), the derivative of the first term on the right-hand side of equation (15) concerning element density ρ i can be calculated as where the third term on the right-hand side of equation (16) is generated by the nonlinear term from equation (1). The derivative of the second term on the right-hand side of equation (15) is given as The adjoint variable method is used to compute the design sensitivity, and the following adjoint equation is introduced: where (15) yields the total derivative off (U, ρ): To derive the derivative of the maximum approximation temperature ζ max of equation (10), the adjoint variable λ ζ is introduced after normalization of ζ max , yielding The derivative of the first term of the right-hand side of equation (20) is given as Substituting equations (17) and (21) into (20), the total derivative ofζ max concerning ρ i is given as where λ ζ is the solution of the following adjoint equation: To obtain the derivatives of the objective and temperature constraint functions, two adjoint equations (18) and (23) should be solved at each optimization iteration.

Continuum shape sensitivity of nonlinear heat conduction problems
The material derivative formulas were derived from the domain and boundary functionals (Pironneau, 1983;Sokolowski & Zolesio, 1992). The material derivative of J ( ) in (11) is given as where κ c = div(n) is the curvature of the boundary 2 , and V n is the normal component of design velocity on the boundaries. To obtain the explicit formula of shape sensitivity, the material derivative of the weak form equation (3) is given as where a (u, ω) and (ω) are of the forms and Substituting equations (26) and (27) into (25) yields An adjoint equation is introduced by replacing ω in equation (28) with a virtual temperature fieldλ, yielding Substituting equations (28) and (29) into (24), the material derivative of the objective function J ( ) is expressed as Therefore, the shape sensitivity analysis, which indicates the effect of shape variation on the objective function, can be obtained in terms of the resulting boundary integrals. Consider the regional temperature constraints in equations (11) and (12), the maximum approximation temperature is written as the normalized form: The material derivative ofζ max (ϕ) is given as where the parameters η, A, and B are defined as To obtain the explicit expression of equation (32) concerning V n , using equation (28), the following adjoint equation is whereλ ζ is the virtual temperature field. Using equations (28) and (34), the following sensitivity expression ofζ max (ϕ) is obtained: The fact that the continuum sensitivities J ( ) andζ max (ϕ) can be expressed as a boundary integral gives a significant advantage in numerical implementation. In Section 3.1, the density-based design sensitivity is given when the element densities are taken as the design variables. This section provides the continuum shape sensitivity for the temperature-constrained topology optimization problems.

Multiconstrained Optimization Algorithms
The nonlinear optimization problems are subject to the regional temperature and volume constraints. Therefore, the MMA algorithm is proposed to solve the topology optimization problems with multiple constraints. The nonlinear heat conduction problems and topology optimization algorithms are given with the density-based and continuum frameworks. The parametrized level set method is adopted to update the shape of the material domain, since the conventional level setbased topology optimization solved by the upwind difference scheme cannot easily handle the multiconstrained optimization problems.

Density-based optimization algorithm
With the SIMP scheme, the thermal conductivity can be interpolated according to the relative element-wise constant density and its derivative is given by where κ min is the thermal conductivity of void phase to avoid singularity of the finite element stiffness matrix, and p is the penalization power. Equation (37) is substituted into (19) and (22) to obtain the sensitivity results for topology optimization.
The standard Newton-Raphson iterative scheme is applied to solve the nonlinear system. The tangent stiffness matrix of the nonlinear system equation (8) is given as The temperature field is updated with a sequence of systems in the following form: The convergent solution is obtained when the termination condition R E − K(U i ) ≤ ε is satisfied, for some ε > 0.

Algorithm 1: Density-based optimization algorithm
Step 1: Initialize and discretize the design domain with finite elements. The material parameters for the nonlinear heat conduction problems in equations (1) and (2) are given. The parameters γ , ρ min , , u ζ , p, κ 0 , κ min , β, and u β for volume constraint, thermal conductivity, and maximum approximation temperature from equations (9), (10), and (36) are chosen for the optimization problems.
Step 2: Find the temperature field u and the adjoint variable λ by solving equations (3) and (18). The derivative off (U, ρ) concerning ρ is then solved by equation (19). For the maximum approximation temperature ζ max , solve the adjoint equation (23) to obtain the adjoint variable λ ζ and calculate the derivative ofζ max concerning ρ using equation (22).
Step 3: The constants for the MMA algorithm are selected as inputs. The constraint functions V and ζ max in equation (9) should be normalized to restrict the constraints to the reasonable ranges that satisfy the input requirements of the MMA algorithm.
Step 4: Solve the topology optimization problems subject to the volume and regional temperature constraints with initial iteration m = 0. Update the design variables ρ i with the MMA algorithm. i = 1, 2, ..., n is the variable index.
Step 5: Calculate the objective function f (U, ρ), volume V, regional maximum temperature ζ max , and derivatives using the updated design variables ρ i . Normalizations of the objective function and constraints are implemented to satisfy the MMA algorithm.
Step 6: Check if a fixed number of design iterations is reached.
If the condition is met, then an optimization result is found. Otherwise, repeat Steps (4)-(6) with m ← m + 1.
The optimization problem is solved with a fixed number of design iterations to show the excellent convergence and stability properties of the proposed algorithm.

Continuum shape-based optimization algorithm
The parametrized level set-based topology optimization is also implemented. The conventional level set solved by the upwind difference scheme cannot directly extend to the multiconstrained topology optimization problems. The parametrized level set method uses the RBFs to represent the level set function (Wang & Wang, 2006;Wei et al., 2018). Therefore, the optimal material distribution is realized by updating the coefficients of RBFs during the optimization process. The coefficients of RBFs are solved by the MMA algorithm, which is different from the 88line code for the parametrized level set method (Wei et al., 2018). To suppress the dependence of the initial topological guess, the strategy of new hole nucleation is introduced to control the frequency of generating holes and the volume changes for several optimization iterations. When the coefficients of RBFs are taken as the design variables, the derivatives of the objective function and constraints concerning RBFs' coefficients should be derived for the MMA algorithm. From equation (30), the derivative of J ( ) concerning α i can be written as where α = {α j } j=1,2,...,nl are the coefficients of RBFs, R( · ) is a fixed basis function which is a radius related to the Euclidean distance, { p j } j=1,2,...,nl are the computational centers and n l is the number of nodes related to the level set function, and δ(φ) is the Dirac function (Wang et al., 2003), and ϑ is given as Using equation (35), the derivative of the regional maximum temperature is written as whereǔ is in the form For the level set representation, the nonlinear thermal conductivity of equation (1) is defined as where φ is the level set function. The proposed algorithms are an iterative method, structured as Algorithms 2 and 3.

Algorithm 2: Optimization algorithm with parametrized level set model
Step 1: Initialize and discretize the design domain with finite elements. The material parameters for the nonlinear heat conduction problems in equations (1), (2), and (45) are given. The parameters γ , , u ζ , κ 0 , β, and u β for volume constraint, thermal conductivity, and maximum approximation temperature in equations (11) and (12) are chosen for the optimization problems.
Step 2: The regions for the temperature constraints are fixed and filled with the thermally conductive material. Calculate the basis functions R( · ) with local support domain size δ s = 5 , where is the mesh size of finite element analysis. R( · ) does not change during the optimization process. Initialize the design variable α by interpolating the level set function.
Step 3: Find the temperature field u and the adjoint variableλ by solving equations (3) and (29). The derivative of J ( ) concerning α is then obtained by equation (41). For the maximum approximation temperatureζ max , solve the adjoint equation (34) to get the adjoint variableλ ζ , and calculate the derivative ofζ max concerning α using equation (43). The derivatives of the volume concerning α can be obtained by settingǔ = 1 in equation (44).
Step 4: Generate new holes in the design domain using the hole nucleation (HN) Algorithm 3 and record the current material volume as V nh . Then, initialize the parameters for the MMA algorithm and calculate the level set function according to the material distribution. The design variable α is updated by interpolating the level set function. Thus, the range of design variable is defined as [α − α , α + α ], where α is a predefined constant for an adaptive range of design variables. The range of α is dynamically updated during the optimization process.
Step 5: Solve the topology optimization problems with the volume and regional temperature constraints with initial iteration m = 0. Update the design variable α using the MMA algorithm and the range of design variable with [α − α , α + α ]. Step 6: Update the level set function and compute material volume V in the design domain. If the volume change V = V nh − V satisfies V ≤ ν in the current iteration step, then implement Step 8; otherwise, continue. ν is admissible volume change between two adjacent hole nucleation generations.
Step 7: Generate new holes in the design domain using the HN algorithm and update V nh using the current material volume. Update the RBFs' coefficients α i and the level set function φ.
Step 8: Implement the finite element analysis to obtain the temperature field u and the adjoint variablesλ andλ ζ . Then, calculate J ( ), V,ζ max , and their derivatives concerning the design variables α i . Update the input parameters for the inputs of the MMA algorithm.
Step 9: Check if a fixed number of design iterations is reached.
If the condition is met, then an optimization result is found. Otherwise, repeat Steps (5)-(9) with m ← m + 1.

Algorithm 3: Hole nucleation (HN) generation
Step 1: Using equations (30) and (35) over the entire design domain, normalize the derivatives of J ( ) andζ max to the range [0, 1] named as d 1 and d 2 . Then, the summation of their derivatives with the definition d t = 0.5(d 1 + d 2 ) is taken as the criterion to generate new holes.
Step 2: The bisection method is used to find the threshold value to remove the material elements with a given volume change δV t according to the value of d t .
Step 3: Convert the positive sign "+" to the negative sign "−" of the level set function corresponding to the removed elements in Step 2. The positive sign represents material regions, and the negative sign denotes the void phase. Step 4: Reinitialize the level set function as the signed distance function, which satisfies |∇φ| = 1. The reinitialization is realized by directly computing the distance from the mesh points to the geometric boundary.
Algorithm 2 implements topology optimization through iterations between the hole nucleation generation and shape optimization for the multiconstrained nonlinear optimization problems. This paper proposes a new hole nucleation generation strategy in Algorithm 3, which combines the derivative of the regional temperature constraint into the derivative of the objective function to handle the multiconstrained topology optimization problems. Then, the shape optimization is realized by updating the coefficients of RBFs with the MMA algorithm. The volume changes generated through the hole nucleation generation and shape optimization are strictly controlled to avoid the local minima to some extent that is due to the drastic volume changes between two successive hole nucleation generations. In Algorithm 2, the local support domain size δ s is selected as 5 according to some practical considerations of the MMA solver (Svanberg, 2002). The variables α j should preferably be scaled such that 0.1 ≤ α max j − α min j ≤ 100, for all j.

Example one: initializations with Fig. 4
In this section, several examples are presented to illustrate the temperature-constrained topology optimization. To compare, the optimal results without and with the regional temperature constraints, the initial design domain, and the boundary conditions are shown in Fig. 4. Figure 4a and b are given as a rectangular domain with dimensions 0.4 m × 0.3 m. The design models are discretized with 80 × 60 mesh elements for the nonlinear finite element analysis. The grey regions subjected to the regional temperature constraints are shown in Fig. 4b. The density-based and continuum shape-based optimization algorithms are tested separately to verify the effectiveness of different optimization approaches. Unless stated otherwise, all the units in this paper are the International System of Units (SI) for describing the physical quantities.

Case one: density-based topology optimization
The material density ρ = 1, the thermal conductivity κ 0 = 1, κ min = 0.001, p = 3, the nonlinear thermal conductivity parameters β = −0.1 and u β = 1 are set for equation (36). The parameter = 15 is given to approximate the maximum temperature in equation (10). ζ 0 does not arise in its derivative equation (22), so it can be neglected. As shown in Fig. 4, the boundary conditions on 0 are u 0 = 0 and the thermal loads are q v = 0.5. The convection boundary condition is not included in the following examples, since the convection formulation usually belongs to the design-dependent problems (Chen & Kikuchi, 2001;Joo et al., 2017;Yoon & Koo, 2019). The design-dependent problems require special consideration for topology optimization. In other words, the positions and boundaries of holes are not known beforehand and change during the iterative progress. Figure 5c is the optimal results without the temperature constraint. Figure 5d is the corresponding temperature distribution of Fig. 5c. The optimization model and constraints (9) without ζ max are adopted. The optimization is run for a fixed number of 120 iterations. Figure 6 is the convergence curves of the volume ratio and maximum temperature for the case of Fig. 5. The temperature value 1.5394 in the third row of Table 1 is the maximum temperature in the regions that hold the same positions represented as the constrained regions in Fig. 4b. In Table 1, the backslash means that the temperature constraint is not considered in the optimization model, as shown in Fig. 4a. Figure 7 is the optimization process with the regional temperature constraints. Figure 8 is the convergence curves of the volume ratio and maximum temperature of Fig. 7. The optimization model is shown in Fig. 4b. The parameter for the temperature constraints in equation (9) is set to u ζ = 1.4. The result shown in Fig. 7e is different from Fig. 5c. The optimal shape and topology are different when the regional temperature constraints are imposed on the optimization model. The maximum temperature in the constraint regions is 1.3996 when the convergence result is obtained, as shown in Fig. 7f and Table 1. Therefore, the temperature can be reduced through the temperature-constrained topology optimization by comparing the results in Table 1.
For the initial design of Fig. 7, the ratio of maximum thermal conductivity to minimum thermal conductivity is about 5, as shown in Fig. 9b. The degree of nonlinearity is moderate. For the optimal result, the temperature and the corresponding thermal conductivity with κ 0 = 1 and β = −0.1 are marked as red lines, as shown in Fig. 9a and b. As shown in Fig. 9b, the thermal conductivity changes with the optimization iterations, since the material distribution changes during the optimization process. Thus, the degree of nonlinearity changes from moderate level to slight level. The Newton-Raphson method is adopted to solve the nonlinear heat conduction problems, which can    handle different levels of degree of nonlinearity. For a highly nonlinear problem, the parameters related to the Newton-Raphson method should be reset to satisfy the nonlinear finite element analysis. Thus, the topology optimization model and the sensitivity analysis for the nonlinear heat conduction problems can also be used for different levels of degree of nonlinearity. Figure 10e is the optimal result when the parameter of the temperature constraints in equation (9) is set to u ζ = 1.0. The shapes of the thin structures shown in Fig. 10e are different from Fig. 7e when the parameter u ζ is set to different values. Figures 7 and 10 show some intermediate results of optimization iterations, which also indicates different convergence processes. The maximum temperature in the temperature regions as shown in Fig. 10f is 0.9969 when the convergence result is obtained with temperature constraint u ζ = 1.0. Figure 11 is the convergence curves of the volume ratio and maximum temperature of Fig. 10.
The density-based topology optimization for the nonlinear heat conduction problems holds different topological distributions when the temperature constraint is imposed on the optimization model, as shown in Figs 5c, 7e, and 10e. The optimal shapes of material distributions are also different due to different temperature constraints, as shown in Figs 7e and 10e. However, the optimal results are not always available in practice. When the constraint parameter u ζ = 0.5 is adopted, the optimization problem leads to a nonoptimal result, as shown in Fig. A1 and Table A1. After 120 iterations, the result is shown in Fig. A1b when Algorithm 1 is performed. The maximum temperatures in the constraint regions cannot satisfy the temperature constraint value u ζ = 0.5 according to the convergent result 0.6828 in Table A1.

Case two: continuum shape-based topology optimization
The examples of this section use the design domain and boundary conditions, as shown in Fig. 4. The optimization Algorithm 2 with the parametrized level set model is adopted to verify the regional temperature-constrained topology optimization for the nonlinear heat conduction problems. For the parametrized level set-based topology optimization, H = 1.5 with the mesh size and ξ = 1e − 8 are set for the Heaviside function and the Dirac delta function in equations (41), (43), and (45) (Wang et al., 2003). The adaptive interval for the design variable α is set to α = 0.001 in Algorithm 2. ν = 0.06V 0 is the volume change between successive hole nucleation generations. An admissible volume δV t = 0.01V 0 is set for each hole nucleation generation in Algorithm 3. The other parameters are the same as those in Section 5.1.1. Figure 12 is the nonlinear optimization process of the parametrized level set-based topology optimization without temperature constraint. Algorithm 2 is performed to obtain the material distribution. The material elements are directly removed to generate new holes according to the derivative information using Algorithm 3 without the term related to the temperature constraint derivative. The volume change is strictly controlled between two successive implementations of hole nucleation generations, which can be seen from the red circles across the red line in Fig. 13a. The volume change between successive hole nucleation generations is controlled below ν = 0.06V 0 , as shown in Fig. 13b. Figure 14 shows the implementation of Algorithm 3. The reinitialization is then implemented, as shown in Fig. 14c, which satisfies the property of the signed distance function. As shown in Figs 5d and 12f, and the corresponding temperatures in Table 1, the optimal results are similar for the two methods. When the regional temperature constraints with u ς = 1.4 are considered with the optimization model (11), the result and the     temperature distribution are shown in Fig. 15e and f. The material distribution of Fig. 15e is similar to the density-based optimal result, as shown in Fig. 7e. The material elements are removed according to the derivatives related to the combination of J ( ) andζ max in Algorithm 3. The process is marked as circles along the red line in Fig. 16a. From the optimization process of Fig. 15b-d, the derivative information ofζ max affects new hole generation. The removed elements around the temperatureconstrained regions are partially due to the derivative d 2 in Algorithm 3 shown as the left and right boundaries in Fig. 15b and c. The volume changes between successive hole nucleation generations are controlled with the parameter ν from Step 6 in Algorithm 2. When the temperature constraints are considered, the hole nucleation generation is different according to the comparison of Figs 14 and 17. The new holes around the geometric boundaries shown in Fig. 17a are due to the term d 2 in Step 1 of Algorithm 3. The volume changes between successive hole nucleation generations are below ν = 0.06V 0 , as shown in Fig. 16b. The frequency of hole nucleation according to the volume change can suppress the local minimum to some extent. The results in Table 1 show that the temperature-constrained topology optimization for the nonlinear heat conduction problems can reduce the temperature to the constraint condition u ζ = 1.4. However, when the temperature constraint is set to u ζ = 1.0, the nonoptimal result is received, as shown in Fig. A3 and Table A1. From Fig. 17, the derivative of the temperature constraint is combined with the derivative of the objective function, which generates new holes around the temperature-constrained regions. The shape optimization is performed several iterations after new hole nucleation generation, as shown in Fig. 16. The topological generation and shape optimization are iteratively implemented to achieve topology optimization with the parametrized level set method. In appendix, the material distribution of Fig. A3f is like that in Fig. 15, but the final temperature in the constrained regions is 1.3562 that does not satisfy the constraint condition with u ζ = 1.0, as shown in Table A1. Therefore, the optimal results using the density-based and continuum shape optimization algorithms are not always available when the inequality temperature constraint ζ max (ϕ) ≤ u ζ is imposed on the optimization model.
The programs are written in Matlab R2021a Prerelease. All computations are performed on a PC with Intel(R) Core(TM) i7-8850H CPU @ 2.6 GHz processor and 64 GB RAM. For the densitybased optimization algorithm, the entire computational cost mainly includes the nonlinear finite element analysis and the MMA algorithm. For Fig. 5 without temperature constraint, the mean computational times are 2.3237, 2.2773, and 0.0463 s concerning the whole loop, nonlinear finite element analysis, and MMA algorithm of each iteration. For the density-based topology optimization with temperature constraints, the mean computational times are 3.3080, 3.2522, and 0.0557 s with u ζ = 1.4, as shown in Fig. 7. For Fig. 10, the mean computational times are 3.3026, 3.2403, and 0.0622 s with u ζ = 1.0. For the continuum shape-based topology optimization, when the temperature constraint is not considered, as shown in Fig. 12, the mean computational times are 12. 6296, 4.2850, 8.0084, and 3.7340 s concerning the whole loop, nonlinear finite element analysis, MMA algorithm, and hole nucleation generation of each iteration, respectively. For the temperature constraint case of Fig. 15, the mean computational times are 15. 8796, 5.9364, 9.5496, and 5.0778 s, respectively. Therefore, the continuum shape-based topology optimization takes more computational times than the density-based topology optimization. In the above numerical cases, the temperature-constrained problems hold more computational cost than the cases without temperature constraints.

Example two: initializations with Fig. 18
To further illustrate the effectiveness of the proposed methods for temperature reduction, the following examples are initialized as a square domain with dimensions 0.3 m × 0.3 m, as shown in Fig. 18. The design models are discretized with 60 × 60 mesh elements for the nonlinear finite element analysis. The grey regions subjected to the temperature constraints are presented in Fig. 18b. The positions of regional temperature constraints are different from those of Fig. 4b. The density-based and continuum shape-based optimization algorithms are tested separately to verify the effectiveness of different optimization approaches. Other parameters for optimization are the same as those of the above examples.

Case one: density-based topology optimization
When the density-based optimization method is adopted, Fig. 19c is the optimal result without temperature constraint       Fig. 18a, and (f) temperature distribution of optimal result 23e.
for the nonlinear heat conduction problems. Figure 19d is the temperature distribution of the optimal result. The optimization model (9) without ζ max is adopted. Figure 20 is the convergence curves of the volume ratio and maximum temperature for the case of Fig. 19. The maximum temperatures are 1.0563 in the temperature constraint regions, as shown in the third row of Table 2. Compared with Fig. 19d, when the temperature constraints are imposed on the design problem with u ζ = 0.5, the maximum temperature is reduced to 0.5008 in the temperature constraint regions, as shown in Fig. 21f and the fourth row of Table 2. Compared with Fig. 19c, the result of multiconstrained topology optimization is different and presented in Fig. 21e. The heat conduction paths and the material distribution in Fig. 21e are different from Fig. 19c to reduce temperatures in the temperature-constrained regions. Therefore, the material topology and shape are affected by the temperature constraints.   Fig. 18b, and (f) temperature distribution of optimal result 25e. Figure 22 is the convergence curves of the volume ratio and maximum temperature of Fig. 21.

Case two: continuum shape-based topology optimization
The initializations of boundary conditions are shown in Fig. 18. Figure 23 is the optimization process of the parametrized level set-based optimization method. The maximum temperature in the temperature-constrained regions is 1.0227, which is given in the fifth row of Table 2. Figure 24a shows the volume and temperature changes during the optimization process. The volume changes between two successive hole generations are controlled below ν = 0.06V 0 , as shown in Fig. 24b. The optimal result with temperature constraints u ς = 0.5 is presented in Fig. 25e. Compared with Fig. 23f, the temperature is reduced to 0.5007, as shown in Fig. 25f and Table 2. Figure 26 shows the optimization process when the temperature constraints are imposed in the structure. Table 2 gives the various cases of topology optimization for the nonlinear heat conduction problems. The parametrized level set-based topology optimization can achieve optimal results with smooth geometric boundaries. The volume change ν between two adjacent hole nucleation generations is set to 0.06V 0 that affects the convergence process. The volume change is generated by the shape optimization with the MMA algorithm. If this parameter is set too small, the successive hole nucleation generations would occur, which means the volume change of shape optimization exceeds ν . Thus, the computational cost would be expensive due to the nonlinear finite element analysis during the generation of hole nucleation. On the contrary, if the volume change ν is too large, the new hole nucleation positions would be affected, which affects the topological distribution of material. Therefore, the parameter should be more reasonable when the shape optimization is performed at least once between two adjacent hole nucleation generations.
To compare with examples in Section 5.1, the volume change ν between two adjacent hole nucleation generations is also set to 0.06V 0 . As shown in Fig. 26, the successive hole nucleation generations are invoked since the volume change generated by shape optimization between two adjacent iterations is larger than the predefined 0.06V 0 . The optimal results, as shown in Figs 23e and 25e, are similar to the results with the densitybased optimization algorithms, as shown in Figs 19c and 21e. That illustrates that the parameter ν with a relatively small value does not affect the optimal result.
For parameter α , if it is set to a large value, the shape optimization by the parametrized level set method may generate large volume change after one iteration, which is not conducive to achieving convergence result, especially with more thin members and details in Fig. 25e. The volume change after one iteration should be less than ν in Step 6 of Algorithm 2, which means it is more reasonable to perform shape optimization at least once after hole nucleation generation. Hence, the parameter α should be selected according to the volume change after shape optimization.
With the density-based algorithm, the mean computational times related to the whole loop, nonlinear finite element analysis, and MMA algorithm of each iteration are 1.6029, 1.5706, and 0.0322 s and 2.1550, 2.1020, and 0.0529 s for the cases of Figs 19 and 21, respectively. For Figs 23 and 25, the mean computational times related to the whole loop, nonlinear finite element analysis, MMA algorithm, and hole nucleation generation of each iteration are 8.8917, 3.0028, 5.5312, and 2.3846 s and 10.8671, 3.9750, 6.1419, and 3.1603 s, respectively. In conclusion, the com-putational cost increases with the mesh density. The continuum shape-based algorithm takes more time than the density-based algorithm with similar optimization parameters. The nonlinear finite element analysis takes up most of the computational time for the density-based algorithm, such as Figs 19 and 21. For the continuum shape-based algorithm, the global MMA algorithm takes more computational time than the nonlinear finite element analysis corresponding to Figs 23 and 25. For numerical examples, Algorithm 1 is the density-based optimization algorithm for the nonlinear heat conduction problems. Algorithm 2 is the continuous shape-based optimization algorithm. Both optimization algorithms are solved by the MMA solver. For different initializations, as shown in Figs 4 and 18, the parameter settings are the same as possible, which ensures the fairness of comparisons of numerical examples.

Conclusions
This paper presents a practical and computationally efficient procedure for solving the temperature-constrained topology optimization of nonlinear heat conduction problems. The maximum approximation temperature measure for the regional temperature constraints is proposed with the density-based and continuum descriptions. Then, the density-based and continuum shape-based optimization algorithms are performed to verify the effectiveness of the temperature constraint formulations for the topology optimization problems. For the continuum shape-based optimization, the material derivative of regional temperature constraint is derived using the adjoint variable method. The numerical examples demonstrate that the proposed method effectively reduces the temperatures in the constrained regions. However, the temperature constraint formulation for topology optimization cannot always obtain the acceptable results as shown in Appendix A, since the material volume that restricts the final temperature range is fixed as the constraint condition. A possible solution to handle the above problem is that minimizing the material volume subjected to the regional temperature constraints is adopted as the optimization model, which does not consider thermal compliance. The density-based topology optimization is more effective to reduce temperatures than the continuous shape-based topology optimization for the nonlinear heat conduction problems, though the smooth boundaries can be obtained without implementing any filter by the parametrized level set-based topology optimization. When the parametrized level set-based optimization model is adopted to perform the continuum shape optimization, the adaptive bounds of the coefficients of RBFs are proposed to control the volume changes during the optimization process. The hole nucleation generation is implemented by considering the derivative of temperature constraint for the multiconstrained optimization problems, which overcomes the dependences of the initial guess for the level set-based topology optimization. Furthermore, the proposed approach can certainly be applied to other problems of structural optimization of nonlinear elasticity involving multiphysics.

Non-Optimal Results with Regional Temperature Constraints
The optimal results of Appendix A are used to compare with the results presented in Section 5. Figures A1 and A2 and Table A1 are the results using the optimization Algorithm 1. This example shows that the optimization algorithm falls into local minima with the regional temperature constraints. Figures A3 and A4 and Table A1 are the optimal results using the parametrized level set-based optimization Algorithms 2 and 3. The iteration steps are 1, 3, 10, 20, 30, and 600, corresponding to results shown in Fig. A3a-f. This example shows that the parametrized level set-based multiconstrained topology optimization also falls into local minima, which cannot reach the optimal result satisfying the regional temperature constraints.