# A topology optimization method based on non-uniform rational basis spline hyper-surfaces for heat conduction problems

Marco Montemurro, Khalil Refai

► To cite this version:

Marco Montemurro, Khalil Refai. A topology optimization method based on non-uniform rational basis spline hyper-surfaces for heat conduction problems. Symmetry, 2021, 13 (5), 10.3390/sym13050888. [https://doi.org/10.3390/sym13050888](https://doi.org/10.3390/sym13050888)

HAL Id: hal-03307581

https://hal.inrae.fr/hal-03307581

Submitted on 29 Jul 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# A Topology Optimization Method Based on Non-Uniform Rational Basis Spline Hyper-Surfaces for Heat Conduction Problems

Marco Montemurro* and Khalil Refai

Arts et Métiers Institute of Technology, Université de Bordeaux, CNRS, INRA, Bordeaux INP, HESAM Université, I2M UMR 5295, F-33405 Talence, France; khalil.refai@ensam.eu

* Correspondence: marco.montemurro@ensam.eu or marco.montemurro@u-bordeaux.fr

**Abstract:** This work deals with heat conduction problems formulation in the framework of a CAD-compatible topology optimization method based on a pseudo-density field as a topology descriptor. In particular, the proposed strategy relies, on the one hand, on the use of CAD-compatible Non-Uniform Rational Basis Spline (NURBS) hyper-surfaces to represent the pseudo-density field and, on the other hand, on the well-known Solid Isotropic Material with Penalization (SIMP) approach. The resulting method is then referred to as NURBS-based SIMP method. In this background, heat conduction problems have been reformulated by taking advantage of the properties of the NURBS entities. The influence of the integer parameters, involved in the definition of the NURBS hyper-surface, on the optimized topology is investigated. Furthermore, symmetry constraints, as well as a manufacturing requirement related to the minimum allowable size, are also integrated into the problem formulation without introducing explicit constraint functions, thanks to the NURBS blending functions properties. Finally, since the topological variable is represented by means of a NURBS entity, the geometrical representation of the boundary of the topology is available at each iteration of the optimization process and its reconstruction becomes a straightforward task. The effectiveness of the NURBS-based SIMP method is shown on 2D and 3D benchmark problems taken from the literature.

**Keywords:** heat conduction; NURBS hyper-surfaces; topology optimization; SIMP method; minimum length scale requirement; additive manufacturing

<figure><img src="image_2.png" /></figure>

check for
updates

**Citation:** Montemurro, M.; Refai, K. A Topology Optimization Method Based on Non-Uniform Rational Basis Spline Hyper-Surfaces for Heat Conduction Problems. *Symmetry* **2021**, *13*, 888. https://doi.org/10.3390/sym13050888

Academic Editors: Mikhail Sheremet
and Sergei D. Odintsov

Received: 17 April 2021

Accepted: 12 May 2021

Published: 17 May 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_3.png" /></figure>

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

## 1. Introduction

The continuous downscaling of semi-conductor electronics in devices like smartphones and laptops, which require increasing power rates that need to be dissipated, calls for major challenges to design dedicated cooling systems [1]. Efficient heat transfer in electronic devices is of paramount importance because it allows operation at higher performance for longer duration [2]. Moreover, the performance of other devices, like heat exchangers, turbine blades, fins, thermoelectric generators, and cooling systems, critically depends on effective heat transfer [3–7]. Improved heat transfer design can reduce the energy consumption of devices during their operative life, thus allowing for a reduction of the operational costs. The use of dedicated optimization strategies in heat transfer problems can be also applied to increase the efficiency of support structures in various additive manufacturing technologies to save material and to allow producing more compact and cost-effective structures [8].

In the last decade, the use of topology optimization (TO) methods has experienced a significant growth in many industrial fields due also to the development of modern additive manufacturing processes. The goal of TO is to determine the optimal distribution of the material within a given design domain under prescribed requirements. Since its introduction in the late 1980s [9], numerous approaches have been proposed in the literature [10,11], as the homogenization method [9,12], Solid Isotropic Material with
Penalization (SIMP) method [10,13], Level Set Method (LSM) [14,15], Moving Morphable Components (MMC) method [16], Evolutionary Structural Optimization (ESO) method [17], and its improved version, i.e., the Bi-directional ESO [18,19] method.

TO has been extensively used in heat transfer problems [20–27]. An exhaustive review on this topic can be found in [28]. Li et al. [20,21] extend the ESO algorithm to shape optimization problems dealing with steady-state heat conduction. Zhuang et al. [22] employed the LSM to deal with TO heat conduction problems under multiple load cases. Yoon [23] proposed a formulation taking into account for forced convective heat transfer. In this background, the heat transfer equation with forced convective heat loss and the Navier-Stokes equation were integrated in the TO problem formulation. As a consequence, four material properties were interpolated in terms of the pseudo-density field: the inverse permeability, the conductivity, the material density, and the specific heat capacity. Iga et al. [27] focused on the influence of design-dependent heat convection coefficients and internal heat generation on the optimized topologies found by means of the homogenization method [29]. The dependency of these coefficients on the design variables was determined through a dedicated surrogate model. Ikonen et al. [24] formulated heat conduction problems in the framework of an interesting TO method based on parametric L-systems and the finite volume method (FVM). Results were compared to those provided by the classical SIMP approach showing equivalent or superior performances. Recently, Yoon and Koo [26] developed a sensitivity analysis for TO of steady-state conductive thermal problems subject to design-dependent thermal loads using density gradients–based boundary detection, whilst Hu et al. [25] presented an interesting application of TO dealing with heat transfer problems in microchannels in order to optimize their performances in terms of both heat dissipation and pressure drop. Bendsøe and Sigmund in [10] discussed the generalization of the well-known 99-line Matlab code to heat conduction problems.

It is noteworthy that the success of the SIMP method is due to its efficiency and compactness [10]. However, two major drawbacks affect such an approach. Firstly, the topological descriptor relies on the mesh of the finite element (FE) model, and it is not possible to obtain an optimized topology compatible with CAD software: accordingly, the results of the TO require a time-consuming CAD reconstruction/reassembly phase in order to integrate the optimized solution within the CAD environment. Secondly, to avoid the well-known checker-board effect, dedicated filtering techniques [10] or projection methods [30,31] must be introduced, which represent a further limitation of this method.

In order to overcome the aforementioned issues, the classical SIMP approach has been recently reformulated in the framework of Non-Uniform Rational Basis Spline hyper-surfaces [32–39]. The resulting method is then referred to as *NURBS-based SIMP method*. As discussed in [32,33], some consequences of outstanding importance result from this approach: (1) the number of design variables is unrelated to the number of elements and a significant reduction of the design variables amount can be obtained with respect to the classical SIMP approach; (2) the optimized topology is unrelated to the quality of the mesh of the FE model; (3) the NURBS formalism allows taking advantage of an implicitly defined filter zone, whose size depends on the NURBS parameters. Moreover, the CAD reconstruction phase is a trivial task [38], and, thanks to the peculiar features of the NURBS entities, it is possible to meet the design requirements on the reassembled geometry.

In this study, only heat conduction problems are considered, wherein the so-called *thermal compliance* is considered as a measure of the overall structural thermal conductivity. Heat conduction problems are formulated in the context of the NURBS-based SIMP method. The contribution of this study is twofold. Firstly, the gradient of thermal responses with respect to the topological variables is derived by exploiting the main properties of NURBS entities. In particular, the formulation benefits from the local support property of the NURBS blending functions [32,33], which establishes an implicit relation among the pseudo-densities of adjacent elements. Secondly, a sensitivity analysis of the optimized topology to the integer parameters of the NURBS entity is carried out. Particularly, these parameters can be set to define a minimum member size requirement without introducing
an explicit constraint into the problem formulation, as discussed in [35]. The effectiveness of the NURBS-based SIMP method is tested on both 2D and 3D test cases taken from the literature.

The paper follows this outline. Section 2 presents the fundamentals of NURBS hyper-surfaces. Section 3 introduces the main concepts at the basis of the NURBS-based SIMP method. The numerical results are presented and discussed in Section 4: the differences obtained in the optimized topologies when using either B-spline or NURBS entities as topology descriptor are highlighted and the influence of the integer parameters involved in the definition of the NURBS entity is also investigated. Finally, meaningful conclusions and prospects are provided in Section 5.

**Notation 1.** Upper-case bold letters are used to indicate tensors and matrices, while lower-case bold letters indicate column vectors.

## 2. Fundamentals of NURBS Hyper-Surfaces

A NURBS hyper-surface is a polynomial-based function, defined over a parametric space (domain), taking values in the NURBS space (co-domain). Therefore, if N is the dimension of the parametric space and M is the dimension of the NURBS space, a NURBS entity is defined as $\mathbf{h}: \mathbb{R}^N \to \mathbb{R}^M$. The mathematical formula of a generic NURBS hyper-surface is

$$ \mathbf{h}(\zeta_1, \dots, \zeta_N) = \sum_{i_1=0}^{n_1} \cdots \sum_{i_N=0}^{n_N} R_{i_1, \dots, i_N}(\zeta_1, \dots, \zeta_N) \mathbf{P}_{i_1, \dots, i_N}, \quad (1) $$

where $n_j (j = 1, \dots, N)$ is the number of control points (CPs) along the $\zeta_j$ parametric direction, $R_{i_1, \dots, i_N}(\zeta_1, \dots, \zeta_N)$ are the piece-wise rational basis functions, which are related to the standard NURBS blending functions $N_{i_k, p_k}(\zeta_k)$, $k = 1, \dots, N$ by means of the relationship

$$ R_{i_1, \dots, i_N}(\zeta_1, \dots, \zeta_N) = \frac{\omega_{i_1, \dots, i_N} \prod_{k=1}^{N} N_{i_k, p_k}(\zeta_k)}{\sum_{j_1=0}^{n_1} \cdots \sum_{j_N=0}^{n_N} [\omega_{j_1, \dots, j_N} \prod_{k=1}^{N} N_{j_k, p_k}(\zeta_k)]} \quad (2) $$

In Equations (1) and (2), $\mathbf{h}(\zeta_1, \dots, \zeta_N)$ is a M-dimension vector-valued rational function, $(\zeta_1, \dots, \zeta_N)$ are scalar dimensionless parameters defined in the interval $[0, 1]$, whilst $\mathbf{P}_{i_1, \dots, i_N}$ are the CPs coordinates. The j-th CP coordinate $(X_{i_1, \dots, i_N}^{(j)})$ is stored in the array $\mathbf{X}^{(j)}$, whose dimensions are $(n_1 + 1) \times \cdots \times (n_N + 1)$. The explicit expression of CPs coordinates in $\mathbb{R}^M$ is:

$$ \mathbf{P}_{i_1, \dots, i_N} = \{X_{i_1, \dots, i_N}^{(1)}, \dots, X_{i_1, \dots, i_N}^{(M)}\}, \quad (3) \\ \mathbf{X}^{(j)} \in \mathbb{R}^{(n_1+1)\times\dots\times(n_N+1)}, \quad j = 1, \dots, M. $$

Curves and surfaces formulae can be easily deduced from Equation (1). The CPs layout is referred to as control polygon for NURBS curves, control net for surfaces and control hyper-net otherwise [32]. The overall number of CPs constituting the hyper-net is:

$$ n_{\text{CP}} := \prod_{i=1}^{N} (n_i + 1). \quad (4) $$

The generic CP does not actually belong to the NURBS entity but it affects its shape by means of its coordinates. A weight $w_{i_1, \dots, i_N}$ is associated to the generic CP. The higher the weight $w_{i_1, \dots, i_N}$, the more the NURBS entity is attracted towards the CP $\mathbf{P}_{i_1, \dots, i_N}$. For each parametric direction $\zeta_k$, $k = 1, \dots, N$, the NURBS blending functions are of degree $p_k$ and can be computed in a recursive way as
$$N_{i_k, 0}(\zeta_k) = \begin{cases} 1, & \text{if } v_{i_k}^{(k)} \le \zeta_k < v_{i_k+1}^{(k)}, \\ 0, & \text{otherwise,} \end{cases} \qquad (5)$$

$$N_{i_k,q}(\zeta_k) = \frac{\zeta_k - v_{i_k}^{(k)}}{v_{i_k+q}^{(k)} - v_{i_k}^{(k)}} N_{i_k,q-1}(\zeta_k) + \frac{v_{i_k+q+1}^{(k)} - \zeta_k}{v_{i_k+q+1}^{(k)} - v_{i_k+1}^{(k)}} N_{i_k+1,q-1}(\zeta_k), \qquad (6)$$
$$q = 1, \dots, p_k,$$

where each blending function is defined on the knot vector

$$\mathbf{v}_{(k)} = \{\underbrace{0, \dots, 0}_{p_k+1}, v_{p_k+1}^{(k)}, \dots, v_{m_k-p_k-1}^{(k)}, \underbrace{1, \dots, 1}_{p_k+1}\}, \qquad (7)$$

whose dimension is $m_k + 1$, with

$$m_k = n_k + p_k + 1. \qquad (8)$$

Each knot vector $\mathbf{v}_{(k)}$ is a non-decreasing sequence of real numbers that can be interpreted as a discrete collection of values of the related dimensionless parameter $\zeta_k$. The NURBS blending functions are characterized by several interesting properties: the interested reader is referred to [40] for a deeper insight into the matter. Here, only the *local support* property is recalled because it is of paramount importance for the NURBS-based SIMP method [32,33]:

$$R_{i_1, \dots, i_N}(\zeta_1, \dots, \zeta_N) \neq 0, \\ \text{if } (\zeta_1, \dots, \zeta_N) \in [v_{i_1}^{(1)}, v_{i_1+p_1+1}^{(1)} \times \dots \times [v_{i_N}^{(N)}, U_{i_N+p_N+1}^{(N)}]. \qquad (9)$$

Equation (9) means that each CP (and the respective weight) affects only a precise zone of the parametric space, which is referred to as local support or influence zone.

### 3. The NURBS-Based SIMP Method

The details of the formulation of the SIMP method in the NURBS hyper-surfaces framework are given in [32,33]. The main features of the NURBS-based SIMP method for TO are briefly recalled here. The mathematical formulation is here provided for 3D steady-state heat conduction problems under the hypothesis that the applied thermal loads and boundary conditions (BCs) do not depend upon the pseudo-density field (of course, the formulation can be generalized by relaxing this hypothesis). Consider the compact space $\mathcal{D} \subset \mathbb{R}^3$ in a Cartesian orthonormal frame $O(x_1, x_2, x_3)$:

$$\mathcal{D} := \{\mathbf{x}^T = \{x_1, x_2, x_3\} \in \mathbb{R}^3 : x_1 \in [0, a_1], x_2 \in [0, a_2], x_3 \in [0, a_3]\}, \qquad (10)$$

where $a_j (j = 1, 2, 3)$ is a reference length defined along $x_j$ axis. The aim of TO is to search for the best distribution of a given “heterogeneous material” satisfying the requirements of the design problem.

Consider the steady-state equation of the FE model in the most general case:

$$\hat{\mathbf{K}}\hat{\mathbf{u}} = \hat{\mathbf{f}}; \quad \hat{\mathbf{u}}, \hat{\mathbf{f}} \in \mathbb{R}^{\hat{\mathbf{N}}_{\text{DOF}}}, \quad \hat{\mathbf{K}} \in \mathbb{R}^{\hat{\mathbf{N}}_{\text{DOF}} \times \hat{\mathbf{N}}_{\text{DOF}}}, \qquad (11)$$

where $\hat{\mathbf{N}}_{\text{DOF}}$ represents the overall number of degrees of freedom (DOFs) before applying the BCs, and $\hat{\mathbf{K}}$ is the non-reduced (singular) conductivity matrix of the FE model, while $\hat{\mathbf{f}}$ and $\hat{\mathbf{u}}$ are the non-reduced vectors of the external generalized nodal thermal forces and temperatures, respectively. Using standard FE notation [41], Equation (11) can be rewritten as:

$$\begin{bmatrix} \mathbf{K} & \mathbf{K}_{\text{BC}} \\ \mathbf{K}_{\text{BC}}^{\text{T}} & \tilde{\mathbf{K}} \end{bmatrix} \begin{pmatrix} \mathbf{u} \\ \mathbf{u}_{\text{BC}} \end{pmatrix} = \begin{pmatrix} \mathbf{f} \\ \mathbf{r} \end{pmatrix}, \qquad (12)$$
with:

$$ 
\begin{aligned}
& \mathbf{u}, \mathbf{f} \in \mathbb{R}^{N_{\mathrm{DOF}}}, \mathbf{u}_{\mathrm{BC}}, \mathbf{r} \in \mathbb{R}^{N_{\mathrm{BC}}}, \mathbf{K} \in \mathbb{R}^{N_{\mathrm{DOF}} \times N_{\mathrm{DOF}}}, \\
& \mathbf{K}_{\mathrm{BC}} \in \mathbb{R}^{N_{\mathrm{DOF}} \times N_{\mathrm{BC}}}, \hat{\mathbf{K}} \in \mathbb{R}^{N_{\mathrm{BC}} \times N_{\mathrm{BC}}}, \hat{N}_{\mathrm{DOF}} = N_{\mathrm{DOF}} + N_{\mathrm{BC}},
\end{aligned}
\quad (13) $$

where $N_{BC}$ represents the number of DOFs where temperature is imposed (Dirichlet's BC), while $N_{\mathrm{DOF}}$ is the number of unknown DOFs. In Equation (12), $\mathbf{u}$ and $\mathbf{u}_{\mathrm{BC}}$ are the unknown and imposed vectors of nodal temperatures, respectively. $\mathbf{f}$ is the vector of generalized external nodal thermal forces, whilst $\mathbf{r}$ is the vector of (unknown) generalized nodal thermal reactions at nodes where BCs are imposed.

$\mathbf{K}$, $\mathbf{K}_{\mathrm{BC}}$, and $\hat{\mathbf{K}}$ are the conductivity matrices of the FE model after applying BCs. Consider, now, the case of zero Dirichlet's BCs and non-zero Neumann's BCs: the thermal compliance of the structure is defined as:

$$ W := \mathbf{f}^T \mathbf{u}. \quad (14) $$

In the SIMP approach, the material domain $\Omega \subseteq \mathcal{D}$ is identified by means of a pseudo-density function $\rho(\mathbf{x}) \in [0, 1]$ for $\mathbf{x} \in \mathcal{D}$: $\rho(\mathbf{x}) = 0$ denotes absence of material, whilst $\rho(\mathbf{x}) = 1$ indicates completely dense bulk material. The density field affects the element conductivity matrix and, accordingly, the global conductivity matrix of the FE model as

$$ 
\begin{aligned}
\mathbf{K} &:= \sum_{e=1}^{N_e} \rho_e^{\alpha} \mathbf{L}_e^T \mathbf{K}_e^0 \mathbf{L}_e = \sum_{e=1}^{N_e} \mathbf{L}_e^T \mathbf{K}_e \mathbf{L}_e, \\
\mathbf{K}_e^0, \mathbf{K}_e &\in \mathbb{R}^{N_{\mathrm{DOF}}^e \times N_{\mathrm{DOF}}^e}, \mathbf{L}_e \in \mathbb{R}^{N_{\mathrm{DOF}}^e \times N_{\mathrm{DOF}}^e},
\end{aligned}
\quad (15) $$

where $\rho_e$ is the fictitious density computed at the centroid of the generic element $e$, whilst $\alpha \ge 1$ is a parameter penalizing the intermediate densities between 0 and 1, in agreement with the classic SIMP approach ($\alpha = 3$ in this study). $N_e$ is the total number of elements and $N_{\mathrm{DOF}}^e$ is the number of DOFs of the generic element. In Equation (15), $\mathbf{K}_e^0$ and $\mathbf{K}_e$ are the non-penalized and the penalized conductivity matrices of element $e$, expressed in the global reference frame of the FE model, whilst $\mathbf{L}_e$ is the connectivity matrix of element $e$.

In the context of the NURBS-based SIMP method, the pseudo-density field for a TO problem of dimension $D$ is represented through a NURBS hyper-surface of dimension $D+1$. Therefore, for a 3D problem a 4D entity is needed and the pseudo-density field reads:

$$ \rho(\zeta_1, \zeta_2, \zeta_3) = \sum_{i_1=0}^{n_1} \sum_{i_2=0}^{n_2} \sum_{i_3=0}^{n_3} R_{i_1,i_2,i_3}(\zeta_1, \zeta_2, \zeta_3) \rho_{i_1,i_2,i_3}. \quad (16) $$

In Equation (16), $n_{\mathrm{CP}} = (n_1+1)(n_2+1)(n_3+1)$ is the total number of CPs, $\rho(\zeta_1, \zeta_2, \zeta_3)$ constitutes the fourth coordinate of the array $\mathbf{h}$ of Equation (1), while $R_{i_1,i_2,i_3}(\zeta_1, \zeta_2, \zeta_3)$ are the NURBS rational basis functions of Equation (2). The dimensionless parameter $\zeta_j$ can be related to the Cartesian coordinates as follows:

$$ \zeta_j = \frac{x_j}{a_j}, \quad j = 1, 2, 3. \quad (17) $$

Among the parameters tuning the shape of the NURBS entity, only the pseudo-density at CPs and the associated weights are identified as design variables and are collected in the vectors $\zeta_1$ and $\zeta_2$, respectively, defined as:

$$ \zeta_1^T := (\rho_{0,0,0}, \dots, \rho_{n_1,n_2,n_3}), \quad \zeta_2^T := (w_{0,0,0}, \dots, w_{n_1,n_2,n_3}). \quad (18) $$

The dimension of these arrays is equal to $n_{\mathrm{CP}}$. Accordingly, in the most general case, the overall number of design variables is $n_{\mathrm{var}} = 2n_{\mathrm{CP}}$. Thus, the classic TO problem of thermal compliance minimization subject to an inequality constraint on the volume can be formulated as:
$$ \min_{\xi_1, \xi_2} \frac{W}{W_{\text{ref}}}, \quad \text{s.t.:} \quad \left\{ \begin{array}{l} \mathbf{K} \mathbf{u} = \mathbf{f}, \\ \frac{V}{V_{\text{ref}}} - \gamma \le 0, \\ \xi_{1k} \in [\rho_{\text{min}}, \rho_{\text{max}}], \quad \xi_{2k} \in [\omega_{\text{min}}, \omega_{\text{max}}], \\ \forall k = 1, \dots, n_{\text{CP}}. \end{array} \right. \qquad (19) $$

In Equation (19), $V_{\text{ref}}$ is a reference volume, $V$ is the volume of the material domain $\Omega$, while $\gamma$ is the fixed volume fraction; $V_e$ is the volume of element $e$, and $\rho_{\text{min}}$ represents the lower bound, imposed to the density field to prevent any singularity for the solution of the equilibrium problem. The objective function is divided by a reference compliance, $W_{\text{ref}}$, to obtain a dimensionless value. The volume of the material domain appearing in Equation (19) is defined as:

$$ V := \sum_{e=1}^{N_e} \rho_e V_e, \qquad (20) $$

where $V_e$ is the volume of element $e$. Moreover, in Equation (19), the linear index $k$ has been introduced for the sake of compactness. The relationship between $k$ and $i_j$, ($j = 1, 2, 3$) is:

$$ k := 1 + i_1 + i_2(n_1 + 1) + i_3(n_1 + 1)(n_2 + 1). \qquad (21) $$

The other parameters involved in the definition of the NURBS entity (i.e., degrees, knot-vector components and number of CPs) are set a-priori at the beginning of the TO analysis and are not optimized.

The computation of the derivatives of both objective and constraint functions with respect to the design variables is needed to solve problem (19) through a deterministic algorithm. This task is achieved by exploiting the local support property of Equation (9). For instance, the general expressions of the derivatives of both the thermal compliance (in the case $\mathbf{u}_{\text{BC}} = 0$) and the volume read

$$ \frac{\partial W}{\partial \tilde{\xi}_{ik}} = -\alpha \sum_{e \in S_k} \frac{w_e}{\rho_e} \frac{\partial \rho_e}{\partial \tilde{\xi}_{ik}}, \quad i = 1, 2, k = 1, \dots, n_{\text{CP}}, \qquad (22) $$

$$ \frac{\partial V}{\partial \tilde{\xi}_{ik}} = \sum_{e \in S_k} V_e \frac{\partial \rho_e}{\partial \tilde{\xi}_{ik}}, \quad i = 1, 2, k = 1, \dots, n_{\text{CP}}, \qquad (23) $$

where $w_e$ is the compliance of the generic element, $S_k$ is the discretized version of the local support of Equation (9), while $\frac{\partial \rho_e}{\partial \tilde{\xi}_{ik}}$ reads

$$ \frac{\partial \rho_e}{\partial \tilde{\xi}_{ik}} = \begin{cases} R_{k'}^e, & \text{if } i = 1, \\ \frac{R_k^e}{\tilde{\xi}_{2k}}(\tilde{\xi}_{1k} - \rho_e), & \text{if } i = 2. \end{cases} \qquad (24) $$

The scalar quantity $R_k^e$, appearing in Equation (24), is the NURBS rational basis function of Equation (2) evaluated at the element centroid. More details on the analytical passages to derive the gradient of each response function are available in [32,33].

## 4. Numerical Results

The effectiveness of the proposed method is illustrated on 2D and 3D benchmark problems taken from the literature. For each case, the pseudo-density field and the optimum topology are shown. The results presented in this section are obtained by means of the code SANTO (SIMP And NURBS for Topology Optimization) developed at the I2M laboratory in Bordeaux [32,33]. SANTO is coded in the Python® environment and exhibits an easily operable code, with a structure adapted to work with any FE code. In this study, the
commercial code ANSYS® is used to build the FE models and assess the responses of the structure, i.e., nodal temperature and thermal compliance.

Moreover, the Globally-Convergent Method of Moving Asymptotes (GC-MMA) algorithm [42] has been used to perform the solution search for the constrained non-linear programming problem (CNLPP) of Equation (19). The parameters governing the behavior of the GC-MMA algorithm are given in Table 1.

<table><caption>**Table 1.** GC-MMA algorithm parameters.</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td><i>move</i></td><td>0.1</td></tr><tr><td><i>albefa</i></td><td>0.1</td></tr><tr><th>Stop Criterion</th><th>Value</th></tr><tr><td>Maximum n. of function evaluations</td><td>100 × <i>n</i><sub>var</sub></td></tr><tr><td>Maximum n. of iterations</td><td>1000</td></tr><tr><td>Tolerance on objective function</td><td>10<sup>-12</sup></td></tr><tr><td>Tolerance on constraints</td><td>10<sup>-12</sup></td></tr><tr><td>Tolerance on input variables change</td><td>10<sup>-12</sup></td></tr><tr><td>Tolerance on Karush-Kuhn-Tucker norm</td><td>10<sup>-6</sup></td></tr></tbody></table>

As far as numerical tests are concerned, the following aspects are considered: (1) the influence of the geometric entity, i.e., B-spline or NURBS, used to describe the pseudo-density field on the optimized topology is studied for 2D and 3D cases; (2) the influence of the integer parameters, i.e., blending functions degree and CPs number, on the optimized topology is investigated only on 2D benchmark problems for the sake of brevity; (3) the effect of the minimum member size requirement on the optimized topology is highlighted for both 2D and 3D cases.

Regarding the design space of the CNLPP of Equation (19), lower and upper bounds of design variables are set as: $\rho_{min} = 10^{-3}$, $\rho_{max} = 1$; $\omega_{min} = 0.5$, $\omega_{max} = 10$. Moreover, the non-trivial knot vectors components in Equation (19) are evenly distributed in the interval [0, 1] for each benchmark problem.

The material thermal conductivity used for all the considered test cases is $k_m = 1 \text{ Wm}^{-1}\text{K}^{-1}$.

### 4.1. 2D Benchmark Problems

Two 2D benchmark problems are considered here: for each problem, the reference volume $V_{ref}$, appearing in the CNLPP formulation of Equation (19), is the overall area of the domain where the structure is embedded. The geometrical parameters and the BCs of each test case are illustrated in Figure 1. Benchmark 1 (BK1-2D), which has been taken from [43], is a plate having the following dimensions: $a_1 = a_2 = 20 \text{ m}$, $t = 1 \text{ m}$. A heating source $s_h = 0.001 \text{ Wm}^{-2}$ is evenly distributed over the whole design domain. As far as BCs are concerned, the temperature is set to zero for nodes located at $x_1 \in [\frac{a_1-a_t}{2}, \frac{a_1+a_t}{2}]$, $x_2 = a_2$ with $a_t = 2 \text{ m}$, while, on the rest of the boundary the adiabatic wall condition ($q = 0$) is imposed. The FE model is made of $N_e = 80 \times 80$ PLANE55 elements (4 nodes with a single DOF per node).

The second benchmark problem (BK2-2D) is characterized by the same geometrical and material properties of BK1-2D, the only difference being in the applied BCs. In this case, a zero temperature is applied on four heat skins located at $(x_1 \in [\frac{a_1-a_t}{2}, \frac{a_1+a_t}{2}], x_2 = ja_2)$ and $(x_1 = ja_1, x_2 \in [\frac{a_2-a_t}{2}, \frac{a_2+a_t}{2}])$, with $j = 0, 1$, whilst a heat flux $q = 0$ is set on the other regions of the boundary. In addition, in this case, the mesh of the FE model is made of $N_e = 80 \times 80$ PLANE55 elements.

For both test cases, the reference volume is $V_{ref} = a_1a_2$ and the volume fraction appearing in Equation (19) has been set as $\gamma = 0.3$. An initial guess characterized by a uniform pseudo-density field $\rho(\zeta_1, \zeta_2) = \gamma$ has been considered for each analysis. The
reference thermal compliance $W_{\text{ref}}$ appearing in Equation (19) corresponds to the thermal
compliance of the initial guess: its value is 5.956 WK for BK1-2D and 1.0389 WK for BK2-2D.

<figure><img src="image_4.png" /><figcaption>Figure 1. Geometry and boundary conditions of benchmark problems (a) BK1-2D and (b) BK2-2D.</figcaption></figure>

4.1.1. BK1-2D: Sensitivity of the Optimized Topology to the B-Spline and NURBS Entities Integer Parameters

An extensive numerical campaign of tests has been performed on BK1-2D: the aim is to study the sensitivity of the optimized topology to the integer parameters involved in the definition of B-spline and NURBS entities. In particular, the CNLPP of Equation (19) is solved by considering the following combinations of blending functions degrees and CPs numbers: (a) $p_j = 2, 3, 4$, ($j = 1, 2$); (b) $n_{CP} = 40 \times 40, 56 \times 56, 68 \times 68$ for both B-spline and NURBS entities. Moreover, in this case, the CNLPP formulation has been enhanced by adding a symmetry constraint with respect to plane $x_1 = \frac{a_1}{2}$.

Results are provided in terms of dimensionless thermal compliance $\frac{W}{W_{\text{ref}}}$ and number of iterations $N_{\text{iter}}$ to achieve convergence for B-spline and NURBS entities in Figures 2 and 3, respectively. For each solution, the requirement on the volume fraction is always satisfied.
<figure><img src="image_5.png" /></figure>

(a) $p_1 = p_2 = 2$, $n_{\text{CP}} = 40 \times 40$, $\frac{W}{W_{\text{ref}}} = 0.0800$, $N_{\text{iter}} = 112$

<figure><img src="image_6.png" /></figure>

(c) $p_1 = p_2 = 2$, $n_{\text{CP}} = 68 \times 68$, $\frac{W}{W_{\text{ref}}} = 0.0685$, $N_{\text{iter}} = 555$

<figure><img src="image_7.png" /></figure>

(e) $p_1 = p_2 = 3$, $n_{\text{CP}} = 56 \times 56$, $\frac{W}{W_{\text{ref}}} = 0.0756$, $N_{\text{iter}} = 415$

<figure><img src="image_8.png" /></figure>

(g) $p_1 = p_2 = 4$, $n_{\text{CP}} = 40 \times 40$, $\frac{W}{W_{\text{ref}}} = 0.0857$, $N_{\text{iter}} = 571$

<figure><img src="image_9.png" /></figure>

(b) $p_1 = p_2 = 2$, $n_{\text{CP}} = 56 \times 56$, $\frac{W}{W_{\text{ref}}} = 0.0726$, $N_{\text{iter}} = 358$

<figure><img src="image_10.png" /></figure>

(d) $p_1 = p_2 = 3$, $n_{\text{CP}} = 40 \times 40$, $\frac{W}{W_{\text{ref}}} = 0.0833$, $N_{\text{iter}} = 309$

<figure><img src="image_11.png" /></figure>

(f) $p_1 = p_2 = 3$, $n_{\text{CP}} = 68 \times 68$, $\frac{W}{W_{\text{ref}}} = 0.0714$, $N_{\text{iter}} = 463$

<figure><img src="image_12.png" /></figure>

(h) $p_1 = p_2 = 4$, $n_{\text{CP}} = 56 \times 56$, $\frac{W}{W_{\text{ref}}} = 0.0779$, $N_{\text{iter}} = 467$

<caption>Figure 2. *Cont.*</caption>
<figure><img src="image_13.png" /><figcaption>**Figure 2.** Benchmark problem BK1-2D: sensitivity of the optimized topology to CP number and basis functions degrees, B-spline solutions of problem (19); the grey-scale bar refers to the pseudo-density field of Equation (16).</figcaption></figure>

<figure><img src="image_14.png" /></figure>

(a) $p_1 = p_2 = 2$, $n_{CP} = 40 \times 40$, $\frac{W}{W_{ref}} = 0.0712$, $N_{iter} = 573$

<figure><img src="image_15.png" /></figure>

(c) $p_1 = p_2 = 2$, $n_{CP} = 68 \times 68$, $\frac{W}{W_{ref}} = 0.0604$, $N_{iter} = 752$

<figure><img src="image_16.png" /></figure>

(e) $p_1 = p_2 = 3$, $n_{CP} = 56 \times 56$, $\frac{W}{W_{ref}} = 0.0677$, $N_{iter} = 707$

<figure><img src="image_17.png" /></figure>

(b) $p_1 = p_2 = 2$, $n_{CP} = 56 \times 56$, $\frac{W}{W_{ref}} = 0.0647$, $N_{iter} = 667$

<figure><img src="image_18.png" /></figure>

(d) $p_1 = p_2 = 3$, $n_{CP} = 40 \times 40$, $\frac{W}{W_{ref}} = 0.0741$, $N_{iter} = 767$

<figure><img src="image_19.png" /></figure>

(f) $p_1 = p_2 = 3$, $n_{CP} = 68 \times 68$, $\frac{W}{W_{ref}} = 0.0637$, $N_{iter} = 701$

<caption>**Figure 3.** *Cont.*</caption>
<figure><img src="image_20.png" /><figcaption>**Figure 3.** Benchmark problem BK1-2D: sensitivity of the optimized topology to CP number and basis functions degrees, NURBS solutions of problem (19); the grey-scale bar refers to the pseudo-density field of Equation (16).</figcaption></figure>

<figure><img src="image_21.png" /></figure>

A synthesis of the results illustrated in Figures 2 and 3 is shown in Figure 4, where the dimensionless thermal compliance is plotted versus the number of CPs and blending functions degrees. Some remarks can be inferred from the analysis of these results.

<figure><img src="image_22.png" /><figcaption>**Figure 4.** Benchmark problem BK1-2D: dimensionless thermal compliance versus CPs number and degrees for B-spline and NURBS solutions</figcaption></figure>
* For B-spline and NURBS solutions, the higher the number of CPs (or the lower the degree) the smaller the objective function value. As explained in [32,33], this is due to the local support size: the higher the CPs number for a given degree (or the lower the degree for a given number of CPs) the smaller the local support size, consequently smaller topological branches appear in optimized topologies. Moreover, the higher the degree (for a given CPs number) the smoother the boundary of the final topology.

* The NURBS local support can be associated to the concept of the filter zone in standard density-based TO algorithms, as stated in Section 3. According to the definition of the local support of Equation (9), the higher the degree (or the smaller the CPs number) the wider the local support; thus, a single control point affects a wider region of the computation domain. Indeed, as discussed in [35], the local support of Equation (9) enforces a minimum length scale in the optimized topology. Consequently, it can be stated that a high number of CPs and small degrees should be considered if minimum member size does not constitute a restriction for the problem at hand. High degrees and/or small CPs number should be considered otherwise.

* The effect of including the weights among the design variables is twofold: on the one hand, weights contribute to improve the final performances (the objective function of a NURBS solution is always lower than the one of a B-spline solution), whilst, on the other hand, they allow for obtaining optimized topologies characterized by a boundary smoother than the B-spline counterpart.

* The constraint on the volume fraction gets a very small negative value (between $-1 \times 10^{-6}$ and 0) for the optimized topologies resulting from problem (19); thus, the local minimizer is located on the boundary between feasible and infeasible regions.

* All the analyses were performed on a work-station with an Intel Xeon E5-2697v2 processor (2.70–3.50 GHz, Santa Clara, CA, USA) and four cores dedicated to the optimization calculations. The highest computational time occurs for the NURBS solution illustrated in Figure 3 (i), which required about 1.5 h to find the local feasible minimizer.

### 4.1.2. BK2-2D: Minimum Member Size Effect on the Optimized Topology

The influence of the minimum member size requirement is investigated on BK2-2D. In particular, the CNLPP formulation of Equation (19) has been enhanced by considering a constraint on the minimum length scale requirement: the minimum dimension of the optimized topology should be greater than or equal to $d_{min} = 0.3$ m. To automatically satisfy the minimum length scale requirement without introducing an explicit constraint in the problem formulation, according to the methodology presented in [35], B-spline and NURBS entities with $p_1 = p_2 = 3$ and $n_{CP} = 68 \times 68$ CPs are used for this analysis. Moreover, the topology is constrained to be symmetric with respect to planes $x_1 = a_1/2$ and $x_2 = a_2/2$.

The optimized topologies solution of problem (19) are shown in Figure 5, for both B-spline and NURBS entities: numerical results are provided in terms of dimensionless thermal compliance $\frac{W}{W_{ref}}$, number of iterations $N_{iter}$ to achieve convergence and measured minimum member size $d_{min}^m$, i.e., the minimum length scale measured after CAD reconstruction of the boundary of the optimized topology [35]. For each solution, the requirement on the volume fraction is always satisfied (this constraint is always almost active, in the sense that it takes a very small negative value between $-1 \times 10^{-6}$ and 0).

As expected, the NURBS solution is characterized by a dimensionless thermal compliance lower than the B-spline counterpart. Moreover, thanks to the geometrical properties of the NURBS blending functions, the requirement on the minimum length scale is always satisfied. In particular, it is noteworthy that the minimum length scale constraint is almost active for the NURBS solution illustrated in Figure 5.
<figure><img src="image_23.png" /><figcaption>**Figure 5.** Benchmark problem BK2-2D: (a) B-spline and (b) NURBS solutions of problem (19) with $p_1 = p_2 = 3$ and $n_{CP} = 68 \times 68$; the grey-scale bar refers to the pseudo-density field of Equation (16).</figcaption></figure>

## 4.2. A 3D Benchmark Problem

The third benchmark (BK1-3D) deals with the TO of the 3D cubic domain illustrated in Figure 6; the geometry of BK1-3D is characterized by the following dimensions: $a_1 = a_2 = a_3 = 20$ m, $a_t = 2$ m. The material properties are the same as those used in test cases BK1-2D and BK2-2D.

<figure><img src="image_24.png" /><figcaption>**Figure 6.** Geometry and boundary conditions of benchmark problem BK1-3D.</figcaption></figure>

The FE model is made of $N_e = 64,000$ SOLID279 elements (eight nodes, one DOF per node). A null temperature is imposed on the nodes belonging to the heat skin located at $(x_1 \in [\frac{a_1-a_t}{2}, \frac{a_1+a_t}{2}], x_2 \in [\frac{a_2-a_t}{2}, \frac{a_2+a_t}{2}], x_3 = a_3)$, while a zero heat flux is imposed on the rest of the boundary. A heating source $s_h = 0.001 \text{ Wm}^{-3}$ is evenly distributed over the whole design domain.

Regarding problem (19), the reference volume and the volume fraction are $V_{ref} = a_1a_2a_3$ and $\gamma = 0.3$, respectively. An initial guess characterized by a uniform pseudo-density field $\rho(\zeta_1, \zeta_2, \zeta_3) = \gamma$ is considered. The reference thermal compliance associated to the starting point is $W_{ref} = 503.32 \text{ WK}$. Furthermore, three design requirements have been added to the CNLPP formulation of Equation (19): a minimum member size requirement $d_{min} = 0.7 \text{ m}$
and a double orthogonal symmetry constraint with respect to planes $x_1 = \frac{a_1}{2}$ and $x_2 = \frac{a_2}{2}$. According to the methodology described in [35], the minimum length scale requirement is fulfilled by choosing a B-spline/NURBS entity characterized by $p_j = 2$, ($j = 1, 2, 3$) and $n_{CP} = 30 \times 30 \times 30$.

The optimized topologies are illustrated in Figure 7, for both B-spline and NURBS entities: numerical results are provided in terms of dimensionless thermal compliance $\frac{W}{W_{ref}}$, number of iterations $N_{iter}$ and measured minimum member size $d_{min}^m$. For each solution, the requirement on the volume fraction is active. The same remarks already done for benchmark problem BK2-2D can be repeated here.

<figure><img src="image_25.png" /><figcaption>**Figure 7.** Benchmark problem BK1-3D: (a) B-spline and (b) NURBS solutions of problem (19) with $p_1 = p_2 = 2$ and $n_{CP} = 30 \times 30 \times 30$.</figcaption></figure>

**5. Conclusions**

In this work, steady-state thermal conduction TO problems have been revisited in the context of a special SIMP approach reformulated in the framework of the NURBS hyper-surfaces theory. Some features of this approach have to be highlighted.

* NURBS hyper-surfaces bring three advantages: (a) unlike the classical SIMP approach, a filter zone does not need to be introduced because the NURBS local support establishes an implicit relationship among the pseudo-density of contiguous mesh elements; (b) when compared to the classical SIMP approach, the number of design variables is reduced; (c) the CAD reconstruction of the boundary of the optimized topology is an easy task.

* A sensitivity analysis of the optimized topology to the NURBS integer parameters has been performed. Some general rules about the choice of the integer parameters can be drawn: the higher the number of CPs (for a given degree) or the lower the degree (for a given number of CPs) the smaller the objective function value, for both B-spline and NURBS solutions.

* The role of NURBS weights has been evaluated. In particular, by keeping the same number of CPs and the same degrees, the objective function of the NURBS solution is lower than that of the B-spline counterpart.

* The minimum-length scale requirement is correctly taken into account, without introducing an explicit optimization constraint, by properly setting the integer parameters
of the NURBS entity. This is one of the most important advantages of the NURBS-
based SIMP approach.

* The topological descriptor is not related to the mesh of the FE model. The FE model is only used to assess the physical responses of the problem at hand. The optimized topology can be easily extracted at the end of the optimization process because it is described by means of a pure geometrical entity, i.e., a CAD-compatible entity.

Relevant paths to take for future works include: (i) the integration of the transient regime within the TO problem formulation; (ii) the implementation of thermo-mechanics problems (by considering both weak and strong couplings); (iii) the extension of the proposed approach to design-dependent boundary conditions in order to deal with more realistic engineering applications. Research is ongoing on all the above aspects.

**Author Contributions:** M.M.: conceptualization, methodology, software, validation, formal analysis, investigation, resources, writing—original draft preparation, writing—review and editing, supervision, funding acquisition. K.R.: investigation, writing—original draft preparation. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work has been funded by Nouvelle-Aquitaine region through OCEAN-ALM project.

**Data Availability Statement:** The raw/processed data required to reproduce these findings cannot be shared at this time as the data also forms part of an ongoing study.

**Conflicts of Interest:** The authors declare no conflict of interest.

## References

1. Garimella, V.; Fleischer, A.; Murthy, J.; Keshavarzi, A.; Prasher, R.; Patel, C.; Bhavnani, S.; Venkatasubramanian, R.; Mahajan, R.; Joshi, Y.; et al. Thermal Challenges in Next-Generation Electronic Systems. *IEEE Trans. Compon. Packag. Technol.* **2008**, *31*, 2577–2592. [CrossRef]

2. Vassighi, A.; Sachdev, M. *Thermal and Power Management of Integrated Circuits*; Springer Science & Business Media: New York, NY, USA, 2006.

3. Dong, X.; Liu, X. Multi-objective optimal design of microchannel cooling heat sink using topology optimization method. *Numer. Heat Transf. Part A Appl.* **2020**, *77*, 90–104. [CrossRef]

4. Zhao, M.; Tian, Y.; Hu, M.; Zhang, F.; Yang, M. Topology optimization of fins for energy storage tank with phase change material. *Numer. Heat Transf. Part A Appl.* **2020**, *77*, 284–301. [CrossRef]

5. Wang, S.; Li, S.; Luo, L.; Zhao, Z.; Du, W.; Sundén, B. A high temperature turbine blade heat transfer multilevel design platform. *Numer. Heat Transf. Part A Appl.* **2021**, *79*, 122–145. [CrossRef]

6. Sun, S.; Liebersbach, P.; Qian, X. 3D topology optimization of heat sinks for liquid cooling. *Appl. Therm. Eng.* **2020**, *178*, 115540. [CrossRef]

7. Zhang, B.; Zhu, J.; Gao, L. Topology optimization design of nanofluid-cooled microchannel heat sink with temperature-dependent fluid properties. *Appl. Therm. Eng.* **2020**, *176*, 115354. [CrossRef]

8. Malekipour, E.; Tovar, A.; El-Mounayri, H. Heat Conduction and Geometry Topology Optimization of Support Structure in Laser-Based Additive Manufacturing. *Mech. Addit. Adv. Manuf.* **2018**, *9*, 17–27.

9. Bendsøe, M.; Kikuchi, N. Generating optimal topologies in structural design using a homogenization method. *Comput. Methods Appl. Mech. Eng.* **1988**, *71*, 197–224. [CrossRef]

10. Bendsøe, M.; Sigmund, O. *Topology Optimization-Theory, Methods and Applications*; Springer: Berlin/Heidelberg, Germany, 2003.

11. Eschenauer, H.; Olhoff, N. Topology Optimization of Continuum Structures: A Review. *Appl. Mech. Rev.* **2001**, *54*, 331–389. [CrossRef]

12. Suzuki, K.; Kikuchi, N. A homogenization method for shape and topology optimization. *Comput. Methods Appl. Mech. Eng.* **2001**, *93*, 291–318. [CrossRef]

13. Bendsøe, M. Optimal shape design as a material distribution problem. *Struct. Optim.* **1989**, *1*, 193–202. [CrossRef]

14. Sethian, J.; Wiegmann, A. Structural Boundary Design via Level Set and Immersed Interface Methods. *J. Comput. Phys.* **2000**, *163*, 489–528. [CrossRef]

15. Wang, M.; Wang, X.; Guo, D. A level set method for structural topology optimization. *Comput. Methods Appl. Mech. Eng.* **2013**, *192*, 227–246. [CrossRef]

16. Guo, X.; Zhang, W.; Zhang, J.; Yuan, J. Explicit structural topology optimization based on moving morphable components (MMC) with curved skeletons. *Comput. Methods Appl. Mech. Eng.* **2016**, *310*, 711–748. [CrossRef]

17. Xie, Y.M.; Steven, G.P. A simple evolutionary procedure for structural optimization. *Comput. Struct.* **1993**, *49*, 885–896. [CrossRef]

18. Huang, X.; Xie, Y.M. Bi-directional evolutionary topology optimization of continuum structures with one or multiple materials. *Comput. Mech.* **2009**, *43*, 393. [CrossRef]
19. Huang, X.; Xie, Y.M. *Evolutionary Topology Optimization of Continuum Structures: Methods and Applications*; John Wiley & Sons, Ltd: Hoboken, NJ, USA, 2010.

20. Li, Q.; Steven, G.P.; Querin, O.M.; Xie, Y.M. Shape and topology design for heat conduction by evolutionary structural optimization. *Int. J. Heat Mass Transf.* **1999**, *42*, 3361–3371. [CrossRef]

21. Li, Q.; Steven, G.P.; Xie, Y.M.; Querin, O.M. Evolutionary topology optimization for temperature reduction of heat conducting fields. *Int. J. Heat Mass Transf.* **2004**, *47*, 5071–5083. [CrossRef]

22. Zhuang, C.G.; Xiong, Z.H.; Han, D. A level set method for topology optimization of heat conduction problem under multiple load cases. *Comput. Method Appl. Mech. Eng.* **2007**, *196*, 1074–1084. [CrossRef]

23. Yoon, G. Topological design of heat dissipating structure with forced convective heat transfer. *J. Mech. Sci. Technol.* **2010**, *24*, 1225–1233. [CrossRef]

24. Ikonen, T.; Marck, G.; Sóbester, A.; Keane, A. Topology optimization of conductive heat transfer problems using parametric L-systems. *Struct. Multidiscip. Optim.* **2018**, *58*, 1899–1916. [CrossRef]

25. Hu, D.; Zhang, Z.; Li, Q. Numerical study on flow and heat transfer characteristics of microchannel designed using topological optimizations method. *Sci. China Tech. Sci.* **2020**, *63*, 105–115. [CrossRef]

26. Yoon, M.; Koo, B. Topology design optimization of conductive thermal problems subject to design-dependent load using density gradients. *Adv. Mech. Eng.* **2019**, *11*, 1–11. [CrossRef]

27. Iga, A.; Nishiwaki, S.; Izui, K.; Yoshimura, M. Topology optimization for thermal conductors considering design-dependent effects, including heat conduction and convection. *Int. J. Heat Mass Transf.* **2009**, *52*, 2721–2732. [CrossRef]

28. Dbouk, T. A review about the engineering design of optimal heat transfer systems using topology optimization. *Appl. Therm. Eng.* **2017**, *112*, 841–854. [CrossRef]

29. Allaire, G. *Shape Optimization by the Homogenization Method*; Springer: New York, NY, USA, 2002.

30. Guest, J.K.; Prévost, J.H.; Belytschko, T. Achieving minimum length scale in topology optimization using nodal design variables and projection functions. *Int. J. Numer. Methods Eng.* **2003**, *61*, 238–254. [CrossRef]

31. Lazarov, B.S.; Wang, F.; Sigmund, O. Length scale and manufacturability in density-based topology optimization. *Arch. Appl. Mech.* **2016**, *86*, 189–218. [CrossRef]

32. Costa, G.; Montemurro, M.; Pailhès, J. A 2D topology optimisation algorithm in NURBS framework with geometric constraints. *Int. J. Mech. Mater. Des.* **2018**, *14*, 669–696. [CrossRef]

33. Costa, G.; Montemurro, M.; Pailhès, J. NURBS Hyper-surfaces for 3D Topology Optimisation Problems. *Mech. Adv. Mater. Struct.* **2021**, *28*, 665–684. [CrossRef]

34. Costa, G.; Montemurro, M.; Pailhès, J.; Perry, N. Maximum length scale requirement in a topology optimisation method based on NURBS hyper-surfaces. *CIRP Ann.* **2019**, *68*, 153–156. [CrossRef]

35. Costa, G.; Montemurro, M.; Pailhès, J. Minimum Length Scale Control in a NURBS-based SIMP Method. *Comput. Methods Appl. Mech. Eng.* **2019**, *354*, 963–989. [CrossRef]

36. Rodriguez, T.; Montemurro, M.; Le Texier, P.; Pailhès, J. Structural Displacement Requirement in a Topology Optimization Algorithm Based on Isogeometric Entities. *J. Optim. Theory Appl.* **2020**, *184*, 250–276. [CrossRef]

37. Bertolino, G.; Montemurro, M.; Perry, N.; Pourroy, F. An Efficient Hybrid Optimization Strategy for Surface Reconstruction. *Computer Graphics Forum* **2021**. [CrossRef]

38. Costa, G.; Montemurro, M. Eigen-frequencies and harmonic responses in topology optimisation: A CAD-compatible algorithm. *Eng. Struct.* **2020**, *214*, 110602. [CrossRef]

39. Montemurro, M.; Bertolino, G.; Roiné, T. A General Multi-Scale Topology Optimisation Method for Lightweight Lattice Structures Obtained through Additive Manufacturing Technology. *Compos. Struct.* **2021**, *258*, 113360. [CrossRef]

40. Piegl, L.; Tiller, W. *The NURBS Book*; Springer: Berlin/Heidelberg, Germany; New York, NY, USA, 1997.

41. Zienkiewicz, O.C.; Taylor, R.L.; Zhu, J.Z. *The Finite Element Method: Its Basis and Fundamentals*; Elsevier: Amsterdam, The Netherlands, 2005.

42. Svanberg, K. A class of globally convergent optimization methods based on conservative convex separable approximations. *SIAM J. Optim.* **2002**, *12*, 555–573. [CrossRef]

43. Alexandersen, J. Topology Optimization for Convection Problems. Master's Thesis, DTU Mekanik, Kongens Lyngby, Denmark, 2011.