Metamaterials can be used to reveal novel physics and develop new applications because of their novel artificial structures. The transformation principle can be used to control processes with governing equations that are form-invariant under coordinate transformations, enabling the control of diffusion processes by changing material parameters such as conductivity. Thermal metamaterials, together with the governing theories, make it possible to actively manipulate macroscopic heat phenomena of artificial systems, which enables people to change the heat phenomena. Such metamaterials or metamaterial-based devices refer to those artificial structures that yield novel functions in controlling heat transfer [1, 2].

The thermal invisibility cloak has a unique thermal insulation function, and its temperature gradient in the thermal stealth area is zero, which can avoid the heat flow very well. Therefore, thermal invisibility cloaks have a great application value in computer chips, spacecraft return modules, satellites, and other equipment.

In 2006, Pendry et al. [3], based on the form invariance of Maxwell’s equation, designed an artificial metamaterial with anisotropic characteristic parameters to make electromagnetic waves propagate along a preset path, to achieve the purpose of electromagnetic invisibility. In 2008, Fan et al. [4] first introduced the concept of an electromagnetic invisibility cloak into the field of thermal science, and theoretically predicted the thermal invisibility cloak. In 2010, Li et al. [5] proposed the principle of using coordinate transformation to design multi-functional stealth cloak. The multi-functional invisibility cloak can not only realize the function of a thermal invisibility cloak, but also realize the function of an electric invisibility cloak. In 2012, Narayana et al. [6] used COMSOL software to simulate the heat flow distribution of two different materials coupled together. They built a prototype thermal cloak in the lab by combining two materials with vastly different thermal conductivities into a spiral-shaped multilayer structure. By measuring its temperature distribution in the temperature field, it is verified that its experimental performance is in complete agreement with the theoretical expectation, thus opening the door for experimental exploration of a thermal invisibility cloak. In 2013, Schittny et al. [7] constructed a new type of heat invisibility cloak using a single copper material. By machining thermal channels of different sizes on the copper sheet (each channel leads to different thermal conductivity due to different sizes), the heat flow in r<R2 is compressed to R1<r<R2. In 2014, Mao et al. [8] designed the expression of the thermal conductivity of the cylindrical heat cloak with non-conformal arbitrary cross-section based on the transformation of thermodynamics and verified that the heat cloak has the function of thermal protection for the stealth area through COMSOL Multiphysics simulation. In 2015, Li et al. [9] introduced the dependence of thermal conductivity on temperature when studying thermal superstructural materials, thus developing the transformation thermal theory (thermal conductivity is independent of temperature) and establishing the nonlinear transformation thermal theory, which then promoted the nonlinear thermal research based on thermal superstructural materials. The so-called nonlinear heat theory focuses on the thermal conductivity of materials, which is no longer a constant independent of temperature, but a physical quantity with specific response to temperature. Based on this response, they designed switching thermal cloaks that can respond intelligently to ambient temperatures and prepared macroscopic thermal diodes. In 2017, Xia et al. [10] deduced the general solution expression of the thermal conductivity of the three-dimensional heat cloak with arbitrary shape and carried out full wave simulation verification. In 2018, Xia et al. [11] deduced the expression of the thermal conductivity of the square heat cloak and designed the square heat cloak with any angle. In 2021, Sha et al. [12] address three long-standing challenges, i.e., transformation optics-induced anisotropic material parameters, the limited shape adaptability of experimental thermal meta devices, and a priori knowledge of back-ground temperatures and thermal functionalities. They took the local thermal conductivity tensors as input, resorted to topology optimization for the freeform designs of topological functional cell (TFCs), and then directly assembled and printed them. They designed and 3D-printed three freeform thermal meta devices (concentrator, rotator, and cloak).

The structural vibro-acoustic problem in a thermal environment [13, 14] is one of the most common problems in aerospace engineering. This paper studies the problem of heat avoidance in thermal environment, aiming at laying a foundation for solving the problem of thermal acoustic vibration. Previous studies have focused on a single cloak and did not propose how to arrange the cloaks periodically to study the manipulation of the heat flow by the periodic plate of the heat cloak. In this paper, two specific ways are given to obtain the transformation equation of the square heat cloak. In addition, the " position parameters " are used in the thermal conductivity of the cloak, to obtain the thermal cloak at any position in practical applications. Based on this, the thermal cloaks are arranged periodically. The regulating effect of the periodic plate of the thermal cloaks on the heat flow is simulated under different boundary conditions by using COMSOL software.

In nature, the heat flux always flows from hot areas to cold areas, which is decided by the heat conduction equation: ρ · c ∂ T ∂ t = ∇ · λ ∇ T + Q (1) where ρ , c , λ represent the density, thermal capacity, and the thermal conductivity of the medium respectively. T is the temperature. Q is the heat resource. If there is no heat resource and at steady state, the Eq. 1 is simplified as: ∇ · λ ∇ T = 0 (2)

Since the Eq. 2 has formal invariance [15], it could be transformed in another space as follows: ∇ ′ · λ ′ ∇ ′ T ′ = 0 (3) where λ ′ , T ′ represent the thermal conductivity coefficients and temperatures in transformed spaces. According to the theory of transformation thermodynamics, the thermal conductivity relationship between the transformation space and the original space [16]: λ ′ = A · λ 0 · A T det A (4) where λ 0 represents the thermal conductivity coefficients in the original space. A is the Jacobian transformation matrix, reflecting the geometric changes from the original space to the transformation space, and its components are: A i j = ∂ x i ′ ∂ x j (5) where x i ′ denotes the three coordinate components x’, y’, z’ in the new coordinate system. x j represents the three coordinate components x, y, z in the original coordinate system.

The square thermal cloak is designed as the Figure 1 shown. A square area with a side length of 2 s 2 in Figure 1A is compressed into a square torus in Figure 1B which is encircled by a square with side length of 2 s 1 and a square with side length of 2 s 2 . Considering that the region has geometric symmetry, the region area could be divided into four triangular areas along its diagonal. The thermal conductivity is calculated according to the transformation equation of four triangular regions respectively [17].

Two approaches will be given to find the thermal conductivities of the square thermal cloak as follows. The first one is a direct method. To illustrate how the transformation equation of each region is obtained, consider the following case Figure 2A: take a point x on a line segment, the length of which is l. When the segment is compressed to a length of l ′ , as Figure 2B shows, the corresponding x becomes x ′ . It is not difficult to get the following formula: x l = x ′ l ′ (6)

Figures 3A, B are the schematic diagram of the area 1 before and after the transformation respectively. Take the left triangle area’s median line, as Figure 4A shows, for example. For the first step shifting the whole segment to the right by s 1 (Figure 4B): x 1 = x + s 1 (7)

Second, compress the line segment into s 2 − s 1 (Figure 4C): x ′ = s 1 + Δ x (8)

According to formula 6, it is obtained: Δ x s 2 − s 1 = x s 2 (9)

Simultaneous Eqs 7–9, it could be obtained: x ′ = s 2 − s 1 s 2 x + s 1 (10)

It can be known from similar triangles: x x ′ = y y ′ , y ′ = s 2 − s 1 s 2 + s 1 x y (11)

From the above example, it is easy to find that the transformation formula of area 1 is as follows: x ′ = x s 2 − s 1 s 2 + s 1 , y ′ = s 2 − s 1 s 2 + s 1 x y (12)

In area 1, the Jacobian matrix A 1 and the thermal conductivity λ 1 ′ of the transformation space are: A 1 = s 2 − s 1 s 2 0 − s 1 y x 2 s 2 − s 1 s 2 + s 1 x , λ 1 ′ = 1 − s 1 x ′ − s 1 y ′ x ′ 2 − s 1 y ′ x ′ 2 x ′ 4 + s 1 2 y ′ 2 x ′ 4 − s 1 x ′ 3 λ 0 (13)

In a similar way, the transformation formula and the Jacobian matrix of areas 2, 3, and 4 could be found as the formulae 14, 15, and 16: x ′ = x s 2 − s 1 s 2 + s 1 y , y ′ = y s 2 − s 1 s 2 + s 1 (14) x ′ = x s 2 − s 1 s 2 − s 1 , y ′ = y s 2 − s 1 s 2 − s 1 x (15) x ′ = x s 2 − s 1 s 2 − s 1 y , y ′ = y s 2 − s 1 s 2 − s 1 (16)

The corresponding Jacobian matrixes of areas 2, 3, and 4 are respectively as the formulae 17, 18, and 19: A 2 = s 2 − s 1 s 2 + s 1 y − s 1 x y 2 0 s 2 − s 1 s 2 (17) A 3 = s 2 − s 1 s 2 0 s 1 y x 2 s 2 − s 1 s 2 − s 1 x (18) A 4 = s 2 − s 1 s 2 − s 1 y s 1 x y 2 0 s 2 − s 1 s 2 (19)

Finally, the thermal conductivities λ 2 ′ , λ 3 ′ , λ 4 ′ of each transformation area 2,3,4 are respectively as follows: λ 2 ′ = y ′ 4 + s 1 2 x ′ 2 y ′ 4 − s 1 y ′ 3 − s 1 x ′ y ′ 2 − s 1 x ′ y ′ 2 1 − s 1 y ′ λ 0 (20) λ 3 ′ = 1 + s 1 x ′ s 1 y ′ x ′ 2 s 1 y ′ x ′ 2 x ′ 4 + s 1 2 y ′ 2 x ′ 4 + s 1 x ′ 3 λ 0 (21) λ 4 ′ = y ′ 4 + s 1 2 x ′ 2 y ′ 4 + s 1 y ′ 3 s 1 x ′ y ′ 2 s 1 x ′ y ′ 2 1 + s 1 y ′ λ 0 (22)

The above is the direct way to get the thermal conductivity of the square thermal cloak.

The second way is using a rotation matrix to figure out the thermal conductivity of the four areas. The areas 2, 3, and 4 can be obtained by rotating the area1 counterclockwise by 90°, 180°, 270°. The rotation process can be implemented by the rotation matrix.

To begin with, let us introduce the rotation matrix P by the flowing example. Consider the situation shown in Figure 5: rotate point A 90° counterclockwise to get point B.

The coordinates of point A are: x = ρ · cos ⁡ θ , y = ρ · sin ⁡ θ (23) where x and y indicate the horizontal and vertical coordinates of point A in the Cartesian coordinate system. ρ, θ are the polar diameter and polar angle of point A in the polar coordinate system.

The coordinates of point B are: s = ρ · cos θ + π 2 = − ρ · sin ⁡ θ = − y , t = ρ · sin θ + π 2 = ρ · cos ⁡ θ = x (24) where s and t indicate the horizontal and vertical coordinates of point B in the Cartesian coordinate system.

The matrix form is used to clarify the relationship between A and B: s t = 0 − 1 1 0 x y (25) P = 0 − 1 1 0 (26) s t = P x y (27)

P is the rotation matrix, which mean a transformation of rotating counterclockwise by 90°. For the area 1, the transformation formula is: x ′ y ′ = s 2 − s 1 s 2 0 0 s 2 − s 1 s 2 + s 1 x x y + s 1 0 (28)

The transformation equation of area 2 can be obtained by rotating the transformation equation of area 1 90° counterclockwise. Multiply both sides of formula 28 by the rotation matrix P: P x ′ y ′ = P s 2 − s 1 s 2 0 0 s 2 − s 1 s 2 + s 1 x P − 1 P x y + P s 1 0 (29)

B and B’ are obtained by rotating A and A’, so we can get: s ′ t ′ = P x ′ y ′ , s t = P x y (30)

Take the formula 30 into the formula 29: s ′ t ′ = P s 2 − s 1 s 2 0 0 s 2 − s 1 s 2 + s 1 t P − 1 s t + 0 s 1 = s 2 − s 1 s 2 + s 1 t s s 2 − s 1 s 2 t + s 1 (31)

Let s ′ , t ′ , s , t in the formula 31 be represented by x ′ , y ′ , x , y respectively, the transformation equation of area 2 is obtained as follows: x ′ = x s 2 − s 1 s 2 + s 1 y , y ′ = y s 2 − s 1 s 2 + s 1 (32)

In a similar way, by multiplying the rotation matrix P, the transformation equations of areas 3 and 4 can be found. Then, the thermal conductivity of each transformation space is calculated by formula 4. The results of the direct method and the rotation matrix method are the same. The two methods’ principles are geometric transformation and matrix rotation respectively. The physical natures are the same, only the mathematical means are different. The cloud map of the thermal conductivity coefficient of the thermal cloak is drawn in Figure 6. Figure 6A is λ x x ′ and Figure 6B is λ x y ′ .

The above thermal cloak is at the origin of the frame. Furthermore, the thermal conductivity of a thermal cloak at any position is calculated by introducing position parameters a, b. Move the coordinate system: m = x + a , n = y + b , m ′ = x ′ + a , n ′ = y ′ + b (33)

The transformation equations of area 1 are: m ′ = m − a s 2 − s 1 s 2 + s 1 + a , n ′ = n − b s 2 − s 1 s 2 + s 1 x − a + b (34)

Let m ′ , n ′ , m , n in the formula 34 be represented by x ′ , y ′ , x , y respectively, the transformation equation of area 1 is obtained as follows: x ′ = x − a s 2 − s 1 s 2 + s 1 + a , y ′ = y − b s 2 − s 1 s 2 + s 1 x − a + b (35)

The corresponding Jacobian matrix and the thermal conductivity of area 1 at any position are: A 1 = s 2 − s 1 s 2 0 − s 1 y − b x − a 2 s 2 − s 1 s 2 + s 1 x − a , λ 1 ′ = 1 − s 1 x ′ − a − s 1 y ′ − b x ′ − a 2 − s 1 y ′ − b x ′ − a 2 x ′ − a 4 + s 1 2 y ′ − b 2 x ′ − a 4 − s 1 x ′ − a 3 λ 0 (36)

In a similar way, in the areas 2, 3, and 4, the thermal conductivity of a thermal cloak at any position can be obtained.

To verify the theoretical results of the above derivation, we use COMSOL to implement the simulation. The square-shaped thermal cloak is built in the heat transfer module. Set s 1 = 3 m , s 2 = 6 m . The thermal conductivity of the original space λ 0 = 1 w / m · k . Set the thermal conductivity calculated above in the material setup module. Two kinds of thermal boundaries are set for the simulation of one thermal cloak cell and the periodic plate.

Two methods have been given to construct the thermal conductivities of the square thermal invisibility cloak. The cloud map of the thermal conductivity coefficient of the thermal cloak has been drawn. Besides, by introducing the position parameters into the calculation of the thermal conductivity coefficient, the cloaks were arranged periodically so that they can manipulate heat flow in a larger area. The heat insulation function of both the single thermal cloak and the thermal cloak periodic plate were tested under different heat boundary conditions using COMSOL Multiphysics.

The results show that both the single thermal and the thermal cloak periodic plate have the function of avoiding heat flow under different heat boundary conditions. The heat fluxes travel around the inner domain with good thermal stealth effect. To prove the necessity of the position parameters, the comparative verification is done. The results show that the heat flow can stably bypass each thermal invisibility cloak only when the position parameters are considered.

In engineering application, vibration reduction and heat insulation have always been the focus of research. Phonon crystal plates can be used to reduce response. How to combine the phonon crystal plate with the thermal cloak periodic plate to play the role of shock absorption and heat insulation at the same time is worthy of further study.