Qinghua Abdouramane Moussa Ali

Abdouramane Moussa
email: amoussaa@irisa.fr

Qinghua Zhang
email: zhang@irisa.fr

Adaptive observer based fault diagnosis applied to differential-algebraic systems ⋆

Keywords: fault diagnosis, nonlinear systems, differential-algebraic equations, adaptive observer, extended Kalman filter

published or not. The documents may come

Adaptive observer based fault diagnosis applied to differential-algebraic systems

INTRODUCTION

Many modern engineering systems can be modeled by an explicit Ordinary Differential Equation (ODE) of the form ẋ = f (x, u, t)

(1) where x and u represent respectively the (vectorial) state and input of the system, and the dot over x denotes the derivative in the time t. This equation has a long-term mathematical history, and a large number of analytical and numerical tools have been developed for its study.

However, in some cases such an explicit state space model for the dynamics of a given system is not available. The system may instead be described by an implicit differential equation of the form F ( ẋ, x, u, t) = 0

(2) This class of systems includes and is broader than state space systems in the ODE form (1). Some of the relations within this implicit model may not involve at all the time derivative ẋ, hence being purely algebraic equations. This motivates calling (2) a Differential-Algebraic Equation (DAE), a descriptor or singular system.

The theory on DAEs goes back a long time ago. The origins of this theory can be traced back to the work of K. [START_REF] Weierstrass | Zur Theorie der bilinearen und quadratischen Formen[END_REF]) and L. Kronecker [START_REF] Kronecker | Algebraische Reduction der Schaaren bilinearer Formen: Algebraische Reduction der Schaaren quadratischer Formen[END_REF]) on parameterized families of bilinear forms. The term algebraic-differential system was used in the circuit context by [START_REF] Brown | Derivative-explicit differential equations for rlc graphs[END_REF]). Some DAEs can be simply regarded as implicitly written ODEs, this is the case of DAEs in which the matrix of partial derivatives ∂F/∂ ẋ has full column rank. In principle, the theory developed in the framework of ODEs ⋆ This work was supported by the ASTech MODIPRO project. can be applied in this case. However, after repeated failure of numerical integration methods on certain DAEs, researchers initiated discussions on treatment of DAEs [START_REF] Petzold | Differential/algebraic equations are not ode's[END_REF], [START_REF] Lewis | A survey of linear singular systems[END_REF]).

In contrast to systems modeled by ODEs, where the theory for fault diagnosis is well-established [START_REF] Ding | Model-Based Fault Diagnosis Techniques -Design Schemes Algorithms and tools[END_REF], [START_REF] Isermann | Fault-Diagnosis System[END_REF] and references therein), for DAE systems such study is not as well developed. The few existing studies on fault diagnosis of DAE systems deal with linear case [START_REF] Benveniste | On the use of descriptor systems for failure detection and isolation[END_REF], [START_REF] Duan | Robust fault detection in descriptor linear systems via generalized unknown input observers[END_REF]), a certain class of nonlinearities [START_REF] Shields | Observer design and detection for nonlinear descriptor systems[END_REF], [START_REF] Zhang | Fault detection and isolation in nonlinear dynamic systems: A combined input-output and local approach[END_REF] or systems with all states measurable [START_REF] Polycarpou | Nonlinear fault diagnosis in differential algebraic systems[END_REF], [START_REF] Vemuri | Fault diagnosis of differential-algebraic systems[END_REF]). This paper considers the problem of nonlinear fault diagnosis for a large class of systems modeled by nonlinear DAEs. The main idea used in this paper for dealing with DAEs is based on a particular method for the discretization of DAEs [START_REF] Milne | Numerical Calculus[END_REF], [START_REF] Petzold | Differential/algebraic equations are not ode's[END_REF]). This discretization method results in explicit discrete time state space equations. The proposed fault diagnosis methodology is then based on an adaptive observer applied to the discrete time state space system. This adaptive observer based approach is an integration and extension of the observer based and parameter estimation approaches. The core of such fault diagnosis methods consists of an adaptive observer which is used both to estimate the monitored faults and to improve the robustness against model uncertainty due to parameter changes.

The paper is organized as follows. Different assumptions on the system and the fault are needed to solve the fault diagnosis problem. Section 2 is devoted to the outline of the approach discussed in this paper. We detail the assumptions made and the motivations of this contribution and explain how adaptive observer based approach is used for the generation of features, and applied to the detection and diagnosis of faults. Section 3 presents a simulation model of a heat exchanger by discarding flow approximations in favor of one-dimensional solutions of the Navier-Stokes in a device with flow ports in the ends such as a pipe. Then, in section 4 the proposed model is used to illustrate the proposed approach.

FAULT DIAGNOSIS OF NONLINEAR DAE SYSTEMS

Consider the DAE system 0 = F c ( ẋ, x, u, θ) (3a) y(t) = h(x, θ) (3b) where x ∈ R n is the state vector, u ∈ R m is the known variables vector (typically the control input vector), y ∈ R p is the output vector, F c and h are known smooth vector-valued functions (linear or nonlinear) respectively in R n and R p . The parameter vector θ represents some characteristics of the system subject to changes caused by faults possibly affecting the system. When the system is operating normally, θ is at its nominal value θ 0 .

With the above formulation, the diagnosis of fault can be achieved by estimating the parameter vector θ. For this purpose, the input u and the output y are the only available signals. The hidden variable vector x has to be dealt with in some manner during the estimation of the parameter vector θ.

The first difficulty of the problem formulated in this section is the implicit form of the state equation, or rather the DAE, in (3). This difficulty is dealt with by adopting a particular discretization of the DAE leading to a discrete time explicit state equation.

After discretizing the continuous time t to the discrete instants (t k = k∆) k=0,1,2••• , the continuous time state vector x(t) is discretized accordingly as x k = x(t k ), and the derivative ẋ(t k ) at time t k is approximated by a backward difference of x(t). The resulting discretized system equations write

F d (x k , x k-1 , u k , θ) F c x k -x k-1 ∆ , x k , u k , θ = 0 (4) y k = h(x k , θ) (5)
Assume that the matrix of partial derivatives ∂F d ∂x k is non singular, equation (4) implicitly defines a function [START_REF] Krantz | The Implicit Function Theorem: History, Theory, and Applications[END_REF])

x k = f d (x k-1 , u k , θ) (6)
This non singular ∂F d ∂x k assumption is equivalent to assume that 1 ∆

∂F c ∂ ẋ + ∂F c ∂x (7)
is non singular, with ∂Fc ∂ ẋ denoting the matrix of partial derivatives of F c with respect to its first argument. Alternatively, to directly transform the DAE (3) to an explicit continuous time state equation, a non singular matrix of partial derivatives ∂F c /∂ ẋ would be required. This condition is different from the one on (7), which concerns two matrices of partial derivatives of F c and is more likely to be satisfied in practice. For example, if some of the equations in (3a) are purely algebraic, or in other words, some components of F c , say grouped in Fc , do not include ẋ, then matrix ∂ Fc /∂ ẋ is clearly zero. As Fc must contain x (otherwise the corresponding equations would not contain any information about the state vector x), the matrix of partial derivatives ∂ Fc /∂x is not generally zero.

The implicitly defined function f d (x, u, θ) can rarely be found analytically and one usually implements an quasi-Newton iterative method [START_REF] Dennis | Quasi-newton methods, motivation and theory[END_REF]) to compute x k from x k-1 , u k and θ, by solving equation ( 4) for x k . If the implicit function value defined by ( 4) is not unique, the value the closest to the value of x k-1 is chosen. This choice should correspond to the true system trajectory if the discretization step size ∆ is small enough and if the system trajectory is continuous in time.

The partial derivatives of x k with respect to x k-1 and θ can be computed by the implicit function theorem

∂x k ∂x k-1 = ∂f d ∂x = ∂F c ∂ ẋ + ∆ ∂F c ∂x -1 ∂F c ∂ ẋ (8) ∂x k ∂θ = ∂f d ∂θ = - 1 ∆ ∂F c ∂ ẋ + ∂F c ∂x -1 ∂F c ∂θ (9)
where the matrices of partial derivatives are evaluated at the point

x k -x k-1 ∆ , x k , u k , θ .
From the relation ( 6), one commonly used method to solve the problem of diagnosis is to augment the state x k with the parameter vector θ and to implement an extended Kalman filter (EKF) [START_REF] Cox | On the estimation of state variables and parameters for noisy dynamic systems[END_REF]). While this approach has proved effective in some applications, at least in the case of state space systems, it has also some well known drawbacks. In particular, the fact of treating equally the state vector x k and the parameter vector θ as if they had similar dynamics may make the tuning of the EKF delicate.

The approach proposed in this section relies on an adaptive observer. Conceptually similar to the EKF applied to the augmented system, the adaptive observer has the advantage of being able to designed in two steps:

• first a state estimator by assuming that the parameter vector θ is known, • then the parameter estimator, which is coupled with the state estimator.

In the case of linear time varying systems, a comparison in some other aspects between the EKF applied to the augmented system and the adaptive observer can be found in the Appendix A of [START_REF] Li | An adaptive observer for joint estimation of states and parameters in both state and output equations[END_REF]).

Adaptive observer

After the previously presented discretization, the continuous time DAE system (3) becomes the discrete time state space system

x k = f d (x k-1 , u k , θ) + w k (10a) y k = h(x k , θ) + v k (10b
) where w k and v k represent errors caused by discretization and modeling/measurement errors ignored in the continuous time model. The purpose of this subsection is to propose an adaptive observer for joint estimation of the state x and the parameter θ in this discrete time state space system. Though adaptive observers with formally proved convergence have been developed through a constructive approach for some particular class of nonlinear systems [START_REF] Xu | Nonlinear system fault diagnosis based on adaptive estimation[END_REF], [START_REF] Zhang | An adaptive observer for sensor fault estimation in a class of uniformly observable nonlinear systems[END_REF], [START_REF] Farza | Adaptive observers for nonlinearly parameterized class of nonlinear systems[END_REF]), the design of such algorithms for wider classes of nonlinear systems remains a difficult task. Here we propose to follow an approach in the spirit of the EKF, in the sense that we apply an adaptive observer originally designed for linear systems to the nonlinear system (10) by linearizing it around the last values of the state and parameter estimates. While the convergence of this adaptive observer applied to linear systems has been proved formally [START_REF] Guyader | Adaptive observer for discrete time linear time varying systems[END_REF], [START_REF] Li | An adaptive observer for joint estimation of states and parameters in both state and output equations[END_REF]), the study on its convergence for general nonlinear systems, like the convergence of the EKF in general, remains an open problem, except in the particular cases studied in [START_REF] Xu | Nonlinear system fault diagnosis based on adaptive estimation[END_REF], [START_REF] Zhang | An adaptive observer for sensor fault estimation in a class of uniformly observable nonlinear systems[END_REF], [START_REF] Farza | Adaptive observers for nonlinearly parameterized class of nonlinear systems[END_REF]).

The main idea for designing this adaptive observer, as originally proposed in [START_REF] Zhang | Adaptive observer for multiple-inputmultiple-output (mimo) linear time-varying systems[END_REF] and for its nonlinear generalizations, consists of a two-step procedure: first design a state estimator by assuming that the parameter vector θ is known, then modify and associate the state estimator with an parameter estimation algorithm. In this design procedure the difference between the dynamics of the state variables (varying and governed by state equations) and of the parameters (typically constant) is explicitly taken into account.

Following this procedure, it is first assumed that the parameter vector θ is known in (10). A classical state estimator for general nonlinear systems of this form is the EKF, that can be formulated as follows:

xk = f d (x k-1 + K k-1 (y k-1 -ŷk-1 ), u k , θ) ŷk = h(x k , θ) ( 11
)
where xk is the state estimate, ŷk is the output estimate, the computation of the Kalman gain K k will be detailed later.

Remark that here the term EKF is in the sense that the Kalman filter is applied to the nonlinear system (10) by linearizing the system, at each time instant t k , around the last state estimate xk-1 . This is a widely used practice, though in general no convergence proof is formally established. The EKF in this sense is to be distinguished from the application of the (extended) Kalman filter to the augmented system by viewing the parameter vector θ as extra states. For the moment θ is simply assumed known.

Now let us consider the estimation of the parameter vector θ. Again inspired by the adaptive observer for linear systems [START_REF] Guyader | Adaptive observer for discrete time linear time varying systems[END_REF], [START_REF] Li | An adaptive observer for joint estimation of states and parameters in both state and output equations[END_REF]), the parameter estimate θk is recursively updated from the output error of the current model:

θk = θk-1 + Θ k (y k -ŷk )
where Θ k is a gain matrix to be specified. At the same time, the parameter vector θ in the state estimation equation (11) should be replaced by θk-1 . After this substitution, the state estimate becomes somehow biased because of the difference between θk-1 and the true θ. It is shown in [START_REF] Guyader | Adaptive observer for discrete time linear time varying systems[END_REF], [START_REF] Li | An adaptive observer for joint estimation of states and parameters in both state and output equations[END_REF]) that, in the case of linear systems, this bias should be compensated by an additional term in the state estimation equation in order to establish the convergence proof of the adaptive observer. Similarly, an additional term will be added to the modified state estimation equation ( 11). The resulting state and parameter estimation equations constitute the adaptive observer for the nonlinear system (10):

θk = θk-1 + Θ k (y k -ŷk ) (12a) xk = f d (x k-1 + K k-1 (y k-1 -h(x k-1 , θk-1 )), u k , θk-1 ) + Υ k ( θk -θk-1 ) (12b) ŷk = h(x k , θk-1 ) (12c)
with the gain matrices Θ k , K k and Υ k computed from the following equations:

P k = F k [(I -K k-1 H k-1 )P k-1 ]F T k + Q k-1 (13)
K k = P k H T k (H k P k H T k + R k ) -1 (14) Υ k = F k (I -K k-1 H k-1 )Υ k-1 + N k-1 (15) Ω k = H k Υ k + L k (16) S k = 1 λ S k-1 - 1 λ S k-1 Ω T k-1 Γ k-1 Ω k-1 S k-1 (17)
Γ k = λR k + Ω k S k Ω T k -1 (18)
Θ k = S k Ω T k Γ k (19) F k ∂f d ∂x (x k-1 ) (20) G k ∂f d ∂θ (x k-1 ) (21)
H k ∂h ∂x xk , θk (22)
L k ∂h ∂θ xk , θk (23)
where λ > 0 is a forgetting factor, Q k ∈ R n×n and R k ∈ R p×p are symmetric positive definite matrices corresponding to the covariance matrices of w k , v k when they are modeled as random noises and in ( 20) and ( 21), xk-1

(x k-1 + K k-1 (y k-1 -h(x k-1 , θk-1 )), u k , θk-1 ).
Notice that, compared to (11), the term Υ k ( θkθk-1 ) has been added to the state estimation equation (12b) to compensate the difference between θ and θk-1 .

Let us remark again that the formulation of this nonlinear adaptive observer has been in the same spirit as the EKF, by applying the adaptive observer initially designed for linear systems [START_REF] Guyader | Adaptive observer for discrete time linear time varying systems[END_REF], [START_REF] Li | An adaptive observer for joint estimation of states and parameters in both state and output equations[END_REF]) to the nonlinear system(10) linearized around the last state and parameter estimates at each time instant t k .

For the convenience of the reader, we recall the outline of the convergence proof in the linear case. Necessary conditions for the existence of the observer are the asymptotical convergence to zero of the state and the parameter estimation errors, respectivily defined by xk = x k -xk and θk = θ -θk in the noise-free case (w k = 0 and v k = 0).

The dynamics of the error xk in the linear case is as follows

xk = F k (I -K k-1 H k-1 ) xk-1 + (G k -F k K k-1 L k ) θk-1 -Υ k ( θk -θk-1 ) (24)
Define the variable z k = xk -Υ k θk . From the above relations, the dynamics of z k is simplified to an exponentially stable system

z k = F k (I -K k-1 H k-1
) zk-1 (25) and from ( 17)-( 19), the homogeneous part of the parameter estimation error governed by

θk = (I -Θ k Ω k ) θk-1 -Θ k H k z k
(26) tends exponentially to zero under some persistent excitation condition [START_REF] Guyader | Adaptive observer for discrete time linear time varying systems[END_REF], [START_REF] Narendra | Stable Adaptive Systems[END_REF]). Thus, by construction, we have

   lim k→+∞ z k = 0 lim k→+∞ θk = 0 ⇒ lim k→+∞ xk = 0
It is also possible to apply the EKF to the augmented system (treating θ as extra states) for joint estimation of x and θ. The proposed adaptive observer has two practical advantages. The gains K k for state estimation and Θ k for parameter estimation can be tuned in two steps in simulation studies. In the first step the parameter θ is assumed known, hence the tuning of K k is like in the case of the classical EKF for state estimation and afterward, in the second step, Θ k is tuned while the tuning of K k is fixed. Another advantage of the adaptive observer is that the recursive computations of K k and Θ k are separated, implying a lower numerical cost, compared to the fully coupled gain matrix computation in the EKF applied to the augmented system.

HEAT EXCHANGER MODEL EQUATIONS

ṁ2 in T 2 in ṁ1 in T 1 in ṁ1 out T 1 out ṁ2 out T 2 out Figure 1. A heat exchanger
To illustrate the fault diagnosis method presented in this paper, the monitoring of a heat exchanger is considered in this section. Heat exchangers are typical devices for heating, cooling, refrigeration and air conditioning. Their task is to transfer energy in the form of heat from one medium (e.g, a gas or liquid) to another. The heat transfer is from a hot stream with entering and exiting temperatures T 1 in and T 1 out and incoming and outgoing mass flow rates ṁ1 in and ṁ1 out , to a cold stream with entering and exiting temperatures T 2 in and T 2 out and incoming and outgoing mass flow rates ṁ2 in and ṁ2 out . A schematic of a heat exchanger with these temperatures and mass flow rates is shown on figure 1. The mathematical model of heat exchanger used in this study has been developed in the CSDL (Complex Systems Design Lab) project1

The actual heat transfer rate is given by the equalities

Q = C h (T 1 in -T 1 out ) (27) = -C c (T 2 in -T 2 out ) (28)
where C h and C c denote respectively the heat capacity rates (product of specific heat capacity and incoming mass flow rate) of the hot and cold stream.

If the heat exchanger effectiveness, specific heat capacities and incoming temperatures are known, the exit temperatures can be obtained as following : 29)

T 1 out = T 1 in -ǫ min{C h , C c } C h (T 1 in -T 2 in ) (
T 2 out = T 2 in + ǫ min{C h , C c } C c (T 1 in -T 2 in ) (30)
where ǫ is the heat exchanger effectiveness coefficient.

A heat exchanger can be described mathematically in terms of Navier-Stokes and heat transfer equations that are partial differential equations. We consider the hot fluid as working fluid with a quasi-one dimensional flow along a cylindrical tube of length L and lateral area A.

Following the finite-volume method, the computational domain is approximated by two cells of volumes v 1 and v 2 , each one characterized at fixed t by a temperature T i and pressure p i (i = 1, 2). Using the well-known ideal gas law ρ = p RT , we obtain the thermodynamic heat exchanger model given by the following DAEs :

                                 0 = v 1 R ṗ1 T 1 -p 1 Ṫ1 T 2 1 + v 2 R ṗ2 T 2 -p 2 Ṫ2 T 2 2 -ṁ1 + ṁ2 0 = v 1 ṗ1 h(T 1 ) RT 1 -1 + p 1 Ṫ1 R ∂h(T1) ∂T T 1 -h(T 1 ) T 2 1 + v 2 ṗ2 h(T 2 ) RT 2 -1 + p 2 Ṫ2 R ∂h(T2) ∂T T 2 -h(T 2 ) T 2 2 -ṁ1 h(T 1 in ) + ṁ2 h(T 2 ) + Q 0 = ṁ2 1 T 1 p 1 -ṁ2 2 T 2 p 2 - AL 2R ( m1 + m2 ) (31)
where h(T ) is the specific enthalpy of a fluid with temperature T , R is the specific gas constant, ṁ1 and ṁ2 are respectivily incoming and outgoing mass flow rate of working fluid.

By assuming known the following variable

• the mass flow rates ṁ1 , ṁ2

• incoming temperatures T 1 in , T 2 in • outgoing temperature T 1 out = T 2
and supposing the pressure p 2 (in the cell 2) measured, the model used for the diagnosis is given by system (32) below 0 = F c ( ẋ, x, u, φ) y = 0 1 0 x + ψ (32)

where the state x ∈ R 3 and the output y ∈ R are

x = p 1 p 2 T 1 , y = p 2 ,
the input vector u contains all known variables i.e. ṁ1,2 , m1,2 , T 1,2 in , T 2 , Q, the scalar values φ and ψ represent respectively the eventual efficiency loss and the pressor sensor bias. The parameters φ and ψ are written in the vectorial form as

θ = φ ψ .
The function F c is defined by

F c ( ẋ, x, u, φ) = f 1 ( ẋ, x, u, φ) f 2 ( ẋ, x, u, φ) f 3 ( ẋ, x, u, φ) with f 1 (z, x, u, φ) = v 1 R z 1 x 3 -x 1 z 3 x 2 3 + v 2 R z 2 T 2 -x 2 Ṫ2 T 2 -ṁ1 + ṁ2 f 2 (z, x, u, φ) = v 1 z 1 h(x 3 ) Rx 3 -1 + x 1 z 3 R ∂h(x3) ∂T x 3 -h(x 3 ) x 2 3 +v 2 z 2 h(T 2 ) RT 2 -1 + x 2 Ṫ2 R ∂h(T2) ∂T T 2 -h(T 2 ) T 2 2 -ṁ1 h(T 1 in ) + ṁ2 h(T 2 ) + φ Q f 3 (z, x, u, φ) = ṁ2 1 x 3 x 1 -ṁ2 2 T 2 x 2 - AL 2R ( m1 + m2 )
When the process operates under normal operating conditions, the nominal parameter value θ = θ 0 1 0 .

In the presence of a a loss of efficiency (i.e 0 < φ < 1) and/or a sensor bias (i.e ψ = 0), the parameters vector θ will deviate from the nominal value θ 0 .

SIMULATION EXAMPLE

In this section, we use the adaptive observer to monitor faults in a heat exchanger modeled by (32). The simulated faults are

• a degradation of the efficiency coefficient,

• a sensor bias.

The sampling period for the discretize time model is 1s and the simulation is performed during 1000s. Centered gaussian white noise is added to the input and output to satisfy the condition of persistent excitation required by the adaptive observer.

The recursively estimated values of φ and ψ are respectively illustrated in figures 2 and 3, where the dotted lines represent the true simulated parameter values, and the solid lines represent the estimated values.

The estimated states are compared with the simulated states in figures 4, 5 and 6, where dotted lines represent the true simulated state variables, and the solid lines represent the estimated values.

CONCLUSION

This paper has dealt with fault diagnosis in DAE systems.

We have focused our study on the diagnosis of faults modeled as parameter changes in a class of dynamic systems modeled by implicit nonlinear DAEs. The adaptive observer technique is used to accomplish the fault diagnosis task. We have shown that the method initially developed for linear systems based on adaptive observers can be extended to general nonlinear systems. The decision for fault diagnosis is based on the time evolution of parameter estimates. Simulation results are produced to illustrate the ability of the proposed approach to detect faults in a heat exchanger between two streams of dry air treated as an ideal gas. Like EKF-based methods, further studies should be made on the robustness of the proposed method to the severity of nonlinearities.

Figure 2 .Figure 3 .Figure 4 .Figure 5 .Figure 6 .

 23456 Figure 2. Graphical representations of the simulated parameter φ in dotted line and its estimate φ in solid line, over time (unit : second).

The CSDL (Complex Systems Design Lab) project funded byFUI (2009FUI ( -2012) ) is for the purpose of developing a collaborative platform for the design of complex systems. http://www.systematic-paris-region.org/fr/projets/csdl

ACKNOWLEDGEMENTS

The authors would like to thank Dassault Aviation and Dassault Systèmes for having provided the CSDL model in Dymola and for their assistance.


Abdouramane Moussa Ali

Qinghua Zhang
email: qinghua.zhang@@inria.fr

Abdouramane Moussa
email: abdouramane.moussaali@@inria.fr

Fault diagnosis for linear time-varying descriptor systems

Keywords: fault diagnosis, discrete-time descriptor dynamics, time-varying models, adaptive observer

In this paper fault diagnosis is studied for linear time varying descriptor systems, which are the discrete time counterpart of the continuous time dynamic systems described by differentialalgebraic equations. This class of systems includes and is broader than the well-known state space systems. The framework of descriptor systems is particularly useful for studying dynamic systems exhibiting time varying singularities. Actuator faults and sensor faults are respectively modeled as parametric changes in the state equation and in the output equation. The main result of this paper consists in extending an adaptive observer, initially designed for state space systems, to descriptor systems. Based on this result, fault diagnosis is performed by estimating the parameters characterizing actuator and sensor faults. Simulations examples are presented to illustrate the proposed method.

INTRODUCTION

Many modern engineering systems can be modeled by an explicit ordinary differential equation (ODE) of the form ẋ(t) = f (x(t), u(t))

(1) where x(t) and u(t) represent respectively the (vectorial) state and input of the system, ẋ(t) denotes the time derivative of x(t), and f (•, •) is some function characterizing the system. Such state space equations have a long-term mathematical history, and a large number of analytical and numerical tools have been developed for their study.

However, in some cases such an explicit state space model for the dynamics of a given system is not available. The system may instead be described by implicit differential equations, known as differential-algebraic equations (DAE), of the form F ( ẋ(t), x(t), u(t)) = 0 (2) where F (•, •, •) is some vector-valued function. If ẋ(t) can be solved for from (2), then the DAE can be converted to an ODE, but this operation is not always possible. It is thus necessary to study DAE systems in some situations.

After linearization around an operating point and discretization in time, the original nonlinear DAE system is approximately described by an implicit discrete time state space equation

E k+1 x(k + 1) = A k x(k) + B k u(k) + µ(k)
(3) where x(k), u(k) and µ(k) are respectively the discrete time state, the input and the modeling errors indexed by k = 0, 1, 2, ..., and the matrix E k may not be of full column rank.

⋆ This work was supported by the ITEA2 MODRIO project.

Systems governed by (3) are known as descriptor systems. This is a general and convenient framework for studying DAE systems which can be appropriately linearized. Some descriptor systems can be simply regarded as implicitly written state-space equations, this is the case of systems in which the matrix E k has full column rank. In principle, the theory developed in the framework of linear state-space equations can be applied in this case. However, even if the matrix E k is invertible for all k, it is preferable to use the descriptor equation because of the possible ill conditioning of E k .

Descriptor systems have attracted considerable attentions in recent decades where great efforts where made to investigate descriptor system theory and applications [START_REF] Nikoukhah | Kalman filtering and riccati equations for descriptor systems[END_REF]), Benveniste et al. (1993), [START_REF] Darouach | Design of observers for descriptor systems[END_REF], [START_REF] Shields | Observer design and detection for nonlinear descriptor systems[END_REF], [START_REF] Polycarpou | Nonlinear fault diagnosis in differential algebraic systems[END_REF], [START_REF] Zhang | Fault detection and isolation in nonlinear dynamic systems: A combined input-output and local approach[END_REF], [START_REF] Vemuri | Fault diagnosis of differential-algebraic systems[END_REF], [START_REF] Koenig | Design of proportional-integral observer for unknown input descriptor systems[END_REF]).

Fault diagnosis (detection and identification of fault) is rarely tackled in the descriptor case, in contrast to the case of systems with classical state space representations, where the theory is well-established [START_REF] Ding | Model-Based Fault Diagnosis Techniques -Design Schemes Algorithms and tools[END_REF], [START_REF] Isermann | Fault-Diagnosis System[END_REF] and references therein). For descriptor system fault diagnosis, most studies are observer-based approaches in the time invariant case [START_REF] Duan | Robust fault detection in descriptor linear systems via generalized unknown input observers[END_REF], [START_REF] Marx | Robust fault-tolerant control for descriptor systems[END_REF]).

In this paper, fault diagnosis will be studied for the general class of time-varying discrete time linear descriptor systems in the form of

E k+1 x(k + 1) = A k x(k) + B k u(k) + Φ k θ + µ(k) (4a) y(k + 1) = C k+1 x(k + 1) + Ψ k θ + ν(k) (4b) k = 0, 1, 2, • • • where x(k) ∈ R n is the state vector, u(k) ∈ R m
is the input vector (possibly including other known variables), y(k) ∈ R p is the output vector, E k+1 , A k , B k , C k are known real matrices of compatible sizes, µ(k) and ν(k) represent modeling/measurement errors and are modeled as independent, zero-mean, Gaussian vector sequences with covariance matrices Q k and R k respectively, the terms Φ k θ and Ψ k θ represent respectively faults possibly affecting the state and the output equation. The matrices Φ k ∈ R n×q and Ψ k ∈ R n×q are assumed known, whereas the parameter vector θ ∈ R q is unknown and subject to changes caused by faults. The fault-free system is characterized by the nominal value θ of the parameter vector.

The development of this paper is based, on the one hand, on the adaptive estimation techniques for explicit state space systems [START_REF] Zhang | Adaptive observer for multiple-inputmultiple-output (mimo) linear time-varying systems[END_REF]), and, on the other hand, on the Kalman filter for time-varying discrete linear descriptor systems. The core of the proposed fault diagnosis methods consists of an adaptive observer which is used both to estimate the monitored faults and to improve the robustness against model uncertainty due to parameter changes.

The paper is organized as follows. The first step in the monitoring of a system modeled by ( 4) is necessarily the state space x reconstruction when the only information available is contained in signals u and y. In Section 2, under some assumptions, a transformation of descriptor systems is introduced which allows us to use a standard Kalman filter for state vector estimation, in the fault-free case. Section 3 is devoted to the outline of the approach discussed in this paper. We detail the motivations of this contribution and explain how adaptive observer based fault diagnosis approach can be extended to general linear descriptor systems. In section 4, a simulation example is used to illustrate the effectiveness of the proposed approach.

KALMAN FILTER FOR FAULT-FREE DESCRIPTOR SYSTEMS

Here we recall some basic knowledge about the Kalman filter of discrete linear descriptor systems described by

E k+1 x(k + 1) = A k x(k) + B k u(k) + µ(k) (5a) y(k + 1) = C k+1 x(k + 1) + ν(k) (5b) k = 0, 1, 2, • • •
Like the classical Kalman filter for state space systems, the purpose of descriptor system Kalman filter is for state estimation from input and output signals. State reconstruction has been frequently studied in continuous time and usual for state space systems case (when E = I) since the early work in [START_REF] Kalman | A new approach to linear filtering and prediction problems[END_REF], [START_REF] Kreindler | On the concepts of controllability and observability of linear systems[END_REF], [START_REF] Silverman | Controllability and observability in time-variable linear systems[END_REF]) and in the time invariant descriptor case. Theoretical work on the observability of linear time-varying continuous time descriptor systems can be found in [START_REF] Campbell | Observability of linear time-varying descriptor systems[END_REF]).

In [START_REF] Nikoukhah | Kalman filtering and riccati equations for descriptor systems[END_REF]), the authors developed a Kalman filter for (5), based on maximum likelihood estimation. The filtered estimate x is defined as the maximum likelihood estimate of x based on (5) together with a priori knowledge of the initial condition x(0). Recursively, we obtain, from (5) and the estimate of x(k) denoted by x(k), the following equation

Y k = H k+1 x(k + 1) + V k , k = 0, 1, • • • (6) where Y k = A k x(k) + B k u(k) y(k + 1) , H k = E k C k , V k = w(k) ν(k) w(k) = -A k x(k) -µ(k) and x(k) = x(k) -x(k) is the state estimate error.
It has been shown in [START_REF] Nikoukhah | Kalman filtering and riccati equations for descriptor systems[END_REF]) (Lemma 2.3 and Lemma 2.4) that the (optimal) estimate x(k+1) of x(k + 1) based on ( 5) is the same as its (optimal) estimate based on (6).

Let R k be the covariance matrix of V k . It is shown in [START_REF] Nikoukhah | Kalman filtering and riccati equations for descriptor systems[END_REF]) that

R k = A k P k A T k + Q k 0 0
R k where P k is a symmetric positive definite matrix, which can be computed at each iteration of the maximum likelihood estimation algorithm, as recalled later in this section.

We shall make the following assumptions on system (5). A1 : The initial state x(0) is Gaussian with known mean x0 and variance P 0 and independent of µ and ν. This can be represented by the following equation x(0) = x0 + η (7) where η is zero-mean, Gaussian and independent of µ and ν, and with covariance P 0 . A2 : The matrix H k is supposed to be full column rank ∀ k ≥ 1 A3 : For all k ≥ 0, the variance matrices Q k and R k are positive definite.

From ( 6), the maximum likelihood estimate of x(k + 1) consists in maximizing the probability density function of V k parameterized by x(k + 1). This estimate is given by

x(k + 1) = arg max α e {-1 2 (Y k -H k+1 α)R -1 k (Y k -H k+1 α)} (8)
The solution of the maximization in (8) gives

x(k + 1) = (H T k+1 R -1 k H k+1 ) -1 H T k+1 R -1 k Y k (9)
The error covariance associated with the estimate x(k + 1) is

P k+1 = (H T k+1 R -1 k H k+1 ) -1
(10) Note that, the inverse matrices in ( 9) and ( 10) are well defined because of assumptions A2 and A3.

Consider the following partition

[L k K k ] = P k+1 H T k+1 R -1 k (11)
The estimate x(k + 1) can be explicitly expressed from ( 6), ( 9) and (11), yielding

x(k + 1) = L k A k x(k) + K k y(k + 1) + L k B k u(k) (12)
and the state estimation error is then governed by the equation

x(k + 1) = L k A k x(k) + L k µ(k) -K k ν(k) (13)
Let Pk be the predicted estimate covariance :

Pk = A k-1 P k A T k-1 + Q k-1 (14)
The following Riccati equation is obtained from (10)

Pk+1 = A k E T k+1 P -1 k E k+1 + C T k+1 R -1 k C k+1 -1 A k + Q k (15)
We know that the stability of the Kalman filter (the optimal estimator) is generally guaranteed if the covariance matrix Pk is bounded from above and from below by two strictly positive definite matrices for all k ≥ 0.

In [START_REF] Nikoukhah | Kalman filtering and riccati equations for descriptor systems[END_REF]), the asymptotic behavior of the error variance P k was examined in the time invariant case (with constant matrices A, C, E, Q, R). It can be shown in this case that this error variance Pk converges exponentially fast to the unique positive definite solution of the algebraic descriptor Riccati equation [START_REF] Nikoukhah | Generalized riccati equations for two-point boundaryvalue descriptor systems[END_REF]):

P = A E T P -1 E + C T R -1 C -1 A T + Q (16)
3. ADAPTIVE OBSERVER Consider the system in (4) subject to actuator and sensor faults. By defining L k and K k as in ( 11), the following explicit state recurrent equation is derived from (4):

x(k + 1) = L k A k x(k) + K k y(k + 1) + L k B k u(k) (17a) + [L k Φ k -K k Ψ k ]θ + L k µ(k) -K k ν(k) y(k + 1) = C k+1 x(k + 1) + Ψ k θ + ν(k) (17b) k = 0, 1, 2, • • •
We assume, in addition of assumptions A1 -A3, that

A4 : E k , A k , B k , C k , Φ k , Ψ k
are all known and bounded A5 : the parameter vector θ is constant In fault diagnosis, the terms Φ k θ and Ψ k θ can be used to model faults. In particular, they can be used to model actuator or sensor faults according to their location in the system model.

In this section, we aim to estimate the hidden variable vector x and the parameter vector θ based on the knowledge of the input u and the output y signals. One commonly used method to solve this problem is to augment the state x(k) with the parameter vector θ and to implement the Kalman filter (12). While this approach has proved effective in some applications, at least in the case of state space systems [START_REF] Cox | On the estimation of state variables and parameters for noisy dynamic systems[END_REF]), it has also some well known drawbacks. In particular,

• treating equally the state vector x(k) and the parameter vector θ as if they had similar dynamics may make the tuning of the Kalman filter delicate. • the computational cost of the Kalman filter increases mainly because of the larger covariance matrix of the augmented state estimation error.

The approach proposed in this section relies on an adaptive observer. We extend the existing results concerning the adaptive observer for fault diagnosis in linear state space systems [START_REF] Zhang | Adaptive observer for multiple-inputmultiple-output (mimo) linear time-varying systems[END_REF]) to the descriptor system.

Conceptually similar to the Kalman filter applied to the augmented system, the adaptive observer explicitly takes into account the difference between state variables and unknown parameters.

For the descriptor system described by (4) we propose the following adaptive observer

x(k + 1) = L k A k x(k) + K k y(k + 1) + L k B k u(k) + [L k Φ k -K k Ψ k ] θ(k) + ω k (18a) ŷ(k + 1) = C k+1 x(k + 1) + Ψ k θ(k) (18b) Υ k+1 = L k A k Υ k + [L k Φ k -K k Ψ k ] (18c) Ω k+1 = C k+1 Υ k+1 + Ψ k (18d) ω k = L k A k Υ k ( θ(k) -θ(k -1)) (18e) θ(k + 1) = θ(k) + Θ k+1 (y(k + 1) -ŷ(k + 1)) (18f) Θ k+1 = S k Ω T k+1 λR k + Ω k+1 S k Ω T k+1 -1 (18g) Γ k+1 = λR k + Ω k+1 S k Ω T k+1 -1 (18h) S k+1 = 1 λ (I -Θ k+1 Ω k+1 ) S k (18i) k = 0, 1, 2, • • •
where λ ∈ (0, 1) is a forgetting factor, θ(k) is the estimate of the parameter vector θ at time k, matrices gain L k and K k are computed as in (11). The term ω k added in (18a) and defined in (18e) may appear unusual. It will help to establish the proof of the convergence of this adaptive observer later in this section.

This adaptive observer is the discrete time counterpart of the continuous time algorithm presented in [START_REF] Li | An adaptive observer for joint estimation of states and parameters in both state and output equations[END_REF]. A simpler discrete time algorithm (no unknown parameter vector in the output equation) has been proposed in [START_REF] Guyader | Adaptive observer for discrete time linear time varying systems[END_REF]). Its generalizations to some particular class of nonlinear systems have been considered in [START_REF] Xu | Nonlinear system fault diagnosis based on adaptive estimation[END_REF], [START_REF] Zhang | An adaptive observer for sensor fault estimation in a class of uniformly observable nonlinear systems[END_REF], [START_REF] Farza | Adaptive observers for nonlinearly parameterized class of nonlinear systems[END_REF]).

The convergence of the adaptive observer (18) can be studied both in the noise-free case and in the noisecorrupted case. As the proofs are similar to those of [START_REF] Guyader | Adaptive observer for discrete time linear time varying systems[END_REF]) and [START_REF] Li | An adaptive observer for joint estimation of states and parameters in both state and output equations[END_REF]), let us formulate here its main lines.

Define the state and parameter estimation errors, respectively given by

x(k) = x(k) -x(k) (19) θ(k) = θ -θ(k) (20)
In the noise-free case (µ(k) = ν(k) = 0), the state estimation error is governed by the following equation derived from (17a) and ( 18a) 21), ( 18c) and (18e), the dynamics of z k is simplified to an exponentially stable system

x(k + 1) = L k A k x(k) + [L k Φ k -K k Ψ k ] θ(k) -ω k (21) Define the variable z(k) = x(k) -Υ k θ(k -1) (22) From (
z(k + 1) = L k A k z(k) (23)
The parameter estimation error is governed by θ

(k + 1) = (I -Θ k+1 Ω k+1 ) θ(k)-Θ k+1 C k+1 z(k + 1) (24)
The key issue is to prove the exponential stability of the homogeneous part of (24) given by

ζ(k + 1) = (I -Θ k+1 Ω k+1 ) ζ(k) ( 25
)
Lemma 3.1. Under conditions

• the matrix L k A k in (18c) is stable
• there exists an integer N and a real α > 0, such that for all k ≥ 0

αI ≤ k+N -1 i=k Ω T i+1 R -1 i Ω i+1 (26)
• the forgetting factor λ ∈ (0, 1)

• the initial gain matrix S 0 is symmetric positive definite then the matrix S k in (18i) is positive definite and bounded from above and from below by two strictly positive definite matrices for all k ≥ 0.

A proof of this lemma can be found in Appendix A. The stability of the matrix L k A k is directly related to the stability of the Kalman filter ( 12), which is assumed in this paper. In the literature the condition ( 26) is referred to as the persistent excitation condition [START_REF] Narendra | Stable Adaptive Systems[END_REF]).

The boundedness of S k ensured by the above lemma allows us to define the Lyapunov as follows

V k = ζ(k) T S -1 k ζ(k) (27) It follows V k+1 = ζ(k + 1) T S -1 k+1 ζ(k + 1) (28a) = λζ(k) T (I -Θ k+1 Ω k+1 ) T S -1 k ζ(k) (28b) = λV k -λζ(k) T Ω T k+1 Θ T k+1 S -1 k ζ(k) (28c) = λV k -λ [Ω k+1 ζ(k)] T Γ k+1 [Ω k+1 ζ(k)] (28d)
where Γ k+1 is the positive definite matrix given by (18h

) i.e [Ω k+1 ζ(k)] T Γ k+1 [Ω k+1 ζ(k)] ≥ 0.
We then conclude that V k+1 ≤ λV k . This means that for any λ ∈ (0 , 1), V k decreases exponentially and that ζ k and thus θk also tends exponentially to zero.

It then follows from ( 22) that, in the noise-free case,

   lim k→+∞ z k = 0 lim k→+∞ θk = 0 ⇒ lim k→+∞ xk = 0
In the noise-corrupted case, instead of ordinary convergence, the state and parameter estimation errors converge in the mean to zeros [START_REF] Guyader | Adaptive observer for discrete time linear time varying systems[END_REF]).

The proposed adaptive observer has two practical advantages. The gains L k and K k for state estimation and Θ k+1 for parameter estimation can be tuned in two steps in simulation studies. In the first step the parameter θ is assumed known, hence the tuning of L k and K k is like in the case of the classical Kalman filter for state estimation and afterward, in the second step, Θ k+1 is tuned while the tuning of L k and K k is fixed. Another advantage of the adaptive observer is that the recursive computations of the gains for state and parameter estimations are separated, implying a lower numerical cost, compared to the fully coupled gain matrix computation in the Kalman filter applied to the augmented system.

NUMERICAL EXAMPLE

In this example, we use the adaptive observer to monitor faults in a system modeled by ( 5) with the following data

E k+1 =   0 0 -1 k k + 1   , A k =    1 - k + 1 k + 0.5 0 1 k + 1    B k = 1 1 , u(k) =   
1 -e -0.018k 0.8 sin(0.03k)

1 + e -0.03k sin(0.1k)

   C k+1 = [0 1], θ = 1, Φ k =   1 -e -k 30 sin(0.1k) 1 -e -k 15 sin(0.5k)   , Ψ k = 0
It is easy to verify that the matrix E k+1 is not invertible and that the matrix E k+1 C k+1 is full column rank for all k ≥ 0. The sampling period for the discrete time model is 1s and the simulation is performed during 1000s. We simulate a degradation of the parameter θ starting from time instant k = 200. Centered Gaussian white noise is added to the output. The simulated output with noise y is plotted in figure 1. The simulated and estimated parameter (fault) are shown in figure 2, where the dotted lines represent the true simulated parameter values, and the solid lines represent the estimated values. The estimated states are compared with the simulated states in figures 3 and 4, where dotted line represents the true simulated state variables, and the solid line represents the estimated values. These results show that, after the transient time of about 100s corresponding to the transient time (unmodeled transient) of θ, the parameter estimates follow closely the evolution of the simulated parameter. It is then possible to detect the simulated faults and estimate their severity.

CONCLUSION

This paper has dealt with fault diagnosis in descriptor systems. We have focused our study on the diagnosis of faults modeled as parameter changes in a class of discretetime descriptor systems. Under realistic assumptions on the system, a Kalman filtering algorithm is derived. Then, the adaptive observer technique is used to accomplish the fault diagnosis task. We have shown that the adaptive observer initially developed for explicit linear state space systems can be extended to general linear descriptor systems. The decision for fault diagnosis is based on the time evolution of parameter estimates. Simulation results are produced to illustrate the ability of the proposed approach to detect faults.

Benveniste, A., Basseville, M., El Ghaoui, L., Nikoukhah, R., and Willsky, Alan, S. (1993). On the use of

M k+1 = λM k + Ω T k+1 R -1 k Ω k+1 M 0 = S -1 0 > 0 (A.1)
Let us first study the upper bound of M k . Due to the stability of L k A k , the matrix Ω k computed from bounded Φ k and Ψ k is also bounded. As 0 < λ < 1, the iterations in (A.1) is exponentially stable. Hence M k driven by the bounded term Ω T k+1 R -1 k Ω k+1 is also bounded. Now consider the lower bound of M k . A recursive evaluation of (A.1) leads to

M k+1 = λ k+1 M 0 + k j=0 λ k-j Ω T j+1 R -1 j Ω j+1 (A.2) M k+1 ≥ λ k+1 M 0 + ⌊ k N ⌋ i=1 iN -1 j=(i-1)N λ k-j Ω T j+1 R -1 j Ω j+1 ≥ λ k+1 M 0 + ⌊ k N ⌋ i=1 λ k-iN +1 iN -1 j=(i-1)N Ω T j+1 R -1 j Ω j+1
where ⌊•⌋ denotes the floor (greatest integer) function.

It follows from ( 26) that

M k+1 ≥ λ k+1 M 0 + α   ⌊ k N ⌋ i=1 λ k-iN +1   I (A.3)
Then, we can conclude that M k is bounded from below by a strictly positive definite matrix. Hence we have proved that M k has both an upper bound and a strictly positive definite lower bound.

From (A.1) the inverse of the matrix M k+1 writes

M -1 k+1 = 1 λ M k + 1 λ Ω T k+1 R -1 k Ω k+1 -1
(A.4) By using the matrix inversion lemma

(A + BCD) -1 = A -1 -A -1 B DA -1 B + C -1 -1 DA -1
it easy to check that S k = M -1

k satisfies S k+1 = 1 λ (I -Θ k+1 Ω k+1 ) S k (A.5)
with the initial condition S 0 .

Therefore, we conclude that S k is bounded from above and from below by two strictly positive definite matrices, for all k ≥ 0.

Figure 1 .Figure 2 .

 12 Figure 1. Graphical representations of the simulated noised output y over time (k).

Figure 3 .

 3 Figure 3. Graphical representations of the simulated statex 1 in dotted line and its estimate in solid line, over time (k).

Figure 4 .

 4 Figure 4. Graphical representations of the simulated statex 2 in dotted line and its estimate in solid line, over time (k).


M Tadie Fogaing

H N Yoshikawa

O Crumeyrolle

I Mutabazi

Heat transfer in the thermo-electro-hdyrodynamic convection under microgravity conditions

Keywords: 44, 24, +f Natural convection -47, 65, -d Magnetohydrodynamics and electrohydrodynamics

This article deals with the thermal convection in a dielectric liquid confined in an infinite-length capacitor with a temperature gradient under microgravity conditions. The dielectrophoretic force resulting from differential polarization of the liquid plays the role of buoyancy force with an effective electric gravity. It induces the convection when the Rayleigh number based on this electric gravity exceeds a critical value. Two-dimensional numerical simulation for a geometry with a large aspect ratio is used to determine the developed convective flow. The Nusselt number is computed for different values of Prandtl number and its dependence on the distance from the critical condition is determined. The viscous and thermal dissipations are found to be identical to those in the Rayleigh-Bénard convection which yields the same heat transfer.

Introduction

Rayleigh-Bénard convection in a fluid layer confined between two horizontal plates with a vertical temperature gradient along the gravitational acceleration g has become the prototype of modern nonlinear physics [START_REF] Manneville | Structures dissipatives chaos et turbulence[END_REF][START_REF] Bergé | Order within chaos: Towards a deterministic approach to turbulence[END_REF][START_REF] Cross | [END_REF]. Similar thermal convections can be induced by using external fields and specific properties of fluids even under microgravity. A magnetic field applied to ferrofluids [4] or an electric field applied to dielectric fluids [5] allow to generate a ponderomotive force which can be regarded as thermal buoyancy in an artificial gravity. In the present work, we are interested in the thermal convection induced by an electric field and a temperature gradient applied to dielectric fluids. The latter convection is often referred to as thermo-electro-hydrodynamic (TEHD) convection. It has been the subject of many investigations since few decades, particularly for the case of dielectric fluids confined between two horizontal parallel plane electrodes [6][7][8][9][10]. The driving force of the convection arises from the dielectrophoretic (DEP) force density F DEP = -1 2 E 2 ∇ǫ, where E is the magnitude of the electric field E and ǫ is the electric permittivity of the fluid. The DEP force is a consequence of differential polarization of the fluid. It is a component of the electrohydrodynamic force [START_REF] Landau | Landau and Lifshitz Course of Theoretical Physics[END_REF][START_REF] Melcher | Continuum Electromechanics[END_REF] and dominates over the other components when the applied electric field is alternating with a high frequency f and when the fluid is incompressible and has no mobile a Present address: Laboratoire J.-A. Dieudonné, UMR 7351 CNRS, Université de Nice-Sophia Antipolis b E-mail: innocent.mutabazi@univ-lehavre.fr boundaries [10]. Under a temperature gradient, the DEP force can be generated through the thermal variation of ǫ:

F DEP = - 1 2 E 2 ∇ǫ (T ) = -αρθg e + ∇ eθE 2 2 , (1)
where θ is the deviation from the reference temperature T 0 : θ = T -T 0 . The introduced electric gravity g e is proportional to the gradient of the electric energy stored in the fluid:

g e = e ρα ∇ ǫ 1 E 2 2 . (2)
In eq. ( 1), the thermal variation of the permittivity has been modeled by a linear relationship ǫ = ǫ 1 (1 -eθ) with the reference permittivity ǫ 1 at T = T 0 . The coefficient e is positive and of the order of 10 -3 -10 -2 K -1 for most dielectric fluids. The first term in eq. ( 1) is the analogue of the Archimedian buoyancy so that a thermal convection will be induced when the Rayleigh number based on the electric gravity g e exceeds a certain critical value. The second term in eq. ( 1) can be lumped with the hydrodynamic pressure term in the momentum equation.

The DEP force can induce convection even in a stable thermal stratification against the Earth's gravity g, i.e., the configuration where the temperature gradient is directed opposite to g, as shown in the linear stability theories [6][7][8]. The nonlinear behavior of the convection just above the critical state has been discussed by Stiles [8] and Yoshikawa et al. [10] by using a weakly nonlinear analysis and DNS, respectively.

Yoshikawa et al. [10] highlighted the role of the perturbation component of the electric gravity. In the unper-turbed state (fig. 1a), the thermoelectric buoyancy exerted on a fluid parcel is completely balanced by the pressure. A small perturbation δθ of the temperature of the parcel results in the buoyancy -αρδθḡ e V, where ḡe is the electric gravity in the unperturbed state and V is the volume of the parcel. It can be directed to the lower temperature zone to destabilize the base state. The temperature perturbation also results in a perturbation electric field δE through the thermal variation of ǫ, which is given by δE ∼ eδθ Ē, according to Gauss's law of electricity. The gradient of the electric energy is varied due to δE by -ǫ 1 ĒδE/D and ǫ 1 ĒδE/D at the upper and lower part of the parcel, where D is the diameter of the parcel: V = πD 3 /6. The thermoelectric buoyancy due to δE can then be estimated as -e θ + ∂ z θ D/2 -ǫ 1 ĒδE/D V/2 for the upper half of the parcel and as -e θ -∂ z θ D/2 ǫ 1 ĒδE/D V/2 for the lower half of the parcel. In the total, the force due to the electric field perturbation is

eǫ 1 ĒδE ∂ θ ∂z V 2 = e 2 ǫ 1 Ē2 δθ ∂ θ ∂z V 2 ( 3
)
The parcel is, in other words, subjected to a gravity g ′ e = -e2 ǫ 1 (∂ z θ) Ē2 /αρ. This perturbation gravity is directed opposite to the temperature gradient so that it has a stabilizing effect. This was confirmed in the linear stability theory and the analysis of the DNS results [10,[START_REF] Yoshikawa | [END_REF].

The present study is the continuation of the latter work and aims to develop, using DNS, the nonlinear properties of the TEHD convection in a two dimensional capacitor with an alternating electric tension and a temperature gradient. The paper is organized as follows: in the next section, we present the flow equations and the numerical solution procedure. The results are given and discussed in sect. 3. The last section is concerned with the conclusion.

" !" " " ! ! " ! ! " ! #!$ #$%&'(!)%&* + ,&-.-/01&/ 23&4
'5$(5$) 1 shows some values of g e,0 evaluated for different dielectric liquids.

e,0 = eǫ 1 V 2 0 B 3 /ραd 3 [log{(1- B/2)/(1+B/2)}] 2 . Table
In the electrohydrodynamic Boussinesq approximation, the equations for TEHD convection read:

∇ • u = 0, (4) ∂u ∂t + u • ∇u = -∇H + ∇ 2 u - L P r θg e , (5)
∂θ ∂t + u • ∇θ = 1 P r ∇ 2 θ, (6)
∇ • [(1 -Bθ) ∇φ] = 0. ( 7
)
where H is the dimensionless Bernoulli function. In the heat conduction equation ( 6), the viscous dissipation and Joule heating have been neglected, following the arguments developed by Yavorskaya et al. [14]. These equations must be solved together with boundary conditions at the electrodes (x = ±1/2) and at the adiabatic walls (z = ±Γ/2):

u = 0, θ = -1 2 , φ = 1 at x = -1 2 , (8)
u = 0, θ = 1 2 , φ = 0 at x = 1 2 , (9)
u = 0, ∂ z θ = 0, ∂ z φ = 0 at z = ± Γ 2 . ( 10
)
When the electric Rayleigh number L is smaller than a critical value, the homogeneous stationary state (u, θ, φ) = (0, θ(x), φ(x)) is established, where The electric gravity in this base state, ḡe , is given by ḡe = ḡe e x with ḡe = -1

θ = -x, φ = log 1+Bx 1+B/2 log 1-B/2 1+B/2 . ( 11
)
(1 + Bx) 3 . (12)
In the present work, we consider a two-dimensional system with a large aspect ratio Γ = 114 for different values of P r in order to simulate the nonlinear behavior of the TEHD convection in an infinite-length capacitor. The thermoelectric parameter will be set at a small value (B = 0.03). The set of partial differential equations ( 4)- (7) with the boundary conditions ( 8)- (10) are solved by the finite element method implemented in a commercial software package (COMSOL Multiphysics 3.5, Comsol AB, Stockholm, Sweden). Numerical grids are made of identical rectangles with sides of ∆x = 0.1 and ∆z = 0.15 so that the fluid domain is divided by 10 and 760 along the x and z directions, respectively. The backward differentiation formula is used for the time integration. The convergence of computation was verified by grid refinements. The initial fields are specified as null for the velocities, the temper-ature, and the electric field: The solved problem corresponds to a situation where the electric potential V 0 and the temperature difference θ are imposed instantaneously on a steady isothermal fluid layer at t = 0.

Results

Velocity, temperature and electric fields have been computed for a given set of parameter values (P r, B, L). All the obtained data exhibit a supercritical bifurcation from the conductive state to a convective regime with a critical electric Rayleigh number L c = 2130. As Yoshikawa et al. [10] have shown, the value of L c is independent from P r and coincides with its value determined by the linear stability theory. The critical modes are stationary and independent on P r.

Hot and cold fluids are, respectively, convected to the cold and hot plates by the thermal buoyancy due to the electric gravity and form convection rolls (fig. 3 a,c). Hot and cold regions are not symmetric with respect to the center line (x = 0) and yield a net enhancement of heat transfer by the convection, similar to the Rayleigh-Bénard convection. Perturbative electric fields are concentrated on the electrodes, where the store of the electrostatic energy can be intensified (fig. 3 e). This intensification is actually found but only by a tiny quantity (fig. 4). This means that the convection does not require much energy supply to be stored in the fluid.

In fig. 3, different quantities computed from the fields are also shown (b,d,f,g). To gain a better insight into the TEHD convection, we have used the equations of kinetic energy and temperature variance averaged over the flow domain:

d dt 1 2 u ′2 = w BG + w P G -Φ , ( 13
)
d dt θ ′2 2 = u ′ θ ′ - 1 P r |∇θ ′ | 2 , ( 14
)
where means the average over the whole fluid domain: = Γ -1 dxdz, and the primes indicate perturbation components. Different terms w BG , w P G and Φ in the energy equation ( 13) are the power input by the basic electric gravity, the power input by the perturbation electric gravity and the viscous dissipation function, respectively:

w BG = - L P r θ ′ u ′ •ḡ e , w P G = - L P r θ + θ ′ u ′ •g ′ e (15)
Φ = 2 ∂u ′ ∂x 2 + ∂u ′ ∂z + ∂w ′ ∂x 2 + 2 ∂w ′ ∂z 2 . ( 16
)
The power input by the basic electric gravity is concentrated at hot and cold cores (fig. 3 f). The power input by the perturbation electric gravity is positive near the electrodes, while it takes negative values over large zones inside the gap (fig. 3 g). The viscous dissipation is primarily due to the shear at the electrodes (fig. 3 b).

In the average over the whole fluid domain, the basic electric gravity provide energy to the flow, while the perturbation electric gravity impedes the convection. Figure 5 shows different terms in eq. ( 13) as a function of the normalized distance δ from the critical electric Rayleigh number: δ = L/L c -1. The viscous dissipation balances with the contributions from the different electric gravity components. The relative importance of the three contributions is independent of the Prandtl number except for small P r (< 1). Furthermore, it converges to the same value as δ → 0 for all the examined values of P r (fig. 6):

! !"!# !"$ !"$# !"% &!"!' &!"!% &!"!$ ! !"!$ !"!% !"!' !"!( )*+,-./0./01*2+34.*5631/70- !" !"# # !" $"# # $ 8.%&'&(9 &.)* &$ 8.:!'&: %. 9 +(!,#
! !"!# !"$ !"$# !"% !"%% !"%& !"%' !"%( !") !")% !")& *+,+! "# $ -+.+,+%+- & "' $ $! '# $!! $! ) !"!$ !"!' !"!( !"$ !"% !"#
w P G / Φ = -0.24 + O (δ) . ( 17
)
The ratio of w BG to Φ is given by w BG / Φ = -w P G / Φ + 1 = 1.24 + O(δ).

Discussions

The heat transfer enhancement is given by the Nusselt number N u:

N u = 1 + 1 Γ Γ/2 -Γ/2 (-∂ x θ ′ + P ru ′ θ ′ ) dz (18)
which compares the convective to conductive heat transfers.

Yoshikawa et al. [10] showed that the Nusselt number in saturated TEHD convection is proportional to the distance as N u -1 = Cδ with the proportionality constant C = 0.78 for P r larger than 1. For small P r, the constant C varies with P r as it does in the Rayleigh-Bénard convection (fig. 7). Seeking a similar correlation to that for the Rayleigh-Bénard convection [15], one can find

N u -1 = 1.28 - 0.0273 P r + 0.0077 P r 2 -1 δ. (19)
For small Prandtl number, the difference in N u between these two convections is small. In the saturated steady state, one can show that N u can be computed by

N u -1 = P r u ′ θ ′ . ( 20
)
Making use of eq. ( 14), one finds that the Nusselt number is given by the averaged value of the thermal dissipation:

N u -1 = |∇θ ′ | 2 .
The latter relationship is also found in the Rayleigh-Bénard convection. For both convections yielding the same heat transfer, the thermal dissipations is identical to each other.

In the limit of small B, the basic electric gravity is uniform: ḡe = -e x , so that w BG = P r -1 Lθ ′ u ′ . Indeed, w BG and u ′ θ ′ are identical in our simulation for B = 0.03 except the factor P r -1 L (fig. 3 d,f). The heat transfer can then also be computed by

N u -1 = P r L w BG = P r L (-w P G + Φ ) . (21)
Since the ratio of w P G to Φ is given by eq. ( 17) in the vicinity of the criticality, we have N u-1 = 1.24 P rL -1 c Φ = P r Φ /1720. One can also derive the similar relationship to eq. ( 21) for the Rayleigh-Bénard convection: N u -1 = P rR -1 Φ [16,17]. It reads in the vicinity of the criticality N u -1 = P r Φ /1708. For both convections yielding the same heat transfer, the viscous dissipations is hence also be identical to each other.

Conclusion

The thermo-electro-hydrodynamic convection in dielectric fluids represents a simple way of realizing thermal convection under microgravity conditions. We have highlighted the effect of the perturbative gravity against the destabilization effect of the basic electric gravity. The heat transfer by the TEHD convection has been quantified for different values of the Prandtl number. This effect becomes important for large values of the thermoelectric coupling parameter.

!!"!#$"#!%!&'($)

$$! %& %& !'$'('#")*'+''''''''''''','''''''''''''-%&! $ %& %& * !"!)./ !"!!..

Fig. 1 .

 1 Fig. 1. Electric field perturbation induced by a small temperature perturbation δθ of a parcel of dielectric fluid.

!

Fig. 2 .

 2 Fig. 2. Geometrical configuration of the problem

Fig. 3 .

 3 Fig. 3. Saturated flow fields and related quantities at the central part of electrodes: (a) velocity vectors u ′ with vorticity intensity |∇ × u ′ |, (b) viscous dissipation function Φ, (c) isotherms θ ′ = cst (solid and broken lines indicate positive and negative θ ′ , respectively) with thermal dissipation |∇θ| 2 , (d) convective heat transfer P r -1 L u ′ θ ′ , (e) equipotentials φ ′ = cst with electric field intensity |E ′ |, (f) power input by the basic electric gravity wBG, (g) power input by the perturbation electric gravity wP G and . The basic electric gravity is directed from right to left.

Fig. 4 .

 4 Fig. 4. Perturbation in electrostatic energy ∆UE = (ǫE 2ǭ Ē2 )/2 (P r = 100, B = 0.03).

Fig. 5 .

 5 Fig. 5. Different energy generation terms computed after the saturation (P r = 100, B = 0.03).

Fig. 6 .

 6 Fig. 6. Relative importance of WP G to Dv computed for saturated convection (B = 0.03).

  This work benefited from a financial support from CNES (French Space Agency), the CPER-Haute Normandie under the program THETE and from the FEDER. M.T.F. thanks the MESR for financial support during her Ph.D preparation. H.N.Y. acknowledges the financial support of the French National Research Agency (ANR), through the program Investissements d'Avenir (ANR-10 LABX-09-01), LABEX EMC 3 .

Fig. 7 .

 7 Fig.7. Proportionality constant C: N u -1 = Cδ, for different Prandtl numbers P r in the vicinity of the criticality (B = 0.03). The broken-line curve is the corresponding constant in the Rayleigh-Bénard convection: C = (0.6994 -0.00472/P r + 0.00832/P r 2 ) -1 (R/Rc -1)[15].

Table 1 .

 1 Electric gravity in some dielectric liquids in a capacitor with d = 10 mm under an electric voltage V0 = 10 kV and a temperature difference 1 K. The working temperature is 25 • C.

	Acetonitrile	0.777	1.38	36	0.155	14.7
	Nitrobenzene	1.198	0.830	34.9	0.188	32.2
	Acetone	0.785	1.43	19.1	0.086	0.756
	Chlorobenzene	1.101	0.985	5.61	0.0157	0.0003

Liquid

ρ (10 3 kg/m 3 ) α (10 -3 K -1 ) Dielectric constant e (K -1 ) ge,0 (m 2 /s)

Flow equationsWe consider a dielectric fluid confined in a capacitor with a temperature gradient and an alternating voltage V (t) = √ 2V 0 sin(2πf t) imposed over the gap between the electrodes. (fig.2). Characteristic scales can be used to introduce non-dimensional control parameters: the viscous relaxation time τ ν for the time, the gap d between the electrodes for the length and ∆θ for the temperature. The