# Target-driven splitting SPH optimization of thermal conductivity
distribution

Bo Zhang<EMAIL_ADDRESS>Chi Zhang<EMAIL_ADDRESS>Xiangyu Hu
<EMAIL_ADDRESS>TUM School of Engineering and Design, Technical University
of Munich,
85748 Garching, Germany  Huawei Technologies Munich Research Center
80992 Munich, Germany

###### Abstract

Efficiently enhancing heat conduction through optimized distribution of a
limited quantity of high thermal conductivity material is paramount in cooling
electronic devices and numerous other applications. This paper introduces a
target-driven all-at-once approach for PDE-constrained optimization and
derives a splitting smoothed particle hydrodynamics (SPH) method for
optimizing the distribution of thermal conductivity in heat conduction
problems. In this method, the optimization iteration of the system is split
into several easily addressed steps. A targeting step is employed to
progressively enforce the direct target, which potentially leads to increased
PDE residuals. Then, these residuals are recovered through an evolution step
of the design variable. After this, a PDE solution step is carried out to
further decrease the PDE residuals, and the system is ready for the next
iteration. Unlike the simulation-based approaches, the present method does not
rely on the adjoint state equation and converged state variable field in each
iteration, and the optimization process is significantly simplified and
accelerated. With the utilization of an implicit SPH splitting operator and a
general numerical regularization formulation, the information propagation is
further accelerated and the numerical stability is greatly enhanced. Typical
examples of heat conduction optimization demonstrate that the current method
yields optimal results comparable to previous methods and exhibits
considerable computational efficiency. Moreover, the optimal results feature
more moderate extreme values, which offers distinct advantages for the easier
selection of appropriate material with high thermal conductivity.

###### keywords:

Heat transfer enhancement, Volume-to-point problem, PDE-constrained
optimization, Target-driven approach, Gradient descent method, Smoothed
Particle Hydrodynamics, Implicit splitting operator

††journal: Elsevier

##  1 Introduction

Electronic devices find wide-ranging applications in various industries,
including aerospace, transportation, network communication, processing, and
manufacturing. The evolution of microelectronics has facilitated the
miniaturization of electronic devices, but this reduction in device size
drastically increases power density and presents considerable challenges for
effective cooling [1]. To guarantee the expected performance and lifespan,
multiple design ideas have been explored to dissipate the generated heat and
lower the operating temperatures [2, 3, 4, 5, 6, 7]. Among these ideas,
inserting high thermal conductivity material [2, 3] has garnered extensive
attention due to its effectiveness and simplicity of implementation. The rapid
innovation in additive manufacturing, specifically in metal 3D printing, has
also opened up new possibilities for implementing this cooling concept [8, 9].

Thus, the primary design challenge lies in finding the optimal distributions
of high thermal conductivity materials to meet specific cooling targets. This
optimization is also known as the volume-to-point (VP) problem [10]. It
involves redistributing a fixed amount (constraint) of material with high
thermal conductivity (design variable) to cool a heat-generating volume within
specified boundaries, with the objective of minimizing the temperature (state
variable) in the given domain. The VP problem is a typical example of PDE-
constrained optimization, characterized by the constraint of physical
principles expressed as partial differential equations (PDEs) [11].

The methods for PDE-constrained optimization can be categorized into
simulation-based and all-at-once approaches, depending on how the PDE
constraint is handled [12, 13]. In simulation-based methods, it is eliminated
by obtaining a converged physical solution using existing solvers. Therefore,
the optimization process computes the gradients of state variables with
respect to the design ones at the hypersurface of the physical solution to
determine the searching direction in each iteration. Typical approaches for
computing gradients, whether done explicitly or implicitly, include the
adjoint technique [14, 15], automatic differentiation (AD) [16] and artificial
neural network (ANN) [17, 18]. Although the simulation-based methods are
conceptually appealing and widely used in solving VP problems [19, 20, 21],
they require repeated and costly solutions of the PDEs, even in the initial
stages when the design variables are still far from their optimal values. All-
at-once methods, which explicitly maintain the PDE constraint as another
optimal target and treat both the state and design variables equally. A clear
advantage is that they do not require repeated PDE solutions but satisfy the
PDE constraint only at the termination of optimization. Some all-at-once
methods, such as the augmented Lagrangian method [22] and sequential quadratic
programming method [23], have been proposed for addressing general PDE-
constrained optimization problems. However, these methods are rare in solving
VP problems due to their inherent complexity and high computational costs.

In early works on simulation-based methods for tackling the VP problem, some
indirect global optimal principles, such as temperature gradient field
homogenization [24, 25, 26], entropy generation minimization (EGM) [27, 28,
29, 30], and entransy dissipation extremum (EDE) [26, 31, 32, 33, 34], have
been proposed and adopted. Although these principles were expected to be
equivalent to the direct target of the lowest average temperature, it has been
argued that they are not equivalent and directly utilizing these principles
for different targets or boundaries may lead to sub-optimal or even
detrimental solutions [32]. Therefore, direct optimization targets, including
the lowest average [35, 15] and minimizing hot spot temperature [18, 36], have
also been employed recently.

In addition to the aforementioned iterative methods based on deterministic
principles, some stochastic approaches, such as the bionic optimization (BO)
method [37, 38], cellular optimization (CA) method [39, 40], simulated
annealing (SA) and genetic algorithm (GA) method [41], have been explored to
address VP problems. Topology optimization, for example the asymptotes
algorithm (MMA) [42] and density-based method [43], has also gained attention
for optimal heat conduction problems. While these stochastic methods offer
novel avenues for achieving optimal solutions, they often deal with a large
number of variables in the spatial discretization of the domain, leading to
reduced efficiency and suitability for large-scale and reliability-sensitive
problems. Therefore, developing an efficient yet straightforward method,
particularly leveraging the direct optimization target and the all-at-once
concept, holds great promise for addressing heat conduction optimization
problems.

With these in mind, this paper introduces a target-driven all-at-once method
for PDE-constrained optimization problems and derives a smoothed particle
hydrodynamics (SPH) method for optimizing the thermal conductivity
distribution to minimize the average temperature. In this method, the
optimization iteration is divided into several easily handled steps. A
targeting step is adopted ot progressively impose the direct target,
potentially resulting in increased PDE residuals. Then, through an evolution
step of the design variables, these residuals are subsequently recovered.
Following this, a PDE solution step is performed to further decrease the PDE
residual and prepare for the next iteration. The novelty of this work can be
summarized in three aspects. Firstly, as the split steps are only weakly
coupled with each other, compared to previous all-at-once approaches, the
present updating of both state and design variables is much easier to handle.
Secondly, by leveraging the splitting-operator SPH method, implicit updating
is achieved without the inversion of large-size matrices. Thirdly, a general
formulation of regularization has been proposed to achieve numerical stability
when evolving the design variable.

In the following sections, Section 2 provides a brief problem description;
Section 3 introduces the SPH method and the numerical scheme; Section 4
elucidates the proposed optimization method, and Section 5 demonstrates
typical examples to validate the effectiveness of the method. Finally, Section
6 concludes the key findings and outlooks the future work.

##  2 Problem description

We consider the optimization of thermal conductivity distribution for two-
dimensional heat conduction problems, specifically addressing typical VP
problems. As illustrated in Fig. 1,

Figure 1: Illustration of the 2D heat conduction problem. $\Omega$ denotes the
thermal domain; $\varGamma_{T}$ and $\varGamma_{q}$ represent constant surface
temperature and heat flux boundaries, respectively, and $\dot{Q}$ indicates
the presence of an internal heat source.

the thermal domain under consideration, denoted as $\Omega$, is subject to
different boundary conditions and may contain internal heat sources. The
steady-state temperature field within the thermal domain can be determined by
solving the transient heat conduction governing equation given by

$\dfrac{\text{d}T}{\text{d}t}=\nabla\cdot\left(k\nabla
T\right)+\dot{Q}\quad\text{in}\quad\Omega,$ (1)

when the left-hand side (LHS) converges to zero. Here, $T$ represents
temperature, $k$ denotes thermal conductivity, and $\dot{Q}$ is the volume
rate of the internal heat source. Note that, for the sake of simplified
analysis, Eq. (1) is reduced from the general heat conduction equation with
$\rho C=1$, where $\rho$ and $C$ represent density and heat capacity,
respectively.

Two typical boundary conditions are considered here, as shown in Fig. 1. The
Dirichlet boundary condition is given as

$T=T_{b}\quad\text{on}\quad\varGamma_{T},$ (2)

where the surface temperature $T_{b}$ is held constant. The Neumann boundary
condition, which maintains a fixed heat flux rate, can be expressed as,

$-k\nabla T\cdot\boldsymbol{n}=q_{b}\quad\text{on}\quad\varGamma_{q},$ (3)

where $\boldsymbol{n}$ indicates the surface normal vector pointing outward.

The objective of the optimization is to obtain a distribution of thermal
conductivity in $\Omega$ by minimizing the average steady temperature
$\overline{T}$ defined as

$\overline{T}=\frac{1}{V}\int_{V}TdV,$ (4)

where $V$ is the total volume of the thermal domain. Additionally, the average
thermal conductivity is constrained to remain constant by

$\int_{V}kdV=k_{0}V,$ (5)

where $k_{0}$ is a reference value throughout the optimization process.

##  3 Numerical scheme

###  3.1 SPH formulation

In this study, the SPH method is employed to solve the temperature and thermal
conductivity fields. SPH is a fully Lagrangian particle method that was
initially proposed for astrophysical applications [44, 45]. Since its
introduction, SPH has demonstrated significant success in simulating a wide
range of scientific problems, including heat transfer problems [46, 47].

In the SPH scheme, the heat conduction governing equation in Eq. (1) can be
discretized at each SPH particle $i$ located at $\boldsymbol{r}_{i}$ with its
neighboring particles $j$ as follows,

$\displaystyle\dfrac{\text{d}T_{i}}{\text{d}t}=2\sum_{j}\overline{k}_{ij}\displaystyle\dfrac{T_{ij}}{r_{ij}}\nabla_{i}W_{ij}V_{j}+\dot{Q}_{i}.$
(6)

Here,
$\nabla_{i}W_{ij}=\nabla_{i}W\left(\left|\boldsymbol{r}_{ij}\right|,h\right)=\frac{\partial
W_{ij}}{\partial r_{ij}}\boldsymbol{e}_{ij}$, where
$\boldsymbol{r}_{ij}=\boldsymbol{r}_{i}-\boldsymbol{r}_{j}$, $h$ is the
smoothing length and the unit vector
$\boldsymbol{e}_{ij}=\frac{\boldsymbol{r}_{ij}}{r_{ij}}$, represents the
derivative of the kernel function. $T_{ij}=T_{i}-T_{j}$ indicates the inter-
particle temperature difference, and $V_{j}$ denotes the volume of neighboring
particles $j$. The quantity $\overline{k}_{ij}=(k_{i}+k_{j})/2$ is the inter-
particle average thermal conductivity.

Near the domain boundary, several layers of dummy particles are introduced to
enforce different boundary conditions. Implementing the Dirichlet boundary
condition is straightforward and involves imposing the temperatures

$T_{w}=2T_{b}-T_{i},$ (7)

at dummy particles implied by the wall boundary condition [48]. To implement
the Neumann boundary condition, the discretization of Eq. (1) is modified into

$\dfrac{\text{d}T}{\text{d}t}=\nabla\cdot\left(k\nabla
T\right)+\dot{Q}+\dot{Q}^{\varGamma q}\quad\text{in}\quad\Omega,$ (8)

following the Ref. [49], where the heat flux in Eq. (3) is replaced by a
volumetric term $\dot{Q}^{\varGamma q}$, which can be discretized as

$\dot{Q}^{\varGamma
q}_{i}=-q_{b}\displaystyle\sum_{j\in\Omega_{j}}\left(\boldsymbol{n}_{i}+\boldsymbol{n}_{j}\right)\cdot\nabla_{i}W_{ij}V_{j}.$
(9)

Here, $\Omega_{j}$ represents the boundary domain defined by the dummy
particles. The unit vectors $\boldsymbol{n}_{i}$ and $\boldsymbol{n}_{j}$ are
normal to the boundary evaluated at the respective positions of particle $i$
and $j$.

###  3.2 Splitting operator based implicit scheme

It is well known that traditional implicit schemes often require large-scale
matrix inversion or iterations across the entire system, which can lead to
significant memory demands and challenges in parallelization. To overcome
these challenges, we employ a splitting operator based implicit scheme to
advance Eq. (6). The implicit solving step is divided into individual
particle-by-particle operations, and each evolves a small system that is easy
to inverse. One commonly used approach for this purpose is the second-order
Strang splitting technique [50], shown as

$\displaystyle S_{i}^{(\Delta t)}=$ $\displaystyle D_{1}^{(\frac{\Delta
t}{2})}\circ D_{2}^{(\frac{\Delta t}{2})}\circ\dots\circ D_{i}^{(\frac{\Delta
t}{2})}\dots\circ D_{N_{t}-1}^{(\frac{\Delta t}{2})}\circ
D_{N_{t}}^{(\frac{\Delta t}{2})}\circ$ (10) $\displaystyle
D_{N_{t}}^{(\frac{\Delta t}{2})}\circ D_{N_{t}-1}^{(\frac{\Delta
t}{2})}\circ\dots\circ D_{i}^{(\frac{\Delta t}{2})}\dots\circ
D_{2}^{(\frac{\Delta t}{2})}\circ D_{1}^{(\frac{\Delta t}{2})}.$

Here, the operator $S_{i}^{(\Delta t)}$ represents the complete step for
advancing the equation. $N_{t}$ refers to the total number of particles, and
$D_{i}$ represents the splitting operator corresponding to particle $i$. The
update of the variable for the entire field involves a forward sweep of all
particles for half a time step, followed by a backward sweep for another half
time step [51].

Within the local implicit formulation, Eq. (6) can be rewritten as

$\displaystyle\dfrac{\text{d}T_{i}}{\text{d}t}=2\sum_{j}\overline{k}_{ij}\displaystyle\dfrac{T_{ij}^{n+1}}{r_{ij}}\nabla_{i}W_{ij}V_{j}+\dot{Q}^{n+1}_{i},$
(11)

where $T_{ij}^{n+1}=T_{ij}^{n}+\text{d}T_{i}-\text{d}T_{j}$. The terms
$\text{d}T_{i}$ and $\text{d}T_{j}$ represent the incremental changes for
particle $i$ and its neighboring particles $j$ at each advancing time step.
For brevity, we introduce the coefficient

$B_{j}=2\overline{k}_{ij}\displaystyle\dfrac{1}{r_{ij}}\nabla_{i}W_{ij}V_{j}\text{d}t,$
(12)

and the residual of Eq. (11), without considering the increment, has the form

$E_{i}=-\displaystyle\sum_{j}B_{j}T_{ij}^{n}-\dot{Q}^{n+1}_{i}\text{d}t.$ (13)

The implicit formulation of Eq. (11) can be further expressed as

$E_{i}=\left(\displaystyle\sum_{j}B_{j}-1\right)\text{d}T_{i}-\displaystyle\sum_{j}B_{j}\text{d}T_{j}.$
(14)

To determine the incremental changes for temperature, we employ the gradient
descent method [52] by reducing the LHS of Eq. (13) following its gradient.
The gradient $\nabla E_{i}$ with respect to the variable
$\left(\text{d}T_{i},\text{d}T_{1},\text{d}T_{2},\cdots,\text{d}T_{N}\right)^{T}$,
where $N$ gives total number of all neighboring particles, can be obtained as

$\nabla E_{i}=\left(\sum_{j}B_{j}-1,-B_{1},-B_{2},\cdots,-B_{N}\right)^{T}.$
(15)

We set

$\left(\text{d}T_{i},\text{d}T_{1},\text{d}T_{2},\cdots,\text{d}T_{N}\right)^{T}=\eta_{i}\nabla
E_{i},$ (16)

where $\eta_{i}$ represents the learning rate [52] for the particle $i$.
Substituting Eqs. (15) and (16) into Eq. (14), the learning rate can be
obtained as

$\eta_{i}=\left(\left(\sum_{j}B_{j}-1\right)^{2}+\sum_{j}\left(B_{j}\right)^{2}\right)^{-1}E_{i}.$
(17)

According to Eqs. (15) and (16), the incremental change in temperature of
particle $i$ and all its neighbors can be obtained and updated as

$\begin{cases}\displaystyle
T_{i}^{n+1}&=T_{i}^{n}-\text{d}T_{i}=T_{i}^{n}+\eta_{i}\left(\sum_{j}B_{j}-1\right)\\\
\displaystyle T_{1}^{n+1}&=T_{1}^{n}-\text{d}T_{1}=T_{1}^{n}-\eta_{i}B_{1}\\\
\displaystyle T_{2}^{n+1}&=T_{2}^{n}-\text{d}T_{2}=T_{2}^{n}-\eta_{i}B_{2}\\\
&\cdots\\\ \displaystyle
T_{N}^{n+1}&=T_{N}^{n}-\text{d}T_{N}=T_{N}^{n}-\eta_{i}B_{N}\\\ \end{cases}.$
(18)

Note that, Eq. (18) involves updating the variables for particle $i$ and its
neighboring particles simultaneously. When a shared-memory parallelization is
employed, conflict may arise when multiple threads attempt to update the
values of a single particle pair simultaneously. To address this issue, we
have implemented a splitting Cell Linked List method [53]. This method
effectively prevents conflicts by ensuring that neighboring particles are
located in the same cell or in adjacent cells that are distributed among the
same threads. Also note that, for an explicit integration of the thermal
diffusion equation, the maximum allowable time step can be defined as

$\Delta t_{d}=0.5\dfrac{\rho Ch^{2}}{k_{max}}.$ (19)

Since the implicit scheme is employed here for obtaining the steady solution
of the Eq. (6), the time step size is chosen as a large value of $10\Delta
t_{d}$ without considering the temporal accuracy.

##  4 Target-driven optimization

###  4.1 Method overview

The target-driven PDE-constrained optimization method proposed here is based
on the principle of residual recovery. The primary objective is to solve the
PDEs while progressively imposing the target function directly defined by the
state variable. As the latter part can potentially lead to modifications or
even an increase in the PDE residuals, to address this issue, the design
variable undergoes an evolution to recover the original residual, before the
PDE solving continues. This process is repeated iteratively until the fields
of both state and design variables converge and reach steady states. At that
point, the optimal distribution of the design variable is obtained while
satisfying both the target function and PDE constraints. Note that, the
residual recovery approach is a typical all-at-once method since the residual
of PDE only converges upon completing the optimization process.

Although the approach proposed above can theoretically be applied to general
PDE-constrained optimization problems, here we apply it to address heat-
conduction based optimizations The flowchart of the proposed method is
illustrated in Fig. 2

Figure 2: The flowchart of the target-driven optimization procedure.

and the detailed steps are given as follow.

  1. Initialization:

     * (a)

The thermal domain is populated with inner and dummy particles for
implementing the SPH method. The thermal conductivity is initialized as a
uniform distribution with $k_{0}$. The initial temperature is randomly
assigned, and then, the average initial temperature $\overline{T}$ is
calculated.

     * (b)

The PDE residual $e_{i}^{*}$, i.e. the LHS of Eq. (6), is calculated for each
particle within the thermal domain. The overall maximum and average residuals,
denoted as $e_{max}^{*}$, and $e_{ave}^{*}$, respectively, are then
determined.

  2. Step 1: Imposing target.

     * (a)

The target is locally imposed on each particle with a strength $\beta$. In
this case, it can be expressed as $T_{i}^{c}=T_{i}^{*}-\beta$ based on the Eq.
(4). The modified PDE residual $e_{i}^{c}$ on each particle is then updated
with the imposed temperature $T_{i}^{c}$.

  3. Step 2: Design variables evolution.

     * (a)

The thermal conductivity $k_{i}$ on each particle and its neighboring
particles $k_{j}$ are evolved with the aim of recovering the residual
$e_{i}^{*}$ from $e_{i}^{c}$, as detailed in Section 4.2.

     * (b)

Re-normalize the thermal conductivity $k_{i}$ on each particle to fulfill the
constraint presented in Eq. (5).

     * (c)

The thermal conductivity distribution is regularized with the coefficient
$\mu$, employing the diffusion analogy numerical regularization method, as
explained in detail in Section 4.3. The maximum thermal conductivity
variation, denoted as $k^{c}_{max}$, is determined by Eq. (24).

  4. Step 3: PDE solution solving.

     * (a)

The PDE solving advances using the updated thermal conductivity values to
obtain the intermediate temperature field. Meanwhile, the average temperature
for the current state, represented as $\overline{T}^{c}$, is calculated. The
maximum and average residual, $e_{max}^{c}$ and $e_{ave}^{c}$ respectively,
are also updated accordingly.

     * (b)

The PDE solving stops advancing when the current maximum residual
$e_{max}^{c}$ and average residual $e_{ave}^{c}$ are both smaller than their
respective values obtained at the last iteration, satisfying the conditions
$e_{max}^{c}\textless e_{max}^{*}$ and $e_{ave}^{c}\textless e_{ave}^{*}$.

  5. Iterations and termination. After updating the new $e_{max}^{*}$ and $e_{ave}^{*}$, the optimization process repeats from Steps 1 to 3 until the maximum residual reaches the specified threshold, and the variations of temperature and thermal conductivity become smaller than the thresholds. Once these criteria are met, the optimization process is considered complete, and the resulting distribution of $k$ is deemed optimal.

Note that the present method divides the optimization process into small,
easily manageable steps, which simplifies and expedites the optimization
process. The magnitude of target strength $\beta$ is chosen as a small
fraction of the average initial temperature, within the range of $0.5\sim 1$
in the present study. In addition, the $\beta$ is adjusted dynamically,
increasing by a factor of 1.05 when the average temperature is lower than in
the previous optimize iteration and decreasing with a decayed factor of 0.8
when the temperature exceeds that of the previous iteration. Actually,
decreasing $\beta$ is quite important for effective convergence in the late
stages of the optimization. Furthermore, since the convexity of these
optimization problems are not able to be established, the current method,
similar to other general optimization approaches, does not guarantee the
global optimal solution.

###  4.2 Evolution of design variable

In Step 3, the residual $e_{i}^{*}$ for particle $i$ in the PDE is calculated
as

$e_{i}^{*}=\sum_{j}\left(k_{i}+k_{j}\right)\dfrac{T_{ij}}{r_{ij}}\nabla_{i}W_{ij}V_{j}+\dot{Q}_{i}.$
(20)

Once the target is imposed on this particle, the PDE residual deviates from
its original value and will be recovered by modifying the design variable $k$
on particle $i$ and its neighboring particles $j$. This process can be
represented by the pseudo-time evolution of following equation

$\frac{\text{d}k_{i}}{\text{d}\tau}=\sum_{j}\left(T^{c}_{i}-T_{j}\right)\left(k_{i}^{m+1}+k_{j}^{m+1}\right)\dfrac{1}{r_{ij}}\nabla_{i}W_{ij}V_{j}+\dot{Q}_{i}+e_{i}^{*}.$
(21)

Here,
$k_{i}^{m+1}+k_{j}^{m+1}=k_{i}^{m}+\text{d}k_{i}+k_{j}^{m}+\text{d}k_{j}$,
where $m$ is the previous time step and $\text{d}k_{i}$ and $\text{d}k_{j}$
represent increments after the new time step. The implicit splitting operator
introduced in Section 3.2 is utilized. Similar to the Eq. (16), a linear
system is formed with respect to
$\left(\text{d}k_{1},\text{d}k_{2},\cdots,\text{d}k_{N-1},\text{d}k_{N}\right)^{T}$.
Note that, the pseudo-time derivative on the LHS is essential for the stable
evolution of $k$. If this term is omitted, the diagonal entries of the matrix
for the linear system become

$\displaystyle\left(\eta_{1}\sum_{j}^{n}B_{1j},\eta_{2}\sum_{j}^{n}B_{2j},\cdots,\eta_{N}\sum_{j}^{n}B_{Nj}\right)^{T}.$
(22)

It is observed that the magnitudes of diagonal entries in the matrix can be
significantly smaller than those of non-diagonal entries, potentially leading
to numerical instability [54]. On the contrast, when the pseudo-time
derivative term is included, the linear system transforms into

$\displaystyle\left(\eta_{1}\sum_{j}^{n}\left(B_{1j}-1\right),\eta_{2}\sum_{j}^{n}\left(B_{2j}-1\right),\cdots,\eta_{N}\sum_{j}^{n}\left(B_{Nj}-1\right)\right)^{T},$
(23)

whose diagonal entries becomes dominant, and therefore stabilize the evolution
of the design variables. In addition, since $k$ is a material property and
should be non-negative, it is clipped at a lower bound of $0.0001$ during each
iteration.

###  4.3 Numerical regularization

After the evolution of the design variable, it is necessary to apply numerical
regularization, serving two essential purposes. One is that, as previously
mentioned in the Ref. [15], the regularization plays a critical role in
maintaining numerical stability and obtaining a smooth solution. Secondly, it
helps prevent over-fitting and avoids finding trivial local optima only. In
this study, we introduce a diffusion analogy approach for regularizing the
distribution of the design variable, i.e. the thermal conductivity $k$ is
treated as the variable again in the pseudo-time SPH discretized diffusion
equation, given as

$\dfrac{\text{d}k_{i}}{\text{d}\tau}=2\sum_{j}\mu\dfrac{k_{ij}}{r_{ij}}\nabla_{i}W_{ij}V_{j},$
(24)

where $k_{ij}=k_{i}-k_{j}$ and $\mu$ is the artificial diffusion coefficient
used to control the rate of regularization. We choose $\mu$ to be general in
the range $1\sim 2$ according to the target strength. The coefficient also
undergoes a similar dynamical adjustment strategy as the target strength
because smaller $\beta$ usually require less regularization to achieve a
smooth field. Note that the pseudo-time derivative term is also used in Eq.
(24) to ensure that the diagonal is dominant.

##  5 Results and discussion

In this section, we present a set of example problems to validate the
effectiveness of the proposed target-driven PDE-constrained optimization
method (referred to as the TD) for optimizing the thermal conductivity
distribution within a heat conduction domain.

###  5.1 The 2/10 heat sinks with uniform internal heat source

The first set of problems involves a square thermal domain measuring 1m on
each side, with a uniform internal heat source. Two heat sinks, each covering
20% of the side length, are symmetrically positioned at the center of opposite
boundaries, maintaining the constant temperature. The remaining boundaries are
set as adiabatic. Various scenarios are considered, including different heat
source intensities, initial thermal conductivity values, and cases with either
identical or non-identical sink temperatures. As shown in Fig. 3,

Figure 3: Illustration of the problem setups. (a) Problem 1: Identical heat
sinks; (b) Problem 2: Non-identical heat sinks; (c) Problem 3: Identical heat
sinks with higher $k_{0}$; (d) Problem 4: Non-identical heat sinks with higher
$k_{0}$.

the detailed setups for these problems are illustrated as follows:

  * 1.

Problems 1 and 2 feature a uniform heat source with an intensity of $\rm
1000W/m^{3}$, while Problems 3 and 4 have a uniform heat source with an
intensity of $\rm 2000W/m^{3}$.

  * 2.

The initial thermal conductivity, denoted as $k_{0}$, is 1W/(m·K) for Problems
1 and 2, while it is 4W/(m·K) for Problems 3 and 4.

  * 3.

Problems 1 and 3 have identical heat sink temperatures, set at 300K and 280K,
respectively. Problems 2 and 4 involve non-identical heat sink temperatures,
with one of the heat sinks reaching a higher temperature of 350K.

The thermal domain is discretized with 100 particles on each side, resulting
in a total of 10,000 particles. Additionally, four layers of dummy particles
are used to enforce boundary conditions. Note that the problems discussed in
the following sections share the same configuration for the SPH
implementation. Reference solutions for Problems 1 and 2 were obtained using
automatic differentiation (AD) and the temperature gradient homogenization
(TGH) method, and can be found in Ref. [18]. Reference solutions for Problems
3 and 4 were achieved through adjoint analysis (AA) and the TGH method, and
are available in Ref. [15].

The steady temperature distributions with uniform thermal conductivity are
presented in Fig. 4.

Figure 4: Temperature distributions with uniform thermal conductivity for
different scenarios. (a) Problem 1; (b) Problem 2; (c) Problem 3; (d) Problem
4.

All the obtained distributions are in good agreement with the reference
results (see Figs. 4 and 10 of Ref. [18]). In these cases, heat generated by
the internal source is conducted throughout the domain and subsequently
dissipated at the heat sinks. A more pronounced temperature gradient is
observed near the heat sinks, indicating a concentration of heat flux and
resulting in higher temperature within the domain. The average temperatures
obtained are slightly higher than that in the reference, which may be
attributed to variations in resolution and discretization employed.

For these problems, the optimization results obtained by the present method
are displayed in Figs. 5 to 8, and the comparisons with previous work are
summarized in Tables 1 and 2. Problems 1 and 3, featuring identical heat
sinks, result in a symmetrical pattern in the optimized results that align
with the reference results. In cases where the heat sinks are identical, the
problem becomes self-adjoint, making the TGH method equivalent to the AA
method when the target is to minimize the global temperature [55]. Therefore,
all methods adopt the direct target and exhibit similar optimal performance,
although the reference result for the Problem 1 suggests slightly superior
outcome using the AD method. Notably, the present method tends to yield a
higher reduction in temperature for these two problems. The optimization
process homogenized the temperature gradient to balance heat flux through the
domain. Optimized temperature distributions in Figs. 5 and 7 reveal an overall
even temperature gradient, matching with the reference results obtained by
other methods (see Fig. 5(b) of Ref. [18] and Fig. 6(a) of Ref. [15]).
Furthermore, the optimization process results in a significant increase in
thermal conductivity near the heat sinks. The distribution of optimized
thermal conductivity in Figs. 5 and 7 exhibits four distinct peaks prominently
located at the edges of the heat sinks, closely matching the reference (see
Fig. 5(a) of Ref. [18] and Fig. 6(b) of Ref. [15]). These peaks effectively
interact with other boundaries, enhancing the efficient dissipation of
generated heat. It’s worth noting, however, that the peak thermal conductivity
obtained by the present method approximates 14 W/(m·K) and 40 W/(m·K) for
Problems 1 and 3, respectively, which are considerable less than the values in
the reference (beyond 20 W/(m·K) and 50 W/(m·K)).

Figure 5: Present optimized results of the Problem 1. (a) Temperature; (b)
Thermal conductivity.

Figure 6: Present optimized results of the Problem 2. (a) Temperature; (b)
Thermal conductivity.

Figure 7: Present optimized results of the Problem 3. (a) Temperature; (b)
Thermal conductivity.

Figure 8: Present optimized results of the Problem 4. (a) Temperature; (b) Thermal conductivity. Table 1: Summary of result comparisons for problems with identical sinks. Problem | Method |  Original $\overline{T}\left(K\right)$ |  Optimized $\overline{T}\left(K\right)$ |  Reduced $\left(\%\right)$
---|---|---|---|---
1 | TD | $584.30$ | $413.88$ | $29.17$
AD [18] | $582.04$ | $412.73$ | $29.09$
TGH [18] | $582.04$ | $413.45$ | $28.97$
3 | TD method | $417.59$ | $339.12$ | $18.79$
AA [15] | $414.36$ | $336.98$ | $18.67$
TGH [15] | $414.36$ | $336.99$ | $18.67$
Table 2: Summary of result comparisons for problems with non-identical sinks. Problem | Method |  Original $\overline{T}\left(K\right)$ |  Optimized $\overline{T}\left(K\right)$ |  Reduced $\left(\%\right)$
---|---|---|---|---
2 | TD | $609.30$ | $438.79$ | $27.98$
AD [18] | $607.04$ | $436.96$ | $28.02$
TGH [18] | $607.04$ | $438.20$ | $27.81$
4 | TD | $452.73$ | $371.66$ | $17.91$
AA [15] | $449.36$ | $368.61$ | $17.97$
TGH [15] | $449.36$ | $371.55$ | $17.32$

The differences in heat sink configurations in Problems 2 and 4 highlight that
the AA and TGH methods are no longer equivalent. The TGH method transitions
into an indirect target approach, while the direct target methods, TD, AD, and
AA, consistently outperform the TGH. The present method exhibits a slightly
lower temperature reduction ratio compared to the other two direct target
methods. Interestingly, when we decrease the artificial regularization
coefficient ($\mu$), it tends to yield a lower temperature, albeit within
certain limits. However, this adjustment results in significantly higher peak
values of the thermal conductivity. Therefore, we have chosen the current
coefficient value for a balance between optimization performance and the peak
value of thermal conductivity. This choice provides greater flexibility in
selecting high thermally conductive but electrically isolated materials for
cooling electronic devices without sacrificing optimal performance. The
optimized temperature distributions in Figs. 6 and 8 exhibit a generally even
gradient, matching with the reference (see Fig. 11(b) of Ref. [18] and Fig.
7(a) of Ref. [15]). Similarly, the optimized thermal conductivity
distributions in Figs. 6 and 8 still feature four peaks around the sinks.
However, the heights of opposite peaks are no longer equal due to temperature
discrepancies, with larger heights observed near the colder sinks. These
features closely align with the reference results (see Fig. 11(a) of Ref. [18]
and Fig. 7(b) of Ref. [15]), and notably, the present method is always
characterized with lower peaks.

All simulations and optimizations were performed on a computer equipped with 2
Intel(R) Xeon(R) CPU E5-2680 v4 processors. For Problems 1 to 4, it takes
approximately 100 seconds to obtain a steady-state temperature field, and
detailed information regarding the optimization duration is presented in Table
3.

Table 3: Summary of the optimization process for Problems 1 to 4. Problem | Steady time(s) | Optimized iteration | Ratio
---|---|---|---
Loop | Step | Time(s)
1 | $102.8$ | $108$ | $160549$ | $264.3$ | $2.6$
2 | $100.6$ | $94$ | $160541$ | $265.4$ | $2.6$
3 | $110.2$ | $106$ | $210906$ | $341.1$ | $3.1$
4 | $122.1$ | $104$ | $169990$ | $279.5$ | $2.3$

Note that, although the actual optimization time for each run may vary and be
depending on the selected parameters, the shown results suggest that the
present method is quite efficient, as it only takes a few times the
computation cost of obtaining a steady-state solution to achieve the optimized
result.

###  5.2 The 2/10 sinks with Gaussian distributed heat source

Problem 5 closely resembles Problem 4; however, it introduces a non-uniform
distribution of the heat source. Four Gaussian heat sources are symmetrically
placed within the domain at four coordinates: (0.25, 0.25), (0.25, 0.75),
(0.75, 0.25), and (0.75, 0.75), as depicted in Fig. 9. The heat sources
$\dot{Q}_{i}$ are determined by

$\dot{Q}_{i}=C_{i}exp\left[-10\left(\left(x-x_{i}\right)^{2}+\left(y-y_{i}\right)^{2}\right)\right].$
(25)

Here $\left(x_{i},y_{i}\right)$ represents the center point of each heat
source, and $C_{i}\rm=3000W/m^{3}$ denotes the intensity. Fig. 9 illustrates
the cumulative heat source intensity across the entire domain. Reference
solutions obtained by the AA and TGH methods are also available in Ref. [15].

Figure 9: Problem 5. (a) Schematic depiction of the problem setup; (b) Overall
heat source intensity; (c) Temperature distribution under uniform thermal
conductivity.

The steady temperature distribution with uniform thermal conductivity is shown
in Fig. 9, and the obtained average temperature is 518.55K, which agrees well
with the reference.

The current optimization results are presented in Fig. 10, and summarizes the
result comparisons in Table 4.

Figure 10: Present optimized results of the Problem 5. (a) Temperature; (b)
Thermal conductivity.

Due to the non-identical heat sink configuration, it’s evident that the TGH
method is no longer equivalent to the AA method, and the reference results
demonstrate that AA still outperforms the latter. By employing a direct
target, the present method further explores the optimal results and yields
more temperature reduction, being aligned with the performance of the AA
method, as it is also based on the direct target. The optimized temperature
contour shown in Fig. 10 reveals a notably lower values, particularly on the
colder heat sink side, indicating an enhanced cooling capacity. The
temperature gradient is more evenly distributed, aligning with the reference
results (their Fig. 10(a)). Additionally, the optimized distribution in Fig.
10 continues to feature the characteristic four peaks similar to the reference
(their Fig. 10(b)). However, the present thermal conductivity increases
continuously around the heat sink region, as opposed to being concentrated at
isolated spots as in the reference. In addition, the simulation time for a
steady solution is 112.4 seconds for this case, while the optimization time is
980.7 seconds, again indicating quite good efficiency.

Table 4: Summary of result comparisons for Problem 5. Problem | Method |  Original $\overline{T}\left(K\right)$ |  Optimized $\overline{T}\left(K\right)$ |  Reduced $\left(\%\right)$
---|---|---|---|---
5 | TD | $518.55$ | $400.33$ | $22.80$
AA [15] | $517.61$ | $417.35$ | $19.37$
TGH [15] | $517.61$ | $422.85$ | $18.31$

###  5.3 The 1/10 heat sinks with uniform internal heat source

Problem 6 pertains to the utilization of smaller, non-identical heat sinks. As
illustrated in Fig. 11,

Figure 11: Problem 6. (a) Schematic depiction of the problem setup. (b)
Temperature distribution under uniform thermal conductivity.

two heat sinks, each occupying $10\%$ of the side length, are positioned at
the center of the top (350K) and the bottom (300K) boundaries. Reference
solutions obtained by the AA and TGH methods are available in Ref. [35]. The
temperature distribution under uniform thermal conductivity as shown in Fig.
11, is consistent with the reference results (their Fig. 6(d)) and the present
average temperature (365.1K) matches that of the reference (363.3K) also.

The present optimization results are depicted in Fig. 12, and Table 5 provides
a summary of comparisons.

Figure 12: Present optimized results of the Problem 6. (a) Temperature; (b)
Thermal conductivity.

It is observed that the present temperature reduction ratio is higher than
that achieved by the TGH method, although slightly lower than the result
obtained using the AA method. As shown in Fig. 12, the present optimized
temperature distribution presents a smoother profile compared to that of the
reference (their Fig. 6(b)), especially in the vicinity of the high-
temperature sink. From the present thermal conductivity distribution in Fig.
12, the most notable feature is the single peak, in contrast to four peaks
obtained in the reference result (their Fig. 7(b)), which refines the mesh
around the small heat sinks. Note that, due to the introduction of
regularization, the present highest value of thermal conductivity is only
about 20W/(m·K), which is significantly lower than the reference result of
about 200W/(m·K). Moreover, the computation time for a steady solution is
137.9 seconds, while the time required for the entire optimization is 890.8
seconds, again showing quite good efficiency.

Table 5: Summary of results comparison for Problem 6. Problem | Method |  Original $\overline{T}\left(K\right)$ |  Optimized $\overline{T}\left(K\right)$ |  Reduced $\left(\%\right)$
---|---|---|---|---
6 | TD | $365.1$ | $333.0$ | $8.79$
AA [35] | $363.3$ | $331.0$ | $8.89$
TGH [35] | $363.3$ | $331.9$ | $8.64$

###  5.4 The 1/10 sinks with the heat flux heater

Problem 7 explores a thermal domain that features two heat sinks on one side
with different operating temperatures, along with a heat flux heater situated
on the opposite side. As depicted in Fig. 13.

Figure 13: Problem 7. (a) Schematic depiction of the problem setup. (b)
Temperature distribution under uniform thermal conductivity.

a heat flux heater with a magnitude of $\rm 2000W/m^{2}$ is positioned at the
central of the upper boundary, covering $10\%$ of the side length, and two
heat sinks maintained at 300K and 350K, respectively, with the same size are
located along the lower boundary. Reference solutions with the target of
minimizing temperature on the flux boundary can be found in Ref. [35]. The
temperature distribution with uniform thermal conductivity is depicted in Fig.
13, and shows good agreement with the reference (their Fig. 10(d)), although
the average temperature on the flux boundary is slightly lower than that in
the reference, possibly due to the mesh refinement near the flux heater in the
latter.

The present results, with the objective of minimizing the average temperature
of the domain, are shown in Fig. 14, and the comparison of results is
summarized in Table 6.

Figure 14: Present optimized results of the Problem 7. (a) Temperature; (b)
Thermal conductivity.

The optimized thermal conductivity contributes to an overall reduction in
temperature. Note that, even though the present objectives is different from
that of the reference, the pattern of the optimized thermal conductivity in
Fig. 14 closely resembles that of the reference (their Fig. 11(b)). It forms a
bridge with high thermal conductivity between the flux heater and the colder
sink. Compared to the reference results (their Fig. 10(b)), the present
optimized temperature distribution also exhibits a uniform gradient
perpendicular to the line connecting the heater and the colder sink, proving
the enhanced cooling capacity of the colder sink.

Table 6: Summary of results comparison for Problem 7.

Problem | Method |  Original $\overline{T}\left(K\right)$ |  Optimized $\overline{T}\left(K\right)$ |  Reduced $\left(\%\right)$
---|---|---|---|---
7 |  TD1 | $470.7$ | $415.3$ | $11.77$
TD2 | $689.5$ | $519.5$ | $24.66$
AA [35] | $693.2$ | $501.6$ | $27.64$

  * 1

The average temperature across the entire thermal domain.

  * 2

The average temperature along the heat flux boundary.

However, due to the different objectives, the present result shows a higher
averaged temperature on the flux boundary compared to the reference. Again,
the computation time for the steady solution is 197.9 seconds, whereas the
optimization requires 970.8, suggesting good efficiency.

##  6 Conclusion and remark

In this paper, a target-driven all-at-once approach for PDE-constrained
optimization is introduced, and it is applied to optimize thermal conduction
problems. By splitting the optimization iteration into small, easily managed
steps and treating both state and design variables in the same way, the need
for deriving complex adjoint equations and obtaining converged state solutions
at each optimization iteration is eliminated. In addition, the mesh-free,
splitting-operator based implicit SPH method is employed as the underlying
numerical method, and a diffusion-analogy regularization approach is developed
to ensure the numerical stability. Typical examples of thermal conduction
problems demonstrate that the present method is able to achieve quite
efficient optimization with the computational cost generally on the same order
as obtaining a single converged PDE solution. Furthermore, the present optimal
results are comparable to those from the previous work, but with the
industrially relevant advantage of lower extreme values. Note that, as the
target-driven concept used in the present method is not restricted to specific
optimization targets, it may be extended for, such as those on the domain
boundary, topology, or other thermal and fluid dynamics applications, which
are also our future work focus.

## Acknowledgments

The first author would like to acknowledge the financial support provided by
the China Scholarship Council (No.202006230071). C. Zhang and X.Y. Hu would
like to express their gratitude to Deutsche Forschungsgemeinschaft(DFG) for
their sponsorship of this research under grant number DFG HU1527/12-4. The
corresponding code of this work is available on GitHub at
https://github.com/Xiangyu-Hu/SPHinXsys.

## References

  * [1] A. L. Moore, L. Shi, Emerging challenges and materials for thermal management of electronics, Materials today 17 (4) (2014) 163–174.
  * [2] M. Almogbel, A. Bejan, Conduction trees with spacings at the tips, International Journal of Heat and Mass Transfer 42 (20) (1999) 3739–3756.
  * [3] Z. Guo, X. Cheng, Z. Xia, Least dissipation principle of heat transport potential capacity and its application in heat conduction optimization, Chinese science bulletin 48 (2003) 406–410.
  * [4] A. Da Silva, G. Lorenzini, A. Bejan, Distribution of heat sources in vertical open channels with natural convection, International Journal of Heat and Mass Transfer 48 (8) (2005) 1462–1469.
  * [5] K. Chen, S. Wang, M. Song, Optimization of heat source distribution for two-dimensional heat conduction using bionic method, International Journal of Heat and Mass Transfer 93 (2016) 108–117.
  * [6] H. Feng, L. Chen, Z. Xie, F. Sun, Constructal design for a disc-shaped area based on minimum flow time of a flow system, International Journal of Heat and Mass Transfer 84 (2015) 433–439.
  * [7] I. A. Ghani, N. A. C. Sidik, N. Kamaruzaman, Hydrothermal performance of microchannel heat sink: The effect of channel design, International Journal of Heat and Mass Transfer 107 (2017) 21–44.
  * [8] E. M. Dede, S. N. Joshi, F. Zhou, Topology optimization, additive layer manufacturing, and experimental testing of an air-cooled heat sink, Journal of Mechanical Design 137 (11) (2015).
  * [9] T. Zegard, G. H. Paulino, Bridging topology optimization and additive manufacturing, Structural and Multidisciplinary Optimization 53 (2016) 175–192.
  * [10] A. Bejan, Constructal-theory network of conducting paths for cooling a heat generating volume, International Journal of Heat and Mass Transfer 40 (4) (1997) 799–816.
  * [11] J. De los Reyes, Numerical pde-constrained optimization. springerbriefs in optimization, Springer, Cham. doi 10 (2015) 978–3.
  * [12] R. Herzog, Lectures notes algorithms and preconditioning in pde-constrained optimization, Technical Report (2010).
  * [13] R. Herzog, K. Kunisch, Algorithms for pde-constrained optimization, GAMM-Mitteilungen 33 (2) (2010) 163–176.
  * [14] K. Ito, K. Kunisch, Lagrange multiplier approach to variational problems and applications, SIAM, 2008.
  * [15] T. Zhao, X. Wu, Z.-Y. Guo, Optimal thermal conductivity design for the volume-to-point heat conduction problem based on adjoint analysis, Case Studies in Thermal Engineering 40 (2022) 102471.
  * [16] L. B. Rall, G. F. Corliss, An introduction to automatic differentiation, Computational Differentiation: Techniques, Applications, and Tools 89 (1996) 1–18.
  * [17] O. I. Abiodun, A. Jantan, A. E. Omolara, K. V. Dada, N. A. Mohamed, H. Arshad, State-of-the-art in artificial neural network applications: A survey, Heliyon 4 (11) (2018).
  * [18] M. Song, K. Chen, X. Zhang, Optimization of the volume-to-point heat conduction problem with automatic differentiation based approach, International Journal of Heat and Mass Transfer 177 (2021) 121552.
  * [19] Q. Chen, H. Zhu, N. Pan, Z.-Y. Guo, An alternative criterion in heat transfer optimization, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 467 (2128) (2011) 1012–1028.
  * [20] W. Qi, K. Guo, H. Liu, B. Liu, C. Liu, Assessment of two different optimization principles applied in heat conduction, Science Bulletin 60 (23) (2015) 2041–2053.
  * [21] X. Chen, T. Zhao, M.-Q. Zhang, Q. Chen, Entropy and entransy in convective heat transfer optimization: A review and perspective, International Journal of Heat and Mass Transfer 137 (2019) 1191–1220.
  * [22] D. P. Bertsekas, Constrained optimization and Lagrange multiplier methods, Academic press, 2014.
  * [23] W. Alt, K. Malanowski, The lagrange-newton method for nonlinear optimal control problems, Computational Optimization and Applications 2 (1993) 77–100.
  * [24] X. Cheng, Z. Li, Z. Guo, Constructs of highly effective heat transport paths by bionic optimization, Science in China Series E: Technological Sciences 46 (2003) 296–302.
  * [25] X. Cheng, X. Xu, X. Liang, Homogenization of temperature field and temperature gradient field, Science in China Series E: Technological Sciences 52 (2009) 2937–2942.
  * [26] Z.-Y. Guo, H.-Y. Zhu, X.-G. Liang, Entransy—a physical quantity describing heat transfer ability, International Journal of Heat and Mass Transfer 50 (13-14) (2007) 2545–2556.
  * [27] A. Bejan, J. Kestin, Entropy generation through heat and fluid flow (1983).
  * [28] A. Bejan, Entropy generation minimization: The new thermodynamics of finite-size devices and finite-time processes, Journal of Applied Physics 79 (3) (1996) 1191–1218.
  * [29] W. Du, P. Wang, L. Song, L. Cheng, Optimization of volume to point conduction problem based on a novel thermal conductivity discretization algorithm, Chinese Journal of Chemical Engineering 23 (7) (2015) 1161–1168.
  * [30] R. Wang, Z. Xie, Y. Yin, L. Chen, Constructal design of elliptical cylinders with heat generating for entropy generation minimization, Entropy 22 (6) (2020) 651.
  * [31] Q. Chen, X.-G. Liang, Z.-Y. Guo, Entransy theory for the optimization of heat transfer–a review and update, International Journal of Heat and Mass Transfer 63 (2013) 65–81.
  * [32] T. Zhao, Y.-C. Hua, Z.-Y. Guo, Irreversibility evaluation for transport processes revisited, International Journal of Heat and Mass Transfer 189 (2022) 122699.
  * [33] T. Zhao, D. Liu, Q. Chen, A collaborative optimization method for heat transfer systems based on the heat current method and entransy dissipation extremum principle, Applied Thermal Engineering 146 (2019) 635–647.
  * [34] S. Wu, K. Zhang, G. Song, J. Zhu, B. Yao, F. Li, Study on the performance of a miniscale channel heat sink with y-shaped unit channels based on entransy analysis, Applied Thermal Engineering 209 (2022) 118295.
  * [35] Z.-X. Tong, M.-J. Li, J.-J. Yan, W.-Q. Tao, Optimizing thermal conductivity distribution for heat conduction problems with different optimization objectives, International Journal of Heat and Mass Transfer 119 (2018) 343–354.
  * [36] J. Zhang, X. Wu, M. Song, K. Chen, An effective method for hot spot temperature optimization in heat conduction problem, Applied Thermal Engineering (2023) 120325\.
  * [37] Z.-Z. Xia, Z.-X. Li, Z. Guo, Heat conduction optimization: high conductivity constructs based on the principle of biological evolution, in: International Heat Transfer Conference Digital Library, Begel House Inc., 2002.
  * [38] Z.-Z. Xia, X.-G. Cheng, Z.-X. Li, Z.-Y. Guo, Bionic optimization of heat transport paths for heat conduction problems, Journal of Enhanced Heat Transfer 11 (2) (2004).
  * [39] R. Boichot, L. Luo, Y. Fan, Tree-network structure generation for heat conduction by cellular automaton, Energy Conversion and Management 50 (2) (2009) 376–386.
  * [40] R. Boichot, L. Luo, A simple cellular automaton algorithm to optimise heat transfer in complex configurations, International Journal of Exergy 7 (1) (2010) 51–64.
  * [41] X. Xu, X. Liang, J. Ren, Optimization of heat conduction using combinatorial optimization algorithms, International journal of heat and mass transfer 50 (9-10) (2007) 1675–1682.
  * [42] F. H. Burger, J. Dirker, J. P. Meyer, Three-dimensional conductive heat transfer topology optimisation in a cubic domain for the volume-to-surface problem, International Journal of Heat and Mass Transfer 67 (2013) 214–224.
  * [43] M. C. E. Manuel, P. T. Lin, Design explorations of heat conductive pathways, International Journal of Heat and Mass Transfer 104 (2017) 835–851.
  * [44] L. B. Lucy, A numerical approach to the testing of the fission hypothesis, The astronomical journal 82 (1977) 1013–1024.
  * [45] R. A. Gingold, J. J. Monaghan, Smoothed particle hydrodynamics: theory and application to non-spherical stars, Monthly notices of the royal astronomical society 181 (3) (1977) 375–389.
  * [46] V. Vishwakarma, A. K. Das, P. Das, Steady state conduction through 2d irregular bodies by smoothed particle hydrodynamics, International Journal of Heat and Mass Transfer 54 (1-3) (2011) 314–325.
  * [47] Y. Xiao, H. Zhan, Y. Gu, Q. Li, Modeling heat transfer during friction stir welding using a meshless particle method, International journal of Heat and mass Transfer 104 (2017) 288–300.
  * [48] S. Adami, X. Y. Hu, N. A. Adams, A generalized wall boundary condition for smoothed particle hydrodynamics, Journal of Computational Physics 231 (21) (2012) 7057–7075.
  * [49] E. M. Ryan, A. M. Tartakovsky, C. Amon, A novel method for modeling neumann and robin boundary conditions in smoothed particle hydrodynamics, Computer Physics Communications 181 (12) (2010) 2008–2023.
  * [50] G. Strang, On the construction and comparison of difference schemes, SIAM journal on numerical analysis 5 (3) (1968) 506–517.
  * [51] K. Nguyen, A. Caboussat, D. Dabdub, Mass conservative, positive definite integrator for atmospheric chemical dynamics, Atmospheric Environment 43 (40) (2009) 6287–6295.
  * [52] C. M. Bishop, N. M. Nasrabadi, Pattern recognition and machine learning, Vol. 4, Springer, 2006.
  * [53] Y. Zhu, C. Zhang, X. Hu, A dynamic relaxation method with operator splitting and random-choice strategy for sph, Journal of Computational Physics 458 (2022) 111105.
  * [54] R. Bagnara, A unified proof for the convergence of jacobi and gauss–seidel methods, SIAM review 37 (1) (1995) 93–97.
  * [55] J. Alexandersen, O. Sigmund, Revisiting the optimal thickness profile of cooling fins: A one-dimensional analytical study using optimality conditions, in: 2021 20th IEEE Intersociety Conference on Thermal and Thermomechanical Phenomena in Electronic Systems (iTherm), IEEE, 2021, pp. 24–30.