Control and design heat flux bending in thermal devices with transformation optics

Guoqiang Xu; Haochun Zhang; Yan Jin; Sen Li; Yao Li

doi:10.1364/OE.25.00A419

1. Introduction

Transformation optics is a tool used for bending the electromagnetic [1,2], optical [3], acoustic [4], and elastodynamic waves [5] in a distort coordinate system. In order to bend waves, the invariance of control equations should be maintained after the coordinate transformation. An electromagnetic cloak, which can map the material characteristic parameters onto the dimensional coordinates by associating electromagnetic transmission with coordinate transformation, is developed upon the independence of Maxwell equations [1,2]. Considering the extensional applications, which are corresponding to different coordinate transformations in electromagnetic field, many representative and peculiar electromagnetic devices have been investigated and established in recent years based on transformation optics, including concentrators [6], electromagnetic cavities [7], beam modulations [8], lenses [9], etc. In addition, a form invariance for the heat diffusion equation is also demonstrated in diffuse fields, both spherical and prolate spheroid cloaks are designed based on transformation optics [10]. Transformation thermotics is developed based on transformation optics by deducing the heat diffusion equation using a Jacobian matrix. For the control of heat, a 2D circular thermal cloak and a thermal concentrator are proposed [11]. Recently, some experiments have proved the characteristics of thermal cloaks and concentrators and indicated the possibility of heat control. Narayana et al. [12] constructed new artificial materials for thermal conduction and experimentally demonstrated the utility of guiding heat flux. Furthermore, they developed a 2D heat shield to reroute heat flux [13]. Schittny et al. [14] fabricated a 2D micro-structured thermal cloak, which is composed of alternate layers of copper and polydimethylsiloxanes (PDMS). They observed the temperature gradient distribution and indicated that the thermal protection performance was transient. Both 2D and 3D bilayer thermal cloaks made of bulk isotropic materials were demonstrated and their protection capabilities were verified experimentally [15]. A novel ultrathin 3D thermal cloak and a controllable thermal cloaking device with active thermoelectric components were developed based on the concept of magnetic cloaks without thermal metamaterials and transformation thermodynamics [16,17].

For the control of heat, some methods require further examination. Yang et al. investigated the equations of arbitrarily shaped objects, and the pushover theory was veriﬁed using numerical results [18]. Huang and associates [19,20] proposed some unconventional thermal cloak ideas in which the cloaked object could feel the flow of heat by controlling heat conduction. The infinite singularity of material parameters could be eliminated in a cylindrical cloak with the help of a new method in which the path trajectory of heat flux is governed using Hamiltonian [21]. Vemuri et al. indicated that the bending of heat flux in multilayered composite materials depends on both the composite rotation angle and the ratio of thermal conductivities [22]. Then, the behaviors of guiding the heat flux bending in composite materials were verified experimentally [23]. Fatih et al. [24] observed that the variability of the thickness as well as the thermal conductivity of interfaces in composites influence the deviation of heat flux bending significantly, and the variation in the thermal conductivity played a larger role than that of the thickness. Because the structure of a thermal cloak with composite materials can force the heat flux to bypass the protected object, the bending of heat flow should be controlled for designing the trajectory of heat flux on the thermal cloak. Present research studies have proposed the variation characters of heat flux bending on 2D thermal materials, which were fabricated by stacking thin sheets [22–24]. With the development and application of novel thermal devices, the design and control of heat flow in 3D cloaks [10,15–17,25,26] became attractive. For this purpose, unlike the present research studies [22–24], we explore the heat transport process in 3D space considering the changes in the geometrical curvature of cloaking systems. In this paper, we deduce the expressions for heat flux bending in 3D space based on transformation thermotics. Both 2D round and 3D spherical cloaks are developed to verify the expressions and simulate the performance of heat flux bending using Fluent.

2. Theoretical method and physical model

2.1 Derivation of heat flux bending in cloaking system

A heat diffusion equation describes the transfer of heat flux from warmer to colder parts. We consider 3D heat diffusion in a domain Ω without a heat source. The heat diffusion equation can be expressed as follows:

ρ (α) c (α) \frac{\partial T}{\partial t} = \nabla \cdot (κ (α) \nabla T) .

where T, ρ, c, and κ denote the transient temperature (K), density (kg∙m⁻³), specific heat capacity (J∙K⁻¹∙kg⁻¹), and thermal conductivity (W∙m⁻¹∙K⁻¹) corresponding to each point α = (x, y, z) in domain Ω of Cartesian coordinates, respectively. ∇ = e₁ ∂/∂x + e₂ ∂/∂y + e₃ ∂/∂z denotes the gradient. Heat transformation is considered under steady state ensuring the continuity of heat flux [4]. Hence, Eq. (1) possesses differential form invariance.

Considering the heat transfer in the series and parallel configurations of composite layers given in [22], the heat conductivity tensor at different positions can be expressed as follows:

κ = (\begin{matrix} κ_{x} & 0 & 0 \\ 0 & κ_{y} & 0 \\ 0 & 0 & κ_{z} \end{matrix}) = (\begin{matrix} \frac{(l_{n - 1} + l_{n}) κ_{n - 1} κ_{n}}{l_{n - 1} κ_{n} + l_{n} κ_{n - 1}} & 0 & 0 \\ 0 & \frac{κ_{n - 1} l_{n - 1} + κ_{n} l_{n}}{l {}_{n - 1}+ l_{n}} & 0 \\ 0 & 0 & \frac{κ_{n - 1} l_{n - 1} + κ_{n} l_{n}}{l {}_{n - 1}+ l_{n}} \end{matrix}), (n \geq 2) .

where l is the thickness of the composite layers and n denotes the number of layers.

To describe the process of heat transfer in a multilayered thermal cloaking scheme, coordinate transformations described by a Jacobian matrix should be performed according to the change in variable α = (x, y, z) in domain Ω to α’ = (x’, y’, z’) based on domain Ω’. The transformation process is illustrated in Fig. 1. Original space Ω is replaced by transformation domain Ω’ using a Jacobian matrix to expand a new invisible region in a transformational Cartesian coordinate system. R₁(α’) and R₂(α’) denote the inside and outside boundaries of the thermal cloak. In order to create a thermal invisible region with multiple layers between R₁(α’) and R₂(α’) as shown in Fig. 1(b), we need to compress the inner region of the radius less than R(α) in Fig. 1(a) along the radial direction considering the geometric transformation in [1,2,10,15–17,25–27]. Where $r = \sqrt{x^{2} + y^{2} + z^{2}}$ , $θ = a r c \cos \frac{z}{r}$ and $a r c \tan \frac{y}{x}$ denote the radial distance, azimuthal, and elevation in the original domain Ω, respectively. r’, θ’ and φ’ are the related direction components of transformation domain Ω’.

Fig. 1 The transport process. (a) Original domain Ω; (b) Transformation domain Ω’.

Download Full Size | PDF

r' = R_{1} (α) + \frac{(R_{2} (α) - R_{1} (α))}{R_{2} (α)} r, θ' = θ, φ' = φ .

The Jacobian matrix used for relating original and transformational Cartesian coordinate system can be expressed as Eq. (4) where ω = (R₂-R₁)/R₂:

\begin{matrix} J = \frac{\partial (x', y', z')}{\partial (x, y, z)} = \frac{\partial (x', y', z')}{\partial (r', θ', φ')} \frac{\partial (r', θ', φ')}{\partial (r, θ, φ)} \frac{\partial (r, θ, φ)}{\partial (x, y, z)} \\ = (\begin{matrix} \frac{\partial x'}{\partial r'} & \frac{\partial x'}{\partial θ'} & \frac{\partial x'}{\partial φ'} \\ \frac{\partial y'}{\partial r'} & \frac{\partial y'}{\partial θ'} & \frac{\partial y'}{\partial φ'} \\ \frac{\partial z'}{\partial r'} & \frac{\partial z'}{\partial θ'} & \frac{\partial x'}{\partial φ'} \end{matrix}) diag (ω, 1, 1) (\begin{matrix} \frac{\partial r}{\partial x} & \frac{\partial r}{\partial y} & \frac{\partial r}{\partial z} \\ \frac{\partial θ}{\partial x} & \frac{\partial θ}{\partial y} & \frac{\partial θ}{\partial z} \\ \frac{\partial φ}{\partial x} & \frac{\partial φ}{\partial y} & \frac{\partial φ}{\partial z} \end{matrix}) . \end{matrix}

Based on the Jacobian matrix, the heat conductivity tensor is associated with the new domain:

κ' = \frac{J κ J'}{\det (J)} .

where J’ is the transpose of Jacobian matrix J and det (J) denotes the determinant of the Jacobian matrix.

Substituting Eqs. (3) and (4) into Eq. (5), we can obtain the distribution of the heat conductivity tensor at different positions in transformation domain Ω’. The values of the matrix on the right side of Eq. (6) represent the varying heat conductivities in different directions:

κ' = (\begin{matrix} κ_{r' r'} & κ_{r' θ'} & κ_{r' φ'} \\ κ_{θ' r'} & κ_{θ' θ'} & κ_{θ' φ'} \\ κ_{φ' r'} & κ_{φ' θ'} & κ_{φ' φ'} \end{matrix}) .

Considering the form invariance of Maxwell equations after coordinate transformation [1], the geometric transformation in Eq. (3) and the Jacobian matrix in Eq. (4) can be also used to design invisible cloak schemes in electromagnetic field [1,2]. The thermal conductivity tensor in Eqs. (5) and (6) can be replaced by permittivity and permeability in the electromagnetic field. A number of novel phenomena and devices [1,2,6–9] corresponding to various coordinate transformations, can be representatively designed and fabricated using the theory of transformation optics in the electromagnetics area. In general, applications of transformation optics can be extended to multi-physical field coupling, which has a great potential value in establishing bifunctional or multifunctional devices [28].

According to the geometric transformation above, the diffusion equation with the heat flux along the x direction of the original domain can be written as [25–27]:

ρ' c' \frac{\partial T}{\partial t} = \frac{1}{{(r' - R_{1})}^{2}} (\begin{matrix} \frac{\partial}{\partial r'} (κ_{r'} {(r' - R_{1})}^{2} \frac{\partial T}{\partial r'}) + \frac{1}{\sin^{2} θ'} \frac{\partial}{\partial φ'} (κ_{φ'} \frac{\partial T}{\partial φ'}) \\ + \frac{1}{\sin θ'} \frac{\partial}{\partial θ'} (κ_{θ'} s i n θ' \frac{\partial T}{\partial θ'}) \end{matrix}) .

With the temperature gradient along the x direction of the original domain Ω as constant [22–24], the transformation Fourier law can be expressed to character the change in the heat flux density along r’, 𝜃’ and 𝜑’ by considering Eqs. (6) and (7) together. The superscripts denote the directions of heat flux and the subscripts are the directions of temperature gradients (∇T_r’, ∇T_𝜃’, ∇T_𝜑’) of the new domain Ω’:

\begin{matrix} q^{r'} = q_{{_{\nabla T}}_{_{r'}}}^{^{r'}} + q_{{_{\nabla T}}_{_{θ'}}}^{^{r'}} + q_{{_{\nabla T}}_{_{φ'}}}^{^{r'}} \\ = - (^{r' - R_{1}) 2} \sin θ' (κ_{x} \cos^{2} φ' \sin^{2} θ' + κ_{y} \sin^{2} φ' \sin^{2} θ' + κ_{z} \cos^{2} θ') \cdot \nabla T_{r'}, \\ - r' (r' - R_{1}) \sin^{2} θ' \cos θ' (κ_{x} \cos^{2} φ' + κ_{y} \sin^{2} φ' - κ_{z}) \cdot \nabla T_{θ'} \end{matrix}

\begin{matrix} q^{θ'} = q_{{_{\nabla T}}_{_{r'}}}^{^{θ'}} + q_{{_{\nabla T}}_{_{θ'}}}^{^{θ'}} + q_{{_{\nabla T}}_{_{φ'}}}^{^{θ'}} \\ = - r' (r' - R_{1}) \sin^{2} θ' \cos θ' (κ_{x} \cos^{2} φ' + κ_{y} \sin^{2} φ' - κ_{z}) \cdot \nabla T_{r'}, \\ - {(r')}^{2} \sin θ' (κ_{x} \cos^{2} φ' \cos^{2} θ' + κ_{y} \sin^{2} φ' \cos^{2} θ' + κ_{z} \sin^{2} θ') \cdot \nabla T_{θ'} \end{matrix}

\begin{matrix} q^{φ'} = q_{{_{\nabla T}}_{_{r'}}}^{^{φ'}} + q_{{_{\nabla T}}_{_{θ'}}}^{^{φ'}} + q_{{_{\nabla T}}_{_{φ'}}}^{^{φ'}} \\ = - r' (r' - R_{1}) \sin^{2} θ' \sin φ' \cos φ' (κ_{y} - κ_{x}) \cdot \nabla T_{r'} . \\ - {(r')}^{2} \cos θ' \sin φ' \cos φ' (κ_{y} - κ_{x}) \cdot \nabla T_{θ'} \end{matrix}

We observed from the above equations that there exist components in the θ’ and 𝜑’ directions of the heat flux density along the radial direction of temperature gradient ∇T_r’ and components in the r’ and 𝜑’ directions along the azimuthal direction of temperature gradient ∇T_𝜃’. Hence, we can derive the angle of heat flux bending in a 3D thermal cloak system space corresponding to the coordinate transformation in Eq. (3):

ϕ_{1} = \tan^{- 1} (\frac{q_{{_{\nabla T}}_{_{r'}}}^{^{θ'}}}{q_{{_{\nabla T}}_{_{r'}}}^{^{r'}}}) = \tan^{- 1} (\frac{r' \sin θ' \cos θ' (κ_{x} \cos^{2} φ' + κ_{y} \sin^{2} φ' - κ_{z})}{(r' - R_{1}) (κ_{x} \cos^{2} φ' \sin^{2} θ' + κ_{y} \sin^{2} φ' \sin^{2} θ' + κ_{z} \cos^{2} θ')}),

ϕ_{2} = \tan^{- 1} (\frac{q_{{_{\nabla T}}_{_{r'}}}^{^{φ'}}}{q_{{_{\nabla T}}_{_{r'}}}^{^{r'}}}) = \tan^{- 1} (\frac{r' \sin θ' \sin φ' \cos φ' (κ_{y} - κ_{x})}{(r' - R_{1}) (κ_{x} \cos^{2} φ' \sin^{2} θ' + κ_{y} \sin^{2} φ' \sin^{2} θ' + κ_{z} \cos^{2} θ')}),

ϕ_{3} = \tan^{- 1} (\frac{q_{{_{\nabla T}}_{_{r'}}}^{^{φ'}}}{q_{{_{\nabla T}}_{_{r'}}}^{^{θ'}}}) = \tan^{- 1} (\frac{\sin φ' \cos φ' (κ_{y} - κ_{x})}{\cos θ' (κ_{x} \cos^{2} φ' + κ_{y} \sin^{2} φ' - κ_{z})}),

ϕ_{4} = \tan^{- 1} (\frac{q_{{_{\nabla T}}_{_{θ'}}}^{^{r'}}}{q_{{_{\nabla T}}_{_{θ'}}}^{^{θ'}}}) = \tan^{- 1} (\frac{(r' - R_{1}) \sin θ' \cos θ' (κ_{x} \cos^{2} φ' + κ_{y} \sin^{2} φ' - κ_{z})}{r' (κ_{x} \cos^{2} φ' \cos^{2} θ' + κ_{y} \sin^{2} φ' \cos^{2} θ' + κ_{z} \sin^{2} θ')}),

ϕ_{5} = \tan^{- 1} (\frac{q_{{_{\nabla T}}_{_{θ'}}}^{^{φ'}}}{q_{{_{\nabla T}}_{_{θ'}}}^{^{θ'}}}) = \tan^{- 1} (\frac{\cot θ' \sin φ' \cos φ' (κ_{y} - κ_{x})}{κ_{x} \cos^{2} φ' \cos^{2} θ' + κ_{y} \sin^{2} φ' \cos^{2} θ' + κ_{z} \sin^{2} θ'}),

ϕ_{6} = \tan^{- 1} (\frac{q_{{_{\nabla T}}_{_{θ'}}}^{^{φ'}}}{q_{{_{\nabla T}}_{_{θ'}}}^{^{r'}}}) = \tan^{- 1} (\frac{r' \sin φ' \cos φ' (κ_{y} - κ_{x})}{(r' - R_{1}) \sin^{2} θ' (κ_{x} \cos^{2} φ' + κ_{y} \sin^{2} φ' - κ_{z})}) .

where ϕ₁–ϕ₃ are the components of heat flux bending in the radial direction of temperature gradient ∇T_r’ and ϕ₄–ϕ₆ signify the components in the ∇T_𝜃’ direction. Because Eq. (2) expresses the heat conductivity in the series and parallel configurations of the composite layers, we can change the form of Eqs. (11)-(16) by dividing both the numerators and denominators by κ_x; thus, we can obtain:

\frac{κ_{y}}{κ_{x}} = \frac{κ_{z}}{κ_{x}} = 1 + \frac{l_{n - 1} l_{n}}{{(l_{n - 1} + l_{n})}^{2}} \cdot \frac{{(κ_{n - 1} - κ_{n})}^{2}}{κ_{n - 1} κ_{n}}, (n \geq 2) .

In order to simplify the expressions, we introduce a new variable ε to indicate the relation between thermal conductivities and the thickness of layer in Eq. (13):

ε = \frac{l_{n - 1} l_{n}}{{(l_{n - 1} + l_{n})}^{2}} \cdot \frac{{(κ_{n - 1} - κ_{n})}^{2}}{κ_{n - 1} κ_{n}}, (n \geq 2) .

With the introduction of variable ε, the angle of heat flux bending can be written as following by substituting Eqs. (17) and (18) into Eqs. (11)-(16):

ϕ_{1} = \tan^{- 1} (\frac{- r' ε \sin θ' \cos θ' \cos^{2} φ'}{(r' - R_{1}) (1 + ε \sin^{2} φ' \sin^{2} θ' + ε \cos^{2} θ')}),

ϕ_{2} = \tan^{- 1} (\frac{r' ε \sin θ' \sin φ' \cos φ'}{(r' - R_{1}) (1 + ε \sin^{2} φ' \sin^{2} θ' + ε \cos^{2} θ')}),

ϕ_{3} = \tan^{- 1} (- \frac{\tan φ'}{\cos θ'}),

ϕ_{4} = \tan^{- 1} (\frac{- (r' - R_{1}) ε \sin θ' \cos θ' \cos^{2} φ'}{r' (1 + ε \sin^{2} φ' \cos^{2} θ' + ε \sin^{2} θ')}),

ϕ_{5} = \tan^{- 1} (\frac{ε \sin φ' \cos φ'}{\tan θ' (1 + ε \cos^{2} φ' \sin^{2} θ' + ε \sin^{2} θ')}), (θ ’ \neq \pm \frac{a π}{2}) .

ϕ_{6} = \tan^{- 1} (- \frac{r' \tan φ'}{(r' - R_{1}) \sin^{2} θ'}), (θ ’ \neq \pm \frac{a π}{2}) .

It is apparent that ϕ₁ and ϕ₂ are limited by the physical dimensions in the r’, θ’, and φ’ directions of the transformational coordinate, the variation in the thermal conductivities, and the thickness of the adjacent material layers. We note that ϕ₁, ϕ₂ and ϕ₄ could be zero simultaneously and ϕ₃ = ± φ’ if and only if $θ ’ \neq \pm \frac{a π}{2}$ (a is even). In other words, the heat flux is only along the ∇T_r’ without any components along ∇T_θ’ at this time. It can be observed from Eq. (21) that ϕ₃ depends only on the geometric changes of azimuths. When $θ ’ \neq \pm \frac{b π}{2}$ and $φ ’ \neq \pm \frac{c π}{2}$ (b is odd and c is integer) i.e. components r’ and θ’ are vertical with ∇T_θ’ = 0, ϕ₁ = 0, ϕ₂ ≠ 0 and ϕ₃, ϕ₄, ϕ₅, ϕ₆ don’t exist. Furthermore, if there’s no geometric changes in the radial direction and $θ ’ \neq \pm \frac{b π}{2}$ i.e. r’ = r in Eq. (3) in 2D domain, the expression of ϕ₂ would be the same as Eq. (9) in [22]. In summary, from the values of ϕ₁–ϕ₆, we can decide the deflection angle of heat flux in different temperature gradient directions, providing a new way to direct the distribution of heat flux density in channeling thermal materials to energy devices design, such as thermal cloaks, concentrators, and diodes.

2.2 Physical model

Inspired by the previous work of Vemuri [22–24] wherein rotating material layers were used, a 2D circle scheme based on coordinate transformation used in [12–15] is built to investigate the heat flux bending in the cloaking region. Meanwhile, a 3D spherical scheme [1,2,10,15–17,25–27] is established considering the space conversion in Eq. (3) to visually express the variations in heat flux bending along different directions by applying a constant temperature gradient along x. The transformed conductivity [25,26] in the cloaking region of 3D spherical scheme can be derived with the Jacobian matrix in Eqs. (4) and (5), where κ₀ denotes the isotropic heat conductivity of the surrounding of the cloak:

κ_{r} = \frac{κ_{0} R_{2}}{R_{2} - R_{1}} {(\frac{r' - R_{1}}{r'})}^{2}, κ_{θ} = κ_{φ} = \frac{κ_{0} R_{2}}{R_{2} - R_{1}} .

Figure 2 presents the 2D and 3D cloaking schemes. For the 2D cloaking scheme shown in Fig. 2(a), we built a 2D square plate of dimension 200 mm × 200 mm as the calculation domain. A copper protected object of radius (r) 25 mm was set at the center, and ten homogeneous thermal material layers of thickness (δ) 3 mm each were filled in the cloaking region. The thermal conductivity (W∙m⁻¹∙K⁻¹) of each layer was κ_2D,₁ = 0.15, κ_2D,₂ = 390.0, κ_2D,₃ = 2.63, κ_2D,₄ = 385.0, κ_2D,₅ = 9.85, κ_2D,₆ = 380, κ_2D,₇ = 18.5, κ_2D,₈ = 375, κ_2D,₉ = 28.5, and κ_2D,₁₀ = 370. For the 3D cloaking scheme shown in Fig. 2(b), we adopted the profile in our previous design [27] based on coordinate transformation in Eq. (3), eight homogeneous thermal material layers of thickness (δ) 3 mm each were filled in the cloaking region, and a copper ball of radius 25 mm was set in the protected region. The thermal conductivities (W∙m⁻¹∙K⁻¹) of the layers were selected as following in order to fulfill the criterion in Eq. (25): κ_3D_,1 = 0.15, κ_3D_,2 = 398, κ_3D_,3 = 2.87, κ_3D_,4 = 394.04, κ_3D_,5 = 11.55, κ_3D_,6 = 385.1, κ_3D_,7 = 18.83, and κ_3D_,8 = 378.73. A composite material of aluminum alloy (92Al–8Mg) and PDMS (polydimethylsiloxane) with an effective thermal conductivity [14] of 90 (W∙m⁻¹∙K⁻¹) was chosen as the background for both the 2D and 3D cloaking schemes. To validate our deductions, Fluent was employed to analyze the heat transfer process. The high-temperature boundary was set at a constant temperature of 393 K, and the low-temperature boundary on the opposite side was set at a constant temperature of 293K. The ambient temperature was adopted as 293 K in order to reduce the effects of convection caused by the temperature difference.

Fig. 2 Models of thermal cloaking schemes. (a) 2D scheme; (b) 3D scheme.

Download Full Size | PDF

3. Results and discussion

According to the geometric transformation and deductions, the heat flux line distribution of the reference cloaking schemes is obtained under certain conditions mentioned above. It can be observed from Figs. 3(a) and 3(b) that the heat flux lines in the background, corresponding to the temperature isotherms in the upper-right inset, are parallel. When the heat flux approached the cloaking structure, the direction of heat flux bent and was not orthogonal to the temperature isotherms in the background. Although there are differences in the thermal conductivity vectors of different layers and directions because of the application of Eq. (2), the bending direction of heat flux lines are diverse, i.e. the components of the heat flux direction on the interface between the adjacent layers are different because of different materials and direction vector components. As the temperature gradient along x is constant, only one component of the temperature gradient exists along r’ in the transformation domain. Hence, ϕ₂ is the single active bending angle. Different from a 2D cloaking scheme, two active temperature gradient components exist in the transformation domain along r’ and θ’ in a 3D cloaking scheme. Therefore, ϕ₁–ϕ₆ calculated using Eqs. (19)-(24) are active. Because of the symmetrical structure of the cloaking schemes, the heat flux bending angle can be visually displayed as in Figs. 3(a) and 3(b) in any direction. We observe that the heat flux path can be controlled by changing 𝜀 in Eq. (18) or the geometric directions toward the heat source. This presents potential applications in fabricating 2D or 3D thermal devices with two or more types of materials, such as cloaks, concentrators [29], rotations [30], and thermal lens [31].

Fig. 3 Heat flux line distribution for cloaking schemes and the upper-right illustrations show the temperature variation corresponding the heat flux line: (a) heat flux line for 2D cloaking scheme; (b) heat flux line for 3D cloaking scheme.

Download Full Size | PDF

3.1 Heat flux bending in 2D scheme

Because the control of heat flux bending can be modified by altering the thermal conductivity ratios and rotation angle as indicated in [22], several interfaces between the adjacent layers are selected to study the variation in bending considering the change in the r’ direction additionally. The bending angle calculated using Eq. (20) was −65.62° when φ’ was −50° on the interface between layers 1 and 2, which corresponds closely to −65.2° observed from the Fluent simulation at the same position. The above deductions can be rationalized with such values. The variations in the bending angle corresponding to the changes in the φ’ direction on the interfaces between layers 1 and 2, layers 5 and 6, layers 9 and 10, and layer 10 and the background with thermal conductivity ratios of 2600, 38.6, 13, and 0.24 are shown in Fig. 4, respectively.

Fig. 4 Variation of heat flux bending at interfaces (between layer 1 and 2, layer 5 and 6, layer 9 and 10, layer 10 and background, respectively) with different radios of thermal conductivities for 2D cloaking scheme.

Download Full Size | PDF

We infer that a larger value of ϕ₂ in each interface corresponds to a larger thermal conductivity ratio if the ratio is equal to or larger than one. In addition, the values of ϕ₂ are equal when the thermal conductivity ratios are the inverse of each other, which can be calculated by using Eq. (20), i.e. the value of ϕ₂ increases with decreasing thermal conductivity ratios if the ratio is less than 1. Because of the influence of the radial position in a cloaking scheme, the values of ϕ₂ in each interface are higher than that without considering the radial direction, which can be observed from Fig. 4. However, the radial position affects ϕ₂ more than the thermal conductivity ratios.

The changes in ϕ₂ in each interface are similar; first, ϕ₂ decreased with increasing φ’ and then increased rapidly when φ’ was in the range of −90° – 0°; if φ’ was in the range of 0° – 90°, ϕ₂ increased rapidly with increasing φ’ then decreased. It is clear from the variation in Fig. 4 that the values of ϕ₂ are antisymmetric, and the axis of symmetry is φ’ = 0. A higher bending angle was observed in the inner interface between layers 1 and 2 because of the highest thermal conductivity ratio and r’/ (r’-R₁). Modifying the heat flux bending angle by controlling both the thermal conductivity ratios and geometrical characteristics (r’ and φ’) enables us to change the heat flux direction in order to avoid/concentrate the heat at a specific position.

3.2 Heat flux bending in 3D scheme

Considering the practical application of thermal materials, a 3D cloaking scheme based on the coordinate transformation in Eq. (3) is developed to investigate the heat flux bending in random directions. Because a constant temperature gradient is adopted along x in the original coordinate, two active components of temperature gradient exist along r’ and θ’ (∇T_r’ and ∇T_θ’), and the component (∇T_𝜑’) along φ’ is zero according to the geometric transformation. We randomly chose θ’ = 52° and φ’ = −26° at the first interface between layers 1 and 2 to evaluate the accuracy of Eqs. (19)-(24). The six bending angles calculated by using Eqs. (19)-(24) are −64.18°, 58.58°, 38.39°, −3.24°, 5.95°, and −64.22°, which approached the values observed from Fluent −64.06°, 58.13°, 38.39°, −3.07°, 5.62°, and −64.18°, respectively. These data support the bending angle in 3D space.

The variations in the six bending angles corresponding to the changes in θ’ and φ’ when r is 28 mm are shown in Fig. 5. Additionally, the heat flux bending along ∇T_r’ is illustrated in Figs. 5(a)–5(c). Because of the effects of θ’ (θ’ ≠ ± 𝜋/2), the higher/lower values were obtained at higher/lower θ’ direction. The values of ϕ₁ are antisymmetric and centered at θ’ = 0°. Meanwhile, with the bending decreasing along θ’, the variation in ϕ₁ increases until θ’ = 0° as shown in Fig. 5(a). The changes in ϕ₂ as illustrated in Fig. 5(b) are similar to that of the 2D scheme in Fig. 4. However, ϕ₂ decreased with decreasing θ’, and the highest bending of ϕ₂ corresponding to φ’ were observed when θ’ was ± 𝜋/2. Furthermore, it was also antisymmetric and centered at θ’ = 0° because of the symmetry of the sine function in Eq. (20). It can be observed from Fig. 5(c) that ϕ₃ is an azimuth of θ’ and φ’ components along ∇T_r’. ϕ₃ is equal to φ’ if θ’ and φ’ axes are coincident, i.e., θ’ = 0°. In Eq. (21), because of the monotonic increase of denominator in the range of −90° < θ’ < 0°, ϕ₃ increases with increasing θ’. Besides, ϕ₃ was the same in the range of 0° < θ’ < 90° because of the symmetry of the cosine function, and it was antisymmetric and centered at 𝜑’ = 0° at a constant θ’ position. The bending variations along ∇T_θ’ is illustrated in Figs. 5(d)–5(f). Both ϕ₄ and ϕ₅ decreased with increasing bending along θ’. Meanwhile, ϕ₄ was 0° and ϕ₅ approached 0° as the tangent value of θ’ was infinite. However, the radial position affects ϕ₄ and not ϕ_5, which can be observed from Eqs. (22) and (23). Also, ϕ₄ was antisymmetric and ϕ₅ was positive symmetrical, and both were centered at θ’ = 0° at a constant φ’. ϕ₆ increased monotonically with φ’ at a constant θ’ as shown in Fig. 5(f), and the bending angle decreased with increased bending in θ’. Furthermore, ϕ₆ was the same in the ranges of −90° < θ’ < 0° and 0° < θ’ < 90°, and it was antisymmetric and centered at 𝜑’ = 0° at a constant θ’.

Fig. 5 Variation in heat flux bending on interface between layer 1 and 2 corresponding to the change of 𝜑 for 3D cloak scheme, (a) variation in ϕ₁; (b) variation in ϕ₂; (c) variation in ϕ₃; (d) variation in ϕ₄; (e) variation in ϕ₅; (f) variation in ϕ₆.

Download Full Size | PDF

The bending angles on the interface between layers 5 and 6 (r = 40 mm) are illustrated in Fig. 6. We infer that ϕ₃ was unchanged as shown in Fig. 6(c) because it is independent of the thermal conductivity and radial component; the trends of ϕ₅ and ϕ₆ presented in Figs. 6(e) and 6(f) were the same as that of Figs. 5(e) and 5(f). Nevertheless, ϕ₅, which is related to the thermal conductivity, reduced significantly, and ϕ₆, which is related to the radial direction, decreased slightly. Moreover, ϕ₁, ϕ₂, and ϕ₄ decreased with increasing r and decreasing thermal conductivity ratio. From the discussion above, we note that the path of heat flux can be designed directly by regulating the bending angles, which is seemed have potential applications in photothermal conversion.

Fig. 6 Variation in heat flux bending on interface between layer 5 and 6 corresponding to the change of 𝜑 for 3D cloak scheme, (a) variation in ϕ₁; (b) variation in ϕ₂; (c) variation in ϕ₃; (d) variation in ϕ₄; (e) variation in ϕ₅; (f) variation in ϕ₆.

Download Full Size | PDF

3.3 Control of heat flux bending in thermal materials

The heat flux bending corresponding to the change in the azimuth of θ’ and φ’ was discussed. In order to realize the control of heat flux path at random positions, we assumed two models. One has two alternated cuboid thermal materials of thickness (δ) 3 mm each with 𝜃’ = −45° without the radial component effect, i.e. R₁ = 0; another is a 3D multilayered cloaking scheme fabricated by using two materials with a constant thermal conductivity ratio of 1000 or 0.001.

The influence of different thermal conductivity ratios on ϕ₁, ϕ₂, and ϕ₄ are indicated in Figs. 7(a)–7(c). Owing to a rotation of −45° of the first model, the denominators in Eqs. (19) and (22) were equal, i.e. ϕ₁ and ϕ₃ were equal, as shown in Figs. 7(a) and 7(c). Generally, all the values of ϕ₁, ϕ₂, and ϕ₄ reduced with decreasing thermal conductivity ratios. In addition, the reduction in bending corresponding to 𝜑’ decreased with the decreasing magnitude of thermal conductivity ratios. The bending angles were zero if the thermal materials were the same. The effects of the second model on ϕ₁, ϕ₂, and ϕ₄ considering the radial component at 𝜃’ = −45° are shown in Figs. 7(d)–7(f). It can be observed that ϕ₁ and ϕ₂ increased with an increase in r’/ (r’-R₁), i.e. a reduction in the r’ component. However, the change in ϕ₆, which was deduced by using Eq. (22), is opposite to that of ϕ₁ and ϕ₂. As the range of thermal conductivity ratio, the increments in the bending angles under certain θ’ and φ’ components reduce until saturation approaching an infinitesimal value. Meanwhile, with reducing r component, the increments in ϕ₁ and ϕ₂ also decrease until saturation approaching an infinitesimal value. However, the change in the increments in ϕ₆ is opposite along its reduction. Combined with the discussion on ϕ₃, ϕ₅ and ϕ₆ above, it is clear that the heat flux bending in 3D thermal materials is achieved, which could be used to design novel thermal materials with a different transformation.

Fig. 7 Variation in heat flux bending upon the change in thermal material ratios without considering variation in radial positions, i.e. $\frac{r ’}{r ’ - R_{1}}$ = 1, R₁ = 0, when θ’ = −45° in 3D space, (a) variation in ϕ₁; (b) variation in ϕ₂; (c) variation in ϕ₄; Variation in heat flux bending upon the change in radial positions when θ’ = −45° and thermal conductivity ratio is 1000 or 0.001 in 3D space, (d) variation in ϕ₁; (e) variation in ϕ₂; (f) variation in ϕ₄.

Download Full Size | PDF

4. Conclusions

In this paper, we have provided the expressions for heat flux bending in both 2D and 3D thermal cloaking schemes based on transformation optics. Furthermore, we studied the influence of geometrical characteristics and anisotropic thermal conductivity on bending angles. From the results obtained, we indicate that there would be only an active component of ∇T_r’ in the 2D domain when applying a constant temperature gradient along the x direction of the original domain, i.e. θ’ = ± π/2, which is in agreement with previous research studies [15–17]. Moreover, we proposed the expressions for bending angles in the 3D transformational domain with a constant temperature gradient along x. Under certain conditions, two active components of ∇T_r’ and ∇T_θ’ with an inactive component of ∇T_𝜑’ = 0 exist with the heat flux density vectors along r’, θ’ and 𝜑’, providing six degrees of freedom of heat flux bending angles. Finally, we observed that the heat flux path can be controlled and designed in an anticipatory function directly by regulating azimuths, radial positions, and thermal conductivity ratios. The heat flux bending between the thermal materials is just like refraction leading to novel phenomena in thermal fields not confined to cloaking or concentrators. Benefitted by the method of channeling heat flux at any positions based on our study, many potential applications will be explored in fabricating new state-of-the-art thermal control devices.

Funding

National Natural Science Foundation of China (NSFC) (51436009, 51536001, 51572058).

References and links

1. J. B. Pendry, D. Schurig, and D. R. Smith, “Controlling electromagnetic fields,” Science 312(5781), 1780–1782 (2006). [CrossRef] [PubMed]

2. D. Schurig, J. B. Pendry, and D. R. Smith, “Calculation of material properties and ray tracing in transformation media,” Opt. Express 14(21), 9794–9804 (2006). [CrossRef] [PubMed]

3. J. Valentine, J. Li, T. Zentgraf, G. Bartal, and X. Zhang, “An optical cloak made of dielectrics,” Nat. Mater. 8(7), 568–571 (2009). [CrossRef] [PubMed]

4. Q. Li and J. S. Vipperman, “Two-dimensional acoustic cloaks of arbitrary shape with layered structure based on transformation acoustics,” Appl. Phys. Lett. 105(10), 101906 (2014). [CrossRef]

5. T. Bückmann, M. Thiel, M. Kadic, R. Schittny, and M. Wegener, “An elasto-mechanical unfeelability cloak made of pentamode metamaterials,” Nat. Commun. 5, 4130 (2014). [CrossRef] [PubMed]

6. J. Yang, M. Huang, C. Yang, Z. Xiao, and J. Peng, “Metamaterial electromagnetic concentrators with arbitrary geometries,” Opt. Express 17(22), 19656–19661 (2009). [CrossRef] [PubMed]

7. V. Ginis, P. Tassin, J. Danckaert, C. M. Soukoulis, and I. Veretennicoff, “Creating electromagnetic cavities using transformation optics,” New J. Phys. 14(3), 033007 (2012). [CrossRef]

8. P. H. Tichit, S. N. Burokur, and A. de Lustrac, “Reducing physical appearance of electromagnetic sources,” Opt. Express 21(4), 5053–5062 (2013). [CrossRef] [PubMed]

9. J. Yi, S. N. Burokur, G. P. Piau, and A. de Lustrac, “Coherent beam control with an all-dielectric transformation optics based lens,” Sci. Rep. 6(1), 18819 (2016). [CrossRef] [PubMed]

10. C. Z. Fan, Y. Gao, and J. P. Huang, “Shaped graded materials with an apparent negative thermal conductivity,” Appl. Phys. Lett. 92(25), 251907 (2008). [CrossRef]

11. S. Guenneau, C. Amra, and D. Veynante, “Transformation thermodynamics: cloaking and concentrating heat flux,” Opt. Express 20(7), 8207–8218 (2012). [CrossRef] [PubMed]

12. S. Narayana and Y. Sato, “Heat flux manipulation with engineered thermal materials,” Phys. Rev. Lett. 108(21), 214303 (2012). [CrossRef] [PubMed]

13. S. Narayana, S. Savo, and Y. Sato, “Transient heat flux shielding using thermal metamaterials,” Appl. Phys. Lett. 102(20), 201904 (2013). [CrossRef]

14. R. Schittny, M. Kadic, S. Guenneau, and M. Wegener, “Experiments on transformation thermodynamics: molding the flow of heat,” Phys. Rev. Lett. 110(19), 195901 (2013). [CrossRef] [PubMed]

15. T. Han, X. Bai, D. Gao, J. T. L. Thong, B. Li, and C. W. Qiu, “Experimental demonstration of a bilayer thermal cloak,” Phys. Rev. Lett. 112(5), 054302 (2014). [CrossRef] [PubMed]

16. H. Xu, X. Shi, F. Gao, H. Sun, and B. Zhang, “Ultrathin three-dimensional thermal cloak,” Phys. Rev. Lett. 112(5), 054301 (2014). [CrossRef] [PubMed]

17. M. N. Dang, H. Xu, Y. Zhang, and B. Zhang, “Active thermal cloak,” Appl. Phys. Lett. 107(12), 016623 (2015).

18. T. Yang, L. Huang, F. Chen, and W. Xu, “Heat flux and temperature field cloaks for arbitrarily shaped objects,” J. Phys. D Appl. Phys. 46(30), 305102 (2013). [CrossRef]

19. Y. Gao and J. P. Huang, “Unconventional thermal cloak hiding an object outside the cloak,” Europhys. Lett. 104(4), 468–477 (2013). [CrossRef]

20. X. Y. Shen and J. P. Huang, “Thermally hiding an object inside a cloak with feeling,” Int. J. Heat Mass Tran. 78(7), 1–6 (2014). [CrossRef]

21. L. Wu, “Cylindrical thermal cloak based on the path design of heat flux,” ASME J. Heat Trans. 137(2), 021301 (2015). [CrossRef]

22. K. P. Vemuri and P. R. Bandaru, “Geometrical considerations in the control and manipulation of conductive heat flux in multilayered thermal metamaterials,” Appl. Phys. Lett. 103(13), 133111 (2013). [CrossRef]

23. T. Yang, K. P. Vemuri, and P. R. Bandaru, “Experimental evidence for the bending of heat flux in a thermal metamaterial,” Appl. Phys. Lett. 105(8), 083908 (2014). [CrossRef]

24. F. M. Canbazoglu, K. P. Vemuri, and P. R. Bandaru, “Estimating interfacial thermal conductivity in metamaterials through heat flux mapping,” Appl. Phys. Lett. 106(14), 143904 (2015). [CrossRef]

25. T. Chen, C. N. Weng, and Y. L. Tsai, “Materials with constant anisotropic conductivity as a thermal cloak or concentrator,” J. Appl. Phys. 117(5), 054904 (2015). [CrossRef]

26. T. Han and C. W. Qiu, “Transformation Laplacian metamaterials: recent advances in manipulating thermal and dc fields,” J. Opt. 18(4), 044003 (2016). [CrossRef]

27. G. Q. Xu and H. C. Zhang, “A concept of heat dissipation coefficient for thermal cloak based on entropy generation approach,” AIP Adv. 6(9), 095107 (2016). [CrossRef]

28. C. Lan, B. Li, and J. Zhou, “Simultaneously concentrated electric and thermal fields using fan-shaped structure,” Opt. Express 23(19), 24475–24483 (2015). [CrossRef] [PubMed]

29. C. Lan, B. Li, and J. Zhou, “Simultaneously concentrated electric and thermal fields using fan-shaped structure,” Opt. Express 23(19), 24475–24483 (2015). [CrossRef] [PubMed]

30. S. Guenneau and C. Amra, “Anisotropic conductivity rotates heat fluxes in transient regimes,” Opt. Express 21(5), 6578–6583 (2013). [CrossRef] [PubMed]

31. Y. Liu, F. Sun, and S. He, “Novel thermal lens for remote heating/cooling designed with transformation optics,” Opt. Express 24(6), 5683–5692 (2016). [CrossRef] [PubMed]

Control and design heat flux bending in thermal devices with transformation optics

Abstract

1. Introduction

2. Theoretical method and physical model

2.1 Derivation of heat flux bending in cloaking system

2.2 Physical model

3. Results and discussion

3.1 Heat flux bending in 2D scheme

3.2 Heat flux bending in 3D scheme

3.3 Control of heat flux bending in thermal materials

4. Conclusions

Funding

References and links

Cited By

Figures (7)

Equations (25)

Optics Express