Keywords: Computer Animation, Modal Analysis, Trees, Real-Time Animation virtuelle, analyse modale, arbres, temps-réel

This scientific report presents a real-time method to animate complex scenes of thousands of trees under a user-controllable wind load. Firstly, modal analysis is applied to extract the main modes of deformation from the mechanical model of a 3D tree. The novelty of our contribution is to precompute a new basis of the modal stress of the tree under wind load. At runtime, this basis allows to replace the modal projection of the external forces by a direct mapping for any directional wind. We show that this approach can be efficiently implemented on a graphics hardware. This modal animation can be simulated at low computation cost even for large scenes containing thousands of trees.

Base de projection du vent pour l'animation temps-réel d'arbres

Vegetation is an important part of the realistic depiction of natural scenes in a virtual environment. While research in plant modeling has produced impressive results, animation techniques remain of visually poor quality, especially for realtime applications.

The current development of computer hardware allows real-time applications to render natural scenes with high quality geometric details. However, the increase in model complexity also generates the need for suitable animation techniques. As the existing accurate simulation methods are still too computationally expensive, developers usually resort to ad-hoc techniques. Such methods can be satisfying for simple models such as grass or palm trees, but fails to provide realistic animation of complex branching structure such as standard trees.

Similarly to [START_REF] Stam | Stochastic dynamics: Simulating the effects of turbulence on flexible structures[END_REF], we rely on a modal analysis to animate 3D models of trees. This technique simulates the deformation of the plant using a precomputed deformation basis. While it provides visually convincing results, it still requires full integration of the dynamic wind forces over the entire plant at runtime. Real-time framerates can be obtained only for very simple models and not for realistically complex trees.

Our contribution is a method that makes the computation of the wind forces integration almost costless and independent of the tree complexity. The only remaining computations at runtime consist in the update of the positions of animation control points, defining the skeleton used to animate the rendered geometry. To this end we define a modal projection basis where the influence of directional wind loads can be precomputed off-line.

After introducing modal animation of trees in section 3, we show in section 4 how the wind projection basis is computed and how it is used to avoid the integration of wind forces at runtime. Finally some implementation issues are discussed in section 5. We show that an entire forest containing thousands of trees can be animated in real time using graphics hardware, producing realistic results for both close and distant trees.

State of the art

Many models defined for real-time animation of trees use mechanical simplification specifically designed for that purpose. For example, in [START_REF] Ota | A hybrid method for real-time animation of trees swaying in wind fields[END_REF] and [START_REF] Akagi | Computer animation of swaying trees based on physical simulation[END_REF] the tree dynamics is modeled using the static response of each branch segment (i.e which does not consider time integration of the dynamics). However both methods focus instead on the wind formulation. Ota To animate forests in real-time, Di Giacomo et al. [START_REF] Di | An interactive forest[END_REF] also use static response to animate the tree's terminal branches. For the rest of the structure they simulate the dynamics of the rotation angles between branches using an explicit Euler method, in a similar way as Sakaguchi and Ohya in [START_REF] Sakaguchi | Modeling and animation of botanical trees for interactive virtual environments[END_REF].

These previous works model the tree's mechanical response expressed in generalized coordinates. However they solve the dynamics locally (i.e segment by segment) and do not fully model the inertia of the subsequent branches that are displaced by a joint's rotation. Simulation techniques of global mechanical structure such as [START_REF] Baraff | Linear-time dynamics using lagrange multipliers[END_REF] provide correct dynamics without this limitation. Weber [START_REF] Weber | Fast simulation of realistic trees[END_REF] developed such a simulation method adapted to trees.

For real-time application only little computation time is commited to the animation of vegetation. The use of parallelizable algorithms has become essential. However animation described by generalized coordinates makes the computation of nodes position hierarchically dependent on each other. Weber proposes a scheme to use multithreading, but it is not yet appropriate for highly parallel computing provided by graphics processing unit (GPU).

Data driven methods are interesting alternatives to direct simulation. Diener et al. [START_REF] Diener | Hierarchical retargetting of 2D motion fields to the animation of 3D plant models[END_REF] use video to extract motion data from real plants and reproduce it on virtual models. James et al. in [START_REF] James | Precomputing interactive dynamic deformable scenes[END_REF] and [START_REF] James | Mesh ensemble motion graphs: Data-driven mesh animation with constraints[END_REF] use precomputed simulation as input to generate realistic response to external impulse and to fluctuating wind respectively. However the main disadvantage of these approaches is the restricted controls at run time.

First introduced to the computer graphic community by Pentland and Williams [START_REF] Pentland | Good vibrations: modal dynamics for graphics and animation[END_REF], modal analysis has proven to be a powerful basis for the simulation of oscillatory behaviors. It has then be used for many applications such as precomputed response in an animation using graphics hardware [START_REF] James | Dyrt: dynamic response textures for real time deformation simulation with graphics hardware[END_REF], real-time collision and manipulation [START_REF] Hauser | Interactive deformation using modal analysis with constraints[END_REF] and even sound generation [START_REF] James | Precomputed acoustic transfer: Output-sensitive, accurate sound generation for geometrically complex vibration sources[END_REF].

Modal analysis provides a basis of deformations where all elements of the basis are independent harmonic oscillators. For trees, Shinya and Fournier [START_REF] Shinya | Stochastic motion-motion under the influence of wind[END_REF] applied this approach to model each of the branches. Stam [START_REF] Stam | Stochastic dynamics: Simulating the effects of turbulence on flexible structures[END_REF] was the first to use modal decomposition of the entire structure to animate plants. The main idea is to generate a stochastic wind and compute the modal dynamics directly in the frequency domain using the specific response of harmonic oscillators. However the main limitation of this modal framework for real-time animation, is the need to compute the modal projection of the wind forces (i.e forces integrations over the whole tree) at runtime.

Our contribution consists in a method that drastically reduces the computations needed to simulate modal dynamics. This is done by precomputing a basis of the wind modal projections. Moreover the dynamics is simulated with a computational cost independent of the tree structure complexity. This makes the implementation of level of detail straightforward as only the position of tree nodes used for rendering has to be computed.

INRIA 3 Modal animation framework

In this section, we briefly describe the mechanical model. We introduce modal analysis and explain how it can be used in a real-time animation framework.

Mechanical modeling and modal analysis have been thoroughly studied. We advise interested readers to refer to specialized literature such as [START_REF] Gradin | Mechanical Vibrations: Theory and Application to Structural Dynamics[END_REF].

Mechanical model

To solve dynamics over time we resort to numerical simulation applied on a discretized model of a tree obtained using the Finite Element Method (FEM). The principle is to decompose the mechanical system in a set of elements separated by nodes. The continuous system's deformation is defined as a specific interpolation between the nodes displacement.

In our case, as in [START_REF] Stam | Stochastic dynamics: Simulating the effects of turbulence on flexible structures[END_REF], the system domain Ω is the lineic representation of the tree, i.e. the set of central axes of the branches. Then the FEM discretizes Ω as a set of segments Ω e (see figure 2). Let u(t) be the vector containing the displacements of all the nodes at time t. The dynamics of a dissipative system is commonly expressed by the equation of motion:

M ü(t) + C u(t) + Ku(t) = f (t) (1)
where f (t) is the external force (see section 4), M and K are the mass and stiffness matrices. The commonly accepted model for the damping matrix C is a linear combination of M and K :

C = αM + βK (2)
More details on the finite element representation and the matrices used in equation ( 1) are given in appendix.

Modal analysis

Modal analysis is the decomposition of the deformations of a mechanical system into a set of special deformations called vibration modes. These modes are the solutions u of the free vibration equations:

M ü(t) + Ku(t) = 0 (3) u = λ i (t) • ϕ i (4)
where the ϕ i are constant vectors called the modal deformations, and the λ i (t) are scalar functions of time. Substituting (4) in ( 3) gives:

M -1 Kϕ i = µ i ϕ i with µ i = - λi (t) λ i (t) (5)
As K, M and ϕ i are constants, µ i are constant scalars. Finally, the eigenvectors and eigenvalues of M -1 K are the solutions of equation [START_REF] Diener | Hierarchical retargetting of 2D motion fields to the animation of 3D plant models[END_REF]. Let Φ be the matrix containing all the modal deformations ϕ i , and q(t) the expression of vector u in eigenspace:

Φq(t) = u(t) (6)
Substituting [START_REF] Di | An interactive forest[END_REF] in equation ( 1) and multiplying on the left by Φ T , we get:

M q(t) + C q(t) + K q(t) = Φ T f (t) (7)
where M = Φ T M Φ, C = Φ T CΦ and K = Φ T KΦ. A classical result of modal analysis is that these matrices are diagonal [START_REF] Gradin | Mechanical Vibrations: Theory and Application to Structural Dynamics[END_REF]. We can then rewrite equation ( 7) as 6n independent equations such that, for each element q i of vector q:

qi (t) + γ i qi (t) + ω 2 i q i (t) = 1 m i f i (t) (8)
with

   ω 2 i = k i m i γ i = α + β * ω 2 i (9)
where f i (t) is the modal force induced by the wind (see section 4), ω i the natural frequency, γ i the damping coefficient of mode i, m i and k i the modal mass and damping respectively (i.e the diagonal elements of M and K).

Animation using the modes of deformation

Modal analysis is well suited to our purpose because it extracts the most representative components of the structure deformations. For the case of plants, the vibration modes associated with low frequencies (ω i ) have the most significant (and largely distributed) influence on the final motion [START_REF] Rodriguez | A scaling law for the effects of architecture and allometry on tree vibration modes[END_REF]. Hence only a small subset of the mode matrix Φ needs to be kept (the columns corresponding to the highest eigenvalues).

Equation [START_REF] Hauser | Interactive deformation using modal analysis with constraints[END_REF] means that each mode behaves as a one-dimensional harmonic oscillator (see figure 3). Each vibration mode can be seen as a deformation of the whole structure oscillating independently from the others.

At each time step, the simulation of modal animation consists in the following procedure:

1. Compute f i (t), the modal projection of the wind. 2. Update the dynamics of the modes using the explicit Euler method:

INRIA

qt+δt i ← f i (t)/m i -γ i qt i -ω 2 i q t i qt+δt i ← qt i + δt qt+δt i q t+δt i ← q t i + δt qt+δt i 3.
Compute the updated displacement:

u t+δt = ∀i q t+δt i ϕ i (10)
Figure 4: The modes of deformation of the walnut model. During animation the final displacement of a tree is a combination of its modal deformation. The modes are sorted in ascending order of their natural frequencies.

Finally, only the evaluation of the modal forces f i (t) requires costly computations. This is the main difficulty that our technique proposes to resolve. To this end we now describe a basis that can be computed off-line and avoids costly computations at runtime.

Wind projection basis

Integrating a general wind load over the whole tree requires a sum over all the elements (branches) at each time step of the animation (as done by [START_REF] Stam | Stochastic dynamics: Simulating the effects of turbulence on flexible structures[END_REF]). Doing so prevents real-time animation of sufficiently complex trees. By making the assumption that the wind load is uniform over the whole tree, a projection basis can be precomputed that allows to animate the tree in real-time under the load of any directional wind.

Note that the wind is considered uniform only over each tree but can vary between disctinct instances. The speed and direction of the wind are left unconstrained and can be freely controlled at runtime.

We start by modeling the drag load induced by the contact of wind on a rigid structure. As in [START_REF] Sellier | A finite element model for investigating effects of aerial architecture on tree oscillations[END_REF][START_REF] De | Effects of wind on plants[END_REF], we model the wind load on any point p ∈ Ω at time t by:

f (p, t) = 1 2 ρC D D(p) v(p, t) -ṗ (v(p, t) -ṗ) (11)
Where v(p, t) is the wind speed vector and D(p) is the diameter of the branch at point p. The air density ρ and the drag coefficient C D , depending on the shape of the branches, are We merge them as one scalar factor C.

From this local definition, one can derive the modal wind load for each mode i by modal projection [START_REF] Gradin | Mechanical Vibrations: Theory and Application to Structural Dynamics[END_REF]:

f i = Ω f (p) • ϕ i (p)dp (12)
where we omit the time parameter for clarity, and where ϕ i (p) is the vector containing the i th modal deformation at p.

To simplify computations, we do not take into account the direction of branches in equation [START_REF] James | Dyrt: dynamic response textures for real time deformation simulation with graphics hardware[END_REF] for the drag coefficient. This simplification is reasonnable because the modes with low frequency are mostly perpendicular to the branches. The influence of the branches parallel to the modal forces is then essentially negligible.

In the case of a freely defined finite element discretization, equation ( 12) is computed as the sum of the elements contribution:

f i = e Ωe f (p) • R e N e (p)ϕ e i dp (13)
where R e is the rotation matrix mapping the local frame of the element e to the global frame. ϕ e i is the vector containing the i th modal deformation of both element's border nodes (in the local frame) and N e (p) is the shape function matrix defining the interpolation between these nodes (see appendix).

To compute the modal forces, equation ( 12) can be discretized over the nodes of the FEM (in the frequency domain) as in [START_REF] Stam | Stochastic dynamics: Simulating the effects of turbulence on flexible structures[END_REF]:

f i = n∈Ω f (n) • ϕ i (n)
For simple structure it is a suitable approach but for complex trees such a modal projection would be prohibitive. Using a simplified structure is not a solution as it would result in a poor approximation of the effective modal stress.

In our method, to reduce the cost of the modal projection of the wind forces at run time, we focus on extracting the time dependent parameters out of the integral. The remaining part is then precomputed.

The first simplification is to consider R e as a constant: the wind load is computed on the tree in its rest position. The orientation changes of the branches is not significant as they remain small for natural trees motion. Compensatory INRIA models were tested but did not add noticeable improvement. With this approximation, the modal deformation term is constant and only the wind force depends on time. To extract runtime parameters, we assume a wind model such that:

1. From the wind force, we only keep the mean aerodynamic stress:

f (p) = C D(p) v(p, t) v(p, t)
The remaining part of the external force, called the aerodynamic can be compensated by an increase of the damping coefficients α and β from equation (2).

2. We define the wind with a variable direction and mean speed, but we assume it to have an uniform value v(t) over the tree.

Using these simplifications, equation ( 12) becomes:

f i = Ω C D(p) v(t) v(t).φ i (p)dp (14) = C v(t) v(t). Ω D(p)φ i (p)dp (15)
It results that if we precompute the wind projection vectors p i defined as:

p i = C Ω D(p)φ i (p)dp
then the modal wind load can be easily computed at runtime using:

f i = v(t) v(t).p i

Implementation issues

The mechanical modeling, modal analysis and the precomputation of the projection basis require the use of complex algorithms that must be performed using specialized tools. We used the software Cast3m 2000 [START_REF]Cast3m[END_REF] specialized in the study of mechanical structures using the finite element method.

GPU implementation

A powerful aspect of the modal representation is that most computations are completely independent. It thus allows to take full advantage of the massive parallel computing power offered by present graphics hardware (GPU). This enables us to animate large forest scenes in real time (see section 6).

The two main parameters of our tree model are the number of nodes of the skeleton and the number of modes kept for animation.

1. As can be seen in equation [START_REF] James | Precomputing interactive dynamic deformable scenes[END_REF] the updated position of each node can be computed independently from the others using only the dynamic modal state (q i ). Thus the nodes positions can be stored in a texture in GPU memory (matrix U of figure 5) and updated using a fragment program and Figure 5: Computation steps of the animation framework. K is the number of tree instances, M the number of modes used for the animation and N the number of control points of the skeleton. The colored rectangles represent textures stored in GPU memory. The input matrices T , P, Φ are constant while the matrices Q, U store the dynamic variables. T contains the tree instance data (position, orientation); P constains the modal parameters used to compute the dynamics; Φ constains the modal deformations of the skeleton control points; Q contains the modal states (q i , qi ); U contains the skeleton's node displacements. Several representations can be used for the wind such as a procedural equation or a flow array.

an offscreen buffer. Animating the tree mesh can then be done by reading this texture in a dedicated vertex program. Any level-of-detail rendering model using a (possibly time-varying) subset of the animation nodes can be used, without the need to update the remaining unused nodes.

2. Adjusting the number of modes used for the computation of updated nodes positions gives a handful tradeoff between animation quality and framerate. For distant views of forests only the first few modes needs to be kept, while adding more modes at close range adds detail to the animation. As for the nodes, each modal state (q i , qi ) can be updated independently from the others, they are stored in a texture in GPU memory (matrix Q of figure 5), and updated by a dedicated fragment shader.

The overall animation process can be summarized in three main steps (see figure 5).

1. Update of the modal state for each mode of each tree (item 1 and 2 of section 3.3). It requires to compute the local wind speed for each tree instance.

2. Deformation of the skeletons from the rest position (item 3 of section 3.3)

3. Rendering each tree with deformed geometry using skinning on the animated skeleton.

INRIA

Error correction

Modal deformation being a linear approximation of the displacement of the branches nodes, it is mostly adapted for low amplitude motion. Notably, branches lengths are not kept constant when submitted to strong wind. However as can be seen on figure [START_REF] Di | An interactive forest[END_REF], the corrected displacement of a node, when a single mode is selected, lies on a smooth trajectory which can be efficiently approximated by a quadratic function (using a non linear reparametrization):

u = ∀i qi • v i + q2 i • w i ( 16
)
with qi = a arctan(b

• q i ) ( 17
)
where vectors v i , w i and scalars a, b are optimized using sampled corrected deformations. The results show that this approach yields a sufficiently good correction for branches lengths to allow the application of a strong wind load.

Figure 6: Correct trajectories obtained by selecting each mode successively (each color corresponds to one mode). These trajectories can be faithfully approximated by second degree polynomial curves.

Results

We show results of our method on several models.

1. A walnut tree (from [START_REF] Sinoquet | Assessment of the three-dimensional architecture of walnut trees using digitising[END_REF]) where the trunk and all the branches were digitized from a real tree (see fig. 8). For this model 25 modes are sufficient to produce detailed animations. We used a version containg 3437 nodes, but it should be possible to reduce this number without impairing the animation quality.

2. An oak model generated using [START_REF] Runions | Modeling Trees with a Space Colonization Algorithm[END_REF] (see figure 7) Table 1: Animation framerates using our method depending on the number of trees (horizontal) and number of control points (vertical) with 25 modes of deformation. These results were obtained on an Intel Xeon 3.20GHz with a GeForce 8800 GTS.

Table 1 shows framerates obtained using our technique. Note how the rendering speed keeps constant when more animation nodes are used. The supplemental video demonstrates the realism of the results and some of the possibilities it offers.

INRIA 7 Limitations and future work

We do not address the problem of animating single leaves as modal analysis is not suitable for this task. We do think that leaves should instead be animated independently using techniques such as [START_REF] Ota | A hybrid method for real-time animation of trees swaying in wind fields[END_REF] or animated textures. In our implementation, the leaves rigidly follow the nodes they are attached to. However the complexity of branches animation already provides very convincing motion.

Mathematically, the main limitation of our method is the assumption of spatial uniformity of the wind over each tree. It does not lower visual realism as the modal dynamics keeps natural oscillatory behaviours. However precomputing a more expressive basis would allow to increase animation precision.

Another restriction is that the attenuation of the wind by the tree is not taken into account (i.e the branches in the 'back' of the tree should receive less wind).

Finally the parallelism of the method coupled with efficient GPU programming enables the animation of thousands of trees. It would require a more aggressive level-of-detail scheme to allow realt-time rendering of even larger forest scenes.

Conclusion and acknowledgement

We showed a new approach to compute the modal projection of the wind load allowing a drastic reduction of computations at runtime. To this end, we introduce a precomputed basis for the projections of directional wind in a modal animation framework.

Our implementation of the method shows that it is possible to efficiently animate and render thousands of trees. We think the presented technique is perfectly adapted to real-time applications such as computer games or simulators.

We would like to thank Adam Runions for providing us many high quality tree models.

A node has thus six degrees of freedom (see fig. 9): the three spatial displacements (u x , u y , u z ), the torsion around the main axis (ψ x ) and the derivatives ψ y and ψ z defined by [START_REF] Gradin | Mechanical Vibrations: Theory and Application to Structural Dynamics[END_REF]:

ψ y = -
du z dx and ψ z = du y dx This interpolation is defined such that, at time t, the displacement δ p (t) of any point p ∈ Ω e is: δ p (t) = N e (p)u e (t) [START_REF] Sellier | A finite element model for investigating effects of aerial architecture on tree oscillations[END_REF] where u e (t) is the 12-dimensional vector containing the displacements of the two element's border nodes, and N (p) the shape function matrix is: where the dots represent zeros, l 1 and l 2 are the coefficients of a linear interpolation and the h i are the cubic Hermite basis functions.

N e (p
In accordance with most models of trees, branches are modeled as cylindrical beams [START_REF] Sellier | A finite element model for investigating effects of aerial architecture on tree oscillations[END_REF]. The local stiffness K e and mass M e matrix are defined as the discretization of the potential energy ν and kinetic energy τ of an element of length l Here the dot notation refers to time derivatives and δ = d dx . In our implementation, we typically use the following values: E = 10 10 P a, G = 2.6E, ρ = 10 3 kg.m -3 . From the beam model used, we have I = πr 4 /16, J = √ 2I, A = πr 2 where r is the average radius of the element. These local matrices are then transformed to express the energies in function of the nodes displacement defined in the global frame of the tree. Finally, they are assembled together in the global matrices M and K of equation [START_REF]Cast3m[END_REF]. A description of this procedure has been given by [START_REF] Stam | Stochastic dynamics: Simulating the effects of turbulence on flexible structures[END_REF], or can be found in specialized literature such as [START_REF] Gradin | Mechanical Vibrations: Theory and Application to Structural Dynamics[END_REF].

INRIA

Figure 1 :

 1 Figure 1: A single tree and a forest animated bythis method.

  et al. propose a realistic stochastic description. Akagi et al. incorporate the influence of the plant in the computation of the wind flow modeled with the Navier-Stokes equation.

Figure 2 :

 2 Figure 2: Finite element discretization.

Figure 3 :

 3 Figure 3: A vibration mode is a deformation of the whole tree behaving as a harmonic oscillator.

Figure 7 :

 7 Figure 7: The oak tree model.

Figure 8 :

 8 Figure 8: Digitalized walnut tree. 1000 2000 4000 2214 40 21 11 3437 40 21 11

0 EA(δu x ) 2 +

 02 EI((δ 2 u y ) 2 + (δ 2 u z ) 2 ) 2 + GJ(δψ x )

  . . . . . l 2 . . . . . . h 1 . . . h 2 . h 3 . . . h 4 . . h 1 . h 2 . . . h 3 . h 4 . . . . l 1 . . . . . l 2 . .

			l 1 
	) =	  	  

RR n 6674

Appendix

In the lineic representation of the tree, each element is defined in its local frame where the x-axis is in the direction of the segment at rest. The deformation of an element is the interpolation of its border nodes displacement: linear for the main axis and Hermitian for the y and z axis.


Lucas Depreux

Morgan Almanza

Fabien Parrain

Martino Lobue

Study of the actuation losses in an electrocaloric cooling device

Keywords: Electrocaloric, Cooling, Electrostatic, Actuator, Viscous losses

, where a planar electrocaloric film is electrostatically actuated to move between the two reservoirs. Relying on a Poiseuille flow model, we examine the trade-off between the actuation time and the viscous losses, and its impact on the electrocaloric cooling power. Our model predicts that the cooling power can be increased up to 325 mW cm -2 .

INTRODUCTION

Electrocaloric materials (ECM) undergo entropy and temperature changes triggered by an electric field under isothermal and adiabatic conditions, respectively. Among emerging technologies, electrocaloric cooling is one of the most promising candidates to replace the traditional vapor compression refrigeration systems with solid-state working substances [START_REF] Moya | Caloric materials near ferroic phase transitions[END_REF]. So far, devices using caloric materials as an active component have hardly been able to reach the specific cooling power and the coefficient of performance (COP) of gas-based refrigerators and heat pumps, so that designing a solid-based heat pump with a high throughout is nowadays one of the key technological challenges.

Here, we focus on a cooling device design that is similar to the one presented by Ma and coworkers [START_REF] Ma | Highly efficient electrocaloric cooling with electrostatic actuation[END_REF]. In their device, they use an electrocaloric-polymer film switching between the hot and cold end ; a COP up to 13, a specific cooling power of 2.8 W g -1 and a working frequency of 0.8 Hz have been reported.

More precisely, in their device, the active film is electrostatically actuated to follow a zipping motion between two planar reservoirs, keeping an S-shape as detailed by Sato et al. [START_REF] Sato | An electrostatically actuated gas valve with an sshaped film element[END_REF]. Both faces of the film are covered with a conductive layer, so that when a voltage is applied between the top or bottom electrode and one of its faces, the film progressively zips to the electrode. Thus, the active substance is also the heat exchanger getting rid of the fluids and pumps needed in traditional regenerator designs, a configuration that is well-suited for light and compact actuators. The working cycle consists of four steps illustrated in Fig. 1: (a) the film heats because of the electric field and delivers the heat to the hot side, (b) the film is moved through electrostatic actuation to the cold side, (c) the film cools as the electric field is removed, and absorbs heat from the cold side, (d) the film is actuated back to the hot side.

The overall transferred heat over a cycle from the cold to the hot side Q c depends on the intrinsic properties of the ECM and on the mass of the film. On the one hand, the cold power results from the product of the transferred heat and the working frequency, which depends on the heat exchange time and the The electrocaloric cooler device is based on a thin electrocaloric film (green), with a layer of metallization on both faces (dark gray lines). These are also used to actuate the film towards the electrodes (black lines) on each reservoir (light gray), using electrostatic forces (U zip ), the speed of the film being u f ilm . On each end, the film is stuck to one of the two electrodes with an insulation layer in-between (yellow), and the film and the electrodes are separated with a spacer (green) of thickness h. When the film is colored in orange, an electric field is applied to trigger the electrocaloric film (the film is heated and transfers Q hot to the hot reservoir), and when the film is blue, the field is removed (the film is cooled and receives Q cold from the cold reservoir). u air is the velocity field of the air, and E loss is the energy lost through one actuation.

actuation time. On the other hand, the efficiency depends on the input work required to transfer heat and on the energy lost during a cycle.

Our main goal here is to find the way towards increasing frequency (i.e. power) without giving up on efficiency. The duration of a full cycle corresponds to the two heat exchange stages (i.e. steps (a) and (c) described above) and the two actuations (i.e. the film switching time when moving from one reservoir to the other, steps (b) and (d)); hence the cycle frequency is

f = 1 2(t exch + t act ) , (1)
where for the sake of simplicity we assume that (b), and (d) stages actuation times t act , as well as (a), and (b) exchange times t exch , are equal two by two.

Using the same actuation principle, we have recently shown that the exchanged heat flux can be significantly enhanced [START_REF] Almanza | Electrostatically actuated thermal switch device for caloric film[END_REF]; therefore the exchange time t exch can be reduced down in the 10 ms range. Here, we shall focus on the cost of reducing t act in terms of actuation energy, and we will also examine the optimal balance between exchange and actuation times at a given frequency.

However, while increasing the force of actuation (i.e the voltage U zip in Fig. 1) will allow to move the film faster, it will also generate additional losses, to the detriment of the efficiency. Accordingly, we need to study the energy-loss versus actuation-time trade-off, in order to determine by how much we can increase the frequency without significantly affecting the efficiency.

As our device is MEMS-sized (Micro-Electro-Mechanical-Systems), this makes the energy losses due to the viscous flow of air the predominant dissipative term; we will therefore study this term, neglecting other known loss sources (e.g. film deformation, Joule heating due to electrodes charging, finite time heat exchange, etc.).

We will start dealing with the problem by using the Hagen-Poiseuille flow model, a quite common first-approximation in MEMS literature [START_REF] Gad-El Hak | MEMS: Introduction and Fundamentals[END_REF][START_REF] Bao | Chapter 3 -air damping[END_REF]. Using this model, we will estimate the energy lost due to actuation, and then compare it to heat values reported in the literature [START_REF] Moya | Caloric materials near ferroic phase transitions[END_REF][START_REF] Basso | Direct measurement of the electrocaloric effect in poly(vinylidene fluoride-trifluoroethylene-chlorotrifluoroethylene) terpolymer films[END_REF]. As the viscous losses heat the device, they are detrimental to the cooling effect. Hence, we use a total cold flux that takes into account the negative contribution of losses. This allows us to illustrate how beneficial increasing the frequency is to the cold flux, and how the losses limit the optimal frequency of the device. Finally, we show that reducing the air density in the device would also reduce the losses, and therefore it would allow further increase to the cooling power.

FLOW MODELING

Here we shall model viscous flow of air between the surface of two parallel planar electrodes of length L, and width w, with an electrocaloric film of thickness d between them. The middle part of the film is elastically bent into an S-shape, the bent thickness being δs. The distance between the electrode surfaces is h + d as shown in Fig. 2.

As d, h, δs L the S-bent can be approximated with a straight vertical interface located in x = s (viz. the dashed magenta line in Fig. 2). This allows to skip all the interface details related to the bent shape focusing on the preponderant volume phenomena associated with the air flow, and considering the Sbent as a moving piston.

When a voltage is applied between one of the two pairs of electrodes, the piston is translated along the x axis with speed depending on the applied electric voltage U zip [START_REF] Sato | An electrostatically actuated gas valve with an sshaped film element[END_REF].

When the film is zipping from one electrode to the other, the S-bent pushes the air in front of it over a distance L, and drags the air behind. In the δs L approximation where the S-bent behaves like a planar interface the gas moving in front and behind can be described as a plane viscous flow.

The most common way of dealing with such an issue is to use the parallel plane Hagen-Poiseuille flow model [START_REF] Bao | Chapter 3 -air damping[END_REF] assuming the flow to be stationary, unidirectional (along x axis), isothermal, incompressible and laminar, the fluid to be Newtonian. where u is the axial velocity of the flow, and the force balance equation is

dp(x) dx = µ d 2 u(x) dy 2 , ( 3
)
where p is the pressure, and µ the viscosity of the air. In order to solve this equation, we need boundary conditions for velocity and pressure. The no-slip condition is generally used for channels; in this case u = 0 at the walls. By integrating equation (3) on y-axis, where the solution is a parabola for the velocity profile, we find the well-known Hagen-Poiseuille equation:

dp(x) dx = 12µ wh 3 D v (4)
with the flow rate D v = whu f ilm . For instance, when the film moves forward (viz. in the positive x axis direction), it compresses the air in front (the air is pushed towards x = L) and generates a depression behind (the air is drained from x = 0). The pressure is p + f in front of the film and p - f behind. Therefore, we can consider two separate flows: the one in front of the film s ≤ x ≤ L and the one behind 0 ≤ x ≤ s, that are connected through the outside pressure p out . By integrating Eq. 4 along x-axis, from s to L and for 0 to s, we have

∆p tot = p + -p -= 12µL wh 3 D v , (5)
showing that the flow generated by the movement of the piston along a distance L is equivalent to the flow in a rectangular pipe of equivalent length. Now considering t act as the time needed to the film to sweep the length L, namely t act = L/u f ilm = (whL)/(D v ), from (5) we get the following expression for the actuation time,

t act = 12µL 2 h 2 ∆p tot (6)
The dissipated power can be written as ∆pD v , whence the viscous energy loss per actuation ∆pD v t act writes:

E visc = 12µwL 3 ht act . (7)
As the air is accelerated in the channel, the viscous losses is not the only energy that the actuator need to provide. It is worth noting that here we focus on the losses associated with the steady state flow taking into account nothing but the air kinetic energy, namely:

K air = 1 2 (ρ air whL)u 2 f ilm χ = 3 5 ρ air whL 3 t 2 act , (8)
where ρ air is the air density, and χ = 6/5 is the integral of the square of the ratio of the velocity over the mean velocity.

LOSSES AND COOLING POWER

One of the the relevant quantities defining the performance of a cooling device is the heat flux through the exchange surface, which is why energies and works are commonly surface normalized. Consequently, from here on we shall switch to surfaceenergies using the plates area A = wL as the normalization factor.

The surface energy lost through viscosity is therefore

e visc = 12µL 2 ht act . (9)
Let us now consider the the worst-case scenario where the inertia of the fluid is entirely dissipated in the air at each actuation. The surface kinetic energy of the air in the channel of a Poiseuille flow is:

k air = 3 5 ρ air hL 2 t 2 act , (10)
the total energy loss is

e tot = e visc + k air = 12µL 2 ht act + 3 5
ρhL 2 t 2 act ( 11
)
Now we shall find the optimal thickness h opt by minimizing e tot (h),

h opt = 20µt act ρ air ( 12
)
which corresponds to the thickness value where the kinetic energy equals the viscous loss. The minimized total energy is then

e tot,min = 12 √ 5 √ µρ air L 2 t 3/2 act . ( 13
)

Electrocaloric effect and cold flux

The two figures of merit of ECM are the adiabatic temperature change associated to the electric field change ∆T adia , and the heat T ∆S that the material can exchange when the field is applied or removed along an isothermal transformation. Polyvinylidene fluoride polymers films are promising electrocaloric materials for which an adiabatic temperature change of 10 K and an isothermal heat exchange of 20 kJ kg -1 have been reported [START_REF] Moya | Caloric materials near ferroic phase transitions[END_REF]. Nevertheless, to give a reasonable estimate of the performance of an ECM mounted into a device under actual working condition, we shall use ∆T adia = 1.1 K and T ∆S = 1.2 kJ kg -1 as measured with a calorimetric setup by Basso and coworkers [START_REF] Basso | Direct measurement of the electrocaloric effect in poly(vinylidene fluoride-trifluoroethylene-chlorotrifluoroethylene) terpolymer films[END_REF]. The electrocaloric material exchanges a fraction of T ∆S that depends on the thermodynamic working cycle. Here we assume that the exchanged heat is

Q c = ρ ec T ∆S 2 dwL, (14)
where ρ ec is the ECM density (1800 kg m -3 ), and dwL the volume of the film. The surface heat received from the cold source at each cycle is then

q c = ρ ec T ∆S 2 d. (15)
Increasing the polymer film thickness is a possible strategy to improve the heat exchanged per cycle. However, even when assuming a perfect thermal interface, the heat flux will be limited by the thermal diffusivity of the film, α. Here we describe the role of ECM diffusivity introducing the phenomenological thermal time constant τ ,

τ = d 2 α , (16)
which we will use to approximate t exch . Neglecting the actuation losses and working with a vanishing actuation time we get the following cooling flux expression

qideal = q c 2t exch = ρ ec T ∆S.α 4d ( 17
)
where reducing the film thickness appears as a straightforward way to overcome the slowness of heat diffusion (decreasing t exch ), and thus to achieve a relevant power increase.

However, in an actual device, the increase of actuation losses associated with the reduction actuation time can hardly be neglected. On the contrary, as shown in Fig. 3 where the actuation looses are compared with the cold flux for different film thickness, they can dramatically increase negatively affecting the overall performance of the device. Actually, for actuation times below 1 ms, the actuation losses exceed the cold flux of the electrocaloric effect. Above 1 ms of actuation time, the cold flux is predominant, and for higher actuation time the viscous loss eventually becomes negligible with respect to the cold flux. Moreover, we see in Fig. 3 that the cold flux is modulated by the film thickness which is related to the exchange time.

Cooling power of the device

A simplified way to assess the detrimental impact of the actuation losses on the cooling device performance, is to estimate the overall power using a total exchanged heat, which is the difference between q c and two times (two actuations per cycle) the aforementioned losses, namely

q tot = q c -2e tot,min , (18)
whence using the cycle frequency we get,

qtot = ρ ec √ α T ∆S 2 √ t exch -24 √ 5 √ µρ air L 2 t 3/2 act 2(t exch + t act ) (19)
This new estimation of the cooling power takes into account both the reduction of heat exchanged with the cold reservoir (due to actuation losses) and the increase of the time period (due to the time of actuation) in a single criterion. Figures 4 and5 show a color map of qtot , in which the parameters are t act (related to the actuation losses) and t exch (related to the cold flux). In both figures the dashed lines represent iso-frequency curves where t exch + t act = f /2. Only positive qtot values are plotted while the white region corresponds to values of the (t act , t exch ) plane where qtot < 0. The used parameters are given in Table 1. These graphs show an interesting trend: for t act < t exch we see that increasing the frequency will also increase the power up to a maximum value. Once the maximum is reached, a further increase of the frequency produces significant losses up to the point where qtot < 0 (viz. white region).

Figure 4. Color map of the cold power per unit of surface function of t exch and tact using a log scale.. The losses were calculated using table 1. The white zone is when the total flux is negative, and thus it is 0 at the border.

In Fig. 4, the maximum of qtot , 183 mW cm -2 , is reached for an exchange time, t exch , of 2 ms, for actuation time, t act , of 1.3 ms. The resulting frequency is 154 Hz. The heat exchanged is 1.47 mJ cm -2 while the actuation losses are equal to 0.14 mJ cm -2 . At maximum power, the actuation time represents 40% of the time period, while the losses represents 19% of the heat transferred. Thus, to evaluate the performance of an electrocaloric cooling device, it is mandatory to study the actuation both in term of actuation time and of actuation losses.

Moreover, Fig. 5 shows that the reduction of air density by a factor 100 allows to increase the maximum to 325 mW cm -2 at a frequency of 486 Hz. In this case t exch and t act are respectively 0.63 ms and 0.40 ms, whereas q c and e tot,min are respectively 0.83 mJ cm -2 and 0.079 mJ cm -2 . When comparing the two graphs, we see that the reduction of air density by a factor 100 shifted the values by more than two decades to the right for t exch and less than one decade to the bottom for t act , while the values of the effective flux were nearly doubled near the maximum.

In spite of the actuation losses, increasing the frequency proves to be effective to increase the cooling power. This optimum could be further improved either by enhancing the electrocaloric material, i.e the T ∆S and the diffusivity α, or by improving the actuation time versus losses trade-off.

To improve the actuation losses, two parameters can be tweaked : L 2 and ρ 1/2 . While reducing the length L is very effective to improve the time versus energy loss trade-off, there is a limit at which we start loosing in effective area for the heat

DISCUSSION

Although the film has a cyclic motion, we assumed that the flow was stationary, and we neglected the transient states where the air has to be accelerated or decelerated. From an energetic angle, this is equivalent to the following statement: e visc k air . This approximation is valid for creeping flows where the Reynolds number is low and if the aspect ratio h/L is also low. However, we earlier sought the minimum of energy where the kinetic energy is balanced with the viscous losses; in this case the time derivative should be considered in the force balance equation.

Moreover, the flow was deemed incompressible. As mentioned by Gad-el-Hak [START_REF] Gad-El Hak | MEMS: Introduction and Fundamentals[END_REF], the general condition of incompressible flow is when the velocity of the flow is way less than the speed of sound. Since we want to actuate films of length of 5 cm in 1 ms, the velocity is 50 m s -1 which is sufficiently lower.

Furthermore, we considered that the kinetic energy is lost at each actuation. The actuation could be severely improved if we could store the inertia of the air during the exchange time and restore it for the next actuation, the same way we could recover the electrical charges from one side to the other when zipping.

Finally, to study the impact of the losses on the cooling power, we chose to model the losses as heat that is exchanged with the reservoirs. Without an experiment, it is however difficult to know how much goes to the cold or hot reservoirs. Thus, we chose to provide an estimation in the worst case scenario, that is when the viscous dissipation ends up heating only the cold reservoir, limiting the cooling power of the device.

CONCLUSION

Increasing the frequency is one of the major goals to enhance electrocaloric coolers so that they can compete with typical gasbased refrigerators. The model describes the major limiting factor towards higher frequency: the trade-off between actuation losses and actuation time. A simple Poiseuille model was used to show which parameters can be used to minimize the energy while reducing the actuation time. With this model, the distance between the hot and cold plate h can be optimized to get the minimum amount of losses. Doing so, we found two key parameters to improve the trade-off : the length of the film L, and the density of the air ρ air . While reducing the length is very effective because the losses diminish to L 2 , it is limited by the fact that the S-bent reduces the effective area for the heat exchange.

What is peculiar about this work lies in the combination of both the actuation and the heat exchange between the ECM and the reservoirs in a single figure of merit, which could be related to the cooling power depending on where the work of actuation will be dissipated.

Ultimately, lowering the pressure ρ air can reduce drastically the losses by reducing the kinetic energy of air that is lost if nothing is done to store it. This would enable to actuate within 1 ms at minimal losses, and thus to increase the frequency and the heat flux of electrocaloric coolers up to 486 Hz and 325 mW cm -2 respectively.

Figure 1 .

 1 Figure 1. A full cycle of an electrocaloric cooler using electrostatic actuation.The electrocaloric cooler device is based on a thin electrocaloric film (green), with a layer of metallization on both faces (dark gray lines). These are also used to actuate the film towards the electrodes (black lines) on each reservoir (light gray), using electrostatic forces (U zip ), the speed of the film being u f ilm . On each end, the film is stuck to one of the two electrodes with an insulation layer in-between (yellow), and the film and the electrodes are separated with a spacer (green) of thickness h. When the film is colored in orange, an electric field is applied to trigger the electrocaloric film (the film is heated and transfers Q hot to the hot reservoir), and when the film is blue, the field is removed (the film is cooled and receives Q cold from the cold reservoir). u air is the velocity field of the air, and E loss is the energy lost through one actuation.

Figure 2 .

 2 Figure 2. Schematic of the electrostatic actuator during the actuation step (b) in Fig. 1. The polymer film is in cyan, the electrical insulation layer is in yellow, u f ilm is the mean velocity of the film in the x direction. The S-bent region is located between s and s + δs, and the film is approximated by the dashed magenta line. The mass conservation under the incompressible assumption gives ∂u(x, y) ∂x = 0 (2)

Figure 3 .

 3 Figure 3. Losses due to actuation (blue line) depending on the actuation time compared to cold flux qc due to the electrocaloric effect for different film thickness/exchange time.

Figure 5 . 5 cm

 55 Figure 5. Colormap of the cold power per unit of surface for p = 10 mbar

Table 1 .

 1 Parameters

used for the graphs. The geometric parameter L is taken from Ma et al.

[START_REF] Ma | Highly efficient electrocaloric cooling with electrostatic actuation[END_REF]

.

ACKNOWLEDGMENTS

This work has benefited from the financial support of the LabEx LaSIPS (ANR-10-LABX-0032-LaSIPS) managed by the French National Research Agency under the "Investissements d'avenir" program (ANR-11-IDEX-0003).

BIBLIOGRAPHY