Topology optimization for thermal structures considering design-dependent convection boundaries based on the bidirectional evolutionary structural optimization method

Based on the bidirectional evolutionary structural optimization (BESO) method, the present article proposes an optimization method for a thermal structure involving design-dependent convective boundaries. Because the BESO method is incapable of keeping track of convection boundaries, virtual elements are introduced to assist in identifying the convection boundaries of the structure. In order to solve the difficult issue of element assignment under a design-dependent convection boundary, label matrixes are employed to modify the heat transfer matrix and the equivalent temperature load vector of elements over topology iterations. Additionally, the optimization objective is set to minimize the maximum temperature of the structure in order to deal with the objective reasonableness, and the p-norm method is then used to fit the objective function to calculate sensitivity. Finally, several cases, including 2D and 3D structures under various heat transfer boundary conditions, are provided to illustrate the effectiveness and good convergence of the proposed method.

1 Introduction

Structural topology optimization, as an effective method of structural weight reduction, has been paid more and more attention in recent years. However, a lightweight design for a thermal structure under various boundary conditions still faces great challenges (Deaton and Grandhi, 2015; Zhu et al., 2015), namely because studies under first-type (Dirichlet boundary) and second-type (Neumann boundary) boundary conditions are more common than those under convection (Robin boundary) boundary conditions (Zhou et al., 2016).

Previous topology optimization of thermal structure considering the Robin (third-type) boundary (i.e., the convection boundary) has mostly been founded on density-based and level set methods. In density-based optimizations, it is necessary to create explicit boundaries by setting a cutoff value for elements' density; some classical manners to do this have been proposed, including the following: the hat function (Iga et al., 2009; Dede et al., 2015), the peak interpolation method (Yin and Ananthasuresh, 2001), the element connection parameter method (Ho Yoon and Young Kim, 2005) and convection interpolation (Bruns, 2007; Alexandersen, 2011). However, the choice of the truncation value of the above methods depends heavily on the experience of the designer, and an inappropriate cutoff value may result in ill-defined boundaries, such as “islands” (i.e., isolated solid materials). Some studies have been based on the explicit bounds of the level set method in order to realize the implementation of convection boundaries (Ahn and Cho, 2010; Coffin and Maute, 2015; Li et al., 2022). However, the level set method has the disadvantage of initialization dependence; thus, the location of convection boundaries will greatly affect the optimization results. More importantly, the above studies did not consider the fact that a structure boundary might be partly convective (Hu et al., 2008; Wang and Qian, 2020; Yan et al., 2021). Hence, it is necessary to identify whether a boundary is convective, adiabatic or partly convective. However, the difficulty here is that the solid or fluid region is updated during the iterations. As a result, the convection boundary changes with the iteratively altered structural boundary, as shown in Fig. 1.

Figure 1 Change in the design-dependent convection boundary with iteration.

Bidirectional evolutionary structural optimization (BESO) methods are also one of the main methods for topology optimization. They have the advantage of good convergence and initialization independence (Radman, 2021) and have been widely used in structural heat transfer (Xia et al., 2018), frequency (Gan and Wang, 2021), acoustics (Pereira et al., 2022) and material microstructure optimization design (Qiao et al., 2019), among others. In addition, BESO-type methods can provide a clear, explicit boundary (Da et al., 2018), which has been employed to solve the transmissible load problem (Fuchs and Moses, 2000; Yang et al., 2005) and the pressure load problem (Picelli et al., 2015, 2017). However, according to the heat transfer equation, the convective boundary is more complex than the transmissible or pressure load, as the convective boundary affects not only the equivalent temperature load but also the heat transfer capacity of the structure itself. Therefore, the current BESO method needs to be further investigated with respect to the design-dependent convection boundary in order to improve its applications to more complicated boundary conditions.

This article is devoted to the topological optimization of thermal structures considering the design-dependent convection boundary based on the BESO method. The remainder of the paper is structured as follows: in Sect. 2, the difficulties involved with thermal topology optimization considering the convection boundary are outlined; in Sect. 3, the corresponding solution method is proposed, including the identification of the convection boundary and the assignment of the convection matrix; in Sect. 4, the selection of the topology optimization objective is discussed, and the corresponding sensitivity is analyzed; and, finally, in Sect. 5, numerical cases with various boundary conditions are analyzed to verify the effectiveness and good convergence of the proposed method.

2 Problem statement

For the heat transfer problem, there are three types of thermal boundary conditions (Palani and Ganesan, 2007): the first-type boundary (S₁), called the Dirichlet condition, corresponds to a given fixed temperature on the system; the second-type boundary (S₂), namely the Neumann condition, corresponds to a given heat flux on the boundary; and the third-type boundary (S₃), the Robin condition, corresponds to a given convection coefficient and the ambient temperature on the boundary. These boundary conditions can be expressed as follows:

$$\begin{array}{}\text{(1)}& \begin{array}{l}{\mathit{q}}_{\mathrm{T}}=\stackrel{\mathrm{‾}}{T}\text{on}{S}_{\mathrm{1}};\\ {\stackrel{\mathrm{‾}}{k}}_x{\displaystyle \frac{\partial {\mathit{q}}_{\mathrm{T}}}{\partial x}}{n}_x+{\stackrel{\mathrm{‾}}{k}}_y{\displaystyle \frac{\partial {\mathit{q}}_{\mathrm{T}}}{\partial y}}{n}_y+{\stackrel{\mathrm{‾}}{k}}_z{\displaystyle \frac{\partial {\mathit{q}}_{\mathrm{T}}}{\partial z}}{n}_z={\stackrel{\mathrm{‾}}{q}}_{\mathrm{f}}\text{on}{S}_{\mathrm{2}};\\ {\stackrel{\mathrm{‾}}{k}}_x{\displaystyle \frac{\partial {\mathit{q}}_{\mathrm{T}}}{\partial x}}{n}_x+{\stackrel{\mathrm{‾}}{k}}_y{\displaystyle \frac{\partial {\mathit{q}}_{\mathrm{T}}}{\partial y}}{n}_y+{\stackrel{\mathrm{‾}}{k}}_z{\displaystyle \frac{\partial {\mathit{q}}_{\mathrm{T}}}{\partial z}}{n}_z={\stackrel{\mathrm{‾}}{h}}_{\mathrm{c}}\left(T_{\mathrm{\infty }}-{\mathit{q}}_{\mathrm{T}}\right)\text{on}{S}_{\mathrm{3}}.\end{array}\end{array}$$

Here, $\stackrel{\mathrm{‾}}{T}$ is temperature given on S₁, q_T is the nodal temperature of the structure and T_∞ is the ambient temperature. ${\stackrel{\mathrm{‾}}{k}}_x$, ${\stackrel{\mathrm{‾}}{k}}_y$ and ${\stackrel{\mathrm{‾}}{k}}_z$ represent the thermal conductivity in the respective x, y and z directions. n_x, n_y and n_z are the cosine of the normal line outside the boundary. ${\stackrel{\mathrm{‾}}{q}}_{\mathrm{f}}$ is the given heat flux on S₂, and ${\stackrel{\mathrm{‾}}{h}}_{\mathrm{c}}$ is the convection coefficient on S₃.

When the boundary conditions are satisfied, the steady-state heat equation of discrete elements is obtained using the finite element method (FEM):

$$\begin{array}{}\text{(2)}& {\mathbf{K}}_{\mathrm{T}}^{\mathrm{e}}{\mathit{q}}_{\mathrm{T}}^{\mathrm{e}}={\mathit{P}}_{\mathrm{T}}^{\mathrm{e}},\end{array}$$

where ${\mathbf{K}}_{\mathrm{T}}^{\mathrm{e}}$ is the element heat transfer matrix, ${\mathit{q}}_{\mathrm{T}}^{\mathrm{e}}$ is the temperature vector of element nodes and ${\mathit{P}}_{\mathrm{T}}^{\mathrm{e}}$ is the equivalent temperature load vector of element nodes. The specific forms of ${\mathbf{K}}_{\mathrm{T}}^{\mathrm{e}}$ and ${\mathit{P}}_{\mathrm{T}}^{\mathrm{e}}$ are as follows:

$$\begin{array}{}\text{(3)}& {\displaystyle }\begin{array}{rl}{\mathbf{K}}_{\mathrm{T}}^{\mathrm{e}}& =\underset{{\mathrm{\Omega }}^{\mathrm{e}}}{\int }[k_x{\left({\displaystyle \frac{\partial \mathit{N}}{\partial x}}\right)}^{\mathrm{T}}\left({\displaystyle \frac{\partial \mathit{N}}{\partial x}}\right)+k_y{\left({\displaystyle \frac{\partial \mathit{N}}{\partial y}}\right)}^{\mathrm{T}}\left({\displaystyle \frac{\partial \mathit{N}}{\partial y}}\right)\\ & +k_z{\left({\displaystyle \frac{\partial \mathit{N}}{\partial z}}\right)}^{\mathrm{T}}\left({\displaystyle \frac{\partial \mathit{N}}{\partial z}}\right)]\mathrm{d}{\mathrm{\Omega }}^{\mathrm{e}}+\underset{S_{{\mathrm{3}}^{\mathrm{e}}}}{\int }{\stackrel{\mathrm{‾}}{h}}_{\mathrm{c}}{\mathit{N}}^{\mathrm{T}}\mathit{N}\mathrm{d}{S}_{{\mathrm{3}}^{\mathrm{e}}};\end{array}\text{(4)}& {\displaystyle }{\mathit{P}}_{\mathrm{T}}^{\mathrm{e}}={\int }_{S_{{\mathrm{2}}^{\mathrm{e}}}}{\stackrel{\mathrm{‾}}{q}}_{\mathrm{f}}{\mathit{N}}^{\mathrm{T}}\mathrm{d}{S}_{{\mathrm{2}}^{\mathrm{e}}}+{\int }_{S_{{\mathrm{3}}^{\mathrm{e}}}}{\stackrel{\mathrm{‾}}{h}}_{\mathrm{c}}{T}_{\mathrm{\infty }}{\mathit{N}}^{\mathrm{T}}\mathrm{d}{S}_{{\mathrm{3}}^{\mathrm{e}}}.\end{array}$$

Here, N is the shape function matrix. It should be noted that the internal heat source intensity is not considered in this article.

According to Eq. (3), the first term on the right of the equation represents the contribution of conductivity, and the second term represents the contribution of convection. According to Eq. (4), the first and the second terms on the right of the equation reflect the equivalent temperature load generated by given heat flow and given convection, respectively. Further, Eqs. (3) and (4) can be expressed as follows:

$$\begin{array}{}\text{(5)}& {\displaystyle }{\mathbf{K}}_{\mathrm{T}}^{\mathrm{e}}={\mathbf{K}}_{\stackrel{\mathrm{‾}}{\mathrm{T}}}^{\mathrm{e}}+{\mathbf{H}}_{\mathrm{h}}^{\mathrm{e}},\text{(6)}& {\displaystyle }{\mathit{P}}_{\mathrm{T}}^{\mathrm{e}}={\mathit{P}}_{\mathrm{f}}^{\mathrm{e}}+{\mathit{P}}_{\mathrm{h}}^{\mathrm{e}}.\end{array}$$

Here,

$$\begin{array}{}\text{(7)}& {\displaystyle }{\mathbf{K}}_{\stackrel{\mathrm{‾}}{\mathrm{T}}}^{\mathrm{e}}=\stackrel{\mathrm{‾}}{k}{{\mathbf{K}}_{\stackrel{\mathrm{‾}}{\mathrm{T}}}^{\mathrm{e}}}^{\mathrm{0}}=\stackrel{\mathrm{‾}}{k}{\int }_{{\mathrm{\Omega }}^{\mathrm{e}}}\mathrm{\nabla }{\mathit{N}}^{\mathrm{T}}\mathrm{\nabla }\mathbf{N}\mathrm{d}{\mathrm{\Omega }}^{\mathrm{e}},\text{(8)}& {\displaystyle }{\mathbf{H}}_{\mathrm{h}}^{\mathrm{e}}={\stackrel{\mathrm{‾}}{h}}_{\mathrm{c}}{{\mathbf{H}}_{\mathrm{h}}^{\mathrm{e}}}_{\mathrm{0}}={\stackrel{\mathrm{‾}}{h}}_{\mathrm{c}}{\int }_{S_{{\mathrm{3}}^{\mathrm{e}}}}{\mathit{N}}^{\mathrm{T}}\mathit{N}\mathrm{d}{S}_{{\mathrm{3}}^{\mathrm{e}}},\text{(9)}& {\displaystyle }{\mathit{P}}_{\mathrm{h}}^{\mathrm{e}}={\stackrel{\mathrm{‾}}{h}}_{\mathrm{c}}{{\mathit{P}}_{\mathrm{h}}^{\mathrm{e}}}^{\mathrm{0}}={\stackrel{\mathrm{‾}}{h}}_{\mathrm{c}}\underset{S_{{\mathrm{3}}^{\mathrm{e}}}}{\int }{T}_{\mathrm{\infty }}{\mathit{N}}^{\mathrm{T}}\mathrm{d}{S}_{{\mathrm{3}}^{\mathrm{e}}}.\end{array}$$

In the above expressions, ${\mathbf{K}}_{\stackrel{\mathrm{‾}}{\mathrm{T}}}^{\mathrm{e}}$ is the elemental heat transfer matrix of thermal conductivity, ${{\mathbf{K}}_{\stackrel{\mathrm{‾}}{\mathrm{T}}}^{\mathrm{e}}}^{\mathrm{0}}$ is the identity matrix form of ${\mathbf{K}}_{\stackrel{\mathrm{‾}}{\mathrm{T}}}^{\mathrm{e}}$, $\stackrel{\mathrm{‾}}{k}$ is the thermal conductivity, ${\mathit{P}}_{\mathrm{T}}^{\mathrm{e}}$ is the element equivalent temperature load, ${\mathit{P}}_{\mathrm{f}}^{\mathrm{e}}$ is the element equivalent temperature load of the heat flux, ${\mathbf{H}}_{\mathrm{h}}^{\mathrm{e}}$ is the element heat transfer matrix of convection, ${{\mathbf{H}}_{\mathrm{h}}^{\mathrm{e}}}^{\mathrm{0}}$ is the identity matrix form of ${\mathbf{H}}_{\mathrm{h}}^{\mathrm{e}}$, ${\mathit{P}}_{\mathrm{h}}^{\mathrm{e}}$ is the element equivalent temperature load of convection and ${{\mathit{P}}_{\mathrm{h}}^{\mathrm{e}}}^{\mathrm{0}}$ is the identity vector form of ${\mathit{P}}_{\mathrm{h}}^{\mathrm{e}}$.

Moreover, the element heat transfer matrix and the equivalent temperature load vector of element nodes are assembled using the FEM to form the heat transfer matrix and the temperature load vector of structure nodes:

$$\begin{array}{}\text{(10)}& {\mathbf{K}}_{\mathrm{T}}{\mathit{q}}_{\mathrm{T}}=\left({\mathbf{K}}_{\stackrel{\mathrm{‾}}{\mathrm{T}}}+{\mathbf{H}}_{\mathrm{T}}\right){\mathit{q}}_{\mathrm{T}}={\mathit{P}}_{\mathrm{T}}={\mathit{P}}_{\mathrm{h}}+{\mathit{P}}_{\mathrm{f}},\end{array}$$

where K_T is the heat transfer matrix, ${\mathbf{K}}_{\stackrel{\mathrm{‾}}{\mathrm{T}}}$ is the heat transfer matrix without convection, H_T is the convection heat transfer matrix, P_T is the equivalent temperature load vector, P_f is the equivalent temperature load generated from the heat flux and P_h the equivalent temperature load generated by convection. When convection is excluded, P_h=0.

In the BESO method, when the value of topological variables is 1, the corresponding elements are solid elements (SEs); when the value is a small quantity, such as 0.001 (not 0 to avoid singularity), the corresponding elements are nonsolid elements with very poor properties (Huang and Xie, 2010). If the design-dependent convection boundary is considered, the boundaries can be classified into two types: convective boundaries and adiabatic boundaries. Therefore, in this article, the nonsolid elements should be further subclassified into fluid elements (FEs) and void elements (VEs). As a result, the convection boundary can be identified from the interfaces between the SEs and the connected FEs (Deshmukh and Warkhedkar, 2013; Chamkha et al., 2017).

However, due to the variation in the FEs, the VEs and the SEs with optimization iteration, it is necessary to take measures to keep track of the iteratively changed convection boundary. Additionally, according to Eq. (10), the assignment of K_T and P_T also needs extra treatments.

3 Convection boundary identification and its assignment

This section is dedicated to realizing the identification of the convection boundary in the design domain during the iteration processes (in Sect. 3.1) and to assigning the corresponding convection contribution matrix (in Sect. 3.2).

3.1 Identification of the convection boundaries

In 2D structures, convection can be classified into two categories: (1) top and bottom convection (TBC), in which the direction is vertical to the structure surface, and (2) side convection (SC), in which the direction is parallel to the structure surface (Alexandersen, 2011). In the TBC case, all sides of an element will be influenced by convection, and the structural boundaries in such problems are all convective. In the SC case, convection affects partial sides of an element, and it is necessary to further confirm whether a boundary of the element is convective. Compared with the optimal design under TBC, the optimal design under BC requires an extra scheme to identify convective boundaries (rather than merely based on the geometry boundary); therefore, the latter design it is more complex and will be investigated in this article. In 3D structures, there is only one type of convection; it is similar to SC in 2D structures, but the difference is that the convection affects partial planes of an element. Hence, the convection boundary identification method for SC proposed later in this article can be applied to both 2D and 3D structures.

As the BESO method cannot keep track of convective boundaries (Qiao et al., 2019), virtual elements are introduced in this article in order to distinguish FEs from nonsolid elements. Using convection acting on the surface of a 2D structure as an example, a schematic diagram is shown in Fig. 2 that outlines the following specific steps:

• i.

  In the design domain, the identifier of the solid elements is labeled “1” and that of the nonsolid elements is labeled “2”.

• ii.

  A layer of virtual elements around the design domain are added, and their identifier is marked as “0”. It is noted that these virtual elements do not participate in FEM calculations.

• iii.

  All elements are then “traversed”. If the identifier of an element is 2 and the identifier of any adjacent element is simultaneously 0, the identifier is changed to 0; if this is not the case, the identifier remains at its original value. To ensure identification, it is necessary to carry out this process enough times. The number of times that this process is carried out is related to the number of elements with an initial identifier of 2. As there is the possibility that only the last element of the traversal is marked as 0 and that the identifiers of all the preceding elements are 2, we can only guarantee the validity of the element identifier conversion by ensuring that the minimum number of iterations is the same as the number of initial elements marked as 2.

• iv.

  All of the virtual elements are deleted. The remaining elements with an identifier value of 2, 1 and 0 are void elements, fluid elements and solid elements, respectively.

In each topology iteration, steps (i)–(iv) are repeated to distinguish the FEs and VEs. The corresponding physical convection boundaries then exist at the interfaces between SEs and FEs.

Figure 2 Identification of the convection boundaries.

3.2 Assignment of the convection matrixes

After the identification of the convection boundary of solid elements, it is necessary to further modify the element heat transfer matrix and temperature load vector influenced by convection.

Because convection only affects part of the boundary of the element, this paper introduces label matrixes to achieve the assignment of the matrixes related to convection. There are two parameters that affect convection: the convection coefficient and the ambient temperature. The label matrix of the convection coefficient is denoted as the G matrix. The label matrix of ambient temperature is denoted as the M matrix. In this article, the elements used in 2D structure are quadrilateral elements, and the elements used in 3D structure are hexahedral elements. Therefore, the number of elements of G matrix and M matrix in 2D and 3D structures is 4 and 6, respectively:

$$\begin{array}{}\text{(11)}& \begin{array}{ccl}\mathit{G}& =& {\left[\begin{array}{cccc}{G}_{\mathrm{1}}& G_{\mathrm{2}}& \mathrm{\cdots }& G_n\end{array}\right]}^{\mathrm{T}};\\ \mathit{M}& =& {\left[\begin{array}{cccc}{M}_{\mathrm{1}}& M_{\mathrm{2}}& \mathrm{\cdots }& M_n\end{array}\right]}^{\mathrm{T}}.\end{array}\end{array}$$

Here, n=4 in a 2D structure and n=6 in a 3D structure.

If a component (such as G₁) of matrix G or M is marked as 1, the convection and the ambient temperature are taken into account for an element boundary. On the contrary, 0 means that the convection and the ambient temperature are not considered for an element boundary. Accordingly, employing the G and M matrixes, Eqs. (8) and (9) can be rewritten as follows:

$$\begin{array}{}\text{(12)}& {\displaystyle }\begin{array}{rl}& {\mathbf{H}}_{\mathrm{h}}^{\mathrm{e}}={\stackrel{\mathrm{‾}}{h}}_{\mathrm{c}}{\mathbf{H}}_{\mathrm{h}}^{{\mathrm{e}}_{\mathrm{0}}}={\stackrel{\mathrm{‾}}{h}}_{\mathrm{c}}{\sum }_{i=\mathrm{1}}^n{G}_i{\int }_{{\mathrm{\Gamma }}_i}{\mathit{N}}^{\mathrm{T}}\mathrm{d}{\mathrm{\Gamma }}_i\\ & \left(n=\mathrm{4}\text{in}\text{2D},n=\mathrm{6}\text{in}\text{3D}\right);\end{array}\text{(13)}& {\displaystyle }\begin{array}{rl}& {\mathit{P}}_{\mathrm{h}}^{\mathrm{e}}={\stackrel{\mathrm{‾}}{h}}_{\mathrm{c}}{\mathit{P}}_{\mathrm{h}}^{{\mathrm{e}}_{\mathrm{0}}}={\stackrel{\mathrm{‾}}{h}}_{\mathrm{c}}{\sum }_{i=\mathrm{1}}^n{T}_{\mathrm{\infty }}{G}_i{M}_i{\int }_{{\mathrm{\Gamma }}_i}{\mathit{N}}^{\mathrm{T}}\mathrm{d}{\mathrm{\Gamma }}_i\\ & \left(n=\mathrm{4}\text{in}\text{2D},n=\mathrm{6}\text{in}\text{3D}\right).\end{array}\end{array}$$

Here, Γ is the element boundary.

For the reader's understanding, using the convection acting on the surface of a 2D structure as an example, the steps related to the assignment of the convection matrix to elements by label matrixes are summarized as follows:

• i.

  If an element is nonsolid element (its identifier is 0 or 2), the G and M matrix of this element is a 0 matrix. If it is a solid element, go to the next step.

• ii.

  Check the identifier of its neighboring elements counterclockwise from the left edge, and then fill the G and M matrixes with 0 or 1 in turn.

• iii.

  The element heat transfer matrix and the element temperature load vector considering convection are obtained by Eqs. (12) and (13).

Due to the design dependence of convective boundaries, steps (i)–(iii) need to be repeated in each iteration.

4 Implementation of topology optimization

4.1 Selection of the optimization objective

In the BESO method, the objective function and constraints of optimization are expressed as follows:

$$\begin{array}{}\text{(14)}& \begin{array}{l}\text{Find:}x=\left(x_{\mathrm{1}},x_{\mathrm{2}},\mathrm{\dots }{x}_{\mathrm{e}}\right)\\ \text{Minimize:}\text{Ob}\\ \text{Subject to:}{q}_{\mathrm{T}}^{\mathrm{T}}{\mathbf{K}}_{\mathrm{T}}={\mathit{P}}_{\mathrm{h}}+{\mathit{P}}_{\mathrm{f}}\\ V^{\ast }-\sum _{i=\mathrm{1}}^N{V}_{\mathrm{e}}{x}_{\mathrm{e}}=\mathrm{0}\\ x_{\mathrm{e}}=(\mathrm{0.001},\mathrm{1}).\end{array}\end{array}$$

Here, Ob is the optimization objective, V_e is the volume of each element, x is the topological variable and x_i is the topological variable of element. For more details on the framework of the BESO method, the reader is referred to Huang and Xie (2008).

The optimization objective in this article is exclusively involved with the thermal field. It is initially assumed that the temperature field and load are related to topological variables; thus, the objective function depends on the temperature field, temperature load and topological variables, namely

$$\begin{array}{}\text{(15)}& \begin{array}{l}{\mathit{q}}_{\mathrm{T}}={\mathit{q}}_{\mathrm{T}}\left(\mathit{x}\right)\\ {\mathit{P}}_{\mathrm{T}}={\mathit{P}}_{\mathrm{T}}\left(\mathit{x}\right)\\ \text{Ob}=\text{Ob}\left({\mathit{q}}_{\mathrm{T}}\right(\mathit{x}),{\mathit{P}}_{\mathrm{T}}(x\left)\right).\end{array}\end{array}$$

As the gradient of the node temperature with respect to the element topology variable is solved as an implicit function, a Lagrangian multiplier is introduced in order to eliminate this term (Deb and Srivastava, 2012; Wah et al., 2000):

$$\begin{array}{}\text{(16)}& min\mathit{L}=\text{Ob}+{\mathit{\lambda }}^{\mathrm{T}}\left({\mathit{q}}_{\mathrm{T}}^{\mathrm{T}}{\mathbf{K}}_{\mathrm{T}}-{\mathit{P}}_{\mathrm{T}}\right).\end{array}$$

For this problem, the gradient of the objective function can be obtained:

$$\begin{array}{}\text{(17)}& \begin{array}{rl}\mathrm{\nabla }\mathit{L}& =\left({\displaystyle \frac{\partial \text{Ob}}{\partial {\mathit{q}}_{\mathrm{T}}}}+{\mathit{\lambda }}^{\mathrm{T}}{\mathbf{K}}_{\mathrm{T}}\right){\displaystyle \frac{\partial {\mathit{q}}_{\mathrm{T}}}{\partial \mathit{x}}}+\left({\displaystyle \frac{\partial \text{Ob}}{\partial {\mathit{P}}_{\mathrm{T}}}}-{\mathit{\lambda }}^{\mathrm{T}}\right){\displaystyle \frac{\partial {\mathit{P}}_{\mathrm{T}}}{\partial \mathit{x}}}\\ & +{\mathit{\lambda }}^{\mathbf{T}}{\displaystyle \frac{\partial {\mathbf{K}}_{\mathrm{T}}}{\partial \mathit{x}}}{\mathit{q}}_{\mathrm{T}}.\end{array}\end{array}$$

In previous works, the optimization goal has usually been set to minimize the thermal compliance of the structure (Picelli et al., 2017; Zhou and Li, 2008) – that is, Ob=P_Tq_T. This is because, under a fixed heat flow, minimizing the objective function is equal to maximizing the structural heat transfer matrix, namely

$$\begin{array}{}\text{(18)}& min{\mathit{P}}_{\mathrm{T}}{\mathit{q}}_{\mathrm{T}}=min{\mathit{P}}_{\mathrm{T}}^{\mathrm{T}}{\mathit{P}}_{\mathrm{T}}/{\mathbf{K}}_{\mathrm{T}}\to min{\mathit{P}}_{\mathrm{T}}{\mathit{q}}_{\mathrm{T}}=max{\mathbf{K}}_{\mathrm{T}}.\end{array}$$

However, according to Eqs. (12) and (13), the existence of convection boundaries means that the structure has both a design-dependent heat transfer matrix and nonconstant external load. According to Eq. (18), minimizing the thermal compliance cannot guarantee that the maximum thermal conduction capability is obtained. Thus, the effectiveness of the objective function in Eq. (18) is open to discussion under a design-dependent convection boundary.

In the present article, the optimization objective is set to minimize the maximum temperature in the structural temperature field, which can directly reflect the thermal conduction capability of structures. However, the maximum temperature is a scalar with no gradient. Thus, the p-norm method (Zhai et al., 2018) with a gradient is used to fit the maximum temperature. The objective function can be expressed as follows:

$$\begin{array}{}\text{(19)}& \text{Ob}=\parallel {\mathit{q}}_{\mathrm{T}}{\parallel }_p\to \text{Ob}=({\mathit{L}}_{\mathrm{T}}^{\mathrm{T}}{\mathit{q}}_{\mathrm{T}}^p{)}^{\mathrm{1}/p},\end{array}$$

where L_T is a unit vector used to calculate the summation and p is the value of the norm. In order to illustrate the advantages of minimizing the maximum temperature as an objective function in topology optimization considering design-dependent convection boundaries, a comparison of the topological results with the maximum temperature and the thermal compliance as the respective objective functions is made in Sect. 5.

4.2 Sensitivity analysis

For the BESO method, sensitivity determines the “birth and death” of elements (Xu et al., 2020), which can be solved by the gradient of the objective function (Huang and Xie, 2009). From Eq. (17), the corresponding Lagrangian adjoint matrix of the objective function is

$$\begin{array}{}\text{(20)}& {\mathit{\lambda }}^{\mathrm{T}}=-\left({\mathit{q}}_{\mathrm{T}}^{\mathrm{1}/p}{)}^{\mathrm{T}}\right({\mathbf{L}}_{\mathrm{T}}^{\mathrm{T}}{\mathit{q}}_{\mathrm{T}}{)}^{-(p-\mathrm{1})/p}/{\mathbf{K}}_{\mathrm{T}}.\end{array}$$

Denoting ${\mathit{q}}_{\mathrm{T}}^{\ast }=-\mathit{\lambda }$, the general equation of sensitivity is

$$\begin{array}{}\text{(21)}& \mathrm{\nabla }\mathit{L}=({\mathit{q}}_{\mathrm{T}}^{\ast }{)}^{\mathrm{T}}{\displaystyle \frac{\partial P_{\mathrm{T}}}{\partial \mathbf{x}}}-({\mathit{q}}_{\mathbf{T}}^{\ast }{)}^{\mathrm{T}}{\displaystyle \frac{\partial {\mathbf{K}}_{\mathrm{T}}}{\partial \mathit{x}}}{\mathit{q}}_{\mathrm{T}}.\end{array}$$

The material interpolation model adopted in this article is the Solid Isotropic Material with Penalization (SIMP) model (Tavakoli and Davami, 2008). Then, compared with the sensitivity without considering design-dependent convection $({\mathit{\alpha }}_{\mathrm{e}}=x_{\mathrm{e}}^{(p_{\mathrm{e}}-\mathrm{1})}\stackrel{\mathrm{‾}}{k}\left|\left({q_{\mathrm{T}}^{\mathrm{e}}}^{\ast }{)}^{\mathrm{T}}{{\mathbf{K}}_{\stackrel{\mathrm{‾}}{\mathrm{T}}}^{\mathrm{e}}}^{\mathrm{0}}\right(q_{\mathrm{T}}^{\mathrm{e}}{)}^{\mathrm{T}}\right|)$ and the sensitivity of thermal compliance (${\mathit{\alpha }}_{\mathrm{e}}=x_{\mathrm{e}}^{(p_{\mathrm{e}}-\mathrm{1})}\stackrel{\mathrm{‾}}{k}\left|\left(q_{\mathrm{T}}^{\mathrm{e}}{)}^{\mathrm{T}}{{\mathbf{K}}_{\stackrel{\mathrm{‾}}{\mathrm{T}}}^{\mathrm{e}}}^{\mathrm{0}}\right(q_{\mathrm{T}}^{\mathrm{e}}{)}^{\mathrm{T}}\right|$), the sensitivity considering design-dependent convection can be obtained as follows:

$$\begin{array}{}\text{(22)}& \begin{array}{rl}{\mathit{\alpha }}_{\mathrm{e}}=x_{\mathrm{e}}^{p_{\mathrm{e}}-\mathrm{1}}|& \stackrel{\mathrm{‾}}{k}\left({{\mathit{q}}_{\mathrm{T}}^{\mathrm{e}}}^{\ast }{)}^{\mathrm{T}}{{\mathbf{K}}_{\stackrel{\mathrm{‾}}{\mathrm{T}}}^{\mathrm{e}}}^{\mathrm{0}}\right({\mathit{q}}_{\mathrm{T}}^{\mathrm{e}}{)}^{\mathrm{T}}+{\stackrel{\mathrm{‾}}{h}}_{\mathrm{c}}\left({{\mathit{q}}_{\mathrm{T}}^{\mathrm{e}}}^{\ast }{)}^{\mathrm{T}}{{\mathbf{H}}_{\mathrm{h}}^{\mathrm{e}}}^{\mathrm{0}}\right({\mathit{q}}_{\mathrm{T}}^{\mathrm{e}}{)}^{\mathrm{T}}\\ & -{\stackrel{\mathrm{‾}}{h}}_{\mathrm{c}}\left({{\mathit{q}}_{\mathrm{T}}^{\mathrm{e}}}^{\ast }{)}^{\mathrm{T}}{{\mathit{P}}_{\mathrm{h}}^{\mathrm{e}}}^{\mathrm{0}}\right|,\end{array}\end{array}$$

where p_e is the power law of the penalty. In Eq. (22), the right side of the equation consists of three parts that represent the sensitivity of the heat transfer coefficient, the convection matrix and the equivalent temperature load vector to the objective function, respectively. The detailed derivation process is presented in Appendix A.

5 Case analysis and discussion

Based on the work in Sects. 3 and 4, some typical cases of thermal structures with convective boundary conditions are analyzed in this section, and the topological optimization flowchart is shown in Fig. 3.

Figure 3 The optimal flowchart for thermal structures with a convective boundary.

The convection boundary condition, also known as the mixed boundary, can be combined with both the first-type boundary and the second-type boundary to solve the temperature field. In order to illustrate all of the types of structural optimization designs under convection conditions, topological optimization problems considering combination boundary conditions including design-dependent convection boundaries are classified into three categories:

1. structural optimization under a combination of convection boundaries and the first-type boundary conditions, i.e., the S₁₃ problem;

2. structural optimization under a combination of convection boundaries and the second-type boundary conditions, i.e., the S₂₃ problem;

3. structural optimization under a combination of convection boundary, first-type boundary and second-type boundary conditions, i.e., the S₁₂₃ problem. The S₁₂₃ problem is a hybrid of the S₁₃ problem and the S₂₃ problem.

In addition, it should be noted that convection, including forced convection and natural convection, is a complex boundary condition. The convection coefficient is influenced by the velocity and properties of the fluid as well as by the shape and properties of solid structures (Yang et al., 2013). Although the convection coefficient maybe varies, a constant convection coefficient is specified during topology iterations, as this article mainly focuses on the method of topology optimization for a thermal structure. The materials properties and topology parameters are listed in Table 1.

Table 1 Parameters for the numerical calculations.

Note that “r_min” denotes the filter radius and “ER” represents the evolutionary rate. Detailed information on r_min and ER are given in Huang and Xie (2010).

5.1 The S₁₃ problem

In this subsection, the thermal structure optimization results under a combination boundary conditions including convection and the temperature constraint are discussed. The design domain is a 0.5 m × 0.5 m flat plate with a thickness of 1 cm, and the number of elements is 50×50. The first boundary condition is applied to the geometric center of the design domain. The convective fluid elements are distributed around the design domain, and the newly generated cavities within the structure during topology optimization have no fluids and are adiabatic boundaries. The initial conditions are shown in Fig. 4.

Figure 4 Initial conditions for the S₁₃ problem.

In this case, as the ambient temperature (100 ^∘C) is higher than the structure temperature (0 ^∘C), the convection boundaries mainly provide the equivalent temperature loads to the structure. With the optimization iterations, the optimal processes based on the present proposed method with the maximum temperature as the objective function and the thermal compliance as the objective function are presented in Fig. 5a and b, respectively.

Figure 5 Optimal process for the S₁₃ problem for different optimization objectives: (a) optimization objective of maximum temperature and (b) optimization objective of thermal compliance.

From Fig. 5a and b, with the iteration proceeding, it can be seen that the value of the objective function increases, which is a typical characteristic of BESO (Huang and Xie, 2008). The final maximum temperature and thermal compliance reached based on the proposed method with the maximum temperature as the objective function are 77.44 ^∘C and 14.83 kJ, respectively, as shown in Fig. 5a. Comparably, the final maximum temperature and thermal compliance reached with thermal compliance as the objective function are 80.12 ^∘C and 14.17 kJ, respectively, as shown in Fig. 5b. These results obtained with the two optimization objectives are very close, which is due to the fact that the final boundaries of convection in this case are the same as the initial boundaries, thereby satisfying the condition of a constant convection boundary in Eq. (18), and taking the thermal compliance as the optimization objective is still valid. However, using the maximum temperature as the optimization objective is still meaningful because the maximum temperature of the optimized configuration for this objective is lower.

Additional numerical instability appeared in iteration 5–7 in Fig. 5a, which was caused by the increase in the convection area. As seen in Fig. 6, in iteration 5, the generated topological cavities are disconnected with the outer surface of the structure, and convection only acts on the same part on the surface as that under the initial conditions. However, in iteration 6, the inner cavities propagate to connect with the outer surface, resulting in the convection areas extending to the boundary of the cavities. The increase in the convection area means that convection provides more heat to the structure, leading to a rise in the temperature of the structure. For this problem, the proposed optimization method in the present article, based on the BESO strategy, can correct the inappropriate element elimination according to the sensitivity information. In the following iteration 7, the internal cavity boundaries are modified to disconnect from the external surface, and the optimization iterations appear stable. The resulting configuration and corresponding temperature field are shown in Fig. 7.

Figure 6 Connectivity between the internal cavities and the outer surface of the structure in iterations 5 to 7: (a) iteration 5, (b) iteration 6 and (c) iteration 7.

Figure 7 (a) Final topology configuration and (b) the corresponding temperature field of elements (obtained by averaging the temperature of element nodes).

In Fig. 7a, the outer surface (convection areas) of the structure remains unchanged, which is reasonable to reduce the convection areas (as mentioned earlier). A similar phenomenon was also pointed out by Wang and Qian (2020). Furthermore, the boundaries of the internal cavities are adiabatic, and the temperature in the cavities' body is 0 ^∘C, as shown in Fig. 7b. Compared with the S₁₂ problem (convection boundary replaced by thermal flux, as shown in Fig. 8a), the similarity of the topological iteration curves (as shown in Fig. 8b) and the final configuration (as shown in Fig. 8c) illustrates that the essence of the convection boundary under the S₁₃ condition is to provide heat to the structure. Moreover, the convergence of the proposed method (32 iterations) is comparable to that of the traditional BESO method (36 iterations).

Figure 8 Initial conditions and optimization results for the S₁₂ problem: (a) initial condition for the S₁₂ problem; (b) optimal process for the S₁₂ problem; (c) final configuration for the S₁₂ problem.

5.2 The S₂₃ problem

In this subsection, the optimization of the thermal structure under a combination of boundary conditions including the convection boundary and the heat flux is discussed. For the initial structure, there is an equivalent temperature load P in the center and periphery convection boundaries; other conditions are the same as the case in Sect. 5.1, as shown in Fig. 9.

Figure 9 Initial conditions for the S₂₃ problem.

In this case, as the ambient temperature (0 ^∘C) is lower than the structure temperature with a heat source (P=100 W) in the central structure, the convection boundary is to functionally dissipate heat. The optimal iterative processes and the structural evolution using the proposed method with the maximum temperature as the objective function and the thermal compliance as the objective function are shown in Fig. 10a and b, respectively.

Figure 10 Optimal process for the S₂₃ problem for different optimization objectives: (a) optimization objective of maximum temperature and (b) optimization objective of thermal compliance.

From Fig. 10a, when taking the maximum temperature as the objective function, the distance between the convective boundary and the heat source becomes closer and closer with proceeding iterations (as shown in Fig. 11), which results in a reduction in the initial structural maximum temperature (92.62 ^∘C) to the final value of 88.42 ^∘C, and the corresponding thermal compliance decreases from 9.26 to 8.84 kJ.

Figure 11 The change in distance between the convection boundary and the heat source (getting closer and closer) as well as the change in the corresponding structural temperature field (based on the average temperature of element nodes) with proceeding iterations.

However, as shown in Fig. 10b, when taking thermal compliance as the optimization objective, the iterative curve continues to oscillate during the volume fraction reduction phase (iterations 1–23). When the volume fraction reaches the target value (iterations 23–68), the overall variation trends in the maximum temperature and thermal compliance of the configuration both increase with proceeding iterations. Finally, the final maximum temperature converges to 98.35 ^∘C, which is 5.73 ^∘C higher than the initial value. The corresponding final thermal compliance converges to 9.84 kJ, which is a 0.54 kJ increase compared with the initial value. This is due to the fact that the convective boundaries are changing iteratively in this case, which no longer satisfies the condition of constant convective boundaries in Eq. (18), causing the minimum of the thermal compliance not to equal the maximum of the heat transfer capability and, thus, an unreasonable final configuration. In contrast, it is obvious that the optimization objective of the maximum temperature proposed in this paper can still lead to a good topological optimal configuration.

As the thermal conductivity of the structure is isotropic (according to Table 1 with ${\stackrel{\mathrm{‾}}{k}}_x={\stackrel{\mathrm{‾}}{k}}_y=\mathrm{100}$), the thermal resistance is the same in all directions and, thus, the optimization configuration is going to be circular (as shown in Fig. 12a). If the thermal conductivity of the structure is orthotropic (i.e., ${\stackrel{\mathrm{‾}}{k}}_x=\mathrm{100}$, ${\stackrel{\mathrm{‾}}{k}}_y=\mathrm{60}$), the thermal resistance of the direction with high conductivity is smaller. In order to reach the best heat dissipation for the whole structure, more structural elements with high thermal conductivity are reserved. Conversely, less structural elements are retained in the direction of low thermal conductivity. Thus, the final configuration tends to be elliptical, as shown in Fig. 12b.

Figure 12 Comparison of configurations for isotropic and orthotropic thermal conductivity: (a) isotropy with ${\stackrel{\mathrm{‾}}{k}}_x=\mathrm{100}$, ${\stackrel{\mathrm{‾}}{k}}_y=\mathrm{100}$; (b) orthotropy with ${\stackrel{\mathrm{‾}}{k}}_x=\mathrm{100}$, ${\stackrel{\mathrm{‾}}{k}}_y=\mathrm{60}$.

As mentioned in Fig. 10, the design dependence of convection boundaries has an influence on the final optimal configuration. If the design-dependent convection boundary becomes a fixed boundary, with the G matrixes of the internal elements of the structure set as a constant zero matrix during iterations (see Fig. 13a), the iterative curves (see Fig. 13b) and topological configurations (see Fig. 13c) resemble those (see Fig. 13e, f) under fixed temperature constraints (see Fig. 13d), which illustrates that the essence of the convection boundary under the S₂₃ condition is to dissipate heat from the structure. Moreover, compared with the optimal process curves, the convergence of the proposed method (36 iterations) is comparable to that of the traditional BESO method (36 iterations).

Figure 13 Comparison of the results for the fixed S₂₃ and S₁₂ problems: (a) initial conditions for the fixed S₂₃ problem; (b) optimal process for the fixed S₁₂ problem; (c) configuration for the fixed S₂₃ problem; (d) initial conditions for the S₁₂ problem; (e) optimal process for the S₁₂ problem; (f) configuration for the S₁₂ problem.

5.3 The S₁₂₃ problem

In this subsection, the optimization of the thermal structure under a combination of conditions including the convection boundary, temperature constraint and heat flux are discussed, i.e., the S₁₂₃ problem. The design domain is a rectangular 4 m ×2 m ×1 m parallelepiped, and the number of the corresponding elements is $\mathrm{40}\times \mathrm{20}\times \mathrm{10}$. The geometric center of the design domain has an equivalent temperature load P. The left plane of the design domain has the convection boundary, and the four corners of the right plane have temperature constraints; the remainder boundaries are adiabatic boundaries. In order to illustrate two cases, namely that the convection boundary (1) acts as a cooling boundary to dissipate heat or (2) provides heat to the structure, the ambient temperature T_∞ is set to 0 and 200 ^∘C, respectively. The initial conditions and topology results are shown in Fig. 14a.

In Fig. 14, when T_∞=0, the topology process (as shown in Fig. 14b) resembles that under S₂₃ conditions, and the proposed method can decrease the temperature of the structure by lessening the distance between the design-dependent convection boundaries and the heat source (Fig. 14c). When T_∞=200, the topological process resembles that under S₁₃ conditions, as shown in Fig. 14d, and the proposed method can decrease the rise in structural temperature through design-dependent convection boundaries in the remaining original convection areas (Fig. 14e). Therefore, it can be concluded that the topological variations in convection boundaries are related to the relative temperature of the convection to structure for the S₁₂₃ problem.

Figure 14 Initial conditions and optimization results for the S₁₂₃ problem under different ambient temperatures: (a) initial conditions; (b) optimal process under T_∞=0; (c) optimal configuration under T_∞=0; (d) optimal process under T_∞=200; (e) optimal process under T_∞=200.

6 Conclusion

The present study emphasizes a topology optimization method under design-dependent convection conditions. The main conclusions are as follows:

1. A topology optimization method for a thermal structure with design-dependent convection boundaries is proposed by introducing auxiliary elements to identify the design-dependent convection boundary, employing label matrixes to assign convection-concerned matrixes and utilizing a more appropriate optimal objective to address the design-dependent loads and heat transfer properties induced by convection changes.

2. The effectiveness of the proposed topological method is illustrated using cases with complex thermal boundary conditions (the S₁₂, S₁₃ and S₁₂₃ problems). For all of the simulated cases, the proposed method can give reasonable topological configurations, and the corresponding convergence is comparable to that of the traditional BESO method.

3. In the S₁₃ problem, as the design-dependent convection boundaries provide heat to the structure, the proposed method can reduce the temperature increase without expanding the convective boundary areas (e.g., the temperature dropped sharply from 83.19 to 70.31 ^∘C from iteration 6 to iteration 7, as shown in Figs. 5a and 6). In the S₂₃ problem, the design-dependent convection boundaries function to dissipate heat from the structure. The proposed method can decrease the structural temperature by lessening the distance between convection boundaries and the heat source (the second-type boundary). In the S₁₂₃ problem, as a comprehensive combination of three types of boundaries, the convection boundaries' performance with respect to providing loads or dissipating heat depends on the relative temperature of the convection to the structure.

4. In the S₁₃ problem, with constant convection boundaries, although taking the thermal compliance as the optimization objective is valid, a lower maximum temperature of the structure can be reached when taking the maximum temperature as the objective (77.44 ^∘C vs. 80.12 ^∘C with the former objective). In the S₂₃ problem, with varying convection boundaries, from the iterative process and calculation results, it is obvious that the optimization objective of the maximum temperature proposed in this paper can still lead to a good topological optimal configuration with both a lower temperature of 88.42 ^∘C and a thermal compliance of 8.84 kJ, whereas the thermal compliance objective is invalid because of the nonequivalence of the thermal compliance minimum and the heat transfer capability maximum under varying convective boundary conditions, which can result in an obviously higher temperature of 98.35 ^∘C and a thermal compliance of 9.84 kJ, as seen from Fig. 10.

Appendix A: The detailed derivation process of the sensitivity considering design-dependent convection

Without convection, the equivalent temperature load is independent of topological variables, and Eq. (21) can be simplified as follows:

$$\begin{array}{}\text{(A1)}& \mathrm{\nabla }L=-({\mathit{q}}_{\mathrm{T}}^{\ast }{)}^{\mathrm{T}}{\displaystyle \frac{\partial {\mathbf{K}}_{\stackrel{\mathrm{‾}}{\mathrm{T}}}}{\partial x}}{\mathit{q}}_{\mathrm{T}}.\end{array}$$

When considering convective heat transfer, P_T is no longer a fixed value, which is related to the change in the topological variable x. In this circumstance, Eq. (21) can be transformed into

$$\begin{array}{}\text{(A2)}& \mathrm{\nabla }L=({\mathit{q}}_{\mathrm{T}}^{\ast }{)}^{\mathrm{T}}{\displaystyle \frac{\partial {\mathit{P}}_{\mathrm{h}}}{\partial \mathit{x}}}-({\mathit{q}}_{\mathrm{T}}^{\ast }{)}^{\mathrm{T}}{\displaystyle \frac{\partial {\mathbf{K}}_{\mathrm{T}}}{\partial \mathit{x}}}{\mathit{q}}_{\mathrm{T}}.\end{array}$$

The corresponding discrete forms of Eqs. (A1) and (A2) are as follows:

$$\begin{array}{}\text{(A3)}& {\displaystyle }\mathrm{\nabla }L=-{\sum }_{i=\mathrm{1}}^n\left({{\mathit{q}}_{\mathrm{T}}^{\mathrm{e}}}^{\ast }{)}^{\mathrm{T}}{\displaystyle \frac{\partial {\mathbf{K}}_{\stackrel{\mathrm{‾}}{\mathrm{T}}}^{\mathrm{e}}}{\partial x_i}}\right({\mathit{q}}_{\mathrm{T}}^{\mathrm{e}}{)}^{\mathrm{T}};\text{(A4)}& {\displaystyle }\mathrm{\nabla }L={\sum }_{i=\mathrm{1}}^n\left(\right({{\mathit{q}}_{\mathrm{T}}^{\mathrm{e}}}^{\ast }{)}^{\mathrm{T}}{\displaystyle \frac{\partial {\mathit{P}}_{{\stackrel{\mathrm{‾}}{h}}_{\mathrm{c}}}^{\mathrm{e}}}{\partial x_i}}-\left({{\mathit{q}}_{\mathrm{T}}^{\mathrm{e}}}^{\ast }{)}^{\mathrm{T}}{\displaystyle \frac{\partial {\mathbf{K}}_{\mathrm{T}}^{\mathrm{e}}}{\partial x_i}}\right({\mathit{q}}_{\mathrm{T}}^{\mathrm{e}}{)}^{\mathrm{T}}).\end{array}$$

Under the SIMP penalization model, from Eq. (7), the heat transfer matrix without convection of elements can be expressed as follows:

$$\begin{array}{}\text{(A5)}& {\mathbf{K}}_{\stackrel{\mathrm{‾}}{\mathrm{T}}}^{\mathrm{e}}\left(x\right)={\mathit{x}}^{p_{\mathrm{e}}}\stackrel{\mathrm{‾}}{k}{{\mathbf{K}}_{\stackrel{\mathrm{‾}}{\mathrm{T}}}^{\mathrm{e}}}^{\mathrm{0}}.\end{array}$$

The sensitivity α_e without convection is then obtained, i.e.,

$$\begin{array}{}\text{(A6)}& {\mathit{\alpha }}_{\mathrm{e}}=\left|{\displaystyle \frac{\mathrm{\nabla }{\mathit{L}}_{\mathrm{e}}}{p_{\mathrm{e}}}}\right|=x_{\mathrm{e}}^{(\mathrm{pe}-\mathrm{1})}\stackrel{\mathrm{‾}}{k}\left|\left({{\mathit{q}}_{\mathrm{T}}^{\mathrm{e}}}^{\ast }{)}^{\mathrm{T}}{{\mathbf{K}}_{\stackrel{\mathrm{‾}}{\mathrm{T}}}^{\mathrm{e}}}^{\mathrm{0}}\right({\mathit{q}}_{\mathrm{T}}^{\mathrm{e}}{)}^{\mathrm{T}}\right|.\end{array}$$

Combined with the SIMP model, Eqs. (5) and (13) are rewritten as follows:

$$\begin{array}{}\text{(A7)}& {\displaystyle }{\mathbf{K}}_{\mathrm{T}}^{\mathrm{e}}\left(x\right)=x_{\mathrm{e}}^{p_{\mathrm{e}}}\stackrel{\mathrm{‾}}{k}{{\mathbf{K}}_{\stackrel{\mathrm{‾}}{\mathrm{T}}}^{\mathrm{e}}}^{\mathrm{0}}+x_{\mathrm{e}}^{p_{\mathrm{e}}}{\stackrel{\mathrm{‾}}{h}}_{\mathrm{c}}{{\mathbf{H}}_{\mathrm{h}}^{\mathrm{e}}}^{\mathrm{0}};\text{(A8)}& {\displaystyle }{\mathit{P}}_{\mathrm{h}}^{\mathrm{e}}\left(x\right)=x^{p_{\mathrm{e}}}{\stackrel{\mathrm{‾}}{h}}_{\mathrm{c}}{{\mathit{P}}_{\mathrm{h}}^{\mathrm{e}}}^{\mathrm{0}}.\end{array}$$

Taking the derivatives of Eqs. (A7) and (A8), we have

$$\begin{array}{}\text{(A9)}& {\displaystyle \frac{\partial {\mathbf{K}}_{\mathrm{T}}^{\mathrm{e}}}{\partial x}}=p_{\mathrm{e}}{\mathit{x}}^{(\mathrm{pe}-\mathrm{1})}(\stackrel{\mathrm{‾}}{k}{{\mathbf{K}}_{\stackrel{\mathrm{‾}}{\mathrm{T}}}^{\mathrm{e}}}^{\mathrm{0}}+{\stackrel{\mathrm{‾}}{h}}_{\mathrm{c}}{{\mathbf{H}}_{\mathrm{h}}^{\mathrm{e}}}^{\mathrm{0}}),\text{(A10)}& {\displaystyle \frac{\partial {\mathit{P}}_{\mathrm{h}}^{\mathrm{e}}}{\partial x}}=p_{\mathrm{e}}{\mathit{x}}^{(\mathrm{pe}-\mathrm{1})}{\stackrel{\mathrm{‾}}{h}}_{\mathrm{c}}{T}_{\mathrm{\infty }}{{\mathit{P}}_{\mathrm{h}}^{\mathrm{e}}}^{\mathrm{0}}.\end{array}$$

The sensitivity α_i considering convection is then derived as follows:

$$\begin{array}{}\text{(A11)}& \begin{array}{rl}{\mathit{\alpha }}_{\mathrm{e}}=x_{\mathrm{e}}^{\mathrm{pe}-\mathrm{1}}|& \stackrel{\mathrm{‾}}{k}\left({{\mathit{q}}_{\mathrm{T}}^{\mathrm{e}}}^{\ast }{)}^{\mathrm{T}}{{\mathbf{K}}_{\stackrel{\mathrm{‾}}{\mathrm{T}}}^{\mathrm{e}}}^{\mathrm{0}}\right({\mathit{q}}_{\mathrm{T}}^{\mathrm{e}}{)}^{\mathrm{T}}+{\stackrel{\mathrm{‾}}{h}}_{\mathrm{c}}\left({{\mathit{q}}_{\mathrm{T}}^{\mathrm{e}}}^{\ast }{)}^{\mathrm{T}}{{\mathbf{H}}_{\mathrm{h}}^{\mathrm{e}}}^{\mathrm{0}}\right({\mathit{q}}_{\mathrm{T}}^{\mathrm{e}}{)}^{\mathrm{T}}\\ & -{\stackrel{\mathrm{‾}}{h}}_{\mathrm{c}}\left({{\mathit{q}}_{\mathrm{T}}^{\mathrm{e}}}^{\ast }{)}^{\mathrm{T}}{{\mathit{P}}_{\mathrm{h}}^{\mathrm{e}}}^{\mathrm{0}}\right|.\end{array}\end{array}$$