A reverse thermal cloak design method based on inverse problem theory

Jun Guo, Zhiguo Qu^*, Xueliang Wang

MOE Key Laboratory of Thermo-Fluid Science and Engineering, School of Energy and Power Engineering, Xi'an Jiaotong University, 710049, China

*Corresponding author:   Fax: +86 029 82668543    Email:  zgqu@mail.xjtu.edu.cn

 

Abstract

A reverse thermal cloak design method is developed based on inverse problem theory. The design of thermal conductivity distribution within a cloak can be regarded as the density-based topology optimization with that penalization parameter equals to 1. The relative density is explored using conjugate gradient method with an adjoint equation (CGMAE). A prescribed objective function is used to evaluate the design result, and the function is no larger than 6.41×10^-6. The maximum deviation between background temperature with and without cloak position is no larger than 0.9% for background with uniform gradient and no larger than 5% for background with non-uniform gradient. The dimension of a cloak can also be designed in case of known material thermal conductivity to avoid complex material synthesis process. Compared with coordinate transformation, the reverse design method is able to remove anisotropy of thermal conductivity and is appropriate for design of irregularly shaped cloak. It is easy to implement programming for proposed method and can boost the flexibility of cloak design. Thus, the proposed method offers a practicable scheme for thermal cloak design to control heat conduction at will.

Table of Content

A reverse thermal cloak design method is developed based on inverse problem theory. Thermal conductivity and dimension of cloak, can be designed numerically.

 

 

 

 

 

Keywords: thermal cloak; reverse design method; inverse problem

---

1. Introduction

In general, cloak is a device that can hide an object from detection by an observer outside the cloak. On the one hand, the external physical field isn’t disturbed by the existence of the cloaked object. On the other hand, the cloaked object inside the cloak is free from the external physical field. The cloak is characterized by manipulating propagation direction of energy rather than reflecting or dissipating the energy, on which the conventional stealth technique primarily depends on. There are sufficient potential use like deception detection¹, energy conversion^2,3, enhancing sensitivity of sensors, noise reduction, etc. In essence, the camouflage effect of a certain cloak is realized by well-designed physical properties to control corresponding physical field. The electromagnetic cloak^4,5 is firstly introduced to carry out concealing objects from visible light and microwave. On these pioneering works’ heels, cloak is expanded to other fields like thermodynamics^6-8, elastodynamics ^9,10, acoustics ^11,12, and particle diffusion ^13,14.

 

Thermal cloak can thermally camouflage an object by means of regulating the heat flux that flows around an “invisible” region ^15,16. For one thing, a thermal detector outside the cloak cannot perceive the object from measuring the variation of background thermal field. For another, there is no heat flux penetrating into the cloaked region, which can shield the cloaked object¹⁷. In order to realize this effect, thermal conductivity distribution of the cloak material in steady state is carefully designed. Apart from material with exquisitely distributed thermal properties, heat sources ^18,19 can also play a role in realizing the thermal cloak. Thermal camouflage is broadened by thermal illusion in conduction^20-23 and convection²⁴, which can potentially transform an actual perception into a pre-controlled perception. Thermal meta-devices can also be extended from macroscale^25,26 to microscale^27,28.

 

At present,one kind of thermal cloak design concept is to calculate material parameters in the case when the cloak’s geometry and dimension are known. However, the satisfactory material usually doesn’t exist and should be synthesized in complex method. A perfect thermal cloak can be theoretically designed by coordinate transformation method ^{29, 30} on the strength of mathematical analysis. Aiming at extending application of thermal cloak, transformation thermodynamics on account of coordinate transformation has been fully developed in the condition of time dependence ³¹ and temperature dependence ³². Nevertheless, the well-designed cloak inevitably contains inhomogeneous and anisotropic properties, the trouble from which practicable application of the cloak suffers. It is a possible solution that the ideal cloak is simplified to a multilayered ³³ or a periodic structure ³⁴ composite. In addition to the simplification, another scheme to solve the problem is a bilayer cloak ^{35, 36} consisting of two layer materials whose thermal conductivities are constant. The design of bilayer cloak depends on no longer coordinate transformation but on scattering cancellation. Further, bilayer thermal cloak ³⁷ and carpet ³⁸ with irregular shapes were proposed. However, since the thermal conductivity of bilayer cloak is strictly determined by the cloak dimension, the design flexibility of bilayer cloak is limited by cloak size.

 

Another kind of thermal cloak design object is to obtain the cloak’s dimension for potential application in the case that the thermal conductivity of material is known. However, coordinate transformation cannot meet this demand, because thermal conductivity is designed by such method through transformation of the space dimension and the design process is irreversible. Given much mathematical analysis required on and of itself for coordinate transformation, the design process can be ever more difficult or even cannot be deployed as the cloak shape is more than regular. Even though the cloak shape can be expressed as analytical form ³⁹, the transformation is possibly complex, let alone the cloak whose shape cannot be denoted by mathematical expression. On that note, the universality of coordinate transformation may be limited by mathematical analysis.

 

In this study, aiming at above shortcomings of coordinate transformation, a reverse design method is proposed based on inverse problem theory^40,41.  Ref [40] and Ref. [41] provided referential experience about estimating thermal conductivity and parameters of thermal conductivity distribution form. Different from Ref. [40], in which thermal conductivity was the function of temperature and parameters of temperature items were calculated, thermal conductivity in this study is the function of radius and parameters of space items are estimated. And meanwhile, different from Ref. [41], in which thermal conductivity was time-dependent and the solution domain was one-dimensional, in this study thermal conductivity is time-independent and the solution domain is two-dimensional. Thermal conductivity distributions of cloaks with regular or irregular geometry are designed numerically by proposed method, and the anisotropy is removed. In addition, the thermal conductivity distribution of a cylindrical cloak is designed as a polynomial form. Such design result can improve the flexibility of engineering realization, and such cloak is approximately applicable to any arbitrary background thermal field with no less than 95% accuracy. Furthermore, the dimension of a cloak is designed in case of known thermal conductivity of material by proposed method. In this way, the cloak can be realized by existing material to avoid complex material synthesis process. Different from the method in Ref. [18], which was developed to calculate heat source intensity of an active cloak, this reverse design method is used to calculate thermal conductivity, parameters of thermal conductivity distribution form and cloak dimension.

 

2. The structure of designed thermal cloak model

Fig. 1(a) shows the schematic diagram for a hollow cylindrical thermal cloak, which is made up of two layer material. An object can be located in the cavity that is a cloaked region. As the cloak is positioned, a thermal detector outside the cloak cannot perceive the object, and the object can be thermally protected. To illustrate the cloak structure more clearly, a cross section is shown in Fig. 1(b). The space is divided into four zones denoted by Z1, Z2, Z3, and Z4 by two concentric annuluses. Z1 (r < R₁) is the cloaked region in which an object that should be thermally camouflaged and protected can be placed. Z2 and Z3 make up the thermal cloak. Z2 (R₁ < r < R₂) is composed of material with thermal conductivity k_h, which is constant and ideally infinite. Thus, there is hardly any temperature gradient in Z2. Temperature field in Z1 can be regarded as constant and the object can be thermally shielded. Z3 (R₂ < r < R₃) is used to recover the background (r > R₃, Z4) thermal field distorted by Z2. The cloak made up of Z2 and Z3 can lead to thermal shield in Z1 and no temperature distortion in Z4. Different from cloak model in Ref. [6], whose thermal conductivity in outer layer was constant, the thermal conductivity of model in Fig. 1 in outer layer is isotropic and inhomogeneous, or is designed as a polynomial form. It should be noted that Fig. 1 just shows a general cloak model. Cloaks of other shapes in this study also own the same two layers as that in Fig. 1.

 

 

     

Fig. 1 Schematic of the thermal cloak. (a) Cloak model. (b) Cross section structure of the cloak

 

3. Construction of background thermal field

Since thermal conductivity of Z2 in Fig. 1 is just ideally infinite, the design is not required. Thermal conductivity of Z3 should be calculated by the reverse design method based on inverse problem theory. The direct heat conduction problem is concerned with the determination of temperature field when the initial and boundary conditions, heat generation and physical parameters are specified. In contrast, inverse heat conduction problem (IHCP) involves the determination of the boundary conditions, heat generation, physical parameters (thermal conductivity in this study), etc., in the knowledge of the temperature field. On the strength of inverse problem theory, thermal conductivity of Z3 can be calculated according to original background thermal field. Thus, the background thermal fields should be constructed firstly.

 

Previously, to keep the design process general, the physical variables should be converted into dimensionless indexes. The computational domain is shown in Fig. 1(b). The dimensionless length X, Y are expressed as

.

(1)

(x, y) represents the Cartesian coordinate. L is the length of the computational domain. X and Y are scaled between 0 and 1. Such (X, Y) can represent all computational domain with different length. The dimensionless temperature T is denoted by

(2)

T^* is the real temperature of the whole domain.T_b^* and T_d^* represent the real temperatures of boundary Ω_p and Ω_d respectively. The dimensionless temperature T ranges from 0 to 1. On that note, the designed Z3 can be applied to background thermal field with any uniform temperature gradient. The form of dimensionless thermal conductivity K(X, Y) is

(3)

K^*(x, y) is the real thermal conductivity of the whole domain. K_b^* represents the thermal conductivity of the background. As a result, only if the designed thermal conductivity of Z3 is multiplied by thermal conductivity of certain background, it can be applied to different background.

 

Owing to that original background thermal field is the known condition, certain background thermal field should be constructed in advance through corresponding boundary conditions. A background thermal field with uniform temperature gradient is taken as an example. Boundary conditions of Ω_p and Ω_d are set as constant temperatures respectively. Other boundaries are kept adiabatic.

 

The original background thermal field with uniform temperature gradient can be constructed by the following governing equation and boundary conditions.

K(X, Y) of the whole domain is set as 1 in Eq. (4a). After Eqs. (4a) to (4d) are solved, the original background thermal field U(X, Y) with uniform temperature gradient can be obtained as shown in Fig. 2(a).

 

In addition, a background thermal field with non-uniform temperature gradient can also be constructed in the same way. The governing equation is Eq. (4a) and the boundary conditions are shown in Eqs. (5a)-5(d). The background thermal field with non-uniform temperature gradient is shown in Fig. 2(b).

 

4. Implementation of reverse design method without restriction for thermal conductivity distribution

Thermal cloak can manipulate heat flux by well-designed thermal conductivity distribution in steady state. Thus, thermal conductivity K(X, Y) is the design variable. In this section, thermal conductivity of a cloak is calculated by proposed method without restriction for thermal conductivity distribution, as long as the cloak function can be achieved. Two cloak design cases, in which the cross sections of cloaks are circular and heart-shaped, are taken as examples based on background thermal field with uniform temperature gradient as shown in Fig. 2(a).

 

The design of cloak can be regarded as the density-based topology optimization. Compared with existing design methodology based on covariance matrix adaptation evolution strategy (CMA-ES) ^42,43, the computational efficiency of density method is higher than that of CMA-ES. It is easier to implement programming for density method. In this study, relative density is calculated by conjugate gradient method with an adjoint equation (CGMAE), which divides the solving process into direct problem, sensitivity problem and adjoint problem. These three problems are all to solve Laplace equation. Thus, only one program can solve three problems. And meanwhile, the algorithm for solving Laplace equation has been fully optimized. So, program based on density method and CGMAE can be simple and accurate. But it is more difficult than CMA-ES to clearly express structural boundaries for density method. Particularly, thermal cloak designed by CMA-ES can be fabricated directly by bulk isotropic material because the design result is discrete. However, our design result is continuous. The designed cloak can be fabricated after processed by effective medium theory.

 

4.1 Design process for calculating thermal conductivity of cloak

The design of cloak can be regarded as the density-based topology optimization with that penalization parameter equals to 1. Thermal conductivity in Z3 varies continuously. The optimization of thermal conductivity is transformed into the optimization of relative density. The relative density can be iteratively calculated as an inverse problem in which a squared residue function is minimized using CGMAE. CGMAE is derived from the perturbation principle and transforms the inverse problem to the solution of three problems, namely, the direct problem, the sensitivity problem, and the adjoint problem. The direct problem is comprised of Eqs. (4a) to (4d) under circumstance of that thermal conductivity is set as infinitely great in Z2 and iteratively updated thermal conductivity in Z3. After solving the direct problem, the estimated temperature can be obtained.

An objective function is defined as

(6)

T_i is the estimated temperature at r = R₃, which is determined from the direct problem on the base of iteratively updated thermal conductivity in Z3. The subscript ‘i’ denotes the serial number of grid points at r = R₃. ‘m’ is the total amount of grid points at r = R_3. The evaluation function of Eq. (6) denotes the sum of difference values between estimated temperature and target temperature at all grid points at r = R₃. The number of m depends on the density of grid division. The denser the grid of the computational domain, the lager ‘m’ is and the larger the evaluation value is. So the evaluation function depends on ‘m’. The evaluation function is dimensionless, because the estimated temperature and target temperature are dimensionless. In order to realize the function of thermal cloak, the estimated background temperature should be identical with the original one. For reducing the amount of calculation, just original background temperature at r = R₃ should be set as the target. If function J is minimized, the background thermal field can recover.

Thermal conductivity of the material occupying Z3 is given by

(7)

where K_min and K_max are minimum and maximum values of thermal conductivity, ρ is the relative density, and p is the penalization parameter. Thermal conductivity in Z3 should be less than thermal conductivity of background, because thermal conductivity in Z2 is infinitely great. Thus K_min equals to 0 and K_max equals to 1. p is set as 1 in order to ensure that thermal conductivity varies continuously. Following constraint is imposed to relative density so as to avoid singularity of material thermal conductivity.

(8)

In iterative process, if ρ is less than 0.1, ρ is set as 0.1. And if ρ is greater than 1, ρ is set as 1. The optimization of thermal conductivity is transformed into the optimization of relative density. The following iterative process is considered for the optimization of ρ.

(9)

The subscript ‘n’ denotes the iteration index. β and P(X, Y) are step size and direction of descent (i.e. search direction) respectively.

So as to obtain P(X, Y), following adjoint problem should be solved.

	(10a)
	(10b)
	(10c)

δ is Dirac function. The gradient of function J is expressed as

(11)

The direction of descent P(X, Y) can be obtained by the following expression.

	(12a)
	(12b)

In order to obtain step size β, following sensitivity problem should be considered.

	(13a)
	(13b)
	(13c)
	(13d)

The step size β is obtained by

(14)

The boundary conditions of adjoint problem and sensitivity problem are constructed according to boundary conditions of background thermal fields in section 3. Above process can also be applied to design cloak based on background thermal field with non-uniform temperature gradient. The difference between cloak design based on background thermal field with uniform and non-uniform temperature gradient is reflected on boundary conditions of direct problem, adjoint problem and sensitivity problem. The process to determine thermal conductivity in Z3 using above steps is described as follows:

Step 1: Calculate the original background thermal field U(X, Y).

Step 2: Choose the initial relative density in Z3.

Step 3: Solve the direct problem.

Step 4: Solve the adjoint problem.

Step 5: Calculate the gradient of function J in Eq. (11).

Step 6: Calculate the direction of descent P(X, Y) in Eq. (12a) and Eq. (12b).

Step 7: Solve the sensitivity problem and calculate the step size β in Eq. (14).

Step 8: Update relative density and thermal conductivity in Z3 by Eq. (9) and Eq. (7).

Step 9: Calculate function J in Eq. (6). If J is small enough, stop the iteration and output thermal conductivity in Z3. If not, go back to step 3 and go on this cycle.

Step 1 to step 9 are used to calculate the thermal conductivity and the complex derivation process is not required to be considered. For coordinate transformation, the mathematical form of the cloak shape should be obtained firstly. Then the thermal conductivity distribution can be calculated. If the cloak shape is irregular, the coordinate transformation requires plenty of mathematical analysis and the process may be complex. The solving process of step 1 to step 9 can be programmed previously into the computer. The discrete form of a cloak shape, whose mathematical expression is not required, should be input into the program and the subsequent solving process can be delivered to the computer.

 

4.2 Design result of the cloak with regular geometry

The design of cloak, whose cross section is circular, is taken as an example as the regular geometry case. As shown in Fig. 1(b), R₁, R₂ and R₃ are set as 0.1, 0.13 and 0.2 respectively. The computational domain is evenly divided into Cartesian coordinate grids of 200×200. Step approximation method is used to mark each zone of the cloak. Fig. 3 is the design result and simulated temperature field showing the effect of the cloak. The iteration is 14 and the convergence time is 33.1 seconds. Fig. 2(a) is the original background thermal field. The designed thermal conductivity in Z3 is shown in Fig. 3(a). It is noted that there is no any anisotropy which appears everywhere in thermal cloaks designed by coordinate transformation method. The thermal conductivity is symmetric. According to Eq. (2), all dimensionless background temperatures are the same, no matter what the real background temperatures are. Thus, the thermal conductivity, which is calculated based on the dimensionless background temperature, could be applied to background with any uniform temperature gradient along Y-axis.

 

Due to the presence of the high thermal conductivity of Z2, contours are distorted but a “blank” region is created in Z1, as shown in Fig. 3(b). The “blank” means there is no temperature gradient and no heat flux infiltrating into Z1. Therefore, the object in Z1 can be thermally protected. In the enlarged view of Fig. 3(b), the isotherms of the background are completely distorted. An observer can be aware of the presence of the ‘blank’ region due to the temperature perturbation of the background. Z2 is like a black cloth covering an object. Although, an observer cannot feel the existence of the object, it can be aware of the black cloth. Thus, the background thermal field should be modified. Furthermore, the designed Z3 is positioned coaxially with Z2 and the tailored temperature field is shown in Fig. 3(c). Isotherms pass around Z1 and return undisturbed to their original trajectory as if the object in Z1 is absent. The background thermal field in Fig. 3(c) is recovered to the original one. In order to facilitate comparison, a black dashed ring with R₃ is added in Fig. 3(b) to indicate the boundary of the cloak. By comparing the enlarged views in Fig. 3(b) and Fig. 3(c), the effect of recovering background thermal field is distinct.

 

The quantitative analysis of the designed thermal cloak is shown in Fig. 3(d). The line with triangle symbols represents temperature distribution at r = R₃ without the cloak positioned, i.e., the original background temperature at r = R₃. The line with circle symbols represents temperature distribution at r = R₃ with the cloak positioned. A higher coincidence of these two lines indicates a more effective thermal cloak. The maximum deviation is 0.049% at 273.6°. Thus, the designed cloak performs effectively. In reality, if thermal conductivity of background is not very high, so long as the thermal conductivity in Z2 is 40 times larger than the background, the temperature deviation in Z1 can be no more than 0.5% compared with the anticipated. If there is no material whose thermal conductivity is 40 times larger than background, a flow field with only tangential velocity ⁶ can be introduced into Z2. Only if the velocity is large enough, the efficient thermal conductivity of Z2 approaches to infinite.

Fig. 3 Design result and simulated temperature field of the thermal cloak designed by the proposed method. The objective function equals to 8.52×10^-7. (a) Designed thermal conductivity in Z3. (b) Temperature field as Z2 is positioned. (c) Temperature field as Z2 and Z3 are positioned. (d) Comparison of temperatures at r = R₃ with and without cloak.

 

The objective function is used to monitor the stabilization characteristic of the present algorithm. Fig. 4(a) shows the converging curve of the design of circular cloak. It is noted that with the increase of iteration step, the objective function decreases sharply at first and then converges to 8.5×10^-7 approximately. There is no significant oscillation in the curve. Thus, the proposed method owns good stability. In iterative process, thermal conductivity in Z3 varies continuously because the penalization parameter equals to 1. Fig. 4(b) indicates the statistical results of thermal conductivity distribution shown in Fig. 3(a). It is noted that thermal conductivity values of discrete elements mainly focus on the interval between 0.3 and 0.5, but vary continuously.

Fig.4 Analysis of stability and continuity for design of circular cloak. (a) Converging curve of the algorithm. (b) Number of elements in different thermal conductivity intervals for circular cloak.

 

In fact, there is really more than one kind of thermal conductivity distribution in Z3 for a certain cloak, and the initial value of relative density can determine the final design result. Above design result in Fig. 3(a) is obtained when the initial value of relative density equals to 1. Three other design results when the initial values of relative density are different are shown in Fig. 5. Fig. 5(a)-5(c) are designed thermal conductivity when the initial value equals to 0.9, 0.7 and 0.5 respectively. It is obvious that the three design results are different, but the maximum deviations are all no larger than 0.1%.

Fig. 5 Design results for different initial value of relative density. (a) Design result when initial value equals to 0.9. (b) Design result when initial value equals to 0.7. (c) Design result when initial value equals to 0.5.

 

4.3 Design result of the cloak with irregular geometry

Without mathematical analysis for geometric shape of cloak, the proposed method can be easily extended to designing cloaks with irregular geometry, which is difficult for coordinate transformation. To validate reliability of proposed method for irregular geometry, a heart-shaped cloak is designed.

 

In the process of reverse design, there is no need for mathematical analysis of the cloak shape, and the computational domain of the cloak with an irregular shape can be discretized into Cartesian coordinate grids. The irregular shape is marked by a series of grid points based on step approximation method. As shown in Fig. 6(a), the computational domain is evenly discretized into Cartesian coordinate grids of 200×200. It should be noted that Fig. 6(a) is just an illustration of the discretized pattern and the amount of grids is not real. The bold dots indicate Z3. The design process is the same as that in part 4.1. The iteration is 16 and the convergence time is 34.9 seconds. Thermal conductivity distribution in Z3 is shown in Fig. 6(b). Anisotropy is removed. Temperature field with cloak positioned is shown in Fig. 6(c). There is no temperature gradient in Z1 and no temperature distortion in the background. It can be applied to backgrounds with any uniform temperature gradient because the design process is dimensionless. The quantitative analysis is shown in Fig. 6(d). The maximum deviation is 0.86% at 92.6°.

Fig. 6 Design results of the heart-shaped cloak. The objective function equals to 6.41×10^-6. (a) Discrete form of the irregular shape. (b) Designed thermal conductivity in Z3. (c) Temperature field with cloak positioned. (d) Comparison of temperatures at outer boundary with and without cloak.

 

Fig. 7(a) shows the converging curve of the design of irregular cloak. It is noted that with the increase of iteration step, the objective function decreases sharply at first and then converges to 6×10^-6 approximately. Fig. 7(b) indicates the statistical results of thermal conductivity distribution shown in Fig. 6(b). It is noted that thermal conductivity values of discrete elements mainly focus on the interval between 0.25 and 0.45, but vary continuously.

Fig.7 Analysis of stability and continuity for design of irregular cloak. (a) Converging curve of the algorithm. (b) Number of elements in different thermal conductivity intervals for irregular cloak

 

The reverse method can also be used to design other thermal meta-devices, like thermal concentrator, thermal rotator and thermal illusion tool. Similarly as the design process of thermal cloak, physical parameters of other thermal meta-devices can be calculated by the reverse design method according to target temperature distribution, which can be obtained on the basis of functions of different thermal meta-devices.

 

5. Design of Cylindrical cloak with thermal conductivity in polynomial distribution

The cloak designed above remove anisotropy and can be applied to backgrounds with any uniform temperature gradient. However, it is unfit for other arbitrary backgrounds whose temperature gradients are not uniform, because the thermal conductivity pattern of its cross section is not symmetrical under all rotations about the centre. Thus, the distribution of thermal conductivity should be simplified. According to our research, if thermal conductivity of the cylindrical cloak only changes along the radius, no matter what kind of distribution form of the thermal conductivity is, the cloak can be applied to any arbitrary background temperature at very limited cost. It turns out to be feasible that the thermal conductivity distribution is set as a polynomial in radius. In essence, thermal conductivity distribution is exerted a constraint condition.

 

5.1 Design process for each parameter of the polynomial

Thermal conductivity form denoted by a polynomial in Z3 is expressed as

(15)

l is the total number of terms. If l is chosen to be 0, in other word, thermal conductivity in Z3 is a constant, the cloak becomes to a bilayer cloak^35,36 when thermal conductivity in Z2 equals to zero. Therefore, the bilayer cloak can be regarded as a special case of the cloak whose thermal conductivity is in polynomial distribution. The direct calculation of thermal conductivity is transformed into computing each coefficient a_s by the proposed method. The design process is similar as above but some changes are needed as follows.

The gradient of function J with respect to each a_s is expressed as

(16)

The direction of descent for each a_s is obtained by

	(17a)
	(17b)

The step size expressed by Eq. (14) is too large for this situation and the difference for each a_s cannot be reflected. Thus, the step size for each a_s is changed to

(18)

The update equation of each a_s is

(19)

 

5.2 Design result of thermal conductivity in polynomial distribution

To verify the validity of thermal conductivity with polynomial distribution, three cases are evaluated under the conditions of l = 0, l = 1 and l = 2 respectively. As shown in Fig. 1(a), R₁, R₂ and R₃ are set as 0.1, 0.13 and 0.2 respectively. Thermal conductivity in Z2 is set as infinitely great. The computational domain is evenly divided into Cartesian coordinate grids of 200×200. Step approximation method is used to mark each zone of the cloak. For l = 0, a₀ equals to 0.41, and the iteration is 21 and the convergence time is 12.5 seconds; for l = 1, a₀ equals to 0.2 and a₁ equals to 1.26, and the iteration is 22 and the convergence time is 12.8 seconds; for l = 2, a₀, a₁ and a₂ equal to 0.16, 1.12 and 2.55 respectively, and the iteration is 22 and the convergence time is 13.4 seconds. Temperature fields of these three cases are shown in Fig. 8(a), 8(b) and 8(c). There is no temperature gradient in Z1 and no distortion on the background thermal field. The quantitative analysis is shown in Fig. 8(d). The maximum deviations for l = 0, l = 1 and l = 2 are 0.75% at 331.8°, 0.27% at 333.8° and 0.21% at 333.8°, respectively. After simulation, it is noted that no matter what the value of l is, the cloak with thermal conductivity denoted by Eq. (13) can perform well, only if every parameter a_s is calculated by the proposed method.

Fig. 8 Design results of cloaks with thermal conductivity in polynomial distribution. (a) Temperature field with cloak of l = 0. The objective function equals to 5.25×10^-6. (b) Temperature field with cloak of l = 1. The objective function equals to 5.69×10^-6. (c) Temperature field with cloak of l = 2. The objective function equals to 5.64×10^-6. (d) Comparison of temperatures at r = R₃ with and without cloak.

 

Above three cases are in the context of background thermal field with uniform temperature gradient. The cloak whose thermal conductivity in Z3 is denoted by Eq. (15) can be applied to other arbitrary background thermal fields as well.

 

The original background thermal field is shown in Fig. 2(b). The cloak under situation of l = 1 is taken as an example. Temperature field with the cloak positioned is shown in Fig. 9(a). There is no temperature gradient in Z1. The background thermal field recovers approximately compared with the original one. The quantitative analysis is shown in Fig. 9(b). The deviation is no larger than 5.00%. The smaller Z3, the better the cloak performs⁴⁴. That is, if the value of R₃ goes down, the cloak can perform better.

 

When thermal conductivity in Z3 is constant, the cloak may seem very simple. However, once the cloak sizes, i.e., R₁, R₂ and R₃, are determined, thermal conductivity in Z3 cannot be changed. If the material with thermal conductivity needed in Z3 cannot be made, the realization of cloak is a failure. Therefore, thermal conductivity in polynomial form can improve the flexibility of engineering realization.

Fig. 9 Temperature field with the cloak applied to background of non-uniform temperature gradient. (a) Temperature field with the cloak positioned. (b) Comparison of temperatures at r = R₃ with and without cloak.

 

The proposed method can be validated by the numerical simulation and experimental study in Ref. [35], in which the bilayer cloak can be deemed to be a special case of the polynomial cloak with l = 0. According to Ref. [35], R₁, R₂ and R₃ were set as 0.133, 0.211 and 0.267 respectively. Thermal conductivity in Z2 was set as 0 and thermal conductivity in Z3 equalled to 4.326. As a contrast, thermal conductivity in Z3 calculated by the proposed method is 4.321. Temperature field with cloak designed by the proposed method is shown in Fig. 10(a). The numerical and experimental results of the bilayer cloak are shown in Fig. 10(b) and Fig. 10(c). From the comparison, it can be reflected that the cloak designed by the proposed method can perform well.

Fig. 10 Comparison of the cloak designed by proposed method and the bilayer cloak. (a) Temperature field with the cloak designed by proposed method. (b) Temperature field with the bilayer cloak. (c) Experimental result of the bilayer cloak.

 

6. Dimension design of cloak in case of known material thermal conductivity

Cloak whose thermal conductivity in Z3 is constant is easily to fabricate. However, such cloak is limited by the cloak size. In other word, satisfactory material whose parameter is designed from cloak dimension is usually difficult to synthesize. The cloak dimension can be calculated at the knowledge of the thermal conductivities of existing materials by means of the proposed method. Thus, a cloak can be realized by existing materials.

 

Different from cloak in Fig. 1(b), Z3 can be substituted by available multilayer materials. Z3 can be divided into several layers, and in each layer, thermal conductivity is known. As in Fig. 11(a), Z3 is divided into two parts named Z3-1 and Z3-2. Thermal conductivities of these two parts are constant K₁ and K₂ respectively, in other words, l = 0 is suitable for each layer. The cloak design process is changed from thermal conductivity to the cloak size R₄. The ideal thermal conductivity of existing multilayer cloak is designed by coordinate transformation, and then is discretized into multilayer whose thermal conductivity is constant, based on effective medium theory. Different from existing ones, existing material is adopted and dimension of each layer is calculated directly by the reverse design method.

Fig. 11 Thermal cloak model with dimension to be determined and dimension design results. (a) Thermal cloak model with dimension to be determined. (b) Temperature field with the cloak positioned. (c) Comparison of temperature at r = R3 with and without cloak.

The gradient of J is expressed as

(20)

δR₄ is a small change of R₄. The direction of descent P is obtained by

	(21a)
	(21b)

The step size β is expressed as

(22)

R₄ is updated by

(23)

The steps determining R₄ are as follows.

Step 1: Calculate the original background thermal field U(X, Y).

Step 2: Choose the initial value of R₄.

Step 3: Solve the direct problem to obtain T (R₄).

Step 4: Add a small number δR₄ to R₄ and solve the direct problem again to obtain T (R₄+δR₄).

Step 5: Calculate the gradient of function J in Eq. (20).

Step 6: Calculate the direction of descent P in Eq. (21a) and Eq. (21b).

Step 7: Calculate the step size β in Eq. (22).

Step 8: Update R₄ by Eq. (23).

Step 9: Calculate function J in Eq. (6). If J is small enough, stop the iteration and output R₄. If not, go back to step 3 and go on this cycle.

 

In numerical simulation, K₁ is set as 0.6 and K₂ is set as 0.3. R₁, R₂ and R₃ are set as 0.1, 0.13 and 0.2 respectively. Thermal conductivities of cloak zone Z3-1 and Z3-2 are known, and the radius R₄ of interface between Z3-1 and Z3-2 should be calculated by proposed method. The calculating domain is evenly divided into Cartesian coordinate grids of 200×200. Step approximation method is used to mark each zone of the cloak. After the design process, R₄ equals to 0.165. The iteration is 5 and the convergence time is 1.8 seconds. Temperature field with the cloak positioned is shown in Fig. 11(b). The quantitative analysis is shown in Fig. 11(c). The maximum deviation is 0.59% at 3.6°. In general, Z3 can be divided into more parts and thermal conductivity of each part is known. Size of each part can be designed by the proposed method. This multilayer cloak is also applicable for arbitrary background thermal fields.

 

7. Conclusions

In this study, a reverse design method based on inverse problem theory is proposed to design thermal cloak. Compared with coordinate transformation, the anisotropy of thermal conductivity can be removed. A cloak with irregular heart shape has been designed by this method, which is difficult for coordinate transformation. Then, thermal conductivity whose distribution is in a polynomial form is designed for a cylindrical cloak with improved flexibility for arbitrary background thermal fields. The recovery extent for background thermal field with non-uniform temperature gradient is no less than 95%. At last, the dimension design of a cloak can also be realized at the knowledge of material thermal conductivity by the proposed method. In this way, the cloak can be designed by existing material to avoid complex material synthesis process. The proposed method can offer a practicable scheme for thermal cloak design to control heat conduction at will. The proposed method is based on Laplace equation, which is the governing equation of heat conduction. There is no any mathematical difference between the thermal field in this study and other Laplace fields. Therefore, this method can be extended to other Laplace fields, like direct current electric field and static magnetic field.

 

Acknowledgement

This work is supported by the Basic Science Center Program for Ordered Energy Conversion of the National Natural Science Foundation of China (No. 51888103), National Natural Science Foundation of China (No. 51906189), Foundation for Innovative Research Groups of the National Natural Science Foundation of China (No. 51721004).

 

Conflict of interest

There are no conflicts to declare.

 

References

1T. Yang, X. Bai, D. Gao, L. Wu, B. Li, J. T. L. Thong and C. Qiu, Adv. Mater. 2015, 27, 7752-7758.

2. J. Wang, J. Shang and J. Huang, Phys. Rev. Appl. 2019, 11, 024053

3. T. Han, J. Zhao, T. Yuan, D. Lei, B. Li and C. Qiu, Energy Environ. Sci. 2013, 6, 3537-3541.

4. X. Ni, Z. Wong, M. Mrejen, Y. Wang and X. Zhang, Science 2015, 349, 1310-1314.

 5. D. Schurig, J. J. Mock, B. J. Justice, S. A. Cummer, J. B. Pendry, A. F. Starr and D. R. Smith, Science 2006, 314, 977-980.

6. Y. Li, K. Zhu, Y. Peng, W. Li , T. Yang, H. Xu, H. Chen, X. Zhu, S. Fan and C. Qiu, Nat. Mater. 2019, 18, 48-54.

7. R. Hu, S. Zhou, Y. Li, D. Lei, X. Luo and C. Qiu, Adv. Mater. 2018, 30, 1707237.

8. R. Hu, S. Huang, M. Wang, X. Luo, J. Shiomi and C. Qiu, Adv. Mater. 2019, 31, 1807849.

9. M. Farhat, S. Guenneau and S. Enoch, Phys. Rev. Lett. 2009, 103, 024301.

10. M. Brun, S. Guenneau and A. B. Movchan, Appl. Phys. Lett. 2009, 94, 061903.

11. S. Zhang, C. Xia and N. Fang, Phys. Rev. Lett. 2011, 106, 024301.

12. B. I. Popa, L. Zigoneanu and S. A. Cummer, Phys. Rev. Lett. 2011, 106, 253901.

13. S. Guenneau and T. M. Puvirajesinghe, J. R. Soc. Interface 2013, 10, 20130106.

14. S. Guenneau, D. Petiteau, M. Zerrad, C. Amra and T. Puvirajesinghe, AIP Adv. 2015, 5, 053404.

15. C. Fan, Y. Gao and J. Huang, Appl. Phys. Lett. 2008, 92, 251907.

16. Y. Li, X. Bai, T. Yang, H. Luo and C. Qiu, Nat. Commun. 2018, 9, 273.

17. J. Huang, ES Energy Environ. 2019, 6, 1-3.

18. J. Guo and Z. Qu, Int. J. Heat Mass Tran. 2018, 127, 1212-1222.

19. D. M. Nguyen, H. Xu, Y. Zhang and B. Zhang, Appl. Phys. Lett. 2015, 107, 121901.

20. T. Han, X. Bai, J. T. L. Thong, B. Li and C. Qiu, Adv. Mat. 2014, 26, 1731-1734.

21. L. Xu, S. Yang and J. Huang, Phys. Rev. E 2019, 99, 022107.

22. L. Xu, S. Yang and J. Huang, Phys. Rev. Appl. 2019, 11, 054071

23. X. Peng and R. Hu, ES Energy Environ. 2019, 6, 39-44.

24. F. Yang, L. Xu and J. Huang, ES Energy Environ. 2019, 6, 45-50.

25. Z. Zhou, X. Shen, C. Fang and J. Huang, ES Energy Environ. 2019, 6, 85-91.

26. S. Huang, J. Zhang, M. Wang, W. Lan, R. Hu and X. Luo, ES Energy Environ. 2019, 6, 51-56.

27. R. Hu and X. Luo, Natl. Sci. Rev. 2019, 6, 1071-1073.

28. Y. Liu, Y. Cheng, R. Hu and X. Luo, Phys. Lett. A 2019, 383, 2296-2301.

29. L. Zhou, S. Huang, M. Wang, R. Hu and X. Luo, Phys. Lett. A 2019, 383, 759-763.

30. J. B. Pendry, D. Schurig and D. R. Smith, Science 2006, 312, 1780-1782.

31. S. Guenneau, C. Amra and D. Veynante, Opt. Express 2012, 20, 8207-8218.

32. Y. Li, X. Shen, Z. Wu, J. Huang, Y. Chen, Y. Ni and J. Huang, Phys. Rev. Lett. 2015, 115, 195503.

33. R. Wang, L. Xu, Q. Ji and J. Huang, J. Appl. Phys 2018, 123, 115117.

34. L. Xu, S. Yang and J. Huang, Phys. Rev. Appl. 2019, 11, 034056.

35. T. Han, X. Bai, D. Gao, J. T. L. Thong, B. Li and C. Qiu, Phys. Rev. Lett. 2014, 112, 054302.

36. H. Xu, X. Shi, F. Gao, H. Sun and B. Zhang, Phys. Rev. Lett. 2014, 112, 054301.

37. Y. Liu, W. Guo and T. Han, Int. J. Heat Mass Tran. 2017, 115, 1-5.

38. J. Qin, W. Luo, P. Yang, B. Wang, T. Deng and T. Han, Int. J. Heat Mass Tran. 2019, 141, 487-490.

39. X. He and L. Wu, Appl. Phys. Lett. 2014, 105, 221904.

40. F. Mohebbi, M. Sellier and T. Rabczuk, Int. J. of Therm. Sci. 2017, 117, 68-76.

41. C. Huang and J. Yan, Int. J. Heat Mass Tran. 1995, 38, 3433-3441.

42. G. Fujii, Y. Akimoto and M. Takahashi, Appl. Phys. Lett. 2018, 112, 061108.

43. G. Fujii and Y. Akimoto, Int. J. Heat Mass Tran. 2019, 137, 1312-1322.

44. F. Gomory, M. Solovyov, J. Souc, C. Navau, J. Prat-Camps and A. Sanchez, Science 2012, 335, 1466-1468.