Materials issues in thermal modeling of thin film electrocaloric solid-state refrigerators

Materials properties affecting EC device operation are discussed based on an analytically tractable model of a layered EC refrigerator. Special attention was paid to thermal and interface thermal resistances. Estimates of the average cooling power of a stacked MEMS-based EC refrigerator were made.


Introduction
The electrocaloric (EC) effect is a reversible temperature change of a material which results from an adiabatic application of an electric field.EC refrigeration is an environment-friendly thermal energy conversion technology.Electric fields required for the EC refrigeration cycle can be supplied much easier and less expensively than the high magnetic fields required for the magnetocaloric (MC) refrigeration.Moreover, electrical energy for EC cooling can be provided by stationary or mobile solar cells, and by electric vehicle batteries.This opens up completely new possibilities for an environment-friendly industrialization of developing countries [1].
An EC refrigerator, i.e. a heat pump, is able to transport thermal energy against a temperature gradient from the load to the EC layer and from here to the heat sink.It consists of a polarizable EC material contacted by electrodes and thermal boundaries to the load and the heat sink.The thermal connections have to be opened and closed appropriately as the layer is heated or cooled.Heat is transferred from the load or to the heat sink either (i) via controlled heat switches as well as uncontrolled thermal rectifiers, or (ii) by pumping a gaseous or liquid heat transfer agent through the solid refrigerant [1].
The EC response of a multilayer ceramic capacitor (MLCC) under real environment and operational conditions, but without thermal connections, was considered in [6].The two-dimensional transient heat transfer model included the EC effect as a source term and accounted for radiative and convective heat losses at the interfaces.Earlier, cooling powers for MLCCs was modelled by introducing thermal capacitances and bulk thermal resistances of the capacitor sublayers [7].Later, this approach was extended to a finite element analysis of six material combinations, electrode-to-EC-layer-thickness ratio and volume of the MLCC terminals [8].A simplified one-dimensional lumped model in [9] assumed an on-state thermal contact resistance of the thermal switch of 10 -6 m 2 K/W.A two-dimensional model of a parallel plate EC refrigerator with a heat transfer fluid circulating between the plates was solved numerically by the finite element method (FEM) in [10].The model accounted for the coupled effects of temperature and fluid velocity fields.During the polarization/depolarization processes, thermal insulation was applied on all the walls.The presence of the cold/ hot heat exchangers during the fluid flow process was modelled by fixing the temperatures on the boundaries.Different EC materials (PVDF-based relaxors, PMN-PT and PLZT) have been tested with this model.
In [11], FEM was used to evaluate a micro-electromechanical system (MEMS)-technology refrigeration system based on the EC effect in a P(VDF-TrFE-CFE) terpolymer.The cooling device included two diaphragm actuators fabricated in the plane of a silicon wafer, which drive a heat transfer fluid back and forth across of 5-layer stacks of terpolymer that are placed between them.Distributions of temperature, pressure, and flow velocity of the heat transfer fluid (Galden HT-70) in the x-direction in the element as well as time-averaged cooling power density and Carnot coefficients of performance were calculated.The hot and cold ends had a constant temperature.The remaining surfaces were considered as thermally insulated bearing significant overestimations.An achievable cooling power density of 3 W/cm 2 was reported at a temperature span of 15 K and a coefficient of performance of 31% Carnot.
In [12], a thermoelectric/EC heat pump was modeled using FEM.The proposed device consists of thermoelectric layers acting as thermal switches laminated with EC layers.Ideal thermal contact between the EC and the thermoelectric layers was assumed, i.e. the contact thermal resistance was not taken into account.Any heat losses from the EC material as well as radiation and convection losses were neglected.Thus, the obtained cooling power density in the order of 1 W/cm 2 is strongly overestimated.
Recently, an EC device was numerically modeled which comprises multi-layers of EC and heat storage material films separated by thermal switches changing the contact thermal conductance by forming a larger gap [13].Calculation were carried out for best-case parameters: electric field E = 180 V/µm, EC temperature change ΔT EC = 24 K, negligible interface thermal resistances, and Fourier number Fo = 1.This results in an average cooling power of 2.4 W/ cm 2 and an average heat transfer efficiency of η = 10% for an operational frequency of f = 100 Hz and a temperature difference between the hot and the cold side of the system of ΔT hc = 20 K. Note, that Fo represents the ratio of the heat conduction rate to the rate of thermal energy storage in a solid [14].Complete heat absorption or rejection steps were realized in EC device prototypes for Fo ~ 8…10 [11,15].
In this work, we derive an analytically tractable model of a layered EC refrigerator, we identify and discuss materials properties affecting EC device operation and we illustrate the impact of these properties on device operation by means of our simplified model.

Heat diffusion equation
The variation of temperature within the EC material is described by the heat transfer equation.For the sake of simplicity, we consider in the following an one-dimensional model.In the case of linear heat flow, the parabolic heat diffusion equation of a solid comprising a heat source is given by [16] where Θ = T − T 0 is the temperature difference to the environment, α = κ/c the thermal diffusivity, κ the thermal conductivity, and c the volumetric specific heat assumed to be constant (c(E,T) ≈ c), V the volume, A the area and d the thickness of the EC layer, and h the heat transfer coefficient at the boundary describing heat loss to the environment.Introducing a bulk thermal resistance R′ th,b = (κA) -1 and a thermal capacity C′ th,b = cA, both per unit length, eq. ( 1) transforms into where Q ˙′ and h′ are the values of Q ˙ and h per unit length, respectively.Since the Biot number Bi = hd/κ (characterizing the ratio of the thermal resistance of the EC materials bulk to that at the boundary) will be less than 0.1 up to values of d in the order of 100 μm, the temperature of the EC element during heat transfer remains nearly constant, i.e. the ∂ 2 Θ/∂x 2 = 0.This enables a lumped system approximation [14] simplifying eq. ( 2) to where R″ th,i = AR th,i = 1/h is the area-specific thermal resistance at the interface, C″ th = cd the thermal capacitance of the EC layer per unit area, and Φ the heat flux through the interface.Eq. ( 3) represent the well-known bolometer equation [17].
The driving force of heat flux is the temperature difference Θ caused by charging or discharging the EC element which itself represents a dielectric capacitor.When charging a capacitor possessing a dielectric permittivity ε under isothermal conditions, the heat dQ delivered by a small change of the electric field dE is given by [18] w w § Since the EC temperature change is usually small, ΔΘ EC << Θ, the corresponding thermal flux Φ(t) = dQ(t)/(Adt) might be written as where <Τ> is the average temperature of the cooling cycle and where the field dependence of the temperature coefficient of dielectric permittivity 1/ε•∂ε/∂T is neglected.Note, that above the temperature T m of maximum dielectric permittivity, 1/ε•∂ε/∂T is negative.The rise time of an output signal of a voltage pulse generator is mainly dependent on the switching speed of the power amplifier, while the fall time is determined by the resistor-capacitor (RC) circuitry at the generator output and the load.Here, we will not consider the physical origin of power amplifier switching.Assuming for the rising pulse with τ e = RC the time constant of the electronic circuit, an approximate solution of eq. ( 3) is given by where τ th = R″ th,i C″ th and where the temperature coefficient of the dielectric permittivity was assumed to be constant.For very short rise times of the electric field, τ e →0, eq. ( 7) reduces to A pulse decay given by transforming in terms of t′ which for τ e → 0 is equivalent to eq. ( 8) with opposite sign.A flexible EC cooling device with electrostatic actuation reported in [15] was driven with a high-voltage circuitry possessing an electrical time constant of about 140 ms.On the other hand, an analysis of eqs.( 7) and ( 10) suggests time constants of less than 50 ms for an efficient device operation.
For sake of simplicity, we describe for a moment the thermodynamic cycles of the EC element by a periodic temperature change with a fundamental angular frequency ω = 2π/τ c , where τ c is the cycle time.Since the parabolic heat diffusion equation ( 1) has the same structure as the equation of an electrical RC-transmission line, multilayer problems are best treated by the matrix methods commonly used in electrical engineering [19].In our case, temperature plays the role of a voltage and thermal flux the role of an electrical current.Each point, we shall describe by two quantities, the temperature Θ and the heat flux Φ. Denoting Θ and Φ as the temperature and thermal flux of a slab at the face z = 0 and Θ′ and Φ′ as their values at the face z = d, the following relations yield [16] Θ′ = LΘ + MΦ, with and where Eq. ( 11) may be considered as matrix equation In eq. ( 14), Θ and Φ must be multiplied by their time dependence exp(iωt) which is omitted here entirely.It should be again included when real or imaginary parts have to be taken at the end of the calculation.Assuming perfect thermal contact between the faces of the slabs a comprising n sublayers, the r th being of thickness d r , thermal conductivity κ r and thermal diffusivity D r with Θ r and Φ r and Θ′ r and Φ′ r at its left-hand and right-hand faces, respectively, recursive application of eq. ( 14) gives for a one-dimensional heat flow across several sublayers Given any two of Θ r , Φ r , Θ′ r and Φ′ r the other two can be found using this method.The multiplication of matrices in eq. ( 15) can be carried out successively.However, deriving explicit formulas for slabs of n layers is very cumbersome [20][21][22].
In the presence of thermal contact resistances between the slabs or at the surfaces, they may be expressed also in matrix notation and included into eq.( 15) [16].For example, in case of a thermal interface resistance R″ th,i1 between the first and second slab we find c cc ) ) 4 4 and, correspondingly, .
The final result of this calculation is are linear relations between the temperature and fluxes Θ 1 , Θ′ n , Φ 1 , Φ′ n at the two surfaces of the composite slab.The surface conditions will provide two more equations so that all four quantities Θ 1 , Θ′ n , Φ 1 , Φ′ n will be determined.Recursive application of this matrix scheme yields for a one-dimensional heat flow across j slabs where M is the heat transfer matrix.Thus, the matrix method allows numerically a very easy analysis of multilayer structure by multiplying numerical matrices.
The relative thickness of the EC layer in prototype EC cooling devices amounts to 0.30 to 0.35 d D [11,15].In accordance with the above made approximation of a thermally thin system, we consider in the following a stack of j films satisfying d j << d D .In this case, we obtain:

Electrocaloric material
According to eq. ( 8), the EC temperature change given by ΔT EC = Θ(t = 0) is defined by the volumetric specific heat c and the temperature coefficient dε/dT of dielectric permittivity.Since the values of c at room temperature approach the high temperature limit, it is expected that they do not vary strongly among various pyroelectric materials [23].Therefore, dε/dT is a key parameter.Relaxor ferroelectrics possess above T m large values of 1/ε•∂ε/∂T as high as a few 10 2 K -1 (Table 1).
Recently, BaZr 0.2 Ti 0.8 O 3 (BZT20) was demonstrated to be a promising material for EC application [5,24].BZT20 is a composition near the invariant critical point, where four different phases, including cubic, tetragonal, orthorhombic, and rhombohedral phases coexist.It is located slightly below the crossover to relaxor behavior [28,32] where an intermediate state between ferroelectric and relaxor behavior is expected.On the other hand, BZT20 thick films show the shift of the dielectric permittivity peak towards higher temperatures with increasing frequency characteristic for relaxor ferroelectrics [24].
EC refrigerants are evaluated by a material criterion which characterizes the efficiency of the physical cooling process and, therefore, is independent on the performance of different thermodynamic cycles [1]: Here, εε 0 •E 2 •tanδ is the non-recoverable electrical loss with tanδ the loss tangent, and cΔT EC the heat transferred from the load to the heat sink within one refrigeration cy- cle with c the specific heat of the cooling element and ΔT EC the EC temperature change.At this point, the tanδ comes into play as an additional material parameter since -following eq. ( 8) -the dielectric permittivity cancels out.It has to been noted that tanδ of eq. ( 20) needs to be investigated in the frequency range of EC device operation (0.1 to 10 Hz).In the presence of a ferroelectric hysteresis, it mainly consists of losses due to minor hysteresis loops described by an apparent loss tangent [33]: where P -(E) and P + (E) are the field dependence of dielectric polarization at the upper and lower branch of the minor hysteresis loop, respectively.Obviously, tanδ depends significantly on materials synthesis.Values of tanδ < 0.07 provide refrigerant efficiencies of Φ > 0.95 exceeding the ones of any other solid state cooling technology.

Thermal interface resistance
At thin film interfaces, the highest value of R " th,i measured at room temperature to date is 1.2•10 -5 m 2 K/W for Bi/Hydrogen-terminated diamond.The thermal interface resistance at room temperature of other metal/dielectric combinations falls within a relatively narrow range, 3.3•10 -8 m 2 K/W < R " th,i < 12•10 -8 m 2 K/W [34].The minimum possible R " th,i ph-ph produced by a purely harmonic process involving two phonons (one phonon on each side of the interface) is given by the phonon radiation limit which in the limit of high enough temperatures, where the Dulong-Petit law holds, is given by Here, k is the Boltzmann constant, f max = kΘ/h the vibrational cutoff frequency of the metal, h the Planck constant, Θ the Debye temperature, and v D , the Debye velocity of a two dimensional dielectric with v l and v t the average longitudinal and transverse sound velocities, respectively.A compilation of all data used for calculations is given in Table 2.In result, eq. ( 22) yields a phonon-phonon contribution to the interface resistance of R " th,i ph-ph = 4.1•10 -10 m 2 K/W.While phonons dominate the heat conduction in dielectrics, electrons dominate it in metals.Hence, an energy transfer occurs between electrons and phonons during heat transport across the metal-dielectric interface.In result, two thermal resistances -volumetric electron-phonon (R " th e-ph ) in the metal and interfacial phonon (R " th,i ph-ph ) -are in series.The former one is given by [47] R R th "e-ph e ph ( ) . / where R e is the electron cooling rate or electron-to-phonon energy transfer per unit volume, R e = c e /τ with c e the electron heat capacity per unit volume and τ e-ph is the relaxation time characterizing electron-phonon energy loss or electron cooling.The value κ ph is the phonon thermal conductivity, with ρ the density, N A the Avogadro constant, M the molar mass, and l ph the phonon mean-free-path of the metal.
Assuming τ e-ph = 2 ps (cf.values of W and Cu in Table 2), we arrive at R " th e-ph = 8.6•10 -10 m 2 K/W.Thus, the interfacial thermal resistance in the order of 10 -9 m 2 K/W calculated by simplified models is much lower than the corresponding experimental data.
From a microscopic point of view, all surfaces in contact have deviations from their idealized geometry.Because of these imperfections, two bodies in contact touch actually only at a few discrete points.The heat flux close to the interface is then constricted in these micro-contact regions [48].
The thermal contact resistance is defined as the ratio between the temperature drop at the interface and the heat flux crossing the interface.Similarly, electrical capacitance C is the ratio of the total charge on a conductor to its electric potential (potential difference to a reference electrode).Both capacitance and thermal conductance are affected by edge effects.Therefore, there is a close relation between R th,i and the interface capacitance C i given by [49] R C where ε b is the bulk dielectric permittivity.The area-related value R″ th,i = AR th,i is then given by with d i the thickness and ε i the dielectric permittivity of the interface layer.With regard to the values of C i and ε b in Pt/ BaSr x Ti 1-x O 3 /Pt structures (cf.Table 2), R″ th,i amounts to 1…5•10 -8 m 2 K/W crossing the lower limit of the thermal interface resistance for thin film metal/ dielectric contacts reported in [34].
In macroscopic cooling systems, the thermal interface resistance is generally negligible since it is in series with other thermal resistances exhibiting much larger dimensions and, thus, much larger resistance.Here, heat is transferred from the load or to the heat sink (i) via thermal switches or (ii) by a gaseous or liquid heat transfer agents [1].A compilation of typical interface resistances of solid-solid and liquid-solid thermal contact resistances is given in Table 3.

Heat transfer fluids
In an active EC regenerator configuration, liquid or gaseous heat transfer agents are pumped through the EC refrigerant.They absorb heat from the bed and transfer it to the cold sink.Thereby, an electric field is applied solely to the refrigerant.The flow channels are usually rectangular ducts requiring a two-dimensional flow analysis.On the other hand, in a one-dimensional model, the determination of a heat transfer coefficient h (the inverse of h represents the thermal interface resistance) describing the heat transfer between the fluid and the solid is sufficient.The heat transfer across a boundary is defined by the Nusselt number Nu -the ratio of convective to conductive heat transfer normal to the considered boundary: with d H the hydraulic diameter.For steady state, laminar flow of an incompressible liquid with constant physical properties in a rectangular duct of constant cross-sectional area characterized by its aspect ratio AR (the ratio of the-smaller-to-the larger duct size), the value of Nu is approximately given by [64]: At Bi < 0.1, the temperature of the EC element during heat transfer remains nearly constant (cf.section 2.1).
Here, a duct high of z = 0.6 mm and AR = 0.15 yield thermal resistances of The thermal resistance of the fluid, R″ th,f , is determined by the inverse of the product of its mass flow rate dm/dt and its specific heat at constant pressure c p [16].The total thermal resistance of the fluid, R″ th,f , for a given fluid velocity v then yields [65]   Summing up, the thermal resistance at the interface to the heat transfer fluid limits the heat transfer of active EC regenerators.

Electrocaloric device and cooling power
A MEMS-based EC element with a stack of three cooling cells is illustrated in Figure 1.Multilayer ceramic thin films with interdigitated electrodes replicating multilayer ceramic capacitors (MLCCs) [7] are positioned on a thin SiN x membrane inside a Si frame produced by silicon wet etching.For thermal insulation, the inclined walls of the silicon die are covered by thin SiO 2 film as adhesion layer and a polydimethylsiloxane (PDMS) film as thermal barrier.The silicon dies are stacked using PDMS interlayers consisting of fluid channels.The bulk thermal resistances of the constituent layers are compiled in Table 5.A contact thermal resistance of the BZT/ Pt interface amounting R″ th,c = 5•10 -8 m 2 K/W was calculated by means of eq. ( 27).For the Pt/SiN x interface, a value of R″ th,c = ~3•10 -9 m 2 K/W obtained in literature for Pt/SiO 2 [58] was used.
The thermodynamic cycle of an EC refrigeration process in such a cooler design includes four steps [11,72]: (i) Adiabatic polarization: the EC material is heated by applying the electric field; (ii) Heat rejection: while the electric field is on, the fluid absorbs the heat from the EC material as it flows from the cold side to the hot side rejecting the heat carried by the fluid to the hot side (heat sink); (iii) Adiabatic depolarization: the EC materials is cooled by turning off the electric field; (iv) Heat absorption: the fluid moves back from the hot side to the cold side, releases heat to the EC material and absorbs heat from the cold side (load).Thus, by periodic cycling, heat is extracted from load and released to the heat sink.Here, the fluid reaches its lowest temperature in the cold side.For a cycle time τ c = 2mC″ th R″ th proportional to the thermal time constant of the EC device, the cycle-averaged cooling power of a single cell thermally isolated at the top is given by [9,72]: Here, TI denotes the thermal insulation layer.Neglecting the small bulk and interface thermal resistances, assuming R″ th,b TI → ∞ -equivalent to zero temperature gradient at the boundary between the cells -and taking into account, that the heat flux is divided between the heat transfer fluid and its PDMS environment, we arrive at ( Taking water as heat transfer fluid, assuming ΔT EC = 5 K and m = 2, the cycle-averaged cooling power of a single cell amount to 180 mW/cm 2 .Following Table 3, the corresponding < q˙″> values for silicone oil and HT-70 amount to 51.5 and 25.7 mW/cm 2 , respectively.This is comparable to an experimental value of 30 mW/cm 2 reported recently for a cooling device using a flexible electrocaloric (EC) polymer film and an electrostatic ac-tuation mechanism [15].The dependence of the averaged cooling power of a single cell on cycle time is depicted in Figure 2.
For a cell stack, only the bottom and top cells should be thermally insulated to the environment by a layer of low thermal conductivity.Here, MEMS-structures consisting of vacuum cavities [69] or highly porous materials find application.At the top and bottom surfaces, an additional heat flux divider appears which lower the cooling power of the cell by a factor of Figure 3 illustrated the dependence of F on the number of stacked cells assuming that a vertical temperature gradient is present only in the top and bottom cells of the stack.
The estimations made in this work illustrate that a solid state EC refrigerator with a cooling power of up to 2 W/cm 2 can be realized using a stack of 10 single MEMS-based cells.

Conclusions
A simplified analytical model of solid state EC refrigerators allows analysis of the materials impact on device operation, enables the indentification of the key material parameters and makes it possible to predict the device performance.For a stack of MEMS-based EC cells in an active regenerator mode, cycle-averaged cooling powers of up to 2 W/cm 2 were estimated.
is the penetration depth of the thermal oscillation.Note, that for dielectric thin and thick films with d < 100 µm we get L = O ≈ 1, M = -d/κ, and N = -κk 2 d for ω = 1 Hz since D is in the order of 10 -6 m 2 /s.The value of M represents the area-specific thermal resistance of the film bulk.

For
calculation, an upper limit of the fluid velocity of v = 0.1 m/s was chosen providing a back and forth fluid motion of 1 cm at a cycle time of τ c = 0.2 s.This results in vc >> h leading to R″ th,f ≈ R″ th,s-f .Values of R″ th,s-f are compiled in Table 4. Also, the thermal penetration depth at ω = 1 Hz, d D,1Hz , is listed.It allows a quick estimation of d D at a given cycle time by

Figure 2 .
Figure 2. Cycle-averaged cooling power of a single EC cell in dependence on cycle time.

Figure 3 .
Figure 3. Efficiency factor versus number of stacked cells.

Table 1 .
Temperature coefficient of dielectric permittivity above T m for BZT20 and selected relaxor ferroelectrics.

Table 2 .
Material parameters used for the calculation of the interface thermal resistance.

Table 3 .
Thermal contact resistances of solid-solid and liquid-solid configurations.

Table 4 .
Physical properties of heat transfer fluids, thermal resistance at the solid-fluid interface R″ th,s-f for laminar flow, and values of d D,1Hz .

Table 5 .
Bulk thermal resistances of the constituent layers of the EC device in Figure1.