Search tips
Search criteria 


Logo of scirepAboutEditorial BoardFor AuthorsScientific Reports
Sci Rep. 2017; 7: 6261.
Published online 2017 July 24. doi:  10.1038/s41598-017-06565-6
PMCID: PMC5524825

Optimization-based design of heat flux manipulation devices with emphasis on fabricability


In this work, we present a new method for the design of heat flux manipulating devices, with emphasis on their fabricability. The design is obtained as solution of a nonlinear optimization problem where the objective function represents the given heat flux manipulation task, and the design variables define the material distribution in the device. In order to facilitate the fabrication of the device, the material at a given point is chosen from a set of predefined metamaterials. Each candidate material is assumed to be a laminate of materials with high conductivity contrast, so it is a metamaterial with a highly anisotropic effective conductivity. Following the discrete material optimization (DMO) approach, the fraction of each material at a given finite element of the mesh is defined as a function of continuous variables, which are ultimately the design variables. This DMO definition forces the fraction of each candidate to tend to either zero or one at the optimal solution. As an application example, we designed an easy-to-make device for heat flux concentration and cloaking.


The use of materials with unprecedented effective properties (the so-called “metamaterials”) to control the electromagnetic flux has led to major innovations in optics, electronics and communications1, 2. Taking advantages of the analogies between electromagnetism and thermodynamics3, very interesting, still academic, applications of metamaterials for heat flux manipulation were recently published, including inversion4, shielding4, 5, concentration4, 6, and cloaking4, 7. Of course, all these examples were effective to prove the potential of thermal metamaterials, but their practical application needs further research.

We identify two main obstacles to the practical use of thermal metamaterials. First, most of them47 were designed on the base of coordinate transformation, an approach inherited from electromagnetism8. Up to the authors’ knowledge, the transformation-based approach has been always applied to design metamaterials occupying geometrically simple domains, under homogeneous external heat flux, in order to accomplish a few specific tasks (shielding, concentration, inversion, cloaking). In other words, it has not been (cannot be?) applied to general problems.

To circumvent such limitations, we apply here an optimization-based methodology9, where the cost function defines an arbitrary heat flux manipulation task, to be accomplished in an arbitrary domain under arbitrary boundary conditions.

The other obstacle for real-life applications is the difficult fabrication of the devices, mainly because they need a precise metamaterial variation inside them. This was circumvented by Vemuri et al.10 by fabricating heat flux manipulating devices using a homogeneous laminate that was arranged in two different orientations at each fourth of the device. By this way, they approached as well as possible the thermal conductivity distribution required to accomplish the given task.

Agreeing with Vemuri et al., we consider the coarse discretization of the metamaterial spatial distribution to be the way to facilitate the fabrication of metamaterial devices. Then, we determine the geometry of the resulting homogeneous subdomains as solution of the problem of minimizing the error in the task accomplishment. If the material for each subdomain is chosen between only two metamaterials, this solution can be obtained by adapting techniques borrowed from topology optimization11. In case of a larger set of candidate materials, analogies can be found with multiphase topology optimization12. Particularly, we will define the fraction of each candidate material as a function of continuous design variables following the discrete material optimization (DMO) approach13, which allows to use efficient gradient-based optimization solvers.

However, the current problem has crucial differences with classical topology optimization problems. First, the objective function in topology optimization is either the material volume or the compliance11, linearly dependent on the design or the state variables, respectively. Meanwhile, the current objective represents the given task and is a highly nonlinear function of both the design and the state variables. Secondly, to consider the material volume (either as the objective function or as a constraint) is imperative for topology optimization but not in this case. Actually, it may be interesting to minimize or to limit the material volume (for instance, that of the most expensive material), but this is out of the scope of the current work.

Finally, we applied the current DMO-based approach for the design of an easy-to-make device for heat flux concentration and cloaking, as alternative to that one with continuously spatially variable metamaterial distribution we previously designed9.

Heat flux Manipulation

Steady thermal field as function of material distribution

Applying the finite element method (FEM), the temperature distribution in the domain Ω is approximated as

T(x) = N(x)T,  ∀ x ∈ Ω, 

where N i is the shape function associated to node i and T i is the temperature at this node. Assuming steady heat conduction, T is the solution of


where K and F are the conductivity matrix and the nodal heat flux vector, respectively, given by

K = ∫ΩBTkBdV

F = ∫∂ΩqqwallNdS

with Bij = ∂Nj/∂xi, such that BT = gradT, k is the effective thermal conductivity tensor, and q wall is the prescribed heat flux at the portion [partial differential]Ωq of the boundary of Ω.

Let Ω be a heterogeneous body, where the material is allowed to have an element-wise variation throughout the finite element mesh Ω = Ω (1)[union or logical sum]Ω (2)[union or logical sum][union or logical sum]Ω (E). Further, the effective material properties at the finite element Ω (e) are assumed to be determined by a finite number of parameters grouped in the vector p(e)=[p1(e),p2(e),,pM(e)], as shown in Fig. 1. So, the effective conductivity at Ω (e) is k (e) = k(p (e)). This makes the conductivity matrix K, Equation (3), a function of the parameters of all the elements of the mesh, grouped in the vector P=[p1(1),p2(1),,pM(E)]. Consequently, the nodal temperature vector T, obtained as solution of Equation (2), is also a function of P, i.e. T = T(P).

Figure 1
Heat flux manipulation problem in the domain Ω where the effective properties at each sub-domain Ω (e) depend on a set of parameters p (e).

Task accomplishment as an optimization problem

At the point x (q) inside the finite element Ω (q), the heat flux is


Now, to manipulate the heat flux inside Ω means to force the heat flux to take prescribed values q¯(q) at a series of points x (q) [set membership] Ω, q = 1, 2,…, Q, see Fig. 1. To accomplish this task, we must find P such that


The search of P is constrained to an admissible design space 𝒟. In general, it will not be possible to accomplish the given task for P [set membership] 𝒟. Then, let us do that as well as possible by solving the following nonlinear optimization problem:


where P plays the role of decision or design variables, and w q is the weight given to the accomplishment of the task at the point x (q), with qwq = 1. If w q = 1/Q = constant, the above objective function is the root mean square error (RMSE) in the accomplishment of the task.

In our previous work9, we designed a heat flux manipulating device allowing a continuous distribution of metamaterials. There, we recognized the need of further constraints to ensure the fabrication of the so-designed devices. These are the so-called manufacturing constraints in topology optimization14.

In this work, following Vemuri et al.10, we propose an alternative way of designing an easy-to-make heat flux manipulating device: it will be made of M distinct materials chosen from a predefined set; any material m = 1, 2, …, M may be a metamaterial itself, having a generally anisotropic conductivity k m.

Parametrization of effective properties

The effective conductivity at the finite element Ω(e) is defined by the mixture law:


where fm(e) is the fraction of material m at Ω(e), such that 0fm(e)1 and ifi(e)=1. For easier fabricability, the material at Ω(e) must be one of the candidates instead of a mixture of them, that is fm(e) should be either one or zero. Assuming pm(e)fm(e) as design variables, they are of integer type, so that the solution of the optimization problem (6) requires integer programming.

But integer programming is unaffordable in presence of a large number of design variables, as it is usually the case for fine enough finite element meshes. In this case, it is convenient to use a gradient-based optimization solver, which requires the design variables to be continuous. To this end, we make use of the “discrete material optimization” (DMO) technique proposed by Stegmann and Lund13, where fm(e) is assumed to be a function of M continuous variables ρi(e), which play now the role of design variables, i.e. pi(e)ρi(e).

For optimal ρi(e), fm(e) must tend to one or zero. Stegmann and Lund13 studied various choices for the function fm(e), comparing them in terms not only of their closeness to one or zero at the optimal solution but also of their effect on the convergence to this solution, recommending the definition to be used in this work:


The above equation defines the fraction of material m at the finite element e as a function of the continuous design variables ρ1(e),,ρM(e). Since this definition automatically yields ifi(e)=1, these are not longer constraints for the continuous optimization problem6, which is now subject to bound constraints only: 0ρi(e)1. Like in topology optimization11, intermediate values of fm(e) are penalized by setting p  3 in Equation (8).

Choice of candidate materials

Metamaterials made of alternating sheets (i.e. laminates) of materials with markedly distinct conductivity exhibit a highly anisotropic effective conductivity, making them a popular choice for heat flux manipulation46, 10.

Following Vemuri et al.10, in the examples to be developed in this work we will use a laminate of alternating, equally-thick sheets of copper and polydimethylsiloxane (PDMS), with isotropic conductivities k copper = 398 Wm−1K−1 and k PDMS = 0.27 Wm−1K−1. Assuming such laminate as an effective thermal medium, its principal conductivities are


where λ and τ are the principal axes in-plane and normal to the sheets, respectively. Above equation is valid if the layers have small thickness and high conductivity contrast15. Here, the conductivity contrast between layers is very high (k copper/k PDMS = 1474). Concerning the layer thickness, being the current macroscale dimensions on the order of centimeters, it should be on the order of 1 mm.

Finally, as candidate materials to build a device for heat flux manipulation, we adopt this copper-PDMS laminate with different orientations θ (angle between λ and x axes).


For the purpose of validating the proposed DMO-based methodology for the design of heat flux manipulating devices, let us reproduce the simple device for heat flux shielding proposed by Vemuri et al.10. Given a plate Ω originally filled with 40% nickel steel (with isotropic conductivity k ns = 10 Wm−1K−1), let us design a device occupying the region Ω free to block the heat flux in the central region Ω fixint (see Fig. 2). Further, let us do it as well as possible by using only two candidate materials: material 1 is the copper-PDMS laminate with θ = 45°, and material 2 is the same laminate with θ = 135°.

Figure 2
Validation example: (a) domain of analysis, (b) metamaterial used to build the device, (c) finite element mesh of Ω; the blue elements belong to the device, and the red ones have heat flux prescribed to be null.

The domain Ω is discretized using 60 × 42 = 2520 bilinear square finite elements, as shown in Fig. 2c, including 840 elements inside the device Ω free (the blue ones in that figure).

Accomplishment of the shielding task

The shielding task is prescribed by setting q¯(q)=0 at the 60 elements in Ω fixint (the red ones in Fig. 2c). Then, the minimization problem6 takes the form:


This is a non-linear optimization problem with 1680 continuous design variables grouped in the vector P=[ρ1(1),ρ2(1),,ρ2(840)], subject to the bound constraints 0  P i  1.


The non-linear constrained optimization problem (10) is solved using the IPOPT interior point algorithm16. Further, we used a density-filtering technique17 to avoid checker-board type instabilities.

Figure 3a and b show the optimal design variables ρ1(e) and ρ2(e) inside the device. Then, using Equation (8) we compute the corresponding fraction f1(e) of copper-PDMS laminate at 45° (i.e. material 1), shown in Fig. 3c; the fraction at 135° is simply f2(e)=1f1(e). The device having such material distribution accomplishes the shielding task with high accuracy: RMSE = 1.536 × 10−4||q 0||, being ||q 0|| = 7.8 kWm−2 the magnitude of the heat flux in the plate without the device.

Figure 3
Validation example: (a) optimal design variable ρ1(e), (b) optimal design variable ρ2(e), (c) fraction of laminate at 45° for optimal ρ1(e) and ρ2(e), (d) actual device, (e) temperature (colormap with isolines every ...

The material distribution in Fig. 3c is mostly free of “grey zones” (those where the material is neither of the candidates but a mixture of them). To completely avoid “grey zones”, recourse can be made to “black-and-white” filters18. In this work, a simple a posteriori “black-and-white” filter is applied: the material at Ω (e) is that having the greatest fraction. The resultant piecewise homogeneous metamaterial distribution in the device is shown in Fig. 3d, using which the temperature and heat flux distribution in the whole domain Ω are those depicted in Fig. 3e. The shielding task is still very well accomplished: RMSE = 9.454 × 10−3||q 0||.

Figure 3d is actually the DMO-designed shielding device, where each quadrant is made of the copper-PDMS laminate with either divergent or convergent orientation upstream and downstream the shielded domain Ω fixint. This design coincides with Vemuri et al.’s10, and gives unsurprisingly the best way of forcing the heat flux off Ω fixint using the two given materials.

Application to a heat concentration and cloaking problem

Now, let us apply the current methodology for the design of a device for heat concentration and cloaking as alternative to that one we designed using an optimization-based continuous metamaterial distribution9 as well as to that one designed by Chen and Lei6 using the transformation-based approach. Figure 4a shows the domain of analysis Ω, originally filled with 40% -nickel steel, with isotropic thermal conductivity k ns = 10 Wm−1K−1. At this point, the plate undergoes a homogeneous heat flux q 0 with magnitude ||q 0|| = 7.14 kWm−2 and direction +x.

Figure 4
Heat flux concentration and cloaking example: (a) domain of analysis, (b) finite element mesh of Ω; the blue elements belong to the device, while the red and the green ones have prescribed heat flux for concentration and cloaking, respectively. ...

Then, let us design a device occupying the region Ω free to accomplish the following simultaneous tasks: 1) to concentrate the heat in Ω fixint, and 2) to keep the heat flux unaltered in Ω fixext. Three candidate materials will be used, all of them consisting of the current copper-PDMS laminate but differently oriented: material 1 with θ = 45°, material 2 with θ = 135°, and material 3 with θ = 0°.

The domain Ω is discretized using 70 × 70 = 4900 bilinear square finite elements, as shown in Fig. 4b, including 1896 elements inside the device Ω free (the blue ones in this figure).

Accomplishment of the concentration and cloaking tasks

The cloaking task q¯(q)=q0 is prescribed at the 140 elements in the stripes Ω(1) [subset or is implied by] Ωfixext and Ω(2) [subset or is implied by] Ωfixext (the green ones in Fig. 4b), while the heat concentration task is forced by setting q¯(q)=qconc for the 80 elements in Ω fixint (the red ones in Fig. 4b); here, we adopted q conc = (R/r)q 0, being R/r = 5 the ratio between the outer and the inner radii in the device. Then, the minimization problem (6) takes the form:


This is a non-linear optimization problem with 5688 continuous design variables grouped in the vector P=[ρ1(1),ρ2(1),ρ3(1),,ρ3(1896)], subject to the bound constraints 0  P i  1. Using IPOPT16 as optimization solver, together with density filtering17 to avoid checkerboard instabilities, we compute first the optimal distribution of the design variables ρ1(e), ρ2(e) and ρ3(e), which are then used to compute the material fractions f1(e), f2(e) and f3(e) according to Equation (8). Finally, assuming that the material at an element is that having the highest fraction, we obtain the device depicted in Fig. 5a, using which the temperature and heat flux in the whole plate are those shown in Fig. 5b. As it can be observed in Fig. 5a, by properly orienting the laminate, the heat flux is guided as directly as possible to Ω fixint. Quantitatively, the RMSE for concentration is 0.284||q conc||. Regarding cloaking, the RMSE is similar: 0.277||q 0||.

Figure 5
Heat flux concentration and cloaking example: (a) DMO-designed device, (b) temperature (colormap with isolines every 4 °C) and heat flux (arrows) in the whole plate, (c) temperature along y = 0 for the current DMO-designed ...

The current easy-to-make DMO-designed device given by Fig. 5a performs almost as well as the hard-to-make device we designed using a continuous metamaterial distribution9 (where the thickness of the alternating copper and PDMS layers and their orientation vary throughout the device), specially for concentration. This is clearly apparent from Fig. 5c, showing the temperature distribution along the line y = 0 for both devices.

Now, after Chen and Lei6, let us consider a concentration efficiency index defined as f = (T c  T b)/(T d  T a), where b and c are points located in the boundary of the heat concentration region Ω fixint, and a and d are points located at the outer boundary of the device Ω free (see Fig. 4a). For the current device, we obtain f = 88.37%. This is slightly better than the efficiency of the device fabricated by Chen and Lei6, made of 100 radial copper-PDMS laminates.

Note that the deal between both tasks can be changed by tuning the weights w q, which were assumed to be constant in this example. Actually, cloaking is not usually the main task but most of the devices found in the literature4, 6 perform it as a collateral effect of their transformation-based design.


Having already demonstrated the potentiality of the optimization-based approach for the design of metamaterial devices for heat flux manipulation9, we faced in this work a crucial deterrent for real-life applications of such devices: their fabricability. To this end, easy-to-make heat flux manipulating devices are designed as assemblages of homogeneous subregions made of materials chosen from a predefined set of candidates. Each candidate may be a metamaterial itself, having a highly anisotropic thermal conductivity.

The optimal geometry of these subregions is determined as solution of a non-linear optimization problem where the desired heat flux manipulating task is the objective function. Then, we used the Discrete Material Optimization (DMO) approach to define the material fraction as a function of continuous design variables that rapidly tend to 0 or 1 in order to ensure that the material at a point is one of the candidates and not a mixture of them.

Making so easy to manipulate the heat flux, we hope this work shows the way for the expansion of the use of thermal metamaterials in industry.


We thank the following institutions for supporting this work: The National Scientific and Technical Research Council (CONICET) from Argentina, through the project PIP 1105 “Computational Simulation of Multiphysics Problems. Application to Metal Solidification and Micromechanical Devices”. The European Research Council (ERC), through the project “Advanced Tools for Computational Design of Engineering Materials (COMP-DES-MAT)” (FP/2007–2013, ERC Grant Agreement 320815).

Author Contributions

Author Contributions

V.D.F. developed the optimization-based methodology for general thermomechanical problems, and I.P. applied it to design heat flux manipulating devices. All the work was done under the supervision of V.D.F., who is the tutor of I.P. (doctoral fellow).


Competing Interests

The authors declare that they have no competing interests.


Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


1. Soukoulis, C. M. & Wegener, M. Past achievements and future challenges in the development of three-dimensional photonic metamaterials (review article). Nature Photonics5, 523–530.
2. Maldovan M. Sound and heat revolutions in phononics. Nature. 2013;503:209–217. doi: 10.1038/nature12608. [PubMed] [Cross Ref]
3. Kadic M, Bückmann T, Schittny R, Wegener M. Metamaterials beyond electromagnetism. Rep. Prog. Phys. 2013;76:126501. doi: 10.1088/0034-4885/76/12/126501. [PubMed] [Cross Ref]
4. Narayana S, Sato Y. Heat flux manipulation with engineered thermal materials. Phys. Rev. Lett. 2012;108:214303. doi: 10.1103/PhysRevLett.108.214303. [PubMed] [Cross Ref]
5. Narayana S, Savo S, Sato Y. Transient heat flux shielding using thermal metamaterials. Appl. Phys. Lett. 2013;102:201904. doi: 10.1063/1.4807744. [Cross Ref]
6. Chen, F. & Lei, D. Y. Experimental realization of extreme heat flux concentration with easy-to-make thermal metamaterials. Sci. Rep. 5 (2015). [PMC free article] [PubMed]
7. Schittny R, Kadic M, Guenneau S, Wegener M. Experiments on transformation thermodynamics: Molding the flow of heat. Phys. Rev. Lett. 2013;110:195901. doi: 10.1103/PhysRevLett.110.195901. [PubMed] [Cross Ref]
8. Pendry JB, Schurig D, Smith DR. Controlling electromagnetic fields. Science. 2006;312:1780–1782. doi: 10.1126/science.1125907. [PubMed] [Cross Ref]
9. Peralta I, Fachinotti V, Ciarbonetti A. Optimization-based design of a heat flux concentrator. Sci. Rep. 2017;7:40591. doi: 10.1038/srep40591. [PMC free article] [PubMed] [Cross Ref]
10. Vemuri KP, Canbazoglu FM, Bandaru PR. Guiding conductive heat flux through thermal metamaterials. Appl. Phys. Lett. 2014;105:193904. doi: 10.1063/1.4901885. [Cross Ref]
11. Bendsøe, M. P. & Sigmund, O. Topology optimization. Theory, methods, and applications (Springer-Verlag, 2003).
12. Sigmund O, Torquato S. Design of materials with extreme thermal expansion using a three-phase topology optimization method. J. Mech. Phys. Solids. 1997;45:1037–1067. doi: 10.1016/S0022-5096(96)00114-7. [Cross Ref]
13. Stegmann J, Lund E. Discrete material optimization of general composite shell structures. Int. J. Numer. Meth. Engng. 2005;62:2009–2027. doi: 10.1002/nme.1259. [Cross Ref]
14. Sigmund O. Manufacturing tolerant topology optimization. Acta Mech. Sin. 2009;25:227–239. doi: 10.1007/s10409-009-0240-z. [Cross Ref]
15. Vemuri KP, Bandaru PR. Geometrical considerations in the control and manipulation of conductive heat flux in multilayered thermal metamaterials. Appl. Phys. Lett. 2013;103:133111. doi: 10.1063/1.4823455. [Cross Ref]
16. Wächter A, Biegler LT. On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming. Math. Program., Ser. A. 2006;106:25–57. doi: 10.1007/s10107-004-0559-y. [Cross Ref]
17. Bruns TE, Tortorelli DA. Topology optimization of non-linear elastic structures and compliant mechanisms. Comput. Methods Appl. Mech. Engrg. 2001;190:3443–3459. doi: 10.1016/S0045-7825(00)00278-4. [Cross Ref]
18. Sigmund O. Morphology-based black and white filters for topology optimization. Struct. Multidisc. Optim. 2007;33:401–424. doi: 10.1007/s00158-006-0087-x. [Cross Ref]

Articles from Scientific Reports are provided here courtesy of Nature Publishing Group