Multiscale Design of Graded Stochastic Cellular Structures for the Heat Transfer Problem

Abstract

Advancesin additive manufacturing technology have expanded the development prospect of structures with complex configurations. Cellular structures have been a hot research topic in recent years for their superior performance and characteristics, such as being lightweight and having high specific strength and good permeability. With a high specific surface area, cellular structures perform noticeably well in heat transfer applications when subjected to a body heat source. In this paper, a scale-separated multiscale design of theVoronoi graded stochastic cellular structure (Voronoi-GSCS) that considers the heat transfer problem is proposed. The design method is composed of three steps: the offline calculation on the microscale, the online optimization on the macroscale, and geometry reconstruction on the full scale. Numerical examples are given to show the effectiveness and superiority of the developed method for designing the Voronoi-GSCS. The results obtained by the solid isotropic material with penalization (SIMP) approach are used for comparison. The size effect analysis was conducted to research the influence of the size of the microstructure on the Voronoi-GSCS. It should be stressed that the smallest struts were larger than the minimum print size of the additive manufacturing so that the Voronoi-GSCS could be conveniently exploited in heat transfer applications.

1. Introduction

Temperature difference is the cause of heat transfer, which exists widely in daily life. In buildings, engines, kilns, and other engineering applications, it is of great significance to realize heat transfer and diffusion. As we know, there are three basic ways of heat transfer: heat radiation, heat conduction, and heat convection. The heat conduction velocity of solid material is greater than that of gaseous material. When subjected to the point of the heat source, heat can be quickly diffused through solid materials around the source. The more solid material there is around the source, the higher the heat transfer efficiency. As a result, a traditional single-scale solid structure has high heat transfer efficiency. However, the opposite is true for a body heat source. Body heat sources are widely distributed in the design domain, and concentrated solid material distribution is not conducive to heat diffusion. More efficient forms of structures need to be considered. With a high specific surface area, the cellular structure has gained great attention. In the heat transfer process of a cellular structure, the heat transfer path is divided into two: one is heat conduction through the solid material, the other is heat conduction and heat conduction through the hole. The porous property increases the heat transfer path of solid materials and speeds up the heat transfer efficiency. Consequently, the cellular structure provides excellent heat transfer properties when subjected to a body heat source. With the development of additive manufacturing, the preparation technology of complex cellular structures is becoming more and more mature. Nowadays, cellular structures have become a research hotspot in the interdisciplinary field of materials and structures [1].

Cellular structures are multiscale structures with two or more scales. The scales are interrelated and influence each other. A change in the microstructure affects the mechanical behavior of the macrostructure, and the macrostructure also affects the configuration of the microstructure. Generally, the scale-separated and scale-related approaches are the two main approaches for designing multiscale structures. Based on the homogenization theory [2,3], the scale-separated approach is of high computational efficiency. However, the size effect and the connectivity between adjacent microstructures should be considered. On the contrary, the scale-related approach solves the size effect and connectivity problem, while too many design variables result in poor efficiency. Until now, various modified approaches were developed based on the above two approaches, such as the full-scale/reduced-scale methods [4,5,6,7,8,9], the ${\mathrm{FE}}^2$ methods [10], and the homogenization methods [11,12,13,14,15].

As a type of structure optimization, topology optimization has more design freedom and can obtain more design space, which is one of the most promising aspects of structure optimization. With the vigorous development in recent decades, several mainstream optimization methods have been gradually formed [16,17,18,19]. In most studies, compliance minimization optimization problems and eigenfrequency maximization optimization problems [20,21,22,23] are the main focuses of structural optimization. The structural stiffness and characteristic frequency are reported in detail. Topology optimization of cellular structures has also made some progress in recent years. The good ductility, overall structural stability, shock absorption, and other superior characteristics have been incisively and vividly reflected, while other properties of the cellular structures should be further addressed in the application field. As we know, efficient heat energy dissipation must be required for artifacts that suffer high temperatures. The cellular structures with a high specific surface area are very suitable for heat transfer-related requirements. However, only a limited number of works concern the heat transfer problem for cellular structures. Liu et al. [24] developed a data-driven M-VCUT optimization method for the heat transfer problem of cellular structures with multiple microstructure prototypes. Montemurro et al. [25] proposed a general approach for the optimization of graded architected cellular structures to predict their effective heat conductivity. Imediegwu et al. [26] presented a multiscale framework for the heat and thermostructural optimization of printable structures using lattice-based microarchitectures. Han et al. [27] developed a novel design method to eliminate constraints in the common ground structure approached for nonuniform lattice structures and applied the method to the design of a heat sink with optimum heat dissipation efficiency. Cheng et al. [28] focused on the heat conductivity of lattice structures with movable features.

The above works focus on the optimization of cellular structures with nonstochastic microstructures, such as box-shaped and x-shaped microstructures. These microstructures are simple in configuration and controlled by fewer parameters. In the optimization process, design variables are few, which is highly favored. However, the nonstochastic microstructures are anisotropic, meaning they are sensitive to boundary conditions. When the loading direction or constraint conditions change, the microstructure configuration will change greatly, resulting in greater instability. For isotropic microstructures, such problems are not encountered in service. Consequently, to overcome the problem and improve the performance of structures, a more isotropic microstructure should be developed [29,30,31,32,33]. In this paper, a multiscale design method for the Voronoi-GSCS for the heat transfer problem is proposed. The design method was based on the scale-separated approach. On the microscale, the offline calculation was conducted to establish the mapping relationship among the relative density, the average aspect ratio of the struts, and the heat conductivity matrix. On the macroscale, the optimization formulation for solving the heat transfer problem was developed by using the density-based method to obtain the relative density distribution. On the full scale, the Voronoi-GSCS was reconstructed based on the mapping relationship and the relative density distribution.

This paper is organized as follows. In Section 2, firstly, the offline calculation is introduced in detail, including the generation and homogenization of the microstructure, the calculation of the relative density of the microstructure, the establishment of the mapping relationship among the relative density, the average aspect ratio of the struts, and the heat conductivity matrix. Then, a modified density-based method is proposed to deal with the heat transfer optimization problem. Finally, the reconstruction process of the Voronoi-GSCS is presented. In Section 3, several examples are given to manifest the effectiveness of the proposed multiscale method for designing the Voronoi-GSCS. Some conclusions are drawn in Section 4.

2. Methodology

In this section, the design process of the Voronoi-GSCS is introduced in detail. The flow chart of the process is shown in Figure 1.

2.1. Offline Calculation on the Microscale

The offline calculation on the microscale was conducted for the modeling and analysis of the stochastic Voronoi microstructure, including the generation and homogenization of the stochastic Voronoi microstructure, the calculation of the microstructural relative density, and the establishment of the microstructural mapping relationship.

In the generation process of the microstructure, a mesh generator [34] was used. Its principle was to regard the nodes and connections within the domain as the truss structure. When the nodes achieved the mechanical equilibrium, a stochastic and uniform node distribution formed the triangle microstructure skeleton, as shown in Figure 2a. Based on the geometrical relationship [35], the Voronoi skeleton was obtained, as shown in Figure 2b. By using the fast marching method [36], the signed distance function $f_d(x,y)$ of the Voronoi skeleton could be calculated, as shown in Figure 2c. The sign distance function is a level set function that gives the distance from a point to a curve. In general, the distance between a point and the boundary of a region is defined as positive inside the boundary, negative outside the boundary, and zero on the boundary. Consequently, the stochastic Voronoi microstructure could be generated by cutting the level set function ${\varphi }_m=f_d(x,y)+t/2$, as shown in Figure 2d.

$$\left\{\begin{array}{cc}\hfill {\varphi }_m& >0\iff \mathrm{solid}\mathrm{material}\hfill \\ \hfill {\varphi }_m& =0\iff \mathrm{material}\mathrm{interface}\hfill \\ \hfill {\varphi }_m& <0\iff \mathrm{void}\mathrm{material},\hfill \end{array}\right.$$

(1)

where $t$ is the width of the strut.

After obtaining the Voronoi microstructure, the homogenization approach was used to calculate the equivalent heat conductivity tensor, $d_{ij}^{\mathrm{H}}$, of the microstructure [37],

$$d_{ij}^{\mathrm{H}}=\frac{1}{|Y_m|}{\int }_{Y_m}({\tilde{T}}_{,p}^i-T_{,q}^i){\mu }_{pq}({\tilde{T}}_{,p}^j-T_{,q}^j)\mathrm{d}Y,$$

(2)

where $|Y_m|$ is the volume of the microstructure, ${\mu }_{pq}$ is the heat conductivity tensor of the material, ${\tilde{T}}_{,p}^i$ represents the prescribed temperature gradient, and T is the generalized temperature field given by

$${\int }_{Y_m}{v}_{,i}{\mu }_{ij}{T}_{,j}^k\mathrm{d}Y={\int }_{Y_m}{v}_{,i}{\mu }_{ij}{\tilde{T}}_{,j}^k\mathrm{d}Y,\forall v\in {\mathbb{T}}_m,$$

(3)

where v is a virtual temperature field.

For a two-dimensional problem, the heat conductivity tensor, $d_{ij}^{\mathrm{H}}$, can be written as a matrix in the following form:

$${\mathbf{d}}_{\mathrm{pre}}^{\mathrm{H}}=\left[\begin{array}{cc}{d}_{xx}& d_{xy}\\ d_{yx}& d_{yy}\end{array}\right],d_{xy}=d_{yx},$$

(4)

where $d_{xx}$ and $d_{yy}$ are the heat conductivity matrices in x and y directions, and $d_{xy}=d_{yx}$ is the coupling term.

In the homogenization process, the heat conductivity of solid material and void material in the microstructure was assumed to be 4 and 0.004, respectively. As can be seen in

Table 1. Results of the homogenized heat conductivity matrix.

Microstructure	${\mathit{d}}_{\mathbf{pre}}^{\mathbf{H}}$	${\mathit{d}}_{\mathbf{post}}^{\mathbf{H}}$
	$\left[\begin{array}{cc}1.2854& 0.003\\ 0.003& 1.2735\end{array}\right]$	$\left[\begin{array}{cc}1.2795& 0\\ 0& 1.2795\end{array}\right]$

, the homogenized heat conductivity matrix was anisotropic with the values of $d_{xx}$ and $d_{yy}$ being very close, and the value of $d_{xy}=d_{yx}$ was much smaller than the value of $d_{xx}$ and $d_{yy}$. A large number of numerical results showed that when the number of struts within the microstructure was large enough, the heat conductivity matrix tended to be isotropic with $d_{xx}=d_{yy}$ and $d_{xy}=d_{yx}=0$. Consequently, the homogenized heat conductivity matrix was further processed as follows:

$${\mathbf{d}}_{\mathrm{post}}^{\mathrm{H}}=\left[\begin{array}{cc}\frac{d_{xx}+d_{yy}}{2}& 0\\ 0& \frac{d_{xx}+d_{yy}}{2}\end{array}\right].$$

(5)

As shown in Figure 3, there are three kinds of elements: solid elements, void elements, and blending elements. A fixed background grid was used to calculate the microstructural relative density,

$${\rho }_m={\int }_{Y_m}H\left({\phi }_m\right)\mathrm{d}\mathsf{\Omega }=\sum _{e=1}^{N_e}{\rho }_m^e,$$

(6)

where the subscript ‘$m$’ denotes the microstructure, H is the Heaviside function, ${\phi }_m$ is the level set function of the microstructure, $N_e$ is the number of elements, and ${\rho }_m^e$ is the relative density of the element calculated by

$${\rho }_m^e={\int }_{Y_m^e}H\left({\phi }_m^e\right)\mathrm{d}Y=\sum _{p=1}^{N_g}H\left({\phi }_m^e\left({\zeta }_p\right)\right)\left|{\mathbf{J}}_e\left({\zeta }_p\right)\right|{\omega }_p,$$

(7)

where ${\zeta }_p$ is the coordinate vector of the Gaussian integration point, ${\mathbf{J}}_e$ is the element Jacobian matrix, ${\omega }_p$ is the weight coefficient, and $N_g$ is the number of integration points.

The key to the offline calculation was the establishment of the mapping relationship among the relative density, ${\rho }_m$, the heat conductivity matrix, $d_{xx}=d_{yy}$, and the average aspect ratio of the struts, $\alpha =\frac{l_{\mathrm{avg}}}{t}$, where $l_{\mathrm{avg}}$ is the average length of struts and t is the width of the struts.

As shown in Figure 4a, a two-term exponential function was used to build the mapping relationship between the relative density, ${\rho }_m$, and the average aspect ratio of the struts, $\alpha$. Verification data were used to check the correctness of the fitting curves.

$$\alpha =a{e}^{b{\rho }_m}+c{e}^{d{\rho }_m},$$

(8)

where $a,b,c$, and $d$ are fitted parameters, as shown in

Table 2. The fitting parameters $a,b,c$, and d in Equation (8).

a	b	c	d
43.67	−13.42	10.08	−2.11

.

It can be seen that the verification data were around the fitting curve, which validated the correctness of the fitting curves. As shown in Figure 4b, the distribution of equivalent heat transfer conductivity had good regularity and showed a good linear shape. In this paper, the relationship between the equivalent heat conductivity and the relative density was given by interpolation.

2.2. Online Optimization on the Macroscale

In this subsection, online optimization was conducted to obtain the relative density distribution of the design domain. The optimization formulation and sensitivity analysis are given and calculated for the heat transfer problem. In the optimization process, an element’s relative density was taken as the design variable:

$$\left(\mathbb{P}\right)\left\{\begin{array}{cc}{\displaystyle \underset{\rho }{min}}\hfill & J_{\mathrm{com}}={\mathit{F}}^{\mathrm{T}}\mathbf{T}\left(\rho \right)\hfill \\ s.t.\hfill & \mathbf{K}\mathit{T}=\mathit{F}\hfill \\ \hfill & V\left(\rho \right)-\overline{V}\le 0\hfill \\ \hfill & 0<{\rho }_{min}\le \rho \le {\rho }_{max}\le 1,\hfill \end{array}\right.$$

(9)

where J is the objection function, $\rho$ is the relative density, $\mathit{F}$ is the external heat load vector, $\mathit{T}$ is the temperature vector, $\mathbf{K}$ is the global heat conductivity matrix, $V={\int }_{\mathsf{\Omega }}\rho d\mathsf{\Omega }$ is the current material volume fraction, $\overline{V}$ is the predefined material volume fraction, and ${\rho }_{\min }$ and ${\rho }_{\max }$ are the lower and upper limit of the relative density, respectively.

The OC optimality criteria [38] were employed in this work to solve the heat transfer problem:

$$\begin{array}{c}\hfill {\rho }_e^{\mathrm{new}}=\left\{\begin{array}{c}\max \left({\rho }_{\min },{\rho }_e-\mathsf{\Delta }\right)\mathrm{if}{\rho }_e{\mathsf{\Lambda }}_e^{\zeta }\le \max \left({\rho }_{\min },{\rho }_e-\mathsf{\Delta }\right)\\ \min \left({\rho }_{\max },{\rho }_e+\mathsf{\Delta }\right)\mathrm{if}{\rho }_e{\mathsf{\Lambda }}_e^{\zeta }\ge \min \left({\rho }_{\max },{\rho }_e-\mathsf{\Delta }\right)\\ {\rho }_e{\mathsf{\Lambda }}_e^{\zeta }\mathrm{otherwise}\end{array}\right.\end{array},$$

(10)

where $\mathsf{\Delta }$ is a positive move limit, $\zeta$ is a damping coefficient, and ${\mathsf{\Lambda }}_e$ is calculated by

$$\begin{array}{c}\hfill {\mathsf{\Lambda }}_e=\frac{-{\displaystyle \frac{\partial J}{\partial {\rho }_e}}}{\overline{\lambda }{\displaystyle \frac{\partial V}{\partial {\rho }_e}}}\end{array},$$

(11)

where $\overline{\lambda }$ is the Lagrange multiplier.

It is stressed that ${\rho }_{\min }$ should be a positive value in the optimization for numerical stability. As can be seen in the mapping relationship, when ${\rho }_{\min }$ was close to zero, the corresponding average aspect ratio would be relatively large and the elastic constants were relatively small. As a result, the optimization process would go through oscillation, which made the reconstruction of the Voronoi-GSCS very difficult. Consequently, in this paper, ${\rho }_{\min }=0.1$ was used.

The finite element method was used to make the design domain discrete and to solve the heat equilibrium equation. The derivative of the objective function concerning an element’s relative density is calculated by

$$\begin{array}{c}\hfill \frac{\partial J}{\partial {\rho }_e}=-{\mathit{T}}^{\mathrm{T}}\frac{\partial \mathbf{K}}{\partial {\rho }_e}\mathit{T}\end{array},$$

(12)

where the derivative of the global heat conductivity matrix concerning an element’s relative density is calculated by

$$\frac{\partial \mathbf{K}}{\partial {\rho }_e}=\sum _{e=1}^{N_e}{\mathbf{C}}_e^{\mathrm{T}}\left({\int }_{Y_e}{\mathbf{B}}_e^{\mathrm{T}}\frac{\partial {\mathbf{d}}_{\mathrm{post}}^{\mathrm{H}}}{\partial {\rho }_e}{\mathbf{B}}_e\mathrm{d}Y\right){\mathbf{C}}_e,$$

(13)

where $N_e$ is the total number of elements, ${\mathbf{C}}_e$ is the expansion matrix, $\mathbf{B}$ is the strain matrix, and ${\displaystyle \frac{\partial {\mathbf{d}}_{\mathrm{post}}^{\mathrm{H}}}{\partial {\rho }_e}}$ can be calculated by,

$${\displaystyle \frac{\partial {\mathbf{d}}_{\mathrm{post}}^{\mathrm{H}}}{\partial {\rho }_e}}=\left[\begin{array}{cc}{\displaystyle \frac{\partial (d_{xx}+d_{yy})}{2\partial {\rho }_e}}& 0\\ 0& {\displaystyle \frac{\partial (d_{xx}+d_{yy})}{2\partial {\rho }_e}}\end{array}\right]$$

(14)

where ${\displaystyle \frac{\partial d_{xx}}{\partial {\rho }_e}}$ and ${\displaystyle \frac{\partial d_{yy}}{\partial {\rho }_e}}$ can be obtained by the mapping relationship.

The derivative of the volume constraint concerning the relative density is calculated by

$$\begin{array}{ccc}\hfill \frac{\partial V}{\partial {\rho }_e}& ={\int }_{Y_e}\mathrm{d}Y\hfill & \hfill =V_e\end{array},$$

(15)

where $V_e$ is the element volume.

2.3. Geometry Reconstruction on the Full-Scale

After obtaining the relative density distribution, the mesh size function $f_h^i\left(x,y\right)$ could be reconstructed based on the mapping relationship between the relative density and the average aspect ratio of the struts. For the relative density distribution, the signed distance function, $f_g^i\left(x,y\right)$, corresponding to the contour of the density, ${\rho }_i$, could be calculated easily. The mesh size function $f_h^i\left(x,y\right)$ could be obtained based on the signed distance function, $f_g^i\left(x,y\right)$:

$$f_h^i\left(x,y\right)=▵h\left[1-H\left(-f_g^i\left(x,y\right)\right)\right]+h_1$$

(16)

where $H\left(x\right)$ is the Heaviside function, $\mathsf{\Delta }h=h_0-h_1$ with $h_0$ and $h_1$ being given as

$$\left\{\begin{array}{cc}\hfill h_0& =\alpha \left({\rho }_{i-1}\right)t\hfill \\ \hfill h_1& =\alpha \left({\rho }_1\right)t\hfill \end{array},\right.$$

(17)

where $t$ is the width of the strut and $\alpha$ is the average aspect ratio of the struts.

After obtaining a series of signed distance functions, $f_h^i\left(x,y\right)$, the final mesh size function, $f_h\left(x,y\right)$, is calculated by

$$f_h\left(x,y\right)=\min \left(f_h^1\left(x,y\right),f_h^2\left(x,y\right),\cdots ,f_h^{N_{\rho }}\left(x,y\right),f_h^{N_{\rho }+1}\left(x,y\right)\right),$$

(18)

where $f_h^{N_{\rho }+1}$ is the mesh size function corresponding to the contour of the relative density, ${\rho }_{max}$.

With the mesh size function, $f_h\left(x,y\right)$, the skeleton of the design domain can be obtained [34,35]. The signed distance function, $f_g(x,y)$, of the Voronoi skeleton can be calculated by using the fast marching method [36]. The Voronoi cellular structure was reconstructed by cutting the level set function, ${\varphi }_g=f_g(x,y)+t/2$,

$$\left\{\begin{array}{cc}\hfill {\varphi }_g& >0\iff \mathrm{solid}\mathrm{material}\hfill \\ \hfill {\varphi }_g& =0\iff \mathrm{material}\mathrm{interface}\hfill \\ \hfill {\varphi }_g& <0\iff \mathrm{void}\mathrm{material}\hfill \end{array}\right.,$$

(19)

where $t$ is the width of the strut.

3. Numerical Examples

In this section, three numerical examples are presented to verify the effectiveness of the proposed multiscale design for the Voronoi-GSCS. All the examples were based on a square design domain with $L=40$, as shown in Figure 5a and Figure 6a,b. In the first benchmark example, a uniform distributed heating with $T=0.025$ was set in the design domain and $T=0$ was imposed in the left middle edge of the domain. A gradient distributed heating, T = 0.01–0.05, was set in the design domain and $T=0$ was imposed in the left middle edge of the domain in the second example. In the third example, a uniform distributed heating, $T=0.025$, was set in the design domain and $T=0$ was imposed in the four corners of the domain. The predefined material volume limit was 0.4. The lower and upper limit of the relative density was 0.1 and 0.9, respectively. All the quantities were dimensionless in the examples.

3.1. 2D Benchmark Example

The first benchmark example was used to demonstrate the generation process and the superiority of the Voronoi-GSCS in the heat transfer problem. The size effect of the multiscale design method for the cellular structure was also surveyed.

As shown in Figure 5b, the relative density distribution was obtained by the modified density method in the online optimization. Based on the relative density distribution and the relationship between the average aspect ratio and the relative density, the average aspect ratio distribution could be obtained, as shown in Figure 5c. It can be seen that the region with a high relative density has a smaller corresponding size function. Conversely, the region with a high relative density has a larger area size function. Then, the Voronoi skeleton could be generated by the mesh generator [34] and the coupling relationship [35], as shown in Figure 5d. The cellular structure skeleton corresponding to the region with a large size function was sparse, and the cellular structure skeleton corresponding to the region with a large size function was dense. The signed distance function of the Voronoi skeleton could be calculated by the fast marching method [36], as shown in Figure 5e. The Voronoi-GSCS was reconstructed by cutting the level set function, as shown in Figure 5f. It can be seen that the microstructures were connected smoothly in the whole design domain.

For comparison, the optimized results obtained by the SIMP methods are given. The SIMP method is based on the work of [39]. The penalty coefficient was 3 and the minimum filter radius was 1.2. The relative density distribution is shown in Figure 7a, which shows the dendritic distribution. The value of relative density became larger when it was close to the left edge with $T=0$. It can be observed in Figure 5b and Figure 7a that the two relative densities showed similar distribution trends. The relative density distribution obtained by the SIMP method had more bifurcation, while the relative density distribution obtained by the method developed in this paper was more concentrated.

As shown in Figure 7b, the heat compliance of the single-scale structure was 4438.30. The heat compliance of the cellular structure was 3737.67, which was $15.79\%$ smaller than the single structure. This indicated that the designed cellular structure had great superiority in heat transfer.

As we know, the size effect should be considered in the scale-separated-based multiscale design method for cellular structures. In the present multiscale method, the size effect was surveyed by adjusting the width of the struts, which determined the aspect ratio of the struts. As shown in Figure 8, the cellular structures obtained by using different widths are given. The widths of the struts for the three cellular structures were 0.08, 0.06, and 0.04, respectively. As can be seen, the three cellular structures were similar to each other. Additionally, the convergence history and objective function value were nearly the same for the full-scale reconstruction of the cellular structure, which was based on the same relative density distribution. The difference was due to the fact that the cellular structure became denser and the time cost became higher when the width of the struts became smaller. In addition, the gradient change in the structure was more obvious with a smaller width. In conclusion, it manifested that the size effect on the cellular structure design in the present multiscale method was not significant.

3.2. Extensive Numerical Examples

In this subsection, two extensive numerical examples are described to verify the effectiveness of the proposed multiscale design method in complicated heat transfer problems. The models are shown in Figure 6. The optimized results obtained by the SIMP method were used for comparison.

As shown in Figure 9, in the gradient-distributed heating example, the relative density distribution of the single-scale structure was asymmetric and in dendritic form. The dendritic distribution was clustered toward the position with $T=0$. Due to gradient-distributed heating, the area of relative density distribution in the high temperature was larger than that in the regions with lower temperatures. The relative density distribution of the cellular structure was also asymmetric and in dendritic form. The dendritic distribution was more concentrated towards the position with $T=0$ and high-temperature regions. As can be seen in the obtained Voronoi-GSCS, in the area of low relative density, the distribution of holes was sparse, and in the area of high relative density, the distribution of holes was dense. The distribution of holes showed a gradient trend.

As shown in Figure 10, in the uniform distributed heating (II), the relative density distribution was symmetric along the diagonal line and around the four positions where the predefined temperature was $T=0$. The obtained relative density distributions of the single-scale structure and cellular structure were nearly the same. The reconstructed Voronoi-GSCS was dense around the four positions and sparse in other regions. The microstructure was well-connected without connection mismatch.

The heat compliances of the extensive numerical examples are given in Figure 11. In the gradient-distributed heating example, the heat compliances for the single-scale structure and cellular structure were 1232.82 and 1028.18, respectively. In the uniform distributed heating (II), the heat compliances for the single-scale structure and cellular structure were 2704.11 and 2573.38, respectively. The heat transfer performance improved by 16.60% and 4.83%, respectively. The advantages of the cellular structure were verified again by the extensive numerical example.

4. Conclusions

In this paper, a scale-separated design of the Voronoi-GSCS for the heat transfer problem was proposed. On the microscale, the offline calculation proceeded. The stochastic Voronoi microstructure was generated and homogenized. The relative density of the microstructure was calculated and the mapping relationship among the relative density, the average aspect ratio of the struts, and the heat conductivity matrix were established. On the macroscale, online optimization was carried out to solve the heat transfer problem. The optimized relative density distribution of the design domain was obtained. On the full scale, the Voronoi-GSCS was reconstructed based on the mapping relationship built on the microscale and the relative density distribution obtained on the macroscale. The body heat source with uniform-distributed heating (I and II) and the gradient-distributed heat transfer problems were considered to demonstrate the effectiveness of the developed multiscale method for designing the Voronoi-GSCS. The results showed that the scale-separated multiscale could ensure good connectivity of the adjacent microstructure. The size effect was surveyed by changing the width of the struts. It could be seen that the size effect has little influence on the obtained graded cellular structure. The heat conductivity of the cellular structure was enhanced considerably compared to the single-scale structure obtained by the SIMP method, showing that the Voronoi-GSCS has great advantages in heat transfer. Additionally, the strut size of the Voronoi-GSCS was larger than the minimum print size of additive manufacturing. Admittedly, much work remains to be performed: (1) the model with arbitrary boundaries should be developed; (2) the graded stochastic cellular structure could be produced with the aid of additive manufacturing, and an experiment could be conducted to verify the properties of the Voronoi-GSCS; (3) the proposed multiscale method could be extended to 3D optimization problems.