Transformation thermodynamics: cloaking and concentrating heat flux

Sebastien Guenneau; Claude Amra; Denis Veynante

doi:10.1364/OE.20.008207

1. Introduction

There is currently a renewed interest in gradient index materials (whose paradigm is the Luneburg lens [1]) in optics and acoustics owing to their links with transformational optics and acoustics. In 2006, two research groups (those of Pendry, Schurig and Smith [2] and Leonhardt [3]) independently proposed a systematic way to control wave trajectories in curvilinear coordinate systems. In order to attract attention of mass media, they designed a cloak that renders any object inside it invisible to electromagnetic radiation. The former team theoretized that a coating consisting of a meta-material whose physical properties are deduced from a coordinate transformation in the Maxwell system displays anisotropy and heterogeneity of permittivity and permeability working as a deformation of the optical space around the object. The physicists consider the blowup of a point, thereby tearing apart the metric space. Though this may seem haphazardous, this can be legitimated by making use of advanced mathematical treatments [4, 5]. The experimental validation [6] of their theoretical considerations was given in 2006 for a copper cylinder invisible to an incident plane wave at 8.5 GHz. In 2009, the group of de Lustrac proposed a non-magnetic cloak, with anisotropic and spatially varying permittivity, which was experimentally demonstrated at microwave frequencies [7]. However, Leonhardt’s team used mathematical tools of complex analysis in order to conformally map the Helmholtz equation, thereby ensuring that the (spatially varying) refractive index be isotropic. These two approaches to cloaking (conformal and non-conformal) markedly enhance our capabilities to manipulate light, even in the intense near field limit [8]. Other research groups looked at how governing equations behave under geometric transforms in other areas of physics, such as the conductivity equation [9], the Schrödinger equation for matter waves [10,11,12,13], the acoustic wave equation [14, 15, 16] and the Navier equation for elastodynamics [17, 18]. One should also note that it is possible to drastically reduce the scattering cross section of an object thanks to a homogeneous coating with close to zero [19], or negative [20] refractive index.

In the present paper, we extend the design of transformation based metamaterials to the area of thermodynamics. We first show that the diffusion equation retains its form under geometric transform and we derive the expression of its heterogeneous anisotropic diffusivity in a general coordinate system. We then apply specific transforms in order to design a cloak and a concentrator for heat flux. We use full wave computations to back up our claims of an unprecedented control of heat flux through a change of metric in the thermic space. We finally propose a multilayered device working as a thermic protection in the homogenization regime.

2. Transformed heat equation

We consider the two-dimensional diffusion equation in a bounded cylindrical domain Ω with a source p

ρ (x) c (x) \frac{\partial u}{\partial t} = \nabla \cdot (κ (x) \nabla u) + p (x, t),

where u represents the distribution of temperature evolving with time t > 0, at each point x = (x,y) in Ω. Moreover, κ is the thermal conductivity (W.m⁻¹.K⁻¹ i.e. watt per meter kelvin in SI units), ρ is the density (kg.m⁻³ i.e. kilogram per cubic meter in SI units) and c the specific heat (or thermal) capacity (J.K⁻¹.kg⁻¹ i.e. joule per kilogram kelvin in SI units). It is customary to let κ go in front of the spatial derivatives when the medium is homogeneous. However, here we consider a heterogeneous medium, hence the spatial derivatives of κ might suffer some discontinuity (derivatives are taken in distributional sense, hence transmission conditions ensuring continuity of the heat flux κ∇u are encompassed in Eq. (1)). In this paper, we consider a source with a time step (Heaviside) variation and a singular (Dirac) spatial variation, that is: p(x,t) = p₀H(t)δ(x−x0), with H the Heaviside function and Delta the Dirac distribution. This means that the source term is constant throughout time t > 0, while it is spatially localized on the line x = x₀.

Upon a change of variable x = (x,y) → x′ = (x′,y′) described by a Jacobian matrix J = ∂(x′,y′)/∂(x,y), this equation takes the form:

ρ (x^{'}) c (x^{'}) det (J) \frac{\partial u}{\partial t} = \nabla \cdot (J^{- T} κ (x^{'}) J^{- 1} det (J) \nabla u) + det (J) p (x^{'}, t) .

We note that Eq. (1) and Eq. (2) have the same structure, except that the transformed conductivity

\underline{\underline{κ^{'}}} = J^{- T} κ J^{- 1} det (J) = κ J^{- T} J^{- 1} det (J) = κ T^{- 1},

is matrix-valued, with T the metric tensor, and the product of density by heat capacity in the left hand side is now multiplied by the determinant of the Jacobian matrix J of the transformation. An elegant way to derive Eq. (2) is to multiply Eq. (1) by a smooth function

ϕ \in C_{0}^{\infty} (Ω)

(that is an infinitely differentiable function with a compact support on Ω) and to further integrate by parts, which leads to the following variational form:

\int_{Ω} ρ c \frac{\partial u}{\partial t} ϕ dxdy + \int_{Ω} (\nabla ϕ \cdot κ \nabla u) dxdy - \int_{\partial Ω} (κ \nabla u \cdot n ϕ - κ \nabla ϕ \cdot n u) dl + < p, ϕ > = 0,

where n is the unit outward normal to the boundary ∂Ω of the integration domain Ω. Moreover, <,> denotes the duality product between the space of Distributions (𝒟′(Ω)) and the space of smooth functions (

𝒟 (Ω) = C_{0}^{\infty} (Ω)

) i.e. a pairing in which one integrates a distribution against a test function, so denoted because it corresponds to a bilinear map from vector spaces to scalars, see [21] and references therein.

We now apply to Eq. (4) the coordinate change x = (x,y) → x′ = (x′,y′) and noting that ∇ = J⁻¹∇′, where ∇′ is the gradient in the new coordinates, we end up with

\begin{array}{l} \int_{Ω} {(ρ c \frac{\partial u}{\partial t} ϕ) det (J)} d x^{'} d y^{'} + \int_{Ω} {(J^{- 1} \nabla^{'} ϕ \cdot κ J^{- 1} \nabla^{'} u) det (J)} d x^{'} d y^{'} \\ - \int_{\partial Ω} {(κ J^{- 1} \nabla^{'} u \cdot n ϕ - κ J^{- 1} \nabla^{'} ϕ \cdot n u) det (J)} d l^{'} + < det (J) p, ϕ > = 0 \end{array}

Upon integration by parts, and noting that J⁻¹∇′ϕ · κJ⁻¹∇′u = (∇′ϕ)^T J⁻^T κJ⁻¹∇′u we obtain the variational form of Eq. (2) which lays the foundation of transformation thermodynamics. Let us now apply the transformed Eq. (2) to the design of two metamaterials of potential pratical interest: an invisibility cloak and a concentrator for heat flux, as shown in Fig. 1.

Fig. 1 Principle of an invisibility cloak and a concentrator for heat flux. (a) Sketch of the transformed metric (dotted lines) for an invisibility cloak of inner and outer boundaries described by radii R₁(θ) and R₂(θ) depending upon the angular variable θ. The metamaterial consists of a coating (in blue) with anisotropic heterogeneous conductivity as described in section 3.1. The invisibility region is shown in red; (b) Sketch of the transformed metric (dotted lines) for a concentrator of inner and outer boundaries described by radii R₁(θ) and R₃(θ) depending upon the angular variable θ. The boundary R₂(θ) (in pink) is virtual and only serves the purpose of the geometric transform. The metamaterial consists of an inner region (in red) and a coating (in blue) filled with anisotropic heterogeneous conductivities as described in section 3.2.

Download Full Size | PDF

3. Transformed based thermic metamaterials

Following the work by Pendry et al. [2], we consider the linear geometric transform:

{\begin{array}{l} r^{'} & = & α (θ) r + β (θ) \\ θ^{'} & = & θ \end{array}

which is simply a radial stretch of polar coordinates, but also depending upon the azimuthal angle [22] in order to model domains of a (smooth) arbitrary shape.

After some elementary algebra, we obtain the following transformation matrix:

T^{- 1} = R (θ) (\begin{matrix} \frac{{(r^{'} - β)}^{2} + c_{12}^{2} . α^{2}}{(r^{'} - β) r^{'}} & - \frac{c_{12} . α}{r^{'} - β} \\ - \frac{c_{12} α}{r^{'} - β} & \frac{r^{'}}{r^{'} - β} \end{matrix}) R {(θ)}^{T},

where c₁₂ = ∂r/∂θ′ and R(θ) denotes the rotation matrix through an angle θ, which amounts to expressing the metric tensor in a Cartesian coordinate basis. The details of the derivation can be found in [22], where a similar analysis was performed for control of transverse time-harmonic electromagnetic waves propagating in cylindrical invisibility cloaks.

One should note that when r′ = β in the transformed coordinates, i.e. r = 0 in the original coordinates, the transformation matrix becomes singular as its first coefficient vanishes and the other three tend to infinity.

3.1. An invisibility cloak for heat

Let us consider two embedded smooth closed curves defined by R₁(θ) < R₂(θ), for 0 ≤ θ < 2π, as shown in Fig. 1 and the transformation Eq. (6) with the parameters $α = \frac{R_{2} (θ) - R_{1} (θ)}{R_{2} (θ)}$ and β = R₁(θ): It maps the field within the domain r ≤ R₂(θ) onto the annular domain R₁(θ) ≤ r′ ≤ R₂(θ). The region r′ ≤ R₁(θ) defines the invisibility zone, while the annulus in the cloak consists of a material with heterogeneous anisotropic conductivity $\underline{\underline{κ^{'}}} = κ T^{- 1}$ given by Eq. (7) with

\begin{array}{l} c_{12} & = \frac{- 1}{{(R_{2} (θ) - R_{1} (θ))}^{2}} {\frac{d R_{1} (θ)}{d θ} (R_{2} (θ) - r^{'}) {(R_{2} (θ) - R_{1} (θ))}^{2} R_{2} (θ) \\ + \frac{d R_{2} (θ)}{d θ} R_{1} (θ) (R_{1} (θ) - r^{'})} \end{array}

We now note that if the cloak is circular, dR₁(θ)/dθ = dR₂(θ)/dθ = 0 so that c₁₂ vanishes. It is enlightening to derive the transformation matrix T⁻¹ in that case as this highlights the algorithm underpining the design of thermal metamaterials.

We start with the computation of the Jacobian matrix of the compound transformation (x,y) → (r,θ) → (r′θ′) → (x′,y′) which can be expressed as:

\begin{array}{l} J_{x x^{'}} = J_{x r} J_{r r^{'}} J_{r^{'} x^{'}} & = R (θ) diag (1, r) diag (α^{- 1}, 1) diag(1, 1 / r^{'}) R (θ^{'}) \\ = R (θ) diag (α^{- 1}, r / r^{'}) R (θ^{'}), \end{array}

where J_xr = ∂(x,y)/∂(r,θ), J_rr_′ = ∂(r,θ)/∂(r′,θ′) and J_r_′_x_′ = ∂(r′,θ′)/∂(x′,y′); Moreover, α = (R₂ − R₁)/R₂ and we note that det(J_xx_′) = α⁻¹r/r′ as the rotation matrices R(θ) and R(θ′) are unimodular.

We then compute the transformation matrix using Eq. (3):

\begin{array}{l} T^{- 1} & = J_{x x^{'}}^{- T} J_{x x^{'}}^{- 1} det (J_{x x^{'}}) \\ = R {(θ^{'})}^{- T} diag (α, r^{'} / r) R {(θ)}^{- T} R {(θ^{'})}^{- 1} diag(α, r^{'} / r) R {(θ)}^{- 1} α^{- 1} r / r^{'} \\ = R (θ^{'}) diag (α r / r^{'}, α^{- 1} r^{'} / r) R {(θ^{'})}^{- 1}, \end{array}

where we have use the fact that the rotation matrix satifies R(θ)⁻¹ = R(θ)^T and θ′ = θ.

We now observe that from Eq. (6)r = (r′ − R₁)/α, hence the transformed conductivity inside the circular coating of the cloak can be expressed as

\underline{\underline{κ^{'}}} = κ T^{- 1} = R (θ^{'}) diag(κ_{r^{'}}^{'}, κ_{θ^{'}}^{'}) R {(θ^{'})}^{T},

where the eigenvalues of the diagonal matrix (principal values of conductivity) are from Eq. (10)

\begin{matrix} κ_{r^{'}}^{'} = \frac{r^{'} - R_{1}}{r^{'}}, & κ_{θ^{'}}^{'} = \frac{r^{'}}{r^{'} - R_{1}}, \end{matrix}

with R₁ and R₂ the interior and the exterior radii of the cloak.

We report in Fig. 2 some finite element computations with COMSOL MULTIPHYSICS for a circular thermal invisibility cloak when we apply a constant temperature normalized to 1 on the left side of a cell and a convective flux boundary condition on the right side i.e. n · κ∇u = h(u_R − u), where u_R is the free temperature on the right side, u = 0 that of the ambient medium, and h = 5W.m⁻².K⁻¹ a convective coefficient. We stress that the temperature is small but does not vanish on the right side of the cell. However, the invisibility region (inner disc) displays a specific protection for heat, as the field is constant there. This comes from the fact that from Eq. (2) taken in distributional sense, we are ensured that u and $n \cdot \underline{\underline{κ^{'}}} \nabla u$ are continuous across the cloak’s inner circular boundary, where the unit outward normal n to the boundary is directed along the radial axis. We thus note that the radial flux κ′_r_′∂u/∂r′ is continuous across this boundary. However, its trace on the outer boundary is null as from Eq. (12) κ′_r_′∂u/∂r′ = (r′ − R₁)/r′∂u/∂r′, which vanishes when r′ = R₁. This means that by continuity, the flux which is simply κ∇u inside the invisibility region (the conductivity is a constant equal to the value outside the cloak) should also vanish there. It is therefore natural to have a constant temperature inside the invisibility region. Moreover, the constant value is allowed to vary with time, which is also in agreement with the four panels in Fig. 2.

Fig. 2 Diffusion of heat from the left on a cloak with R₁ = 2.10⁻⁴m and R₂ = 3.10⁻⁴m. The temperature is normalized throughout time on the left side of the cell. Snapshots of temperature distribution at t = 0.001s (a), t = 0.005s (b), t = 0.02s (c), t = 0.05s (d). Streamlines of thermal flux are also represented with white color in panel (d). The mesh formed by streamlines and isothermal values illustrates the deformation of the transformed thermal space: the central disc (‘invisibility region’) is a hole in the metric, which is curved smoothly around it.

Download Full Size | PDF

For t > 0.005s, temperature is actually not vanishing in the inner disc, see panel (c) in Fig. 2 for t = 0.02s. This is an effect due to the non-harmonic regime. Interestingly, its value increases monotonically with time until t = 0.05s. We checked that when t > 0.05s, the normalized temperature no longer increases inside the invisibility region and does not exceed T/2, with T = 1 the temperature on the left side of the cell. Moreover, the isovalues of the field showing the diffusion of temperature within the metamaterial cloak clearly demonstrate that temperature inside the coating is smoothly detoured around the invisibility region, while going unperturbed elsewhere. Such a cloak could be used in thermal protection.

We note that on the inner boundary of the cloak, i.e. when r′ = R₁, the conductivity becomes infinitely anisotropic: the radial coefficient κ′_r_′ vanishes while the azimuthal coefficient κ′_θ_′ tends to infinity: this means that the heat flux is detoured around the invisibility region. Moreover, the left-hand member of Eq. (2) vanishes i.e. the transformed heat equation is thermo-static on the inner boundary of the cloak. Indeed, using again r = (r′ − R₁)/α, we note that the transformed density and heat capacity in Eq. (2) are such that

ρ^{'} c^{'} = det (J_{r r^{'}}) ρ c = \frac{r}{r^{'}} \frac{ρ c}{α} = \frac{r^{'} - R_{1}}{r^{'}} {(\frac{R_{2}}{R_{2} - R_{1}})}^{2} (ρ c),

which vanishes for r′ = R₁. Physically, this means it should take an infinite amount of time for heat to cross the inner boundary of the cloak. This is consistent with the fact that 0 ≤ det(J) ≤ R₂/(R₂ − R₁) (as R₁ ≤ r′ ≤ R₂), which can be seen as a time-scaling coefficient according to Eq. (2), corresponding to the blow up of the point r = 0 onto the disc r′ ≤ R₂. However, the boundary condition at r′ = R₁ can only be satisfied approximately by the finite element package, and one can see that temperature increases steadily inside the invisibility region until it reaches a plateau for t > 0.05s. The temperature is then constant inside the invisibility region.

3.2. A concentrator for heat flux

Let us consider three embedded smooth closed curves defined by R₁(θ) < R₂(θ) < R₃(θ), for 0 ≤ θ < 2π, as shown in Fig. 1 and the transformation Eq. (6) with $α = \frac{R_{1} (θ)}{R_{2} (θ)}$ and β = 0 for 0 ≤ r ≤ R₂(θ) and with $α = \frac{R_{3} (θ) - R_{1} (θ)}{R_{3} (θ) - R_{2} (θ)}$ and $β = R_{3} (θ) \frac{R_{1} (θ) - R_{2} (θ)}{R_{3} (θ) - R_{2} (θ)}$ for R₂(θ) ≤ r ≤ R₃(θ). This transformation maps the field in the region 0 ≤ r ≤ R₂(θ) onto 0 ≤ r′ ≤ R₁(θ′) (i.e. compression of thermal space) and the field in the region R₂(θ) ≤ r ≤ R₃(θ) onto R₁(θ′) ≤ r′ ≤ R₃(θ′) (i.e. extension of thermal space). Importantly, the compression and extension compensate each other for 0 < r′ ≤ R₃(θ′) and the transformation should be the identity from R₃(θ) < r to R₃(θ′) < r′. The resulting heat concentrator (in transformed variables) then consists of two parts:

0 ≤ r′ ≤ R₁(θ′) (inner core): $c_{12} = \frac{\partial r}{\partial θ^{'}} = - \frac{r}{R_{1} {(θ)}^{2}} (R_{2} (θ) \frac{\partial R_{1} (θ)}{\partial θ} - R_{1} (θ) \frac{\partial R_{2} (θ)}{\partial θ}),$
R₁(θ′) ≤ r′ ≤ R₃(θ′) (coating): $\begin{array}{l} c_{12} & = \frac{- 1}{{(R_{2} (θ) - R_{1} (θ))}^{2}} {\frac{\partial R_{1} (θ)}{\partial θ} (R_{3} (θ) - r) (R_{3} (θ) - R_{2} (θ)) \\ + \frac{\partial R_{2} (θ)}{\partial θ} (R_{1} (θ) - R_{3} (θ)) (R_{3} (θ) - r) - \frac{\partial R_{3} (θ)}{\partial θ} (R_{2} (θ) - R_{1} (θ)) (R_{1} (θ) - r)} . \end{array}$

We now note that if the concentrator is circular, dR₁(θ)/dθ = dR₂(θ)/dθ = dR₃(θ)/dθ = 0 so that c₁₂ vanishes in Eq. (14) and Eq. (15). The derivation of the T matrix in that case follows mutatis mutandis that of the previous section with (α,β) = (R₁/R₂, 0) if 0 ≤ r ≤ R₂ and (α,β) = ((R₃ − R₁)/(R₃ − R₂), R₃(R₁ − R₂)/(R₃ − R₂)) if R₂ ≤ r ≤ R₃.

This leads us to T⁻¹ = R(θ′)Diag(κ′_r_′, κ′_θ_′)R(θ′)^T where the cloak is described by the following parameters

\begin{array}{l} κ_{r^{'}}^{'} & = 1, & κ_{θ^{'}}^{'} = 1, & if 0 \leq r^{'} \leq R_{2} \\ κ_{r^{'}}^{'} & = \frac{r^{'} + R_{3} \frac{R_{2} - R_{1}}{R_{3} - R_{2}}}{r^{'}}, & κ_{θ^{'}}^{'} = \frac{r^{'}}{r^{'} + R_{3} \frac{R_{2} - R_{1}}{R_{3} - R_{2}}}, & if R_{2} \leq r^{'} \leq R_{3} \end{array}

where R₁ and R₂ are the interior and the exterior radii of the cloak. Importantly, κ′_r_′ and κ′_θ_′ can no longer vanish or become infinite: this is due to the fact that we do not tear apart the metric when we design a concentrator, unlike for an invisibility cloak. There is a one-to-one correspondance between the boundary of the circle of radius r = R₂ and that of radius r′ = R₁ under the geometric transform (isomorphism).

Moreover, the left-hand side of Eq. (2) has the following transformed density and heat capacity:

ρ^{'} c^{'} = det (J) ρ c = {\begin{array}{l} {(\frac{R_{2}}{R_{1}})}^{2} ρ c & , if 0 \leq r^{'} \leq R_{1}, \\ \frac{r^{'} + R_{3} \frac{R_{2} - R_{1}}{R_{3} - R_{2}}}{r^{'}} {(\frac{R_{3} - R_{2}}{R_{3} - R_{1}})}^{2} ρ c & , if R_{1} \leq r^{'} \leq R_{3}, \end{array}

which is a product that cannot vanish. Here again, det(J) can also be seen as a time-scaling coefficient according to Eq. (2).

We report in Fig. 3 some finite element computations for a circular thermal concentrator when we apply a constant heat source on the left side of a cell. We note that temperature isovalues are bent towards the center of the concentrator: heat diffuses faster in the left part of the device and slower in the right part. Such a concentrator could be used to focus heat in a tiny region, what might have potential applications in photovoltaics.

Fig. 3 Diffusion of heat from the left on a concentrator with R₁ = 1.10⁻⁴m, R₂ = 2.10⁻⁴m and R₃ = 3.10⁻⁴m. The temperature is normalized throughout time on the left side of the cell. Snapshots of heat distribution at t = 0.002s (a), t = 0.005s (b), t = 0.01s (c) and t = 0.02s (d). Streamlines of thermal flux are also represented with white color in panel (d). The mesh formed by streamlines and isothermal values illustrates the deformation of the transformed thermal space: it is squeezed into the central disc.

Download Full Size | PDF

4. A broadband multilayered cloak via homogenization of the heat equation

We now wish to propose a practical design of a broadband cloak consisting of a large number of homogeneous layers with piecewise constant diffusivity (ratio of the conductivity by the product of density by specific heat capacity, which is in SI units of square meter per second). For this, we need to introduce the notion of anisotropic diffusivity.

4.1. Reduced set of parameters for anisotropic diffusivity

We first note that the coordinate transformation $r^{'} = R_{1} + r \frac{R_{2} - R_{1}}{R_{2}}$ can also compress the region r < R₂ into the shell R₁ < r < R₂, provided that the cloak is described by the following reduced anisotropic conductivity $\underline{\underline{κ^{'}}}$ (written in its diagonal basis):

\begin{matrix} κ_{r}^{″} = {(\frac{R_{2}}{R_{2} - R_{1}})}^{2} {(\frac{r^{'} - R_{1}}{r^{'}})}^{2}, & κ_{θ}^{"} = {(\frac{R_{2}}{R_{2} - R_{1}})}^{2}, \end{matrix}

where R₁ and R₂ are the interior and the exterior radii of the cloak. This simplification amounts to multiplying κ_r_′ and κ_θ_′ in Eq. (12) by the spatially varying determinant in Eq. (13). By doing so, the determinant is now equal to 1 in Eq. (2), and the rank-2 tensor

\underline{\underline{K}} = \underline{\underline{κ^{″}}} / (ρ c)

is in SI units of diffusivity (m².s⁻¹).

However, such a simplified set of parameters leads to an approximation of Eq. (2) as

\nabla \cdot (\underline{\underline{κ^{″}}} \nabla u) = det (J) \nabla \cdot (\underline{\underline{κ^{'}}} \nabla u) + \underline{\underline{κ^{'}}} \nabla \cdot (det (J) \nabla u) .

An approximation of this type is legitimate provided that the perturbation

\underline{\underline{κ^{'}}} \nabla \cdot (det (J) \nabla u)

is small enough, that is when det(J) ≪ 1. This condition breaks down near the outer boundary of the cloak, which induces an impedance mismatch in the diffusivity between the cloak and the background. However, we shall see in the sequel that the thermal protection is preserved.

4.2. Homogenization model for the heat equation

We still face the obstacle of an anisotropic conductivity inside the cloak. This problem can be handled using techniques of homogenization. Let us consider that the temperature field is solution of the following governing equation (outside the source)

ρ_{η} c_{η} \frac{\partial u_{η}}{\partial t} = \nabla \cdot (κ_{η} \nabla u_{η}),

inside the heterogeneous isotropic cloak Ω_f, where

κ_{η} = κ (\frac{r}{η})

and

ρ_{η} c_{η} = ρ (\frac{r}{η}) c (\frac{r}{η})

. When heat penetrates the structured cloak Ω_f, it undergoes fast periodic variations. To filter these variations, we consider an asymptotic expansion of u solution of Eq. (20) in terms of a macroscopic (or slow) variable x = (r,θ) and a microscopic (or fast) variable

x_{η} = (\frac{r}{η}, θ)

, where η is a small positive real parameter. The homogenization of this parabolic equation can be derived by considering the ansatz

u_{η} (x) = \sum_{i = 0}^{\infty} η^{i} u^{(i)} (x, x_{η}),

where for every x ∈ Ω_f, u⁽ⁱ⁾(x,·) is 1-periodic along r. Note that we evenly divide Ω_f (R₁ ≤ r ≤ R₂, 0 ≤ θ < 2π) into a large number of thin curved layers of radial length (R₂ − R₁)/η, but the time variable t in Eq. (20) is not rescaled.

Rescaling the differential operator in Eq. (20) accordingly as $\nabla = \nabla_{x} + e_{r} \frac{1}{η} \frac{\partial}{\partial r}$ , and collecting the terms sitting in front of the same powers of η (see e.g. [21]), we obtain

< ρ c > \frac{\partial u_{0}}{\partial t} = \nabla \cdot (\underline{\underline{κ_{hom}}} \nabla u_{0}),

which is the homogenized heat equation where u₀ is the leading order term in the ansatz of Eq. (21). Here,

< f > = \int_{0}^{1} f (r) d r

is the arithmetic mean over a unit cell along the radial axis and

\underline{\underline{κ_{hom}}}

is a homogenized rank-2 diagonal tensor which has the physical dimensions of a homogenized anisotropic conductivity

\underline{\underline{κ_{hom}}} = Diag(κ_{r}, κ_{θ})

given by

\underline{\underline{κ_{hom}}} = Diag (< κ^{- 1} >^{- 1}, < κ >) .

We note that if the cloak consists of an alternation of two homogeneous isotropic layers of thicknesses d_A and d_B and conductivities κ_A, κ_B, densities ρ_A, ρ_B and heat capacities c_A, c_B, we obtain the following effective parameters

\begin{matrix} \frac{1}{κ_{r}} = \frac{1}{1 + η} (\frac{1}{κ_{A}} + \frac{η}{κ_{B}}), & κ_{θ} = \frac{κ_{A} + η κ_{B}}{1 + η}, < ρ c > = \frac{ρ_{A} c_{A} + η ρ_{B} c_{B}}{1 + η}, \end{matrix}

which are respectively described by one harmonic and two arithmetic means, where η = d_B/d_A is the ratio of thicknesses for layers A and B and d_A + d_B = 1.

To mimic the spatially varying anisotropic diffusivity $\underline{\underline{K}} = \underline{\underline{κ^{″}}} / (ρ c)$ where $\underline{\underline{κ^{″}}}$ is given in Eq. (18), we proceed in two steps, following [23]: In a first step, we normalise the homogenized conductivity $\underline{\underline{κ_{hom}}}$ by <ρc>, what we call anisotropic homogeneous diffusivity $\underline{\underline{K_{hom}}} = \underline{\underline{κ_{hom}}} / < ρ c >$ . This approach enables us to consider layers of diffusivities K_A = κ_A/(ρ_Ac_A) and K_B = κ_B/(ρ_Bc_B), what simplifies the numerical implementation. We then approximate the ideal cloak by a multi-layered cloak with M such anisotropic homogeneous concentric layers. In a second step, we approximate each layer i, i = 1,..,M by N thin isotropic layers through the homogenization process described above. This means the overall number NM of isotropic layers can be fairly large, in our case N = 2 and M = 10. A mathematical justification for this can be found in the context of acoustic and Schrödinger equations in [10] by applying the homogenization technique to $κ_{η} = κ (r, \frac{r}{η})$ and $ρ_{η} c_{η} = ρ (r, \frac{r}{η}) c (r, \frac{r}{η})$ , so that the homogenized coefficients in Eq. (22) are still spatially varying on the slow scale r, but this leads to mathematical technicalities beyond the scope of the present paper.

4.3. Ilustrative numerical results

We show in Fig. 4 four snapshots of heat diffusing through a multilayered cloak working via homogenization. These plots have the same features as in Fig. 2. It should be noted that the diffusivity varies over a large range of values. However, this is within an easy reach with realistic materials, as one can rescale diffusivities in caption of Fig. 2 according to, say, Polyvinyl Chloride (PVC) for the bulk material (K = 8 10⁻⁸m².s⁻¹), in which case silver could be the inner most layer of the cloak (K = 1.65 10⁻⁴m².s⁻¹). Snapshots of heat distribution at t = 0.01s (a) and t = 0.05s (b) show that normalized temperature on the left side of the cell is non vanishing inside the inner disc (invisibility region), but it never exceeds a half of the maximum value of normalized temperature (even for t > 0.05s). The fact that the temperature inside the inner disc is still a constant for a given time as in Fig. 2, whereas in Fig. 4 the cloak no longer displays singular values of the diffusivity, requires a more sophisticated explanation (i.e. it is not legitimate to invoke a thermostatic problem as in section 3.1). In order to understand why this is the case, one needs to invoke a maximum principle for parabolic equations [25] that lies outside the scope of the present study.

Fig. 4 Diffusion of heat from the left on a multilayered cloak of inner radius R₁ = 1.5 10⁻⁴m and outer radius R₂ = 3 10⁻⁴m, consisting of 20 homogeneous layers of equal thickness and of respective diffusivities (in unit of 10⁻⁵.m².s⁻¹) 1680.70, 0.25, 80.75, 0.25, 29.39, 0.25, 16.37, 0.25, 10.99, 0.25, 8.18, 0.25, 6.50, 0.25, 5.40, 0.25, 4.63, 0.25, 4.06,0.25 in a bulk material of diffusivity 10⁻⁵m².s⁻¹ (inner disc and homogeneous region outside the concentric layers). Importantly, one layer in two has the same diffusivity. Snapshots of heat distribution at t = 0.01s (a) and t = 0.05s (b) show that normalized temperature is non vanishing inside the inner disc (invisibility region), but it never exceeds a half of the maximum value of temperature (even for t > 0.05s). We note the similarities with Fig. 2.

Download Full Size | PDF

5. Conclusion

In conclusion, we have studied analytically and numerically the extension of transformation optics to the domain of thermodynamics. We have proposed to design an invisibility cloak, which reduces the temperature inside an arbitrary region, and a concentrator which focuses heat flux in a given region. We stress that the efficiency of the thermal protection with the invisibility cloak depends upon the position of its center, because the maximum value of the constant field inside the invisibility region corresponds to the value of the temperature at this point before the geometric transform is applied. We have finally proposed a design of multi-layered cloak consisting of homogeneous isotropic layers of materials with realistic diffusivities. Applications may adress the field of microelectronics, for which reason the diameters of metamaterials are 400 to 600 micrometers and the time scale considered is 20 to 50 milliseconds in our numerical examples. However the key parameter is the ratio of time to length, so that the scale can be directly modified for most general thermodynamics applications. Notice here that we assumed the materials not to be temperature dependent, but this can be directly included within our approach. Moreover, analysis of cloaking for other types of heat sources, such as instantaneous ones, or harmonic ones, is also a possible extension of this work.

Finally, one should note that there is a strong analogy between the heat Eq. (1) and the time-dependent Schrödinger equation in the context of quantum mechanics. Some time-independent quantum cloaks for matter waves have been proposed by the groups of Zhang and Greenleaf using coordinate transforms [10, 11], and space foldings have also been investigated in this context [12]. Such results would underpin cloaking in the context of time-harmonic heat equation. However, in the present case, the transformed heat Eq. (2) is rather analogous to the time dependent transformed Schrödinger equation

det (J) (i \bar{h} \frac{\partial Ψ}{\partial t}) = - \frac{{\bar{h}}^{2}}{2} \nabla \cdot (J^{- T} m^{- 1} (x^{'}) J^{- 1} det (J) \nabla Ψ) + det (J) V Ψ,

which has not been studied thus far in the cloaking literature. Here, what plays the role of the transformed conductivity in Eq. (2) is

{\underline{\underline{m^{'}}}}^{- 1} = J^{- T} m^{- 1} (x^{'}) J^{- 1} det (J)

, which is a transformed (anisotropic heterogeneous) mass, whereas the coefficient in front of the source term in Eq. (2) is now replaced by det(J)V, which is a transformed (heterogeneous) potential. Moreover, det(J) in the left-hand side is a time-scaling parameter, h̄ is the Planck constant and Ψ is the wave function. In order to facilitate the physical interpretation of Eq. (25), one can e.g. divide throughout by det(J), and consider a reduced set of parameters for the transformed mass.

Acknowledgments

S.G. is thankful for an ERC funding through ANAMORPHISM project. The authors are also grateful for stimulating discussions with M. Zerrad and for insightful comments from the anonymous reviewer on the general form of the heat equation and on useful analogies with the Schrödinger equation, which helped improving the manuscript.

References and links

1. R. K. Luneburg, Mathematical Theory of Optics (University of California Press, 1964).

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

3. U. Leonhardt, “Optical conformal mapping,” Science 312, 1777–1780 (2006). [CrossRef] [PubMed]

4. A. Greenleaf, Y. Kurylev, M. Lassas, and G. Uhlmann, “Full-wave invisibility of active devices at all frequencies,” Commun. Math. Phys. 275(3) 749–789 (2007). [CrossRef]

5. R. V. Kohn, H. Shen, M. S. Vogelius, and M. I. Weinstein, “Cloaking via change of variables in electric impedance tomography,” Inverse Probl. 24015016 (2008). [CrossRef]

6. D. Schurig, J. J. Mock, B. J. Justice, S. A. Cummer, J. B. Pendry, A. F. Starr, and D. R. Smith, “Metamaterial electromagnetic cloak at microwave frequencies,” Science 314, 977–980 (2006). [CrossRef] [PubMed]

7. B. Kanté, D. Germain, and A. de Lustrac, “Experimental demonstration of a nonmagnetic metamaterial cloak at microwave frequencies,” Phys. Rev. B 80, 201104(R) (2009). [CrossRef]

8. F. Zolla, S. Guenneau, A. Nicolet, and J. B. Pendry, “Electromagnetic analysis of cylindrical invisibility cloaks and the mirage effect,” Opt. Lett. 32, 1069–1071 (2007). [CrossRef] [PubMed]

9. A. Greenleaf, M. Lassas, and G. Uhlmann, “On nonuniqueness for Calderon’s inverse problem,” Math. Res. Lett. 10, 685–693 (2003).

10. A. Greenleaf, Y. Kurylev, M. Lassas, and G. Uhlmann, “Isotropic transformation optics: approximate acoustic and quantum cloaking,” New J. Phys. 10, 115024 (2008). [CrossRef]

11. S. Zhang, D. A. Genov, C. Sun, and X. Zhang, “Cloaking of matter waves,” Phys. Rev. Lett. 100, 123002 (2008). [CrossRef] [PubMed]

12. A. Diatta and S. Guenneau, “Non singular cloaks allow mimesis,” J. Opt. 13, 024012 (2011). [CrossRef]

13. A. Greenleaf, Y. Kurylev, M. Lassas, and G. Uhlmann, “Schrödinger’s Hat: Electromagnetic, acoustic and quantum amplifiers via transformation optics,” (preprint:arXiv:1107.4685v1).

14. S. A. Cummer and D. Schurig, “One path to acoustic cloaking,” New J. Phys. 9, 45–45 (2007). [CrossRef]

15. H. Chen and C. T. Chan, “Acoustic cloaking in three dimensions using acoustic metamaterials,” Appl. Phys. Lett. 91, 183518 (2007).

16. A. Norris, “Acoustic cloaking theory,” Proc. R. Soc. London 464, 2411–2434 (2008). [CrossRef]

17. G. W. Milton, M. Briane, and J. R. Willis, “On cloaking for elasticity and physical equations with a transformation invariant form,” New J. Phys. 8, 248–268 (2006). [CrossRef]

18. M. Brun, S. Guenneau, and A.B. Movchan, “Achieving control of in-plane elastic waves,” Appl. Phys. Lett. 94, 061903 (2009). [CrossRef]

19. A. Alu and N. Engheta, “Achieving transparency with plasmonic and metamaterial coatings,” Phys. Rev. E 95, 016623 (2005). [CrossRef]

20. G. W. Milton and N. A. Nicorovici, “On the cloaking effects associated with anomalous localised resonance,” Proc. R. Soc. London, Ser. A462,3027–3059 (2006). [CrossRef]

21. V. V. Jikov, S. M. Kozlov, and O. A. Oleinik, Homogenization of Differential Operators and Integral Functionals (Springer-Verlag, New-York, 1994). [CrossRef]

22. A. Nicolet, F. Zolla, and S. Guenneau, “Electromagnetic analysis of cylindrical cloaks of arbitrary cross-section,” Opt. Letters 33, 1584–1586 (2008). [CrossRef]

23. Y. Huang, Y. Feng, and T. Jiang, “Electromagnetic cloaking by layered structure of homogeneous isotropic materials,” Opt. Express 15, 11133–11141 (2007). [CrossRef] [PubMed]

24. M. Rahm, D. Schurig, D. A. Roberts, S. A. Cummer, D. R. Smith, and J. B. Pendry, “Design of electromagnetic cloaks and concentrators using form-invariant coordinate transformations of Maxwell’s equations,” Photonics Nanostruct. Fundam. Appl. 6, 87–95 (2008). [CrossRef]

25. L. Nirenberg, “A strong maximum principle for parabolic equations,” Commun. Pure Appl. Math. 6, 167–177 (1953). [CrossRef]

Transformation thermodynamics: cloaking and concentrating heat flux

Abstract

1. Introduction

2. Transformed heat equation

3. Transformed based thermic metamaterials

3.1. An invisibility cloak for heat

3.2. A concentrator for heat flux

4. A broadband multilayered cloak via homogenization of the heat equation

4.1. Reduced set of parameters for anisotropic diffusivity

4.2. Homogenization model for the heat equation

4.3. Ilustrative numerical results

5. Conclusion

Acknowledgments

References and links

Cited By

Figures (4)

Equations (25)

Optics Express