# Design of metamaterial-based heat manipulators by isogeometric shape
optimization

Chintan Jansari  Stéphane P.A. Bordas<EMAIL_ADDRESS>Elena Atroshchenko  Institute of Computational Engineering, Faculty of
Sciences, Technology and Medicine, University of Luxembourg, Luxembourg City,
Luxembourg.  Clyde Visiting Fellow, Department of Mechanical Engineering, The
University of Utah, Salt Lake City, Utah, United States.  School of Civil and
Environmental Engineering, University of New South Wales, Sydney, Australia.

###### Abstract

There has been a growing interest in controlled heat flux manipulation to
increase the efficiency of thermal apparatus. Heat manipulators control and
manipulate heat flow. A key to the effective performance of these heat
manipulators is their thermal design. Such designs can be achieved by a
periodic assembly of unit cells (known as metamaterials or meta-structure),
whose geometry and material properties can be optimized for a specific
objective. In this work, we focus on thermal metamaterial-based heat
manipulators such as thermal concentrator (which concentrates the heat flux in
a specified region of the domain). The main scope of the current work is to
optimize the shape of the heat manipulators using Particle Swarm Optimization
(PSO) method. The geometry is defined using NURBS basis functions due to the
higher smoothness and continuity and the thermal boundary value problem is
solved using Isogeometric Analysis (IGA). Often, nodes as design variables (as
in Lagrange finite element method) generate the serrate shapes of boundaries
which need to be smoothened later. For the NURBS-based boundary with the
control points as design variables, the required smoothness can be predefined
through knot vectors and smoothening in the post-processing can be avoided.
The optimized shape generated by PSO is compared with the other shape
exploited in the literature. The effects of the number of design variables,
the thermal conductivity of the materials used, as well as some of the
geometry parameters on the optimum shapes are also demonstrated.

###### keywords:

Shape optimization , Heat manipulators , Thermal concentrator, Thermal
metamaterials , Particle swarm optimization , Isogeometric analysis.

##  1 Introduction

The modern world of technological and engineering advancements is concerned
about global warming and energy. Renewable energy is central to the advances.
Even though much emphasis is given to producing affordable and clean energy,
improving heat transfer is central to improving the energetical efficiency of
engineering systems. Consequently, new perspectives to enhance the efficiency
and accuracy of thermal appliances by optimizing the heat transfer process are
emerging. One such perspective focuses on the development of heat
manipulators, which are devices that control heat flow. However, the concept
of controlled heat flow is not thoroughly investigated compared to other forms
of energy transport such as electric and photonic currents. The ability to
manage heat flow can lead to the development of thermal equivalents of devices
such as electric transistors, resistors, rectifiers and diodes.

On the other hand, the invention of artificial metamaterials lights up a new
spark in material sciences and, subsequently, in the field of heat transfer.
Because of their well-designed artificial structures, thermal metamaterials
offer outstanding heat transfer capabilities that surpass those of natural
materials. Recognizing these properties, thermal metamaterials were explored
for creating several heat manipulators. The idea of a metamaterial-based
thermal cloak was formulated [1] and later experimentally investigated [2]. A
thermal cloak is a device that minimizes the temperature disturbance caused by
an item and makes it unidentifiable through temperature measurements [2, 3, 4,
5, 6, 7, 8, 9]. Similarly, other heat manipulators such as thermal
concentrator (that concentrates the heat flux in the specified region of the
domain) [2, 3, 4, 10, 8], thermal camouflage (that creates multiple images of
an item away from its actual position) [5, 11], heat flux inverter (that
inverts the direction of the heat flux in the specific region of the domain)
[2] have also been developed. All these heat manipulators work on macro-scale
heat manipulation offered by conductive thermal metamaterials.

As the name suggests, conductive thermal metamaterials guide the flux flow
over a path of interest by engineering the requisite thermal conductivity. The
required thermal conductivity is constrained by the spatial distribution of
member materials. Hence, acknowledging the pivotal role played by the spatial
configuration in the efficacy of a metamaterial-based heat manipulator, we
present a structural optimization method to optimize the relative shape of the
regions in a heat manipulator and their associated thermal conductivity.

Structural optimization aims at identifying a possible material distribution
to achieve the defined objective. The objective could be reducing stress or
concentrating heat flux that can be defined in the quantifiable form. Mainly
two types of structural optimization, shape and topology optimization, are
prevalent. In shape optimization, a structure cannot change its topology, but
it can modify the shapes of topological features. Meanwhile, in topological
topology optimization, the aim is to find the optimal material distribution in
a design domain and optimize the topology.

Once the optimization problem is defined, we need an appropriate algorithm to
find the solution. Optimization algorithms are divided into two types:
gradient-based and gradient-free. Gradient-based algorithms can provide faster
convergence than gradient-free algorithms. But they may get stuck in the local
minima, which depends on the initial guess of the design variables. On the
other hand, gradient-free algorithms, though comparatively slower, are easy to
implement and do not require the optimization problem to be differentiable. We
use one such algorithm called the particle swarm optimization (PSO). PSO was
proposed by Kennedy and Eberhart [12]. The PSO algorithm works on a simple
rule and that makes the process to search for a new solution pool faster. In
comparison with another class of gradient-free algorithm - genetic algorithm
(GA), PSO does not use crossover or mutation operations and works on real
numbers instead of coding operators.

One more aspect of the present work is the application of isogeometric
analysis (IGA) [13] to analyze the thermal boundary value problem and describe
the boundary of the domains and subdomains. As we use IGA, both geometry and
temperature fields are interpolated using non-uniform rational B-splines
(NURBS) basis functions. The NURBS shape functions represent all conic
geometries exactly. In addition to that, as CAD geometries are also based on
NURBS shape functions, IGA creates a seamless integration of modelling and
analysis and is suitable for a fully integrated design-analysis-optimization
model. We can use the control points for the design variables instead of nodes
as in Lagrange finite element methods (FEM). Often, the nodes as design
variables generate the serrate shapes of boundary which needs to be smoothened
later to avoid irregularities. However, we can predefine the required
smoothness through knot vectors in the case of NURBS-based boundaries. In
addition, we can easily provide higher continuity and smoothness across the
elements in IGA compared to FEM. Only the geometries with predefined
smoothness in the boundaries are explored during optimization. With this
predefined higher smoothness through IGA, we can avoid serrate shapes and
other irregularities in the boundaries. Several articles have been published
that exploit the spline-based bases for optimization problems e.g topology
optimization [14, 15, 16, 17, 18, 19, 20], shape optimization [21, 22, 23, 24,
25, 26, 27].

Limited work has been published on optimization for heat manipulators [9, 28,
29]. Existing work has focused on topology optimization with a stochastic
evolution strategy. It is based on immersed-boundary level set method in
combination with finite elements. However, to the best of our knowledge, there
has not been any work done on shape optimization for heat manipulators. The
salient features of this work are:

  * 1.

The shape optimization method, which requires design variables only on the
boundary, is investigated as a tool to obtain better geometries for
metamaterial-based heat manipulators. IGA is utilized for the thermal boundary
value problem, while a gradient-free PSO algorithm is utilized for
optimization.

  * 2.

The use of control points as design variables, as well as the NURBS
approximation for geometry and solution fields, enables easy control of
smoothness and provides accurate geometrical representation.

  * 3.

In addition, the effect of several factors (such as geometric parameters,
material conductivities, boundary conditions, etc.) on the optimized shapes is
studied.

The paper is organized as follows: Section 2 exhibits the governing equations,
the weak formulation, NURBS approximation and matrix formulation of the
thermal boundary problem using IGA. The optimization problem and PSO algorithm
are presented in Section 3. In this work, we present two examples of heat
manipulators, but the proposed method is generic and can be applied to other
heat manipulators as well. Section 4 and Section 5 show detailed examples of a
thermal concentrator and a thermal cloaked concentrator respectively. Section
6 presents the main conclusions of the proposed work.

##  2 Boundary value problem

Consider a metamaterial-based heat manipulator embedded in a plate made of
isotropic material as shown in Figure 1. The whole domain under consideration
can be denoted as $\Omega=\bigcup\limits_{m}\Omega_{m}\in\mathbb{R}^{2}$
bounded by $\Gamma$=$\partial\Omega$. The boundary $\Gamma$ is uniquely
decomposed into two parts $\Gamma_{D}$ & $\Gamma_{N}$, where Dirichlet and
Neumann boundary conditions are applied, respectively. Furthermore, several
internal boundaries in the form of material interfaces are present, which are
collectively denoted by $\Gamma_{I}$,
$\Gamma_{I}=\bigcup\limits_{k}\Gamma_{I_{k}}$. We assume that there is no
internal heat generation and conduction is the only mode of heat transfer
present. Let $T$ be the temperature field over $\Omega$, the steady state heat
conduction equation can be written as,

$\nabla\cdot\left(\bm{\kappa}\nabla T\right)=0\quad\text{in}\quad\Omega,$ (1)

with the boundary conditions

$\displaystyle T$ $\displaystyle=T_{D}\quad\textrm{on}\quad\Gamma_{D},$ (2a)
$\displaystyle\nabla T\cdot\bm{n}$
$\displaystyle=Q_{N}\quad\textrm{on}\quad\Gamma_{N},$ (2b)

where $\bm{\kappa}$ be the thermal conductivity matrix (for isotropic
material, $\bm{\kappa}=\kappa\mathrm{\mathbf{I}_{2}}$ with
$\mathrm{\mathbf{I}_{2}}$ be an identity matrix of $\mathbb{R}^{2}$), $Q_{N}$
the externally applied flux on boundary $\Gamma_{N}$, $T_{D}$ the prescribed
temperature on boundary $\Gamma_{D}$, $\bm{n}$ the unit normal on boundary and
$\nabla=\small{\left({\dfrac{\textstyle\partial}{\textstyle\partial
x},\dfrac{\textstyle\partial}{\textstyle\partial y}}\right)}$.

Figure 1: Domain description of the boundary value problem

Following the standard Galerkin approach, the weak form of the given heat
conduction boundary value problem (Equations (1)-(2)) is written as follows:
Find $T^{h}\in\mathscr{T}^{h}$ such that $\forall
S^{h}\in\mathscr{S}^{h}_{0}$,

$a(T^{h},S^{h})=\ell(S^{h}),$ (3)

with $\mathscr{T}^{h}$ and $\mathscr{S}^{h}_{0}$ be the trial and test spaces,
respectively,

$\displaystyle\mathscr{T}^{h}\subseteq\mathscr{T}=\left\\{T\in\mathbb{H}^{1}({\Omega}),T=T_{D}\hskip
4.26773pt\textrm{on}\ {\Gamma}_{D}\right\\},$ (4a)
$\displaystyle\mathscr{S}^{h}_{0}\subseteq\mathscr{S}_{0}=\left\\{S\in\mathbb{H}^{1}({\Omega}),S=0\hskip
4.26773pt\textrm{on}\ {\Gamma}_{D}\right\\},$ (4b)

and

$a(T^{h},S^{h})=\int_{\Omega}(\nabla S^{h})^{\rm T}\bm{\kappa}\nabla
T^{h}d\Omega,$ (5) $\ell{(S^{h})}=\int_{\Gamma_{N}}(S^{h})^{\rm
T}Q_{N}d\Gamma.$ (6)

In the present study, we assume that the temperature and normal flux is
continuous along each material interface. If the patches connected at
interface $\Gamma_{I}$ are denoted by 1 and 2, then the corresponding
interface conditions are written as,

$\displaystyle\left\llbracket T\right\rrbracket$ $\displaystyle=0\quad$
$\displaystyle\textrm{on}\quad\Gamma_{I},$ (7a)
$\displaystyle\bm{n}^{1}\cdot\bm{\kappa}^{1}\nabla T^{1}$
$\displaystyle=-\bm{n}^{2}\cdot\bm{\kappa}^{2}\nabla T^{2}\quad$
$\displaystyle\textrm{on}\quad\Gamma_{I},$ (7b)

where $\bm{n}^{1}$, $\bm{n}^{2}$; ${\kappa}^{1}$, ${\kappa}^{2}$ and $T^{1}$,
$T^{2}$ are unit normals, conductivity matrices, and temperatures associated
with patch 1 and 2, respectively. The jump operator
$\llbracket\cdot\rrbracket$ is defined in the next paragraph.

The continuity conditions shown in Equation (7) are applied in the given
boundary value problem by modifying the weak form using Nitsche’s method.
Nitsche’s method is a method between the Lagrange multiplier method and the
penalty method. In the weak formulation, it replaces the Lagrange multipliers
by their physical representation, normal flux. In addition to that, it keeps
the coercivity of the bilinear form intact and the variational form
consistent. Nitsche’s method has been successfully applied for patch coupling
in [30, 31]. When Nitsche’s method is applied to couple patches, the linear
form (Equation (6)) on the right side remains the same, while the bilinear
form (Equation (5)) is altered as follows,

$a(T^{h},S^{h})=\int_{\Omega}(\nabla S^{h})^{\rm T}\bm{\kappa}\nabla
T^{h}d\Omega-\int_{\Gamma_{I}}\left(\bm{n}\cdot\\{\bm{\kappa}\nabla
S^{h}\\}\right)^{\rm T}\llbracket T^{h}\rrbracket~{}d\Gamma\\\
-\int_{\Gamma_{I}}\llbracket S^{h}\rrbracket^{\rm
T}\left(\bm{n}\cdot\\{\bm{\kappa}\nabla
T^{h}\\}\right)d\Gamma+\int_{\Gamma_{I}}\beta~{}\llbracket
S^{h}\rrbracket^{\rm T}\llbracket S^{h}\rrbracket~{}d\Gamma,$ (8)

where $\bm{n}$ is the normal at $\Gamma_{I}$ for any one patch from the
patches connected at $\Gamma_{I}$ ($\bm{n}=\bm{n}^{1}=-~{}\bm{n}^{2}$),
$\beta$ is the stabilization parameter. The jump operator
$\llbracket\cdot\rrbracket$ described in Equation (7) and average operator
$\\{\cdot\\}$ in Equation (8) can be defined as,

$\llbracket\theta\rrbracket=\theta^{1}-\theta^{2},\quad\\{\theta\\}=\gamma\theta^{1}+(1-\gamma)\theta^{2},$
(9)

where $\theta$ is a property of interest, superscript denotes the patch (with
which property $\theta$ relates to) and $\gamma$ is the averaging parameter
($0<\gamma<1$).

Figure 2: Parametrization of a point from parametric domain to a point in
physical domain using NURBS basis functions.

Let $\bm{x}\in\Omega$ and $\bm{\xi}$ be the corresponding knot value. We
parametrize the domain using $n$ NURBS, $N_{i}$ (as shown in Figure 2)

$\bm{x}=\sum_{i=1}^{n}\mathbf{P}_{i}N_{i}(\bm{\xi}),$ (10)

where $\mathbf{P}_{i}$ is the $i^{th}$ control point.

As we are using isogeometric analysis with standard Galerkin approach, the
test and trial function both are approximated with the same NURBS shape
functions as geometry. The trial and test function approximation can be
written as,

$\displaystyle T^{h}(\bm{\xi})=\sum_{i=1}^{n}T_{i}N_{i}(\bm{\xi}),$ (11a)
$\displaystyle S^{h}(\bm{\xi})=\sum_{i=1}^{n}S_{i}N_{i}(\bm{\xi}),$ (11b)

where $T_{i}$, $S_{i}$ are the temperature and arbitrary temperature
perturbation at $i^{th}$ control point.

Figure 3: The matching pairs of control points from different patches
connected at interface by Nitsche method

By substituting Equation (11) in Equation (3), we obtain the linear system

$\mathbf{K}\mathbf{T}=\mathbf{Q},$ (12)

where $\mathbf{T}$ is the vector of unknown temperature $T_{I}$ at all control
points. The global stiffness matrix $\mathbf{K}$ and the global flux vector
$\mathbf{F}$ are written as,

$\mathbf{K}=\mathbf{K}^{b}+\mathbf{K}^{n}+(\mathbf{K}^{n})^{\rm
T}+\mathbf{K}^{s},$ (13)
$\mathbf{F}=\int_{\Gamma_{N}}\mathbf{N}^{\textrm{T}}Q_{N}~{}d\Gamma,$ (14)

where $\mathbf{K}^{b}$ is the bulk stiffness matrix defined as below,

$\mathbf{K}^{b}=\sum_{m}\int_{\Omega_{m}}(\mathbf{B}^{m})^{\textrm{T}}\bm{\kappa}^{m}\mathbf{B}^{m}~{}d\Omega.$
(15)

$\mathbf{K}^{n}$ and $\mathbf{K}^{s}$ are the interfacial stiffness matrices,
as these matrices are used to couple the adjacent patches with conditions
given in Equation (7). A point to note, before defining $\mathbf{K}^{n}$ and
$\mathbf{K}^{s}$, is that the connecting patches must have matching control
points at the interface to apply Nitsche’s method as shown in Figure 3.
Following the notation of Equations (7)-(9), $\mathbf{K}^{n}$ and
$\mathbf{K}^{s}$ are given by the following equations.

$\mathbf{K}^{n}=\begin{bmatrix}-\gamma\displaystyle\int_{\Gamma_{I}}(\mathbf{N}^{1})^{\textrm{T}}\bm{n}\bm{\kappa}^{1}\mathbf{B}^{1}~{}d\Gamma&-(1-\gamma)\displaystyle\int_{\Gamma_{I}}(\mathbf{N}^{1})^{\textrm{T}}\bm{n}\bm{\kappa}^{2}\mathbf{B}^{2}~{}d\Gamma\\\\[1.00006pt]
\gamma\displaystyle\int_{\Gamma_{I}}(\mathbf{N}^{2})^{\textrm{T}}\bm{n}\bm{\kappa}^{1}\mathbf{B}^{1}~{}d\Gamma&(1-\gamma)\displaystyle\int_{\Gamma_{I}}(\mathbf{N}^{2})^{\textrm{T}}\bm{n}\bm{\kappa}^{2}\mathbf{B}^{2}~{}d\Gamma\end{bmatrix}$
(16)
$\mathbf{K}^{s}=\begin{bmatrix}\beta\displaystyle\int_{\Gamma_{I}}(\mathbf{N}^{1})^{\textrm{T}}\mathbf{N}^{1}~{}d\Gamma&-\beta\displaystyle\int_{\Gamma_{I}}(\mathbf{N}^{1})^{\textrm{T}}\mathbf{N}^{2}~{}d\Gamma\\\\[1.00006pt]
-\beta\displaystyle\int_{\Gamma_{I}}(\mathbf{N}^{2})^{\textrm{T}}\mathbf{N}^{1}~{}d\Gamma&\beta\displaystyle\int_{\Gamma_{I}}(\mathbf{N}^{2})^{\textrm{T}}\mathbf{N}^{2}~{}d\Gamma\end{bmatrix}$
(17)

Matrices $\bm{\kappa}^{m}$, $\mathbf{B}^{m}$ and $\mathbf{N}^{m}$ in Equation
(15)-(17) are the conductivity matrix, shape function derivative matrix and
shape function vector respectively (for a given patch $m$). The matrices
$\mathbf{B}^{m}$ and $\mathbf{N}^{m}$ are given as below,

$\mathbf{B}^{m}=\begin{bmatrix}{N_{1,x}^{m}}&{N_{2,x}^{m}}&...&{N_{I,x}^{m}}&...\\\\[1.00006pt]
{N_{1,y}^{m}}&{N_{2,y}^{m}}&...&{N_{I,y}^{m}}&...\end{bmatrix},\quad\mathbf{N}^{m}=\begin{bmatrix}{N_{1}^{m}}&{N_{2}^{m}}&...&{N_{I}^{m}}&...\end{bmatrix}.$
(18)

The verification is performed for the first example to check the accuracy of
Nitsche’s method (refer Section 4.4.1). For the current work, the
stabilization parameter $\beta=1\times 10^{12}$. With the large value of the
stabilization parameter, Nitsche’s method behaves close to the penalty method.
Here the objective is to force the temperature continuity in a stricter sense
and normal flux continuity in a weaker sense. In the literature [30, 31], it
is also reported that the large stabilization parameter might cause ill-
conditioning of the system, but we did not face any conditioning issue for our
boundary value problem. For the averaging parameter, a standard value
$\gamma=0.5$ is taken, which gives equal weights to fluxes on both sides of
the interface.

##  3 Optimization

###  3.1 Optimization problem description

In a standard shape optimization problem, the design variables are the
parameters that control the geometry and try to optimize a function called
objective function. In our case, suppose $\bm{x}$ is the vector of the
$N_{\mathrm{var}}$ design variables and $f_{\mathrm{obj}}$ is the objective
function, then the shape optimization problem for a heat manipulator can be
defined in a mathematical form as,

$\min_{\bm{x}\in\mathbb{R}^{N_{\mathrm{var}}}}~{}f_{\mathrm{obj}}(\bm{x}),$
(20)

with

$\displaystyle
f_{\mathrm{obj}}:\mathbb{R}^{N_{\mathrm{var}}}\rightarrow\mathbb{R},$ (21a)
$\displaystyle f_{\mathrm{obj}}:\bm{x}\mapsto f_{\mathrm{obj}}(\bm{x}),$ (21b)

such that the following constraints are satisfied,

Equality constraints: $\displaystyle h_{i}(\bm{x})=~{}0\quad$ $\displaystyle
i=1,2,...,n_{h},$ (22) Inequality constraints: $\displaystyle
g_{j}(\bm{x})\leq~{}0\quad$ $\displaystyle j=1,2,...,n_{g},$ (23) Box
constraints: $\displaystyle
x_{k,\rm{min}}~{}\leq~{}x_{k}\leq~{}x_{k,\rm{max}}\quad$ $\displaystyle
k=1,2,...,N_{\mathrm{var}},$ (24)

where $n_{h}$ and $n_{g}$ are the number of equality constraints and
inequality constraints respectively. $x_{k,\rm{min}}$ and $x_{k,\rm{max}}$ are
lower and upper bounds of the design variable $x_{k}$.

###  3.2 Particle swarm optimization

Kennedy and Eberhart, in their article [59], proposed Particle Swarm
Optimization (PSO), an optimization concept based on the swarming behaviour of
birds flock and fish school. The PSO is a non-gradient based optimization
method and can search very large spaces of candidate solutions for the
optimization of continuous nonlinear functions.

In this method, each candidate solution is considered as a particle. The
algorithm starts with a swarm of initial particles, then advances towards the
optimum solution by updating the positions and velocities of these particles
according to their fitness to objective function. In our minimization-
optimization problem, the position of a particle is a vector of design
variable values, and the velocity is the correction applied to the position at
the each iteration. Furthermore, the fitness measures how small the objective
function is for a given particle. The position and velocity of a particle are
updated based on its individual best position as well as the position of the
best particle (the particle with the smallest objective function).

Suppose, $\omega$ is a non-negative number which defines the contribution of
old velocity on the updated velocity; $\mathrm{pbest}_{k}^{i}$ is the best
position of the $k^{th}$ particle at the $i^{th}$-iteration, which has the
smallest value of objective function among all positions utilized before by
the $k^{th}$ particle; $\mathrm{gbest}^{i}$ is the position of the best
particle that has the smallest objective function up to the current iteration
from the whole swarm; $\alpha_{1}$ and $\alpha_{2}$ are the acceleration
coefficients and can be indicated as self adjustment weight and social
adjustment weight respectively. The self adjustment weight defines the maximum
correction in the direction of the individual best position, while social
adjustment weight defines the maximum correction in the direction of the
position of best particle; $\mathrm{rand}_{1}$ and $\mathrm{rand}_{2}$ are two
random numbers between 0 and 1. Then, by using the velocity and position at
the $i^{th}$-iteration, the equations to find the updated velocity and
position of the $k^{th}$ particle at the $(i+1)^{th}$-iteration can be written
as,

$v_{k}^{i+1}=\omega
v_{k}^{i}+\alpha_{1}~{}\mathrm{rand}_{1}~{}(\mathrm{pbest}_{k}^{i}-x_{k}^{i})\\\
+\alpha_{2}~{}\mathrm{rand}_{2}~{}(\mathrm{gbest}^{i}-x_{k}^{i}),$ (25a)
$x_{k}^{i+1}=x_{k}^{i}+v_{k}^{i+1}.$ (25b)

The detailed methodology of the PSO is presented in the form of a flow chart
in Figure 4. At first, we need to define the values of the parameters utilized
in the Equation (25b) such as $\alpha_{1}$, $\alpha_{2}$, and $\omega$. The
next step is to initialize the position $\bm{x}$, such that it satisfies the
equality and inequality constraints, and velocity $\bm{v}$. Then, the geometry
and mathematical model are constructed by using $\bm{x}$. Moreover, all
constraints are being checked. For box constraints, if any $x_{k}$ is out of a
bound defined by Equation (24), it is set equal to the bound. However, if
equality constraints or inequality constraints are not satisfied, then
$f_{\mathrm{obj}}$ is set to a penalty number to avoid that specific particle
as a solution.

As we proceed further, using the geometry and mathematical model, the boundary
value problem (In our case, heat conduction problem as described in section
3.1) is being solved and corresponding objective function $f_{\mathrm{obj}}$
is evaluated. Afterwards, $\mathrm{pbest}$ for each particle and
$\mathrm{gbest}$ for the whole swarm are calculated based on
$f_{\mathrm{obj}}$ values of all particles. At the end of each iteration, the
stopping criteria are checked. If the criteria are not satisfied then the
velocity and position of each particle are updated by Equation (25b) and the
procedure will be repeated with the new values. This optimization loop
continues until the stopping criteria are satisfied. Once, the stopping
criteria are satisfied, the loop ends and $\mathrm{gbest}$ is taken as the
optimum solution.

Start Define the values of parameters $\alpha_{1}$, $\alpha_{2}$, $\omega$,
parameters related to constraints and stopping criteria Initialize the design
variables $\bm{x}$ satisfying the equality and inequality constraints, and
velocity $\bm{v}$ Are the box constraints satisfied?
$x_{k,\rm{min}}~{}\leq~{}x_{k}\leq~{}x_{k,\rm{max}}$; $k=1,2,...,n_{x}$ Create
a geometry using $\bm{x}$, then prepare a mathematical model Are the equality
and inequality constraints satisfied? Solve heat conduction boundary value
problem (as described in section ) using IGA Evaluate the objective function
$f_{\mathrm{obj}}(\bm{x})$ Find $\mathrm{pbest}_{k}^{i}$ and
$\mathrm{gbest}^{i}$ Update position and velocity $\displaystyle v_{k}^{i+1}$
$\displaystyle=\omega
v_{k}^{i}+\alpha_{1}~{}\mathrm{rand}_{1}~{}(\mathrm{pbest}_{k}^{i}-x_{k}^{i})$
$\displaystyle\quad+\alpha_{2}~{}\mathrm{rand}_{2}~{}(\mathrm{gbest}^{i}-x_{k}^{i})$
$\displaystyle x_{k}^{i+1}$ $\displaystyle=x_{k}^{i}+v_{k}^{i+1}$ Is the
stopping conditions satisfied? $\mathrm{gbest}^{i}$ is the optimal solution
Stop $i=i+1$ $x_{k}=x_{k,\rm{min}},\quad\text{if}~{}x_{k}<x_{k,\rm{min}}$
$x_{k}=x_{k,\rm{max}}\quad\text{if}~{}x_{k}>x_{k,\rm{max}}$
$f_{\mathrm{obj}}(\bm{x})$= penalty value YESNOYESYESNONO

Figure 4: Flow chart of particle swarm optimization algorithm

In the present work, we utilize the inbuilt MATLAB function “particleswarm” to
implement PSO algorithm. We consider tolerance of $1\times 10^{-6}$ in the
change of objective function $f_{\mathrm{obj}}$ value and 15 consecutive
stalled iterations (improvement in $f_{\mathrm{obj}}$ value is less than
$1\times 10^{-6}$) as stopping criteria. We take MATLAB default values of
$\alpha_{1}$, $\alpha_{1}$=1.49 and $\alpha_{2}$, $\alpha_{2}$=1.49 for our
case as well. Furthermore, the value of $\omega$ is adaptive between [0.1,1.1]
as per the inbuilt scheme. The particular value of $\omega$ is decreased or
increased based on stalled iteration count.

##  4 Thermal Concentrator

###  4.1 Problem definition

(a) A base material plate under application of constant flux.

(b) A thermal concentrator embedded in the plate.

(c) A metamaterial based concentrator in the base material plate.

Figure 5: Schematic design of (a) The base material domain $\Omega$ (nickel
steel plate) under constant heat flux applied by high temperature source on
the left side and low temperature sink on the right side;
$\Omega_{\mathrm{in}}$ is the region of focus where we want to concentrate the
flux, $\Gamma_{\mathrm{flux}}$ is the mid-section along which we measure the
flux (b) A thermal concentrator embedded in the base material plate to
increase the flux concentration in $\Omega_{\mathrm{in}}$,
$\Omega_{\mathrm{design}}$ is the area of the concentrator where we optimize
the shape, $\Omega_{\mathrm{out}}$ is the outside area of remaining base
material,
$\Omega=\Omega_{\mathrm{in}}\cup\Omega_{\mathrm{design}}\cup\Omega_{\mathrm{out}}$
(c) An embedded matermaterial-based thermal concentrator created by putting
alternative sectors of copper and PDMS.

Figure 6: NURBS parameterization of the domain $\Omega$. The detailed view of
connectivity between adjacent patches as well as the design control points
which are the means to manipulate the interfaces between copper and PDMS
sectors, and eventually the shape of the concentrator.

In this example, we optimize a thermal concentrator, whose main objective is
to concentrate the heat flux in the specified region of the domain. The
schematics of the problem is given in Figure 5. The geometry and problem
details are referred from [32]. We consider a square base material-nickel
steel plate ($\kappa_{\rm ni}=10$ W/mK) with side length 140 mm. To increase
the flux concentration in the circular core region of 10 mm radius at center,
an annular sector-shaped concentrator of outer radius 50 mm is introduced
around the core. The geometry is constructed by placing, side by side, annulus
sectors made of two materials, one with higher conductivity and another with
lower conductivity. The number of sectors employed in the geometry is denoted
by $N_{\mathrm{sec}}$. To ascertain that each sector is adjacent to the sector
of other material from both side, $N_{\mathrm{sec}}$ is always taken as an
even number, $N_{\mathrm{sec}}=2\ell$ with $\ell=2,3,4,...$. Now, the sectors
are arranged such that the center axis of one of sectors is overlaps with the
$x$-axis. Moreover, the center axes of all sectors are uniformly distributed
angle-wise on the annular region. Altogether, the concentrator appeared as a
uniform series connection of two alternating materials in the azimuthal
direction. Hence, the effective medium behaves as an anisotropic medium and
alters the heat flux direction to follow a particular path of interest. In our
case, copper and polydimethylsiloxane (PDMS), with thermal conductivity
$\kappa_{\rm copper}=398$ W/mK and $\kappa_{\rm PDMS}=0.27$ W/mK, are
exploited to design the concentrator. Initially, the core material is taken as
same as background base material nickel steel. The left and right side of the
plate are fixed at $300$ K and $200$ K temperature respectively.

###  4.2 Objective function

The objective of the concentrator is to concentrate the heat flux in inner
core $\Omega_{\mathrm{in}}$. The concentrated flux is measured along the mid-
section $\Gamma_{\mathrm{flux}}$ of $\Omega_{\mathrm{in}}$. Mathematically,
the concentrator function can be written as,

$\Psi_{\mathrm{flux}}=\dfrac{1}{\overline{\Psi}_{\mathrm{flux}}}\int_{\Gamma_{\mathrm{flux}}}-\kappa_{\mathrm{in}}\nabla
T\cdot\bm{n}~{}d\Gamma,$ (26)

where $\kappa_{\mathrm{in}}$ is the conductivity of core, and
$\overline{\Psi}_{\mathrm{flux}}$ is the normalisation value given as,

$\overline{\Psi}_{\mathrm{flux}}=\int_{\Gamma_{\mathrm{flux}}}-\overline{\kappa}_{\mathrm{in}}\nabla\overline{T}\cdot\mathbf{n}~{}d\Gamma,$
(27)

Here, all the properties with overline represent the reference case when whole
domain is filled with base material.

Now, since we are solving a minimization problem for optimization, we define
the objective function as

$f_{\mathrm{obj}}=\left|\dfrac{1}{\Psi_{\mathrm{flux}}}\right|.$ (28)

###  4.3 NURBS parameterization, design variables and constraints

For the NURBS approximation of the geometry, each sector, inner core and
outside base material are considered as a separate patch as shown in Figure 6.
As mentioned in Section 2, Nitsche’s method is implemented to provide
continuity of temperature and normal flux across patch interfaces. All
external and internal boundaries are kept unchanged during optimization,
except the interfaces between each pair of sectors (shown as white curves in
the detailed view in Figure 6). In fact, only some of the control points of
these interfaces are used for shape manipulation, which are called design
control points (shown as white dots in the detailed view in Figure 6). The
knot vectors and the weight of control points are kept fixed as well.

In the present work, the radial positions of design control points are fixed
and distributed uniformly between the inner and outer radius. And the
circumferential positions on the predefined circular paths by radial positions
are interpreted as design variables. The number of design variables is denoted
by $N_{\mathrm{var}}$. To reduce $N_{\mathrm{var}}$ and subsequent
computational burden on the optimization problem, symmetry along center axis
is imposed for each sector. Therefore during optimization, the interface
between only one pair of sectors is designed. All other interfaces will follow
through the imposed symmetry. Now, due to overlap of the center axis of one
sector on $x$-axis as mentioned earlier, the symmetry along $x$-axis is also
preserved. The position of design control points (calculated using design
variables), the predefined knot vector of the interface and the imposed
symmetry will decide the overall shape of sector patches.

We consider the sector, whose center lies along the negative $x$-axis, as the
first sector and assigned it the material with higher conductivity, copper.
Then, each sector is assigned a material alternatively as shown in Figure 6.
Now, to simplify the explanation of the results in the next section, the
possible configurations are divided into four types A,B,C and D according to
$N_{\mathrm{sec}}$ as shown in Table 1. The axes of symmetry present in each
configuration are also mentioned in the table.

Table 1: Different possible configurations of the concentrator due to imposed
symmetry and $N_{\mathrm{sec}}$.

Type |  A |  B |  C |  D
---|---|---|---|---
$N_{\mathrm{sec}}$ |  $8\ell-4$ |  $8\ell-2$ |  $8\ell$ |  $8\ell+2$
Axes of symmetry |  $x$-axis, $y$-axis and sector center axis |  $x$-axis and sector center axis |  $x$-axis, $y$-axis and sector center axis |  $x$-axis and sector center axis
Sector center-axes along $x$-axis |  copper sectors |  one copper and one PDMS sector |  copper sectors |  one copper and one PDMS sector
Sector center-axes along $x$-axis |  PDMS sectors |  - |  copper sectors |  -
Schematic |  |  |  |

For the optimization, box constraints and inequality constraints are applied.
The box constraints restrict the position of design control points between two
center-axes, and avoid the singularities occurring due to the intersection of
two interfaces. An inequality constraint ensuring the minimum area of each
NURBS patch is imposed too. This constraint helps to evade very thin patches
and the numerical instability created by them. In addition, another inequality
constraint in terms of non-negativity of the Jacobain of NURBS
parameterization is imposed, which guarantees only physically feasible
geometry (avoid self-overlap or self-intersection).

###  4.4 Results and discussion

####  4.4.1 Verification of Nitsche’s method

(a) Temperature $T$ distribution with $\llbracket T\rrbracket=5.318\times
10^{-10}$ at $\Gamma_{I}$

(b) Flux flow with $\dfrac{\bm{n}\cdot\llbracket\bm{\kappa}\nabla
T\rrbracket}{|\bm{n}\cdot\bm{\kappa}^{1}\nabla T^{1}|}=4.306\times 10^{-3}$ at
$\Gamma_{I}$

Figure 7: Verification of applied Nitsche’s method through checking the
continuity conditions after boundary value problem solution of initial
geometry with $N_{\mathrm{sec}}=4$.(a) Temperature distribution with total
jump in temperature value $\llbracket T\rrbracket=5.318\times 10^{-10}$ at
$\Gamma_{I}$ (b) Flux flow with total jump in the relative normal flux value
$\dfrac{\bm{n}\cdot\llbracket\bm{\kappa}\nabla
T\rrbracket}{|\bm{n}\cdot\bm{\kappa}^{1}\nabla T^{1}|}=4.306\times 10^{-3}$ at
$\Gamma_{I}$.

In this section, we verified the Nitsch’e method utilized to apply the
interface continuity condition. We consider the initial configuration (used
for optimization) of the concentrator with straight radial edges of sectors
with $N_{\mathrm{sec}}=4$. The boundary value problem is solved and the
results are presented in Figure 7. The temperature along the mid-section
(along the $x$-axis) is shown in blue. From Figure 7(a), it is evident that
the temperature is continuous across all patch interfaces. Now, from Figure
7(b), we can see that the flux streamlines flow from one patch to another
without any jump. The numerical values of total jump of temperature and
relative normal flux (considering all the interfaces) are measured
$5.318\times 10^{-10}$ and $4.306\times 10^{-3}$. From these results, we can
verify that the Nitsche’s method impose the required continuity condition.

####  4.4.2 A convergence study

$10^{2}$$10^{3}$$10^{4}$$10^{5}$$10^{-2}$$10^{-1}$Degrees of freedomRelative
error in $f_{obj}$$N_{\mathrm{sec}}$=4, S1$N_{\mathrm{sec}}$=12,
S1$N_{\mathrm{sec}}$=4, S2$N_{\mathrm{sec}}$=12, S2$N_{\mathrm{sec}}$=4,
S3$N_{\mathrm{sec}}$=12, S3

(a) Type A, $N_{\mathrm{var}}=3$

$10^{3}$$10^{4}$$10^{5}$$10^{-3}$$10^{-2}$Degrees of freedomRelative error in
$f_{obj}$$N_{\mathrm{sec}}$=4, S1$N_{\mathrm{sec}}$=12,
S1$N_{\mathrm{sec}}$=4, S2$N_{\mathrm{sec}}$=12, S2$N_{\mathrm{sec}}$=4,
S3$N_{\mathrm{sec}}$=12, S3

(b) Type A, $N_{\mathrm{var}}=4$

$10^{3}$$10^{4}$$10^{5}$$10^{-2}$Degrees of freedomRelative error in
$f_{obj}$$N_{\mathrm{sec}}$=6, S1$N_{\mathrm{sec}}$=14,
S1$N_{\mathrm{sec}}$=6, S2$N_{\mathrm{sec}}$=14, S2$N_{\mathrm{sec}}$=6,
S3$N_{\mathrm{sec}}$=14, S3

(c) Type B, $N_{\mathrm{var}}=3$

$10^{3}$$10^{4}$$10^{5}$$10^{-3}$$10^{-2}$Degrees of freedomRelative error in
$f_{obj}$$N_{\mathrm{sec}}$=6, S1$N_{\mathrm{sec}}$=14,
S1$N_{\mathrm{sec}}$=6, S2$N_{\mathrm{sec}}$=14, S2$N_{\mathrm{sec}}$=6,
S3$N_{\mathrm{sec}}$=14, S3

(d) Type B, $N_{\mathrm{var}}=4$

$10^{3}$$10^{4}$$10^{5}$$10^{-3}$Degrees of freedomRelative error in
$f_{obj}$$N_{\mathrm{sec}}$=8, S1$N_{\mathrm{sec}}$=16,
S1$N_{\mathrm{sec}}$=8, S2$N_{\mathrm{sec}}$=16, S2$N_{\mathrm{sec}}$=8,
S3$N_{\mathrm{sec}}$=16, S3

(e) Type C, $N_{\mathrm{var}}=3$

$10^{3}$$10^{4}$$10^{5}$$10^{-4}$$10^{-3}$$10^{-2}$Degrees of freedomRelative
error in $f_{obj}$$N_{\mathrm{sec}}$=8, S1$N_{\mathrm{sec}}$=16,
S1$N_{\mathrm{sec}}$=8, S2$N_{\mathrm{sec}}$=16, S2$N_{\mathrm{sec}}$=8,
S3$N_{\mathrm{sec}}$=16, S3

(f) Type C, $N_{\mathrm{var}}=4$

$10^{3}$$10^{4}$$10^{5}$$10^{-2}$$10^{-1}$Degrees of freedomRelative error in
$f_{obj}$$N_{\mathrm{sec}}$=10, S1$N_{\mathrm{sec}}$=18,
S1$N_{\mathrm{sec}}$=10, S2$N_{\mathrm{sec}}$=18, S2$N_{\mathrm{sec}}$=10,
S3$N_{\mathrm{sec}}$=18, S3

(g) Type D, $N_{\mathrm{var}}=3$

$10^{3}$$10^{4}$$10^{5}$$10^{-2}$Degrees of freedomRelative error in
$f_{obj}$$N_{\mathrm{sec}}$=10, S1$N_{\mathrm{sec}}$=18,
S1$N_{\mathrm{sec}}$=10, S2$N_{\mathrm{sec}}$=18, S2$N_{\mathrm{sec}}$=10,
S3$N_{\mathrm{sec}}$=18, S3

(h) Type D, $N_{\mathrm{var}}=4$

Figure 8: Convergence of the objective function $f_{\mathrm{obj}}$ for thermal
concentrator problem for Type-A ($N_{\mathrm{sec}}=$ 4, 12), Type-B
($N_{\mathrm{sec}}=$ 6, 14), Type-C ($N_{\mathrm{sec}}=$ 8, 16), Type-D
($N_{\mathrm{sec}}=$ 10, 18) configurations, and $N_{\mathrm{var}}=$ 3, 4.
$4\%$ relative error (shown by pink horizontal line) is taken as the tolerance
value to select the minimum mesh size for optimization. Samples 1, 2, and 3
are denoted by S1, S2 and S3. Type-A requires a finer mesh compared to Type-B
and Type-D. Type-C needs a comparatively coarse mesh.

$10^{3}$$10^{4}$$10^{-2}$$10^{-1}$Degrees of freedomRelative error in
$f_{obj}$$N_{\mathrm{sec}}$=4, S1$N_{\mathrm{sec}}$=12,
S1$N_{\mathrm{sec}}$=4, S2$N_{\mathrm{sec}}$=12, S2$N_{\mathrm{sec}}$=4,
S3$N_{\mathrm{sec}}$=12, S3

(a) Type A, $N_{\mathrm{var}}=5$

$10^{3}$$10^{4}$$10^{5}$$10^{-3}$$10^{-2}$Degrees of freedomRelative error in
$f_{obj}$$N_{\mathrm{sec}}$=4, S1$N_{\mathrm{sec}}$=12,
S1$N_{\mathrm{sec}}$=4, S2$N_{\mathrm{sec}}$=12, S2$N_{\mathrm{sec}}$=4,
S3$N_{\mathrm{sec}}$=12, S3

(b) Type A, $N_{\mathrm{var}}=6$

$10^{3}$$10^{4}$$10^{5}$$10^{-2}$Degrees of freedomRelative error in
$f_{obj}$$N_{\mathrm{sec}}$=6, S1$N_{\mathrm{sec}}$=14,
S1$N_{\mathrm{sec}}$=6, S2$N_{\mathrm{sec}}$=14, S2$N_{\mathrm{sec}}$=6,
S3$N_{\mathrm{sec}}$=14, S3

(c) Type B, $N_{\mathrm{var}}=5$

$10^{3}$$10^{4}$$10^{5}$$10^{-3}$$10^{-2}$Degrees of freedomRelative error in
$f_{obj}$$N_{\mathrm{sec}}$=6, S1$N_{\mathrm{sec}}$=14,
S1$N_{\mathrm{sec}}$=6, S2$N_{\mathrm{sec}}$=14, S2$N_{\mathrm{sec}}$=6,
S3$N_{\mathrm{sec}}$=14, S3

(d) Type B, $N_{\mathrm{var}}=6$

$10^{3}$$10^{4}$$10^{5}$$10^{-4}$$10^{-3}$Degrees of freedomRelative error in
$f_{obj}$$N_{\mathrm{sec}}$=8, S1$N_{\mathrm{sec}}$=16,
S1$N_{\mathrm{sec}}$=8, S2$N_{\mathrm{sec}}$=16, S2$N_{\mathrm{sec}}$=8,
S3$N_{\mathrm{sec}}$=16, S3

(e) Type C, $N_{\mathrm{var}}=5$

$10^{4}$$10^{5}$$10^{-4}$$10^{-3}$Degrees of freedomRelative error in
$f_{obj}$$N_{\mathrm{sec}}$=8, S1$N_{\mathrm{sec}}$=16,
S1$N_{\mathrm{sec}}$=8, S2$N_{\mathrm{sec}}$=16, S2$N_{\mathrm{sec}}$=8,
S3$N_{\mathrm{sec}}$=16, S3

(f) Type C, $N_{\mathrm{var}}=6$

$10^{4}$$10^{5}$$10^{-2}$Degrees of freedomRelative error in
$f_{obj}$$N_{\mathrm{sec}}$=10, S1$N_{\mathrm{sec}}$=18,
S1$N_{\mathrm{sec}}$=10, S2$N_{\mathrm{sec}}$=18, S2$N_{\mathrm{sec}}$=10,
S3$N_{\mathrm{sec}}$=18, S3

(g) Type D, $N_{\mathrm{var}}=5$

$10^{4}$$10^{5}$$10^{-3}$$10^{-2}$Degrees of freedomRelative error in
$f_{obj}$$N_{\mathrm{sec}}$=10, S1$N_{\mathrm{sec}}$=18,
S1$N_{\mathrm{sec}}$=10, S2$N_{\mathrm{sec}}$=18, S2$N_{\mathrm{sec}}$=10,
S3$N_{\mathrm{sec}}$=18, S3

(h) Type D, $N_{\mathrm{var}}=6$

Figure 9: Convergence of the objective function $f_{\mathrm{obj}}$ for thermal
concentrator problem for Type-A ($N_{\mathrm{sec}}=$ 4, 12), Type-B
($N_{\mathrm{sec}}=$ 6, 14), Type-C ($N_{\mathrm{sec}}=$ 8, 16), Type-D
($N_{\mathrm{sec}}=$ 10, 18) configurations, and $N_{\mathrm{var}}=$ 5, 6.
$4\%$ relative error (shown by pink horizontal line) is taken as the tolerance
value to select the minimum mesh size for optimization. Samples 1, 2, and 3
are denoted by S1, S2 and S3. Type-A requires a finer mesh compared to Type-B
and Type-D. Type-C needs a comparatively coarse mesh.

At first, we study the convergence of the objective function with refinement
of mesh. In this study, the optimization (with chosen $N_{\mathrm{sec}}$ and
$N_{\mathrm{var}}$) is run with a coarse mesh to find the optimized shape.
Then, the mesh is refined without changing the shape to investigate the effect
of refinement on the objective function value. This study assists us in
finding a minimum mesh size that can ensure sufficient accuracy for the
objective function during optimization. In Figures 8-9, the relative error in
the objective function value with respect to the degrees of freedom is plotted
for several numbers of sectors as well as for numbers of design variables. The
relative error is measured with reference to the most refined mesh for each
case. $N_{\mathrm{sec}}$= $4,6,8,10,12,14,16$ and 18, while
$N_{\mathrm{var}}$= $3,4,5$ and 6. For $N_{\mathrm{var}}=3,5$, $C^{0}$
continuity is applied between the interface NURBS elements, while for
$N_{\mathrm{var}}=4,6$, $C^{1}$ continuity is applied. We assume that a mesh
with a relative error within $4\%$ is sufficiently fine and will subsequently
use this discretization for optimization.

From Figures 8-9, it is apparent that the Type-A configuration requires a
finer mesh for convergence compared to Type-B and Type-D. On the other hand,
Type-C needs a comparatively coarse mesh for convergence. Overall, the
relative error tolerance can be reached by a mesh with less than $10^{4}$
degrees of freedom.

####  4.4.3 Shape optimization

(a) $N_{\mathrm{var}}=3$

(b) $N_{\mathrm{var}}=4$

(c) $N_{\mathrm{var}}=5$

(d) $N_{\mathrm{var}}=6$

Figure 10: Variation of concentrator function $\Psi_{\mathrm{flux}}$ of the
optimized shape with respect to number of sectors for different number of
design variables. $N_{\mathrm{var}}=$3, 4, 5, 6. Samples 1, 2, and 3 are
denoted by S1, S2 and S3. The un-optimized geometry concentrates around 3
times more flux than a base material flat plate, while the optimized
geometries concentrates up to 6 times more flux. Type-A configuration
concentrates more flux compared to other configurations. There is no
significant gain by using larger $N_{\mathrm{var}}$.

(a) Type-A, $N_{\mathrm{sec}}=4$, $N_{\mathrm{var}}=3$

(b) Type-A, $N_{\mathrm{sec}}=4$, $N_{\mathrm{var}}=4$

(c) Type-A, $N_{\mathrm{sec}}=4$, $N_{\mathrm{var}}=5$

(d) Type-A, $N_{\mathrm{sec}}=4$, $N_{\mathrm{var}}=6$

(e) Type-B, $N_{\mathrm{sec}}=6$, $N_{\mathrm{var}}=3$

(f) Type-B, $N_{\mathrm{sec}}=6$, $N_{\mathrm{var}}=4$

(g) Type-B, $N_{\mathrm{sec}}=6$, $N_{\mathrm{var}}=5$

(h) Type-B, $N_{\mathrm{sec}}=6$, $N_{\mathrm{var}}=6$

(i) Type-C, $N_{\mathrm{sec}}=8$, $N_{\mathrm{var}}=3$

(j) Type-C, $N_{\mathrm{sec}}=8$, $N_{\mathrm{var}}=4$

(k) Type-C, $N_{\mathrm{sec}}=8$, $N_{\mathrm{var}}=5$

(l) Type-C, $N_{\mathrm{sec}}=8$, $N_{\mathrm{var}}=6$

(m) Type-D, $N_{\mathrm{sec}}=10$, $N_{\mathrm{var}}=3$

(n) Type-D, $N_{\mathrm{sec}}=10$, $N_{\mathrm{var}}=4$

(o) Type-D, $N_{\mathrm{sec}}=10$, $N_{\mathrm{var}}=5$

(p) Type-D, $N_{\mathrm{sec}}=10$, $N_{\mathrm{var}}=6$

Figure 11: Flux flow in the optimized shape for different number of sectors
and number of design variables. $N_{\mathrm{sec}}=$4 (Type-A), 6 (Type-B), 8
(Type-C), 10 (Type-D) and $N_{\mathrm{var}}=$3, 4, 5, 6. Type-A configuration
concentrates more flux compared to other configurations. Increasing
$N_{\mathrm{var}}$ only brush up the details of the optimized interfaces.

Figure 10 and Table 2 show the concentrator function $\Psi_{\mathrm{flux}}$
variation with respect to number of sectors for $N_{\mathrm{var}}=3,4,5,6$.
Here, we show $\Psi_{\mathrm{flux}}$-value for a geometry with straight lines
as interfaces (i.e. geometry without any optimization, same as utilized in
[32]) as well as the optimized geometries obtained by three sample
optimization runs. For each sample, PSO algorithm goes through and compares
several solutions (completely different from another sample but in the design
space) created by random numbers before eventually leading to an optimized
shape at the end. We can consider each sample as a separate path from initial
guess to final optimized solution. As can be seen from Figure 10,
$f_{\mathrm{obj}}$ is almost constant for all three runs of optimization,
which corroborates the reproducibility of the present optimization process and
sufficiency of the exploited mesh for calculating $f_{\mathrm{obj}}$
accurately. Figure 10 highlights the benefit achieved in
$\Psi_{\mathrm{flux}}$-value with optimization. We can notice that Type-A
configuration has bigger $\Psi_{\mathrm{flux}}$-value compared to all other
configurations. In particular, $N_{\mathrm{sec}}=4$ performs better than any
other number of sectors. For $N_{\mathrm{sec}}=4$, the un-optimized geometry
concentrates around 3 times more flux than a base material flat plate, while
the optimized geometry concentrates approximately 6 times more flux (refer
Table 2).

Now, we analyse the effect of $N_{\mathrm{var}}$ on objective function. To
understand the effect of $N_{\mathrm{var}}$, streamline plots are plotted for
Type A to D configuration ($N_{\mathrm{sec}}=4,6,8,10$) in Figure 11. From the
figure, we can observe that the optimized shapes are more or less identical
for different $N_{\mathrm{var}}$, and there is no significant gain by using
larger $N_{\mathrm{var}}$. As large $N_{\mathrm{var}}$ is the indication of
more degrees of freedom to manipulate the shape, increasing $N_{\mathrm{var}}$
only brushes up the details of the optimized interface. Especially
$N_{\mathrm{var}}=4,6$, which have $C^{1}$ continuity between interface NURBS
elements, smoothens the interface between the sectors. Nevertheless, the
improvement provided in the objective function value is very limited.

####  4.4.4 Shape optimization combined with conductivity optimization

(a) $N_{\mathrm{var}}=3$

(b) $N_{\mathrm{var}}=4$

(c) $N_{\mathrm{var}}=5$

(d) $N_{\mathrm{var}}=6$

Figure 12: Variation of concentrator function $\Psi_{\mathrm{flux}}$ of the
optimized shape with respect to number of sectors for three case (I) without
any optimization (II) only shape optimization, and (III) shape optimization
combined with conductivity optimization. $N_{\mathrm{var}}=$3, 4, 5, 6. For
case-III, The concentrated flux is approximately 10 to 13 times as large as
for a base material plate, 2 to 3 times as large as for case-II & 3 to 6 times
as large as for case-I.

(a) Type-A, $N_{\mathrm{sec}}=4$, shape opt

(b) Type-B, $N_{\mathrm{sec}}=6$, shape opt

(c) Type-C, $N_{\mathrm{sec}}=8$, shape opt

(d) Type-D, $N_{\mathrm{sec}}=10$, shape opt

(e) Type-A, $N_{\mathrm{sec}}=4$, shape + $\kappa$ opt

(f) Type-B, $N_{\mathrm{sec}}=6$, shape + $\kappa$ opt

(g) Type-C, $N_{\mathrm{sec}}=8$, shape + $\kappa$ opt

(h) Type-D, $N_{\mathrm{sec}}=10$, shape + $\kappa$ opt

(i) Type-A, $N_{\mathrm{sec}}=12$, shape opt

(j) Type-B, $N_{\mathrm{sec}}=14$, shape opt

(k) Type-C, $N_{\mathrm{sec}}=16$, shape opt

(l) Type-D, $N_{\mathrm{sec}}=18$, shape opt

(m) Type-A, $N_{\mathrm{sec}}=12$, shape + $\kappa$ opt

(n) Type-B, $N_{\mathrm{sec}}=14$, shape + $\kappa$ opt

(o) Type-C, $N_{\mathrm{sec}}=16$, shape + $\kappa$ opt

(p) Type-D, $N_{\mathrm{sec}}=18$, shape + $\kappa$ opt

Figure 13: Flux flow in the optimized shape for different number of sectors for two cases (II) only shape optimization, and (III) shape + conductivity optimization. $N_{\mathrm{sec}}=$4, 12 (Type-A), 6, 14 (Type-B), 8, 16 (Type-C), 10, 20 (Type-D) and $N_{\mathrm{var}}=$6\. Base material is shown in grey color. Copper in the core and PDMS sectors on $y$-axis (inline with $\Gamma_{\mathrm{flux}}$) provide benefits in concentrating more flux. Table 2: Objective function $f_{\mathrm{obj}}$ and concentrator function $\Psi_{\mathrm{flux}}$ values for case-I: without any optimization, case-II: shape optimization, and case-III: shape optimization combined with conductivity optimization. $N_{\mathrm{sec}}=$ 4, 6, 8, 10, 12, 14, 16, 18, 20. $N_{\mathrm{var}}=$ 3, 4, 5, 6. Number of sectors |  Number of design variables | w/o opt | Shape opt |  Shape + $\kappa$ opt
---|---|---|---|---
| $f_{\mathrm{obj}}$ |  $\Psi_{\mathrm{flux}}$ | $f_{\mathrm{obj}}$ |  $\Psi_{\mathrm{flux}}$ | $f_{\mathrm{obj}}$ |  $\Psi_{\mathrm{flux}}$
$N_{\mathrm{sec}}=4$ |  $N_{\mathrm{var}}$=3 |  2.391$\times 10^{-1}$ |  4.182 |  1.407$\times 10^{-1}$ |  7.105 |  7.138$\times 10^{-2}$ |  14.010
$N_{\mathrm{var}}$=4 |  1.421$\times 10^{-1}$ |  7.036 |  7.148$\times 10^{-2}$ |  13.991
$N_{\mathrm{var}}$=5 |  1.433$\times 10^{-1}$ |  6.977 |  7.143$\times 10^{-2}$ |  14.001
$N_{\mathrm{var}}$=6 |  1.407$\times 10^{-1}$ |  7.106 |  7.136$\times 10^{-2}$ |  14.013
$N_{\mathrm{sec}}=6$ |  $N_{\mathrm{var}}$=3 |  3.051$\times 10^{-1}$ |  3.277 |  2.660$\times 10^{-1}$ |  3.760 |  8.448$\times 10^{-2}$ |  11.838
$N_{\mathrm{var}}$=4 |  2.698$\times 10^{-1}$ |  3.706 |  8.512$\times 10^{-2}$ |  11.748
$N_{\mathrm{var}}$=5 |  2.688$\times 10^{-1}$ |  3.721 |  8.422$\times 10^{-2}$ |  11.874
$N_{\mathrm{var}}$=6 |  2.660$\times 10^{-1}$ |  3.760 |  8.420$\times 10^{-2}$ |  11.876
$N_{\mathrm{sec}}=8$ |  $N_{\mathrm{var}}$=3 |  4.503$\times 10^{-1}$ |  2.221 |  4.448$\times 10^{-1}$ |  2.248 |  7.207$\times 10^{-2}$ |  13.875
$N_{\mathrm{var}}$=4 |  4.444$\times 10^{-1}$ |  2.250 |  7.357$\times 10^{-2}$ |  13.593
$N_{\mathrm{var}}$=5 |  4.437$\times 10^{-1}$ |  2.254 |  7.337$\times 10^{-2}$ |  13.629
$N_{\mathrm{var}}$=6 |  4.443$\times 10^{-1}$ |  2.251 |  7.307$\times 10^{-2}$ |  13.685
$N_{\mathrm{sec}}=10$ |  $N_{\mathrm{var}}$=3 |  2.592$\times 10^{-1}$ |  3.858 |  2.414$\times 10^{-1}$ |  4.143 |  7.887$\times 10^{-2}$ |  12.679
$N_{\mathrm{var}}$=4 |  2.458$\times 10^{-1}$ |  4.068 |  7.946$\times 10^{-2}$ |  12.585
$N_{\mathrm{var}}$=5 |  2.451$\times 10^{-1}$ |  4.080 |  7.871$\times 10^{-2}$ |  12.704
$N_{\mathrm{var}}$=6 |  2.414$\times 10^{-1}$ |  4.143 |  7.858$\times 10^{-2}$ |  12.726
$N_{\mathrm{sec}}=12$ |  $N_{\mathrm{var}}$=3 |  2.348$\times 10^{-1}$ |  4.260 |  1.636$\times 10^{-1}$ |  6.113 |  7.393$\times 10^{-2}$ |  13.526
$N_{\mathrm{var}}$=4 |  1.642$\times 10^{-1}$ |  6.091 |  7.397$\times 10^{-2}$ |  13.519
$N_{\mathrm{var}}$=5 |  1.639$\times 10^{-1}$ |  6.102 |  7.395$\times 10^{-2}$ |  13.522
$N_{\mathrm{var}}$=6 |  1.634$\times 10^{-1}$ |  6.122 |  7.393$\times 10^{-2}$ |  13.526
$N_{\mathrm{sec}}=14$ |  $N_{\mathrm{var}}$=3 |  2.401$\times 10^{-1}$ |  4.165 |  2.331$\times 10^{-1}$ |  4.290 |  7.698$\times 10^{-2}$ |  12.990
$N_{\mathrm{var}}$=4 |  2.371$\times 10^{-1}$ |  4.217 |  7.751$\times 10^{-2}$ |  12.902
$N_{\mathrm{var}}$=5 |  2.364$\times 10^{-1}$ |  4.229 |  7.688$\times 10^{-2}$ |  13.008
$N_{\mathrm{var}}$=6 |  2.331$\times 10^{-1}$ |  4.291 |  7.669$\times 10^{-2}$ |  13.039
$N_{\mathrm{sec}}=16$ |  $N_{\mathrm{var}}$=3 |  2.994$\times 10^{-1}$ |  3.340 |  2.979$\times 10^{-1}$ |  3.357 |  7.265$\times 10^{-2}$ |  13.764
$N_{\mathrm{var}}$=4 |  2.950$\times 10^{-1}$ |  3.390 |  7.437$\times 10^{-2}$ |  13.446
$N_{\mathrm{var}}$=5 |  2.957$\times 10^{-1}$ |  3.382 |  7.392$\times 10^{-2}$ |  13.528
$N_{\mathrm{var}}$=6 |  2.960$\times 10^{-1}$ |  3.379 |  7.372$\times 10^{-2}$ |  13.565
$N_{\mathrm{sec}}=18$ |  $N_{\mathrm{var}}$=3 |  2.465$\times 10^{-1}$ |  4.057 |  2.279$\times 10^{-1}$ |  4.388 |  7.621$\times 10^{-2}$ |  13.121
$N_{\mathrm{var}}$=4 |  2.316$\times 10^{-1}$ |  4.318 |  7.663$\times 10^{-2}$ |  13.049
$N_{\mathrm{var}}$=5 |  2.313$\times 10^{-1}$ |  4.323 |  7.603$\times 10^{-2}$ |  13.152
$N_{\mathrm{var}}$=6 |  2.279$\times 10^{-1}$ |  4.389 |  7.593$\times 10^{-2}$ |  13.171
$N_{\mathrm{sec}}=20$ |  $N_{\mathrm{var}}$=3 |  2.285$\times 10^{-1}$ |  4.376 |  1.733$\times 10^{-1}$ |  5.772 |  7.371$\times 10^{-2}$ |  13.566
$N_{\mathrm{var}}$=4 |  1.736$\times 10^{-1}$ |  5.760 |  7.375$\times 10^{-2}$ |  13.559
$N_{\mathrm{var}}$=5 |  1.731$\times 10^{-1}$ |  5.777 |  7.373$\times 10^{-2}$ |  13.563
$N_{\mathrm{var}}$=6 |  1.729$\times 10^{-1}$ |  5.785 |  7.371$\times 10^{-2}$ |  13.566

(a) Type A, $N_{\mathrm{sec}}=4$, shape optimization

(b) Type A, $N_{\mathrm{sec}}=4$, shape + $\kappa$ optimization

(c) Type A, $N_{\mathrm{sec}}=12$, shape optimization

(d) Type A, $N_{\mathrm{sec}}=12$, shape + $\kappa$ optimization

Figure 14: Detailed view of flux concentration in the optimized shape at the
junction for Type-A configuration for two cases (II) shape optimization and
(III) shape optimization combined with conductivity optimization
$N_{\mathrm{sec}}$=4, 12.

Next, we conducted an observation by allowing the thermal conductivity to vary
with shape optimization. As our optimization uses non-gradient based
algorithm, it is easy to implement the conductivity as a design variable apart
from regular design variables. In this case, the materials in the concentrator
and core region are allowed to change their conductivity patch-wise. The upper
and lower limits for the conductivity box constraints are provided as
$\kappa_{\rm copper}=398$ W/mK and $\kappa_{\rm PDMS}=0.27$ W/mK,
respectively.

The results are presented in Figure 12 and Table 2. Here three different cases
are compared: (I) without any optimization, (II) only shape optimization, and
(III) shape optimization combined with conductivity optimization. It is
observed that for each $N_{\mathrm{sec}}$ and $N_{\mathrm{var}}$ value, case-
III provides the best resultss. The concentrated flux at
$\Gamma_{\mathrm{flux}}$ is approximately 10 to 13 times as large as for a
base material plate and 2 to 3 times as large as for case-II & 3 to 6 times as
large as for case-I.

In order to understand the material distribution and corresponding optimized
shape, the streamline plots for case-II and case-III are compared in Figure
13. In the figure, the color of a patch represents its thermal conductivity
obtained by optimization. The color scale is defined over a logarithmic scale,
where the orange color denotes a material with higher conductivity and green
denotes a material with lower conductivity. However, the base material is
separately kept in grey color, independent of the color scale.

It is evident from Figure 13 that, for all configurations, the core region
reaches the upper limit of conductivity range. Subsequently, the imposition of
continuity conditions at interfaces as well as the symmetry constraints allow
the copper sectors to cover a large periphery along outer radius. This large
opening works as a gate to attract more flux, before eventually driving the
same flux to the core region. Therefore, the large value of
$\Psi_{\mathrm{flux}}$ for case III compared to case II is justified. Thus, in
order to maximize the performance of the concentrator, the material with the
highest conductivity should be associated to the inner core.

Apart from that, if we focus on the sector-wise material distribution in the
concentrator region, it remains the same as before for Type-A configurations
(Figures 13(a), 13(e), 13(i), 13(m)). On the other hand, for Type-C
configuration (Figures 13(c), 13(g), 13(k), 13(o)), the materials in the
concentrator region are exchanged. In other words, the sectors inline with
$\Gamma_{\mathrm{flux}}$ should be of the lowest conductivity. Therefore, for
Type-B (Figures 13(b), 13(f), 13(j), 13(n)) and Type-D configurations (Figures
13(d), 13(h), 13(l), 13(p))(where it is not possible to arrange the lowest
conductivity along the sectors inline with $\Gamma_{\mathrm{flux}}$) for some
instances the sectors swap materials, for some instances they do not.

In case-II, one more point to note is that the flux concentration at the very
small interface created between the sector with lower conductivity and inner
core for Type-A configuration as shown in Figure 14. The fundamental reason
behind it is when the flux flows from a higher $\kappa$ patch to a lower
$\kappa$ patch, the flux diverges before crossing the interface due to the
normal flux continuity condition. The same effect occurs when the flux flows
from the copper sector to the inner core which has comparatively low $\kappa$
nickel steel material). Eventually the diverged flux gets accumulated at the
junction of three materials (shown in the detailed view in Figure 14) before
entering in the inner core. As this flux concentration can cause numerical
burden, we have provided a limit on the minimum size of the interface by box
constraints. However in case-III, the inner core takes the highest possible
$\kappa$-value and avoids the flux concentration. The smooth flow of flux from
the copper sector to the inner core is evident in Figure 14.

####  4.4.5 Application of inclined flux on optimized shape

(a) $\phi=-3\pi/8$

(b) $\phi=-\pi/4$

(c) $\phi=-\pi/8$

(d) $\phi=0$

(e) $\phi=\pi/8$

(f) $\phi=\pi/4$

(g) $\phi=3\pi/8$

Figure 15: Variation of concentrator function $\Psi_{\mathrm{flux}}$ of the
optimized shape with respect to number of sectors for constant flux applied at
angle $\phi$ = $-3\pi/8$, $-\pi/4$, $-\pi/8$, $0$, $\pi/8$, $\pi/4$ and
$3\pi/8$ for three case (I) without any optimization (II) only shape
optimization, and (III) shape optimization combined with conductivity
optimization. $N_{\mathrm{var}}=$4\. The optimized geometries for case-III
perform better than case-I and II.

In this study, we test the optimized geometries from case-I, II and III under
application of constant flux at different directions. The constant flux of
magnitude $Q_{N}=7.143\times 10^{3}$ W is applied as the Neumann boundary
conditions on all four sides of the plate. Also, the upper left corner is
fixed at 300 K to make the boundary value problem well-posed. We consider
seven different values for the angle $\phi$ of the applied flux with the
$x$-axis ($\phi$=$-3\pi/8$, $-\pi/4$, $-\pi/8$, $0$, $\pi/8$, $\pi/4$ and
$3\pi/8$). Since, as shown earlier, there is no apparent variation of
$\Psi_{\mathrm{flux}}$ with respect to $N_{\mathrm{var}}$, we only consider
one value $N_{\mathrm{var}}$, $N_{\mathrm{var}}$=4. The results are presented
in Figure 15.

From Figure 15, it is evident that the optimized geometries achieved for case-
III, which is combined shape and conductivity optimization, perform better
than case-I and II. On the other hand, we do not observe any obvious
superiority between case-I and II optimized geometries. Since the geometries
for case-II are optimized for the specific boundary conditions, they will not
necessarily perform well for other cases. However, the results for case-III
are understandable considering relatively wider opening for copper sectors to
attract more flux towards the core. One more thing to note for case-III
geometries is that the effect of configuration (A/B/C/D) of sectors on
$\Psi_{\mathrm{flux}}$ diminishes for larger $N_{\mathrm{sec}}$, and all
configurations perform equally well. This comes from the fact that for larger
$N_{\mathrm{sec}}$, almost the whole outer and inner peripheries are covered
by copper and a very small part remains for PDMS. Nevertheless their thin
structure, PDMS sectors are needed inline with $\Gamma_{\mathrm{flux}}$ to
provide a barrier to the flux in $\Omega_{\mathrm{design}}$, and guide it
towards $\Omega_{\mathrm{in}}$ and eventually through
$\Gamma_{\mathrm{flux}}$.

##  5 Thermal cloak-concentrator

###  5.1 Problem definition

In this example, we optimize the thermal cloak-concentrator for simultaneously
concentrating flux and cloaking the inner core. The example is referred from
[29] and the problem details are kept unchanged to compare the results in the
later stage. However, we apply different tools from the tools applied in [29].
Instead of topology optimization with covariance matrix adaptation evolution
strategy (CMA-ES), we utilize shape optimization with PSO algorithm.
Additionally, we exploit isogeometric analysis (IGA) instead of Lagrange
finite element method (FEM) to solve the boundary value problem. The
schematics of the problem is similar to the last example as shown in Figure 5,
however the dimensions are different. The dimensions of the current geometry
are shown in Figure 16(a). All dimensions are given relative to the outer
radius of the heat manipulator, $R_{D}$. Here, we take $R_{D}$=10. We consider
iron ($\kappa_{\rm iron}=67$ W/mK) as the base material and copper & PDMS
(with slightly different values of conductivity from the last example) as the
materials for thermal cloak-concentrator. The conductivity of copper and PDMS
is taken as $\kappa_{\rm copper}=386$ W/mK and $\kappa_{\rm PDMS}=0.15$ W/mK
respectively. The left boundary works as a sink with $0^{0}$C constant
temperature, and the right boundary works as a source with $1^{0}$C.

(a) A metamaterial based cloak-concentrator in the base material plate

(b) An annulus ring-sector shaped heat manipulator in the base material plate

Figure 16: Schematic design of (a) A thermal cloak-concentrator embedded in
the base material plate to increase the flux concentration in
$\Omega_{\mathrm{in}}$, $\Omega_{\mathrm{design}}$ is the area of the cloak-
concentrator, where we want to optimize the shape, $\Omega_{\mathrm{out}}$ is
the outside area of remaining base material,
$\Omega=\Omega_{\mathrm{in}}\cup\Omega_{\mathrm{design}}\cup\Omega_{\mathrm{out}}$;
For our case, we assume $R_{D}$=10. The base material is iron. (b) An embedded
matermaterial-based thermal cloak-concentrator made of block type of structure
of copper and PDMS.

###  5.2 Objective function

The objective of a thermal cloak is to avoid the disturbance caused by an
obstacle/inner region, and mimic the temperature distribution in the outer
base material as if there was no obstacle/inner region. Here, we will try to
cloak the inner core, where we are also concentrating the flux, by a single
heat manipulator. The combined objective function is defined as,

$f_{\mathrm{obj}}=\underbrace{\left|\dfrac{1}{\Psi_{\mathrm{flux}}}\right|^{4}}_{f_{\mathrm{obj,flux}}}+\underbrace{\left|\Psi_{\mathrm{cloak}}\right|}_{f_{\mathrm{obj,cloak}}},$
(29)

where the first part is related to the concentrator’s objective and second is
related to the cloak’s objective. The concentrator function
$\Psi_{\mathrm{flux}}$ is same as defined in Equation (26), while the cloak
function $\Psi_{\mathrm{cloak}}$ is defined as below,

$\Psi_{\mathrm{cloak}}=\dfrac{1}{\widetilde{\Psi}_{\mathrm{cloak}}}\int_{\Omega_{\mathrm{out}}}|T-\overline{T}|^{2}~{}d\Omega,$
(30)

with $\widetilde{\Psi}_{\mathrm{cloak}}$ be the normalisation value given as,

$\widetilde{\Psi}_{\mathrm{cloak}}=\int_{\Omega_{\mathrm{out}}}|\widetilde{T}-\overline{T}|^{2}~{}d\Omega,$
(31)

where $T$, $\overline{T}$ are same as described in Section 4.2.
$\widetilde{T}$ is the temperature field when $\Omega_{\mathrm{design}}$ is
filled with PDMS material.

###  5.3 NURBS parameterization, design variables and constraints

Most of the aspects of NURBS approximations, design variables and constraints
are similar to the last example as mentioned Section 4.3. The geometry of the
cloak-concentrator is taken as a structure made of several concentric annular
rings connected at the periphery. Besides, each ring is made of sectors.
Overall, the geometry of a cloak-concentrator is a structure made of sector-
shaped blocks in annulus area as shown in Figure 16(b). Each block is
considered as a different patch, and continuity conditions at interfaces are
imposed by Nitsche’s method. The design variables, constraints and symmetry
conditions are kept unchanged from the last example. However, all the radius
between the inner and outer radius are considered as the extra design
variables in addition to the usual design variable.

###  5.4 Results and discussion

####  5.4.1 Conductivity optimization

At first, we need an initial topology for the cloak-concentrator to start the
optimization. Since the new features cannot be created in the geometry during
shape optimization, it is crucial to take an appropriate topology that can
ensure sufficient scope for shape optimization. To get the initial topology,
we run solely conductivity optimization. For the optimization, we impose
$x$-axis and $y$-axis symmetry. The results of the conductivity optimization
will exhibit suitable conductivity values for all blocks, which will assist us
to define initial material distribution in cloak-concentrator. We consider two
initial geometries: one with $N_{\mathrm{ring}}=3$ & $N_{\mathrm{sec}}=12$,
another with $N_{\mathrm{ring}}=5$ & $N_{\mathrm{sec}}=20$.

(a) Material distribution for $N_{\mathrm{ring}}$=3, $N_{\mathrm{sec}}$=12

(b) $|T-\overline{T}|$ for $N_{\mathrm{ring}}$=3, $N_{\mathrm{sec}}$=12,
$\Psi_{\mathrm{cloak}}=1.425\times 10^{-5}$

(c) Flux flow for $N_{\mathrm{ring}}$=3, $N_{\mathrm{sec}}$=12,
$\Psi_{\mathrm{flux}}=7.587$

(d) Material distribution for $N_{\mathrm{ring}}$=5, $N_{\mathrm{sec}}$=20

(e) $|T-\overline{T}|$ for $N_{\mathrm{ring}}$=5, $N_{\mathrm{sec}}$=20,
$\Psi_{\mathrm{cloak}}=7.904\times 10^{-6}$

(f) Flux flow for $N_{\mathrm{ring}}$=5, $N_{\mathrm{sec}}$=20,
$\Psi_{\mathrm{flux}}=8.500$

Figure 17: optimized block-wise $\kappa$-distribution of thermal cloak-
concentrator; temperature disturbance and flux concentration after adding
thermal cloak-concentrator in the base material plate. (a,b,c)
$N_{\mathrm{ring}}$=3, $N_{\mathrm{sec}}$=12 $f_{\mathrm{obj}}=3.161\times
10^{-4}$ (d,e,f) $N_{\mathrm{ring}}$=5, $N_{\mathrm{sec}}$=20,
$f_{\mathrm{obj}}=1.995\times 10^{-4}$. The optimized conductivity
distribution primarily creates fan-type of shapes of low-conductivity material
above and below inner core.

(a) $N_{\mathrm{ring}}$=3, $N_{\mathrm{sec}}$=12. Conductivity values for all
materials: $f_{\mathrm{obj}}=4.632\times 10^{-4}$ with
$\Psi_{\mathrm{flux}}=7.688$ and $\Psi_{\mathrm{claok}}=1.769\times 10^{-4}$

(b) $N_{\mathrm{ring}}$=5, $N_{\mathrm{sec}}$=20. Conductivity values for all
materials: $f_{\mathrm{obj}}=5.367\times 10^{-4}$ with
$\Psi_{\mathrm{flux}}=8.212$ and $\Psi_{\mathrm{claok}}=3.168\times 10^{-4}$

Figure 18: Block-wise initial material distribution for shape optimization. To
avoid large number of materials involved into the geometry, the blocks with
close conductivity values (from conductivity optimization) are taken as a
single material with a round-off conductivity value.

The results obtained from conductivity optimization are shown in Figure 17.
Figure 17 shows the block-wise conductivity distribution, temperature
disturbance $T-\overline{T}$ and flux flow. For the geometry with
$N_{\mathrm{ring}}=3$ & $N_{\mathrm{sec}}=12$, the optimization produces
$f_{\mathrm{obj}}=3.161\times 10^{-4}$ with $\Psi_{\mathrm{flux}}=7.587$ and
$\Psi_{\mathrm{cloak}}=1.425\times 10^{-5}$. Similarly for the geometry with
$N_{\mathrm{ring}}=5$ & $N_{\mathrm{sec}}=20$, $f_{\mathrm{obj}}=1.995\times
10^{-4}$ with $\Psi_{\mathrm{flux}}=8.500$ and
$\Psi_{\mathrm{cloak}}=7.904\times 10^{-6}$. The optimized conductivity
distribution primarily creates fan-type of shapes of low-conductivity material
above and below inner core, similar to the optimized topology mentioned in
[29]. In addition to that, similar to the dam structure (a structure made of
low $\kappa$ material obstructing the flux as it enters in cloak-concentrator
along the $x$-axis) observed in [29], we also get low-conductivity material
distribution in the direction of the incoming flux for $N_{\mathrm{ring}}=5$ &
$N_{\mathrm{sec}}=20$ geometry. As mentioned in [29], narrow structures of
high $\kappa$ material appear between fan shape and base material in the
optimal geometries. These narrow structures diminish the temperature
disturbance produced by fan shapes. In our material distribution, the same
cloaking effect is provided by the blocks with conductivity value close to the
base material. After analysing the block-wise conductivity distribution, it is
observed that some of the blocks have very close conductivity values. Now, to
avoid increasing the number of materials involved in the geometry, the blocks
with close conductivity values are taken as a single material with a round-off
conductivity value. The updated material distribution and corresponding
conductivity values are shown in Figure 18. However, the new conductivity
values, which are round-off values of optimized values, changes
$f_{\mathrm{obj}}$, $\Psi_{\mathrm{flux}}$ and $\Psi_{\mathrm{cloak}}$ to some
extent. The altered values are shown in the caption. Lastly, the updated
material distribution of block-type of geometry will be used as the initial
topology for the shape optimization problem in the next stage.

(a) Optimized shape for $N_{\mathrm{var}}$=23

(b) $|T-\overline{T}|$ for $N_{\mathrm{var}}$=23,
$\Psi_{\mathrm{cloak}}=4.749\times 10^{-6}$

(c) Flux flow for $N_{\mathrm{var}}$=23, $\Psi_{\mathrm{flux}}=8.077$.

(d) Optimized shape for $N_{\mathrm{var}}$=50

(e) $|T-\overline{T}|$ for $N_{\mathrm{var}}$=50,
$\Psi_{\mathrm{cloak}}=7.451\times 10^{-6}$

(f) Flux flow for $N_{\mathrm{var}}$=50, $\Psi_{\mathrm{flux}}=8.012$

Figure 19: optimized shape of thermal cloak-concentrator; temperature
disturbance ($T-\overline{T}$) and flux flow after adding thermal cloak-
concentrator in the base material plate for $N_{\mathrm{ring}}$=3,
$N_{\mathrm{sec}}$=12 (a,b,c) $N_{\mathrm{var}}=23$,
$f_{\mathrm{obj}}=2.397\times 10^{-4}$ (d,e,f) $N_{\mathrm{var}}=50$,
$f_{\mathrm{obj}}=2.501\times 10^{-4}$. For both cases, the concentrate flux
is around 8 times as large as for the flat plate. The normalized temperature
difference in the outer region is in the order of $10^{-6}$.

(a) Opt. shape $N_{\mathrm{ring}}$=5, $N_{\mathrm{sec}}$=20,
$N_{\mathrm{var}}$=59

(b) $|T-\overline{T}|$ for $N_{\mathrm{ring}}$=5, $N_{\mathrm{sec}}$=20,
$N_{\mathrm{var}}$=59, $\Psi_{\mathrm{cloak}}=7.363\times 10^{-6}$

(c) Flux flow for $N_{\mathrm{ring}}$=5, $N_{\mathrm{sec}}$=20,
$N_{\mathrm{var}}$=59, $\Psi_{\mathrm{flux}}=8.316$

(d) Opt. shape $N_{\mathrm{ring}}$=5, $N_{\mathrm{sec}}$=20,
$N_{\mathrm{var}}$=132

(e) $|T-\overline{T}|$ for $N_{\mathrm{ring}}$=5, $N_{\mathrm{sec}}$=20,
$N_{\mathrm{var}}$=132, $\Psi_{\mathrm{cloak}}=1.232\times 10^{-5}$

(f) Flux flow for $N_{\mathrm{ring}}$=5, $N_{\mathrm{sec}}$=20,
$N_{\mathrm{var}}$=132, $\Psi_{\mathrm{flux}}=8.318$

Figure 20: optimized shape of thermal cloak-concentrator; temperature disturbance ($T-\overline{T}$) and flux flow after adding thermal cloak-concentrator in the base material plate for $N_{\mathrm{ring}}$=5, $N_{\mathrm{sec}}$=20 (a,b,c) $N_{\mathrm{var}}=59$, $f_{\mathrm{obj}}=2.165\times 10^{-4}$ (d,e,f) $N_{\mathrm{var}}=132$, $f_{\mathrm{obj}}=2.212\times 10^{-4}$. For both cases, the concentrate flux is around 8 times as large as for the flat plate. The normalized temperature difference in the outer region is in the order of $10^{-5}$-$10^{-6}$. Table 3: Comparison between literature, conductivity optimization and shape optimization results for thermal cloak concentrator. Optimization method |  Number of design variables |  $N_{\mathrm{ring}}=3$, $N_{\mathrm{sec}}=12$ |  $N_{\mathrm{ring}}=5$, $N_{\mathrm{sec}}=20$
---|---|---|---
|  |  $\Psi_{\mathrm{cloak}}$ |  $\Psi_{\mathrm{flux}}$ |  $\Psi_{\mathrm{cloak}}$ |  $\Psi_{\mathrm{flux}}$
Literature results [29] |  2850 |  7.64$\times 10^{-5}$ |  7.16 |  7.64$\times 10^{-5}$ |  7.16
Conductivity optimization |  |  1.425$\times 10^{-5}$ |  7.587 |  7.904$\times 10^{-6}$ |  8.500
Initial geometry for shape optimization |  |  1.769$\times 10^{-4}$ |  7.688 |  3.168$\times 10^{-4}$ |  8.212
Shape optimization |  23 |  4.749$\times 10^{-6}$ |  8.077 |  |
|  50 |  7.451$\times 10^{-6}$ |  8.012 |  |
|  59 |  |  |  7.363$\times 10^{-6}$ |  8.316
|  132 |  |  |  1.232$\times 10^{-5}$ |  8.318

####  5.4.2 Shape optimization

We now perform shape optimization for the geometries shown in Figure 18. The
results of the shape optimization are presented in Figures 19-20. For each
geometry we consider two different values of $N_{\mathrm{var}}$. The optimized
shapes for all cases generate negligible disturbance in the temperature
profile ($T-\overline{T}$) as shown in Figures 19-20. Also, we can observe
that the material distribution in the inner rings guide the flux towards the
inner core, and the material distribution in the outermost ring suppresses all
the temperature disturbance created by the inner region. Thus, we can say that
the shape and conductivity values of the patches on the outermost ring are
very critical to the cloaking performance. In Table 3, results from
literature, conductivity optimization, and shape optimization are compared.
From the comparison, we can say that, the shape optimization process can
produce slightly better results than the topology optimization results from
literature [29] with much fewer design variables. However, this advantage also
comes with the complexity of using more than two materials. We also use high
conductivity copper in the inner core, instead of base material iron as in
[29]. However, the increment $N_{\mathrm{var}}$ does not provide any specific
advantages.

For the first type of geometry with $N_{\mathrm{ring}}=3$ &
$N_{\mathrm{sec}}=12$, $N_{\mathrm{var}}=23$ and $N_{\mathrm{var}}=50$, while
for the second type of geometry with $N_{\mathrm{ring}}=5$ &
$N_{\mathrm{sec}}=20$, $N_{\mathrm{var}}=59$ and $N_{\mathrm{var}}=132$. The
optimized shape effectively reduces the disturbance caused by inner core as
well as concentrate the flux. For all four cases, the concentrate flux is
around 7 times as large as for the flat plate. The normalized temperature
difference in the outer region is very small (mostly in the order of
$10^{-6}$). The exact values of $\Psi_{\mathrm{flux}}$ and
$\Psi_{\mathrm{cloak}}$ are presented in Figures 19-20. The optimized shape
results are $10^{-1}$ order better (both in $\Psi_{\mathrm{flux}}$ and
$\Psi_{\mathrm{cloak}}$) than the results from literature.

##  6 Conclusions

In the present article, we investigated the shape optimization method for
metam-
aterial-based heat manipulators. The approach works as a tool to find better
geometries with good manufacturability and higher efficiency. The proposed
method utilized a gradient-free Particle Swarm Optimization (PSO) algorithm.
The geometry and solution fields are approximated using NURBS basis functions.
It enables easy control of shape and provides good geometrical accuracy
without meshing, remeshing, or any special post-processing approaches. Also,
Nitsche’s method is used to impose interface continuity conditions. We
presented examples of a thermal concentrator and a thermal cloak-concentrator
to demonstrate the efficiency of the proposed method.

For the thermal concentrator, we analyzed sector-type geometries made of two
materials. We studied the effect of geometric parameters such as the number of
sectors, the number of design variables, the imposed symmetry conditions, the
conductivity and spatial distribution of member materials, and the direction
of incident flux on the optimized shape and performance of the thermal
concentrator. The shape optimization indicates that,

  * 1.

With applied symmetry across $x$-axis and center axes of sectors, the
optimized geometries can accumulate heat flux 2 to 7 times as large as for un-
optimized geometries with straight sector edges.

  * 2.

Type A ($N_{\mathrm{sec}}=8\ell-4$, $\ell=1,2,...$) configurations perform
better than other types of configurations, as the low $\kappa$ sectors for
Type A configurations lie inline with $\Gamma_{\mathrm{flux}}$ (across which
the concentrated flux is measured) in $\Omega_{\mathrm{design}}$. The low
$\kappa$ material obstructs the flux flow circumventing the core and guides it
towards the core, and eventually, through $\Gamma_{\mathrm{flux}}$.

  * 3.

By associating the high $\kappa$ material to the core, the openings (along the
outer radius) for the high $\kappa$ sectors increase significantly. That
provides a wider area to draw more flux towards the core. Also, the heat flux
transfer from $\Omega_{\mathrm{design}}$ to $\Omega_{\mathrm{in}}$ becomes
smoother.

  * 4.

The flux concentration capacity of the concentrator can be improved up to 2 to
3 times by placing high $\kappa$ material in the core and low $\kappa$
material inline with $\Gamma_{\mathrm{flux}}$ in $\Omega_{\mathrm{design}}$.

For the thermal cloak-concentrator, we analyzed the sector-shaped block-type
geometries. We studied the effect of geometric parameters such as the number
of sectors and the number of rings, the conductivity and spatial distribution
of member materials. The shape optimization indicates that,

  * 1.

With the optimized geometries, the thermal cloak-concentrator can collect 7 to
8 times as much flux as the un-optimized geometries. It also cloaks the core
by diminishing the normalized temperature disturbance up to a negligible value
of order $10^{-5}$ -$10^{-6}$.

  * 2.

With more member materials, the results in terms of the objective function are
around $10^{-1}$ order better than the literature results (which are based on
topology optimization) with fewer design variables.

The proposed method is applicable to any heat manipulator including thermal
cloaks, active cloaks, heat inverters, heat illusions, and so on. Since it
uses a gradient-free optimization algorithm, it can also handle non-
differentiable problems. However, the limitation of the method is that it can
not generate a new topology. Given that, the initial topology is crucial for
the optimization. Furthermore, it can be extended to 3D applications, which is
the subject of future studies.

## Acknowledgements

Stéphane P.A. Bordas acknowledges the funding from the European Union’s
Horizon 2020 research and innovation programme under grant agreement No 811099
TWINNING Project DRIVEN for the University of Luxembourg.

## References

  * Chen et al. [2008] T. Chen, C.-N. Weng, J.-S. Chen,  Cloak for curvilinearly anisotropic media in conduction,  Applied Physics Letters 93 (2008) 114103.
  * Narayana and Sato [2012] S. Narayana, Y. Sato,  Heat flux manipulation with engineered thermal materials,  Phys. Rev. Lett. 108 (2012) 214303.
  * Guenneau et al. [2012] S. Guenneau, C. Amra, D. Veynante,  Transformation thermodynamics: cloaking and concentrating heat flux,  Opt. Express 20 (2012) 8207–8218.
  * Schittny et al. [2013] R. Schittny, M. Kadic, S. Guenneau, M. Wegener,  Experiments on transformation thermodynamics: Molding the flow of heat,  Phys. Rev. Lett. 110 (2013) 195901.
  * Han et al. [2014a] T. Han, X. Bai, D. Gao, J. T. L. Thong, B. Li, C.-W. Qiu,  Experimental demonstration of a bilayer thermal cloak,  Phys. Rev. Lett. 112 (2014a) 054302.
  * Han et al. [2014b] T. Han, X. Bai, J. T. L. Thong, B. Li, C.-W. Qiu,  Full control and manipulation of heat signatures: Cloaking, camouflage and thermal metamaterials,  Advanced Materials 26 (2014b) 1731–1734.
  * Sklan et al. [2016] S. R. Sklan, X. Bai, B. Li, X. Zhang,  Detecting Thermal Cloaks via Transient Effects,  Scientific Reports 6 (2016) 32915.
  * Li et al. [2019] J. Li, Y. Li, T. Li, W. Wang, L. Li, C.-W. Qiu,  Doublet thermal metadevice,  Phys. Rev. Applied 11 (2019) 044021.
  * Fujii and Akimoto [2019] G. Fujii, Y. Akimoto,  Topology-optimized thermal carpet cloak expressed by an immersed-boundary level-set method via a covariance matrix adaptation evolution strategy,  International Journal of Heat and Mass Transfer 137 (2019) 1312–1322.
  * Shen et al. [2016] X. Shen, C. Jiang, Y. Li, J. Huang,  Thermal metamaterial for convergent transfer of conductive heat with high efficiency,  Applied Physics Letters 109 (2016) 201906.
  * Peng et al. [2020] Y.-G. Peng, Y. Li, P.-C. Cao, X.-F. Zhu, C.-W. Qiu,  3d printed meta-helmet for wide-angle thermal camouflages,  Advanced Functional Materials 30 (2020) 2002061.
  * Kennedy and Eberhart [1995] J. Kennedy, R. Eberhart,  Particle swarm optimization,  in: Proceedings of ICNN’95 - International Conference on Neural Networks, volume 4, pp. 1942–1948 vol.4.
  * Hughes et al. [2005] T. Hughes, J. Cottrell, Y. Bazilevs,  Isogeometric analysis: Cad, finite elements, nurbs, exact geometry and mesh refinement,  Computer Methods in Applied Mechanics and Engineering 194 (2005) 4135–4195.
  * Hassani et al. [2012] B. Hassani, M. Khanzadi, S. M. Tavakkoli,  An isogeometrical approach to structural topology optimization by optimality criteria,  Structural and Multidisciplinary Optimization 45 (2012) 223–233.
  * Shojaee et al. [2012] S. Shojaee, M. Mohamadianb, N. a. Valizadeh,  Composition of isogeometric analysis with level set method for structural topology optimization,  International Journal of Optimization in Civil Engineering 2 (2012).
  * Dedè et al. [2012] L. Dedè, M. J. Borden, T. J. Hughes,  Isogeometric Analysis for Topology Optimization with a Phase Field Model,  Archives of Computational Methods in Engineering 19 (2012) 427–465.
  * Tavakkoli et al. [2013] S. M. Tavakkoli, B. Hassani, H. a. Ghasemnejad,  Isogeometric topology optimization of structures by using mma,  International Journal of Optimization in Civil Engineering 3 (2013).
  * Tavakkoli and Hassani [2014] S. M. Tavakkoli, B. a. Hassani,  Isogeometric topology optimization by using optimality criteria and implicit function,  International Journal of Optimization in Civil Engineering 4 (2014).
  * Seo et al. [2010] Y.-D. Seo, H.-J. Kim, S.-K. Youn,  Isogeometric topology optimization using trimmed spline surfaces,  Computer Methods in Applied Mechanics and Engineering 199 (2010) 3270–3296.
  * Qian [2013] X. Qian,  Topology optimization in b-spline space,  Computer Methods in Applied Mechanics and Engineering 265 (2013) 15–35.
  * Espath et al. [2011] L. F. R. Espath, R. V. Linn, A. M. Awruch,  Shape optimization of shell structures based on nurbs description using automatic differentiation,  International Journal for Numerical Methods in Engineering 88 (2011) 613–636.
  * Koo et al. [2013] B. Koo, M. Yoon, S. Cho,  Isogeometric shape design sensitivity analysis using transformed basis functions for kronecker delta property,  Computer Methods in Applied Mechanics and Engineering 253 (2013) 505–516.
  * Park et al. [2013] B. U. Park, Y. D. Seo, O. Sigmund, S. K. Youn,  Shape optimization of the stokes flow problem based on isogeometric analysis,  Structural and Multidisciplinary Optimization 48 (2013) 965–977.
  * Blanchard et al. [2014] L. Blanchard, R. Duvigneau, A. V. Vuong, B. Simeon,  Shape Gradient for Isogeometric Structural Design,  Journal of Optimization Theory and Applications 161 (2014) 361–367.
  * Kiendl et al. [2014] J. Kiendl, R. Schmidt, R. Wüchner, K.-U. Bletzinger,  Isogeometric shape optimization of shells using semi-analytical sensitivity analysis and sensitivity weighting,  Computer Methods in Applied Mechanics and Engineering 274 (2014) 148–167.
  * Lian et al. [2016] H. Lian, P. Kerfriden, S. P. A. Bordas,  Implementation of regularized isogeometric boundary element methods for gradient-based shape optimization in two-dimensional linear elasticity,  International Journal for Numerical Methods in Engineering 106 (2016) 972–1017.
  * Lian et al. [2017] H. Lian, P. Kerfriden, S. Bordas,  Shape optimization directly from cad: An isogeometric boundary element approach using t-splines,  Computer Methods in Applied Mechanics and Engineering 317 (2017) 1–41.
  * Fujii and Akimoto [2019] G. Fujii, Y. Akimoto,  Optimizing the structural topology of bifunctional invisible cloak manipulating heat flux and direct current,  Applied Physics Letters 115 (2019) 174101.
  * Fujii and Akimoto [2020] G. Fujii, Y. Akimoto,  Cloaking a concentrator in thermal conduction via topology optimization,  International Journal of Heat and Mass Transfer (2020).
  * Nguyen et al. [2014] V. P. Nguyen, P. Kerfriden, M. Brino, S. P. Bordas, E. Bonisoli,  Nitsche’s method for two and three dimensional NURBS patch coupling,  Computational Mechanics (2014).
  * Hu et al. [2018] Q. Hu, F. Chouly, P. Hu, G. Cheng, S. P. Bordas,  Skew-symmetric nitsche’s formulation in isogeometric analysis: Dirichlet and symmetry conditions, patch coupling and frictionless contact,  Computer Methods in Applied Mechanics and Engineering 341 (2018) 188–220.
  * Chen and Yuan Lei [2015] F. Chen, D. Yuan Lei,  Experimental Realization of Extreme Heat Flux Concentration with Easy-to-Make Thermal Metamaterials,  Scientific Reports (2015).