Transformed Fourier and Fick equations for the control of heat and mass diffusion

We review recent advances in the control of di ff usion processes in thermodynamics and life sciences through geometric transforms in the Fourier and Fick equations, which govern heat and mass di ff usion, respectively. We propose to further encom-pass transport properties in the transformed equations, whereby the temperature is governed by a three-dimensional, time-dependent, anisotropic heterogeneous convection-di ff usion equation, which is a parabolic partial di ff erential equation combining the di ff usion equation and the advection equation. We perform two dimensional finite element computations for cloaks, concentrators and rotators of a complex shape in the transient regime. We precise that in contrast to invisibility cloaks for waves, the temperature (or mass concentration) inside a di ff usion cloak crucially depends upon time, its distance from the source, and the di ff usivity of the invisibility region. However, heat (or mass) di ff usion outside cloaks, concentrators and rotators is una ff ected by their presence, whatever their shape or posi-tion. Finally, we propose simplified designs of layered cylindrical and spherical di ff usion cloaks that might foster experimental e ff orts in thermal and biochem-ical metamaterials. C 2015 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution License.


I. INTRODUCTION
There is currently a keen interest in the control of heat flux using thermal metamaterials in steady, 1-7 transient [8][9][10][11][12] and periodic 13 regimes. In the present paper, we discuss the functionality of diffusion cloaks for heat and mass, as well as concentrators and rotators, 14 via transformed Fourier 8,15 and Fick [16][17][18] equations.
Such designs are based upon the extension of metamaterials designed using tools of transformation optics to the fields of thermodynamics and biochemistry. Let us note that cloaking for diffusion processes 19,20 is more subtle than for waves [21][22][23][24][25][26][27][28] wherein the field vanishes in the invisibility region irrespective of its material constituent, time and its distance of the source. This appears to be in sharp contrast with thermal cloaks wherein temperature inside the invisibility region (or mass concentration in the context of diffusion of chemical species) appears to depend on its diffusivity, upon time, and the distance from the source.
More precisely, in a way similar to the English physicist Pendry and his American colleagues Schurig and Smith who proposed in 2006 to design an invisibility cloak by mapping Maxwell's equations on a curvilinear space with a hole in it (where an object can be hidden), one can make coordinate changes in governing equations of thermodynamics and biochemistry 8,16 with convection diffusion phenomena.
The plan of the paper is as follows: In section II, we apply a change of coordinates to the diffusion-convection equation and show that the diffusivity becomes in general heterogeneous (i.e. spatially varying) and anisotropic (i.e. matrix valued). In section III, we apply specific mappings in order to make a hole in the transformed coordinates so as to design a cloak (with a blow-up of a point), a concentrator (with a compression of a region) and a rotator (with a rotation of axes) for heat flux or for mass concentration flux. We use finite element computations to validate our hypothesis of control of heat or concentration flux in the transformed thermal or chemical spaces. We also notably analyze the potential protection offered by diffusion cloaks, which could have potential applications in microelectronics, see Fig. 1 for an artistic view. In section IV, we propose a multilayered spherical device which approximates the functionality of a diffusion cloak in the homogenization regime. 13,29 In section V, we look at so-called diffusion carpet cloaks and we conclude the paper in section VI.

II. TRANSFORMED CONVECTION-DIFFUSION EQUATION
Let us consider the convection-diffusion equation which is a parabolic partial differential equation combining the diffusion equation and the advection equation. This equation describes physical phenomena where particles or energy (or other physical quantities) are transferred inside a physical system due to two processes: diffusion and convection. Let us assume that the diffusion coefficient and the convection velocity are constant and there are no sources or sinks in a bounded domain Ω (the source lies outside Ω). The convection-diffusion equation is then expressed as where u denotes the distribution of temperature in thermodynamics (or mass concentration in biochemistry) at each point x = (x, y, z) within the domain Ω. In the transient regime, u also changes with time t > 0 (this is markedly different from the thermostatic case, where the left-hand side of (1) vanishes). We note that the coefficient κ is the thermal conductivity (W.m −1 .K −1 i.e. watt per meter kelvin in SI units), ρ is the density (kg.m −3 i.e. kilogram per cubic meter in SI units) and c the specific heat (or thermal) capacity (J.K −1 .kg −1 i.e. joule per kilogram kelvin in SI units).
The bulk velocity v has the unit of length by time. Similar coefficients hold in the mass diffusion case and they generally depend depend upon space and time variables. The diffusion flux −κ∇u measures the amount of substance that will flow through a small volume during a small time interval (mol.m −3 .s −1 ).
One can see at this point that if the medium has an anisotropic conductivity κ, heat will diffuse in awkward directions. Moreover, if the medium is heterogeneous (e.g. when κ is piecewise constant) the spatial derivatives of κ might suffer some jump, hence one could observe strong heat gradient at interfaces. In fact, transmission conditions ensuring continuity of heat u and heat flux κ∇u are encompassed in Eq. (1) since derivatives are taken in distributional sense. Let us finally stress that in Eq. (1) the source term will be for our numerical purpose a time step (Heaviside) variation and a singular (Dirac) spatial variation, that is: p(x,t) = p 0 H(t)δ(x − x0), with H the Heaviside function and Delta the Dirac distribution. Physically speaking, the source term is constant throughout time t > 0 and it is spatially localized on the line x = x 0 .
Let us now make the change of variable (x, y, z) → (x ′ , y ′ , z ′ ) in Eq. (1). We find that the transformed convection-diffusion equation has the form: where J = ∂(x ′ , y ′ , z ′ )/∂(x, y, z) is the Jacobian matrix. One can see that Eq. (1) and (2) share the same structure. This is better seen if we introduce the transformed diffusivity and velocity as Clearly, κ ′ is a matrix-valued field, whereas v ′ is a (heterogeneous) vector field. It is interesting to note that T is the metric tensor. There are few ways to derive Eq. (2), but a simple one is to multiply (1) by a test function φ which is smooth with a compact support on Ω (i.e. it vanishes on the boundary of Ω). One can then integrate the resulting equation over Ω. After integration by parts, we are led to a weak form of Eq. (1): where <, > denotes the duality product between the space of distributions and the space of smooth functions (in which one integrates a distribution multiplied by a test function), see Ref. 29. It is easy to make the coordinate change (x, y, z) → (x ′ , y ′ , z ′ ) in (4). For this, one notes that ∇ = J −1 ∇ ′ , where ∇ ′ is the gradient in the new coordinates and also one should keep in mind there is a change in the infinitesimal volume dxdy dz proportional to det(J):  Equation (5) is the weak form of (2). In order to see this, one should integrate by parts, and note that We further note that Fick's equation which is used in the context of mass diffusion is encompassed by the convection-diffusion equation (take the product ρc to be equal to 1 and assume that the velocity field v vanishes. Thus, what we discuss in the sequel applies straightforwardly to mass diffusion.

III. DIFFUSION INVISIBILITY CLOAK, CONCENTRATOR AND ROTATOR OF AN ARBITRARY SHAPE
We now wish to implement the transformed diffusion-convection Eq. (2) in order to design three heterogeneous anisotropic media with specific functionalities: an invisibility cloak, a concentrator and a rotator for heat flux (or mass diffusion). In order to describe boundaries of general shapes, a finite Fourier expansion may be used. Throughout this section, we consider a linear transform: where α and β are θ dependent coefficients, with 0 < θ ≤ 2π. In this way, we can design some cloak, concentrator and rotator of a complex shape. For the sake of illustration, we shall consider three boundaries with up to three terms in the Fourier expansion:

A. Diffusion cloak
One can see in Figure 2 that heat flux (resp. mass in the context of biochemistry 16 or chemical engineering 17 ) is smoothly detoured around the core of the cloak. Moreover, the cloak is itself transparent for heat since the isothermal lines are unperturbed outside this metamaterial. More precisely, we consider the geometric transform (7) with parameters FIG. 2. Diffusion cloak: (Left panel) Diffusion of heat (or mass) from the left on a cloak of complex shape of inner radius R 1 (θ) = 0.4R(1 + 0.2sin(3θ)) and outer radius R 3 (θ) = R(1 + 0.2(sin(3θ) + cos(4θ))), with R = 0.4 10 −4 m. The cloak is placed close to the source (left), further away (middle) and removed (right). Snapshots of heat distribution at t = 0.001s (a,b,c), t = 0.07s (d,e,f) and t = 0.21s (g,h,i) show that isovalues of temperature (resp. mass) are nearly unperturbed outside the cloak. One notes the maximum temperature (resp. mass) in the core of the cloak is achieved for t ≥ 21s and depends upon the distance between the center of the cloak to the source.
We obtain the following Jacobian and transformation matrices (R(θ) is the matrix of rotation through an angle θ): with det(J) = (r − β)/(α 2 r) and The numerical results performed with the finite element package COMSOL are shown in Figure 2. It is interesting to examine the thermal cloak parameter values at its non-reflecting outer boundary R 2 (θ). This can be done through the analysis of the entries of the inverse of the metric tensor T in the polar basis. We first notice that the off-diagonal terms (T −1 ) r θ = (T −1 ) θr are generally nonzero unlike for the circular case when T −1 is diagonal. In general we observe that −1 < (T −1 ) r θ < 1, which reflects the rotation of the tensor T −1 with respect to its eigenbasis. We further note that (T −1 ) r r also varies with θ, unlike for circular thermal cloaks. Last, 0 < (T −1 ) r r < 1 and (T −1 ) θθ > 1.5, in agreement with the fact that the cloak ought to exhibit a strong azimuthal anisotropy for heat to flow around the inner core.
It is illuminating to compare these numerical results carried out in the intense near-field limit, to those obtained from an analytical approach, in order to better understand how this thermal cloak works. The temperature (resp. mass) field u is given in the transformed coordinates by with ρ(ρ ′ , θ ′ ) and θ(θ ′ ) given by the inverse map of the map defined by Eq. (9). We used the software Matlab to produce the map of isothermal values, see Fig. 3. Comparing this ideal thermal (resp. mass diffusion) cloaking with that of Fig. 2, we notice that the temperature (resp. mass) inside the invisibility region (inner core) of the cloak has the same value at any time step as the temperature (resp. mass) at the point that we blow up in the original coordinate system. This form of protection is markedly different from what is achieved with invisibility cloaks for waves. In the latter case, the field vanishes inside the invisibility region, whereas in the former case the field is uniform but not zero therein.

B. Diffusion concentrator
We now propose to extend the design of circular concentrators for heat 8 to more complex shapes. For the sake of simplicity, we consider a concentrator with the same boundaries as in (8). One can see in Figure 4 that heat (resp. mass concentration) isovalues are smoothly squeezed (hence the flux is enhanced) within the concentrator. Moreover, the concentrator is itself invisible for heat (resp. mass) as the isothermal lines are unperturbed outside this metamaterial. More precisely, we consider the geometric transform (7) with parameters: We obtain the same Jacobian and transformation matrices as in (10). However, the expression for entry c 22 has changed: • 0 ≤ r ≤ R 1 (θ) : . Snapshots of heat (resp. mass) distribution at t = 0.002s (a,d,g), t = 0.07s (b,e,h) and t = 0.21s (c,f,i) show that isovalues of temperature (resp. mass) are nearly unperturbed outside the cloak. One notes the maximum temperature (resp. mass) in the core of the cloak is achieved for t ≥ 21s and depends upon the distance between the center of the cloak to the source.

C. Diffusion rotator
We now extend the design of circular cylindrical rotators for heat 15 to more complex geometries. For the sake of simplicity, we consider a rotator with the same boundaries as the cloak and concentrator, see (8). One can see in Figure 5 that heat (resp. mass) flux is smoothly rotated within the rotator (through an angle π/4). Moreover, the rotator is itself invisible for heat (resp. mass) as the isothermal values are unperturbed outside this metamaterial. It could be used to enhance thermal (or mass) exchanges.
In order to design the rotator, we consider the geometric transform: where θ 0 = π/4 denotes the tilt. We obtain the following Jacobian and transformation matrices: with det(J) = r −β α 2 r . Also, we have

D. Three-dimensional diffusion cloak of a complex shape
A nontrivial question to ask is whether one can design thermal or mass diffusion cloaks of non-spherical shapes. The extension of the previous section to the general three-dimensional case requires to describe the inner and outer boundaries of the cloak with varying radii We further assume there is some revolution for the design with an arbitrary cross-section described by two functions R 1 (φ) and R 2 (φ) giving an angle dependent distance from the origin. These functions correspond respectively to the interior and exterior boundary of the cloak. We shall only assume that these two boundaries can be represented by a differentiable function. Their finite Fourier expansions are thought in the form R j (φ) = a j 0,0 +  p n=1 a j 0, n cos(nφ), j = 1, 2, where p can be a small integer, and a j 0, n = 0 for n p for computational easiness. The geometric transformation which maps the field within the full domain ρ ≤ R 2 (φ) onto the annular domain R 1 (φ) ≤ ρ ′ ≤ R 2 (φ) can be expressed as: 30 where 0 ≤ ρ ≤ R 2 (φ). Note that the transformation maps the field for ρ > R 2 (φ) onto itself through the identity transformation. This change of coordinates is characterized by the transformation of the differentials through the Jacobian: After some elementary algebra, we find that with det(J) = r −β α 2 r . Also, we have for Elsewhere, T −1 reduces to the identity matrix (c 11 = 1 and c 13 = 0 for ρ ′ > R 2 (φ ′ )).

IV. A MULTILAYERED CLOAK WITH SIMPLIFIED ISOTROPIC PARAMETERS
It is interesting to simply the design of diffusion cloaks by considering layered cloaks with homogeneous isotropic layers. As detailed in the introduction, a few papers have already proposed such a route based upon well-known formulae derived notably in the homogenization literature. 29 A. Two-dimensional layered diffusion cloak Figure 6 shows snapshots of heat (or concentration in the context of mass diffusion) distribution for a six-layer cloak and a two-layer cloak. At t = 0.02s (a,b,e,f) and t = 0.12s (c,d,g,h) isovalues of temperature (or concentration in the context of mass diffusion) are nearly unperturbed for the six-layer cloak, but not quite so for the two-layer cloak. However, normalized temperature (or concentration) is nonvanishing inside the inner disc of both types of cloaks (and slightly higher for the six-layer one) and depends upon the distance between the center of cloak to the source. The left-hand side has a constant temperature (or concentration) of 1 and the right hand side of 0.

B. Three-dimensional layered diffusion cloak
We now propose a design of a layered spherical cloak for diffusion of heat (resp. mass) based on the homogenization approach. We show in Figure 7 snapshots of heat (resp. mass concentration) for a twenty-layer diffusion cloak. Heat (resp. mass concentration) is clearly shown to be constant inside the invisibility region at every time-step.

V. DIFFUSION INVISIBILITY CARPET: MAPPING A CURVED SURFACE ON A FLAT SURFACE
Another interesting application of transformational thermodynamics is the concept of carpet cloaks, which is inspired by Li and Pendry's original proposal of invisibility carpets in electromagnetics. 31 For this, let us consider the linear geometric transform: where y ′ is a stretched vertical coordinate. It is easily seen that this linear geometric transform maps the segment [a, b] within the horizontal axis y = 0 onto the curve y ′ = y 1 (x), and leaves the curve y = y 2 (x) unchanged. Importantly, there is a one-to-one relationship between the segment and y 1 . The curves y 1 and y 2 are assumed to be differentiable, and this ensures that the carpet won't display any singularity on its inner boundary.

A. Two-dimensional carpets
The linear transform (24) is expressed in a Cartesian basis as: where α = ( y 2 − y 1 )/y 1 and from the chain rule This leads to the inverse symmetric tensor T −1 We note that the transformed product of density and specific heat has the factor det(J) = α −1 = y 1 /( y 2 − y 1 ) which is spatially varying, but non-vanishing. When the carpet's boundaries are piecewise linear, which is the case in Fig. 8, one gets piecewise constant, anisotropic, coefficients. It should be noted that the temperature (resp. mass concentration) needs to be uniform all throughout the lower boundary of the computational domain for the carpet to work properly. Moreover, when the protruding object is narrower, the required anisotropy in the carpet is larger.