Design and Analysis of a Temperature-Sensitive Thermal Meta-Regulator Possessing Different Heat Distribution Modes

Abstract

The control and regulation of thermal fields is of great significance in solving various thermal management problems in human life. Benefitting from the emerging space transformation technique and thermal meta-material, thermal meta-structures with unique thermal control capabilities have been rapidly developed in recent years. However, the exploration of the functional diversity of thermal meta-materials and structures is still inadequate; most related works are still limited to the single-field control effect and lack sensitivity to external environment changes. For the designed functional structures, observation and analysis of energy fluctuations and irreversible heat loss during the regulation process of the diffusive thermal field are also scare. Therefore, in this current work, we design a thermal meta-regulator (based on the space transformation theory) that is capable of differently distributing thermal energy according to the heat input direction and switching field control pattern with the change of ambient temperature. In addition to the common indicator of temperature, we also introduce the local entropy production rate and the total entropy production in the thermo-dynamic category to carry out entropy analysis of the energy processes involved in the thermal meta-regulator, making a multi-angle evaluation of the structural performance. Furthermore, we use the statistical response surface method to explore the comprehensive/interaction effect of multiple influencing factors on the thermal meta-regulator; the derived regression equations can be used to accurately predict the structural effects under different design schemes and temperature conditions. Our work further enriches the diversity and flexibility of thermal field manipulation manners and the demonstrated functions are also expected to be realized in other physical fields.

1. Introduction

For decades, arbitrary control of energy in physical fields [1,2,3], such as the electromagnetic, optical, acoustics, and thermal fields, has historically been the focus of research. Since flexible and diverse manipulation of thermal energy is essential for various aspects of human life and macroscopic thermal control significantly influences the thermal rectification effect of many thermal management problems, achieving the on-demand domination of thermal energy has long been an urgent research direction. This topic has acquired breakthroughs in recent years, given the promotion and extension of space transformation technology in the thermal field [4,5,6,7,8,9], which originated from the concept of transformation optics [10]. Transformation optics is a mathematical tool for geometric interpretation of field governing equations. Benefitting from field governing equations that remain form-invariant during space transformation processes, the flux and energy distributions in one initial space can be mapped into another new space with only appropriate modifications to the constitutive parameters, thus building strong relationships between physical phenomena and material responses. The transformation optics theory first appeared in the optical field, then quickly spread to other wave fields such as the acoustics [11] and elastic fields [12], due to the formal similarity of the fields’ governing equations. Although the fundamental difference between the wave equation and the diffusion equation made the promotion of the transformation technique in the thermal field start relatively late, spatial mapping is still feasible to deal with heat transfer problems after the form invariance of the heat transfer equation has been rigorously confirmed [13]. For the exotic distributions of thermal properties resulting from transformation processes, which are difficult to achieve with a single natural material, thermal metamaterials as man-made composites with special structures can well match and exhibit these special characteristics [5,6]. On this basis, many novel thermal meta-structures have flourished, including thermal cloaks [4,14,15,16,17,18,19,20,21], thermal concentrators [20,21,22,23,24,25,26,27,28], thermal rotators [20,22,29], thermal camouflages [30], thermal illusions [31,32], and thermal messages [33,34], etc. Some studies have also built multi-functional meta-structures that can passively or actively adapt their functions to external environmental conditions or can be simultaneously applied to multiple physics fields [35,36,37,38,39,40,41]. However, the way that existing structures control and regulate thermal fields still needs to be improved; most of them only show the same response to different thermal environments, which may reduce the suitability of the thermal meta-structures to more diverse and complex application conditions. Additionally, the inherent diffusivity of the thermal field makes it worthwhile to observe and analyze energy fluctuations as well as irreversible losses during manipulation processes; this has not received sufficient attention so far. Although equilibrium and stability have been emphasized in classical science, irreversible processes are increasingly present in modern science. In this respect, entropy can visualize the thermo-dynamic properties of one system and is suitable for characterizing irreversible phenomena [28,42]; the degree of energy fluctuations and irreversible loss in the heat transfer process can be observed through the change in entropy production.

Therefore, in this current work, a temperature-dependent thermal meta-regulator with different heat manipulation modes is put forward by combing hybrid space transformation mappings and nonlinear responses. Under the comprehensive action of compression, regional rotation, and exponential function, the proposed thermal meta-regulator can show different thermal field control effects for opposite heat flow incoming directions and switch the working state according to the environmental temperature changes. In order to make a more complete observation and evaluation of the energy regulation process in the diffusive thermal field, the local entropy production rate and the total entropy production of the thermo-dynamics discipline are utilized, apart from the general temperature index. Aiming at critical design and environmental variables, after investigating the impact of each single factor on the structural function, the statistical response surface method is further applied to explore the comprehensive and interaction effect of multiple factors. The obtained regression equations are useful tools for predicting structural performance under specific design and ambient conditions. The findings can help to enrich the diversity and flexibility of thermal field regulation modes and the demonstrated functional structure can be potentially realized in various physical fields.

2. Model Establishment

To begin with, we first consider the following governing equation for the two-dimensional heat conduction process without inner sources:

(1)$\nabla ·\left(\kappa \nabla T\right)=I$

where κ denotes the thermal conductivity, T denotes the temperature distribution in two-dimensional Cartesian coordinate system (x, y). Different values of the state parameter I lead to different system states: I = 0 represents steady state, $I=\rho c\partial T/\partial t$ represents transient state, where ρ is mass density, c is specific heat capacity, and t is time. Based on the strictly proven form-invariance characteristic of heat transfer equations during the space transformation from virtual space (x′, y′) to real space (x, y) [13], the transformed form of the above formula can be written as:

(2)$\nabla \prime ·\left({\kappa }^{\prime }\nabla {\prime T}^{\prime }\right)=I\prime$

where the upper quotation mark denotes the parameters in virtual space. Still, I’ = 0 points to steady state and $I\prime =\rho \prime c\prime \partial T/\partial t$ points to transient state. In contrast, the space transformation for the transient heat transfer is more complex than the steady case. The application of transformation technique to transient heat transfer problems involves more issues, including the difference between the expected effect and the actual realization, the characteristics of the required thermal properties, and the applicability range of the corresponding structures; different specific transformation forms will have different trade-offs in these aspects [9,43]. Here, we use Jacobian transformation matrix J to establish the association between different systems and obtain new distributions of the key thermal properties after transformation:

(3)$J=\frac{\partial (x\prime ,y\prime )}{\partial (x,y)}$

(4)$\kappa \prime =\frac{J\kappa J^T}{\mathit{det}J}$

(5)$\rho \prime c\prime =\frac{\rho c}{detJ}$

Typically, thermal conductivity of the host is always set to a scalar value, thereby facilitating subsequent hybrid spatial mappings to handle the distribution of thermal conductivity in the thermal field without generating redundant terms [15,24,44]. We choose a constant κ₀ as the κ in original space, then separately replace ρ and c in original space with a constant ρ₀ and a constant c₀ on the same principle. J^T and detJ represent the transpose and determinant of the Jacobian transformation matrix, respectively.

Inspired by the hybrid space transformation theory [45], here we combine basic spatial mapping methods of compression and rotation to establish the relationship between coordinate systems and between thermal physical parameters in different spaces. In a two-dimensional Cartesian coordinate system, the compression/stretching of the x-axis is practically equivalent to the stretching/compression of the y-axis [45], so we compress in the x-axis direction alone and rotate the two main axes around the coordinate origin:

(6)$\begin{array}{cc}{\kappa }^{\prime }=\frac{J{\kappa }_0{J}^T}{\mathit{det}J}={\kappa }_0& ·\frac{R\left(\theta \right)·\frac{\left[\begin{array}{cc}n& 0\\ 0& 1\end{array}\right]·\left[\begin{array}{cc}n& 0\\ 0& 1\end{array}\right]}{n}·{R\left(\theta \right)}^T}{detR\left(\theta \right)}={\kappa }_0·\left[\begin{array}{cc}cos\theta & sin\theta \\ -sin\theta & cos\theta \end{array}\right]·\left[\begin{array}{cc}n& 0\\ 0& \frac{1}{n}\end{array}\right]·\left[\begin{array}{cc}cos\theta & -sin\theta \\ sin\theta & cos\theta \end{array}\right]\hfill \\ & ={\kappa }_0·\left[\begin{array}{cc}n{cos}^2\theta +\frac{1}{n}{sin}^2\theta & \frac{1}{n}sin\theta cos\theta -nsin\theta cos\theta \\ \frac{1}{n}sin\theta cos\theta -nsin\theta cos\theta & \frac{1}{n}{cos}^2\theta +n{sin}^2\theta \end{array}\right]\hfill \\ & ={\kappa }_0·\left[\begin{array}{cc}{\kappa \prime }_{xx}& {\kappa \prime }_{xy}\\ {\kappa \prime }_{yx}& {\kappa \prime }_{yy}\end{array}\right]\hfill \end{array}$

(7)${\rho }^{\prime }{c}^{\prime }=\frac{{\rho }_0{c}_0}{\mathit{det}J}=\frac{{\rho }_0{c}_0}{n}$

where 0 < n <1 is the compression ratio of the x-axis and θ represents the rotation angle of the main axis. By embedding this Cartesian coordinate system into a square functional domain and making the origin of the coordinate system coincide with the center of the functional domain, θ can be considered as a function of the y coordinate value [45]:

(8)$\theta =\pm \frac{1}{L}\xi y\pi$

where L is the side length of the square functional domain and ζ is a cycle parameter representing the periodic variation of the thermal conductivity distribution. The plus and minus signs, respectively, represent clockwise or counterclockwise of the main axis, which will result in different transformed thermal conductivity tensor and different heat distributions. Several classical works have rigorously demonstrated that nonlinear features including nonlinear thermal physical properties and nonlinear background media can be smoothly incorporated into the space transformation technique [35,36,37,38,46]. So, here, we integrate two forms of θ relying on temperature change:

(9)$\theta =\left\{\begin{array}{ccc}\frac{1}{L}\xi y\pi \left(\frac{2}{1+exp\left[\beta \left({T\mathrm{r}}_{}-{T\mathrm{c}}_{}\right)\right]}-1\right)& ,& T_{\mathrm{r}}\ne T_{\mathrm{c}}\\ \frac{1}{L}\xi y\pi & ,& T_{\mathrm{r}}=T_{\mathrm{c}}\end{array}\right.$

where Tr denotes the temperature at the center point of the square functional domain, which we take as a reference to compare with the preset critical temperature value Tc. As a result, the rotation direction of the main axis in the functional domain depends directly on the difference between Tr and Tc, which lays the foundation for the function switchability of the structure. The specific value of the constant coefficient β [36,37,46] is determined based on a concise mathematical analysis of Equations (8) and (9). In order to highlight the effect of β, we assume the remaining relevant parameters in both equations in a uniform and fundamental way: L = 1 m and ζ = 1. The temperature difference term Tr-Tc is also set to 1 K to amplify the sensitivity of the structure to ambient temperature variations. The structure that is valid for this small temperature difference is naturally effective for the much larger temperature difference hereafter in the settings. The variation range of the independent variable y is (0,1) m; the following Figure 1 shows the degree of convergence between the two equations for different values of β based on the objective θ. Taking the variation curve from Equation (8) as the reference, it is easy to observe that the larger the value of β, the more the variation curve from Equation (9) tends to be similar with Equation (8). When β increases to a certain extent, as shown by the dotted lines in the figure, the curves in this range have good approximations to the reference line, especially for small y values. Thus, combining the structural dimensions and environmental conditions during the subsequent simulations, we choose β = 3.5 K⁻¹ in the following work to ensure the function switching sensitivity of the thermal meta-regulator structure.

Furthermore, given the diffusivity of the thermal field and complementarity between heat transfer and thermo-dynamics disciplines, concepts of local entropy production rate and total entropy production are also applied to observe and analyze the function, as well as related energy process of the designed structure in a more comprehensive way [28,47]. The expression for the local entropy production rate ${{\dot{S}}^{″}}_g$ is as follows:

(10)$\begin{array}{c}{{\dot{S}}^{″}}_g=\frac{\partial }{\partial x^{\prime }}\left(\frac{q{\prime }_x}{T}\right)+\frac{\partial }{\partial y^{\prime }}\left(\frac{q{\prime }_y}{T}\right)=\frac{\partial }{\partial x^{\prime }}\left(\frac{-1}{T}\left[\begin{array}{cc}\kappa {\prime }_{xx}& \kappa {\prime }_{xy}\end{array}\right]\frac{\partial T}{\partial x^{\prime }}\right)+\frac{\partial }{\partial y^{\prime }}\left(\frac{-1}{T}\left[\begin{array}{cc}\kappa {\prime }_{yx}& \kappa {\prime }_{yy}\end{array}\right]\frac{\partial T}{\partial y^{\prime }}\right)\\ =\frac{-1}{T^2}\left[\begin{array}{cc}{\left(\frac{\partial T}{\partial x\prime }\right)}^2& {\left(\frac{\partial T}{\partial y\prime }\right)}^2\end{array}\right]\cdot \left[\begin{array}{cc}\kappa {\prime }_{xx}& \kappa {\prime }_{xy}\\ \kappa {\prime }_{yx}& \kappa {\prime }_{yy}\end{array}\right]=\left[\begin{array}{cc}{{\dot{S}}^{″}}_{gx}& {{\dot{S}}^{″}}_{gy}\end{array}\right]\end{array}$

where the matrix form resulting from the anisotropic transformed thermal conductivity does not exist in the original space where thermal conductivity is the constant κ₀. It can be seen from the above formula that the local entropy production rate ${{\dot{S}}^{″}}_g$ inside the distorted space varies with the temperature gradient and the thermal properties at each location, thus allowing visualization of the energy fluctuations and irreversible losses during heat transfer. Through the area integration of the local entropy production rate, the total entropy production $\dot{S_g}$ of in the target system can be written as:

(11)${\dot{S}}_g=\sum {\iint }_A{{\dot{S}}^{″}}_{g}dA$

where A denotes the area of integral region and can be defined and divided as needed.

To validate the function of the designed thermal meta-regulator structure, two sets of structural models are built in Figure 2 and a series of finite element simulations are carried out using the commercial software package COMSOL Multiphysics. For the thermal meta-regulator model in Figure 2a, the entire simulation domain is in a rectangular background that holds a constant thermal conductivity κ₀ = 19.6 W/(m∙K) and a constant product of mass density and specific heat capacity ρ₀c₀ = 13.44 × 10³ KJ/(m³∙K). Data of thermal properties here correspond to the nickel–steel material (50% Ni), which is one of the main materials for the construction of thermal meta-structures/devices due to its advantaged thermal stability and plasticity [28]. The background is split into two square areas with a side length of H = 0.2 m; these two square areas are connected to the same hot end, while the two cold ends are placed on the two other parallel sides of the square areas. The critical temperature value Tc = 100 K, the value of which is chosen by reference to the classical application of two-dimensional thermal meta-material in electronic devices [48,49] and the average maximum tolerable temperature of most electronic components. The temperatures at the hot and cold ends are chosen to be 90 K–40 K and 160 K–110 K for visual access to two thermal environments with the same temperature gradient distribution but below and above the critical temperature value, respectively. Other upper and lower boundaries of background are adiabatic. Then, in each divided square, a smaller square functional region with side length L = 0.14 m is dialed out from the same geometric center to which the transformed thermal conductivity κ’, as well as the transformed product ρ′c′, calculated by Equations (6)–(9), are assigned, where the compression ratio of x-axis n = 0.2 and cycle parameter ζ = 1. As a control, the bare plate model in Figure 2b has completely the same geometric dimensions, background materials, and boundary conditions as the thermal meta-regulator model, except for the absence of internal functional regions. Nevertheless, the overall area of the bare plate model is still divided into the layout of the thermal meta-regulator model; we define all local areas of both structural models for convenience and clarity for the follow-up research.

3. Results and Discussion

The steady-state temperature clouds of the thermal meta-regulator in different temperature environments are displayed in Figure 3. Within the temperature environment in Figure 3a, the temperature of the functional region center Tr is below the preset critical value Tc and the heat input from the right side in the functional region I_1,dis moves away from the upper and lower boundaries and gathers more towards the center. Hence, the heat accumulation near the center is significantly enhanced, forming a zone that is around 12.5% hotter than the background space in the same vertical direction, while, for heat transferred from left to right in the functional region I_2,dis, thermal energy is distributed in the opposite form. More heat is extruded to the top and bottom edges passing through the functional region, leading to a temperature difference of about 12.5%, but with a colder region at the center. When the temperatures at the hot end and cold ends rise so that Tr exceeds Tc, as shown in Figure 3b, the thermal conductivity distributions within two functional regions change so that the overall temperature distribution in the thermal meta-regulator is completely reversed. The heat comes from the right side in the functional region I_1,dis and is drawn more to the top and bottom rather than the middle path, while the heat entering the functional region I_2,dis from the left side is driven mainly in the center direction to notably increase the temperature there. Either case is clearly different from traditional thermal conductive media under the same ambient conditions. The unified and isotropic thermal properties of traditional thermal conductive media make it impossible to distinguish the direction of the incoming flow and can only assign the same temperature value to points on the same vertical line in the heat transfer path. Furthermore, based on the coordinate system in Figure 3a, we extract the temperature values on the control line y = 0 m in both the thermal meta-regulator model and the bare plate model for more quantitative comparisons, as shown in Figure 3c. Regardless of the temperature environment in which the bare plate is placed, the temperature on the control line rises then falls linearly from left to right, with temperature gradients in the same absolute value. However, in the thermal meta-regulator model, the temperature values on the control lines passing through the functional regions I_1,dis and I_2,dis are, respectively, higher and lower than those in the bare plate model when the applied low temperature field makes Tr lower than the preset Tc, while the applied high temperature field resulting in Tr higher than the preset Tc brings the exact opposite situation. The temperature difference between the two models is prominent inside the functional regions and increases then decreases along the positive x-axis, but significantly reduces in the background space. These all imply the peculiar field regulation ability of the designed thermal meta-regulator to different input heat directions in different thermal environments.

Additionally, the above settings are set to specifically show the special control effect on the thermal field directed by the concept of the designed thermal meta-regulator. This concept is universal; different math values of the threshold as well as ambient temperature may lead to different derivation results from the mathematical model, but the essential function and characteristics of the structure will not be affected. When further exploring the potential applications of the designed thermal meta-regulator, we can decide the temperature threshold value according to the actual environmental conditions or requirements at that time and still obtain a similar thermal response. For example, Figure 4 shows temperature distributions of the thermal meta-regulator model at Tc = 373 K and temperature settings of the hot end/cold ends are 363 K–313 K and 433 K–383 K. Obviously, the structure still exhibits essentially the same functions and characteristics as described above with the new threshold value and correspondingly different environmental conditions.

A series of transient simulations are also performed to observe and analyze the detailed process of function alternation of the thermal meta-regulator model stimulated by ambient temperature changes. Keeping other settings unchanged, we reset the initial temperature of the system to 40 K and reset the boundary conditions of the hot and cold ends to time-dependent piecewise functions. For the hot end there is:

(12)$T=\left\{\begin{array}{c}90\mathrm{K},0\mathrm{s}<t\le 15\mathrm{s}\\ 160\mathrm{K},15\mathrm{s}<t\le 40\mathrm{s}\end{array}\right.$

For the cold end there is:

(13)$T=\left\{\begin{array}{c}40\mathrm{K},0\mathrm{s}<t\le 15\mathrm{s}\\ 110\mathrm{K},15\mathrm{s}<t\le 40\mathrm{s}\end{array}\right.$

where each piecewise function is smoothed by the continuous function.

Figure 5 shows the distributions of the temperature as well as the components of the local entropy production rate at different times. The boundary condition settings for the current transient simulations ensure that the temperature of the functional region center Tr is below the preset critical temperature value Tc in the first 15 s. Therefore, at t = 13 s, the heat input from the right and left into the square functional regions is individually concentrated towards the middle of the region I_1,dis and the sides of the region I_2,dis, respectively, resulting in temperature rises at corresponding locations. Meanwhile, the local entropy production rate component ${{\dot{S}}^{″}}_{gx}$ is distributed in a parabola form, with the opening to the right in both functional regions I_1,dis and I_2,dis. Larger ${{\dot{S}}^{″}}_{gx}$ values exist at the heat transfer front, corresponding to the temperature profile. In background spaces II_1,dis and II_2,dis, some high ${{\dot{S}}^{″}}_{gx}$ values appear near the upper and lower right corners of two functional regions at this moment. These locations correspond to the demarcation between the special thermal conductivity distribution as well as peculiar thermal distribution in the functional regions and the ordinary thermal conductivity, as well as common thermal distribution in the background spaces. After more than 15 s of simulation time, at t = 17 s, the boundary temperatures changes for Tr above the preset Tc and the newly input heat to each square functional region begins to be regulated in the opposite way to the earlier distribution pattern; the heat incoming from the right (left) begins to be squeezed to the sides (center) of the region I_1,dis (I_2,dis). The new heat input in both functional regions also propagates in a parabolic form and produces the local entropy production rate component ${{\dot{S}}^{″}}_{gx}$, but the opening is to the left. At this point, different distribution patterns of ${{\dot{S}}^{″}}_{gx}$ coexist, and large ${{\dot{S}}^{″}}_{gx}$ values can be seen at the junction of these two distribution shapes. The background spaces around corners of the functional regions near the hot end hold significantly high ${{\dot{S}}^{″}}_{gx}$ due to the alternating thermal energy distribution and the violent field fluctuations inside the functional regions. The previous ${{\dot{S}}^{″}}_{gx}$ values near the upper right and lower right corners of the functional regions obviously decrease, as the energy distribution in these local areas tends to equalize before the arrival of new thermal energy. Finally, at t = 28 s, the heat transfer process under the new ambient condition is basically completed. A central low-temperature zone and a central high-temperature zone appear, respectively, inside the functional regions I_1,dis and I_2,dis, hinting at a complete functional transition of the thermal meta-regulator model. Correspondingly, the distribution of ${{\dot{S}}^{″}}_{gx}$ in the two functional regions is shaped as a parabola with the opening purely to the left. The high ${{\dot{S}}^{″}}_{gx}$ values in the background space disappear; the significant decrease in ${{\dot{S}}^{″}}_{gx}$ value, together with the basically uniform color allocation of the whole system, indicate that the system is very close to equilibrium state. Moreover, in Figure 5(c1–c3), the changing patterns of ${{\dot{S}}^{″}}_{gy}$ are very similar to those of ${{\dot{S}}^{″}}_{gx}$; some subtle discrepancies are caused by the differences in thermal conductivity components and temperature gradients in different directions.

It also should be noted that, based on the current model where the change of thermal conductivity in the functional regions with the external environment is directly controlled by the mathematical terms, the experimental composition used in future corresponding experimental studies for the practical realization of this nonlinear response will also have an impact on the functional transition between the different working states of the thermal meta-regulator. This is to be expected based on some previous experimental research, e.g., chemical variations from phase change materials (PCMs) and physical movements from shape memory materials (SMMs) or temperature-dependent mechanical modules that may cause additional energy loss to the field control structures and associated systems [37,38,46,50,51,52]. Therefore, as an example, the distributions of temperature and the local entropy production rate during the function alternation procedure of the thermal meta-regulator model containing latent heat are additionally studied. Under the same environment and boundary conditions, the change in thermal conductivity distribution of the functional regions in different environments is simulated as the effect of phase change and the latent heat ∆H = 243 kJ/kg of n-Octadecane, one of the PCMs commonly used in experimental studies to construct thermal meta-devices [50,51], is supplemented to the settings of thermal physical parameters. Observing the results in Figure 6, it can be found that, although the introduction of latent heat makes the thermal meta-regulator model take significantly longer to complete the functional adjustment with ambient changes, the distributions and variation patterns of the temperature and the local entropy production rate components are still similar to the discussion above. For the two thermal field regulation modes that basically formed before and after the function switching occurs, the temperature and the local entropy production rate components are essentially the same as in Figure 5; the main difference is during the function switching process. Since the change of thermal conductivity in the functional regions is at the expense of releasing latent heat, the thermal meta-regulator model only gradually demonstrates the new manipulation state for the heat inflow after a relatively long period of time when the environmental temperature changes. The energy fluctuations and irreversible losses within the thermal meta-regulator model and the whole domain also rise significantly, thus the overall functional switching process has more instability compared with the case in Figure 5. Nevertheless, the distributions of the temperature and the local entropy production rate components during this the process still maintain the pattern and characteristics of the corresponding state in Figure 5. These results further corroborate the effectiveness of the constructed thermal meta-regulator model.

According to theoretical derivations, the thermal conductivity distribution that directly determines the operating behavior of the thermal meta-regulator under each environmental condition is strongly correlated with various parameters, including the compression ratio of x-axis n, the temperature difference between the hot and cold ends △T, the cycle parameter ζ, and the area of the square functional regions that can be quantified by the side length L. Hence, it is necessary to investigate the influence of these key factors on the structural performance of the thermal meta-regulator in terms of field regulation effect and energy quality. From the results available above, we can find the high symmetry of thermal field pattern not only in the functional regions I_1,dis and I_2,dis, but also for the whole system in different temperature environments. So, we take the absolute value of the temperature difference between points A and B |△T_t| and the total entropy production ${\dot{S}}_g$ of the whole system as the target parameters. The former remains constant in both high- and low-temperature environments, while keeping the boundary temperature difference fixed and other parameters unchanged, so only variation laws in low-temperature environments where Tr is lower than the preset Tc are studied. The latter is a thermo-dynamic indicator of energy fluctuations and irreversible loss contained in the internal physical process inside a system. Taking the initial steady-state simulation scheme as the benchmark, we choose the above four strongly correlated parameters as independent variables and two target parameters as dependent variables in both the thermal meta-regulator and the bare plate models. Results are shown in Figure 7 and Figure 8.

In all cases in Figure 7, |△T_t| is always equal to 0 in the bare plate model, independent of the variation of the four independent variables, but possesses much higher value in the thermal meta-regulator model, indicating that all design schemes of the thermal meta-regulator model are effective. For the thermal meta-regulator model, the increase in the spatial distortion degree in the functional regions based on the decreasing compression ratio n will lead to a stronger regulation of the thermal field, so the smaller n leads to the larger |△T_t|. The change in the boundary temperature difference △T (including that in Figure 8 below) comes from fixing the hot-end temperature and then adjusting the cold-end temperature. The almost linear positive correlation between |△T_t| and △T suggests that the thermal meta-regulator model has greater application potential in environments with greater variability. The cycle parameter ζ denotes the repeatability of temperature profile shape in the functional regions. Referring to the initial scheme with ζ = 1, ζ values less than 1 imply a consistent but smaller amplitude of the thermal distribution trend, indicating the declined control degree of the thermal field in functional regions, so it is logical that |△T_t| falls with the decreasing ζ when ζ < 1. Instead, multiple identical heat distribution shapes are repeated in the functional regions when ζ > 1; the larger the ζ value, the higher the repetition frequency. Under this circumstance, locations of heat accumulation or dispersion will no longer be fixed and will gradually extend to the entire functional regions as the repetition frequency increases. As a result, |△T_t| drops significantly once ζ exceeds 1 and, although |△T_t| rises to some extent when ζ becomes an odd integer, it falls more violently as ζ increases to the next even integer, illustrating that the effect of the thermal meta-regulator model is relatively unsatisfactory in all scenarios where ζ > 1. As for the side length L, the larger L represents the larger area occupied by the functional regions, so the thermal energy will be further away from or closer to the center of the functional regions under the action of the thermal meta-regulator model. Hence, the larger L is, the stronger the thermal field control ability of the thermal meta-regulator.

Figure 8 shows the effect of each design variable on total entropy production ${\dot{S}}_g$ in two models. Here, we analyze the correlations between total entropy production ${\dot{S}}_g$ and different design variables under both high and low environmental conditions. Overall, the changing pattern of ${\dot{S}}_g$ with each variable is consistent for different temperature environments, but the ${\dot{S}}_g$ value in hot ambient is always lower than in cold ambient, which should be attributed to the lower local entropy production resulted from larger temperature values at spatial points in high-temperature environments for the same input energy. Then, for the bare plate model reflecting the traditional materials, the functional regions do not exist (let alone the spatial distortions and the peculiar distribution of thermal energy), so the parameters n, ζ, and L do not affect the thermo-dynamic characteristics of the bare plate system. However, for the thermal meta-regulator model, the growing n value causes the reduction in the spatial compression of the functional regions and subsequent energy fluctuations within the entire system, so ${\dot{S}}_g$ decreases exponentially with the increase in n. Furthermore, the increase in input thermal energy will intensify the rearrangement of the thermal field by the thermal meta-regulator model, so the increase in the boundary temperature difference △T will lead to an exponential growth in ${\dot{S}}_g$. It should be noted that, during the gradual increase in △T, the total entropy production ${\dot{S}}_g$ of the bare plate model also promotes exponentially rather than remaining constant, which can be inferred from the definition of the entropy concept as well as the increased energy and decreased temperature at most locations within the domain due to the increase in △T. Under the definition of repeatability of the temperature profile shape in the functional regions, ζ directly characterizes the spatial distortion degree in the functional regions; higher ζ values necessarily induce more energy disturbances and monotonic growth of ${\dot{S}}_g$ in the thermal meta-regulator model. In addition, the area expansion of the functional regions also leads to the broadening of the spatial deformation scope of the whole system, so that ${\dot{S}}_g$ increases with the increase in the side length L of functional regions, although this increasing trend deteriorates as L approaches the dimensionality of the background space.

Generally, a qualified thermal field control structure should produce as few energy fluctuations and losses as possible while high-efficiently controlling the thermal field, that is, achieving the high |△T_t| value while reducing the ${\dot{S}}_g$ value. In the previous analysis, we have discussed the single impacts of each key variable separately on |△T_t| and ${\dot{S}}_g$, which are not always consistent, so synergistic effects of variables on these two target parameters need to be further explored to directly predict and evaluate the performance of the thermal meta-regulator under the combined influences. Therefore, based on the temperature difference between points A and B |△T_t| in the thermal meta-regulator model and the total entropy production difference between the thermal meta-regulator model and the bare plate model |△${\dot{S}}_g$|, the response surface method (RSM) is selected to make comprehensive predictions of the structural performance and to assess the interaction between different key variables. Maintaining the original geometric dimensions and boundary conditions, we still choose the compression ratio of x-axis n, the temperature difference between the hot end and cold ends △T, the cycle parameter ζ, and the side length of the square functional regions L as variables and, respectively, denote them with symbols A, B, C, and D for convenience, as Table 1 shows.

Taking the critical temperature value Tc = 100 K, the variation range of n and △T is set to 0.1–0.5 and 10 K–50 K, respectively, where the acquisition of different temperature differences is the same as the settings in Figure 7 and Figure 8. Previous analysis shows that the ζ values greater than 1 are not conducive to the performance improvement of the thermal meta-regulator and the reduction of the system energy fluctuations, so here we only consider the case where ζ values are no more than 1 and only vary in the range of 0.2–1. L is adjusted within the range of 0.02 m–0.18 m. Afterwards, |△T_t|^0.2 and |△${\dot{S}}_g$|^0.2 are further processed as target parameters to refine the regression model and reduce unnecessary regression coefficients and a series of simulation schemes are designed for these two target parameters in different environments through arranging and combining different variables. Results of “analysis of variance” (ANOVA) for the obtained quadratic models are presented in Tables S1–S4 in the Supplementary Material, which strongly illustrate the rationality and reliability of each quadratic model. Therefore, on this basis, we can obtain the following regression equations for different target parameters in different environments, which are written in terms of coding factors. In the low-temperature environment with Tc = 100 K and the fixed hot-end temperature 90 K:

(14)$\begin{array}{c}{\left|∆T_{\mathrm{t}}\right|}^{0.2}=1.40-0.085A+0.20B+0.17C+0.20D-0.014AB-0.045AC-0.017AD+0.026BC+0.032BD\\ +0.061CD-0.021A^2-0.064B^2-0.12C^2-0.13D^2\end{array}$

(15)$\begin{array}{c}{\left|∆{\dot{S}}_g\right|}^{0.2}=1.89-0.26A+0.28B+0.20C+0.15D-0.064AB-0.093AC+0.16CD+0.054A^2\\ -0.074B^2-0.23D^2\end{array}$

In the high-temperature environment with Tc = 100 K and the fixed hot-end temperature 160 K:

(16)$\begin{array}{c}{\left|∆T_{\mathrm{t}}\right|}^{0.2}=1.40-0.085A+0.20B+0.17C+0.20D-0.014AB-0.045AC-0.017AD+0.026BC\\ +0.032BD+0.061CD-0.021A^2-0.064B^2-0.12C^2-0.13D^2\end{array}$

(17)$\begin{array}{c}{\left|∆{\dot{S}}_g\right|}^{0.2}=1.68-0.22A+0.25B+0.16C+0.14D-0.047AB-0.083AC+0.12CD+0.042A^2-0.071B^2\\ -0.21D^2\end{array}$

The same form of Equations (14) and (16) in turn confirms the high symmetry of the heat distribution within the thermal meta-regulator in different temperature environments.

Figure 9 shows the normal probability of internally studentized residuals for |△T_t|^0.2 and |△${\dot{S}}_g$|^0.2 in two temperature environments; the similarity between Figure 9a,c still comes from the consistency between Equations (14) and (16). In all cases, most internally studentized residuals are concentrated around the diagonal, which further implies that the regression models in Equations (14)–(17) are valid and informative enough to show the combined effects of selected variables on the regulation effect as well as the involved energy situation of the thermal meta-regulator model. It can be observed in the fitting equations that the influences of the single- and multi-variable items on the target parameters are not uniform, that is, the combination effects of different variable cannot be inferred from the influence law of each single variable on the target parameters. Therefore, it is necessary to visually explore the interactions between multiple variables. Based on the significant multi-variate items, we show a series of two-dimensional contour plots to illustrate the interaction characteristics of the different variables in each significant combination.

Benefiting from the same form of expression, the laws of action between the variables within the significant items discussed are applicable to the two environments chosen. Figure 10 shows the contour plots of interactions between each two variables within the significant items on |△T_t|^0.2. In Figure 10a–c, for the variable combinations containing the compression ratio of x-axis n, the smaller the n, the larger the values of other variables and the higher |△T_t|^0.2 value is obtained. The |△T_t|^0.2 value increases linearly with the decrease in n and the increase in other variables and the gradient of |△T_t|^0.2 along the variation direction of n is relatively small, which implies that the significance and the impact of compression ratio of x-axis on the field regulation effect of the thermal meta-regulator are weaker compared with other variables. In the paired combinations of the remaining three variables in Figure 10d–f, |△T_t|^0.2 has a monotonic positive correlation with each variable. Contour plots are almost symmetrically distributed along the diagonal from the upper right corner to the lower left and the gradient of |△T_t|^0.2 along the variation direction of each variable is essentially the same, which means that the three factors have almost the same degree of influence on the field regulation effect of the thermal meta-regulator. Thus, the high boundary temperature difference △T, large cycle parameter ζ, and the more proportional functional regions (represented by the larger side length L) together point to a stronger regulation function of the thermal meta-regulator.

Contour plots in Figure 11 show the interactions between variables within the significant items on |△${\dot{S}}_g$|^0.2, which are basically consistent under different environmental conditions. In the combinations of n-△T and n-ζ, |△${\dot{S}}_g$|^0.2 decreases monotonously with the increase in n and the decrease in △T or ζ. From the gradients of |△${\dot{S}}_g$|^0.2 along the variation directions of multiple variables, the impact of compression ratio of x-axis on the energy fluctuations and losses are similar with that of the boundary temperature difference but stronger than the cycle parameter. Furthermore, in the significant item consisting of ζ and L, |△${\dot{S}}_g$|^0.2 has a monotonic positive correlation with both variables, but the degree of influence of the circulation parameter on the energy fluctuations and losses is slightly more pronounced than that of the proportion of functional regions. Within the given ranges, either combine the side length of functional regions and the cycle parameter above and below moderate levels, respectively, or confine the side length of functional regions to vary only below moderate levels for facilitating control of the energy fluctuations and losses to a lesser extent.

Figure 12 illustrates the relationships between the predicted and the actual values of each target parameter in different environments; the similarity between Figure 12a,c still comes from the consistency between Equations (14) and (16). The high concentration of data near the diagonal means that the gap between the predicted and actual schemes is relatively small, so that the obtained regression equations can be used to predict the performance of the thermal meta-distributor when each variable takes on a certain value within a given range. In the following, Equations (14)–(17) are tested using random schemes to verify the accuracy and effectiveness for the performance prediction of the thermal meta-distributor based on the two target parameters. The variable values within each scheme and the comparisons between the predicted and actual results are presented in Table 2. The maximum error between the predicted and actual results is 2.88%, which strongly demonstrates the accurate predictive power of the resulting regression equations.

In the follow-up work, we will continue to conduct relevant experimental work as a validation and extension of current research. The nonlinear anisotropic thermal conductivity that exists in the theoretical model is key to realizing the special structural function. One solution based on our experience is to construct the composite configuration of the structure itself using material combinations with high thermal conductivity differences (e.g., copper and polydimethylsiloxane (PDMS)) and the equivalent medium theory (EMT) [36,38,43], then introduce temperature-dependent mechanical modules [52] to achieve the further temperature response. The alternative is to directly apply specific chemical material such as n-Octadecane in PCMs [51] (as the results in Figure 6 and related analysis) and SMMs [37,38,46], which has high requirements for the nonlinearity and deformability of material properties. Furthermore, based on the two-dimensional model used in theoretical analysis, a thickness of about 1 mm will be assigned to the out-of-plane direction to realize the practical processing of the solid thermal meta-regulator structure. It has been confirmed that this degree of out-of-plane thickness does not interfere with the performance of thermal meta-structures derived from planar designs [52], which also allows us to adopt the simplified two-dimensional design in previous sections. Naturally, multi-angle and multi-factor performance analysis and optimization of the actual experimental system will be carried out either.

Furthermore, in terms of the fundamental design strategies, although the theoretical method in this current work simultaneously involving hybrid spatial mapping and nonlinear responses has shown apparent uniqueness and high flexibility in thermal field manipulation, designing thermal meta-devices in a more concise and subtle way for breaking through the limitations of uncommon thermal properties thus maximizing structural feasibility is still a challenge to overcome in the future. Some scholars have introduced scattering cancellation to obtain homogeneous and isotropic distribution of thermal properties [8], but the advantages of this method are mainly confirmed in the thermal cloaking effect. For the multi-functional thermal meta-regulator that is expected currently, in addition to exploring the possibility of applying scattering cancellation strategy, introducing different mathematical techniques to improve the transformation thermotics method is an inevitable trend. All these efforts may contribute to promoting industrial applications of regulation devices for energy control and utilization.

4. Conclusions

In this work, we design a thermal meta-regulator structure based on space transformation technology and nonlinear response settings. The proposed structure can allocate the heat input from different directions in opposite ways and adjust regulation modes according to the ambient temperature changes, so that the center of the structure can be approximately 12.5% warmer or cooler than the same position of the bare plate. Then, the involved energy fluctuations and irreversible losses during regulation process are observed and the impacts of pivotal design and environmental parameters on structural regulation performance are investigated from both heat transfer and thermo-dynamics perspectives. Higher spatial deformation degree n = 1/11, larger environmental temperature difference △T = 50 K, moderate cycle parameter ζ = 1, and more functional area proportion L = 0.18 m not only enhance the structural allocation of thermal field but also aggravate the energy distribution and loss in the system. Finally, we introduce the statistical response surface analysis method to figure out the comprehensive and interaction effect of multiple influencing factors; the derived regression equations can be used to precisely predict structural performance under specific design and ambient conditions, whose error is no more than 2.88%. It is also worth noting that, benefitting from the proven high flexibility of the space transformation method [53,54], the established structure can be further extended to richer regulation forms and arbitrary shape profiles. Findings in this current work improve the flexibility and diversity of thermal field regulation modes; the proposed novel functional structure, as well as the design philosophy, is expected to be realized in other physical fields.

Author Contributions

Investigation, J.Z.; software, Y.C.; writing—original draft, Y.L.; writing—review and editing, H.Z. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Data is unavailable due to privacy.

Conflicts of Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Nomenclature

T	temperature
x	x-axis in Cartesian coordinates
y	y-axis in Cartesian coordinates
c	specific heat capacity
t	time
J	Jacobin matrix
n	compression ratio of x-axis
L	side length of functional region
$\dot{S_g^{″}}$	local entropy generation rate
${\dot{S}}_g$	total entropy generation
A	area of integral region
H	side length of computational domain
Greek Letters
∇	gradient operator
κ	thermal conductivity
ρ	density
θ	rotation angle of main axis
ζ	cycle parameter
β	scale coefficient
Superscripts
‘	parameter in transformation space
T	transposition of matrix
Subscripts
0	parameter value in original space
c	preset value
x	parameter component in the x-axis direction
y	parameter component in the y-axis direction

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. The degree of convergence between Equations (8) and (9) for different values of β based on the objective θ. Within the range marked by the dotted lines, β = 3.5 K−1 is chosen based on the structural dimensions and environmental conditions hereafter in the settings.

Figure 2. Model schematics of: (a) the designed thermal meta-regulator and (b) the bare plate. The rectangular background of the entire simulation domain is split into two square areas, which are connected to the same hot end, while the two cold ends are placed on the two other parallel sides. Other upper and lower boundaries of background are adiabatic. Except for the smaller square functional regions possessing transformed thermal properties present only in (a), the geometric dimensions, background materials, and boundary conditions of the two structural models are identical.

Figure 3. At the critical temperature value Tc = 100 K, steady-state temperature distributions of the thermal meta-regulator model in: (a) low temperature environment, the temperatures at high and low temperature ends are 90 K and 40 K, respectively; (b) high temperature environment, the temperatures at high and low temperature ends are 160 K and 110 K, respectively. Based on the coordinate system in (a), temperature values on the control line y = 0 m in two models under different environment conditions are shown in (c). A and B are selected research points for subsequent quantitative analysis, symbols I and II corresponding to the different regions divided in Figure 2.

Figure 4. At the critical temperature value Tc = 373 K, steady-state temperature distributions of the thermal meta-regulator model in: (a) low temperature environment, the temperatures at high and low temperature ends are 363 K and 313 K, respectively; (b) high temperature environment, the temperatures at high and low temperature ends are 433 K and 383 K, respectively. The thermal meta-regulator model still exhibits similar functions and characteristics.

Figure 5. Different distributions of the thermal meta-regulator model at different times: (a1–a3) temperature and (b1–c3) the local entropy production rate components. The initial temperature of the system is reset to 40 K and the boundary conditions of the hot and cold ends are reset to time-dependent piecewise functions. Other settings remain the same as in the steady-state simulations.

Figure 6. Different distributions of the thermal meta-regulator model containing latent heat at different times: (a1–a3) temperature and (b1–c3) the local entropy production rate components. Other settings remain the same as in Figure 5.

Figure 7. Effects of different design variables on temperature difference between points A and B |△Tt| in two models: (a) the compression ratio of x-axis n; (b) the temperature difference between the hot and cold ends △T; (c) the cycle parameter ζ; (d) the side length of the functional regions L. Based on the high symmetry of the thermal field pattern, only variation laws in low-temperature environments are studied, where the critical temperature value Tc = 100 K and the hot-end temperature is fixed at 90 K.

Figure 8. Effects of different design variables on total entropy production [Forumla omitted. See PDF.] in two models: (a) the compression ratio of x-axis n; (b) the temperature difference between the hot and cold ends △T; (c) the cycle parameter ζ; (d) the side length of the functional regions L. Variation laws in both high- and low-temperature environments are studied, where the critical temperature value Tc = 100 K and the hot-end temperature is fixed at 160 K and 90 K, respectively.

Figure 9. Normal probability of internally studentized residuals for |△Tt|0.2 and |△[Forumla omitted. See PDF.]|0.2 in different environments. Maps (a–d) are from regression equation Equation (14) and Equation (17), respectively. The concentrations of internally studentized residuals around the diagonals indicate that the regression models obtained are reasonably valid.

Figure 10. Contour plots of interactions between each two variables within the significant items on |△Tt|0.2: (a) interactions between n and △T; (b) interactions between n and ζ; (c) interactions between n and L; (d) interactions between △T and ζ; (e) interactions between △T and L; (f) interactions between ζ and L. The laws of action are applicable to both high and low temperature environments (Tr > Tc and Tr < Tc). Numbers in squares on the contour lines represent the promotion of the interaction between factors on the target value, the larger the number, the higher the target value.

Figure 11. Contour plots of interactions between variables within the significant items on |△[Forumla omitted. See PDF.]|0.2. In the low-temperature environment (Tr < Tc): (a1) interactions between n and △T; (a2) interactions between n and ζ; (a3) interactions between ζ and L. In the high-temperature environment (Tr > Tc): (b1) interactions between n and △T; (b2) interactions between n and ζ; (b3) interactions between ζ and L. Numbers in squares on the contour lines represent the promotion of the interaction between factors on the target value, the larger the number, the higher the target value.

Figure 12. Relationships between the predicted and the actual values of each target parameter in different environments. Maps (a–d) are from regression equation Equations (14)–(17), respectively. Colorful squares represent residuals and the high concentration of data near the diagonal means that the obtained regression equations have good predictability.

Table 1

Changing range and levels of key variables.

Symbol	Variable		Levels
		−1	0	1
A	compression ratio of x-axis n	0.1	0.3	0.5
B	boundary temperature difference △T (K)	10	30	50
C	cycle parameter ζ	0.2	0.6	1
D	side length of functional regions L (m)	0.02	0.1	0.18

Table 2

Comparison of the predicted and actual results based on different target parameters in different environments.

Tests		1	2	3	4	5
Design parameters	A(n)	0.20	0.20	0.40	0.40	0.35
	B(△T/K)	40	20	40	35	45
	C(ζ)	0.3	0.7	0.8	0.5	0.9
	D(L/m)	0.4	0.16	0.08	0.12	0.16
|△Tt|0.2(in both high- and low-temperature environments,Tr > Tc and Tr < Tc)	Predicted	1.098688	1.448663	1.420697	1.393189	1.692707
	Actual	1.131322	1.458638	1.443399	1.372943	1.659722
	Error	2.88%	0.68%	1.57%	−1.47%	−1.99%
|△${\dot{S}}_g$|0.2(in the low-temperature environment, Tr < Tc)	Predicted	1.839126	1.924902	1.880846	1.802404	2.186365
	Actual	1.836368	1.892262	1.913513	1.777581	2.204273
	Error	−0.15%	−1.72%	1.71%	−1.40%	0.81%
|△${\dot{S}}_g$|0.2(in the high-temperature environment, Tr > Tc)	Predicted	1.622161	1.714084	1.673369	1.623723	1.933115
	Actual	1.624022	1.680720	1.696791	1.600287	1.947644
	Error	0.11%	−1.99%	1.38%	−1.46%	0.75%

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/en16155807/s1, Table S1. Results of ANOVA for the quadratic model of |ΔT_t|^0.2 (in the low temperature environment, Tr < Tc). Table S2. Results of ANOVA for the quadratic model of |Δ${\dot{S}}_g$|^0.2 (in the low temperature environment, Tr < Tc). Table S3. Results of ANOVA for the quadratic model of |ΔT_t|^0.2 (in the high temperature environment, Tr > Tc). Table S4. Results of ANOVA for the quadratic model of |Δ${\dot{S}}_g$|^0.2 (in the high temperature environment, Tr > Tc).

References

1. Al-Thabaiti, S.; Mostafa, M.; Ahmed, A.; Salama, R. Synthesis of Copper/Chromium Metal Organic Frameworks- Derivatives as an Advanced Electrode Material for High-Performance Supercapacitors. Ceram. Int.; 2023; 49, pp. 5119-5129. [DOI: https://dx.doi.org/10.1016/j.ceramint.2022.10.029]

2. Alshorifi, F.T.; Tobbala, D.E.; El-Bahy, S.M.; Nassan, M.A.; Salama, R.S. The Role of Phosphotungstic Acid in Enhancing the Catalytic Performance of UiO-66 (Zr) and Its Applications as an Efficient Solid Acid Catalyst for Coumarins and Dihydropyrimidinones Synthesis. Catal. Commun.; 2022; 169, 106479. [DOI: https://dx.doi.org/10.1016/j.catcom.2022.106479]

3. Altass, H.M.; Ahmed, S.A.; Salama, R.S.; Moussa, Z.; Jassas, R.S.; Alsantali, R.I.; Al-Rooqi, M.M.; Ibrahim, A.A.; Khder, M.A.; Morad, M. et al. Low Temperature CO Oxidation Over Highly Active Gold Nanoparticles Supported on Reduced Graphene Oxide@Mg-BTC Nanocomposite. Catal. Lett.; 2023; 153, pp. 876-886. [DOI: https://dx.doi.org/10.1007/s10562-022-04026-y]

4. Fan, C.; Gao, Y.; Huang, J. Shaped Graded Materials with an Apparent Negative Thermal Conductivity. Appl. Phys. Lett.; 2008; 92, 251907. [DOI: https://dx.doi.org/10.1063/1.2951600]

5. Li, Y.; Li, W.; Han, T.; Zheng, X.; Li, J.; Li, B.; Fan, S.; Qiu, C. Transforming Heat Transfer with Thermal Metamaterials and Devices. Nat. Rev. Mater.; 2021; 6, pp. 488-507. [DOI: https://dx.doi.org/10.1038/s41578-021-00283-2]

6. Yang, S.; Wang, J.; Dai, G.; Yang, F.; Huang, J. Controlling Macroscopic Heat Transfer with Thermal Metamaterials: Theory, Experiment and Application. Phys. Rep.-Rev. Sect. Phys. Lett.; 2020; 908, pp. 1-65. [DOI: https://dx.doi.org/10.1016/j.physrep.2020.12.006]

7. Huang, J. Theoretical Thermotics: Transformation Thermotics and Extended Theories for Thermal Metamaterials; Springer: Berlin, Germany, 2020.

8. Xu, L.; Huang, J. Transformation Thermotics and Extended Theories: Inside and Outside Metamaterials; Springer: Berlin, Germany, 2023.

9. Sklan, S.R.; Bai, X.; Li, B.; Zhang, X. Detecting Thermal Cloaks via Transient Effects. Sci. Rep.; 2016; 6, 32915. [DOI: https://dx.doi.org/10.1038/srep32915]

10. Pendry, J.B.; Schurig, D.; Smith, D.R. Controlling Electromagnetic Fields. Science; 2006; 312, pp. 1780-1782. [DOI: https://dx.doi.org/10.1126/science.1125907]

11. Zhang, S.; Xia, C.; Fang, N. Broadband Acoustic Cloak for Ultrasound Waves. Phys. Rev. Lett.; 2011; 106, 024301. [DOI: https://dx.doi.org/10.1103/PhysRevLett.106.024301]

12. Bueckmann, T.; Thiel, M.; Kadic, M.; Schittny, R.; Wegener, M. An Elasto-Mechanical Unfeelability Cloak Made of Pentamode Metamaterials. Nat. Commun.; 2014; 5, 4130. [DOI: https://dx.doi.org/10.1038/ncomms5130] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24942191]

13. Guenneau, S.; Amra, C.; Veynante, D. Transformation Thermodynamics: Cloaking and Concentrating Heat Flux. Opt. Express; 2012; 20, pp. 8207-8218. [DOI: https://dx.doi.org/10.1364/OE.20.008207]

14. Li, J.Y.; Gao, Y.; Huang, J.P. A Bifunctional Cloak Using Transformation Media. J. Appl. Phys.; 2010; 108, 1780. [DOI: https://dx.doi.org/10.1063/1.3490226]

15. Han, T.; Yuan, T.; Li, B.; Qiu, C. Homogeneous Thermal Cloak with Constant Conductivity and Tunable Heat Localization. Sci. Rep.; 2013; 3, 1593. [DOI: https://dx.doi.org/10.1038/srep01593]

16. Xu, H.; Shi, X.; Gao, F.; Sun, H.; Zhang, B. Ultrathin Three-Dimensional Thermal Cloak. Phys. Rev. Lett.; 2014; 112, 054301. [DOI: https://dx.doi.org/10.1103/PhysRevLett.112.054301] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24580599]

17. Hu, R.; Xie, B.; Hu, J.; Chen, Q.; Luo, X. Carpet Thermal Cloak Realization Based on the Refraction Law of Heat Flux. EPL; 2015; 111, 54003. [DOI: https://dx.doi.org/10.1209/0295-5075/111/54003]

18. Nguyen, D.M.; Xu, H.; Zhang, Y.; Zhang, B. Active Thermal Cloak. Appl. Phys. Lett.; 2015; 107, 121901. [DOI: https://dx.doi.org/10.1063/1.4930989]

19. Guo, J.; Qu, Z.G. Thermal Cloak with Adaptive Heat Source to Proactively Manipulate Temperature Field in Heat Conduction Process. Int. J. Heat Mass Transf.; 2018; 127, pp. 1212-1222. [DOI: https://dx.doi.org/10.1016/j.ijheatmasstransfer.2018.07.035]

20. Sha, W.; Xiao, M.; Zhang, J.; Ren, X.; Zhu, Z.; Zhang, Y.; Xu, G.; Li, H.; Liu, X.; Chen, X. et al. Robustly Printable Freeform Thermal Metamaterials. Nat. Commun.; 2021; 12, 7228. [DOI: https://dx.doi.org/10.1038/s41467-021-27543-7]

21. Wang, Z.; Zhu, Z.; Liu, T.; Hu, R. Inverse Design of Thermal Metamaterials with Holey Engineering Strategy. J. Appl. Phys.; 2022; 132, 145102. [DOI: https://dx.doi.org/10.1063/5.0108743]

22. Narayana, S.; Sato, Y. Heat Flux Manipulation with Engineered Thermal Materials. Phys. Rev. Lett.; 2012; 108, 214303. [DOI: https://dx.doi.org/10.1103/PhysRevLett.108.214303]

23. Dede, E.M.; Nomura, T.; Schmalenberg, P.; Lee, J.S. Heat Flux Cloaking, Focusing, and Reversal in Ultra-Thin Composites Considering Conduction-Convection Effects. Appl. Phys. Lett.; 2013; 103, 063501. [DOI: https://dx.doi.org/10.1063/1.4816775]

24. Han, T.; Zhao, J.; Yuan, T.; Lei, D.Y.; Li, B.; Qiu, C. Theoretical Realization of an Ultra-Efficient Thermal-Energy Harvesting Cell Made of Natural Materials. Energy Environ. Sci.; 2013; 6, pp. 3537-3541. [DOI: https://dx.doi.org/10.1039/c3ee41512k]

25. Hu, R.; Wei, X.; Hu, J.; Luo, X. Local Heating Realization by Reverse Thermal Cloak. Sci. Rep.; 2014; 4, 3600. [DOI: https://dx.doi.org/10.1038/srep03600]

26. Chen, F.; Yuan Lei, D. Experimental Realization of Extreme Heat Flux Concentration with Easy-to-Make Thermal Metamaterials. Sci. Rep.; 2015; 5, 11552. [DOI: https://dx.doi.org/10.1038/srep11552] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26109080]

27. Han, T.; Bai, X.; Liu, D.; Gao, D.; Li, B.; Thong, J.T.; Qiu, C. Manipulating Steady Heat Conduction by Sensu-Shaped Thermal Metamaterials. Sci. Rep.; 2015; 5, 10242. [DOI: https://dx.doi.org/10.1038/srep10242] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25974383]

28. Xu, G.; Zhang, H.; Zou, Q.; Jin, Y.; Xie, M. Forecast of Thermal Harvesting Performance under Multi-Parameter Interaction with Response Surface Methodology. Int. J. Heat Mass Transf.; 2017; 115, pp. 682-693. [DOI: https://dx.doi.org/10.1016/j.ijheatmasstransfer.2017.07.085]

29. Guenneau, S.; Amra, C. Anisotropic Conductivity Rotates Heat Fluxes in Transient Regimes. Opt. Express; 2013; 21, pp. 6578-6583. [DOI: https://dx.doi.org/10.1364/OE.21.006578]

30. Hu, R.; Xi, W.; Liu, Y.; Tang, K.; Qiu, C.W. Thermal Camouflaging Metamaterials. Mater. Today; 2021; 45, pp. 120-141. [DOI: https://dx.doi.org/10.1016/j.mattod.2020.11.013]

31. Hu, R.; Zhou, S.; Li, Y.; Lei, D.; Luo, X.; Qiu, C. Illusion Thermotics. Adv. Mater.; 2018; 30, 1707237. [DOI: https://dx.doi.org/10.1002/adma.201707237]

32. Zhou, S.; Hu, R.; Luo, X. Thermal Illusion with Twinborn-like Heat Signatures. Int. J. Heat Mass Transf.; 2018; 127, pp. 607-613. [DOI: https://dx.doi.org/10.1016/j.ijheatmasstransfer.2018.07.034]

33. Hu, R.; Huang, S.; Wang, M.; Luo, X.; Shiomi, J.; Qiu, C. Encrypted Thermal Printing with Regionalization Transformation. Adv. Mater.; 2019; 31, 1807849. [DOI: https://dx.doi.org/10.1002/adma.201807849]

34. Hu, R.; Huang, S.; Wang, M.; Zhou, L.; Peng, X.; Luo, X. Binary Thermal Encoding by Energy Shielding and Harvesting Units. Phys. Rev. Appl.; 2018; 10, 054032. [DOI: https://dx.doi.org/10.1103/PhysRevApplied.10.054032]

35. Sklan, S.R.; Li, B. A Unified Approach to Nonlinear Transformation Materials. Sci. Rep.; 2018; 8, 4436. [DOI: https://dx.doi.org/10.1038/s41598-018-22215-x] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29535315]

36. Li, Y.; Shen, X.; Huang, J.; Ni, Y. Temperature-Dependent Transformation Thermotics for Unsteady States: Switchable Concentrator for Transient Heat Flow. Phys. Lett. A; 2016; 380, pp. 1641-1647. [DOI: https://dx.doi.org/10.1016/j.physleta.2016.02.040]

37. Shen, X.; Li, Y.; Jiang, C.; Ni, Y.; Huang, J. Thermal Cloak-Concentrator. Appl. Phys. Lett.; 2016; 109, 031907. [DOI: https://dx.doi.org/10.1063/1.4959251]

38. Shen, X.; Li, Y.; Jiang, C.; Huang, J. Temperature Trapping: Energy-Free Maintenance of Constant Temperatures as Ambient Temperature Gradients Change. Phys. Rev. Lett.; 2016; 117, 055501. [DOI: https://dx.doi.org/10.1103/PhysRevLett.117.055501] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27517778]

39. Ma, Y.; Liu, Y.; Raza, M.; Wang, Y.; He, S. Experimental Demonstration of a Multiphysics Cloak: Manipulating Heat Flux and Electric Current Simultaneously. Phys. Rev. Lett.; 2014; 113, 205501. [DOI: https://dx.doi.org/10.1103/PhysRevLett.113.205501]

40. Moccia, M.; Castaldi, G.; Savo, S.; Sato, Y.; Galdi, V. Independent Manipulation of Heat and Electrical Current via Bifunctional Metamaterials. Phys. Rev. X; 2014; 4, 021025. [DOI: https://dx.doi.org/10.1103/PhysRevX.4.021025]

41. Lan, C.; Li, B.; Zhou, J. Simultaneously Concentrated Electric and Thermal Fields Using Fan-Shaped Structure. Opt. Express; 2015; 23, pp. 24475-24483. [DOI: https://dx.doi.org/10.1364/OE.23.024475]

42. Torabi, M.; Aziz, A. Entropy Generation in a Hollow Cylinder with Temperature Dependent Thermal Conductivity and Internal Heat Generation with Convective–Radiative Surface Cooling. Int. Commun. Heat Mass Transf.; 2012; 39, pp. 1487-1495. [DOI: https://dx.doi.org/10.1016/j.icheatmasstransfer.2012.10.009]

43. Schittny, R.; Kadic, M.; Guenneau, S.; Wegener, M. Experiments on Transformation Thermodynamics: Molding the Flow of Heat. Phys. Rev. Lett.; 2013; 110, 195901. [DOI: https://dx.doi.org/10.1103/PhysRevLett.110.195901] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23705719]

44. Chen, H.; Chan, C.T. Transformation Media That Rotate Electromagnetic Fields. Appl. Phys. Lett.; 2007; 90, 241105. [DOI: https://dx.doi.org/10.1063/1.2748302]

45. Shen, X.; Jiang, C.; Li, Y.; Huang, J. Thermal Metamaterial for Convergent Transfer of Conductive Heat with High Efficiency. Appl. Phys. Lett.; 2016; 109, 201906. [DOI: https://dx.doi.org/10.1063/1.4967986]

46. Li, Y.; Shen, X.; Wu, Z.; Huang, J.; Chen, Y.; Ni, Y.; Huang, J. Temperature-Dependent Transformation Thermotics: From Switchable Thermal Cloaks to Macroscopic Thermal Diodes. Phys. Rev. Lett.; 2015; 115, 195503. [DOI: https://dx.doi.org/10.1103/PhysRevLett.115.195503]

47. Xu, G.; Zhang, H.; Zou, Q.; Jin, Y. Predicting and Analyzing Interaction of the Thermal Cloaking Performance through Response Surface Method. Int. J. Heat Mass Transf.; 2017; 109, pp. 746-754. [DOI: https://dx.doi.org/10.1016/j.ijheatmasstransfer.2017.02.056]

48. Dede, E.M.; Schmalenberg, P.; Nomura, T.; Ishigaki, M. Design of Anisotropic Thermal Conductivity in Multilayer Printed Circuit Boards. IEEE Trans. Compon. Packag. Manuf. Technol.; 2015; 5, pp. 1763-1774. [DOI: https://dx.doi.org/10.1109/TCPMT.2015.2473103]

49. Dede, E.M.; Schmalenberg, P.; Wang, C.-M.; Zhou, F.; Nomura, T. Collection of Low-Grade Waste Heat for Enhanced Energy Harvesting. AIP Adv.; 2016; 6, 4950861. [DOI: https://dx.doi.org/10.1063/1.4950861]

50. Velez, C.; Khayet, M.; Dezarate, J.M.O. Temperature-Dependent Thermal Properties of Solid/Liquid Phase Change Even-Numbered n-Alkanes: N-Hexadecane, n-Octadecane and n-Eicosane. Appl. Energy; 2015; 143, pp. 383-394. [DOI: https://dx.doi.org/10.1016/j.apenergy.2015.01.054]

51. Kang, S.; Cha, J.; Seo, K.; Kim, S.; Cha, Y.; Lee, H.; Park, J.; Choi, W. Temperature-Responsive Thermal Metamaterials Enabled by Modular Design of Thermally Tunable Unit Cells. Int. J. Heat Mass Transf.; 2019; 130, pp. 469-482. [DOI: https://dx.doi.org/10.1016/j.ijheatmasstransfer.2018.10.127]

52. Li, Y.; Zhang, H.; Chen, Y.; Zhang, D.; Huang, Z.; Wang, H. Realization and Analysis of an Intelligent Flux Transfer Regulator by Allocating Thermal and DC Electric Fields. Int. J. Heat Mass Transf.; 2021; 179, 121677. [DOI: https://dx.doi.org/10.1016/j.ijheatmasstransfer.2021.121677]

53. Xu, G.; Zhou, X. Manipulating Cell: Flexibly Manipulating Thermal and DC Fields in Arbitrary Domain. Opt. Express; 2019; 27, pp. 30819-30829. [DOI: https://dx.doi.org/10.1364/OE.27.030819] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31684325]

54. Zhang, X.; He, X.; Wu, L. Three-Dimensional Thermal Illusion Devices with Arbitrary Shape. Chin. Phys. B; 2019; 28, 064403. [DOI: https://dx.doi.org/10.1088/1674-1056/28/6/064403]