Designing thermal energy harvesting devices with natural materials through optimized microstructures
Qingxiang Ji, Xueyan Chen, Jun Liang, Vincent Laude, Sébastien Guenneau, Guodong Fang, Muamer Kadic
DDesigning thermal energy harvesting devices with natural materialsthrough optimized microstructures
Qingxiang Ji a , b , Xueyan Chen a , b , Jun Liang c , Vincent Laude b , S´ebastien Guenneau d , Guodong Fang a , ∗ ,Muamer Kadic b a National Key Laboratory of Science and Technology on Advanced Composites in Special Environments, Harbin Institute ofTechnology, Harbin, 150001, China b Institute FEMTO-ST, CNRS, University Bourgogne Franche-Comt´e, 25000 Besan¸con, France c Institute of Advanced Structure Technology, Beijing Institute of Technology, Beijing 100081, China d UMI 2004 Abraham de Moivre-CNRS, Imperial College, London SW7 2AZ, UK
Abstract
Metamaterial thermal energy devices obtained from transformation optics have recently attracted wide atten-tion due to their vast potential in energy storage, thermal harvesting or heat manipulation. However, thesedevices usually require inhomogeneous and extreme material parameters which are di ffi cult to realize in large-scale applications. Here, we demonstrate a general process to design thermal harvesting devices with avail-able natural materials through optimized composite microstructures. We apply two-scale homogenizationtheory to obtain e ff ective properties of the microstructures. Optimal Latin hypercube technique, combinedwith a genetic algorithm, is then implemented on the microstructures to achieve optimized design parame-ters. The optimized microstructures can accurately approximate the behavior of transformed materials. Wedesign such devices and numerically characterize good thermal-energy harvesting performances. To validatethe wide-range application of our approach, we illustrate other types of microstructures that mimic well theconstitutive parameters. The approach we propose can be used to design novel thermal harvesting devicesavailable with existing technology, and can also act as a beneficial vehicle to explore other transformationoptics enabled designs. Keywords:
Transformation thermodynamics; Thermal energy harvesting; Microstructures; Optimization
1. Introduction
Transformation optics was first proposed to per-form cloaking on electromagnetic waves, based onform-invariance of the governing equations after co-ordinates transformation [1]. Since then the con-cept was promoted in various physical fields [2–10],creating miscellaneous devices, like invisible cloaks[11], carpets [12], invisible sensors [13], illusion de-vices [14, 15] or hyperlenses [16–18].Over the past decades, the technique of thermalenergy harvesting and management, which can col- ∗ Corresponding authorE-mail address: [email protected] lect and store heat energy from the ambient sur-roundings, has attracted renewed interest [19, 20].Transformation thermodynamics, a counterpart oftransformation-optics, has been proposed to guideheat flux in thermal management and generate novelthermal meta-devices, such as heat cloaks, thermalenergy harvesting devices, thermal sensors, etc [21–26]. Thermal energy harvesting devices, which canfocus and harvest ambient thermal energy withoutsevere perturbations to the heat profile outside thedevices, have vast potential in improving the energy-conversion e ffi ciency of existing technologies [27].A challenge for thermal meta-devices is that theyoften require inhomogeneous and anisotropic consti- Preprint submitted to Elsevier August 21, 2020 a r X i v : . [ phy s i c s . a pp - ph ] J u l utive parameters which are di ffi cult to realize es-pecially for large-scale applications. As the het-erogenous constitutive profile is position-dependentand continuous, some form of discretization is re-quired. Moreover, anisotropic materials can be ap-proximated by structures of thin and alternating lay-ers based on e ff ective medium theory [8, 28]. Ther-mal devices following this principle were fabricatedand experimentally characterized [23, 29]. In theseworks, step-wise approximation of the ideal pa-rameters were made, which sacrificed the perfor-mances to reduce fabrication di ffi culty. Other re-searchers realized thermal cloaks with bulk isotropicmaterials using the scattering-cancellation approach[30–33], which was then extended to the design ofthermal concentrating devices. Following this ap-proach, a class of solar-shaped thermal harvestingdevices were demonstrated, which could manipulateand concentrate the heat flux using natural materialswithout singularities [34, 35]. However, these meta-devices are always shaped as regular profiles (cylin-der, sphere or ellipse), as the scattering-cancellationmethod is non-trivial for irregular or complex shapes.Except for the mentioned strategies, do we havemethods that can both o ff er good performance andfabrication simplicity? The answer is positive. Re-searchers [36, 37] demonstrated that the constitu-tive medium can be approximated by microstructuresin the context of acoustic metamaterials. Thanksto similarities between the governing equations foracoustic and thermal fields, a similar approach canbe implemented for thermal metamaterials. Then thekey problem is to determine the e ff ective property ofthe built microstructure and find the proper param-eters of the microstructure to mimic desired trans-formed parameters. Ji et al. [38] achieved thermalconcentration using fiber reinforced microstructuresbased on simplified e ff ective medium theory. Pomot et al. [39] realized acoustic cloaking by microstruc-tures with several types of perforations, combinedwith a genetic algorithm. We note that the partialdi ff erential equations (PDEs) therein have the samestructure to the heat equation in the static limit, inwhich case one simply has to replace the density bythe conductivity and all the asymptotic analysis car-ries through. Moreover, PDEs in Ref. [39] are sup-plied with homogeneous Neumann boundary condi- tions that stand for rigid walls in acoustics, and insu-lating walls in thermodynamics. It is thus temptingto implement homogenization and e ff ective mediumtheories similar to those already used in acousticmetamaterials, however bearing in mind that in thedynamic regime the acoustic wave (elliptic) PDE andthe heat di ff usion (parabolic) PDE are very di ff erentin nature.In our work we focus on the thermal field andestablish a general design road map to obtain re-alizable thermal harvesting devices utilizing micro-structures. We apply two-scale homogenizationtheory to determine the equivalent thermal proper-ties and employ the Optimal Latin-Hypercube tech-nique to obtain desired design parameters. Be-yond this problem, the process is applicable to othermicrostructures and to wave problems (such as inacoustics). We stress that because of their complex-ity, the microstructures we investigate here would bechallenging to achieve otherwise. Finally we numer-ically characterize such thermal harvesting deviceswith natural materials and verify their harvesting ef-ficiency by finite element simulations.
2. Methods
We recall the heat conduction equation withoutheat sources ∇ · ( k ∇ T ) − ρ c ∂ T ∂ t = , (1)where ∇ is the gradient, T is the temperature, ρ c isthe product of density by heat capacity and k denotesthe heat conductivity. Following the theory of trans-formation optics [1] and thermodynamics [10], thegoverning equation will remain unchanged under acoordinate transformation when the transformed pa-rameters satisfy k (cid:48) = JkJ T det( J ) and ( ρ c) (cid:48) = ρ cdet(J) , (2)where J is the Jacobian of the geometric transfor-mation, J T its transpose and det( J ) its determinant.We first consider steady-state case where parameters ρ c vanish. The transformed conductivity k (cid:48) derived2y Eq. (2) is usually anisotropic and space depen-dent, which is di ffi cult to achieve in practice. To re-move the singularities and make the proposed devicesimpler to realize, we choose the following nonlineartransformation that maps domain Ω ( r ) onto domain Ω (cid:48) ( r (cid:48) ): r (cid:48) = r (cid:32) rr (cid:33) C , (3)which yields constant relative radial ( k (cid:48) r ) and tangen-tial heat conductivity ( k (cid:48) θ ) in cylindrical coordinatesas k (cid:48) r = C and k (cid:48) θ = , (4)where C is a constant with C >
1. Here, r and r arethe inner radius and the outer radius of the designedshell region, respectively (see details in the supple-mental material).In stark contrast to the material parameters ob-tained from rigorous transformation optics, such adevice is homogeneous in materials composition andits performance is only determined by the thermal-conduction anisotropy (characterized by constant C ;see details in the supplemental material). We em-ploy this geometric transformation as it avoids theneed for extreme spatially-varying parameters andthus the transformed medium will be much easier toimplement in practice. Moreover, it is enough forthe proof-of-concept demonstration. We stress thatthe method we describe is indeed also applicable tothe design of rigorous transformation-based thermaldevices.For clarity, we outline the design process for thegeneral case with a specific example. Assume thatwe create a cylindrically symmetric thermal harvest-ing device with inner radius r = .
01 m and outerradius r = .
04 m (see Fig. 1). The proposed de-vice can harvest thermal energy from the surround-ings and concentrate it into the inner domain r (cid:48) < r .Heat energy density in the inner domain is thus sig-nificantly increased. In our design, the thermal con-ductivity of the background is k b =
132 Wm − K − .We set C =
2, which implies that for the shell regionwe have k r =
264 Wm − K − , k θ =
66 Wm − K − .Now we turn to the realization of shell region pa-rameters by two natural materials A and B throughoptimized microstructures. Fig. 1.
Schematic representation of a possible realization of thethermal harvesting device. The shell displays some periodicityalong the azimuthal direction, which is a hallmark of a concen-trator. Such a device can harvest thermal energy from its sur-roundings. Heat flows are concentrated into the inner domain r (cid:48) < r thanks to the designed shell region r ≤ r (cid:48) ≤ r . The material distributions in Eq. (4) need to bemapped onto a microstructure exhibiting prescribedconstitutive parameters. We build a medium withidentical elementary cells repeating periodically inspace (see Fig. 2 ). Generally, the well-establishede ff ective medium theory plays a dominant role indetermining e ff ective properties and is easy and di-rect for simple geometries. When designing meta-devices with complex-shaped structures, however, itwould be far from trivial to evaluate the equivalentproperties. Here we apply instead two-scale homog-enization theory [40] to determine e ff ective parame-ters with an asymptotic approach.We consider a two-dimensional periodic mediumwith square elementary cells [0 , η ] of side-length η (cid:28)
1. The solution T η of the steady-state heat equa-tion with fast oscillating parameters A η = k ( x /η, y /η ) ∇ · (cid:16) k η ∇ T η (cid:17) = η tends to zero, tothe solution T hom of the homogenized heat equation ∇ · ([ k hom ] ∇ T hom ) = . (6)3 ig. 2. Schematic illustration of 2D periodic lattices definedby 2 parameters (area fractions f and f ). We choose this mi-crostructure to mimic the desired medium of the shell regionthrough homogenization and optimization. The e ff ective property of the periodic medium isgiven by[ k hom ] = (cid:32) (cid:104) k (cid:105) − (cid:104) k ∂ x V (cid:105) −(cid:104) k ∂ x V (cid:105)−(cid:104) k ∂ y V (cid:105) (cid:104) k (cid:105) − (cid:104) k ∂ y V (cid:105) (cid:33) , (7)where ∂ x : = ∂/∂ x , ∂ y : = ∂/∂ y and < . > denotesthe mean operator over the periodic cell. V ( x , y ) and V ( x , y ) are solutions defined up to an additive con-stant of auxiliary problems of thermostatic type onthe periodic cell [40]: (cid:40) ∇ · (cid:2) k ( x , y ) ∇ ( V − x ) (cid:3) = ∇ · (cid:2) k ( x , y ) ∇ ( V − y ) (cid:3) = . (8)We solve the auxiliary problems in weak form byusing COMSOL Multiphysics, which sets up the fi-nite element problem with periodic conditions im-posed to the field on opposite ends of the elemen-tary cell. We note in passing V and V are uniquesolutions of Eq. (8) up to additive constants, butthese constants do not a ff ect the homogenized con-ductivity, as one can see that in Eq. (7) only thepartial derivatives of V and V are involved. Thepotentials V ( x , y ) and V ( x , y ) of an illustrative case( f = . , f = . ff ective conductivity as[ k hom ] = (cid:32) .
01 3 . − − . − . (cid:33) = (cid:32) k k (cid:33) . (9)It is noticed that in Eq. (9) the o ff -diagonal com-ponents are almost negligible and originate fromnumerical errors. The resulting spurious artificial Fig. 3.
Potentials V ( x , y ) and V ( x , y ) of the illustrative casewhere f = . , f = . anisotropy can be safely ignored. We emphasizethat with the finite element method, the e ff ective ten-sor can be obtained for any periodic composite andthat the same technique was implemented before inacoustics [39] and electromagnetism [40], in whichcase the e ff ective density and permittivity tensors canbe deduced from the same annex problems as in ourthermal case.The e ff ective properties of the medium can betuned by several design parameters such as geome-try (here we use area fractions f and f ) and ma-terial properties (thermal conductivities k A and k B ).We restrict our attention here to predefined mate-rial properties and set the geometrical parameters asvariables. Our goal is to find the set of geometri-cal parameters which properly mimic the homoge-nized medium. Therefore, the geometry of elemen-tary cells should be tuned to obtain desired equiv-alent properties. Considering the heavy workloadof the trial-and-error method, we implement an Op-timal Latin hypercube technique to solve the prob-lem. We note here that the two-scale homogenizationtechnique has been already used to design a thermalconcentrator similar to the one shown in Fig. 1 andto show that a concentrator consisting of concentriclayers would require some complex valued conduc-tivities with sign-shifting imaginary parts [42]. In thepresent case, we investigate doubly periodic designsand thus we do not face such pitfalls. Optimal Latin Hypercube Sampling (OLHS) tech-nique is applied to optimize the spatial positions of4ontrol points. The aim of this process is to designa matrix where the sample points spread as evenlyas possible within the design region. This methodis e ffi cient and robust due to its enhanced stochas-tic evolutionary algorithm and significant reductionin matrix calculations to evaluate new / modified de-signs during searching [43]. In this work, we pre-define the two materials as air ( k A = . k B = − K − and select thegeometrical (area fraction) parameters f and f asdesign variables. The ranges of these variables aredefined as 0 . < f < . . < f < . ff ective thermal conductivity of the peri-odic medium for each sample point using the two-scale homogenization theory. Details of the gener-ated sample points and calculated results are listed inTable S1. We emphasize the importance of choos-ing proper ranges of the variables. If, for instance,we use much larger ranges for f and f , the derivedvalues can be less accurate unless a more refined dis-cretization (more sample points) is chosen.We then create an approximation surrogate modelfrom the obtained one hundred space samples [44].The surrogate model is built by the Elliptical Ba-sis Function Neural Network technique which estab-lishes a relation between design targets ( k and k )and variables ( f and f ), as shown in Fig. 4.We use two estimators to evaluate the reliability ofthe surrogate model, the coe ffi cient of determination( R ) and the root mean square error (RMSE). Theseare defined as R = − (cid:80) mi = ( y i − ˆ y i ) (cid:80) mi = ( y i − ¯ y i ) , (10)RMSE = (cid:114) (cid:80) mi = ( y i − ˆ y i ) m , (11)where y i and ˆ y i are respectively the real value andthe predicted value of the objective function over thesame sample points, ¯ y i is the mean value of all objec-tive functions and m is the total number of samplingpoints. The closer R is to 1 and RMSE is to 0, themore accurate the model. In Fig. 4 (c-d), we can ob-serve that the predicted values are in good agreementwith the actual values. The calculated R and RMSE are listed in Table1. The average error and the maximum error amongall samples are also shown. In the constructed sur-rogate models, R is larger than 0.99997 and RMSEare smaller than 0.00012. The maximum error re-mains close to 0, which demonstrates that the surro-gate models are accurate. Table 1
Accuracy of the constructed surrogate model.
Error type k k RMSE 1 . − .
88 10 − Average 6 .
45 10 − .
52 10 − Maximum 5 .
47 10 − .
88 10 − R . E ( k , k , α, β ) = α | k − k r | + β | k − k θ | , (12)where α and β are weighting factors for the two ob-jective sub-functions such that α + β =
1. Thefunction E ( k , k , α, β ) measures the overall di ff er-ence between obtained and objective values. We de-fine α = β = .
5, i.e. equal weights for the di-agonal tensor elements k and k . The principle isto obtain a global minimum within the discrete so-lution space. We solve this inverse problem by us-ing a Non-Dominated Sorting Genetic Algorithm ap-proach [46]. The one hundred random structures,corresponding to the parameter space samples ob-tained by the OLHS method, form the first genera-tion. During the search process, the population sizeand the number of generations are defined as 12 and200, respectively. New generations are created us-ing crossover and mutation processes. We set themutation distribution index and crossover distribu-tion index as 20 and 10, respectively. The crossoverprobability is set as 0.9. For the sake of clarity, wedo not detail the well-established genetic algorithmapproach. In short, the Non-Dominated Sorting Ge-netic Algorithm performs well enough in our caseand enables us to determine e ffi ciently the desireddesign parameters.5 ig. 4. Panels (a) and (b) show the relation between design targets (heat conductivities k and k ) and design variables (geometricalparameters f and f ). The results are obtained basing on the chosen sample points and using two-scale homogenization theory.Panels (c) and (d) show good agreement between actual values and predicted values of the heat conductivities, demonstrating thereliability of the surrogate model. The indices k and k denote heat conductivity of the elementary cell in the x and y directions,respectively. Following this approach, we obtain the desired setof parameters as f = . f = . C =
2. We then implement the microstructure andcalculate the corresponding heat conductivities. Theresult are listed in Table.2. It can be seen that theobtained parameters closely mimic the desired trans-formed medium.
3. Scheme validation and discussion
We now turn to the theoretical recipe for realiz-ing thermal harvesting based on optimized compos-
Table 2
Comparison of predicted value and targeted value for the de-rived set of parameters.
Predicted Targeted Relativedi ff erence k (Wm − K − ) 263.97 264 0.01% k (Wm − K − ) 66.01 66 0.02%ite microstructures. Note that the heat conductivityin Eq. (4) is expressed in ( r , θ ) polar coordinateswhile the microstructure is designed in ( x , y ) Carte-sian coordinates. To build a cylindrical thermal har-6esting device, we discretize the shell region intonumerous units and transplant microstructures withmatched geometrical parameters into each unit. Wedesign in this paper a discrete thermal harvesting de-vice with 15 radial layers and 45 tangential sectors.Similarly to what has been done in Ref. [39], wecould have considered an increasing number of lay-ers and sectors to very accurately approximate theidealized thermal concentrator parameters. However,we focus here on a practically implementable design.The materials constituting the microstructure are nat-urally available materials: copper (material A) andair (material B). It is understood that the device couldbe built from other materials. We use copper and airhere considering their applicability for realistic ex-periments.Numerical calculations were conducted, wheretemperatures at the left and right boundary are re-spectively imposed as 1 K and 0 K, for easiness innormalization. We apply Neumann (perfect insula-tor) conditions at other boundaries. As indicated inthe scheme (Fig. 5), iso-thermal lines are signif-icantly compressed to the inner domain ( r ≤ r ).Hence, heat flux density in the inner domain is en-larged, implying that more heat energy is concen-trated into the central region. In addition, iso-thermallines in the background are uniform with little per-turbations. That is, the heat energy is harvested andconcentrated into the central region without muchperturbation of the external thermal field. We empha-size that similar computations would also hold fortime-harmonic acoustic and electromagnetic equa-tions.We further conduct a quantitative analysis of thethermal harvesting behavior of the proposed scheme.Two measurement lines are defined. An horizon-tal line ( y =
0) is selected to reveal perturbationsof the external thermal field whereas a vertical line( x = r ) is chosen to illustrate the heat harvesting ef-ficiency. We also build an additional contrast platethat occupies the same area as the harvesting devicebut that is composed only of a homogeneous back-ground medium. The following index η is defined toevaluate the energy harvesting e ffi ciency as η = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) T x = r − T x = − r T x = r − T x = − r (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . (13) Fig. 5.
Temperature fields in the steady state for (a) the idealcase and (b) the proposed thermal harvesting device. In bothcases, heat flows are concentrated to the inner domain withoutmuch perturbations to the external thermal field.
An index M V is introduced to characterize the pertur-bations to external fields ( i.e. the thermal neutralityof the concentrator) M V = (cid:82) Ω | T ( x , y , z ) − T r ( x , y , z ) | d Ω (cid:82) Ω d Ω , (14)where Ω denotes the probe domain of external ther-mal fields and T r represents the temperature field ofthe homogeneous medium. The index M V reveals all7 ig. 6. Comparison of temperature profiles between the harvesting device and the contrast plate along the measurement lines (a) y = x = r . All results are obtained in the steady-state and the insets denote positions of measurement lines. We define ∆ T = | T ( x , y , z ) − T r ( x , y , z ) | to quantitatively demonstrate the heat harvesting performances. Significant heat-focusing in the innerdomain ( − r < x < r ) and minor fluctuations in external domain (along x = r ) are both oberved, demonstrating good heatconcentrating performances. perturbations to the external heat profile [47].It is apparent in Fig. 6(a) that thermal energy issignificantly concentrated into the inner domain, asa tight focusing with a local heat-intensity increaseis observed. We list thermal gradients inside the in-ner domain (characterized by ∆ T in = (cid:12)(cid:12)(cid:12) T x = r − T x = − r (cid:12)(cid:12)(cid:12) )and harvesting e ffi ciencies for di ff erent cases in Tab.3. The results indicate that thermal gradients in thecentral region are almost twice as large as in thecontrast plate. The concentration e ffi ciency is sig-nificantly lifted. Both results reach nearly theoreti-cal values, demonstrating that the heat concentratingscheme is both e ff ective and accurate. Table 3
Thermal harvesting performances for di ff erent cases. Index Ideal Proposed Bare plate ∆ T in (K / m) 0.4 0.39 0.2 η M V (K) 0 0.0025 0.0003It can be observed in Fig. 6(b) that temperaturesalong the vertical measurement line are almost uni-form, indicating that the external thermal field is onlyslightly influenced. We notice some minor perturba-tions in Tab. 3 which are mainly due to the discretiza-tion process and to numerical errors. Thus far, all that has been achieved for heat canbe directly translated to airborne acoustics (withrigid inclusions) and electromagnetism (with per-fectly conducting inclusions, in the case when themagnetic field is polarized perpendicular to the xy-plane), since governing equations are identical in thestatic limit, and so results in Fig. 6 and 7 hold forthe acoustic and electromagnetic counterparts of thethermal concentrator. We would like to investigatenow the di ff usive nature of heat conduction. We de-fine same boundary conditions as the steady statecase and focus here on the evolution of thermal har-vesting performance over time. Temperature distri-butions at di ff erent time steps t are shown in Fig. (7),where thermal gradients and perturbations are alsoplotted as a function of time. It is observed that per-turbations of the external thermal profile increase andthen decrease after a certain lapse of time, whereasthe thermal gradient of the object increases graduallytoward its maximum value. Generally, the proposeddesign performs well for harvesting thermal energy,with significant heat flux concentrated in the innerdomain and little perturbation to the outside thermalfield, once the permanent regime has been reached.We stress that the design process proves feasible asthe created device works e ffi ciently and convergesto almost the same harvesting performance as in the8deal case. Fig. 7.
Thermal harvesting performances of the concentratorin the transient regime. (a) Temperature fields are illustrated atdi ff erent times t . The performances get gradually better after acertain lapse of time. (b) Thermal gradients of the inner domain( ∆ T in = (cid:12)(cid:12)(cid:12) T x = r − T x = − r (cid:12)(cid:12)(cid:12) ) and perturbations to external thermalfield ( M V ) are plotted as functions of time. The quantitativeresults agree well with those revealed in (a). The heat concentration e ffi ciency is directly deter-mined by the constant C [35], which in turn can beapproached by the material properties and the geom-etry of the microstructures. To gain more insightsinto the underlying mechanism of the microstructure,we derive and show in Fig. 8 the relation between de-sign variables (area fraction f and f ) and the heatconcentration e ffi ciency. For comparison, we furtherassume two other materials C and D to substitute for material B, with heat conductivity k C < k B < k D . Itis observed that higher concentration e ffi ciency re-quires larger f and smaller f which directly in-dicates larger geometrical anisotropy. Besides, wenotice that smaller design variables are needed toachieve a given e ffi ciency for a larger conductivity ofthe second material (for constant material A). Largermaterial anisotropy also allows for a larger maximumconcentration e ffi ciency. Fig. 8.
Relation between design variables (area fractions f and f ) and the concentration e ffi ciency. The inset shows thestudied elementary cell where one material is pre-defined as k A = k B =
400 Wm − K − ), we also considered material C( k C =
300 Wm − K − ) and material D ( k D =
600 Wm − K − ) fora comparison. Note that the microstructure we employed in Fig.2 is not the unique choice. One can use other elemen-tary cells as long as they are appropriate to obtain therequired anisotropy. To validate the wide-range ap-plication of our approach, we further present severalother elementary cells and build corresponding de-vices. Steady state simulations are conducted withsame boundaries with the aforementioned case. InFig. 9(a), we set material properties as design vari-ables and recover the widely used solar-shaped de-vice. In Fig. 9(b), we obtain the same design asin Fig. 1 but based on a di ff erent elementary cellobtained by a sub-lattice translation. Besides, weshow in Fig. 9(c) a split ring element cell for which9t seems unlikely that the e ff ective medium theorycould be easily applied. Good harvesting perfor-mance is again achieved without much perturbationsto the external field. As a note, the symmetry of thesplit ring element cell is reduced compared to theother cases, but the e ff ective tensor remains diago-nal. The strength of our approach becomes apparentwhen dealing with such more complex geometries,where the application of e ff ective medium theory isfar from trivial.We stress that using the optimization approach inthis paper high-performance thermal-harvesting de-signs can be obtained that meet given external con-straints, such as maximum or minimum values formaterial properties and geometrical parameters. Wecan search for the best solution among materialsproperties, microstructures and geometrical param-eters within the full available range, providing a sub-stantial flexibility and many degrees of freedom inpractical applications.We considered in this paper two-dimensional (2D)harvesting devices as fabrication and measurementare presumably simpler than in 3D. The proposed op-timization method, however, is also well suited to thedesign of 3D thermal harvesting devices. A majordi ff erence is that more design variables are requiredin 3D, and thus more computational resources wouldbe required to undertake such a study, for which thetheoretical part of the present work would apply mu-tatis mutandis.
4. Conclusion
We proposed a general method to design thermal-energy harvesting devices from naturally availablematerials. We designed composite microstructuresand calculated their e ff ective conductivity by thetwo-scale homogenization technique. We then im-plemented the Optimal Latin hypercube method toobtain design parameters that can best mimic thetransformed medium. We built a harvesting devicemodel based on the obtained optimized microstruc-tures and demonstrated good thermal harvesting per-formance, thereby validating the e ff ectiveness of ourdesign method.The optimization method adds great flexibility tothe constraints that can be imposed on devices, which Fig. 9.
Additional optimized microstructures and correspond-ing thermal harvesting performance. (a) Laminar structures de-fined by two material parameters ( k and k ) if L = L and twoadditional geometrical parameters if L (cid:44) L . The optimizedparameters are k = . k = . L = L . (b) Elementary cell defined by two geometrical pa-rameters f = . , f = . f = . f = . f = . f = . f = . ff erent microstructures and obtain desired ther-mal harvesting performances. allows one to find the best possible solutions withdi ff erent geometrical structures and component ma-terials. We stress that the flexibility and simplicity ofthe method is a good addition to existing heat manip-10lation techniques, including other thermal function-alities, i.e. cloaking and illusion. It also paves a pathfor novel optical, acoustic and electromagnetic de-vices based on form-invariant governing equations. Acknowledgments
This work was supported by the EIPHI Gradu-ate School [grant number ANR-17-EURE-0002]; theFrench Investissements d’Avenir program, projectISITEBFC [grant number ANR-15-IDEX-03]; andthe National Natural Science Foundation of China[grant numbers 11732002 and 11672089].
References [1] J. B. Pendry, D. Schurig, D. R. Smith, Controlling electro-magnetic fields, Science 312 (5781) (2006) 1780–1782.[2] M. Kadic, T. B¨uckmann, R. Schittny, M. Wegener, Meta-materials beyond electromagnetism, Reports on Progressin Physics 76 (12) (2013) 126501.[3] G. W. Milton, M. Briane, J. R. Willis, On cloaking forelasticity and physical equations with a transformation in-variant form, New Journal of Physics 8 (10) (2006) 248.[4] T. B¨uckmann, M. Thiel, M. Kadic, R. Schittny, M. We-gener, An elasto-mechanical unfeelability cloak madeof pentamode metamaterials, Nature Communications 5(2014) 4130.[5] Y. Achaoui, A. Diatta, M. Kadic, S. Guenneau, Cloakingin-plane elastic waves with swiss rolls, Materials 13 (2)(2020) 449.[6] M. Kadic, M. Wegener, A. Nicolet, F. Zolla, S. Guenneau,A. Diatta, Elastodynamic behavior of mechanical cloaksdesigned by direct lattice transformations, Wave Motion92 (2020) 102419.[7] R. Schittny, M. Kadic, T. B¨uckmann, M. Wegener, In-visibility cloaking in a di ff usive light scattering medium,Science 345 (6195) (2014) 427–429.[8] D. Schurig, J. Mock, B. Justice, S. A. Cummer, J. B.Pendry, A. Starr, D. R. Smith, Metamaterial electromag-netic cloak at microwave frequencies, Science 314 (5801)(2006) 977–980.[9] S. A. Cummer, D. Schurig, One path to acoustic cloaking,New Journal of Physics 9 (3) (2007) 45.[10] S. Guenneau, C. Amra, D. Veynante, Transformation ther-modynamics: cloaking and concentrating heat flux, Op-tics Express 20 (7) (2012) 8207–8218.[11] M. Farhat, S. Guenneau, S. Enoch, Ultrabroadband elasticcloaking in thin plates, Physical Review Letters 103 (2)(2009) 024301.[12] J. Li, J. B. Pendry, Hiding under the carpet: a new strat-egy for cloaking, Physical Review Letters 101 (20) (2008)203901. [13] T. Yang, X. Bai, D. Gao, L. Wu, B. Li, J. T. Thong,C.-W. Qiu, Invisible sensors: Simultaneous sensing andcamouflaging in multiphysical fields, Advanced Materi-als 27 (47) (2015) 7752–7758.[14] M. Liu, Z. Lei Mei, X. Ma, T. J. Cui, Dc illusion and itsexperimental verification, Applied Physics Letters 101 (5)(2012) 051905.[15] R. Hu, S. Zhou, Y. Li, D.-Y. Lei, X. Luo, C.-W. Qiu,Illusion thermotics, Advanced Materials 30 (22) (2018)1707237.[16] Z. Liu, H. Lee, Y. Xiong, C. Sun, X. Zhang, Far-field opti-cal hyperlens magnifying sub-di ff raction-limited objects,Science 315 (5819) (2007) 1686–1686.[17] J. Li, L. Fok, X. Yin, G. Bartal, X. Zhang, Experimentaldemonstration of an acoustic magnifying hyperlens, Na-ture Materials 8 (12) (2009) 931–934.[18] M. Kadic, S. Guenneau, S. Enoch, S. A. Ramakrishna,Plasmonic space folding: Focusing surface plasmons vianegative refraction in complementary media, ACS nano5 (9) (2011) 6819–6825.[19] J. Freeman, I. Guarracino, S. A. Kalogirou, C. N.Markides, A small-scale solar organic rankine cycle com-bined heat and power system with integrated thermal en-ergy storage, Applied thermal engineering 127 (2017)1543–1554.[20] M. Amin, N. Putra, E. A. Kosasih, E. Prawiro, R. A. Lu-anto, T. Mahlia, Thermal properties of beeswax / graphenephase change material as energy storage for buildingapplications, Applied Thermal Engineering 112 (2017)273–280.[21] C. Fan, Y. Gao, J. Huang, Shaped graded materialswith an apparent negative thermal conductivity, AppliedPhysics Letters 92 (25) (2008) 251907.[22] L. Zeng, R. Song, Experimental observation of heat trans-parency, Applied Physics Letters 104 (20) (2014) 201905.[23] R. Schittny, M. Kadic, S. Guenneau, M. Wegener, Ex-periments on transformation thermodynamics: moldingthe flow of heat, Physical Review Letters 110 (19) (2013)195901.[24] T. Han, X. Bai, J. T. Thong, B. Li, C.-W. Qiu, Full con-trol and manipulation of heat signatures: Cloaking, cam-ouflage and thermal metamaterials, Advanced Materials26 (11) (2014) 1731–1734.[25] X. Peng, R. Hu, Three-dimensional illusion thermoticswith separated thermal illusions, ES Energy & Environ-ment 6 (2019) 39–44.[26] S. Huang, J. Zhang, M. Wang, R. Hu, X. Luo, Macroscalethermal diode-like black box with high transient rectifica-tion ratio, ES Energy & Environment 6 (2019) 51–6.[27] G. Xu, H. Zhang, Y. Jin, Achieving arbitrarily polygonalthermal harvesting devices with homogeneous parametersthrough linear mapping function, Energy Conversion andManagement 165 (2018) 253–262.[28] H. Chen, C. T. Chan, P. Sheng, Transformation optics andmetamaterials, Nature Materials 9 (5) (2010) 387.[29] S. Narayana, Y. Sato, Heat flux manipulation with ngineered thermal materials, Physical Review Letters108 (21) (2012) 214303.[30] T. Han, X. Bai, D. Gao, J. T. Thong, B. Li, C.-W. Qiu,Experimental demonstration of a bilayer thermal cloak,Physical Review Letters 112 (5) (2014) 054302.[31] X. Zhang, X. He, L. Wu, Ellipsoidal bifunctional thermal-electric transparent device, Composite Structures 234(2020) 111717.[32] H. Xu, X. Shi, F. Gao, H. Sun, B. Zhang, Ultrathinthree-dimensional thermal cloak, Physical Review Letters112 (5) (2014) 054301.[33] T. Han, P. Yang, Y. Li, D. Lei, B. Li, K. Hippal-gaonkar, C.-W. Qiu, Full-parameter omnidirectional ther-mal metadevices of anisotropic geometry, Advanced Ma-terials 30 (49) (2018) 1804019.[34] F. Chen, D. Y. Lei, Experimental realization of extremeheat flux concentration with easy-to-make thermal meta-materials, Scientific Reports 5 (2015) 11552.[35] T. Han, J. Zhao, T. Yuan, D. Y. Lei, B. Li, C.-W. Qiu, The-oretical realization of an ultra-e ffi cient thermal-energyharvesting cell made of natural materials, Energy & Envi-ronmental Science 6 (12) (2013) 3537–3541.[36] J. Mei, Z. Liu, W. Wen, P. Sheng, E ff ective dynamic massdensity of composites, Physical Review B 76 (13) (2007)134205.[37] D. Torrent, J. S´anchez-Dehesa, Acoustic cloaking in twodimensions: a feasible approach, New Journal of Physics10 (6) (2008) 063015.[38] Q. Ji, G. Fang, J. Liang, Achieving thermal concentra-tion based on fiber reinforced composite microstructuresdesign, Journal of Physics D: Applied Physics 51 (31)(2018) 315304.[39] L. Pomot, C. Payan, M. Remillieux, S. Guenneau, Acous-tic cloaking: Geometric transform, homogenization and agenetic algorithm, Wave Motion 92 (2020) 102413.[40] F. Zolla, S. Guenneau, Artificial ferro-magneticanisotropy: homogenization of 3d finite photoniccrystals, in: IUTAM Symposium on Asymptotics, Singu-larities and Homogenisation in Problems of Mechanics,Springer, 2003, pp. 375–384.[41] G. Allaire, Homogenization and two-scale convergence,SIAM Journal on Mathematical Analysis 23 (6) (1992)1482–1518.[42] D. Petiteau, S. Guenneau, M. Bellieud, M. Zerrad,C. Amra, Thermal concentrator homogenized with solar-shaped mantle, arXiv preprint arXiv:1508.05081 (2015).[43] J.-S. Park, Optimal latin-hypercube designs for computerexperiments, Journal of Statistical Planning and Inference39 (1) (1994) 95–111.[44] M.-W. Mak, S.-Y. Kung, Estimation of elliptical basisfunction parameters by the em algorithm with applicationto speaker verification, IEEE Transactions on Neural Net-works 11 (4) (2000) 961–969.[45] E. Cherkaev, Inverse homogenization for evaluation of ef-fective properties of a mixture, Inverse Problems 17 (4)(2001) 1203. [46] K. Deb, A. Pratap, S. Agarwal, T. Meyarivan, A fast andelitist multiobjective genetic algorithm: Nsga-ii, IEEETransactions on Evolutionary Computation 6 (2) (2002)182–197.[47] Q. Ji, X. Chen, G. Fang, J. Liang, X. Yan, V. Laude,M. Kadic, Thermal cloaking of complex objects with theneutral inclusion and the coordinate transformation meth-ods, AIP Advances 9 (4) (2019) 045029.ective dynamic massdensity of composites, Physical Review B 76 (13) (2007)134205.[37] D. Torrent, J. S´anchez-Dehesa, Acoustic cloaking in twodimensions: a feasible approach, New Journal of Physics10 (6) (2008) 063015.[38] Q. Ji, G. Fang, J. Liang, Achieving thermal concentra-tion based on fiber reinforced composite microstructuresdesign, Journal of Physics D: Applied Physics 51 (31)(2018) 315304.[39] L. Pomot, C. Payan, M. Remillieux, S. Guenneau, Acous-tic cloaking: Geometric transform, homogenization and agenetic algorithm, Wave Motion 92 (2020) 102413.[40] F. Zolla, S. Guenneau, Artificial ferro-magneticanisotropy: homogenization of 3d finite photoniccrystals, in: IUTAM Symposium on Asymptotics, Singu-larities and Homogenisation in Problems of Mechanics,Springer, 2003, pp. 375–384.[41] G. Allaire, Homogenization and two-scale convergence,SIAM Journal on Mathematical Analysis 23 (6) (1992)1482–1518.[42] D. Petiteau, S. Guenneau, M. Bellieud, M. Zerrad,C. Amra, Thermal concentrator homogenized with solar-shaped mantle, arXiv preprint arXiv:1508.05081 (2015).[43] J.-S. Park, Optimal latin-hypercube designs for computerexperiments, Journal of Statistical Planning and Inference39 (1) (1994) 95–111.[44] M.-W. Mak, S.-Y. Kung, Estimation of elliptical basisfunction parameters by the em algorithm with applicationto speaker verification, IEEE Transactions on Neural Net-works 11 (4) (2000) 961–969.[45] E. Cherkaev, Inverse homogenization for evaluation of ef-fective properties of a mixture, Inverse Problems 17 (4)(2001) 1203. [46] K. Deb, A. Pratap, S. Agarwal, T. Meyarivan, A fast andelitist multiobjective genetic algorithm: Nsga-ii, IEEETransactions on Evolutionary Computation 6 (2) (2002)182–197.[47] Q. Ji, X. Chen, G. Fang, J. Liang, X. Yan, V. Laude,M. Kadic, Thermal cloaking of complex objects with theneutral inclusion and the coordinate transformation meth-ods, AIP Advances 9 (4) (2019) 045029.