Control of heat flux using computationally designed metamaterials

EngOpt 2016 - 5th International Conference on Engineering Optimization Iguassu Falls, Brazil, 19-23 June 2016. Control of Heat Flux Using Computationally Designed Metamaterials V´ıctor D. Fachinotti, Ignacio Peralta, Alfredo E. Huespe, ´ Angel A. Ciarbonetti Centro de Investigaci´on de M´etodos Computacionales (CIMEC), UNL, CONICET, Predio “Dr. Alberto Cassano”, Colectora Ruta Nacional 168 s/n, 3000, Santa Fe, Argentina, e-mail: [email protected]. Abstract To gain control over the diffusive heat flux in a given domain, one has to design metamaterials with a specific distribution of the generally anisotropic thermal conductivity throughout the domain. Until now, the appropriate conductivity distribution was usually determined using transformation thermodynamics. By this way, only a few particular cases of heat flux control in simple domains having simple boundary conditions were studied. As a more general approach, we propose to define the heat control problem as an optimization problem where we minimize the error in guiding the heat flux in a given way, taking as design variables the parameters that define the variable microstructure of the metamaterial. Anisotropic conductivity is introduced by using a metamaterial made of layers of two materials with highly different conductivities, the thickness of the layers and their orientation throughout the domain are the current design variables. As an application example we design a device that thermally shields the region it encloses, while it keeps unchanged the flux outside it. Keywords: Design of metamaterials, heat flux manipulation, optimization, thermal shield, thermal cloaking. 1 Introduction The control of the electromagnetic flux using metamaterials led to major innovations in electronics and communications [1]. Taking advantage of the analogies between electromagnetism and thermodynamics, some researchers developed materials with unprecedented thermal properties (the thermal “metamaterials”) for heat flux manipulation, for instance the heat flux inverter by Narayana and Sato [2]. Compared to the advances in electromagnetism, the design of thermal metamaterials is an emerging research and development area. In a first approach, metamaterials can be empirically designed (e.g., the thermal shield of Narayana and Sato [2]). More sophisticated thermal metamaterials can be designed using the transformation thermodynamics concept (e.g., the inverter and the concentrator proposed by Narayana and Sato [2] or the cloaking device of Schittny et al. [3], inherited from from electromagnetism [4]). A straightforward example of the application of ideas from electromagnetism in thermal problems is the heat flux inverter of Narayana and Sato [2], derived from the device to rotate electromagnetic fields proposed by Chen and Chan [5] to rotate electromagnetic fields. The transformation-based approach has been applied to specific heat control problems. For general problems (i.e., having arbitrary prescribed magnitude and direction of the heat flux, geometry of the manipulating device, geometry and boundary conditions of the domain where the device is embedded), we propose a new, optimization-based approach for the design of thermal metamaterials. We solve a nonlinear constrained optimization problem where the objective function to minimize is the error in the accomplishment of the given heat manipulation task, and the design variables characterize the spatial distribution of the metamaterial throughout the manipulating device. We show the capability of the present method by designing a device for thermal shield and cloaking, as alternative to the transformation-based design of Schittny et al.[3] 2 Heat conduction as a function of the microstructure Let us consider the domain Ω in Figure 1, with boundary ∂Ω divided in two non-overlapping portions: ∂Ωq(where the heat flux qwall is prescribed) and ∂ΩT(where the temperature Twall is prescribed). In steady state, the heat flux conduction in Ω is governed by the equation −div (kgrad T) + s= 0 in Ω,(1) 1 and the boundary conditions: T=Twall in ∂ΩT,(2) −kgrad T·n=qwall in ∂Ωq,(3) where Tis the temperature, sis the internal heat source, kis the (effective) thermal conductivity tensor, and nis the unit vector normal to and pointing outwards ∂Ω. Using the finite element method (FEM), the temperature field in Ω is approximated as follows: T(x) = Nj(x)Tj∀x∈Ω,(4) where Njis the shape function associated to the node jof the finite elements mesh (discretized Ω) and Tiis the temperature at node j(unknown). Using a standard (Galerkin) FEM, the nodal temperature Tjis the solution of the algebraic system of equations KijTj=Fi,(5) where Kij and Fiare the components of the global conductivity matrix and the nodal heat flux vector respectively, given by Kij =ZΩ grad Ni·kgrad NjdV, (6) Fi=ZΩ sNidV+Z∂Ωq qwallNidS. (7) The system of equations (5) is the FEM version of the heat conduction (1) subject to the boundary conditions (2) and (3). This is a classical FEM problem, whose details can be found for instance in the book of Zienkiewicz and Taylor of the basics of FEM [6]. Microstructure characterized by p(m) Domain W (macroscale) x q wall q T T wall ¶W ¶W m l = D D (m) Fig. 1. Thermal problem in a macroscopic domain Ω where the effective properties depend on a quantitatively characterized microstructure. Let the microstructure vary throughout Ω, being sampled at a series of points x(µ)∈Ω (µ= 1, . . . , N). Further, let the microstructure at any x(µ)be characterized by nparameters p(µ) i, grouped in the vector p(µ). Then, the effective conductivity kat x(µ)is at last a function of p(µ), i.e: k(x(µ)) = k(p(µ)). The microstructure throughout Ω is characterized by the vector P= [p(1),...,p(N)]. Then, the global conductivity matrix K(equation (6)) is a function of P, and so they are the nodal temperatures Tj (solution of equation (5)) and the temperature field T(approximated by equation (4) for FEM). 2 3 Design of the microstructure to control the heat flux To design the microstructure in Ω consists of finding Psuch that Ω has a given response. In this case, we aim to enforce the heat flux to take prescribed values ¯ q(q)at a series of points x(q),q= 1, . . . , Q. The heat flux at x(q)is given by q(x(q))=[−kgrad T]x(q)=−k(p(q)) [ grad Nj]x(q)Tj(P)≡q(q)(P) Then, we have to find Psuch that q(q)(P) = ¯ q(q),for q= 1, . . . , Q. (8) Let us look for Pwithin a space Pof admissible solutions. Generally, the task (8) cannot be exactly satisfied by any P∈ P. So, let us accomplish this task as well as possible by solving the nonlinear constrained optimization problem min P∈P 1 Q Q X q=1   q(q)(P)−¯ q(q)   2,(9) that is, by minimizing the root mean square error in the accomplishment of the task (8). 4 Application to thermal shield and cloak Let us apply the proposed methodology to design a device for thermal shielding and cloaking as alternative to that designed by Schittny et al. [3] based on transformation thermodynamics. This device, embedded in a plate with prescribed heat flux, is designed to thermally shield the region it encloses without altering the outer flux (this is thermal cloaking). The Ω domain is the entire plate, a rectangle of sides W= 18 cm and H= 12 cm subject to the following boundary conditions: T=Tmax = 80◦C for x=−W/2, T=Tmin = 25◦C for x=W/2, and q·n= 0 for y=±H/2 (see Figure 2). The heat flux normal to the plate is neglected. The plate is made of a homogeneous and isotropic material with thermal conductivity k0= 85 W/mK. Without the device, the heat flux in the plate is q0=k0(Tmax −Tmin)/W 0=26.0 kW/m2 0. The considered device is the ring Ωfree ⊂Ω with inner radius r= 2.5 cm and exterior R= 5 cm, see Figure 2. This ring is designed to thermally shield the region Ωfixint, made of copper with isotropic conductivity kcopper = 394 W/mK. A further design requirement on the device is to keep the heat flux outside it (i.e., in the remainder portion of the plate, Ωfixext) unaltered. The domain Ω is discretized using a mesh of 75 ×50 bilinear rectangular finite elements, as shown in Figure 3. Each blue element, belonging to the device Ωfree, has a microstructure sampling point. In the other elements, the material is prescribed. In the red elements, the heat flux is prescribed (as it will be explained in Section 4.2). 4.1 Metamaterial characterization The current problem of thermal shielding and cloaking is similar to that solved by Schittny et al. [3]. They used a composite consisting of a copper plate with holes filled by polydimethylsiloxane (PDMS), as shown in Figure 4. Both materials have highly different isotropic conductivities: kcooper = 394 W/mK and kPDMS = 0.15 W/mK. Varying the shape and density of the holes throughout the plate, they achieved the distribution of effective anisotropic conductivity allowing the desired heat flow. In this work, following [2], the anisotropy in the effective conductivity is achieved by using a stacked composite made of alternating sheets of materials A and B with different isotropic conductivities. Like Schittny et al. [3], we adopted A=copper and B=PDMS. As shown in Figure 3(b), the representative 3 Adiabatic fijint fixext libre Wfree Wfixint r R T=Tmax W Hx y T=Tmin Adiabatic Wfixint W Wfree Fig. 2. Domain of analysis for the thermal shielding and cloaking problem. volume element (RVE) of the microstructure of such composite at the sampling point x(µ)∈Ωfree is a unit square characterized by the vector p(µ)of components p(µ) 1=d(µ) A(thickness of sheet A) and p(µ) 2=θ(µ)(orientation of the sheets); the thickness of material B is dB=lµ−dA, where lµ= 1 is the thickness of the RVE. The effective thermal conductivities at x(µ)in the direction of the sheets (λ) and normal to the sheets (τ) are kλλ(x(µ)) = d(µ) AkA+d(µ) BkB lµ =kλλ(d(µ) A), kττ (x(µ)) = lµ d(µ) A/kA+d(µ) B/kB =kττ (d(µ) A), which are arranged in the matrix k0(x(µ)) = "kλλ(d(µ) A) 0 0kττ (d(µ) A)#=k0(d(µ) A). Now, the matrix of tensorial components of the effective conductivity referred to the fixed Cartesian frame x-yat the point x(µ)can be computed as k(x(µ)) = [R(θ(µ))]Tk0(d(µ) A)R(θ(µ)) = k(d(µ) A, θ(µ)) = k(p(µ)),(10) where Ris the rotation matrix R(θ) = cos θsin θ −sin θcos θ. The equation (10) explicitly defines the effective conductivity at a point as a function of the microstructure at that point. 4.2 Design of the microstructure in the shielding-cloaking device Let us design the device Ωfree such that it does its best to thermally shield Ωfixint without altering the heat flux in Ωfixext. This amounts to determine the metamaterial distribution in Ωfree by solving the 4 Representative volume element (RVE) W(1) W(2) W(3) Wfree q b) Material B Material A l t lm lm dB (m) dA (m) (m) a) Discretized domain =1 Fig. 3. (a) Finite element mesh of the analyzed domain Ω; the blue elements belong to the heat flux manipulating device Ωfree, and the red ones have prescribed heat flux. (b) Representative volume element (RVE) of the microstructure at an element in the device. rR Wfree Wfixint Wfixext Fig. 4. Drilled cooper plate from the experiment of thermal shielding and cloaking proposed by Schittny et al. [3] optimization problem min P 1 QX q  q(q)(P)−¯ q(q)   2,(11) where the sum is extended to the Q= 444 red elements in Figure 3, those lying in Ω(1) ⊂Ωfixext, Ω(2) ⊂Ωfixext, and Ωfixint;¯ q(q)=q0at the center of elements in Ω(1) and Ω(2) (this is the cloaking task), and ¯ q(q)=0at the center of elements in Ωfixint (this is the shield task). The vector of design variables Pis the set of vectors p(µ)characterizing the microstructure at the N= 1014 elements of Ωfree, with P2µ−1=p(µ) 1=d(µ) Aand P2µ=p(µ) 2=θ(µ),µ= 1,2,...N. For the chosen metamaterial, the following box constraints arise: 0≤P2µ−1=d(µ) A≤1, 0≤P2µ=θ(µ)≤π, with µ= 1, . . . , N. Here, this nonlinear constrained optimization problem was solved using the IPOPT interior point algorithm [7]. To avoid “checkerboard” instabilities at the optimal solution, we used the density filter technique proposed by Sigmund [8]. 5 4.3 Results The optimal solutions for dA(that is in fact the fraction of copper since the RVE was assumed to be a unit square), dB(the fraction of PDMS) and θ(the orientation of the sheets) in the device are plotted in Figure 5. The fraction of copper was found to increase radially in the whole ring. This fraction also increased radially in the solution of Schittny et al. [3] (see Figure 4), but in a concentric way. Figure 6 shows the effective conductivity field in the device Ωfree, generally anisotropic, as result of the optimal metamaterial distribution. For this conductivity distribution in the device, the temperature field in the whole plate is that shown in Figure 7. There, it can be seen that the cloaking task was closely accomplished (the isotherms are mostly vertical and evenly separated outside the device) and so it was the shield task (the temperature range inside Ωfixint is not larger than 1◦C). Quantitatively, the root mean square error in the accomplishment of both tasks is 69.5 W/m2= 2.7×10−3kq0k, which is considered to be highly satisfactory. dA 1 0.88 0.66 0.44 0.22 0.12 dB 0.88 0.66 0.44 0.22 0 q 2.68 2.02 1.37 0.71 0.05 Fig. 5. Variation of the fraction of copper (dA), the fraction of PDMS (dB), and the orientation of copper and PDMS sheets in the shielding-cloaking device. kxy 177.1 91.2 0 -91.2 -187.6 kxx 388.4 291.1 194.1 97.0 0.24 kyy 393.9 293.5 195.6 97.8 2.7 Fig. 6. Cartesian components of the effective thermal conductivity in the shielding-cloaking device, in W/mK. 5 Conclusions We presented a novel method for designing metamaterials to control the diffusive heat flux in ways that were inconceivable using ordinary materials. This method consists in solving an optimization problem where the objective function to be minimized is the error in the accomplishment of a given heat flux control task, and the design variables define the microstructure in a heat flux manipulating device. Its potentiality was proved by designing a highly efficient device for thermal shielding and cloaking. Such solution may be useful for protecting sensitive regions in an electrical circuit or in a chip from excessive heating. This and other real problems are the target of future works. 6 80.0 66.3 52.5 38.8 25.0 Wfixext Wfixint Wfree Fig. 7. Temperature distribution in the plate, in ◦C. The difference between isotherm is 1◦C. Acknowledgments We thank the following institutions for supporting this work: •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). •The National Scientific and Technical Research Council (CONICET) from Argentina, through the project PIP 1105 “Simulaci´on Computacional de Problemas en Multif´ısica. Aplicaci´on a Solidificaci´on de Metales y Dispositivos Microelectromec´anicos” [Computational Simulation of Multiphysics Problems. Application to Metal Solidification and Micromechanical Devices]. References [1] M. Maldovan. Sound and heat revolutions in phononics. Nature, 503:209–217, 2013. [2] S. Narayana and Y. Sato. Heat flux manipulation with engineered thermal materials. Phys. Rev. Lett., 108(21):214303, 2012. [3] 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. [4] J. B. Pendry, D. Schurig, and D. R. Smith. Controlling electromagnetic fields. Science, 312(5781): 1780–1782, 2006. [5] H. Chen and C. T. Chan. Transformation media that rotate electromagnetic fields. Appl. Phys. Lett., 90:241105, 2007. [6] O. C. Zienkiewicz and R. L. Taylor. The finite element method, volume 1: The basis. ButterworthHeinemann, 2000. [7] A. W¨achter and L. T. Biegler. On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming. Math. Program., Ser. A, 106:25–57, 2006. [8] O. Sigmund. Morphology-based black and white filters for topology optimization. Struct. Multidisc. Optim., 33:401–424, 2007. 7