Romeo Tatsambon Fomena

Christophe Collewet Cemagref
email: christophe.collewet@inria.fr

Fluid Flows Control using Visual Servoing

Keywords:

This paper introduces for the first time the control of a fluid flow using visual servoing. The control consists to regulate a flow around a desired velocities profile. The originality here is the possibility to use optical flow measurements, computed from the observed flow, in a visual servo control scheme. Compared to existing approaches that uses a limited set of measurements to control a fluid flow, our approach is shown to be a major improvement in terms of the flow state estimation and control.

INTRODUCTION

The goal of this paper is to show the benefit of using visionbased control in a closed-loop scheme to control a fluid flow, that is to modify the current state of a flow in order to reach a desired state. Closed-loop vision-based control, also called visual servoing, is now a well established technique in the robotics community (see [START_REF] Espiau | A new approach to visual servoing in robotics[END_REF]). More precisely, to achieve a visual servoing task, a set of visual features s(t) is selected from the image of the scene. A control law is then designed so that this set of visual features s(t) reaches a desired value s * corresponding to a desired state of the system. The dynamic of the error vector e(t) = s(t)s * is then given by

ė(t) = ∂e(t) ∂t + L e (t)u(t), ( 1
)
where u is the system control inputs, L e (t) is the jacobian matrix encoding the time variation of the visual features to the variation of the control signal acting on the system and ∂e(t)/∂t expresses the variation of the error vector due to the free motion of the visual features. The control principle is to regulate the error vector e to zero. The control law is built from (1) using the knowledge of the approximations L e (t) and ∂e(t) ∂t . Visual servoing has shown impressive results in numerous complex contexts such as underwater, medical and aerial (helicopters, blimp) robotic as shown in [START_REF] Bonin-Font | Visual navigation for mobile robots: a survey[END_REF].

In this work, we aim at stabilizing a flow, i.e. maintaining a flow in a desired state whatever the external perturbations. Fluid flows control is of great economical interest in aerospace. For instance reducing the drag of an airplane while enhancing lift can reduce fuel consumption. This can be done by preventing transition to a turbulent air flow around the plane wings.

Flows control can be achieved in two different ways: passive or active control. Passive control provides a permanent action on the system. Most often it consists in optimizing shapes or in choosing suitable surfacing (see [START_REF] Choi | Direct numerical simulation of turbulent flow over riblets[END_REF]). The main problem with this approach is that the control is, of course, inefficient when the system pa-rameters change. Conversely, in active control an external energy is required to act on the system, such as techniques based on blowing and suction as described in [START_REF] Joshi | A systems theory approach to the feedback stabilization on infinitesimal and finite amplitude disturbances in plane Poiseuille flow[END_REF]. However, most often, only open-loop control is used. As a matter of fact, designing a closed-loop control law requires the use of sensors that can be at the same time non-intrusive, accurate and adapted to the time and space scale of the phenomenon under monitoring. Unfortunately, such sensors are hardly available in the real context of control applications. The most commonly used measurement, obtained from MEMS, is the shear stress at a limited set of measurement points (see [START_REF] Joshi | A systems theory approach to the feedback stabilization on infinitesimal and finite amplitude disturbances in plane Poiseuille flow[END_REF]).

To deal with these issues, we propose a closed-loop visionbased control approach. As far as we know, the proposed approach has never been used for flows control issues. By using vision, the complete flow can be fully observed in real time (see [START_REF] Champagnat | How to calculate dense piv vector fields at video rate[END_REF]). This full observation is shown to be of great improvements for flows control in comparison with existing control approaches that use only a limited set of walls shear stress measurements.

In this paper, we focus on 2D plane Poiseuille flow which is well-known in fluid mechanics. The control framework for this flow is presented in Section 2. We then compare, in Section 3, the existing control approaches and the proposed vision-based control approach. In Section 4, the proposed theoretical results are validated in simulation using synthetic data sets of spatio-temporal variations of flow velocities.

2D PLANE POISEUILLE FLOW MODELING FOR CONTROL DESIGN

In this section we first present the basics of Poiseuille flow, then we introduce the control principle for this flow and finally we present the reduced linearized model used to derive the existing shear stress based control laws.

Basics of the Poiseuille flow

Poiseuille flow is a flow in an infinite length channel due to a pressure gradient.

       ∂V ∂t + (∇V)V = -∇P + 1 R e ∇ 2 V ∇ • V = 0 V(x, y = ±1, t) = 0 ( 2
)
where P is the pressure; V is the flow velocity; V(x, y = ±1, t) = 0 represents the no slip boundary condition ; R e is a dimensionless number called the Reynolds number.

Since the Poiseuille flow is simple, the analytical solution (V bx , V by , P b ) of (2) in the steady state case, i.e.

dV dt = ∂V ∂t + V • ∇V = 0, can be found: (V bx , V by , P b ) = (1 -y 2 , 0, - 2 R e x). (3)
This solution is illustrated in Fig. 1.

Control principle of the Poiseuille flow

A 2D Plane flow Poiseuille flow can be controlled via boundaries. Boundary control consists in modifying the upper and lower boundary conditions (see [START_REF] Joshi | A systems theory approach to the feedback stabilization on infinitesimal and finite amplitude disturbances in plane Poiseuille flow[END_REF]). The boundary control on the upper and the lower channels can be theoretically represented by functions χ u and χ l that allows mass conservation in the controlled system as pictured in Fig. 2. Note that in the absence of control the red dashed curves (see Fig. 2) are aligned with the lower and upper boundary lines as expected.

The existing active control approaches, presented in [START_REF] Joshi | A systems theory approach to the feedback stabilization on infinitesimal and finite amplitude disturbances in plane Poiseuille flow[END_REF], [START_REF] Bewley | Optimal and robust control and estimation of linear paths to transition[END_REF] and [START_REF] Mckernan | Control of plane Poiseuille flow: a theoretical and computational investigation[END_REF], focus on the stabilization of the flow, i.e. maintaining the flow in a desired state, typically the steady state, whatever the external perturbations. In this case the flow is initially in the steady state but in an unstable equilibrium, i.e. a small velocities pertubation value V p (x, y, t) destabilizes the non-controlled fluid flow.

Reduced linearized model of Poiseuille flow

The modeling, required to derive the control law, consists first of all in linearizing the Navier-Stokes equation (NSE) around the steady state solution. Then the continuous linearized model of the NSE of the flow is first represented in the Fourier domain at a specific wavenumber α n = 2π/L c containing the instability: V p (x, y, t) = V n p (y, t)e jαn . Then the obtained spectral representation V n p (y, t) is discretised in the y direction by pressure High Low pressure

M V(M, t) x y x = 0 x = L c y = 1 y = -1 χ l (x, y 1 , t) χ u (x, y M , t)
Fig. 2. Boundaries control.

using a number of M Gauss-Lobatto collocation points

y k = cos ((k -1)π/(M -1)) , 1 ≤ k ≤ M .
Gauss-Lobatto collocation points are chosen to improve the spectral accuracy. The obtained reduced linearized model is given by (see [START_REF] Mckernan | Control of plane Poiseuille flow: a theoretical and computational investigation[END_REF]):

   ṗn (t) = A n p n (t) + B n u(t) z(t) = C ⊤n p n (t) p n (0) = p n 0 ( 4
)
where p n (t) is the state vector, A n is the state matrix, u(t) is the system input, like blowing or suction actions on the channel boundaries, B n is the input matrix, C n is the output matrix and z(t) is the vector of shear stress measurements on the upper and lower boundaries.

In the classical case (see [START_REF] Joshi | A systems theory approach to the feedback stabilization on infinitesimal and finite amplitude disturbances in plane Poiseuille flow[END_REF], [START_REF] Mckernan | Control of plane Poiseuille flow: a theoretical and computational investigation[END_REF], [START_REF] Bewley | Optimal and robust control and estimation of linear paths to transition[END_REF]) of Reynolds number Re = 10 000 and wavenumber α n = 1, the reduced linearized model of Poiseuille flow ( 4) is unstable as shown in [START_REF] Orszag | Accurate solution of the Orr-Sommerfeld stability equation[END_REF]. This instability can be seen through the poles of the state matrix A n . These poles are illustrated on Fig. 3, where we can see the unstable mode λ = 0.00373967 ± i0.23752649. Control theory for linear systems can be used to cancel this instability.

In the next section, we compare the existing control approaches with our vision-based approach.

COMPARISON OF THE CONTROL APPROACHES

Existing control approaches

Using the classical output feedback control u(t) = -k ⊤ z(t), the controller generally fails to suppress unobservable high transient modes which could trigger transition to turbulence as shown in [START_REF] Joshi | A systems theory approach to the feedback stabilization on infinitesimal and finite amplitude disturbances in plane Poiseuille flow[END_REF]. The currently existing approaches focus on the state feedback LQR (see [START_REF] Mckernan | Control of plane Poiseuille flow: a theoretical and computational investigation[END_REF]) and the output feedback LQG regulation (see [START_REF] Bewley | Optimal and robust control and estimation of linear paths to transition[END_REF]). While the LQR approach can not be used in practice since it requires the knowledge of the true state value p n (t), the LQG approach is built on an observer. This observer provides an estimated value p n (t) of the state vector from the shear stress measurements z(t) using a Linear Quadratic Estimation (LQE) scheme. The control signal for the output feedback LQG regulator is thus given by

u(t) = -k ⊤ p n (t),
(5) where vector k is the LQR optimal gain.

In the next section we focus on state estimation using an observer in existing approaches and using optical flow measurements in our approach.

Imaginary axis

Poles Zeros

Unstable mode

Fig. 3. Poles and zeros of the reduced linearized system.

Flow State estimation

Before presenting the vision-based state estimation approach we briefly recall the LQE of the state vector using a limited number of shear stress measurements.

Existing LQE state estimation.

In the LQE framework it is supposed that the reduced linearized system (4) has process disturbances ε p and measurements noise ε z . In addition ε p and ε z are uncorrelated Gaussian white noise with covariance matrices Ξ p and Ξ z respectively. The reduced linearized system (4) is thus rewritten as

         ṗn (t) = A n p n (t) + B n u(t) + ε p (t) z(t) = C ⊤n p n (t) + ε z (t) E ε p ε ′ p = Ξ p , E {ε z ε ′ z } = Ξ z E ε z ε ′ p = 0, E {ε p ε ′ z } = 0 p n (0) = p n 0 , (6)
where x ′ denotes the conjugate transpose of vector x.

From the plant model given in (6), the state vector LQE is given by ˙

p n (t) = A n p n (t) + B n u(t) + L(z(t) -C n⊤ p n (t)) p n (0) = unknown, (7)
where p n (t) is the estimated state vector, and L the LQE optimal gain. From ( 6) and ( 7), the estimation error dynamic is given by

˙ p n (t) -ṗn (t) = (A n -LC n⊤ )( p n (t) -p n (t)) + Lεz(t) -εp(t) p n (0) -p n (0) = unknown. (8)
A well-known key issues in the LQE is the choice of the initial value p n (t = 0) that has to be fed to the observer (see 7)). Indeed this initial value highly influences the convergence time after which p n (t) equals the true value of the state vector p n (t). Modeling the initial condition p n (t = 0) with known physical statistics about the studied flow and using a time-varying estimator gain L(t) can reduce the time to which p n (t) equals the true value of the state vector p n (t) as proposed in [START_REF] Hoepffner | State estimation in wall-bounded flow systems. part 1. perturbed laminar flows[END_REF]. But this solution is still not satisfactory since it needs additional a-priori parameters in the model of the initial condition. In addition, it is also well-known that the above presented LQE approach produces a noisy estimation when shear measurements are noisy (see ( 8)).

The initialization and the noise issues in the LQE could provoke turbulence in the controlled flow and divergence of the control as shown in [START_REF] Mckernan | Control of plane Poiseuille flow: a theoretical and computational investigation[END_REF].

3.2.2. Vision-based state estimation in the general case of noisy measurements. In this section, we overcome the above mentioned issues related to state vector estimation Our contribution relies on the fact that from the knowledge of the optical flow it is possible to estimate the flow perturbation velocities field V p (x, y, t) and therefore to compute the estimation p n (t) of the state vector. The initial value p n (t = 0) is therefore no longer of concerned in our approach. This estimation is used thereafter in a closed-loop LQR scheme: the obtained closed-loop LQR is thus a visual servoing control scheme. Indeed using s(t) = p n (t), around the steady state, the error dynamic given by ( 1) is identified as the first equation in (4).

The optical flow can be obtained from the fronto-parallel projection on the 2D flow (see Fig. 4) using the technique described in [START_REF] Heitz | Variational fluid flow measurements from image sequences: synopsis and perspectives[END_REF]. The obtained image velocties field is exactly 1 d V p (x, y, t) which is proportional to the 2D flow perturbation velocity V p (x, y, t) up to the scale factor d > 0 (distance between the flow plane and the image plane). In the following we assume that the distance d is obtained by a calibration technique, in which case we can measure the flow velocity field V p (x, y, t).

Let ε(x, y, t) be a 2D independent identically distributed white Gaussian process. Let V p be the M × N pixels size image of the pertubation velocity obtained from noise-free optical flow measurements. In the practical case where optical flow measurements provide noisy velocity perturbation V p (x, y, t) = V p (x, y, t) + ε(x, y, t), these noisy measurements are used to compute estimation of the state vector p n (t) as shown in Appendix A. The estimated state vector is given by

p n (t) = p n (t) + 1 N e n (t), (9)
where vector e n (t) is related to the measurements noise matrix ε(x, y, t). It is clear from (9), that the larger the value of N the closer p n (t) is from p n (t).

Now we compare the behavior of the system closed with the two control approaches.

3.3

Comparison of the behavior of the closed-loop systems 3.3.1. Behavior of the system closed by the LQG control.

Combining the estimation error dynamic given in (8) and the LQR state dynamic (plugging (5) into the first equation of ( 6)), we obtain

   ṗn (t) ˙ p n (t) -ṗn (t) = M lqg p n (t) p n (t) -p n (t) + Mn εp(t) εz(t) p n (0) -p n (0) = unknown (10)
where

M lqg = A n -B n k ⊤ -B n k ⊤ 0 A n
-LC ⊤n and M n = I 0 -I L , with I an identity matrix. Equation ( 10) shows that the true state dynamic ṗn (t) depends on p n (t) and on p n (0). Therefore, a poor initialization of the observer could cause the divergence of the non-linear system describing the fluid flow. In addition the noise in the measurements propagates (see ( 10)) in the control law (see ( 5)); and this is not suitable for actuators in practice.

Behavior of the system closed by the vision-based control.

In the case of measurements noise in the optical flow, the reduced linearized system (4) controlled by the vision-based approach is given by

     ṗn (t) = A n p n (t) + B n u(t) u(t) = -k ⊤ p n (t) p n (t) = p n (t) + 1 N e n (t), ( 11
)
where k is the LQR optimal gain. The closed-loop reduced linearized system (11) can be rewritten as

ṗn (t) = (A n -B n k ⊤ )p n (t) - 1 N B n k ⊤ e n (t). ( 12
)
The initial value p n (t = 0) is therefore no longer of concerned in our approach. In addition for a large number of pixels N , the reduced linearized system dynamic equation ( 12) is less affected by measurements noise since 1 N B n k ⊤ e n (t) tends to 0. This is another great improvement over the LQG control scheme that is always noise dependent when noisy shear stress values are used in the LQE as shown in (10).

RESULTS

In this section the validations of our approach do not use real optical flow measurements. We use synthetic data sets of spatio-temporal variations of the perturbation velocities V p (x, y, t) obtained from the Poiseuille flow reduced linearized model presented in (4). This reduced linearized model is used to build a linear simulator of 2D channel Poiseuille flow with the following classical characteristics as used in related works [START_REF] Joshi | A systems theory approach to the feedback stabilization on infinitesimal and finite amplitude disturbances in plane Poiseuille flow[END_REF], [START_REF] Bewley | Optimal and robust control and estimation of linear paths to transition[END_REF] and [START_REF] Mckernan | Control of plane Poiseuille flow: a theoretical and computational investigation[END_REF]: the Reynolds number is R e = 10 000, the length of the channel is L c = 4π, and the reduced model Fourier wavenumber that contains instabilities is α n = 1. Using Matlab codes provided in [START_REF] Mckernan | Control of plane Poiseuille flow: a theoretical and computational investigation[END_REF], matrices A n , B n , C n given in (4) are computed.

The optimal gain k is computed from the weighting matrix R set as a scaled identity matrix R = rI with r = 200, penalizing thus high control values; and from the weighting matrix Q developed to compute the kinetic energy density E(t) = p n′ (t)Qp n (t) in [START_REF] Mckernan | Control of plane Poiseuille flow: a theoretical and computational investigation[END_REF] .

We first present results concerning the vision-based control approach in the ideal case where there is no measurements noise. Fig. 5 shows the different steps in the control of the perturbed flow with N x = 252. where we can see an increase due to the pertubation growth in the case where the flow is not controlled; and then a decrease also towards 0 once the control law is applied. At this step, we can see that the final velocities profile given in Fig. 5(d) is very similar to the desired velocities profile in Fig. 5(a). Therefore, the vision-based approach performs as expected.

Comparison of the estimation methods

In this section we show that the vision-based state estimation provides better results than the LQE. We consider a perturbed flow that is not controlled. Results are given in Fig. 6 in terms of the square norm of the state vector instead of the more relevant 2M (normally greater than 40 for a more accurate reduced model ( 4)) components of the state vector for the sake of clarity and readability. Fig. 6(a) presents the ideal case where there is no measurements noise and no initialization error. From this figure we can see that both estimations perfectly correspond to the ground truth value of the state vector. Fig. 6(c) highlights the initialization and the asymptotic convergence issues in the LQE; these issues are not of concerned in the visionbased approach which provides the ground truth value of the state vector. This result confirmes that the optical flow based estimation performs better than the LQE from shear stress measurements in any case.

The strong robustness to noise of the vision-based state estimation is presented on the right column figures where the standard deviation (STD) σ of on the optical flow noise is purposely set to a value 10 times higher than the STD σ ss on the shear stress noise. These figures present an average over a large number of realizations of the stochastic noises. Using only N = 252 pixels on the streamwise direction, Fig 6(b) shows a poor vision-based state estimation compared to Fig. 6(d) where N = 501 pixels are used. Note that N = 501 is far less than the number of pixels available in real situations where the images size can be at least 1280 × 960 (N = 1280). Due to a large number of flow velocities provided by the optical flow, the new approach is very robust to noisy measurements.

Behavior of the closed-loop systems

The behavior of the closed-loop system is shown to be better with the visio-based control than with the LQG control. Results are presented in Fig. 7. Fig. 7(a) depicts the behavior of the control signal in the ideal case (no measurements noise, no initialization error). Fig. 7(b) depicts the behavior of the control signal when the initial value is set as p n (0) = 0 by default since the value of p n (0) is unknown. In this case we can see that the value of the control signal is 100 times higher than the ideal control signal case which includes the vision-based approach (compare the highest control signal values in (2006). In addition, as expected, the control signal (see Fig. 7(b)) takes more time to converge to 0 (3000 iterations compared to the VB-LQG approach). This leads to an energy consumption far much higher for the LQG control than for the vision-based control.

The figures in the second row present the control signals in presence of measurements noise. Fig. 7(c) pictures the case of LQG control where we can see that the control signal does not converge to zero. Although the noise STD has been set to a small value, σ ss = 0.03, the control signal is very noisy, which is not suitable for actuators. Finally, Fig. 7(d) illustrates the robustness of the vision-based control where the STD in the optical flow noise is 10 times higher than the STD in the shear stress noise: we can see from this last figure that the control converges, and we can also see that the larger the sample of flow particles velocities used the lesser the noise in the control signal.

CONCLUSION

In this paper we have proposed a new approach for fluid flows control. This approach is based on optical flow techniques, which are non-intrusive, robust to measurements noise and thus suitable for fluid flows control. Theoretical proofs has been presented in order to show the improvements provided by the vision-based approach over the commonly used output feedback LQG control in terms of the state vector estimation and flows control. Indeed the output feedback LQG limitations concern the initialization of the observer and measurements noise. The initialization issue is not of concerned in the vision-based approach which provides lower control signal values than the LQG scheme in the case of an initialization error. In addition the vision-based approach has been shown to be robust to measurements noise since a large number of flow velocities is available in real practical situations. Future works will focus on a real implementation of the proposed approach.

1 N ε hn k (t), (A.8) where ε hn k (t) = ε n k (t) -fu(y k ) fu(y1) ε n 1 (t) -f l (y k ) f l (yM ) ε n M (t) . The homogeneous coefficients vector is thus given by       0 V hn py (y 2 , t) . . .

V hn py (y M -1 , t) 0       =      0 V hn py (y 2 , t) . . . V hn py (y M -1 , t) 0      + 1 N      0 ε hn 2 (t) .
. .

ε hn M -1 (t) 0     
. (A.9)

In the last step, the result of projection of V h py (y k , t) onto Chebychev polynomials leads to Plugging (A.9) into (A.10) leads to

a n (t) = a n (t) + 1 N Σ -1   ε hn 3 (t) . . . ε hn M -2 (t)   , (A.11) since a n (t) = Σ -1   
V hn py (y 3 , t) . . . a n (t) q n u (t) q n l (t) and e n (t) =

        Σ -1   ε hn 3 (t) . . . ε hn M -2 (t)   1 fu(y 1 ) ε n 1 (t) 1 f l (y M ) ε n M (t)        
, it is easy to express the estimated coefficient p n (t) in terms of the ideal one and the noise as given in (9).

FigFig. 1 .

 1 Fig. 1. Steady state velocities profile of the 2D plane Poiseuille flow: h c is the channel half height, L c is the flow streamwise period length and M is a point in the flow. state velocities profile of the plane Poiseuille flow in a streamwise period of the infinite length channel. The non dimensionalized Navier-Stokes equation of this flow is given by    

Fig. 4 .

 4 Fig. 4. Visualization of the flow using a laser sheet which role is to enlighten the particles seeded in the fluid. by using optical flow measurements as shown in Fig. 4. Optical flow is the apparent velocity vector field representing the motion of photometric pattern (pixels luminance) in successive image sequences. Optical flow techniques can be used to estimate instantaneous velocities of a fluid flow from any image sequences as detailed in Heitz et al. (2010).

  Fig. 5(a) pictures the desired image of the flow corresponding to the steady state velocities profile; Fig. 5(b) shows the image of the flow

Fig. 5 .

 5 Fig. 5. Vision-based control (the color in the figures represents the vorticity map): (a) desired image of the flow, (b) initial image of the flow, the perturbation has grown k = 750, (c) controlled flow at k = 1047, (d) controlled flow at k = 1500, (e) control law versus frame iteration and (f) kinetic energy density versus frame iteration. just before the application of the vision-based control law where we can see that the flow has become turbulent. Figs. 5(c) and 5(d) show different steps of the controlled flow at arbitrary selected iteration numbers k = 1047 and k = 1500 respectively: the control at each selected instant is represented by green vertical arrows on the upper and the lower channel boundaries. The control law converges since it tends towards 0 as shown on Fig. 5(e). Moreover, Fig. 5(f) depicts the kinetic energy density of the flow perturbationwhere we can see an increase due to the pertubation growth in the case where the flow is not controlled; and then a decrease also towards 0 once the control law is applied. At this step, we can see that the final velocities profile given in Fig.5(d) is very similar to the desired velocities profile in Fig.5(a). Therefore, the vision-based approach performs as expected.

Fig. 6 .

 6 Fig. 6. Comparison of state vector estimations using shear stress and optical flow: (a) ideal case, (b) measurements noise, (c) LQE initialization, (d) measurements noise with a large number of pixels N .

Fig. 7 .

 7 Fig. 7. Comparison LQG and the vision-based control approaches: (a) ideal case, (b) initialization error in the LQG regulator, (c) measurements noise in the LQG regulator, (d) measurements noise in the visionbased control. Fig. 7(b) and 7(a)). This higher control signal value could lead to an unsuitable state trajectory which can cause the real non-linear system to diverge as shown in McKernan(2006). In addition, as expected, the control signal (see Fig.7(b)) takes more time to converge to 0 (3000 iterations compared to the VB-LQG approach). This leads to an energy consumption far much higher for the LQG control than for the vision-based control.

  where Σ = (Σ i (y k )) i=5,...,M, k=3,...,M -2 is the combination of Chebychev basis evaluated at Gauss-Lobatto Collocation points (see[START_REF] Mckernan | Control of plane Poiseuille flow: a theoretical and computational investigation[END_REF]).

Appendix A. VISION-BASED STATE ESTIMATION

In this worth mentioning that only the normal component V py (x, y, t) of the perturbation veloticy V p (x, y, t) is used to compute the state vector p n (t) given in (4) (see [START_REF] Mckernan | Control of plane Poiseuille flow: a theoretical and computational investigation[END_REF]). Here we show the three steps for estimation of p n (t) (given by ( 9)) from the optical flow field measurements V py (y, x, t) (transpose of V py (x, y, t)) in the general case where measurements are corrupted by noise:

V py (y k , x i , t) = V py (y k , x i , t) + ε ki (t) (A.1) where ε ki (t) is an idd Gaussian noise on the veloctiy value at each pixel with standard deviation (STD) σ of and 0 mean; k = 1, ..., M ; i = 1, ..., N .

We first compute the Fourier series coefficients of the measured image velocities V py at the wavenumber α n = 1. Indeed by multiplying V py = V py (y k , x i , t) k=1,...,M, i=1,...,N by the transpose of the vector 1 N e -jαnx1 • • • e -jαnxM , we obtain the Fourier series of the noisy velocities

Since the upper boundary condition is a sinusoidal function, its Fourier series coefficient is given by V n py (y 1 , t) = q n u (t)f u (y 1 ), (A.3) from which we deduce that

Similarly the coefficient q n l (t) for the lower boundary can be expressed as

In the second step we compute the homogeneous coefficients vector V hn py . These homogeneous coefficient measurements V hn py (y k , t) are given by the expression V h py (y k , t) = V py (y k , t) -q n u (t)f u (y) -q n l (t)f l (y). (A.7) Using (A.2), (A.5) and (A.6), it is easy to rewritte expression (A.7) as V hn py (y k , t) = V hn py (y k , t) +


Chao-Sheng Tang
email: tangchaosheng@nju.edu.cn

Yu-Jun Cui
email: yujun.cui@enpc.fr

Bin Shi
email: shibin@nju.edu.cn

Anh-Minh Tang
email: tang@cermes.enpc.fr

Chun Liu

Desiccation and cracking behaviour of clay layer from slurry state under wetting-drying cycles

Keywords: wetting-drying cycle, clay layer, desiccation crack, aggregate formation, image processing, crack pattern

Laboratory tests were conducted to investigate the effect of wetting-drying (W-D) cycles on the initiation and evolution of cracks in clay layer. Four identical slurry specimens were prepared and subjected to five subsequent W-D cycles. The water evaporation, surface cracks evolution and structure evolution during the W-D cycles were monitored. The effect of W-D cycles on the geometric characteristics of crack patterns was analyzed by image processing. The results show that the desiccation and cracking behaviour was significantly affected by the applied W-D cycles: the measured cracking water content θ c , surface crack ratio R sc and final thickness h f of the specimen increased significantly in the first three W-D cycles and then tended to reach equilibrium; the formed crack patterns after the second W-D cycle were more irregular than that after the first W-D cycle; the increase of surface cracks was accompanied by the decrease of pore volume shrinkage during drying. In addition, it was found that the applied W-D cycles resulted in significant rearrangement of specimen structure: the initially homogeneous and non-aggregated structure was converted to a clear aggregated-structure with obvious inter-aggregate pores after the second W-D cycle; the specimen volume generally increased with increasing cycles due to the aggregation and increased porosity. The image analysis results show that the geometric characteristics of crack pattern were significantly influenced by the W-D cycles, but this influence was reduced after the third cycle. This is consistent with the observations over the experiment, and indicates that the image processing can be used for quantitatively analyzing the W-D cycle dependence of clay desiccation cracking behaviour.

Introduction

The formation of desiccation cracks on soil surface due to loss of water is a common natural phenomenon, and can significantly affect the soil performance in various geotechnical, agricultural and environmental applications. For example, a cracked soil is more compressible than an intact one at the same water content and the overall mechanical strength is weakened due to the presence of cracks [START_REF] Morris | Cracking in drying soils[END_REF].

The size (width, length and depth), tortuosity, spatial distribution and connectivity of cracks govern the rate and the velocity at which solutes and microorganisms are transported in the soil, and thus control the dispersal of substances in soil [START_REF] Horgan | An empirical stochastic model for the geometry of two-demensional crack growth in soil (with discussion)[END_REF]. Most importantly, the soil hydraulic properties are directly controlled by the desiccation crack networks [START_REF] Chertkov | Tortuosity of crack networks in swelling clay soils[END_REF][START_REF] Chertkov | Using surface cracks spacing to predict crack network geometry in swelling soils[END_REF].

Many previous studies have indicated that the hydraulic conductivity of cracked soils is several orders of magnitude greater than that of intact soils [START_REF] Boynton | Hydraulic conductivity tests on compacted clay[END_REF][START_REF] Albrecht | Effect of desiccation on compacted natural clays[END_REF]. This issue is therefore a major concern in design and construction of low permeability structures as clay buffers and barriers for nuclear waste isolation, liners and covers for landfill, etc.

Over the past decades, a number of field studies and laboratory experiments have been undertaken to investigate the initiation and propagation of desiccation cracks in soils [START_REF] Corte | Experimental research on desiccation cracks in soil[END_REF][START_REF] Kleppe | Desiccation cracking of soil barriers. Hydraulic barriers in soil and rock[END_REF][START_REF] Morris | Cracking in drying soils[END_REF][START_REF] Konrad | An idealized framework for the analysis of cohesive soils undergoing desiccation[END_REF][START_REF] Miller | Experimental analysis of desiccation crack propagation in clay liners[END_REF][START_REF] Velde | Structure of surface cracks in soil and muds[END_REF][START_REF] Nahlawi | Laboratory experiments on desiccation cracking of thin soil layers[END_REF][START_REF] Tang | Influencing factors of geometrical structure of surface shrinkage cracks in clayey soils[END_REF][START_REF] Tang | Experimental evidence on the temperature dependence of desiccation cracking behavior of clayey soils[END_REF]. However, these investigations have been largely qualitative and most are limited in the description of desiccation cracking phenomena.

More recently, techniques for quantifying the main features of the crack patterns have evolved from direct field measurement to more sophisticated analysis by image processing [START_REF] Miller | Experimental analysis of desiccation crack propagation in clay liners[END_REF][START_REF] Velde | Structure of surface cracks in soil and muds[END_REF]Vogel et al., 2005). Image analysis has proved to be a powerful tool by which several geometric and morphologic parameters such as crack width, length, area, angle and their distribution characteristics can be determined effectively. In addition, some modelling and theoretical studies on desiccation cracking have also been conducted [START_REF] Abu-Hejleh | Desiccation Theory for Soft Cohesive Soils[END_REF][START_REF] Ayad | Desiccation of a sensitive clay: application of the model CRACK[END_REF][START_REF] Konrad | An idealized framework for the analysis of cohesive soils undergoing desiccation[END_REF][START_REF] Chertkov | Modeling the crack network of swelling clay soils[END_REF][START_REF] Chertkov | Using surface cracks spacing to predict crack network geometry in swelling soils[END_REF][START_REF] Chertkov | Modelling cracking stages of saturated soils as they dry and shrink[END_REF][START_REF] Deng | Numerical simulation of crack formation process in clays during drying and wetting[END_REF][START_REF] Péron | Desiccation cracking of soils[END_REF]Péron et al., 2009a).

However, as soil is a highly complex material, the desiccation cracking behavior is governed by a large number of factors including mineral composition, clay content, relative humidity, temperature, layer thickness, boundary conditions etc. [START_REF] Fang | Introduction to environmental geotechnology[END_REF][START_REF] Albrecht | Effect of desiccation on compacted natural clays[END_REF][START_REF] Nahlawi | Laboratory experiments on desiccation cracking of thin soil layers[END_REF][START_REF] Rodríguez | Experimental and numerical analysis of desiccation of a mining waste[END_REF][START_REF] Tang | Developing law and morphological analysis of shrinkage cracks of clayey soil under different temperatures[END_REF][START_REF] Tang | Influencing factors of geometrical structure of surface shrinkage cracks in clayey soils[END_REF][START_REF] Tang | Experimental evidence on the temperature dependence of desiccation cracking behavior of clayey soils[END_REF]. It is therefore difficult to propose a rational model to describe this complex phenomenon with a reasonable number of parameters. The essential mechanism of desiccation cracking is still not well understood today and the prediction of cracks initiation and the associated crack network propagation also faces several challenges.

It is recognised that the soil in-situ is subject to diurnal changes and seasonal rainy and sunny weather, and undergoes periodical wetting-drying (W-D) cycles. Several studies have been performed to investigate the effect of wetting and drying on soil physical properties [START_REF] Rao | The impact of cyclic wetting and drying on the swelling behavior of stabilized expansive soils[END_REF][START_REF] Alonso | Expansive bentonite-sand mixtures in cyclic con-trolled suction drying and wetting[END_REF][START_REF] Nowamooz | Experimental and numerical studies of the hydromechanical behaviour of a natural unsaturated swelling soil[END_REF][START_REF] Tang | Effects of the maximum soil aggregates size and cyclic wetting-drying on the stiffness of a lime-treated clayey soil[END_REF]. However, the relationship between W-D cycles and desiccation cracking behaviour is still not clearly understood. In this study, the effect of W-D cycles on cracking behaviour was investigated by subjecting a group of initially saturated clay layers to five W-D cycles. The initiation and evolution of cracks on the specimen surface during the cycles were monitored. The geometric characteristics of surface crack patterns are described and quantified through image processing and the mechanisms involved are discussed.

Material and methods

Material

The Romainville clay was used in this investigation. The physical properties are presented in Table 1. This clay is a lagoonal-marine deposit which is part of the Paris Basin Tertiary (Oligocene) formations. It is widely distributed over the eastern region of Paris, and has been considered responsible for the large amount of damages to buildings due to the swelling-shrinkage and cracking phenomenon. Various studies have been undertaken to analyse the hydro-mechanical behaviour of the clay under weather effects such as W-D cycles [START_REF] Audiguier | Caractérisation au laboratoire de la sensibilité au retrait-gonflement des sols argileux[END_REF][START_REF] Laribi | Assessing shrink/swell properties of two argillaceous soils from the Paris Basin: a comparison of cation exchange determination methods[END_REF][START_REF] Geremew | Analysis of the behaviour of a natural expansive soil under cyclic drying and wetting[END_REF]Tang et al.,2009a;Tang et al., 2009b). However, little attention has been paid to the effect of W-D cycles on the desiccation cracking behaviour.

Test methods

The air-dried Romainville clay was crushed and sieved at 2 mm. Saturated slurry specimens were prepared by mixing the crushed powder with distilled water at a water content of approximately 170 % (g g -1 ) (Note that all the given water contents of the tested specimens in the next sections is 'gravimetric' water content, unless otherwise specified). A desired quantity of slurry was then poured into glass cups (117 mm in diameter). To remove entrapped air bubbles in the slurry, these cups were placed on a vibration device for 5 minutes. Finally, the cups were sealed with plastic membrane and left for at least 3 days. The final settled layer thickness was about 10 mm. Four parallel specimens were prepared and dried at constant room temperature (25 ± 1 °C, 50 ± 5 % of relative humidity) until the weight of specimen was stabilized, i.e., the first W-D cycle was completed. The subsequent wetting was started by pouring distilled water directly into the glass cups. During this wetting process, sufficient water was provided to ensure full submergence of the specimen and no mixing was applied. The glass cups were again sealed with plastic membrane to prevent water evaporation. After three days, the specimens were exposed to room conditions to be dried again. This procedure was repeated and finally a total of five W-D cycles was applied.

The schematic set-up used in this study is illustrated in Fig. 1. In order to measure the variation of water content during drying, the specimen was placed on a balance (with an accuracy of 0.01 g) to monitor the weight. A digital camera was fixed above the specimens to capture the surface image during wetting and drying cycles in different time steps. Fig. 2 presents a typical crack pattern taken during the first drying path. With the application of image processing, the surface crack ratio (R sc ) or crack intensity factor CIF [START_REF] Miller | Experimental analysis of desiccation crack propagation in clay liners[END_REF], defined as the ratio of the area of cracks to the total initial area of specimen, was determined to quantify the cracking extent during drying. At the end of each W-D cycle, other geometric parameters -the number of intersections per unit area N int , number of crack segments per unit area N seg , mean crack length L av , mean crack width W av and mean clod area A av -were also determined from the final crack patterns. Note that the crack length was determined by calculating the length of the mid-axis of crack segment between two intersections.

For this purpose, a skeletonising operation [START_REF] Gonzalez | Digital Image Processing Second Edition[END_REF] was initially performed by repeatedly removing pixels from the boundaries of cracks until they are reduced to single-pixel-wide skeletons. The crack width was sampled in orthogonal direction and the shortest distance between the crack boundaries was calculated. The clods are surrounded by cracks, and their areas are defined as the pixel number of the regions. All these parameters were obtained by applying the software CIAS which was developed by our research group. More details can be found in [START_REF] Liu | The analysis method of Morphological parameters of rock and soil crack based on image processing and recognition[END_REF] and [START_REF] Tang | Influencing factors of geometrical structure of surface shrinkage cracks in clayey soils[END_REF]. Eight clods were selected from each specimen in the end of each drying period to determine their thickness using calliper, and the mean value is regarded as the final thickness of specimen layer.

In order to investigate the volume shrinkage behaviour of specimens during drying, four other identical specimens were prepared following the same procedure described above. During the first drying path, small clods were taken from the specimens in different time steps to determine their water contents θ and density ρ. The clod volume used for the density determination was measured by immersing it in a non-wetting hydrocarbon liquid. Then the void ratio e and the corresponding degree of saturation S r of specimens can be determined using the following equations:

1 ) 01 . 0 1 ( - + = ρ θ ρ s e (1) e S s r θρ = (2)
where ρ s is the density of the solid phase of the Romainville clay, equal to 2.79 Mg m -3 (Table 1).

Results

Evaporation, shrinkage and cracking process

The measured water content θ at various times t during the first drying path for the four separate specimens are shown in Fig. 3 (desiccation curve). Two distinct evaporation stages can be indentified: a constant evaporation stage during which water content decreases linearly with time; and a subsequent falling evaporation stage during which water loss slows down gradually until the residual water content of about 4.3 % is reached.

The determined shrinkage curve (e versus θ) is presented in Fig. 4, from which two transition points between linear and exponential phases can be determined -the air entry (AE) point and the shrinkage limit (SL) point [START_REF] Stirk | Some aspects of soil shrinkage and the effect of cracking upon water entry into the soil[END_REF][START_REF] Boivin | Relationship between clay content, clay type, and shrinkage properties of soil samples[END_REF][START_REF] Groenevelt | Analysis of soil shrinkage data[END_REF][START_REF] Cornelis | Measuring and modeling the soil shrinkage characteristic curve[END_REF][START_REF] Chertkov | The reference shrinkage curve of clay soil[END_REF]. The water contents at the point of AE and SL are approximately 18% and 12%, respectively.

After desiccation cracks appear on the surface of the specimen, the surface crack ratio R sc at different water content θ was determined by image processing, and the result is also plotted in Fig. 4 as the cracking curve (R sc versus θ). This indicates that the water content at the onset of cracking ranges from 38 to 43%, which is much higher than the AE value. Therefore the specimens are still fully saturated as cracking occurs. During the initial stage, the R sc increases rapidly with decreasing water content and decreasing pore volume. However, with further drying, the increment of R sc slows down once the water content is close to AE point (θ = 18%). After the SL (θ = 12%) is reached, the R sc approaches a relative steady value of 14.6-15.8 % for the four separate specimens.

Cracking water content, final surface crack ratio and layer thickness

During drying, when the first crack is observed on the specimen surface, the corresponding water content is herein defined as the cracking water content θ c . Fig. 5 shows the determined θ c of the four specimens during each drying path. It indicates that θ c increases rapidly during the first three drying paths, and does not change significantly during the subsequent fourth and fifth drying paths. For instance, the mean value of θ c of the four specimens during the first drying path is 40.5%; increasing to 76.1% during the third drying path, whereas the mean value of θ c during the fifth drying path is 80.2%, only 1.2% higher than the fourth drying path.

After each W-D cycle, the final surface crack ratio R sc and the mean final thickness h f were also determined and are presented in Fig. 6 and Fig. 7 respectively. Similar to θ c , the final mean value of R sc and h f of the four specimens also increase quickly in the first three W-D cycles, while the increment rate slows down in the subsequent cycles.

Geometric parameters of crack pattern

Quantitative analysis of the crack pattern is important when studying clay cracking. Indeed, the geometric parameters are helpful in evaluating the hydro-mechanical properties of the clay-water system. They are also related to the variation of the stress-strain state and may provide a way to investigate the essential mechanisms of cracking. In addition, the crack propagation behavior and the geometric characteristics of the crack pattern reflect the material plasticity and mineral compositions. [START_REF] Perrier | Computer construction of fractal soil structures: simulation of their hydraulic and shrinkage properties[END_REF] indicated that if the real structure features of crack pattern can be determined, the soil response to wetting and drying can be predicted.

Table 2 summarizes the mean values of the geometric parameters which were determined from the final crack patterns after each W-D cycle, the corresponding standard deviation (SD) are also presented. It is found that the mean values of all the parameters in the first three cycles are much more variable than that in the subsequent cycles. For N int and N seg , they reach the maximum value after the second cycle and decrease after the third cycle. After that, N int and N seg increase slightly with further increase of W-D cycles. Moreover, the ratio of N seg /N int is greater than 1.5 and less than 2 for all the cycles, indicating that, three or four crack segments sharing one intersection, is the dominant outline of the crack pattern. For L av , the maximum and the minimum values are 8.43 mm and 7.57 mm respectively after the first and second cycle. It is between 8.07 and 8.22 mm after the third cycle and no obvious trend can be observed. For W av , the minimum value is 1.01 mm that is observed after the first cycle; it reaches a maximum value at 1.86 mm after the third cycle and decreases slightly in the subsequent cycles. Based on the obtained values of L av and W av , it can be deduced that most of the formed cracks in the first drying cycle are slim and threadlike. For A av , the maximum value is 1.18 cm 2 observed after the first cycle. It generally decreases with increasing cycles and finally stabilizes at 0.65 cm 2 , indicating that the clods become smaller due to fragmentation. These results indicate that the geometric and morphologic characteristics of crack patterns are significantly influenced by W-D cycles, but the influence intensity declines after the third cycle.

Moreover, the corresponding SD of the determined parameters generally decreases with increasing cycles, indicating that the initial densities, sizes or dimensions of the crack elements show significant variability but become more homogeneous under the effects of W-D cycles. 2009b), and can be explained using the maximum stress release criterion and crack propagation criterion [START_REF] Lachenbruch | Mechanics of thermal contraction cracks and ice-wedge polygons in permafrost[END_REF][START_REF] Morris | Cracking in drying soils[END_REF].

Visual observations and discussion

Crack pattern and structure evolution after W-D cycles

After the second W-D cycle (Fig. 8 (b)), the shapes of the separated clods are more irregular and the crack segments are more jagged than that observed in Fig. 8 2. In the tests of [START_REF] Yesiller | Desiccation and cracking behavior of three compacted landfill liner soils[END_REF] and [START_REF] Tang | Influencing factors of geometrical structure of surface shrinkage cracks in clayey soils[END_REF], most cracks after the second and third W-D cycles were found to be situated at the same locations as that after the first W-D cycle. However, this was not evident in this study when comparing Figs. 8 (a), (b) and (c). This can be attributed to the heterogeneous aggregated-structure formed during the wetting process. Many weak zones were created, which significantly control the start positions of cracks during the subsequent drying period [START_REF] Yong | Soil properties and behavior[END_REF].

Based on the results of surface crack ratio R sc and layer thickness h f shown in Figs. 6 and7, the final total volume of the dried clods (without the cracks) as well as the increment of porosity after each W-D cycle can be calculated. It is found that, after the fifth cycle, the porosity of the clods increases by approximately 21%. This is consistent with the findings of the previous studies on remolded and natural soils [START_REF] Sartori | Changes in soil porosity and surface shrinkage in a remolded, saline clay soil treated with compost[END_REF][START_REF] Pires | Micromorphological analysis to characterize structure modifications of soil samples submitted to wetting and drying cycles[END_REF], and can be explained by the progressive increase in pore volume and the average diameter of the pores with the increase of W-D cycles [START_REF] Geremew | Analysis of the behaviour of a natural expansive soil under cyclic drying and wetting[END_REF].

Crack pattern evolution and aggregate formation during wetting

A time series of 6 images is shown in Fig. 9 to illustrate the evolution of surface cracks and aggregate formation process during the second wetting path. For reasons of simplicity and clarity, only the white marked area in Fig. 8 Notably, the collapse or destruction firstly occurred on the clods edges, especially on the corner positions (circled zones in Fig. 9 (b)). This is because rapid hydration firstly occurred at these less stable positions, and the previous cracks provided free spaces for clay swelling.

Upon further wetting, the original desiccation cracks tend to be increasingly narrow due to the filling of the exfoliated aggregates from the clods as well as the volume swelling of the clods. Fig. 9 (d) shows that the original desiccation cracks were fully closed after about 2.5 min wetting. However, there were more and more new micro-cracks induced by wetting. These new cracks are mainly resulted from the differential swelling pressures and the inner stresses during the wetting process of the clods. Fig. 9 (d) also shows that the specimen surface is not as flat and smooth as its initial state, but is 'heaved' along the trace of the original desiccation crack segments, as quadrangle marked in Figs. 9 (c)-(e). This is because the high swelling potential of the clay minerals was constrained so that the clods jostle against each other.

Comparison between Fig. 9 (e) and Fig. 9 (f) shows that the surface morphologies are similar to each other. The edges of the aggregates are clear and the wetting-induced new cracks can still be identified even after a long time (72 hours) has elapsed. However, in the test of [START_REF] Tang | Influencing factors of geometrical structure of surface shrinkage cracks in clayey soils[END_REF], the aggregates breakdown into individual fine particles and the new cracks were completely closed after a short wetting time of 2 hours. This can be attributed to the different clay fractions in the studied materials. For Romainville clay, the clay fraction is 79 %, which is 57% higher than Xiashu clayey soil used by [START_REF] Tang | Influencing factors of geometrical structure of surface shrinkage cracks in clayey soils[END_REF]. The dispersion of aggregates was therefore prevented by the strong bonds between clay particles.

Pictures of the specimen subjected to the third wetting path are presented in Fig. 10.

The observations are slightly different from the second wetting process (Fig. 9). After water was poured into the glass cup, clods degraded rapidly and some small aggregate plates floated on the water surface (Fig. 10 (a)). About 1.5 min later, the original desiccation cracks were fully closed (Fig. 10 (b)); after a period of 72 h wetting, all the clods were broken down and a large amount of individual aggregates can be observed (Fig. 10 (c)). However, no new micro-cracks or "jostle effects" were observed, as in the second wetting path (Fig. 9).

During the fourth and fifth wetting paths, no significant difference was observed as compared to the third wetting path, and as a result, no pictures are presented.

Effect of W-D cycles on cracking water content

Fig. 5 demonstrates that the W-D cycles result in significant increase in cracking water content θ c . Moreover, near 90% of the increase occurs in the first three cycles.

Previous studies indicate that cracking occurs when the tensile stress induced by the developed suction exceeds the tensile strength of the specimen [START_REF] Corte | Experimental research on desiccation cracks in soil[END_REF][START_REF] Morris | Cracking in drying soils[END_REF][START_REF] Miller | Experimental analysis of desiccation crack propagation in clay liners[END_REF]Péron et al., 2009b). The effects of W-D cycles on this cracking behaviour can be explained by considering the specimen fabric, as follows: i) during W-D cycles, volume shrinkage and cracking would result in irreversible fabric changes, decrease of specimen integrity and increase of weak zones in specimen [START_REF] Yong | Soil properties and behavior[END_REF]. As a result, the specimen tensile strength was reduced, promoting specimen cracking at a higher θ c upon the subsequent drying;

ii) the multiple W-D cycles gave rise to an increase in material heterogeneity as discussed in the above sections. For instance, after the first W-D cycle, the initial homogenous structure was changed to aggregated-structure during the subsequent W-D cycles (Fig. 8). Due to the non-uniform distribution characteristics of aggregate sizes and inter-aggregate pore sizes, the tensile stresses developed during subsequent drying cycles would be also non-uniform and can easily be concentrated at defects with lower tensile strength. [START_REF] Weinberger | Initiation and growth of cracks during desiccation of stratified muddy sediments[END_REF] and [START_REF] Tang | Influencing factors of geometrical structure of surface shrinkage cracks in clayey soils[END_REF] have observed that cracks first initiated at surface defects, and emphasized that the surface defects can trigger or promote crack initiation due to shrinkage distortion of the surrounding medium and tensile stress concentration.

Note that because of the difficulty in drying tests, there has been no experimental data concerning the actual value of the local water content at the moment of crack initiation. The commonly measured water content in the laboratory is a mean value and can not completely reflect the moisture distribution characteristics in localised positions. To reduce the test error, the simplest way is to prepare an initially relatively homogenous specimen with a limited thickness -thus thin layer (10 mm thick) slurry specimens were prepared in this study. Generally, the physical, hydraulic and mechanical properties of the specimen vary permanently with moisture changes during drying. In order to better understand the intrinsic mechanism of desiccation cracking , future studies should also take into account some dynamic factors such as evaporation rate, suction development rate, volume shrinkage rate and tensile strength development rate, in addition to the critical water content θ c .

Equilibrium state during multiple W-D cycles

Based on all the results shown in Figs. 5 (θ c ), 6 (R sc ), 7 (h f ) and Table 2, it can be concluded that the cracking behaviour reached a relative equilibrium state after the specimens were subjected to the third cycle. This conclusion is also supported by the observations shown Fig. 8: no significant difference can be observed for the crack patterns after the third, fourth and fifth W-D cycle. Presumably the rearrangement of specimen fabric diminishes and finally ceases after certain a number of W-D cycles.

Previous studies on this subject also confirmed that there is an equilibrium state during multiple W-D cycles, but the corresponding numbers of cycles depend on material nature especially the clay fraction. Generally, clayey soils with high plasticity need more W-D cycles to reach the equilibrium state than silty or sandy soils (Al-Wahab and EI-Kedrah, 1995;[START_REF] Yesiller | Desiccation and cracking behavior of three compacted landfill liner soils[END_REF][START_REF] Omidi | Effect of desiccation cracking on the hydraulic conductivity of a compacted clay liner[END_REF].

Conclusions

Desiccation and cracking behaviour of clay layers from a slurry state upon five wetting-drying (W-D) cycles were investigated through laboratory experiments. The process of water evaporation, surface cracks evolution, structure evolution and volume shrinkage behaviour were monitored and have been discussed here. The geometric characteristics of crack pattern after each W-D cycle were quantitatively analyzed by image processing. The following conclusions can be drawn:

(1) during the first drying path, the water evaporation process is composed of two stages: a constant rate zone and a subsequent decreasing rate zone; the increase of surface crack ratio R sc during drying was accompanied by progressive pore volume shrinkage, and the R sc reached stabilization as the water content reached the shrinkage limit; the final crack pattern was dominated by polygonal clods and smooth crack networks;

(2) during the second wetting path, the poured water resulted in rapid collapse of clods, and the desiccation cracks formed in the previous drying path were quickly closed; meanwhile, a large number of new micro-cracks induced by wetting appeared on the specimen surface that divided the clods into smaller aggregates, and a typical aggregated-structure was developed and significantly intensified the specimen heterogeneity; moreover, the second wetting path led to significant rearrangement of clay particles and modification of the pore network; these processes are generally irreversible and drastically influence the desiccation cracking behaviour during the subsequent drying path;

(3) after the second drying path, since the specimen homogeneity decreased, the shapes of the clods were more irregular and the crack segments were more jagged than that observed after the first drying path; the clods showed clear aggregated-structure with a large amount of inter-aggregate pores; the bonds between the formed aggregates were broken quickly in the third wetting path and no new micro-cracks were induced;

(4) an equilibrium state can be reached after a certain numbers of W-D cycles; for the tested Romainville clay, it was observed that the rearrangement of specimen fabric or structure reached an apparent state of equilibrium after the third W-D cycle: before that, the measured cracking water content θ c , final surface crack ratio R sc and layer thickness h f increased significantly with increasing cycles, but after, the effects of W-D cycles on the cracking behaviour and the cracking pattern became insignificant;

(5) image processing provided useful information on the geometric characteristics of crack patterns. The effect of W-D cycle on desiccation and cracking behaviour can therefore be quantitatively described through these specific parameters.

Generally, field soils are very complex and conditioned by a large number of variables. Cracking is a 3D phenomenon and the patterns are usually composed of irregular polygonal shapes. Quantification of a 3D crack pattern is a much more challenging task than the 2D process described in this investigation. However, it is believed that the findings of this laboratory experiment are helpful for better understanding the interrelationship between the sequence of seasonal processes and the evolution of soils or mud properties. In addition, as multiple W-D cycles can result in possible changes in clay structure and cracking behaviour, the hydraulic and mechanical properties of this type of materials must depend on the number of W-D cycles. Further attention should therefore be paid to the engineering fields where swelling clays are involved.

Tables

Figure captions

Fig. 8

 8 Fig. 8 presents the typical crack pattern after each W-D cycle. It indicates that the specimen surface was split to separate clods by the crack networks. After the first W-D cycle (Fig.8 (a)), the shapes of the clods are relatively regular and most of the clods are close to quadrangles or pentagons; the crack segments are smooth and generally perpendicular to each other. This is consistent with the observations of Vogel et al. (2005) and Peron et al. (2009b), and can be explained using the maximum

  Figs. 8 (c), (d) and (e) indicate that the geometric and morphologic characteristics of the crack patterns are similar to each other after the third, fourth and fifth W-D

  (a) is presented. Fig. 9 (a) shows the initial morphology of desiccation crack pattern taken from Fig. 8 (a). Figs. 9 (b) and (c) indicate that injecting water led to immediate collapse of the clods: the clods were broken down into several small aggregates. Meanwhile, some new micro-cracks induced by wetting appeared on the specimen's surface.

Fig. 1

 1 Fig. 1 Schematic drawing of set-up

  Balance

Fig. 1 Fig. 3 Fig. 4 Fig. 7

 1347 Fig. 1 Schematic drawing of set-up

Table 1

 1 Physical properties of Romainville clay

	Physical properties	Values
	Density of the solid phase	2.79 Mg m -3
	Liquid limit	77 %
	Plastic limit	40 %
	Plasticity index	37 %
	USUC a classification	CH
	Clay (<2µm)	79 %
	Clay composition	Illite and smectite
	Specific surface area	340 m 2 /g
	a-Unified Soil Classification System

Table 2

 2 Mean values of crack quantitative parameters and the corresponding standard deviation (SD) after each wetting-drying cycle

	After each	N int	SD of	N seg	SD of	L av	SD of	W av	SD of	A av	SD of
	wetting-drying cycle	(cm -2 )	N int	(cm -2 )	N seg	(mm)	L av	(mm)	L av	(cm 2 )	L av
	First	1.09	0.17	1.78	0.12	8.43	4.50	1.01	0.63	1.18	1.20
	Second	1.52	0.12	2.45	0.19	7.57	3.75	1.52	0.49	0.78	0.56
	Third	1.37	0.08	2.13	0.11	8.12	3.42	1.86	0.42	0.69	0.53
	Fourth	1.43	0.05	2.16	0.08	8.07	3.31	1.74	0.38	0.65	0.53
	Fifth	1.44	0.07	2.22	0.05	8.22	2.75	1.73	0.37	0.65	0.49

Acknowledgements

This work was supported by the National Natural Science of China (Grant No. 41072211), Research Fund for the Doctoral Program of Higher Education of China (Grant No. 20090091120037), State Key Program of National Natural Science of China (Grant No. 40730739) and the College Graduate Student Innovation Program of Jiangsu Province (CX09B_011Z). It was performed within the project ANR-RGCU «Analyse du Retrait-Gonflement et de ses Incidences sur les Constructions » -ARGIC. All the tests were performed at UR Navier/CERMES, Ecole des Ponts-ParisTech.

Thanks to the two anonymous reviewers for their constructive comments that helped us to improve the manuscript.


Romeo Tatsambon Fomena

Christophe Collewet Cemagref
email: christophe.collewet@cemagref.fr

VISION-BASED CONTROL OF 2D PLANE POISEUILLE FLOW

This paper introduces a new approach for active fluid flows control: the vision-based approach. By using vision, dense flow velocity maps can be estimated and used in an observer-free closed-loop scheme to control a flow. This new approach is validated on a 2D plane Poiseuille flow and is proven to outperform the existing Poiseuille flow control appraches which use a limited number of available shear stress measurements.

INTRODUCTION

A significant part of the work carried out in the field of flows control has been dedicated to the control of the transition from laminar to turbulent states. Delaying, accelerating or modifying this transition can be of great environmental and economical interests for industrial applications. For instance, a drag reduction can be obtained while enhancing the lift, leading consequently to limit fuel consumption for aircrafts.

Flows control can be achieved in two different ways: passive or active control. Passive control provides a permanent action on a system. Most often it consists in optimizing shapes or in choosing suitable surfacing. Conversely, in active control an external energy is required to act on the system like, for example, techniques based on blowing and suction [START_REF] Joshi | A systems theory approach to the feedback stabilization on infinitesimal and finite amplitude disturbances in plane Poiseuille flow[END_REF]). This approach can be seen as an optimal problem where one has to apply an optimal control law based on a certain cost (minimization of the drag, minimization of the actuators power, etc.) [START_REF] Bewley | Optimal and robust control and estimation of linear paths to transition[END_REF]). However, very often, only open-loop control and forcing are used [START_REF] Protas | Drag force in the open-loop control of the cylinder wake in the laminar regime[END_REF]). As a matter of fact, designing a closedloop control law requires the use of sensors that can be at the same time non-intrusive, accurate and adapted to the time and space scale of the phenomenon under monitoring. Unfortunately, such sensors are hardly available in the real context of control applications. The most commonly used measurement, obtained from MEMS, is the shear stress at a limited set of measurement points of the wall [START_REF] Joshi | A systems theory approach to the feedback stabilization on infinitesimal and finite amplitude disturbances in plane Poiseuille flow[END_REF]).

To deal with these issues, this paper proposes a visionbased control approach. As far as we know, such an approach has never been used for flows control issues. Vision-based control is now a well established technique in the robotics and automatic control communities. This technique has shown impressive results in numerous complex contexts such as underwater, medical and aeral robotic [START_REF] Bonin-Font | Visual navigation for mobile robots: a survey[END_REF]). This technique consists in using feedback information provided by a vision sensor to control a dynamic system [START_REF] Espiau | A new approach to visual servoing in robotics[END_REF]). By using vision, the complete flow can therefore be fully observed in real time [START_REF] Champagnat | How to calculate dense piv vector fields at video rate[END_REF]). This full observation is shown here to be of great improvements for flows control in comparison with existing control approaches that use only a limited set of walls shear stress observations (measurements).

In this paper we use the vision-based control approach to regulate the 2D plane Poiseuille flow around its steady state. In the next section we recall the modeling for the stabilization of this flow. In Section 3 we clearly define the unstable Poiseuille flow problem and we present the existing control solutions which include the Linear Quadratic Gaussian (LQG) regulator. The limitations of these solutions are clearly highlighted in the same section. In Section 4 the vision-based control approach is presented and theoretically proven to overcome the aformentioned limitations in Section 5. Since the LQG approach is the standard effective approach for flows control, In Section 6 we compare the vision-based and the LQG approaches. Simulations are given in Section 7 to validate the proposed approach.

MODELING OF 2D PLANE POISEUILLE FLOW FOR CONTROL

In this section we first present the basics of plane Poiseuille flow, then we recall the boundary control principle for this flow and finally we present the reduced linearized model used to derive the currently existing shear stress based control laws.

Basics

       ∂ V ∂t + (V • ∇)V = -∇P + 1 R e ∇ 2 V ∇ • V = 0 V(x, y = ±1,t) = 0 (1)
where P is the pressure; V is the flow velocity; V(x, y = ±1,t) = 0 represents the no slip boundary condition and R e is the Reynolds number.

Since Poiseuille flow is simple, the analytical solution (V bx ,V by , P b ) of (1) in the steady state case, i.e.

dV dt = ∂ V ∂t + (V • ∇)V = 0, can be found: (V bx ,V by , P b ) = (1 -y 2 , 0, - 2 R e x).
(

) 2
This solution is illustrated on Fig. 1.

Control principle

A perturbed plane Poiseuille flow can be controlled via boundaries. Boundary control consists in modifying boundaries conditions as shown in [START_REF] Joshi | A systems theory approach to the feedback stabilization on infinitesimal and finite amplitude disturbances in plane Poiseuille flow[END_REF]).

From a physical point of view, boundary control can be interpreted as a geometric alteration of the boundary as pictured on Fig. 2. The boundary control on the upper and the lower channels can be theoretically modeled by functions χ u and χ l that ensure the mass conservation in the controlled system as pictured on Fig. 2. Note that in the absence of control, i.e. when χ u = χ l = 0, the red dashed curves (see Fig. 2) are aligned with the lower and upper boundary lines as expected.

Reduced order model of the unstable 2D plane Poiseuille flow

For a pratical implementation of flow control methods, the infinite dimension of a flow prompts the need of a reduced flow model. This section aims at deriving the reduced model of a controlled Poiseuille flow.

Concerning Poiseuille flow, most of the works focus on temporal instabilities caused by a perturbation velocity V p (x, y,t). In order to keep permanent such instabilities in the infinite channel when the flow is not controlled, a periodic boundary finite length channel is assumed [START_REF] Joshi | A systems theory approach to the feedback stabilization on infinitesimal and finite amplitude disturbances in plane Poiseuille flow[END_REF]). That is why the perturbation velocity V p (x, y,t) can be expanded in a Fourier series

V p (x, y,t) = +∞ ∑ n=-∞ V n p (y,t)e iα n ( 3
)
where α n = nα 0 , with α 0 = 2π L x the fundamental wavenumber and L x is the streamwise period length (see Fig. 2).

The modeling, required to derive shear stress based control laws, consists first of all in linearizing the NSE around the steady state solution (2). Then the continuous linearized model of the NSE is reduced by approximation of the perturbation velocity V p (x, y,t) at a specifically selected wavenumber α n of the Fourier series (3); and by decomposition of the Fourier series coefficient V n p (y,t) through the evaluation of combinations of Chebychev polynomials Φ m at Gauss-Lobatto collocation points y k as follows

V n p (y k ,t) = M ∑ m=1 p n m (t)Φ m (y k ) with 1 ≤ k ≤ M. (4)
Finally, the null boundary conditions of the closed-loop control system is obtained by setting the upper and lower boundaries to the values of the control inputs χ u and χ l respectively. The details of the derivation of the reduced linearized model are given in [START_REF] Mckernan | Control of plane Poiseuille flow: a theoretical and computational investigation[END_REF]).

All computation done, the reduced linearized model can be presented in the following canonical expression

   ṗn (t) = A n p n (t) + B n u(t) z(t) = C ⊤n p n (t) p n (0) = p n 0 ( 5
)
where p n (t) = (p n m ) is the state vector, A n is the state matrix, u(t) = (u u (t), u l (t)) is the system control inputs on the upper and lower channel boundaries, B n is the input matrix, C n is the output matrix and z(t) is the vector of shear stress measurements on the upper and lower boundaries.

A classical problem and the existing solutions

In this section we present the unstable flow problem related to the 2D plane Poiseuille flow. Then we present existing control solutions based on shear stress measurements and we also highlight their main limitations.

Unstable flow

In the case where the Reynolds number is set to Re = 10 000, the reduced linearized model ( 5) is unstable only for the wavenumber α n = 1 [START_REF] Orszag | Accurate solution of the Orr-Sommerfeld stability equation[END_REF]). This instability can be seen through the poles of the state matrix A n . These poles are illustrated on Fig. 3, where we can see the unstable mode λ = 0.00373967 ± i0.23752649. Note that a mode is a complex conjugate pair of eigenvalues. In this case the flow is initially in the steady state but in an unstable equilibrium, i.e. a small disturbance velocity value V p (x, y,t) destabilizes the uncontrolled fluid flow.

Existing control solutions

Using the classical output feedback control u(t) = -αz(t) in ( 5) with α a scalar gain z(t) is a single shear stress measurement, the unstable 2D Poiseuille flow can be stabilized as shown in [START_REF] Joshi | A systems theory approach to the feedback stabilization on infinitesimal and finite amplitude disturbances in plane Poiseuille flow[END_REF]). However this simple proportional controller generally fails to suppress unobservable high transient modes which could trigger transition to turbulence as explained in [START_REF] Joshi | A systems theory approach to the feedback stabilization on infinitesimal and finite amplitude disturbances in plane Poiseuille flow[END_REF]. By assuming the availability of the value of the state vector p n (t), a state feedback Linear Quadradic Regulator (LQR), can easily stabilize the unstable pertubed Poseuille flow [START_REF] Bewley | Optimal and robust control and estimation of linear paths to transition[END_REF], [START_REF] Mckernan | Control of plane Poiseuille flow: a theoretical and computational investigation[END_REF]). However, in practice the LQR approach can not be used since it requires the unknown value of the true state value p n (t). That is why the output feedback LQG regulation is generally considered as an effective approach for flow control issues [START_REF] Bewley | Optimal and robust control and estimation of linear paths to transition[END_REF]Liu (1998), McKernan (2006)). This approach uses an estimated value p n (t) of the state vector: the value of p n (t) is obtained from the shear stress measurements z(t) using an observer built from a Linear Quadratic Estimation (LQE) scheme. The control signal for the output feedback LQG regulator is thus given by

u(t) = -k ⊤ p n (z(t)), (6)
where vector k is the LQR optimal gain. We recall that the LQR signal minimizes the following cost function

∞ 0 p n′ (t)Qp n (t) + u ′ (t)Ru(t) dt (7)
where Q and R are positive semidefinite weighting matrices. Matrix Q is used to limit the kinetic energy density E n (t) = p n′ (t)Qp n (t) developed in [START_REF] Mckernan | Control of plane Poiseuille flow: a theoretical and computational investigation[END_REF]). Matrix R is used to limit energy comsumption by actuators. In order to maintain wall symmetry, matrix R is set as a scaled identity matrix (McKernan ( 2006)), i.e. R = r 2 I, where r is a real positive parameter. This last control law will be refer to as shear stress based LQG (SSB-LQG) control.

The obsever used by the SSB-LQG control suffers from the initialization issue as we will show in Section 6.1. In addition in the case where the limited shear stress measurements z(t) are noisy, the estimate p n (t) is also noisy. Both of these problems are solved the proposed vision-based control approach, which is introduced in the next section.

FUNDAMENTALS OF VISION-BASED CON-TROL

To achieve a closed-loop vision-based control task, a set of visual features s(t) is selected from the image of the scene. Indeed, only a part of the image (provided by the sensor) is used to define a diffeomorphic map between the observed scene and a judiciously selected set of features in the image. A control law is then designed so that the visual features s(t) reaches a desired value s * corresponding to a desired state of the system. The control principle is thus to regulate the error vector e(t) = s(t)-s * to zero.

To design the control law, the dynamic of the error vector must be known. This dynamic is given by

ė(t) = ∂ e(t) ∂t + L e (t)u(t), (8)
where u(t) is the system control inputs, L e (t) is the jacobian matrix that encodes the time variation of the visual features with respect to the variation of the control signal acting on the system [START_REF] Espiau | A new approach to visual servoing in robotics[END_REF]), and ∂ e(t)/∂t expresses the variation of the error vector due to the free motion of the visual features.

A key point in vision-based control is that this control technique belongs to the class of sensor-based control of dynamic systems: the control law is computed in the sensor frame [START_REF] Samson | Robot control: the Task Function approach[END_REF]). Consequently, this approach corresponds clearly to an observer-free feedback control. In the following we apply the vision-based approach to the regulation of plane Poiseuille flow.

VISION-BASED CONTROL OF FLOWS

In the particular case of flows control, a control law is designed from visual features obtained from the vision system sensing the flow. Of course a great advantage of such a sensor is that it is non-intrusive. This sensor is also an extremely rich and dense source of information on the flow. Indeed a large spectrum of visual features s(t) could be selected from the image, such as coordinates of singular points in the flow or the vorticity map. Nevertheless, to directly compare our vision-based approach to the SSB-LQG approach we choose the same data, i.e. s(t) = p n (t). However, as shown in the next section, the main difference between both approaches is the way to estimate p n (t). We first show in our approach how p n (t) is estimated from visual measurements, and then we present the control law.

Visual sensing of flows

A laser sheet is used to enlighten the particules for which the velocity has to be computed (see Fig. 4). Consequently, from this visualization process, it is possible to compute dense flow velocity maps from optical flow techniques. Optical flow can be defined as the apparent velocity vector field representing the motion of photometric pattern (pixels brightness) in successive image sequences [START_REF] Heitz | Variational fluid flow measurements from image sequences: synopsis and perspectives[END_REF]). We first present the perspective projection of a flow particle, then we show how to estimate a flow particle velocity from its image velocity, and finally we present the computation of the state vector from the velocity of a fow particle.

Perspective projection of a flow particle:

Let o M = ( o M x (t), o M y (t), o M z (t)) be the spacetime coordinates of the flow particle M (see Fig. 5) expressed in the flow frame F o . The perspective projection of M is obtained in three steps.

• The first step consists to express M in the camera frame. Let c M = ( c M x (t), c M y (t), c M z (t)) be the space-time coordinates of M expressed in the camera frame F c . The relationship between c M(t) and o M(t) is given by where

c M(t) = c R o o M(t) + c t o (9)
( c R o , c t o )
is the rigid constant kinematic link between the camera and the flow frames (see Fig. 5). This rigid link is also known as extrinsic camera parameters.

• In the second step, the perspective projection c m(t) = ( c m x (t), c m y (t)) of point c M(t) obtained from the well-known pinhole camera model is given by:

c m x (t) = c M x (t) c M z (t) , c m y (t) = c M y (t) c M z (t) . ( 10
)
• Finally in the last step the perspective coordinates vector c m(t) is expressed in the sensor space (i.e. in pixel unit) as p m(t) = ( p m x (t), p m y (t)), by the relation

p m x (t) = f l x c m x (t) + u 0 , p m y (t) = f l y c m y (t) + v 0 ( 11
)
where l x (respectively l y ) is the pixel size (in meter) in the x (respectively) y direction, f is the focal length and (u 0 , v 0 ) is the vector coordinates of the principal point of the camera. Note that vector ( f x = f /l x , f x = f /l x , u 0 , v 0 ) represents the intrinsic camera parameters. Both the intrinsic and the extrinsic camera parameters can be determined using the calibration method described in [START_REF] Tsai | A versatile camera calibration technique for highaccuracy 3D machine vision metrology using off-the-shelf tv cameras and lenses[END_REF]).

Estimation of the flow velocity particle from its image velocity:

Now we show in three steps how to compute the velocity of a flow particle M from its perspective image p m. The first step consists in expressing the relationship between the flow particle and its perspective image velocities. From (11), it is easy to show that the relationship between the image velocities p ṁ and the flow veloc-

ities o Ṁ(t) = o Ṁx (t) o Ṁy (t) o Ṁz (t)
⊤ is given by

p ṁ(t) = ∂ p m ∂ c M c R o o Ṁ ( 12
)
where the jacobian matrix ∂ p m ∂ c M is given by

∂ c m ∂ c M =   f u M z (t) 0 - -f u M x (t) M 2 z (t) 0 f v M z (t) - f v M y (t) M 2 z (t)   . (13)
In the 2D case, since o M z (t) is constant, we have o Ṁz (t) = 0 and ( 12) can be rewritten as

p ṁx (t) p ṁy (t) = ∂ p m ∂ c M c R o | 12 o Ṁx (t) o Ṁy (t) (14)
where R o | 12 is the restriction of the orientation matrix c R o to its two first columns. The second step consists to estimate the image velocity p ṁ. This estimation can be done using optical flow techniques [START_REF] Heitz | Variational fluid flow measurements from image sequences: synopsis and perspectives[END_REF]). In brief, this technique exploits the brightness conservation equation at pixel p m(t) as follows

0 = dI m (t) dt = ∂ I m (t) ∂t + ∇I ⊤ m p ṁ, ( 15
)
where I m (t) is the brightness of pixel p m(t) and ∇I m is the brightness spatial gradient. There are many methods which can be used to compute p ṁ(t). All of these methods use the assumption that p ṁ is also the image velocities of all the pixels in the local neighbourhood of pixel p m(t). The most advanced methods take also into account an a priori information such as a local parametric motion model. A comprehensive review of some NSE consitent optical flow methods is available in [START_REF] Heitz | Variational fluid flow measurements from image sequences: synopsis and perspectives[END_REF]).

Finally, assuming a perfectly calibrated camera, the flow perturbation velocities V p (x, y,t) can be computed from optical flow measurements by inverting (14).

Estimation of the state vector:

This last step shows how to compute p n (t) from the estimation V p (x, y,t). This is simply done by projecting the perturbation velocities V p (x, y,t) onto Fourier and Chebychev bases. Due to a lack of space, the details of this projection are omitted in this paper.

It is worth mentioning that an important contribution of our vision-based state estimation method is the fact that the initial value p n (t = 0) is no longer of concerned in our approach as shown in Section 6.2. In addition, because of dense flow velocity maps, the vision-based approach provides less noisy estimations of the state vector as shown also in Section 6.2.

Closed-loop vision-based control of flows

As mentioned at the very begining of Section 5, we choose the same data s(t) = p n (t) as in the SSB-LQG approach so that our vision-based approach can be directly compare to the SSB-LQG approach. In that case, it becomes easy to express the dynamic of the error of the visual features e(t) (see ( 8)) around the steady state solution (the equilibrium point). In the ideal case where there is no state disturbances and there is no image noise, from the ideal state dynamic equation given in (5), it is straightforward that the ideal image error dynamic given by ( 8) is such that ∂ e(t) ∂t = A n e(t), L e (t) = B n with e(t) = p n (t).

(16)

More precisely ( 16) can be interpreted as follows: the term

∂ e(t)
∂t describes the instationary aspect of the flow image velocity map due to the motion of the uncontrolled flow, and the term L e (t) encodes the spatial inhomogeneity aspect of the flow image velocity map due to unsteady actuation (for instance blowing and suction actions). Therefore, around the desired state, (8) can be rewritten as

ṗn (t) = A n p n (t) + B n u(t) (17)
and the simple state-feedback control law

u(t) = -k ⊤ p n (t), (18)
can be used. We will refer to this vision-based control law as the vision-based LQG (VB-LQG) control law since p n (t) is obtained from visual measurements p ṁ(t) instead of shear stress measurements as used in (6). In the next section, we compare the SSB-LQG approach with the proposed VB-LQG method.

COMPARISON BETWEEN THE SSB-LQG AND THE VB-LQG APPROACHES

The major difference between the SSB-LQG and the VB-LQG approaches is the estimation of the state vector p n (t): the SSB-LQG regulator uses an observer built from the LQE approach whereas the VB-LQG approach relies on an observer-free estimation method. We now highlight the influence of these estimation methods on the closed-loop system.

Closed-loop system with the SSB-LQG control

In the LQE framework it is supposed that the reduced linearized system (5) has process disturbances ε p and measurements noise ε z which are uncorrelated white Gaussian noise. The estimation error dynamic is given by

˙ p n (t) -ṗn (t) = (A n -LC n⊤ )( p n (t) -p n (t)) + Lε z (t) -ε p (t) p n (0) -p n (0) = unknown, (19)
where L is the estimator optimal gain. The resulting closedloop system is given by [START_REF] Mckernan | Control of plane Poiseuille flow: a theoretical and computational investigation[END_REF]

):      ṗn (t) ˙ p n (t) -ṗn (t) = M lqg p n (t) p n (t) -p n (t) + M n ε p (t) ε z (t) p n (0) -p n (0) = unknown. (20)
where

M lqg = A n -B n k ⊤ -B n k ⊤ 0 A n -LC ⊤n , M n =
I 0 -I L with I an identity matrix. Equation ( 20) clearly shows that the true state dynamic ṗn (t) (around the desired state) depends on the estimation error p n (t)p n (t). Therefore, since this error highly depends on the initial unknown estimation error p n (0)p n (0), a poor initialization of the observer could cause drive the system to a turbulent state as shown in [START_REF] Mckernan | Control of plane Poiseuille flow: a theoretical and computational investigation[END_REF]. In addition the noise in the measurements propagates in the control law (6), this is not suitable at all for the lifetime of the actuators.

Closed-loop system with the VB-LQG control

We first present the estimation of the state vector p n (t) in the case of noisy optical flow measurements and then we show system closed by the VB-LQG control. Let V p (x, y,t) be the M × N pixels size image of the pertubation velocity map obtained from noise-free optical flow measurements. From ( 14), in the practical case where optical flow measurements p ṁ are corrupted by a Gaussian noise process, the computed flow perturbation velocities V p are also affected by Gaussian noise. Without lost of generality the noisy perturbation velocity map can be written as V p (x, y,t) = V p (x, y,t) + ε(x, y,t). This noisy pertubation velocity map V p (x, y,t) is used to compute the estimation of the state vector p n (t) as shown in [START_REF] Fomena | Fluid flows control using visual servoing[END_REF]):

p n (t) = p n (t) + 1 N e n (t), (21)
where vector e n (t) the projection of the measurements noise matrix ε(x, y,t) onto Fourier and Chebychev basis. The resulting vision-based closed-loop system, obtained from ( 5) and ( 21), is given by

   ṗn (t) = A n p n (t) + B n u(t) u(t) = -k ⊤ p n (t) p n (t) = p n (t) + 1 N e n (t), ( 22
)
where k is the LQR optimal gain. The closed-loop reduced linearized system ( 22) can be rewritten as

ṗn (t) = (A n -B n k ⊤ )p n (t) - 1 N B n k ⊤ e n (t). (23)
The initial value p n (t = 0) is therefore no longer of concerned in our approach. In addition for a large number of points N in the streamwise direction (typically N=1024 for a standard image), the reduced linearized system dynamic equation ( 23) is less affected by measurements noise since 1 N B n k ⊤ e n (t) tends to 0. This is another great improvement over the SSB-LQG control scheme that is always noise dependent when noisy shear stress values are used in the LQE as shown in (20).

RESULTS

Validations of our approach is done using synthetic data sets of spatio-temporal variations of the perturbation velocities V p (x, y,t) obtained from the Poiseuille flow reduced model presented in (5). This reduced model is used to build a linear simulator of 2D channel Poiseuille flow with the following classical characteristics as used in related works [START_REF] Joshi | A systems theory approach to the feedback stabilization on infinitesimal and finite amplitude disturbances in plane Poiseuille flow[END_REF], [START_REF] Bewley | Optimal and robust control and estimation of linear paths to transition[END_REF]Liu (1998), McKernan (2006)): R e = 10 000, L x = 4π and α n = 1. Matrices A, B, C given in (5) have been computed for M = 250 using Matlab codes provided in [START_REF] Mckernan | Control of plane Poiseuille flow: a theoretical and computational investigation[END_REF]).

Comparison of the state estimations

Results are given in Fig. 6. Fig. 6(a) presents the ideal case where we can see that both estimations perfectly correspond to the ground truth value of the state vector. Fig. 6(b) highlights the initialization and the asymptotic convergence issues in the LQE; these issues are not of concerned in the vision-based approach which provides the ground truth value of the state vector. In order to show the robustness to noise of the vision-based approach, the standard deviation (STD) on the optical flow noise has been purposely set to a value 10 times higher than the STD on the shear stress noise. As shown on Fig. 6(c), due to a large number of flow velocities provided by the optical flow, the new approach is robust to noisy measurements.

Comparison of the behavior of the closedloop systems

Results are presented in Fig. 7. Fig. 7(a) depicts the behavior of the control signal in the ideal case (no measurements noise, no initialization error). Fig. 7(b) depicts the behavior of the control signal when the initial value is set as p n (0) = 0 by default since the value of p n (0) is unknown. In this case we can see that the value of the control signal can be 100 times higher than the ideal control signal case which includes the LQG-VB approach (compare the highest control signal values in Fig. 7(b) and 7(a)). This higher control signal value could lead to an unsuitable state trajectory which can cause the real non-linear system to diverge as shown in [START_REF] Mckernan | Control of plane Poiseuille flow: a theoretical and computational investigation[END_REF]). In addition, as expected, the control signal (see Fig. 7(b)) takes more time to converge to 0 (3000 iterations compared to the VB-LQG approach). This leads to an energy consumption far much higher for the SSB-LQG control than for the VB-LQG control. The figures on the second row present the measurements noise case. Fig. 7(c) pictures the case of the SSB-LQG control, where we can see that the control signal does not converge to 0. Although the noise STD has been set to a small value, σ ss = 0.03, the control signal is very noisy, which is not suitable for actuators. Finally,Fig. 7(d) illustrates the robustness of the vision-based control where the STD in the optical flow noise is 10 times higher than the STD in the shear stress noise: we can see from this last figure that the control converges, and we can also see that the larger the sample of flow velocities used the lesser the noise in the control signal.

CONCLUSION

In this paper we have proposed a vision-based approach for fluid flow estimation and control which builds on optical flow techniques, which is non-intrusive, robust to measurements noise and suitable for fluid flows control. Theoretical proofs have been given in order to validate the improvements provided by the vision-based approach over the shear stress based LQG control in terms of the state vector estimation and flows control. Simulation results have been presented to validate the proposed approach. Future works will focus on the extension of the vision-based approach to the control of 3D flows.

Fig. 1 Figure 1 .

 11 Fig. 1 illustrates the steady state velocities profile of the 2D plane Poiseuille flow where we can see that the x-direction is associated to the streamwise direction and the y-direction is associated to the normal direction. A 2D plane Poiseuille flow is a flow in a infinite length channel due to a pressure gradient. The non-dimensionalised Navier-Stokes equations (NSE) of this flow are given by

Figure 2 .

 2 Figure 2. Boundaries control of a 2D plane Poiseuille flow, L x is the streamwise channel length.

Figure 3 .

 3 Figure 3. Poles and zeros of the reduced linearized system for Re = 10 000, α n = 1.

Figure 4 .

 4 Figure 4. Fronto-parallel visualization of a 2D flow using a laser sheet which role is to enlighten the particles seeded in the fluid: p m is the perspective image of the fluid particle M .

Figure 5 .

 5 Figure 5. Perspective projection of a flow particle.

Figure 6 .Figure 7 .

 67 Figure 6. Comparison of state vector estimations using shear stress and optical flow: (a) ideal case, (b) LQE poor initialization, (c) measurements noise with a large number of pixels N.