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

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.

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,4,5,6,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 aproaches 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,21,22,23,24,25,26,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,33,34,35,36,37,38,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⟶{\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}\dots \sum _{i_N=0}^{n_N}{R}_{i_1,\dots ,i_N}({\zeta }_1,\dots ,{\zeta }_N){\mathbf{P}}_{i_1,\dots ,i_N},$$

(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}\left({\zeta }_k\right)$, $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}\left({\zeta }_k\right)}{{\sum }_{j_1=0}^{n_1}\dots {\sum }_{j_N=0}^{n_N}\left[{\omega }_{j_1,\dots ,j_N}{\prod }_{k=1}^N{N}_{j_k,p_k}\left({\zeta }_k\right)\right]}.$$

(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 $\left[0,1\right]$, whilst ${\mathbf{P}}_{i_1,\dots ,i_N}$ are the CPs coordinates. The j-th CP coordinate ($X_{i_1,\dots ,i_N}^{\left(j\right)}$) is stored in the array ${\mathbf{X}}^{\left(j\right)}$, whose dimensions are $(n_1+1)\times \dots \times (n_N+1)$. The explicit expression of CPs coordinates in ${\mathbb{R}}^M$ is:

$$\begin{array}{c}{\mathbf{P}}_{i_1,\dots ,i_N}=\{X_{i_1,\dots ,i_N}^{\left(1\right)},\dots ,X_{i_1,\dots ,i_N}^{\left(M\right)}\},\\ \\ {\mathbf{X}}^{\left(j\right)}\in {\mathbb{R}}^{(n_1+1)\times \dots \times (n_N+1)},j=1,\dots ,M.\end{array}$$

(3)

Curves and surfaces formulæ 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_{\mathrm{CP}}≔\prod _{i=1}^N(n_i+1).$$

(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}\left({\zeta }_k\right)=\left\{\begin{array}{c}1,\mathrm{if}{v}_{i_k}^{\left(k\right)}\le {\zeta }_k<v_{i_k+1}^{\left(k\right)},\hfill \\ 0,\mathrm{otherwise},\hfill \end{array}\right.$$

(5)

$$\begin{array}{c}{N}_{i_k,q}\left({\zeta }_k\right)=\frac{{\zeta }_k-v_{i_k}^{\left(k\right)}}{v_{i_k+q}^{\left(k\right)}-v_{i_k}^{\left(k\right)}}{N}_{i_k,q-1}\left({\zeta }_k\right)+\frac{v_{i_k+q+1}^{\left(k\right)}-{\zeta }_k}{v_{i_k+q+1}^{\left(k\right)}-v_{i_k+1}^{\left(k\right)}}{N}_{i_k+1,q-1}\left({\zeta }_k\right),\\ q=1,\dots ,p_k,\end{array}$$

(6)

where each blending function is defined on the knot vector

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

(7)

whose dimension is $m_k+1$, with

$$m_k=n_k+p_k+1.$$

(8)

Each knot vector ${\mathbf{v}}_{\left(k\right)}$ 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]:

$$\begin{array}{c}{R}_{i_1,\dots ,i_N}({\zeta }_1,\dots ,{\zeta }_N)\ne 0,\\ \mathrm{if}({\zeta }_1,\dots ,{\zeta }_N)\in \left[v_{i_1}^{\left(1\right)},v_{i_1+p_1+1}^{\left(1\right)}\right[\times \dots \times \left[v_{i_N}^{\left(N\right)},U_{i_N+p_N+1}^{\left(N\right)}\right[.\end{array}$$

(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 $\mathrm{O}(x_1,x_2,x_3)$:

$$\mathcal{D}≔\{{\mathbf{x}}^{\mathrm{T}}=\left\{x_1,x_2,x_3\right\}\in {\mathbb{R}}^3:x_1\in [0,a_1],x_2\in [0,a_2],x_3\in [0,a_3]\},$$

(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:

$$\widehat{\mathbf{K}}\widehat{\mathbf{u}}=\widehat{\mathbf{f}};\widehat{\mathbf{u}},\widehat{\mathbf{f}}\in {\mathbb{R}}^{{\widehat{N}}_{\mathrm{DOF}}},\widehat{\mathbf{K}}\in {\mathbb{R}}^{{\widehat{N}}_{\mathrm{DOF}}\times {\widehat{N}}_{\mathrm{DOF}}},$$

(11)

where ${\widehat{N}}_{\mathrm{DOF}}$ represents the overall number of degrees of freedom (DOFs) before applying the BCs, and $\widehat{\mathbf{K}}$ is the non-reduced (singular) conductivity matrix of the FE model, while $\widehat{\mathbf{f}}$ and $\widehat{\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:

$$\left[\begin{array}{cc}\mathbf{K}& {\mathbf{K}}_{\mathrm{BC}}\\ {\mathbf{K}}_{\mathrm{BC}}^{\mathrm{T}}& \tilde{\mathbf{K}}\end{array}\right]\left(\begin{array}{c}\mathbf{u}\\ {\mathbf{u}}_{\mathrm{BC}}\end{array}\right)=\left(\begin{array}{c}\mathbf{f}\\ \mathbf{r}\end{array}\right),$$

(12)

with:

$$\begin{array}{c}\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}}},\hfill \\ {\mathbf{K}}_{\mathrm{BC}}\in {\mathbb{R}}^{N_{\mathrm{DOF}}\times N_{\mathrm{BC}}},\widehat{\mathbf{K}}\in {\mathbb{R}}^{N_{\mathrm{BC}}\times N_{\mathrm{BC}}},{\widehat{N}}_{\mathrm{DOF}}=N_{\mathrm{DOF}}+N_{\mathrm{BC}},\hfill \end{array}$$

(13)

where $N_{\mathrm{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 $\widehat{\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}.$$

(14)

In the SIMP approach, the material domain $\Omega \subseteq \mathcal{D}$ is identified by means of a pseudo-density function $\rho \left(\mathbf{x}\right)\in [0,1]$ for $\mathbf{x}\in \mathcal{D}$: $\rho \left(\mathbf{x}\right)=0$ denotes absence of material, whilst $\rho \left(\mathbf{x}\right)=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{array}{cc}\hfill & \mathbf{K}≔\sum _{e=1}^{N_e}{\rho }_e^{\alpha }{\mathbf{L}}_e^{\mathrm{T}}{\mathbf{K}}_e^0{\mathbf{L}}_e=\sum _{e=1}^{N_e}{\mathbf{L}}_e^{\mathrm{T}}{\mathbf{K}}_e{\mathbf{L}}_e,\hfill \\ \hfill & {\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}}},\hfill \end{array}$$

(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}.$$

(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},j=1,2,3.$$

(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 ${\mathbf{\xi }}_1$ and ${\mathbf{\xi }}_2$, respectively, defined as:

$${\mathbf{\xi }}_1^{\mathrm{T}}≔({\rho }_{0,0,0},\dots ,{\rho }_{n_1,n_2,n_3}),{\mathbf{\xi }}_2^{\mathrm{T}}≔(w_{0,0,0},\dots ,w_{n_1,n_2,n_3}).$$

(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:

$$\begin{array}{ccc}{\displaystyle \underset{{\mathbf{\xi }}_1,{\mathbf{\xi }}_2}{min}}{\displaystyle \frac{W}{W_{\mathrm{ref}}}},\hfill & \mathrm{s}.\mathrm{t}.:& \hfill \left\{\begin{array}{c}\mathbf{K}\mathbf{u}=\mathbf{f},\hfill \\ {\displaystyle \frac{V}{V_{\mathrm{ref}}}}-\gamma \le 0,\hfill \\ {\mathbf{\xi }}_{1k}\in [{\rho }_{\min },{\rho }_{\max }],{\mathbf{\xi }}_{2k}\in [{\omega }_{\min },{\omega }_{\max }],\hfill \\ \forall k=1,\dots ,n_{\mathrm{CP}}.\hfill \end{array}\right.\end{array}$$

(19)

In Equation (19), $V_{\mathrm{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 }_{\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_{\mathrm{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,$$

(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).$$

(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}}_{\mathrm{BC}}=0$) and the volume read

$$\frac{\partial W}{\partial {\xi }_{ik}}=-\alpha \sum _{e\in S_k}\frac{w_e}{{\rho }_e}\frac{\partial {\rho }_e}{\partial {\xi }_{ik}},i=1,2,k=1,\dots ,n_{\mathrm{CP}},$$

(22)

$$\frac{\partial V}{\partial {\xi }_{ik}}=\sum _{e\in S_k}{V}_e\frac{\partial {\rho }_e}{\partial {\xi }_{ik}},i=1,2,k=1,\dots ,n_{\mathrm{CP}},$$

(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 {\xi }_{ik}}$ reads

$$\frac{\partial {\rho }_e}{\partial {\xi }_{ik}}=\left\{\begin{array}{cc}{R}_k^e,\hfill & \mathrm{if}i=1,\hfill \\ {\displaystyle \frac{R_k^e}{{\xi }_{2k}}}\left({\xi }_{1k}-{\rho }_e\right),\hfill & \mathrm{if}i=2.\hfill \end{array}\right.$$

(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 the GC-MMA algorithm are given in

Table 1. GC-MMA algorithm parameters.

Parameter	Value
move	0.1
albefa	0.1
Stop Criterion	Value
Maximum n. of function evaluations	$100\times n_{\mathrm{var}}$
Maximum n. of iterations	1000
Tolerance on objective function	${10}^{-12}$
Tolerance on constraints	${10}^{-12}$
Tolerance on input variables change	${10}^{-12}$
Tolerance on Karush–Kuhn–Tucker norm	${10}^{-6}$

.

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 $\left[0,1\right]$ for each benchmark problem.

The material thermal conductivity used for all the considered test cases is $k_m=1$ Wm⁻¹K⁻¹.

4.1. 2D Benchmark Problems

Two 2D benchmark problems are considered here: for each problem, the reference volume $V_{\mathrm{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$ m, $t=1$ m. A heating source $s_h=0.001$ Wm⁻² 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$ m, while, on the rest of the boundary the adiabatic wall condition ($\mathbf{q}=\mathbf{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 $\left(x_1\in [\frac{a_1-a_t}{2},\frac{a_1+a_t}{2}],x_2=j{a}_2\right)$ and $\left(x_1=j{a}_1,x_2\in [\frac{a_2-a_t}{2},\frac{a_2+a_t}{2}]\right)$, with $j=0,1$, whilst a heat flux $\mathbf{q}=\mathbf{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_{\mathrm{ref}}=a_1{a}_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_{\mathrm{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.

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_{\mathrm{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 ${\displaystyle \frac{W}{W_{\mathrm{ref}}}}$ and number of iterations $N_{\mathrm{iter}}$ to achieve convergence for B-spline and NURBS entities in Figure 2 and Figure 3, respectively. For each solution, the requirement on the volume fraction is always satisfied.

A synthesis of the results illustrated in Figure 2 and Figure 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.

• 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_{\mathrm{CP}}=68\times 68$ CPs are used for this analysis. Moreover, the topology is constrained to be symmetric with respect to planes $x_1=\frac{a_1}{2}$ and $x_2=\frac{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 ${\displaystyle \frac{W}{W_{\mathrm{ref}}}}$, number of iterations $N_{\mathrm{iter}}$ to achieve convergence and measured minimum member size $d_{\min }^{\mathrm{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.

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.

The FE model is made of $N_e=\mathrm{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 $\left(x_1\in \left[\frac{a_1-a_t}{2},\frac{a_1+a_t}{2}\right],x_2\in \left[\frac{a_2-a_t}{2},\frac{a_2+a_t}{2}\right],x_3=a_3\right)$, while a zero heat flux is imposed on the rest of the boundary. A heating source $s_h=0.001$ Wm⁻³ is evenly distributed over the whole design domain.

Regarding problem (19), the reference volume and the volume fraction are $V_{\mathrm{ref}}=a_1{a}_2{a}_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_{\mathrm{ref}}=503.32$ WK. Furthermore, three design requirements have been added to the CNLPP formulation of Equation (19): a minimum member size requirement $d_{\min }=0.7$ 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_{\mathrm{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 ${\displaystyle \frac{W}{W_{\mathrm{ref}}}}$, number of iterations $N_{\mathrm{iter}}$ and measured minimum member size $d_{\min }^{\mathrm{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.

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.