Introduction

In recent years, thermal metamaterials have attracted much attention from researchers due to their extraordinary ability on the manipulation of heat flow1,2,3,4, which benefits from the spatial distribution of specific thermal conductivities in a given thermal metadevice5,6,7,8. According to Fourier’s law9, thermal conductivity κ is a tensor matrix (e.g., \(\kappa = \left[ {\begin{array}{*{20}{c}} {\kappa _{xx}} & {\kappa _{xy}} \\ {\kappa _{yx}} & {\kappa _{yy}} \end{array}} \right],\left\{ {\begin{array}{*{20}{l}} {\kappa _{xx}\kappa _{yy} \,>\, \kappa _{xy}\kappa _{yx}} \\ {\kappa _{xx,yy} \,>\, 0,\;\kappa _{xy} = \kappa _{yx}} \end{array}} \right.\) for the 2D case10,11) that characterizes the general linear heat conduction behavior between heat flux and temperature gradients, i.e., \({{{\mathbf{q}}}} = \left[ {\begin{array}{*{20}{c}} {q_x} \\ {q_y} \end{array}} \right] = - \kappa \cdot \nabla {{{\mathbf{T}}}} = - \left[ {\begin{array}{*{20}{c}} {\kappa _{xx}} & {\kappa _{xy}} \\ {\kappa _{yx}} & {\kappa _{yy}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\partial x/\partial T} \\ {\partial T/\partial y} \end{array}} \right]\). For isotropic materials, κxy = κyx = 0 and κxx = κyy, which make identical thermal conductive behaviors in all directions and present some limitations on manipulating heat flow. In contrast, anisotropic materials can flexibly manipulate heat flow. Unfortunately, the vast majority of natural materials are isotropic, which naturally motivates the use of mixtures of natural materials to obtain anisotropic effective thermal conductivities (ETCs).

For anisotropic ETCs, they can be divided into two categories by whether κxy is equal to 0. When κxy = 0, they have been frequently studied in orthotropic composite materials12,13. Considering achieving an arbitrary ETC with κxy = 0, the combination of natural materials based on Effective Medium Theory (EMT) for predicting properties of mixtures are widely employed, such as Wiener bounds13,14 (WB), Maxwell–Euken model14,15, Bruggeman model16, and Asymptotic Homogenization Method (AHM)17. As shown in Fig. 1a, for two given materials (e.g., M1 and M4), it has been demonstrated that WB achieved by layered structures determine the lower and upper bounds of ETCs, and the orange area in between represents the entire ETCs which can be achieved by proper mixtures of M1 and M4. To design a mixture configuration with the targeted anisotropic ETC inside the orange area of Fig. 1a, we can directly resort to designing mixtures of M1 and M4 along the black dashed curve based on EMT in Fig. 1a, which is efficiency-limited and case-dependent. Alternatively, we can start with a pre-determined structure shape and then find proper materials, like the layered structures with M2 and M5 lying on the WB along the blue curves in Fig. 1a. However, the selected M2 and M5 may not exist in nature, and even if we can find two proper naturally-occurring thermal conductive materials, their contact thermal resistance may be an issue. Using the aforementioned strategies, the mixtures with targeted ETC can only be designed along the curves as seen in Fig. 1a, leaving ample space rarely explored. For anisotropic ETCs with κxy ≠ 0, as illustrated in Fig. 1c, we take κxx, κyy and κxy (κyx) as three orthogonal axes to construct a 3D coordinate system Φ(κxx, κyy, κxy) \({\mathbb{R}}^3\) and use Ω(Ma, Mb) to denote the full-parameter anisotropic space of thermal conductivity with just two given materials Ma and Mb, e.g., Ω(M1, M4). For ETCs with κxy ≠ 0 in Ω(M1, M4), generally they can be obtained by the rotation of the mixtures with ETCs in Ω(M1, M4, κxy = 0), i.e., the orange area in Fig. 1a (see Supplementary Note 1 for details about the rotation of the mixtures). Therefore, using the same strategies mentioned above, the mixture with a targeted ETC in Ω(M1, M4) can only be designed along some surfaces, inevitably leaving more ample space rarely explored. Here, we propose a universal and efficient strategy by employing structural topology optimization18,19,20,21 to obtain the optimized material layout of a mixture with desired ETC, which is dubbed topological functional cell (TFC) hereinafter (seen in Fig. 1a–c), and ETCs could easily traverse the full-parameter anisotropic space of thermal conductivity by the judiciously designed TFCs.

Fig. 1: Schematic of full-parameter anisotropic space of thermal conductivity and four thermal metadevices.
figure 1

a Full-parameter anisotropic space with κxy = 0 and the purple dash-dotted line indicates isotropic ETCs. A targeted ETC inside Ω(M1, M4, κxy = 0) can be achieved by the mixture of materials M1 and M4 based on EMT (along the black dashed curve). Alternatively, we can resort to selecting proper materials, i.e., M2 and M5, based on WB (along the blue curves). In contrast, all ETC s inside the orange region can be easily achieved via TFCs. The legend on the right illustrates different mixture configurations. The subscript denotes the materials whose thermal conductivities are given at the top of the figure. b Illustration of the full-parameter anisotropic space with κxy = 0 achievable by TFCs with different materials. c Full-parameter anisotropic space with the rotation of TFCs (RTFCs). The orange liked-ellipsoid space is the full-parameter anisotropic one Ω(M1, M2) of thermal conductivity. The inset shows a side view of Ω(M1, M2). d Schematic of ideal temperature field for four thermal metadevices, including thermal connector, reflector, concentrator and cloak, which is designed based on the full-parameter anisotropic space. Black arrows indicate the direction of heat flow, and the colors indicate the distribution of temperatures.

In this paper, we successfully traverse the full-parameter anisotropic space Ω(Copper, PDMS) of thermal conductivity by TFCs with different ETCs. Then, the thermal properties of several typical TFCs inside Ω(Copper, PDMS) are validated numerically and experimentally. We propose the regionalized scattering cancellation method to design a series of thermal metadevices, including thermal connector, reflector, concentrator and cloak. The temperature field distributions of the four thermal metadevices are schematically shown in Fig. 1d. Finally, we fabricate them by 3D printing and experimentally characterize their thermal functionalities. This study offers an avenue for traversing the full-parameter anisotropic space with topology-optimized thermal metamaterials under two given natural materials, hence enabling highly robust manipulation of heat flow and achieving the powerful design capability for thermal metadevices.

Results

Traversal of full-parameter anisotropic space and manipulation of heat flow

Naturally-occurring materials usually have specific and isotropic thermal conductivities22, forming a discrete database of thermal conductivity values only located on the purple dash-dotted line in Fig. 1a–c. In contrast, the ETC values of mixtures with naturally-occurring materials can occupy a large space bounded by WB, i.e., the full-parameter anisotropic space of thermal conductivity. Here, we start with consideration of how much space can be occupied by the ETC of a mixture with two given materials Ma and Mb. For simplicity, we consider the 2D-case thermal conductivity tensor \({{{\mathbf{\kappa }}}}{{{\mathrm{ = }}}}\left[ {\begin{array}{*{20}{c}} {\kappa _{xx}} & {\kappa _{xy}} \\ {\kappa _{yx}} & {\kappa _{yy}} \end{array}} \right]\). It is known that the ETC of a mixture based on WB, i.e., Series and Parallel models, can be quantified as \({{{\mathbf{\kappa }}}}_s = \left[ {\begin{array}{*{20}{c}} {\kappa _y} & 0 \\ 0 & {\kappa _x} \end{array}} \right]\) or \( {{{\mathbf{\kappa }}}}_p = \left[ {\begin{array}{*{20}{c}} {\kappa _x} & 0 \\ 0 & {\kappa _y} \end{array}} \right]\), where κx = faκa + fbκb and κy = 1/(faa + fbb), κa and κb are respectively the thermal conductivities of two given materials, and fa and fb are their corresponding volume fractions with fa + fb = 1. When fa increases from 0 to 1, the ETC of this mixture will vary along the upper or lower bounds. Thus, as seen in Fig. 1a, the ETCs of mixtures with M1 and M4 occupy the orange area Ω(M1, M4, κxy = 0). Then, considering the rotation of the mixtures, the 3D orange surface-bounded space Ω(M1, M4) in Fig. 1c will be occupied by the ETCs of mixtures with M1 and M4, i.e., the full-parameter anisotropic space. Specifically, the corresponding relationship between the ETC \(\left[ {\begin{array}{*{20}{c}} {\kappa _{xx}} & {\kappa _{xy}} \\ {\kappa _{yx}} & {\kappa _{yy}} \end{array}} \right]\) in Ω(M1, M4) and the ETC \(\left[ {\begin{array}{*{20}{c}} {\kappa _x} & 0 \\ 0 & {\kappa _y} \end{array}} \right]\) in Ω(M1, M4, κxy = 0) can be expressed as

$$\left[ {\begin{array}{*{20}{c}} {\kappa _{xx}} & {\kappa _{xy}} \\ {\kappa _{yx}} & {\kappa _{yy}} \end{array}} \right] = {{{\mathbf{R}}}}^{{{\mathrm{T}}}}(\theta )\left[ {\begin{array}{*{20}{c}} {\kappa _x} & 0 \\ 0 & {\kappa _y} \end{array}} \right]{{{\mathbf{R}}}}(\theta ),\;\kappa _{x,y} = \left( {\kappa _{xx} + \kappa _{yy} \pm \sqrt {\left( {\kappa _{xx} - \kappa _{yy}} \right)^2 + 4\kappa _{xy}^2} } \right)/2$$
(1)

where R(θ) is the rotated matrix, and θ is the rotated angle. Thus, the ETC \(\left[ {\begin{array}{*{20}{c}} {\kappa _{xx}} & {\kappa _{xy}} \\ {\kappa _{yx}} & {\kappa _{yy}} \end{array}} \right]\) in Ω(M1, M4) and the ETC \(\left[ {\begin{array}{*{20}{c}} {\kappa _x} & 0 \\ 0 & {\kappa _y} \end{array}} \right]\) in Ω(M1, M4, κxy = 0) can be convertible by rotating the equivalent structure with an angle θ. Note that the symmetry axis of Ω(M1, M4) is the purple dash-dotted line with κxx = κyy and κxy = 0 in Fig. 1c, composed of the points with isotropic ETCs. Space Ω(M1, M4) is not generated by the rotation of the plane space Ω(M1, M4, κxy = 0) around the symmetry axis (purple dash-dotted line in Fig. 1c) because the cross-section of Ω(M1, M4) perpendicular to the symmetry axis is proved to be an ellipse mathematically, as shown in the inset of Fig. 1c (see Supplementary Note 1 for the rigorous derivation). As shown in Fig. 1b, we can obtain different full-parameter anisotropic space when changing the thermal conductivity of one material. It is concluded intuitively that when keeping one material with the high thermal conductivity unchanged, the other material with a lower thermal conductivity will create a larger full-parameter space. Additionally, the lower thermal conductivity of the material will narrow the gap between the Ω and physical limit of Second law of thermodynamics10,11 (see Supplementary Fig. 2 for more details).

As mentioned above, the full-parameter anisotropic space Ω(Ma, Mb) represents all the 2D-case ETCs that can be achieved by mixtures with materials Ma and Mb. It is known that layered structures have been frequently adopted to design thermal metamaterials due to their easy realization, such as thermal cloak1,5,6,7,23, concentrator5, rotator5, reflection24 and illusion25. However, the ETC value of the layered structures with rotation is only located on the shell boundary ∂Ω(Ma, Mb). We can anticipate that if all ETCs in Ω(Ma, Mb) can be achieved robustly, the thermal metamaterials can be freely designed, even if involving a variety of inhomogeneous anisotropic ETCs, such as those induced by coordinate transformation method4,26. Here, we employ topology optimization to design TFCs with different ETCs, so as to traverse the full-parameter anisotropic space (see Methods for detailed information about topology optimization). Taking copper and polydimethylsiloxane (PDMS) as typical materials, the TFCs with these two materials are designed. As shown in Fig. 2a, the blue points denote the ETCs of TFCs, i.e., the discrete anisotropic ETCs in Ω(Copper, PDMS, κxy = 0). They can almost fill half the Ω(Copper, PDMS, κxy = 0). According to Supplementary Note 1 and Eq. (1), when rotating these TFCs from 0 to 180 degrees, all the ETC values within Ω(Copper, PDMS) can be obtained. Therefore, the whole full-parameter anisotropic space of thermal conductivity is traversed by the ETCs of TFCs. Schematically, Fig. 2b shows three typical TFCs with anisotropic ETCs, i.e., κ1, κ2 and κ3, which are directly obtained via topology optimization. Their topology optimization iteration processes are shown in Fig. 2c and Supplementary Fig. 3, respectively. It is found that the TFCs become distinct along with the optimization iteration process, and their ETCs are close to the target ones.

Fig. 2: Traversal of full-parameter anisotropic space Ω(Copper, PDMS) by topology optimization and manipulation of heat flow.
figure 2

a Illustration of traversing full-parameter anisotropic space with TFCs. The ETCs of TFCs correspond to the blue points, and each TFC is composed of copper (yellow) and PDMS (green). b Three typical ETC tensors and the corresponding TFCs. The size of each cell is 5 mm × 5 mm. c Topology optimization process for κ1 with several intermediate designs. d Heat flux vectors across interface l(j)(j = 1) of two materials with different κ from Area i (i = 1) to Area i + 1 (i = 1). Different colors represent different areas.

Because of the powerful design ability of topology optimization, the TFCs with desired ETCs within the full-parameter anisotropic space can be designed robustly. This computational way adequately excavates the design ability of structural materials, thereby conveniently traversing the full-parameter anisotropic space. Compared with the previously designed ways, topology optimization provides an inverse design paradigm, which is universal, elegant, efficient and exact, and can be used for designing TFCs with the desired κ within Ω. According to Fourier’s law equation \({{{\mathbf{q}}}} = - {{{\mathbf{\kappa }}}} \cdot \nabla {{{\mathbf{T}}}}\), this design paradigm enables freely manipulating heat flow and conveniently designing thermal metamaterials with different functionalities. To manipulate heat flow, we first deduce the relationship between the heat flux vectors across the interface of two materials with different κ. As shown in Fig. 2d, the heat flux changes its direction when it conducts across the interface l(j) (l(j) denotes the j-th interface) from Area i (i = 1) to Area i + 1 (i = 1). It is known that the component of the heat flux vector along the normal direction of l(j) maintains the same without interfacial thermal resistance, namely

$$\begin{array}{lll}\left( {{{{\mathbf{q}}}}^{(in)}} \right)_{ \bot l^{(j)}} \!\!\!\!&=\!\!\!& \left( {{{{\mathbf{q}}}}^{(out)}} \right)_{ \bot l^{(j)}}\\ \left( {{{{\mathbf{q}}}}^{(in,out)}} \right)_{ \bot l^{(j)}} \!\!\!\!&=\!\!\!& \left| {\left( {{{{\mathbf{q}}}}^{(in,out)}} \right)_{l^{(j)}}} \right| \cdot \left( {\sin \alpha ^{(in,out)}} \right)_{l^{(j)}}\\ \left( {\alpha ^{(in)}} \right)_{l^{(j)}} \!\!\!\!&=\!\!\! &\left( \theta \right)_{l^{(j)}} - \arccos \left( {\frac{{\left( {q_x^{(in,out)}} \right)_{l^{(j)}}}}{{\left| {\left( {{{{\mathbf{q}}}}^{(in,out)}} \right)_{l^{(j)}}} \right|}}} \right)\\ \left( {\alpha ^{(out)}} \right)_{l^{(j)}} \!\!\!\!&=\!\!\! &\left( \theta \right)_{l^{(j)}} + \arccos \left( {\frac{{\left( {q_x^{(in,out)}} \right)_{l^{(j)}}}}{{\left| {\left( {{{{\mathbf{q}}}}^{(in,out)}} \right)_{l^{(j)}}} \right|}}} \right)\end{array}$$
(2)

where \(\left( {{{{\mathbf{q}}}}^{(in,out)}} \right)_{ \bot l^{(j)}}\)denotes the normal component of heat flux vector \({{{\mathbf{q}}}}^{(in,out)} = \left[ {\begin{array}{*{20}{c}} {q_x^{(in,out)}} \\ {q_y^{(in,out)}} \end{array}} \right]\) across the interface l(j), and \(\left( \theta \right)_{l^{(j)}}\) denotes the angle between the interface l(j) and the horizontal line. By substituting Eq. (2) into Fourier’s law equation, we can obtain

$$\left( {{{{\mathbf{q}}}}^{(in,out)}} \right)_{l^{(j)}} = - {{{\mathbf{\kappa }}}}^{(i,i + 1)} \cdot \nabla \left( {{{{\mathbf{T}}}}^{(i,i + 1)}} \right)_{l^{(j)}}$$
(3)

where \({{{\mathbf{\kappa }}}}^{(i)} = \left[ {\begin{array}{*{20}{c}} {\kappa _{xx}^{(i)}} & {\kappa _{xy}^{(i)}} \\ {\kappa _{yx}^{(i)}} & {\kappa _{yy}^{(i)}} \end{array}} \right]\) and T(i) represent the thermal conductivity tensor and temperature distribution of Area i, respectively. Obviously, the heat flux q(in,out), thermal conductivity tensor \({{{\mathbf{\kappa }}}}^{(i,i + 1)} = \left[ {\begin{array}{*{20}{c}} {\kappa _{xx}^{(i,i + 1)}} & {\kappa _{xy}^{(i,i + 1)}} \\ {\kappa _{yx}^{(i,i + 1)}} & {\kappa _{yy}^{(i,i + 1)}} \end{array}} \right]\), and the temperature distribution T(i,i+1) of different areas are correlated by Eq. (3). Therefore, under a given temperature gradient, q(in,out) can be controlled by adjusting \({{{\mathbf{\kappa }}}}^{(i,i + 1)} = \left[ {\begin{array}{*{20}{c}} {\kappa _{xx}^{(i,i + 1)}} & {\kappa _{xy}^{(i,i + 1)}} \\ {\kappa _{yx}^{(i,i + 1)}} & {\kappa _{yy}^{(i,i + 1)}} \end{array}} \right]\), which can be selected from the full-parameter anisotropic space of thermal conductivity.

Simulated and experimental verifications of thermal properties of TFCs

Here, the thermal properties of the three TFCs in Fig. 2b are validated by numerical simulations and experiments under the linear temperature gradient. We first consider the simulated temperature field distribution of a pure continuous plate, whose thermal conductivity is assigned with the target ETCs, i.e., κ1, κ2 and κ3, as shown in Fig. 3a–c, where materials with different anisotropies present different perturbations to the temperature field under the same thermally boundary conditions. They are regarded as references and called as theoretical simulations hereinafter. Then, the simulated temperature field distribution of the three periodic TFCs, called as structural simulations hereinafter, are respectively shown in Fig. 3d–f (see Methods for details about numerically theoretical and structural simulations). From Fig. 3a–f, the corresponding isotherms are similar under the above two simulated cases, which indicates that each temperature field under structural simulations is close to the corresponding reference one under theoretical simulations.

Fig. 3: Simulated and experimental results of three TFCs.
figure 3

ac Theoretical simulations of pure materials with target ETC tensors κ1, κ2 and κ3. The white curves are isotherms. df Structural simulations of TFCs composed of two materials (copper and PDMS). Temperature observation along the red dashed lines is performed for comparison. gi Samples of TFCs by 3D printing in experiments, not filled with PDMS and covered by polyvinylchloride (PVC). j, k Experimentally measured temperature field distribution for 10 × 10 periodic TFCs in gi. mo Measured temperature values along the observation lines under three cases. The markers are the calculated points based on the temperature data gathered by the infrared (IR) camera. The insets show the TFCs and the target ETCs.

Next, as shown in Fig. 3g–i, we further experimentally validate the thermal behaviors of the three TFCs after fabricating their samples by 3D printing, called as structural experiments hereinafter (see Methods for details of an experimental setting). As shown in Fig. 3j–l, each experimentally steady-state temperature field distribution is close to the corresponding reference one as well. To quantitatively compare the temperature field distribution, we calculate the temperature value of the points on the red dashed lines in Fig. 3d–f under the above three cases, i.e., theoretical simulation, structural simulation and structural experiment. The temperature profiles under three cases are plotted in Fig. 3m–o, and their differences are small quantitatively. Therefore, it is validated that the ETCs of TFCs are close to the target ones both theoretically and experimentally, implying the robust accuracy and traversal of our design paradigm.

Applications of full-parameter anisotropic space

We design the topology-optimized thermal metadevices to show the significance of traversing full-parameter anisotropic space. Here, we consider two types of thermal metadevices, i.e., thermal connector and reflector, which divert heat flow stably and can be used to bypass non-embedded obstacles, reverse bending and achieve thermal periscope. Besides, we further consider thermal cloak and concentrator that make their interior object from being thermally detected by restoring their exterior background temperature profile. Inspired by the scattering cancellation method6, we propose a regionalized scattering cancellation method to design thermal metadevices based on Eqs. (2) and (3). We divide the whole design area into several subareas, solve the thermal conduction equation with particular boundary constraints between the subareas, and obtain the relationship between the thermal conductivity tensors in each subarea. Taking the thermal connector as an example, as shown in Fig. 4a, we divide the design area into four areas by three interfaces, and the isotherms in each area are parallel. The shape of the thermal connector is determined by the geometry parameters β1 and β2, and the temperature gradient components of different areas are described as

$$\left\{ {\begin{array}{*{20}{c}} {\frac{{\partial T^{(1)}}}{{\partial x}} = \frac{{\partial T^{(2)}}}{{\partial x}}{{{\mathrm{ = }}}}\frac{{\partial T^{(3)}}}{{\partial x}} = \frac{{\partial T^{(4)}}}{{\partial x}} = - \phi } \\ {\frac{{\partial T^{(1)}}}{{\partial y}} = \frac{{\partial T^{(2)}}}{{\partial y}} = \frac{{\partial T^{(3)}}}{{\partial y}} = \frac{{\partial T^{(4)}}}{{\partial y}} = 0{{{\mathrm{ }}}}} \end{array}} \right.$$
(4)

where ϕ is a constant to be determined. By substituting Eq. (4) into Eq. (2), we obtain a general solution as follows (see Supplementary Notes 2 and 3 for the details):

$$\kappa _b = \kappa _{xx}^{(2)},\kappa _{yx}^{(2)} = \kappa _{yx}^{(3)},\kappa _b = \kappa _{xx}^{(3)},\cot \beta _1 = \frac{{\kappa _{xx}^{(2)}}}{{\kappa _{yx}^{(2)}}},\cot \beta _2 = \frac{{\kappa _{xx}^{(3)}}}{{\kappa _{yx}^{(3)}}}$$
(5)

where κb is the thermal conductivity of background material (Area 1 in Fig. 4a). When cotβ1 = 1 and cotβ2 = 1, Eq. (5) can be simplified to \(\kappa _b = \kappa _{xx}^{(2)} = \kappa _{yx}^{(2)} = \kappa _{xx}^{(3)} = \kappa _{yx}^{(3)}\). With consideration of the reciprocity and nonnegativity of κ, we obtain \({{{\mathbf{\kappa }}}}^{(1)} = {{{\mathbf{\kappa }}}}^{(4)} = \left[ {\begin{array}{*{20}{c}} {13} & 0 \\ 0 & {13} \end{array}} \right]\) and \({{{\mathbf{\kappa }}}}^{(2)} = {{{\mathbf{\kappa }}}}^{(3)} = \left[ {\begin{array}{*{20}{c}} {13} & {13} \\ {13} & {26} \end{array}} \right]\). Then, the theoretical simulation with the pure continuous materials for the thermal connector is performed. The result is displayed in Fig. 4b, where the thermal isotherms are parallel in the whole area, and the thermal connection functionality is theoretically achieved. Next, we need to design the corresponding TFCs with the desired anisotropic κ to fill each area. Topology optimization is employed by taking the required κ as the target, and the corresponding TFCs are obtained. The structurally simulated result of the thermal connector is shown in Fig. 4c, which looks similar to the theoretically simulated one (see Methods for details of theoretical simulations). Then, we fabricate the sample of topology-optimized thermal connector (seen in Fig. 4d) by 3D printing for experiment, and its experimentally steady-state temperature field distribution under the same boundary condition is shown in Fig. 4e (see Methods for details of experiments). For theoretical simulation, structural simulation and structural experiment, the corresponding temperature fields are similar, and thermal connection functionality is achieved. It is worth noting that the theoretical simulation result is so perfect and can be considered as the reference for structural simulation and structural experiment. We further plot the temperature profile on the red dashed line in Fig. 4c to evaluate the deviation between the simulated and experimental results. As shown in Fig. 4f, the results of theoretical simulation (black dash-dotted curve), structural simulation (red curve) and structural experiment (blue markers) agree well. The temperature values on the observation line are approximately equal, which is also a good verification of thermal connection functionality. Therefore, the TFCs with the desired ETCs from the full-parameter anisotropic space can robustly manipulate heat flow and achieve thermal connector numerically and experimentally.

Fig. 4: Design of thermal connector, and simulated and experimental verifications of its functionality.
figure 4

a Design of thermal connector based on regionalized scattering cancellation method. The orange area represents the background, and the blue and green areas denote thermal metamaterials. The white curves are isotherms. b Temperature field under theoretical simulation with target ETCs. c Temperature field under structural simulation of TFCs composed of two materials (copper and PDMS). d Samples of thermal metamaterials by 3D printing (without the PVC cover and PDMS fillings). f Experimentally measured temperature field distribution by IR camera. e Measured temperature values along the observation line marked in c under three cases. The blue markers denote the calculated temperature values from the data gathered by the IR camera.

Moreover, we also design other kinds of thermal metadevices, i.e., thermal reflector, concentrator and cloak, to further show the significance of TFCs in manipulating heat flow (see Supplementary Notes 26 for the details). For each thermal metadevice, we obtain the desired anisotropic κ by regionalized scattering cancellation method and provide a feasible solution in Supplementary Table 1. We achieve these thermal metadevices based on the computationally TFCs with required ETCs within Ω(Copper, PDMS). Figure 5a–f shows the theoretical and structural simulation results of thermal metadevices, respectively. Besides, we fabricate the samples (Fig. 5g–i) of the thermal metamaterials by 3D printing for experiments, and their experimentally steady-state temperature fields are shown in Fig. 5j–l. The temperature profiles along the observation red dashed lines are also plotted in Fig. 5m–o. From the above results, it is further illustrated that the proposed regionalized scattering cancellation method is effective and the full-parameter anisotropic space can be easily and robustly traversed by TFCs, which can be employed to manipulate heat flow freely and achieve thermal metadevices experimentally.

Fig. 5: Simulated and experimental results of thermal reflector, concentrator and cloak.
figure 5

ac Temperature fields under theoretical simulation with the target ETCs for the thermal reflector, concentrator and cloak. The white curves are isotherm. df Temperature fields under structural simulation of TFCs composed of two materials (copper and PDMS) for reflector, concentrator and cloak. Temperature observation along the red dashed lines is performed for comparison. gi Samples of thermal metamaterials by 3D printing in the experiments, not filled with PDMS and covered by PVC. jl Experimentally measured temperature fields by IR camera. mo Measured temperature value along the observation lines marked in d–f for thermal reflector, concentrator and cloak under three cases. The blue markers denote the calculated temperature value from the data gathered by the IR camera.

Discussion

We have demonstrated that the ETCs of rotated TFCs can traverse the full-parameter anisotropic space. Considering the manufacturing accuracy of 3D printing, the minimum feature size of metal structures should be in the order of magnitude of 10−4 m to ensure that the designed thermal metamaterials can be fabricated27, which can resort to increasing the filter radius28,29 during topology optimization of TFCs (see Supplementary Fig. 5a–c) and using the materials with lower thermal conductivity (see Supplementary Fig. 5b, d). Besides, although smaller TFCs can achieve better thermal performance, the size of TFCs cannot be infinitely small considering the manufacturing accuracy of 3D printing. Thus, the size of TFCs depends on the dimension of thermal metamaterials and the desired thermal performance of each TFC. For example, in designing a thermal connector, the length of |AB| in Fig. 4a is 50 mm. It is numerically and experimentally verified that the TFCs with 10 × 10 periodic arrangements can achieve the desired thermal performance well. Then, the size of each TFC is set as 5 mm × 5 mm.

In this work, thermal metamaterials with regular shape are designed in steady states. Complex and transient-state thermal metamaterials can also be designed based on the traversed full-parameter anisotropic space of thermal conductivity. Specifically, for designing thermal metamaterials in transient states, follow-up studies can be carried out by taking heat capacity and physical density of materials into account. For designing complex thermal metamaterials, such as fast heat sink setup30,31, enhanced energy harvesting medevices32, and multilayer printed circuit boards33, the combination of the proposed method in this work and some advanced multiscale topology optimization (MTO) strategies can be investigated, such as data-driven MTO34 and kriging-assisted MTO35.

In conclusion, we believe that discussions about manufacturing technology limitations and the design of complex thermal metamaterials will not alter the overall conclusions and insight gained from the presented work. This paper provides a universal and efficient inverse design of mixtures with two solid-state materials, whose ETCs can traverse their full-parameter anisotropic space of thermal conductivity, i.e., their entire 2D-case ETC set. Each ETC in the copper-PDMS full-parameter anisotropic space can be conveniently achieved by a topology-optimized TFC. The ETCs of three typical TFCs are validated numerically and experimentally. Then, we propose the regionalized scattering cancellation method and design the topology-optimized thermal metadevices based on the full-parameter anisotropic space. Their thermal functionalities are validated numerically and experimentally. As a powerful way to design mixtures with any ETC in the full-parameter anisotropic space, our work not only shows a promising strategy to traverse the full-parameter anisotropic space of thermal conductivity, but also exhibits a profound advance in conveniently designing thermal metamaterials to freely manipulate the heat flow, which may trigger the counterpart exploration methods in other physical fields, like acoustic36, mechanics37 and electromagnetics26.

Methods

Multiscale topology optimization

This study employs multiscale topology optimization to design thermal metamaterials. At the macroscale, the κ distribution of each subarea is computationally obtained via the regionalized scattering cancellation method. Then, at the microscale, the TFC with a target κ will be designed by topology optimization. The density-based topology optimization method SIMP18,19 (solid isotropic material with penalization) is employed to optimize TFCs, as it is effective and convenient for obtaining the optimized material distribution. In finite element analysis (FEA)38, the four-node rectangle elements are employed and the element thermal conductivity is defined by the interpolation \(\kappa (\rho _e){{{\mathrm{ = }}}}\kappa _{PDMS} + \rho _e^p(\kappa _{copper} - \kappa _{PDMS})\), where κcopper and κPDMS are the thermal conductivities of the conductor copper and insulator PDMS, respectively. ρe is a design variable, and p is the penalization factor that can help force the design variable as 0 and 1 in the final optimized structure. The element thermal conductivity matrix, global thermal conductivity matrix and global thermal load matrix are respectively expressed as \({{{\mathbf{k}}}}_e{{{\mathrm{ = }}}}\kappa (\rho _e){\int}_{V_e} {\left[ {\left( {\frac{{\partial {{{\mathbf{N}}}}}}{{\partial x}}} \right)^T\left( {\frac{{\partial {{{\mathbf{N}}}}}}{{\partial x}}} \right) + \left( {\frac{{\partial {{{\mathbf{N}}}}}}{{\partial y}}} \right)^T\left( {\frac{{\partial {{{\mathbf{N}}}}}}{{\partial y}}} \right)} \right]dV_e}\), \({{{\mathbf{K}}}}(\rho _e) = \mathop {\sum}\nolimits_{e = 1}^N {{{{\mathbf{k}}}}_e(\rho _e)}\) and \({{{\mathbf{Q}}}} = \mathop {\sum}\nolimits_{e = 1}^N {{{{\mathbf{Nq}}}}_e^0}\), where N is the shape function and Ve is the volume of a finite element. The thermal conduction governing partial differential equations is solved by using FEA, i.e., \({{{\mathbf{K}}}}(\rho _e){{{\mathbf{T}}}} = {{{\mathbf{Q}}}}\), where T is the global temperature matrix. To obtain a structure with the target ETC \(\kappa _{lm}^i(l,m = 1,2)\) effectively and conveniently, we set the minimization of the volume fraction of copper as the objective function C and the relative error between the desired and target ETC as a constraint function G. Hence, the topology optimization model can be depicted as

$$\begin{array}{l}\mathop {{\min }}\limits_{\rho _e} {{{\mathrm{ }}}}C{{{\mathrm{ = }}}}\frac{1}{{\left| V \right|}}\mathop {\sum}\limits_{e = 1}^N {\rho _e} \\ s.t.{{{\mathrm{ : }}}}{{{\mathbf{K}}}}(\rho _e){{{\mathbf{T}}}} = {{{\mathbf{Q}}}}\\ {{{\mathrm{ }}}}\kappa _{lm}^H{{{\mathrm{ = }}}}\frac{1}{{\left| V \right|}}\mathop {\sum}\limits_{e = 1}^N {(\Delta {{{\mathbf{T}}}}_e^{(l)})^T{{{\mathbf{k}}}}_e\Delta {{{\mathbf{T}}}}_e^{(m)}} ,(l,m = 1,2)\\ {{{\mathrm{ }}}}G = g\left( {(\kappa _{lm}^H - \kappa _{lm}^i)^2} \right) \le \varepsilon \\ {{{\mathrm{ }}}}0 \le \rho _e \le 1,{{{\mathrm{ }}}}e = 1,2...N\end{array}$$
(6)

where the numerical homogenization method39,40 is adopted to calculate the ETC tensor \(\kappa _{lm}^H\). V is the volume of the TFC. N is the total number of finite elements. \({\Delta} {\bf{T}}_{e}={\bf{T}}_{e}^{0}-{\bf{T}}_{e}\) is the temperature difference vector, where \({{{\mathbf{T}}}}_e^0\) is the nodal temperature vector under the uniform test heat flow \({{{\mathbf{q}}}}_e^0\) (e.g., {1, 0}T and {0, 1}T in the 2D case). Te is the induced nodal temperature field resulting from the FEA of a TFC, and ε is a small positive number. The heuristic formulas of the constraint function G are set as

$$\begin{array}{l}G = (\kappa _{11}^H - \kappa _{11}^i)^2/a + (\kappa _{22}^H - \kappa _{22}^i)^2/b + (\kappa _{12}^H - \kappa _{12}^i)^2 + (\kappa _{21}^H - \kappa _{21}^i)^2,\quad \quad \,\,\,{{{\mathrm{for }}}}\,\,\kappa _{12}^i = 0\\ G = (\kappa _{11}^H - \kappa _{11}^i)^2/a + (\kappa _{22}^H - \kappa _{22}^i)^2/b + (\kappa _{12}^H - \kappa _{12}^i)^2/c + (\kappa _{21}^H - \kappa _{21}^i)^2/d,\,{{{\mathrm{for }}}}\,\,\kappa _{12}^i \,\ne\, 0\end{array}$$
(7)

where a, b, c and d are dimensionless and take the absolute values of \(\kappa _{11}^i\), \(\kappa _{12}^i\), \(\kappa _{21}^i\) and \(\kappa _{22}^i\), respectively. A filter28 is used to avoid checkerboards and ensure manufacturability in each optimization iteration. The gradient-based method of moving asymptotes (MMA)41,42 is applied to update the design variables. The sensitivity analysis of C and G with respect to ρe is conducted by the adjoint method43. Finally, we will reprocess the structure obtained by topology optimization to get a smooth one. The designed TFCs in partitioned areas by the regionalized scattering cancellation method are then assembled to generate thermal metamaterials, which are further fabricated by 3D printing.

In topology optimization for obtaining the TFC with a specific ETC, the design domain is set as 5 mm × 5 mm and further meshed into 200 × 200 finite elements. We set that the penalization factor p = 5, the small positive number \(\varepsilon = 1e^{ - 4}\), and the filter radius of 3 is used. The thermal conductivities of copper and PDMS are set as 400 Wm−1K−1 and 0.16 Wm−1K−1, respectively. Schematically, Supplementary Fig. 3 shows the topology optimization process of some TFCs. The clear structures and small values of the constraint function G illustrate the effectiveness of the topology optimization model in Eq. (6).

Numerical simulation

For consistency with topology optimization, the thermal conductivities of copper and PDMS in the simulation are also set as 400 Wm−1K−1 and 0.16 Wm−1K−1, respectively. The numerical simulation is performed via the software COMSOL Multiphysics 5.5. The interface between softwares COMSOL Multiphysics 5.5 and Matlab R2017a is used to ensure the consistency of simulated and topological structure models. We conduct the linear heat flow to test the performance of TFCs and thermal metamaterials. Then, we calculate the steady-state temperature field distributions.

For numerical simulations of TFCs shown in the inset of Fig. 3m–o, each structure for simulation is constructed by 10 × 10 periodic arrangements of the cell obtained by topology optimization. Then, the total size of each structure for simulation is 50 mm × 50 mm, and the temperatures of the left and right boundaries of each structure are maintained at 423 K and 273 K, respectively. The upper and lower boundaries of each structure are thermally insulated. The temperature configuration is sustained during the entire simulation process, and the simulated temperature fields are shown in Fig. 3d–f. For comparison, we set a pure plate with the size of 50 mm × 50 mm, filled with an anisotropic material. The κ of the anisotropic material is set as the targeted one of the TFC. The same simulation boundary conditions are set. Numerically simulated results of the pure plate are shown in Fig. 3a–c.

For numerical simulations of thermal metamaterials, based on the specific design domains, thermal metamaterials are generated by shaping the 10 × 10 periodic arrangement of the TFCs into the desired shape and then assembling them. We set the background as the high thermal conductivity silica gel (13 Wm−1K−1). Besides, the left boundary of each thermal metadevice is set as the high-temperature wall (i.e., 363 K for connector and reflector, and 423 K for cloak and concentrator), while the right boundary (bottom boundary in thermal reflector) is set as the low-temperature wall (273 K). The other boundaries in the simulation are thermal insulations. Numerically simulated temperature fields of thermal metadevices are shown in Figs. 4c and 5d–f. For comparison, we set the corresponding area as pure anisotropic materials, and the κ of the anisotropic material is set as the target one (see Supplementary Table 1 for the κ of corresponding area). The same simulation boundary conditions are set. Numerically simulated results of thermal metadevices are shown in Figs. 4b and 5a–c, which are regarded as ideal thermal behaviors for the results in Figs. 4c and 5d–f.

Experiments

In experiments, the samples of TFCs and thermal metamaterials are fabricated by the following steps. We first shape the 10 × 10 periodic arrangement of the TFCs into the desired shape and assemble them. Then the STL models for 3D printing are generated after stretching 2D porous copper structures into 5 mm thickness. Finally, we fabricate the thermal metamaterials by 3D printing (seen in Figs. 3g–i, 4d and 5g–i), and fill PDMS into the copper structure in a mold. After deaerating, heating and trimming, the samples are ready. The whole thermal metadevices are covered by the 0.1 mm thick polyvinylchloride (PVC) and placed on the insulating foam.

For the experimental thermal system of the TFCs, the high temperature of 423 K is imposed on the left side of the fabricated sample through a temperature-control heating strip, while the cooling strip is imposed on the right side, and the temperature of the right side is fixed at 273 K.

For the experimental thermal system of thermal metamaterials, we use the high thermal conductivity silica gel (13 Wm−1K−1) with 5 mm thick as the background material, which can be cut to the desired shape. To decrease the effect of thermal resistance, we fill thermally conductive silicone grease (13 Wm−1K−1) into the seams between the samples and background material. The high temperature (i.e., 363 K for connector and reflector, and 423 K for cloak and concentrator) is imposed on the left side of the system through a temperature-control heating strip, while the cooling strip is imposed on the right side (bottom side in thermal reflector) of the system and the temperature is fixed at 273 K.

The temperature fields of all the thermal systems under the above thermal boundaries are captured by the infrared (IR) camera, and the temperature values of points on the observation lines are calculated by the collected data from the IR camera. The final temperature fields of TFCs and thermal metadevices are respectively shown in Figs. 3j–l, 4e and 5j–l (see Supplementary Videos 1 and 2 for the recorded and accelerated experimental dynamic temperature field during the whole process). The calculated temperature values of points on the observation lines for TFCs and thermal metadevices are respectively shown in Figs. 3m–o, 4f and 5m–o.