Christophe Giraud-Audine

Frédéric Giraud

Michel Amberg

Betty Lemaire-Semail

Generalized Modal Analysis for Closed Loop Piezoelectric Devices

Stress in piezoelectric material can be controlled by imposing the electrical field. Thanks to a feedback, this electrical field can be a function of some strain related measurement so as to confer to the piezoelectric device a closed loop macroscopic behaviour. In this paper we address the modelling of such system by extending the modal decomposition methods to account for the closed loop. To do so the boundary conditions are modified to include the electrical feedback circuit, hence allowing a closed-loop modal analysis. A case study is used to illustrate the theory and to validate it. The main advantage of the method is that design issue such as coupling factor of the device and closed loop stability are simultaneously captured.

Introduction

Piezoelectric devices are intrinsically suited for closed loop structure due to their reversibility: they can be actuators and sensors. It is therefore sensible to modify the dynamic of such a device by feeding back the strain state using the direct piezoelectric effect. Then, applying a voltage depending on this measured strain state, it is therefore possible for instance to confer a modified compliance to the system or increase dissipation in an assigned frequency band. This is the basic idea to many applications where frequency control or damping is required especially in the case of collocated structure. In [START_REF] Moheimani | On the feedback structure of wideband piezoelectric shunt damping systems[END_REF] or [START_REF] Fairbairn | A new approach to active q control of an atomic force microscope micro-cantilever operating in tapping mode[END_REF] for instance, a feedback model is used to modify the damping of structures with different strategies. As a prerequisite to apply the methods though, the frequency response of the mechanical structure must be identified. In [START_REF] Preumont | The damping of a truss structure with a piezoelectric transducer[END_REF], Preumont et al. address the problem of controlling the vibrations of large trusses. A quasi-static model of the actuator was considered, because in this case the dynamic of the piezoelectric actuator was much faster than the one of the controlled structure. This assumption cannot hold in some cases e.g MEMS or energy harvesters. Indeed, in this later cases, the overall dynamic must be accounted for since the resonance is the key issue [START_REF] Dutoit | Design considerations for mems-scale piezoelectric mechanical vibration energy harvesters[END_REF]Kim 2005, Lavrik, Sepaniak and[START_REF] Lavrik | Cantilever transducers as a platform for chemical and biological sensors[END_REF]. Moreover, it is also a well known fact that the coupling factor, is inherent to the overall feedback structure and the design of the actuator [START_REF] Preumont | On the damping of a piezoelectric truss[END_REF]. Therefore, the designer should consider a comprehensive approach in order to take into account the dynamic of the piezoelectric actuator in the case of the closed loop structure. Of course, the model is crucial in this context. Literature on the subject of piezoelectric actuators modelling is abundant, concentrating on fine quasistatic modeling [START_REF] Smits | The constituent equations of piezoelectric heterogeneous bimorphs[END_REF]Choi 1991, Goldfarb and[START_REF] Goldfarb | Modeling piezoelectric stack actuators for control of micromanipulation[END_REF] or dynamic operation near a specific resonance [START_REF] Mason | An electromechanical representation of a piezoelectric crystal used as a transducer[END_REF], Ballas 2007). When addressing broadband operations, one has to consider several modes [START_REF] Meirovitch | Fundamentals of Vibration[END_REF]. Therefore, modal decomposition is a natural and widely used tool [START_REF] Erturk | Assumed-modes modeling of piezoelectric energy harvesters: Euler-bernoulli, rayleigh, and timoshenko models with axial deformations[END_REF]Inman 2009, Ducarne, Thomas andDeü 2012). The approach consists in finding the modal shapes for a given electrical condition (open circuit or short circuit, the later being the most used) where the electric field or charge can be eliminated. These modal shapes are used to form a orthonormal basis on which any response can be decomposed (Meirovitch 2003, Geradin and[START_REF] Geradin | Mechanical Vibrations: Theory and Application to Structural Dynamics[END_REF]. Moreover, in some cases, the modal shape have analytical expressions that are helpful to deduce design rules (Ducarne et al. 2012). Even for more complex geometry, modal shapes of known simpler structures can be used to model the dynamic [START_REF] Erturk | Assumed-modes modeling of piezoelectric energy harvesters: Euler-bernoulli, rayleigh, and timoshenko models with axial deformations[END_REF]. However, since the initial analysis is performed in open loop, many modes may be necessary to model the closed loop response since the resonant frequencies may be modified. Thus, the insight can be lost, and in practice the accuracy of the design is dependent on the number of modes that are considered. In this paper, we show that the modal analysis method can be generalized to take into account the closed loop. The remainder of the paper continues as follows. Section 2 recalls the modeling of a clamped-free beam equipped with a sensor to provide feedback possibilities. Section 3 introduces the extension of the modal analysis in closed loop. We show that the boundary conditions can be modified to incorporate the effect of the feedback as moments dependent on the measurement. By properly dividing the contribution of the feedback induced moment into dissipative or generative (active power) and non-dissipative (reactive power) contributions, we prove that modal projection on a closed loop modal basis is possible. Section 4 presents experimental results to illustrate the approach and validate the prediction of a model based on the theory developed.

Modal analysis of a bender

Assumptions and notations

The schematic of Fig 1 depicts a simplified model of a slender beam used in this paper. The bender is mechanically excited at one end, and a feedback is used to modify its dynamic behaviour. In the following, the overall feedback circuit will be modelized as an amplifier and a series impedance so as to take into account simple filtering circuits. Regarding the mechanical modelling, the problem is supposed to be two dimensional in the sense that the different field are not dependent on x 2 . The bender is constituted by two piezoelectric layers polarized to operate on a 31 coupling. The end at x 1 = 0 is clamped to the moving frame while a proof mass M is clamped at the free end at x 1 = l. The moving frame is translated along the direction e 3 , its position relative to the reference frame is y(t). The bender is supposed to be mechanically unloaded, apart from the effect of inertia. The electrodes are connected as shown. In practice, as will be explained later, voltage V 2 is fixed while the voltage V 1 is varied. We consider small displacements and since the electrodes impose equipotentials, it will also be assumed that the electrical field is mainly along the e 3 direction. In the following, the voltage V 1 (t) is generated from the measure of the free tip of the bender w ,1 (L, t). Van-Voigt notation for tensors will be used in this paper. x ⊤ is the transpose of vector x. Matrices will be indicated by brackets e.g c E is the matrix of compliances at fixed electrical field. f ,ij denotes the partial derivative with respect to space variables x i ,x j , a repeated index denotes repeated partial derivative with respect to one of the space variable. Finally ḟ denote time derivative of the function f . 2.2. Model 2.2.1. Piezoelectric law for the following discussion, it is convenient to consider displacement and voltage to formulate the problem thus the relevant piezoelectric equations used throughout this paper are [START_REF] Ikeda | Fundamentals of Piezoelectricity[END_REF] :

T = c E S -[e] t E D = [e] S + ǫ S E (1)
with [e] the piezoelectric coefficient matrix and ǫ S the permittivity at fixed strain matrix. T, S are the stress and strain tensor, E, D are the electric and electrical displacement fields. The corresponding thermodynamic potential then writes [START_REF] Tiersten | Linear Piezoelectric Plate Vibrations Elements of the Linear Theory of Piezoelectricity and the Vibrations Piezoelectric Plates[END_REF] :

dG 2 (S, E) = TdS -DdE (2)
2.2.2. Kinematic of the model The bonding of the layers is supposed to be perfect, thus a continuous displacement field is considered. In this study shear strains are not considered the bender being thin and the frequencies being low. Thus a Euler-Bernouilli displacement field is used, given by (W. Weaver, Timoshenko and Young 1990, Geradin and Rixen 2014) :

u(x 1 , x 3 , t) =      u 1 (x 1 , x 3 , t) = u(x 1 , t) -x 3 w ,1 (x 1 , t) u 2 (x 1 , x 3 , t) = 0 u 3 (x 1 , x 3 , t) = w(x 1 , t) + y(t) (3)
Which results in the strain field :

S = {u ,1 (x 1 , t) -x 3 w ,1 (x 1 , t), 0, 0, 0, 0, 0} ⊤ (4)
2.2.3. Electrostatic solution Considering the electrical displacement D, the condition divD = 0 must hold. It is assumed that the layers are thin along e 3 compared to the dimensions in the other directions. Moreover, the potentials are imposed on the electrodes located on the horizontal sides of the bender, hence it is assumed that E 3 ≫ E 1 and E 3 ≫ E 2 . Let d ∈ {a, b} denotes the considered sub domain relative to the layers, then for each layer the electric field is given by E = -grad(φ d ). Using the (4) and (1), and neglecting the influence of the electric field along e 1 and e 2 , this condition reduces to :

D d3,3 = 0 ⇔ -e 31 w ,1 (x 1 , t) -ε 33 φ d (x 1 , x 3 , t) ,33 = 0 (5)
Integrating two times with respect to x 3 gives the expression of the electrical potential φ d in one layer of the bender which should write (Nadal, Giraud-Audine, Giraud, Amberg and Lemaire-Semail 2014) :

φ d (x 1 , x 3 , t) = - e 31 ǫ S 33 w ,11 (x 1 , t) x 2 3 2 + A d (x 1 , t)x 3 + B d (x 1 , t) (6)
The functions A d (x 1 , t),B d (x 1 , t), are functions to be determined from the equipotential conditions at the electrodes :

layer(a) φ a (-h) = 0 φ a (0) = V 1 layer(b) φ b (0) = V 1 φ b (h) = V 2 (7)

Dynamic equations

The augmented Lagrangian of the piezoelectric energy harvester is given by :

L = 1 2 ρ ut (x 1 , x 3 , t) u(x 1 , x 3 , t) -G 2 (x 1 , x 3 , t)dΩ+ 1 2 M u(l, 0, t) t u(l, 0, t) + V 1 Q 1 + V 2 Q 2 + F y (8)
where the four last contributions correspond respectively to the kinetic energy of the proof mass M, the energy supplied by the two generators, and the mechanical energy due to the external force applied for the mechanical excitation.

Applying variational calculus to the Lagrangian (8) leads to the equations of the problem [START_REF] Hammond | Energy Methods in Electromagnetism[END_REF]. The internal mechanical equilibrium involves the following equations for the extension and the flexion displacement field respectively :

-Mü(x 1 , t) + Ku ,11 (x 1 , t) = 0 (9a) -M ( ẅ(x 1 , t) + ÿ(x 1 , t)) + J ẅ,11 (x 1 , t) -Dw ,1111 (x 1 , t) = 0 (9b)
where :

• M is the mass per unit length

• J is the rotational inertia per unit length

• K is the compressive rigidity

• D is equivalent the flexion rigidity

The mechanical boundary conditions at the free end are :

Ku ,1 (l, t) = -N ex V 2 (t) (10a) Dw ,11 (l, t) = N f l V 1 - V 2 2 (10b) Dw ,111 (l, t) -J ẅ,1 (l, t) = M ẅ(l, t) + M ÿ(t) (10c)
where the electromechanical conversion coefficient for extension and bending are respectively :

N ex = e 31 a (11a)
N f l = e 31 ah (11b)
At the other end of the bender, the displacement are supposed to be perfectly imposed, hence :

u(0, t) = 0 (12a) w(0, t) = 0 (12b) w ,1 (0, t) = 0 (12c)
The electrical charges can be deduced from :

Q 1 (t) = C S (2V 1 -V 2 ) -N f l [w ,1 (x, t)] l 0 (13a) 2Q 2 (t) + Q 1 (t) = C S V 2 -N ex [u(x, t)] l 0 (13b) with C S = ǫ S 33 al h (13c)
C S is the clamped capacitance.

Theoretical study

The study will consider the case when V 2 is constant, which is necessary in practice to apply a compressive prestress in order to avoid excessive tensile stress during bending. Therefore, the displacement field u(x, t) is static since the extension and flexion fields are decoupled as can be deduced from the previous equations. Besides,without loss of generality, the frequency range considered is supposed to be sufficiently low to neglect the rotational inertia effects in (9b) and thus they will not be considered.

The equilibrium (9b) with the boundary conditions (10b), ( 10c), ( 12b) and (12c) are used for the modal analysis of the bender. To do so, the system is supposed to be free, hence every independent sources of the system are cancelled. This implies that for the mechanical side, y(t) must be set to zero. At a considered resonant angular frequency ω, the solution to (9b) is :

w(x 1 , t) = ψ(x 1 )η(t) (14)
where ψ(x 1 ) is the modal shape written as a linear combination of the Duncan functions (see appendix Appendix A) :

ψ(x 1 ) = As 1 β x 1 l + Bc 1 β x 1 l + Cs 2 β x 1 l + Dc 2 β x 1 l (15)
and η(t) is the vibration :

η(t) = η sin (ωt + α) (16)
Moreover the dispersion condition is :

Dβ 4 -Mω 2 l 4 = 0 (17)
The problem can be simplified further thanks to (12b) and (12c) :

A = B = 0
We proceed in the following by modifying the boundary condition (10b) to include the feedback effect.

Electrical equations

The electrical circuit that filters and amplifies the measurement signal is now introduced. Assuming a linear circuit, a Thévenin equivalent circuit is considered. We restrict the study to case when the voltage source is proportional to the measured signal ‡ :

v c (t) = G c w(l, t) (18)
where G c is the sensor gain possibly followed by an amplification. For the purpose of the study, we substitute 10b) and (13a). It reflects the fact that the bending moment is actually controlled by disturbing the middle plane voltage from the value it would normally have in a pure extension case. Since in the following harmonic oscillation will be studied, the various variables are expressed using complex vectors. Then, the displacement w(x 1 , t) writes § :

V 1 = V 2 2 + v in (
w(x 1 , t) → w(x 1 , t) = ψ(x 1 )η(t) = ψ(x 1 )ηe jωt
with η = ηe jα . Furthermore, equations will be written in a rotating frame so that the e jωt can be dropped, and for shorthand : W ,1... = ψ ,1... (L)η. ‡ A more general case involves a frequency dependent measure voltage. However, since there are no fundamental differences in the outline of the demonstration, this simpler case is considered. § x denote a complex number, and j = √ -1.

Let Z = N (jω) P (jω) be the equivalent series impedance of the Thévenin equivalent generator, the equations of the circuit are :

i = v c -v Z (19a) jC S ωv = i + jωN f l W ,1 (19b)
The second equation is deduced by time derivation of (13a). According to (13a) and (10b), v and v c can be expressed in terms of displacement or derivatives of displacement.

Writing that (19b) and (10b) can be combined, resulting in the following equation :

P (jω) = p k (jω) k and N(jω) = n k (jω) k k ∈ N + , equations (19a),
  deg(P ) k=0 C S p k (jω) k+1 + deg(N ) k=0 n k (jω) k   DW ,11 + κ 2 f l deg(P ) k=0 p k (jω) k+1 W ,1 -κ f l G c deg(N ) k=0 n k (jw) k W = 0 (20)
For the sake of clarity, the condition n 0 = 0 is assumed. In this case, the previous result can be reformulated as :

DW ,11 = R(jω)W ,11 + S(jω)W ,1 + T (jω)W (21)
where R(jω) = -

deg(P ) k=0 C s p k n 0 (jω) k+1 - deg(N ) k=1 n k n 0 (jω) k , S(jω) = - κ 2 f l
n 0 P (jω) and T (jω) = κ f l n 0 N(jω) are polynomials. This defines the new boundary condition for the bending moment including the closed loop. For the following discussion, rewrite the previous equation :

DW ,11 = R d (jω)W ,11 + S d (jω)W ,1 + T d (jω)W + jω R q (jω)W ,11 + S q (jω)W ,1 + T q (jω)W (22)
where R d (jω), R q (jω), S d (jω), S q (jω), T d (jω) and T q (jω) are real even polynomials of ω.

Only the non-dissipative contribution of this condition must be considered for modal analysis. Multiplying both sides of (20) by the the conjugate of the rotation speed at the end of the tip (jωW ,1 ) * , gives the complex power transmitted to the bender by the feedback :

-jωDW ,11 W * ,1 = -jω R(jω)W ,11 + S(jω)W ,1 + T (jω)W W * ,1
Two different contributions should be distinguished :

• imaginary terms, corresponding to real projection on the real axis of the closed loop bending moment complex vector (therefore leading or lagging in quadrature with the rotation speed), which are conservative ;

• real terms, corresponding to imaginary part of the closed loop bending moment complex vector, which are either dissipative or supplying power.

Since W , W ,1 and W ,11 are in phase, it can be concluded that real parts of R, S and T (denoted R d , S d and T d in the sequel) will contribute to the reactive power, while imaginary parts (denoted R q , S q , T q ) will be responsible for active power .

Secular equation

For the modal analysis, only the conservative contribution of the system must be considered. Thus, according to the previous discussion, only the R d , S d and T d are relevant, and thus, after cancelling the η that appear on both sides of the equality sign, the boundary conditions ( 10b) and (10c) are rewritten as :

Dψ ,11 (l) = R d (ω)ψ ,11 (l) + S d (ω)ψ ,1 (l) + T d (ω)ψ(l) (23a) Dψ ,111 (l) = -Mω 2 ψ(l) (23b)
Using ( 17) to rewrite these conditions in β, a general secular equation can be obtained (see Appendix B). Solving this equation for β = β k gives the resonant frequencies ω k using ( 17), and once replaced in ( 14), the mode shape ψ k can be obtained up to a multiplicative constant .

In the usual modal analysis method, the modes are used to decompose the solution to a given excitation. These mode shapes which already partly include the electrical feedback are now compared to the classical mode shapes that would arise from considering short circuit conditions.

Properties of the closed-loop modal shapes

The closed loop (C.L.) modal shapes can be used to describe the solution to the forced vibrations problem. Indeed, let the θ j (x 1 ) be the normalized open loop (O.L.) modal shapes of the bender, obtained by replacing the feedback circuit by a short circuit which verify :

-Mν 2 j θ 2 j + Dθ j,1111 = 0 (24)

where ν k are the corresponding modal angular frequencies, for the boundary conditions :

θ j (0) = 0 (25) θ j,1 (0) = 0 (26) Dθ j,11 (l) = 0 (27) Dθ j,111 (l) = -Mν 2 j θ(l) (28) (29)
The normalization considered is defined by [START_REF] Erturk | An experimentally validated bimorph cantilever model for piezoelectric energy harvesting from base excitations[END_REF] :

l 0 Mθ 2 j (x 1 )dx 1 + Mθ 2 j (l) = 1 (30)
Note that the same normalization will be applied to the ψ k . The O.L. modal shapes verify :

l 0 Mθ j (x 1 )θ k (x 1 )dx 1 + Mθ k θ j = δ jk (31)
where δ jk = 1 if j = k and δ jk = 0 if j = k. Therefore, the C.L. modal shapes ψ k (x 1 ) can be expressed as :

ψ k (x 1 ) = ∞ j=1 α kj θ j (x 1 ) α kj ∈ R (32)
Since the ψ k (x 1 ) depend on the wave numbers {β k } which form and increasing sequence, it is clear from the expression (15) that :

ψ k (x 1 ) = ψ j (x 1 ) if k = j (33)
It follows that the ψ k are independent vectors that can be related to a free basis, and thus can be used also as a basis, that is any response of the bender can be written :

w(x 1 , t) = ∞ k=1 ψ k (x 1 )η k (t) (34)
The standard procedure to solve the vibration problem involves the projection of the dynamic equilibrium verified by a modal shape ψ k on a modal shape ψ j is :

l 0 [Mψ k (x 1 )η k (t) + Dψ k,1111 (x 1 )η k (t)] ψ j (t)dx 1 = 0 (35)
Using part integration, and the boundary conditions (23a) and (23b), the previous equation for harmonic vibrations (i.e ηk (t) = -ω 2 k η k (t)) is rewritten as follows :

-

ω 2 k l 0 Mψ k (x 1 )ψ j (x 1 )dx 1 + Mψ k (l)ψ j (l) -[R d (ω k )ψ k,11 (l) + S d (ω k )ψ k,1 (l) + T d (ω k )ψ k (l)] ψ j,1 (l) + D l 0 ψ k,11 (x 1 )ψ j,11 (x 1 )dx 1 = 0 (36)
Considering the special case where j = k, and using the proposed normalization given by (30), the potential energy for the modal shape, including the feedback, verifies :

-[R d (ω j )ψ j,11 (l) + S d (ω j )ψ j,1 (l) + T d (ω j )ψ j (l)] ψ j,1 (l) + D l 0 ψ j,11 (x 1 ) 2 dx 1 = ω 2 j (37)
The projection of the dynamic equilibrium of mode ψ j on mode ψ k leads to an equation similar to (36) where the k and j indexes are swapped. The difference of these equations would give :

(ω 2 j -ω 2 k ) l 0 Mψ k (x 1 )ψ j (x 1 )dx 1 + Mψ k (l)ψ j (l) + [R d (ω j )ψ j,11 (l) + S d (ω j )ψ j,1 (l) + T d (ω j )ψ j (l)] ψ k,1 (l) -[R d (ω k )ψ k,11 (l) + S d (ω k )ψ k,1 (l) + T d (ω k )ψ k (l)] ψ j,1 (l) = 0 (38)
Therefore, the C.L. modal shapes are generally not orthogonal in closed-loop for the functional (35) . Actually, in C.L. the functional should be modified to account for the electrical equilibrium of the bender and the circuit. Fortunately, this tedious task can be avoided to solve the forced response due to the feature discussed below.

3.3.1. Response to the forced vibration Indeed, the problem can still be decomposed.

Assuming that the solution can be written as in (34) the projection of the solution on mode ψ j writes :

l 0 ∞ k=1 Mψ k (x 1 )η k (t) + Dψ k,1111 (x 1 )η k (t) + mÿ ψ j (x 1 )dx 1 = 0 (39)
Isolating the contribution of the j th mode and performing part integration gives :

l 0 ∞ k =j Mψ k (x 1 )η k (t) + Dψ k,1111 (x 1 )η k (t) ψ j (x 1 )dx 1 + M l 0 ψ j (x 1 ) 2 dx 1 ηj (t) + D l 0 ψ j,11 (x 1 ) 2 dx 1 η j (t) -[R d (ω)ψ j,11 (l) + S d (ω)ψ j,1 (l) + T d (ω)ψ j (l)] ψ j,1 (l)η j (t) + ÿ(t) l 0 Mψ j (x 1 )dx 1 + Mψ j (l) 2 ηj (t) + M ÿ(t)ψ j (l) = 0 (40)
First integral in this expression is nil because modal shape verify (35), thus it can be deduced using ( 30) and ( 37) that the equation simplifies to :

ηj (t) + ω 2 j η j (t) = -φ j (t) (41)
with the modal inertia forces :

φ j (t) = l 0 Mψ j (x 1 )dx 1 + Mψ j (l) ÿ(t) (42)
These equations allow to calculate the modal shapes, the resonant frequencies as well as the coupling factors. However, to be able to predict the vibration amplitude, it is necessary to include the active effects of the electrical circuit.

This depends on the feedback law : for instance, it can be easily deduced from (38) if the feedback is proportional to the derivative of the displacement w ,1 (l, t), orthogonality property is retrieved for Z = 1 21) and ( 22) show that the odd powers of the polynomials R(jω), S(jω) and T (jω) (that is R q , S q and T q ) can either dissipate or provide power. The procedure is the same, that is assuming the solution (34), the dynamic equation is projected onto a mode ψ j . However, the general boundary condition for moment is now considered :

Dw ,11 (l, t) = R d w ,11 (l, t) + S d w ,1 (l, t) + T d w(l, t)
R q ẇ,11 (l, t) + S q ẇ,1 (l, t) + T q ẇ(l, t) (43) where R d , S d , T d , R q , S q , and T q are the linear differential operators with respect to time corresponding to the ω polynomials defined in ( 22). The modal equations are then :

ηj (t) + 2ξω j ηj (t) + ω 2 j η j (t) = -φ j (t) - ∞ k=1 φ jk (t) (44)
In the previous equation the modal damping ξ has been introduced [START_REF] Geradin | Mechanical Vibrations: Theory and Application to Structural Dynamics[END_REF]. A cross coupling between modes appears with the contributions φ jk (t) given by :

φ jk (t) = [R q ηk (t)ψ k,11 (l) + S q ηk (t)ψ k,1 (l) + T q ηk (t)ψ k (l)] ψ j (l) (45)

Experimental validation

Experimental setup

A picture of the experimental set up is reproduced in Fig. 2. The co-fired multilayer bender is a Noliac CMBP05 (NOLIAC 2011). The material is NCE57. The main specifications are summed up in tab. 1, and were used for the model that is discussed hereafter. The bender is clamped on one side fixed on a vibrating pot. On the other side, a magnet is fixed (weight 0.8 g) and for some test a tungsten mass can be fixed. The magnet is used to measure the end tip displacement thanks to a Hall sensor. This sensor was calibrated using a laser sensor (Polytec OFV 505), the sensitivity was estimated to 8.3 mV/ñm. The feedback is realized by amplifying the displacement measure using operational amplifier to select the gain by combining resistors then a power amplifier (NF HSA 4051) is used to post amplify with gains of 20 and 40, resulting in supply voltages varying from 0 to 250 V. An accelerometer is fixed on the moving frame. In order to establish the frequency responses of displacement, a dynamic signal analyser (Stanford Research SR 785) is used. Finally, electrical measurements (voltage, current, acceleration and displacement sensors output) are realised using an oscilloscope (Tektronik TDS 3014).

4.2. Case study : effect of a position feedback and output series resistance on the closed-loop resonance 4.2.1. Model For practical reasons, a series resistance was implemented. Indeed, sensor noise are amplified in the closed loop exciting higher resonances, and can even become destabilising. To prevent this, the resistance was initially introduced to constitute in combination with the capacitance C S low pass filter.

To address the effects of this resistance, since Z = R s , the following polynomials are considered :

N = R s (46a) P = 1 (46b)
Then the condition (43)is : with τ = R s C S . So, for this simple case :

w ,11 (l, t) = N f l D G c w(l, t) -τ ẇ,11 (l, t) - N 2 f l DR ẇ,1 (l, t) (47)
     R d w ,11 (l, t) = 0 S d w ,1 (l, t) = 0 T d w(l, t) = N f l G c w(l, t)
and

     R q ẇ,11 (l, t) = -τ D ẇ,11 (l, t) S q ẇ,1 (l, t) = - N 2 f l
Rs ẇ,1 (l, t) T q ẇ(l, t) = 0 4.2.2. Mechanical resonance As far as the effect of the feedback on resonance is concerned, only T d w(l, t) has to be considered. For the sake of generality, the variable change ζ = x l is introduced. The plot of β against the normalized gain k =

N f l Gcl 2 D
(fig. 3) then depicts the effect of this feedback on the modal shapes (see fig 4 for the case of the first mode), and thus on the resonance frequencies for various mass ratio µ = M M : • k > 0 reduces β, that is it decreases resonance frequencies, thus it can be assimilated to a softening in the sense that it acts as if the rigidity of the bender was less than in open-loop. The effect is rather drastic as for k = 1 the first mode will become unstable due to buckling.

• On the contrary, k < 0 increases β, although rather progressively, so a hardening is induced by the closed loop. Asymptotic analysis shows that for the first mode lim k→∞ β = π

The theoretical (solid and dashed lines) and measured (cross and dots) wave number are compared for some values of the normalized gain k on fig. 5. Theoretical -40, -26, -20, -13, 0, 13, 20, 26, 40}. Taking into account the gain of the sensor and the flexion rigidity, the reduced gain were k = {-1.13, -0.74, -0.57, -0.37, 0, 0.37, 0.57, 0.74, 1.13}.The experimental values of β are deduced from the measured resonance frequencies. Two values of the tip mass to bender mass ratio µ were considered (µ = {0.06, 0.54}). The experimental results are in good agreement with the theoretical prediction. The effects of the closed loop on the dissipation are now examined. Due to the variation of β with respect to k, the modal shape ψ j are modified, and therefore so are the modal forces φ j (42). This is depicted on fig. 6 (top) for the two first modes. Moreover, they affect the modal damping as can be seen in relations ( 44) and (45). In this case, assuming that the cross coupling of the modes is negligible near the mode resonance since the considered mode dominates, it is possible to express the closed loop damping factor as follows :

ξ ′ j = ξ j -1 2 (r q ψ j,11 (l) + s q ψ j,1 (l)) ψ j,1 (l)

where :

r q = - τ D l 3 and s q = Rκ 2 f l l 2
The curves of the damping versus k for the case µ = 0.06 are presented on fig. 6 and it can be seen that the feedback has a clear impact especially on the first mode. Note that to obtain these curves. The O.L. damping was estimated form the O.L. quality factor on the frequency response. For both modes, the damping vanishes as the gain tend to low values, which indicate a possible instability. The predicted resonance curve of the tip vibration amplitude and phase shift at imposed acceleration are depicted on fig. 7. The model prediction are in good agreement with the experiment as for the resonance amplitudes and frequencies. Still, some discrepancies can be found as the experimental resonance curves are non symmetric, especially for low values of damping. Tests showed that the resonance curves are different for increasing and decreasing frequency sweeping 4.2.3. Current at resonance Finally, fig. 9 compares the measured and predicted currents delivered by the amplifier to the piezoelectric device at the mechanical resonance frequencies. The theoretical current is obtained by replacing the expression of the derivative of the displacement into (13a) and performing a time derivation :

i 1 (t) = Q1 (t) = 2C S V1 -κ f l ẇ,1 (L, t) motional current
Again, the results are in good agreement. Additionally, the motional current is represented to illustrate the effect of variation of β due to k that results in an inflexion of the motional current curve when the resonant frequency increases.

Conclusion

This paper proposes an alternative modal decomposition so as to study the effects of feedback on the dynamic of a piezoelectric device. The method is general and consists in the following steps :

• the boundary conditions of the dynamic equation of the system are modified to include directly the feedback as written in (21). Two contributions can be distinguished : reactive (non dissipative) terms depending on even time derivatives of the displacement or spatial derivative of the displacement, active terms depending on odd time derivative of the same functions.

• a modal analysis, consisting in solving the partial derivative equations when the dissipative terms are cancelled. This results in a set of modal shapes (closed loop modal shapes) that presents some properties similar to the classical modal decomposition : the equation of motion can still be projected to form a infinite system of independent equations depending solely on time.

• finally, the dissipative contribution are included as external sources and the solution is then written as a decomposition on the closed loop modal basis. The method has been applied to a simple set up that allows to derive the solutions and has been validated against experimental results. The benefit of the approach is to consider the closed loop system as a whole system and represent the electromechanical interaction including the control. The classical modal decomposition, obtained by considering the short circuited resonance, can also be used to obtain this but the truncation of the projected solution is still the key to the accuracy of the model near a resonance of the closed loop system. In contrast, the proposed method can efficiently model the dynamic with a few modes. Moreover, some issues such as stability or energy consumption can be directly addressed. Therefore, the expected benefits of the method are in the field of design or control of devices where the decoupling of the mechanical structure dynamic and the piezoelectric device is not verified e.g energy harvesters such as the one used here, or atomic force microscope tips.

Figure 1 .

 1 Figure 1. Schematic of the system

Figure 2 .

 2 Figure 2. Experimental setup : (1) bender with a magnet and proof mass fixed at the end; (2) Hall sensor; (3) shaker; (4) accelerometer

Figure 3 .Figure 4 .

 34 Figure 3. Theoretical effect of the normalised feedback gain k on the wave number β

Figure 5 .

 5 Figure 5. Theoretical (th.) and measured (exp.) wave number β vs normalised gain k

Figure 6 .

 6 Figure 6. Theoretical dependence of the modal forces (top) and of the closed loop damping coefficients (bottom) on the normalised gain k for different values of µ

Figure 7 .

 7 Figure 7. Comparison of the theoretical (plain lines) and measured (dotted line) of the dynamic mechanical gain (top) and phase (bottom) on the normalised gain k (µ = 0.06)

Figure 8 .Figure 9 .

 89 Figure 8. Measured and predicted mechanical gains at mechanical resonance frequencies (top) and comparison of the model predictions when including the closed loop damping introduced by the electrical circuit (dashed) or not (plain) (bottom graph)

Table 1 .

 1 Characteristics of NCE 57 piezoelectric material

	ǫ T 33 ǫ 0	1800	-	Relative dielectric constant
	tanδ 170.10 -4	-	Dielectric loss factor
	-d 31 170.10 -12 C/N k 31 0.33 -	Piezoelectric charge constant Electromechanical coupling factor
	ρ	7.7.10 3 kg.m 3	density
	s E 11	17.10 -12 m 2 /N	elastic compliances
	Q	70	-	Mechanical quality factor

Acknowledgments

This work has been carried out within the framework of the Stimtac project of IRCICA (Institut de Recherche sur les composants logiciels et matriels pour l'information et la communication avance), and the Mint project of Inria. This work is also supported by SEEDS and the FUI Touchit.

Appendix A. Duncan functions

Duncan functions are defined as :

with i ∈ {1, 2} They verify the property :

Expression of the secular equation

Defining :

the secular equation is obtained after some calculus :


C Giraud-Audine
email: christophe.giraud-audine@ensam.eu

F Giraud

M Amberg

Influence of feedback on modal shapes and its consequences on broadband energy harvesting

Keywords: Piezoelectric, broadband energy harvesting, modal analysis, I

We study the effect of feedback on an energy harvester for broadband applications. Using modal analysis we show that it has an influence on modal shape as a side effect of the resonance shifting. We show that the piezoelectric coupling can thus be reduced. Experimental results are in agreement and also reveal that non-linear effects are promoted in some cases.

INTRODUCTION

Energy harvesting from ambient vibrations is considered as one solution to provide energy to autonomous system requiring low power or working intermittently [START_REF] Dutoit | Design considerations for MEMS scale piezoelectric mechanical vibration energy harvester[END_REF]. Resonance is the optimal condition for energy harvesting because elastic stress is compensated by inertia forces. However, in practice the resonance condition is rarely verified. First, the mechanical source usually generates vibrations in a bandwidth rather than a single frequency during operation. Second, the parameters of the harvester can vary resulting in a shift of its own resonance frequency, hence the ongoing research on broadband harvesters [START_REF] Twiefel | Survey on broadband techniques for vibration energy harvesting[END_REF].

Another possibility is to use the inverse piezoelectric effect to compensate the elastic stress or the inertia [START_REF] Roundy | Toward self-tuning adaptive vibration-based microgenerators[END_REF][START_REF] Lallart | Frequency Self-tuning Scheme for Broadband Vibration Energy Harvesting[END_REF]. The resonance can then be shifted thus enabling to harvest energy on a larger bandwidth. The contribution of this paper is to show that proportional feedback using the measurement of the tip displacement can modify both the frequencies and the modal shape. However, we show that the coupling factor of the harvester will also be impacted.

The paper is structured as follow: in next section, the model is briefly presented. In the third section, modal bimorph bender with feedback is analyzed. Finally the model is compared to experimental tests.

II. MODEL

A. Assumption and notations

The schematic of Fig 1 depicts a simplified model of a slender bender used for energy harvesting. The problem is supposed to be bidimensionnal that is the different fields are Figure 1. schematic of the studied system not dependent on x 2 (direction orthogonal to the plane of the schematic).The bender is constituted by two piezoelectric (a and b) layers polarized to operate in 31 mode. The end at x 1 =0 is clamped to the moving frame while a proof mass M is clamped at the free end (x 1 =l). The moving frame is translated along the direction e 3 , its position relative to the reference frame is y(t). The bender is supposed to be mechanically unloaded, apart from the effect of inertia.

The electrodes are connected as shown. The voltage V 2 is constant while the voltage v is varied accordingly to some measurement of the end tip. We consider small displacements and since electric potentials are imposed on the electrodes, it will also be assumed that the electrical field is mainly along the e 3 direction.

B. Piezoelectric

For the following discussion, it is convenient to consider displacement and voltage to formulate the problem thus the relevant piezoelectric equation used throughout this paper are :

E ε eS D E e S c T S Τ E     (1)
c E is the rigidity matrix at constant electric field, e is the piezoelectric matrix, and  S is the permittivity at constant strain. The corresponding thermodynamic potential is [START_REF] Ballas | Piezoelectric multilayer beam bending actuators static and dynamic behavior and aspects of sensor integration[END_REF] :

E D S T E S, d - d ) ( dG 2    (2)
This work has been carried out within the framework of the Mint project of Inria (France) and the StimTac Project of the Institut de Recherche sur les Composantans logIciels et matériels pour la Communication Avancée (IRCICA)

C. Kinematic

The bonding of the layers is supposed to be perfect, thus a continuous displacement field is considered. In this study shear strain is not considered, the bender being thin. The deformation considered consists in the superposition of an extension strain and a flexion. Thus, following the classical Bernouilli beam theory [START_REF] Ballas | Piezoelectric multilayer beam bending actuators static and dynamic behavior and aspects of sensor integration[END_REF][START_REF] Géradin | Mechanical vibrations: theory and application to structural dynamics[END_REF], the displacement is given by:

      ) ( ) , ( ), , ( , ) , ( ) , , ( 1 1 1
3 1 3 1 t y t x w t x w x t x u t x x u (3)
Which results in a strain field:

     0 , 0 , 0 , 0 , 0 , 0 ), , ( ) , ( ) , , ( 1 11 , 3 1 1 , 3 1 t x w x t x u t x x S (4)

D. Electrostatic

Considering now the electrical displacement D, the condition divD=0 must hold. According to the piezoelectric law (1) and the strain (4), and neglecting the component of the electric field along e 1 , the electrical potential in one the piezoelectric of the bender is given by:

) , ( ) , ( 2 ) , ( ) , , ( 1 3 1 2 3 1 11 , 33 31 3 1 t x B x t x A x t x w e t x x i i s i       (5)
The functions A i (x 1 ,t), B i (x 1 ,t), where i {a,b} denote the two piezoelectric domains, are functions to be determined from the potential on the electrodes.

E. Dynamic Equations

The augmented Lagrangian of the piezoelectric energy harvester is given by:

  2 2 2 1 2 2 2 2 1 2 1 1 1 Q V Q Q v V Fy M d G L l x l x                      u u u u      (6)
where the four last contributions corresponds to the mechanical energy due to the external force, the kinetic energy of the proof mass and the energy supplied by the generators. Applying variation calculus, the following equations for the extension and the flexion displacement field respectively are obtained:

0 ) ( ) , ( ) , ( ) , ( 0 ) , ( ) , ( 1 1111 , 1 11 , 1 1 11 , 1         t y t x Dw t x w J t x w t x Ku t x u            (7)
where,  is the mass per unit length, J is the rotational inertia per unit length, K is the compressive rigidity, D is equivalent the flexural rigidity. w ,i denotes the derivative with respect to x i , dots are time derivatives. The mechanical boundary conditions at the free end are:

             ) ( ) , ( ) , ( ) , ( ) , ( ) , ( ) ( ) , (

III. THEORETICAL STUDY

A. Secular Equation

In practice V 2 is constant to apply a pre-stress during bending. Therefore, u(x 1 ,t) is constant and will not be considered in the following. Besides, the rotational inertia effects in [START_REF] Erturk | An experimentally validated bimorph cantilever model for piezoelectric energy harvesting from base excitations[END_REF] will not be considered in this study. To isolate the effect of feedback, M is not considered. The equilibrium equation ( 7) with boundary conditions (8), and (9) are used to form the set of conditions that will lead to the modal analysis of the bender. Being an external excitation, y(t) is set to zero [START_REF] Géradin | Mechanical vibrations: theory and application to structural dynamics[END_REF]. Measurement of the displacement is realized by a sensing ceramic followed by an integrator. The resulting voltage is proportional to the deformation at the end of the bender. To realize a feedback, the applied voltage is proportional to the measurement voltage:

B. Resonant Frequencies and Modal Shapes

Eq. ( 11) means that the values of  can be modified by the gain k as shown on Fig. 2. The curves present asymptotes for  such that tantanh between the modes. A softening or a hardening effect depending on the sign of k can be obtained. Indeed, considering the first mode which is of special interest for energy harvesting, the wave number will tend toward 2.635 for an infinitely negative, so the ratio between the clamped free resonance frequency and the closed loop resonance frequency is limited to 1.588.

A positive gain is more effective. Indeed, a variation form 0 to 1 would theoretically decrease  from 1.875 to 0. For k=1 a bifurcation occurs. The eigenvectors can be determined from ( 8) to (10), and the solutions  k to (11) up to a multiplicative constant. In the case at hand they write:

                     l x l x l x l x A x W k k k k k k k k k k         sin sinh cosh cos sinh sin cos cosh ) ( 1 (13)
The modes are orthogonal and can be normalized. For increasing positive gain, the curvature becomes more pronounced, but the mode amplitude decreases. This is confirmed by Fig. 4 where the modal mass are represented vs the feedback gain.

C. Modal Force

As mentioned, the normalized modes form an orthonormal base. Eq [START_REF] Erturk | An experimentally validated bimorph cantilever model for piezoelectric energy harvesting from base excitations[END_REF] can be solved by projection on each mode. The forces are then transformed to modal forces. The mechanical load y  μ in ( 7) is uniform, thus, the corresponding modal force is proportional to:

  1 1 ) ( dx x W k k  (15)
The amplitude are likely to be larger for high participation factor. Fig. 4 shows that for softening gains, the participation factor collapses, while for hardening gain, it remains relatively constant. This is similar to the evolution of the modal mass, presented on the same graph.

D. Modal Electromechanical Coupling

Using the piezoelectric equations (1), the strain field (4) and the electrical potential (5), the charges on the electrodes are obtained by integrating D on their surface. One obtains for Q 1 : /   versus the feedback gain and the wave number. This quantity is proportional to the modal electromechanical coupling [START_REF] Erturk | An experimentally validated bimorph cantilever model for piezoelectric energy harvesting from base excitations[END_REF]. On the softening side, the coupling first slightly increasing (due to the increased curvature) then rapidly collapses (as the modal mass decreases). On the hardening side, the effect is similar, because the moment blocks the free end rotation. This in turn degrades the electromechanical coupling as it depends on the difference of slopes at both ends.

  l fl s t x w V h al t Q 0 1 1 , 2 33 1

B. Results

The result for the magnitude (in dB) between the vibrating pot's acceleration and the speed of the bender's tip is presented on Fig. 6 for various gains. The variation range of the gain was imposed by the voltage limitation for the hardening side (negatives gains) and instability of the closed loop on the softening side (positive gains) due to the bifurcation of the wave number as k approaches 1. As a consequence, the closed loop gain's sign changes with the inversion on the free end slope inducing a positive feedback. The resonant frequency can be shifted from 225 Hz to 410 Hz. The short circuit resonance is obtained by setting G=0 (f cc =329 Hz), therefore the percentage of variation (f cl -f cc )/f cc is -32% to 25%. For high positive values of k, non linearities become prominent because the feedback compensates the linear elasticity of the material.

The shapes were measured for some frequencies and compared to the predicted ones using (14) on Fig. 7. Since the amplitude is irrelevant, the graphs are scaled. The agreement between measurement and prediction confirm that the modal analysis is correct. To further assert the model, the measured and theoretical wave numbers are compared on Fig. 8. Again, a good agreement can be observed.

V. CONCLUSION

An analysis of a bender using a feedback for broadband harvesting has been proposed. The shifting of the resonant frequency is indeed possible within a range that was identified. However, it was shown that the potential benefit is in fact limited even though losses have not been considered in the study. Regarding softening, the mode are less easily excited by a uniform load that would result from the acceleration, this suggests a different shape of the bender for instance a varying section. Regarding hardening, as the bending moment tends to block the rotation at the free end, the coupling is degraded. Tests have confirmed these conclusions, although the electrical coupling could not be directly investigated due to the difficulty to measure the motional current. This point will be addressed in future works.

Figure 2 .Figure 3 .

 23 Figure 2. Wave number  versus global feedback gain k

Fig. 3

 3 Fig.3represents mode shapes according to the wave number and the corresponding feedback gain. For negative gain, the shapes tend toward a clamped-blocked rotation shape.

Figure 4 .Figure 5 .

 45 Figure 4. Participation factor and modal vs gain k

  can be decomposed into a charge involving the blocked capacitance, and a charge induced by the strain in the material Q m which is potentially harvested. On Fig.5the graphs represent the reduced induced charge

Figure 6 .Figure 7 .

 67 Figure 6. resonance for various experimental feedback gain


R Chabiniok

E Sammut

M Hadjicharalambous

L Asner

D Nordsletten

R Razavi

N Smith

Steps towards quantification of the cardiological stress exam

In this work we aim to advance the translation of modelbased myocardial contractility estimation to the clinical problem of quantitative assessment of the dobutamine stress exam. In particular, we address the question of limited spatial resolution of the observations obtained from cine MRI during the stress test, in which typically only a small number of cine MRI slices are acquired. Due to the relative risk during the dobutamine infusion, a safe acquisition protocol with a healthy volunteer under the infusion of a beta-blocker is applied in order to get a better insight into the contractility estimation using such a type of clinical data. The estimator is compared for three types of observations, namely the processed short axis cine stack contiguously covering the ventricles, the short axis stack limited to only 3 slices and the combination of 3 short and 3 long axis slices. A decrease of contractilities in AHA regions under the beta-blocker infusion was estimated for each observation. The corrected model (by using the estimated parameters) was then compared with the displacements extracted from 3D tagged MRI.

Introduction

A dobutamine stress test is used clinically to identify regional wall motion abnormalities to guide management on therapeutic options, such as revascularisation [START_REF] Jahnke | Prognostic value of cardiac magnetic resonance stress tests: Adenosine stress perfusion and dobutamine stress wall motion imaging[END_REF]. It is performed at increasing doses of dobutamine, an inotropic drug, with simultaneous imaging by echocardiography or magnetic resonance imaging (MRI) at each dosage level. Significant risks such as arrhythmia and cardiac arrest at higher doses, and the unpleasant sensation of the pharmacologically induced stress requires a rapid data acquisition which should not exceed three minutes at each dobutamine dosage. This limits the spatio-temporal resolution of the obtained image data. In particular in cardiac MRI -the modality used in this paper -a typical clinical protocol includes an acquisition of 3 slices in short axis placed in basal, mid-and apical third of left ventricle (LV), see Fig. 1 (middle), and 3 standard long axis images of LV (the so-called 2-, 3-and 4-chamber views). In addition to the limited spatial coverage, the demanding breath-holding under the dobutamine infusion often compromises the quality of image data.

A positive test is defined by a regional wall motion abnormality during stress and is currently clinically assessed entirely qualitatively. An estimation of some constitutive parameters in a biomechanical model is a way of assessing mechanical properties of myocardium [START_REF] Chapelle | Fundamental principles of data assimilation underlying the Verdandi library: applications to biophysical model personalization within euHeart[END_REF][START_REF] Xi | The estimation of patient-specific cardiac diastolic functions from clinical measurements[END_REF][START_REF] Wang | Modelling passive diastolic mechanics with quantitative MRI of cardiac structure and function[END_REF], and in [START_REF] Chabiniok | Estimation of tissue contractility from cardiac cine-MRI using a biomechanical heart model[END_REF] quantitative values directly targeting regional myocardial contractility were estimated using cine MRI as observations. Such an approach has a potential to provide a higher reproducibility of the exam, lower inter-observer variability and possibly may even allow a decrease in the dose of dobutamine as the quantitative measure may increase the sensitivity. However, in comparison to [START_REF] Chabiniok | Estimation of tissue contractility from cardiac cine-MRI using a biomechanical heart model[END_REF] -where a full short axis cine stack covering whole ventricles was used for contractility estimation -in the stress test we need to deal with the limited spatio-temporal coverage. The presented study focuses on exploring the feasibility of translating the methodology of [START_REF] Chabiniok | Estimation of tissue contractility from cardiac cine-MRI using a biomechanical heart model[END_REF] into the clinical setup of the dobutamine stress exam with respect to the challenge of limited coverage.

Due to the associated risks in using dobutamine, it would be impossible to perform trial acquisitions in a healthy volunteer or significantly prolong the patient exam to acquire additional short axis slices under the dobutamine stimulation. Therefore, to get an insight into various image sampling modes, we estimate a change of contractility in a healthy volunteer under the safe pharmacological stimulus of the beta-blocker esmolol. This allows safe acquisition of high quality image data of a full cine stack, which can be then downsampled in space or time, in accordance to the real clinical data of the dobutamine stress test. As opposed to the positively inotropic drug dobutamine, the negatively inotropic beta-blocker globally reduces the contractility. This is reflected in the decreased stroke volume (see Fig. 1, right) and we expect that the estimated contractility values should follow a similar trend.

After describing the experimental data together with image processing in Section 2, and the modeling and parameter estimation framework (Section 3), we present in Section 4 the results of contractility estimation at baseline and under beta-blocker stimulation by using processed cine MRI of various resolution. We discuss our results and give some perspective remarks in Section 5 and conclude the paper in Section 6.

Data acquistion

Cardiac MRI was performed on 1.5T Philips Achieva system. The acquisition of a healthy volunteer dataset was performed at baseline and under the infusion of beta-blocker esmolol (50-200 µg/kg/min to achieve 10-20% decrease in heart rate) and the following data were acquired in each part of the study:

-Cine bSSFP sequence in retrospective ECG gating with spatio-temporal resolution 2x2x8mm and 40 time frames / cardiac cycle, FOV 350x350 mm. The short axis cine stack covering contiguously whole ventricles and standard long axis views of left ventricle (2-, 3-, and 4-chamber view) were taken. -3D tagged MRI of LV in prospective ECG triggering with acquired spatial resolution 3.4x7.7x7.7 mm (for 3 orientations of tag planes) reconstructed into one resulting 3D tagged image interpolated to 1x1x1 mm voxel size and temporal resolution ∼33 ms for both baseline and esmolol scan.

Data processing

Image processing consists first of spatial registration of all the image sequences (rigid registration using the Image Registration Toolkit IRTK3 ). Then, we select 3 slices from the short axis stack such that the inter-slice spacing corresponds to that one in real stress exam. Non-rigid image registration based motion tracking [START_REF] Shi | A comprehensive cardiac motion estimation framework using both untagged and 3D tagged MR images based on non-rigid registration[END_REF] is consequently employed to track the LV endo-and epicardial surfaces in cine MRI data. By motion tracking of various types (and combinations) of cine images, we generate three types of observations (deforming LV endo-and epicardial surfaces) which are used in sequel for the contractility estimation. These will be denoted as: Additionally, we extract full tissue displacements from 3D tagged MRI by using the tracking algorithm [START_REF] Shi | A comprehensive cardiac motion estimation framework using both untagged and 3D tagged MR images based on non-rigid registration[END_REF] tuned to the specific spatio-temporal resolution of the tagged data. The full displacements will not be used in this study as observations, however the error between the simulated and extracted displacement will be assessed. More details about image processing of such particular type of MRI dataset can be found in [START_REF] Chabiniok | Manifold learning for cardiac modeling and estimation framework[END_REF].

Modeling and parameter estimation framework

The model and data-asimilation procedure used in this work are described in [START_REF] Chapelle | An energy-preserving muscle tissue model: formulation and compatible discretizations[END_REF], [START_REF] Chabiniok | Estimation of tissue contractility from cardiac cine-MRI using a biomechanical heart model[END_REF], respectively, with more particular details of the estimator setup in [START_REF] Chabiniok | Manifold learning for cardiac modeling and estimation framework[END_REF]. In a nutshell, we use a continuum mechanics based model of biventricular anatomy. The model consists of an active contraction component (Bestel-Clément-Sorine model) and the passive tissue is represented by the visco-elastic material with the hyperelastic potential given by W e = κ 1 (J 1 -3) + κ 2 (J 2 -3) + κ(J -1ln J) (Ciarlet-Geymonat). The model is activated by an analytically-prescribed physiological electrical activation wave.

The data-assimilation is based on reduced-order unscented Kalman filtering [START_REF] Moireau | Reduced-order unscented Kalman filtering with application to parameter identification in large-dimensional systems[END_REF]. The discrepancy between the model and processed image data is evaluated by means of signed distances between the corresponding LV endo-and epicardial surfaces [START_REF] Moireau | Filtering for distributed mechanical systems using position measurements: Perspectives in medical imaging[END_REF]. The parts of the cardiac surfaces on which the model-data discrepancy is evaluated is given by the coverage of the heart by the observations (see red surfaces in Fig. 2). Anatomically-created AHA regions are encompassed by the full s.a. and 3-s.a.+3-l.a. observation surfaces. In the 3-slice s.a. case the mid-cavity AHA segments are nearly fully encompassed by those surfaces, however the top half of the basal segments and the apical segments are not included in the imaging data (Fig. 2,right).

The model is pre-calibrated manually at the baseline so that the global volume indicators correspond to the clinical data. The Windkessel models representing the circulations are calibrated so that the simulated end-systolic and end-diastolic aortic pressures correspond to the measured values at baseline, and physiological values for the pulmonary circulation are assumed. This precalibration is used as an initialization for the contractility estimation both for the baseline and esmolol cases. The estimation of the contractility parameter is performed on the standard AHA subdivision for all 3 types of observations. Fig. 3 (top) shows consistent estimates of myocardium contractility in AHA regions of a healthy volunteer at baseline using all types of observations. Under the beta-blocker stimulation, a decrease of contractility was detected in all segments when the full s.a. stack was used, and we can notice that the contractility decrease is lowest in the septal wall (particularly the apical septal part). When using the 3-s.a.+3-l.a. standard clinical protocol the drop of contractility is smaller but from a qualitative point of view it is similar to the full s.a. stack case -relatively homogeneous except for the septal wall. Finally, the 3-s.a. slices observations detect the lowest contractility drop and in particular "do not see" any contractility change in the regions not covered by the observation surfaces, i.e. the apical regions. A global decrease of contractility is expected under the beta-blocker infusion in a healthy volunteer. However, the results from Fig. 3 cannot be compared to any ground truth as this is not available in vivo. To support our result, we employ the displacements extracted from the 3D tagged MRI and assess the error of radial component of displacement in each node of the tetrahedral mesh within the AHA regions over the systole and average it over the AHA region i:

error i = 1 |n∈AHAi||t∈systole| t∈systole node n∈AHAi |(d t n,model -d t n,3Dtag
) • u r n |, with d denoting the incremental displacement between two consecutive time frames of tagged MRI and u r n the vector of radial direction. Fig. 4 shows the improvement of this radial error for the esmolol case after inputting the estimated contractility parameters into the heart model with respect to the initial model parameterization.

Discussion

The presented paper addresses a specific issue of applying model-based contractility estimation [START_REF] Chabiniok | Estimation of tissue contractility from cardiac cine-MRI using a biomechanical heart model[END_REF] in the quantification of dobutamine stress tests, used routinely to assess cardiac patients for consideration of revascularisation, namely the decreased spatial sampling of cine MRI. In this work we used esmolol, a negatively inotropic agent, to evaluate the optimal coverage and resolution of image data, and tested the whole data acquisition-processing-modeling and parameter estimation framework using a high quality volunteer dataset with pharmacologically modulated contractility.

As shown in Fig. 3 (top, estimation at baseline), if the model is sufficiently pre-calibrated by means of a low discrepancy with the observations, the estimator does not need to do much correction and the final regionally estimated parameters are consistent in between all types of observations. However, if the initialization is further apart from the data -as in our esmolol case -the estimator needs to correct more. As expected, the estimated quantities better correct the part of the model encompassed by the observation surfaces, and the regions further from those surfaces (e.g. apical region in the 3-s.a. observations) will be less sensitive and the estimator will perform worse. Evaluating Grammian matrix [START_REF] Moireau | Joint state and parameter estimation for distributed mechanical systems[END_REF] would provide an insight into the sensitivity of parameters in each region, and will enrich our work in future.

In this study we have focused on detecting the decrease of regional contractility and on the associated correction of the model displacements. In the absence of ground truth, we evaluated the error with respect to the displacements extracted from the 3D tagged MRI -the image data not used in the modeling-parameter estimation framework applied in this paper.

Visually, a simple qualitative assessment of the patterns of relative contractility decrease and the relative error correction for the esmolol case (Figs. 3 and4) supports the idea that the model is being corrected by decreasing the contractility. These figures show that in the septal segments (particularly the apical septal part) the estimator was not able to detect a contractility drop sufficiently and the error in the displacements in the model was not corrected. The apical septal regions may suffer either from a low quality of observations (suboptimal quality of segmentation) or from a more significant modeling error. Both may be caused by the complex geometry and fibre directions at the insertion place of RV. Adding the observations obtained by segmenting the RV surfaces might improve this issue.

In our tests, both the full cine stack and the 3-s.a.+3-l.a. protocol showed qualitatively similar decrease of contractility with esmolol. This suggests that the usual clinical protocol of 3-s.a.+3-l.a. could be as suitable for estimating the relative contractility change as the full s.a. cine stack. The 3-s.a.+3-l.a. protocol may however have two drawbacks: Firstly, the quantitative values of the detected contractility drop are significantly lower than if the full cine stack is used. Secondly, an important prerequisite for the motion tracking combining s.a. and l.a. views [START_REF] Shi | A comprehensive cardiac motion estimation framework using both untagged and 3D tagged MR images based on non-rigid registration[END_REF] is that those images are very well spatially aligned (to correct the non-reproducibility of patient's breath-holding). This spatial registration was efficient for our volunteer scan in which the full short axis cine stack provided sufficient information from the whole heart. When using patient data the spatial registration of the 3-s.a. slices with 3-l.a. slices would be a very challenging task.

The acquisition time for the full short axis cine stack in a spatio-temporal resolution as was used in our estimation trials is comparable to the acquisition of 3-s.a.+3-l.a. slices in standard temporal resolution. The presented work shows that a full-coverage short axis cine stack with a reduced temporal resolution may be an interesting alternative if quantification of results becomes a priority. We are aware that this hypothesis would need to be assessed by a study including a higher number of subjects, and it is one of our current targets.

We have calibrated the aortic pressure to the measured diastolic and systolic values at baseline, and we kept the pressure for the case with drug. It is reasoned by the fact that esmolol acts as a cardiac specific beta-blocker and therefore has a lower effect on the periphery circulation. Even if the pressure was not taken at each level of drug, the estimation would still provide a relative change in contractilities between the LV regions. Although we would be missing the scaling factor for the contractilities, the heterogeneity between the regions might provide the information to diagnose regional defects in contractility.

Long running times of the estimation (convergence of parameters) could be addressed by the preconditioning for instance by a reduced model [START_REF] Caruel | Dimensional reductions of a cardiac model for effective validation and calibration[END_REF] or statistical methods [START_REF] Chabiniok | Manifold learning for cardiac modeling and estimation framework[END_REF]. The state estimator [START_REF] Moireau | Joint state and parameter estimation for distributed mechanical systems[END_REF] was not used in this work to allow for a higher sensitivity to the parameters, but would most likely accelerate the convergence. The optimal setup of the estimator is our ongoing work. A coarser subdivision e.g. into 6 radial regions may be considered when only a limited spatial coverage by images is available (e.g. the 3-s.a. slices in our trial). This result of estimation might be of a similar clinical value for a significantly lower computational cost.

We did not consider tagged MRI as observations, although it was shown that they are superior to cine MRI for estimating the tissue contractility [START_REF] Imperiale | Constitutive parameter estimation using tagged-MRI data[END_REF]. In a very limited image acquisition time and the current practice of qualitative assessment of the dobutamine stress exam, it is unlikely that the tagged data would be acquired in near future in such an exam.

Conclusion

In this paper, we have investigated the types of data acquisition able to estimate a change of myocardial contractility in the clinical setup of dobutamine stress test. By carefully designing experimental data acquisition and adjusting the modeling-parameter estimation framework, the approach can provide an insight into the optimal imaging protocol to most accurately answer the clinical question. Simultaneously to the "computer-assisted diagnosis" we are therefore targeting the topic of "model-assisted optimal imaging".

The problem of estimating the contractilities in dobutamine stress tests has a number of unresolved subtasks and this work is providing initial steps to some of them, in particular the issue of coverage of myocardium by the image data, and the level of spatial discretization for the estimated parameters depending on the image data. Further development of the modeling and data-assimilation components, extending the number of subjects, challenging the image processing part by the image quality realistic in the dobutamine stress examinations and starting the estimation in the real dobutamine cases -with an accordingly modified scanning protocol -is a natural continuation of the presented work.

Fig. 1 .

 1 Fig. 1. Usual coverage of whole ventricles by short axis full cine stack (left) as oppose to the 3 short axis slices typically used in the stress exam (middle). The LV volume plot on the right demonstrates a decreased ejection during infusion of beta-blocker esmolol.

1 .Remark 1 .

 11 full s.a. stack: the short axis cine stack covering whole ventricles with the temporal resolution downsampled to 20 time frames / cardiac cycle; 2. 3-s.a.+3-l.a.: the selected 3 slices in short axis with 3 long axis views in original temporal sampling 40 time frames / cycle (standard clinical protocol); 3. 3-slice s.a.: the selected 3 slices of the full stack (40 time frames / cycle). The temporal resolution for full s.a. stack used in the parameter estimation is such that the MRI acquisition time corresponds to the standard clinical protocol 3-s.a.+3-l.a.

Fig. 2 .

 2 Fig. 2. From left to right: "pseudo 4-chamber" view reconstructed from full short axis cine stack consisting of 8mm-slices without gap with the contours of geometrical model; p4-ch view from 3-s.a. slices (inter-slice spacing 16mm); standard AHA segments. The red parts of LV endo-and epicaridal surfaces are used in the estimation to measure the data-model discrepancy.

Fig. 3 .

 3 Fig. 3. Contractility (in kPa) estimated using various type of observation surfaces at baseline (top) and during esmolol infusion (mid row); relative drop of contractility under the infusion of beta-blocker for each type of observations (bottom).

Fig. 4 .

 4 Fig. 4. Evaluation of radial error in the model (the initial parametrization and after the contractility estimation) with respect to radial displacements extracted from 3D tagged MRI: absolute error (top); relative drop of error per segment (bottom).

http://www.doc.ic.ac.uk/∼dr/software

Acknowledgments

The authors acknowledge the support of Engineering and Physical Sciences Research Council EP/H046410/1, British heart foundation grant NH/11/5/29058 and Cardiovascular Healthcare Technology Cooperative. In addition, the author are thankful to P. Moireau and D. Chapelle (Inria, France) for providing the HeartLab software library, used in this work for all modeling and estimation computations, and to D. Rueckert and W. Shi (Imperial College London, Ixico) for providing the IRTK based motion tracking and valuable discussions. This research was supported by the National Institute for Health Research (NIHR)

Biomedical Research Centre at Guy's and St Thomas' NHS Foundation Trust and King's College London. The views expressed are those of the author(s) and not necessarily those of the NHS, the NIHR or the Department of Health.