El Kamal

Yves Omari

Le Guer

Kamal El Omari
email: kamal.elomari@univ-pau.fr

Yves Le Guer
email: yves.leguer@univ-pau.fr

Alternate rotating walls for thermal chaotic mixing

Keywords: Chaotic mixing, Laminar mixing, Heat transfer, High Prandtl number, Thermal eigenmodes, Unstructured finite volume method

published or not. The documents may come

Nomenclature

A c surface of a mesh element (m 2 ) A T σ composite mixing indicator c p heat capacity (J.kg -1 .K -1 ) d vector joining the centers of two cells E energy (J) k thermal conductivity (W.m -1 .K

Introduction

The topic of heat transfer and mixing plays an important role in many fields as various as earth and life sciences, chemical engineering and material science. This subject is particularly difficult to solve when highly viscous fluids are concerned. A variety of industrial processes involve the heating of highly viscous liquids; they include material processing (molten polymers for example), petroleum engineering, food engineering or environmental engineering. The motivation for this study is to analyze the heating and mixing performance of a chaotic mixer with a rather simple design. Chaotic mixing is obtained by considering the chaotic advection phenomenon [START_REF] Aref | Stirring by chaotic advection[END_REF][START_REF] Aref | The development of chaotic advection[END_REF][START_REF] Ottino | The kinematics of mixing: stretching, chaos, and transport[END_REF] that occurs in laminar flow. Today it is well recognized that chaotic mixing is the most efficient mixing process when highly viscous fluids are involved [START_REF] Raynal | Energy saving in chaotic laminar mixing[END_REF]. Chaotic mixing is also recommended for the mixing of delicate fluids that do not resist the high strains encountered in turbulent flows. Despite the simplicity of their velocity fields, chaotic advection flows are able to create very complex patterns of the advected scalar with highly stretched and folded structures. The more numerous are very thin striations produced by the flow, the more efficient the diffusion and the faster the homogenization of the scalar.

When heat transfer is considered, temperature is the scalar and heat conduction is the diffusion source. In the present study we focus only on the relationship between chaotic mixing and the enhancement of heat transfer in a particular flow. Despite its obvious industrial relevance, there are only a limited number of works that consider this problem. In fact, many more works in the literature concern the relationship between chaotic mixing and chemical reaction more than the coupling of chaotic mixing with heat transfer. When heat transfer enhancement is involved, mainly two classes of flow geometries are encountered in the production of chaotic mixing: those that use rotating elements as eccentric cylinders [START_REF] Ghosh | Heat-transfer enhancement due to slender recirculation and chaotic transport between counter-rotating eccentric cylinders[END_REF][START_REF] Ganesan | Chaotic heat transfer enhancement in rotating eccentric annular-flow systems[END_REF][START_REF] Lefevre | Chaotic advection and heat transfer enhancement in Stokes flows[END_REF][START_REF] Mota | Optimization of heat-transfer rate into time-periodic two-dimensional Stokes flows[END_REF] and those that use multiple pipe bends [START_REF] Peerhossaini | Heat exchanger design based on chaotic advection[END_REF][START_REF] Mokrani | The effects of chaotic advection on heat transfer[END_REF][START_REF] Chagny | Chaotic heat transfer for heat exchanger design and comparison with a regular regime for a large range of Reynolds numbers[END_REF][START_REF] Acharya | Analysis of heat transfer enhancement in coiled-tube heat exchangers[END_REF][START_REF] Lemenand | A thermal model for prediction of the Nusselt number in a pipe with chaotic flow[END_REF][START_REF] Guer | A mapping tool using anisotropic unstructured meshes to study mixing in periodic flows[END_REF][START_REF] Boesinger | Experimental study of reactive chaotic flows in tubular reactors[END_REF][START_REF] Gibout | Coupling of a mapping method and a genetic algorithm to optimize mixing efficiency in periodic chaotic flows[END_REF][START_REF] Kumar | Coiled flow inverter as a heat exchanger[END_REF][START_REF] Kumar | Laminar convective heat transfer in chaotic configuration[END_REF][START_REF] Yamagishi | Chaotic analysis of mixing enhancement in steady laminar flows through multiple pipe bends[END_REF]. For the latter geometries, chaotic trajectories are determined by continuous modification of pipe wall design; the third space dimension provides the required additional degree of freedom and no displacement of the wall is imposed. As a consequence, the number of operating parameters that determines the options for effective stirring protocols is reduced, and the flow structure is given only by the chosen pressure gradient (i.e., the Reynolds number for a particular fluid). Another interesting 3D geometry is the continuous chaotic flow generated inside the Rotated Arc Mixer of Lester et al. [START_REF] Lester | Low Reynolds number scalar transport enhancement in viscous and non-Newtonian fluids[END_REF]. The mixer consists of an inner stationary cylinder lacking internal structure, with apertures cut through its wall, wrapped by a concentric outer rotating cylinder. This geometry offers a large, tunable, optimization parameter space for the enhancement of heat and/or mass transfer. For all of the cases discussed above, thermal chaotic mixing is used to attain multiple objectives: to encourage heat transfer between the heated walls and the neighboring fluid and to achieve temperature homogeneity for the whole fluid domain. For this last condition, we need to both prevent the formation of a large, unmixed KAM island, where the fluid remains cold, and also avoid hot spots in the vicinity of heated walls, which can degrade the fluid.

One of the main differences between thermal chaotic mixing and reactive chaotic mixing is due to the large contrast in transfer rates between mass and heat. Effectively, it is common to notice a difference of two decades or more between the values of the molecular and thermal diffusivities. Thermal diffusivities are higher than molecular diffusivities, which has consequences on the solutions obtained for the advection-diffusion problem. Typically, the patterns of the strange eigenmodes [START_REF] Pierrehumbert | Tracer microstructure in the large-eddy dominated regime[END_REF][START_REF] Liu | Strange eigenmodes and decay of variance in the mixing of diffusive tracers[END_REF] obtained for the scalar fields will differ. These patterns appear repeatedly for periodic velocity fields and repeat themselves every period with an exponential decay of the scalar contrast. More often, asymptotic transport is controlled by the slowest eigenmodes, but that is not always the case [START_REF] Pikovsky | Persistent patterns in deterministic mixing flows[END_REF][START_REF] Popovych | Abnormal mixing of passive scalars in chaotic flows[END_REF].

Thus, the main issue of the present work is to analyze the heating and mixing performances in a realistic mixer of simple geometry and, in particular, to study the effect of wall rotation on the enhancement of heat transport in the whole fluid domain. This crucial effect of wall on the mixing efficiency has been recently highlighted by Gouillart et al. [START_REF] Gouillart | Walls Inhibit Chaotic Mixing[END_REF][START_REF] Gouillart | Slow decay of concentration variance due to no-slip walls in chaotic mixing[END_REF] for the homogenization of concentration in a 2D closed flow environment.

This article is organized as follows: in section 2, we present the problem, the geometrical description of the mixer, the fluid properties and the flow parameters. We also describe the different types of stirring protocols used to study the performance of the mixer for heat transport. The section ends with a presentation of the governing equations. In section 3, the numerical methods used to perform the simulations are described; in section 4, the mixing and energy indicators defined to quantify the efficiency of heat transport and homogenization are presented. All the results are detailed and discussed in section 5. At first, different flow topologies resulting from the choice of the stirring protocol imposed on the walls are shown, and their incidences on mixing and heat transport are given. The results are also analyzed with the help of temperature probability distribution functions (PDFs), which are very useful in globally characterizing the homogenization of temperature in time. We show the existence of recurrent patterns of the temperature fields and explain them in terms of thermal strange eigenmodes. We also analyze the temporal evolution of the temperature gradients in the fluid domain. We conclude the section with a study of the effect of rod eccentricity. Finally, in section 6, concluding remarks are drawn and some perspectives for future work are given.

Problem statement

Geometrical description: the two-rod rotating mixer

A sketch of the studied mixer is presented in Figure 1. It is composed of two circular rods of equal radii, which are maintained vertically inside a cylindrical tank (a bounded domain). The tank and the rods are heated and can rotate around their respective revolution axis. This geometry is similar to the two-roll-mills studied in the literature by Price et al. [START_REF] Price | Numerical and experimental characterization of a family of two-roll-mill flows[END_REF] ; however, in our case, the outer cylinder has the ability to rotate. The rotation of the outer cylinder induces important consequences for the flow and mixing properties that will be explained later.

The geometry of the flow domain is characterized by two dimensionless parameters: the ratio of the radii of the rod and the cylindrical tank (R 1 or R 2 )/R 3 , and the dimensionless eccentricity e = ε/R 3 . The geometrical parameters are set to (R 1 or R 2 )/R 3 = 1/5 and a range between 0.35 and 0.60 for e (R 1 = R 2 = 10 mm , R 3 = 50 mm and an eccentricity of ε = 25 mm for the standard case study). In section 5.4, we will investigate the influence of variable eccentricity on thermal mixing.

The rods and the tank rotate alternately or with a continuous modulation of velocities. They can also have different directions of rotation. The classical eccentric cylinder geometry studied in the past [START_REF] Ghosh | Heat-transfer enhancement due to slender recirculation and chaotic transport between counter-rotating eccentric cylinders[END_REF][START_REF] Ganesan | Chaotic heat transfer enhancement in rotating eccentric annular-flow systems[END_REF][START_REF] Lefevre | Chaotic advection and heat transfer enhancement in Stokes flows[END_REF][START_REF] Mota | Optimization of heat-transfer rate into time-periodic two-dimensional Stokes flows[END_REF] produced a large recirculation bubble, which could not be avoided even with a careful choice of stirring protocols. In contrast, the two-rod mixer [START_REF] Omari | A numerical study of thermal chaotic mixing in a two rod rotating mixer[END_REF] studied here is suitable to obtain full chaotic flow without KAM regions.

Fluid properties and flow parameters

The fluid considered in this study has the thermophysical properties listed in Table 1 This fluid is Newtonian and has a high Prandtl number of P r = 10 4 . Considering these properties, the thermal diffusivity a is equal to 1.515 × 10 -7 m 2 .s -1 . The maximum angular wall velocities are fixed to: Ω 1 = Ω 2 = 30 rpm for the rods and Ω 3 = 6 rpm for the outer tank. The tangential wall velocity is then the same, and equal to U = 31.41 mm.s -1 . A characteristic Reynolds number for this flow can be evaluated as:

Re = ρ U • 2 • (R 3 -R 1 ) µ = 1.66 (1)
The walls of the rods and tank are kept at a constant hot temperature of T hot . The initial uniform temperature of the fluid (before heating) is set to a cold temperature T cold . We define the dimensionless fluid temperature as:

T * = T -T cold T hot -T cold (2)
Therefore, T * cold = 0 and T * hot = 1 and the maximum temperature difference between the walls and the fluid is always 1. The heated walls play the role of a variable heat source that is continuously dissipated in the 2D modulated flow field. This is a specific situation, not usually encountered when dealing with the scalar dissipation of concentration in a 2D velocity field. The Dirichlet boundary condition imposed to the walls implies that during the mixing of the fluid in the tank, the parietal fluxes will change along the walls as well as with time depending on the local flow conditions.

The other important dimensionless parameter for this non-isothermal mixing problem is the Péclet number, Pe = Re Pr. The Péclet number can be seen as the ratio of thermal diffusion time τ td over advection time τ ad for the scalar temperature. The limit Pe = 0 corresponds to the pure diffusion case. In our case, Pe is large (1.66 10 4 and then τ td is 16, 600 times greater τ ad ); hence, the need to speed-up the mixing of temperature scalars is clearly demonstrated.

Stirring protocols

Chaotic mixing is produced here by varying in time the angular speed of the circular rods and the tank using a sine-squared waveform. The stirring protocols studied can be defined by two parameters: the respective direction of rotation between the rods and cylindrical tank, and the duration of the time periodic modulation. Other parameters, such as the ratio between the maximum velocities of the rod and cylindrical tank, and the amplitude of the angular velocity modulation, are not considered in this work. The three possible stirring configurations corresponding to the different flow topologies are specified in Table 2. The symbol (+) indicates a counter-clockwise direction and (-), a clockwise direction. For chaotic flows, two types of temporal modulation are considered for the stirring protocols: a continuous one (a sine-squared modulation of wall velocities) and a non-continuous one (the rods are stopped together for half a period, during which time the outer cylinder is rotating; for the next half-period the contrary is considered). These stirring protocols are illustrated in Figure 2. The maximum angular velocity is the same for the two types of modulation. The rods and the cylinder tank follow sine-squared modulations as defined by equations 3 (here for a modulation period τ = 30 s): Non-continuous case:

Stirring configuration Rod 1 Rod 2 Tank 1 (+) (-) (+) 2 (+) (+) (+) 3 (-) (-) (+)
Ω 1,2 = 30 -60 sin 2 ( πt T ) if Ω 1,2 < 0 ⇒ Ω 1,2 = 0 Ω 3 = 6 -12 sin 2 ( πt T + π 2 ) if Ω 3 < 0 ⇒ Ω 3 = 0 Continuous case: Ω 1,2 = 30 -15 sin 2 ( πt T ) Ω 3 = 6 -3 sin 2 ( πt T + π 2 ) (3)
It was common practice in the chaotic mixing process design to introduce a temporal modulation to the velocity of movement of mixer walls, which generally enhanced the efficiency of mixing. But often, the alternate rotation was not investigated.

In this work, the modulation period was varied from 15 to 100 s. For a modulation period of 15 s, the Strouhal number is typically on the order of St = 0.34. The Strouhal number is defined as:

St = 2 • (R 3 -R 1 ) U moy • τ (4)
with U moy ≈ 0.5 U . Considering this value, we can notice that this flow does not satisfy the quasi-steady hypothesis. As a consequence, in our modeling of the flow we need to account for inertial effects.

Governing Equations

The evolution of the scalar temperature transported in an incompressible laminar flow of a Newtonian fluid is described by the advection-diffusion equation written in its integral form:

∂ ∂t V ρc p T dV + S ρc p T U • n dS = S k ∇T • n dS ( 5
)
where k is the thermal conductivity. The integration is over a volume V surrounded by a surface S, which is oriented by the outward unit normal vector n. As indicated before, the fluid is initially at a uniform cold temperature T cold , and during the process the walls are maintained at a uniform hot temperature T hot (Dirichlet condition). The velocity field of the flow U is obtained by resolution of the unsteady Navier-Stokes equations:

S U • n dS = 0 (6) ∂ ∂t V ρ U dV + S ρ U U • n dS = V -∇p dV + S = τ • n dS (7)
where = τ is the viscous stress tensor. Here, the full Navier-Stokes equations are solved in order to take into account the inertial effects that are induced by the accelerated or decelerated movement imposed on the walls. The evolution of the fluid velocity and temperature fields are obtained by a CFD method, described in the next section.

Numerical methods

Two-dimensional continuity and Navier-Stokes equations, as well as the energy conservation equation, are solved by means of a inhouse code called Tamaris. This code has an unstructured finite volume framework that is applied to hybrid meshes. Variable values ( U , p and T) are stored at cell centers in a collocated arrangement. Cell shapes can be triangular or quadrilateral. To describe the discretization practice used in the code, we can write the above equations (5 and 7) in a unique generic form:

∂ ∂t V ρφ dV + S ρφ U • n dS = S Γ ∇φ • n dS + V S φ dV (8)
where S φ is a source term. Spatial schemes approximating diffusive and convective fluxes are both second-order accurate. The discretization of the diffusion term is performed by approximating the surface integrals by a sum over all the cell faces f (see Figure 3):

S Γ ∇φ • n dS = f Γ f A f ( -→ ∇φ) f • n f (9)
In unstructured meshes, the normal gradient ( -→ ∇φ) f • n f is decomposed into an implicit contribution that uses the values of f at the centers of the two cells sharing the face f (first term in the RHS of Eq. 10 ) and a non-orthogonality correction term treated explicitly by a deferred approach in order to preserve second order accuracy:

( -→ ∇φ) f • n f = φ nb -φ c || d|| + -→ ∇φ • n - d || d|| (10)
d is the vector joining the centers of the two cells (see Figure 3). The average gradient -→ ∇φ is interpolated from the gradients of these neighboring cells. The gradients of variables at the cell centers are computed by a second order least-square method that uses weighted values of φ at all the neighboring cells. Convection terms are also transformed in a sum over the faces f composing the surface S:

S ρφ U • n dS = f (ρφA) f U f • n f (11)
where face values φ f require an appropriate interpolation. In the momentum equations, the blending difference scheme is used [START_REF] Ferziger | Computational methods for fluid dynamics[END_REF]. This scheme combines a second order accurate centered scheme with a low amount of upwind differencing to enforce the stability of the global scheme:

φ f = φ upwind f + β(φ centered f -φ upwind f ) m-1 (12)
where m -1 stands for the previous iteration. Here the weighting factor β is fixed at 0.95.

In the energy equation, special care was given to the face value interpolation in order to limit the numerical diffusion introduced, while preserving a strict boundedness. For this purpose the Normalized Variable Diagram (NVD) Gamma high resolution scheme of Jasak et al. [START_REF] Jasak | High resolution NVD differencing scheme for arbitrarily unstructured meshes[END_REF] is used. This scheme introduces a smooth transition between centered and upwind differencing, depending on the local variation of the variable. This property is important for our study since we need to evaluate the efficiency of the heat advection mechanism in the absence of the numerical diffusion effect that may overshoot the thermal conduction role. Pressurevelocity coupling is ensured by the SIMPLE algorithm [START_REF] Patankar | Numerical heat transfer and fluid flow Series in computational methods in mechanics and thermal sciences[END_REF], while the mass fluxes at the cell faces are evaluated using the Rhie-Chow interpolation [START_REF] Rhie | Numerical study of the turbulent flow past an airfoil with trailing edge separation[END_REF] to avoid pressure checkerboarding. Unsteady computations are advanced in time by the implicit three time step Gear's scheme that ensures second-order accuracy. At each iteration, the discretization technique presented above leads to the creation of a linear system with a non-symmetric sparse matrix for each variable. These linear systems are solved by means of an ILU preconditioned GMRES solver using the implementation of the IML++ library [START_REF] Dongarra | A sparse matrix library in C++ for high performance architectures[END_REF]. After a grid size-dependence study, a mesh of 9944 computational cells is adopted. This mesh was generated by the open-source software Gmsh [START_REF] Geuzaine | Gmsh: a three-dimensional finite element mesh generator with built-in pre-and post-processing facilities[END_REF] and is shown in Figure 4, where the use of regular quadrilateral cells near the walls to enhance the resolution of the boundary layers is noticeable. More details about the sensitivity of the results to the grid size as well as code validation are given in our recent paper [START_REF] Omari | A numerical study of thermal chaotic mixing in a two rod rotating mixer[END_REF].

Mixing and energy indicators

As mentioned by Finn et al. [START_REF] Finn | Mixing measures for a two-dimensional chaotic Stokes flow[END_REF], a large variety of computational mixing measures exist. Many different criteria by which the quality of the mixing (with heat transfer or chemical reaction) may be evaluated will depend upon the application of interest. In order to characterize the efficiency of heat transfer by chaotic mixing, we introduce different mixing indicators, which are monitored throughout the evolution of the process. The first one is the mean temperature T * m of the fluid, which represents the energy extracted from the walls:

T * m = 1 c A c c A c T * c ( 13
)
where the summation is made over all 2D computational cells c of area A c . Indeed, this mean temperature evolution can be seen as an indicator of the ratio of the total energy supplied to the fluid from time 0 to t:

E(t) = V tank ρc p (T (t) -T cold )dV = ρc p V tank T (t)dV -T cold V tank = ρc p V tank (T m (t) -T cold ) (14)
As a consequence, the mean temperature is asymptotically bounded by the fixed temperature imposed on the walls (i.e., T hot or 1 for non-dimensional temperature). The second mixing efficiency indicator is the standard deviation σ of the fluid temperature, which accounts for the level of homogenization of the scalar temperature inside the 2D tank:

σ = 1 c A c c A c (T * c -T * m ) 2 1 2 (15)
Efficient thermal mixing requires good values for both of the aforementioned indicators (i.e., T * near 1 and σ near 0). In order to compare several stirring protocols on the basis of these two indicators, it is useful to introduce a new indicator that combines the two effects:

A T σ = 1 t f inal t f inal 0 T * m σ dt (16)
Then, for an efficient stirring protocol, A T σ must tend toward a high value. While the T * m and σ are unsteady indicators that evolve during the mixing process, A T σ is a single value that characterizes the whole process or at least during a period of t f inal . This duration is chosen in such a way that, for all the studied stirring processes, T * m exceeds 0.95.

Results and discussion

Flow topologies and the effect of stirring protocols

In order to analyze the flow topologies involved in this mixer, we first present the flows obtained without modulation of the magnitude of the angular velocity. For the three stirring configurations discussed (cf. Table 2), we plotted the steady-state streamline topologies with and without the rotation of the outer cylinder.

From the plots of Figure 5 (a), (b) and (c) showing the steady streamline flow topologies, we can observe respectively:

• (a) Two stagnation points inside the fluid located towards the outer cylindrical wall on the vertical axis, seven separatrices, four parabolic points on the wall of the outer cylinder, counter-rotating closed streamlines around the rods, and four small vortices.

• (b) One stagnation hyperbolic point located at the centre of the tank, four separatrices, four parabolic points on the wall of the outer cylinder, co-rotating closed streamlines around the rods, and two large vortices. This configuration corresponds to the one shown in the paper of Price et al. [START_REF] Price | Numerical and experimental characterization of a family of two-roll-mill flows[END_REF] concerning a study of a family of two-roll-mill flows.

• (c) Closed streamlines along the wall of the cylindrical tank as in case (b), a hyperbolic stagnation point at the center of the tank, two large vortices, two separatrices around the vortices, two parabolic points on the wall of each rod and located on the vertical axis of symmetry of the rod, and two separatrices that join the wall rods.

In the lower half of Figure 5, the outer cylinder rotates and closed streamlines are always present in a layer near the wall regardless of the stirring configuration. This layer is thicker for stirring configuration 2 (Fig. 5 (e)), corresponding to co-rotation of the rods and cylinder.

When the rods and outer cylinder rotate steadily at the same time (Fig. 5 (d), (e) and (f)), the stagnation points observed on the walls in cases (a), (b) and (c) of the same figure disappear. We know from dynamical system theory that chaotic mixing will appear from the break-up of the separatrix, which dynamically separates regions of confined closed streamlines. A separatrix defines a heteroclinic trajectory that relies on two fixed (stagnation) points in the dynamical system terminology. Through the introduction of time-periodic perturbations, this heteroclinic connection is broken into two manifolds (one stable and the other unstable) that intersect transversely and give rise to what is called a heteroclinic tangle. The breakup of the separatrix occurs in our case when the rods and tank alternate their rotation. As a consequence, each stagnation point located on a wall is interesting because it defines the starting point for separatrix emanation. We will see later that the presence of these parietal fixed points play a key role in the extraction of heat from the walls. The fixed wall temperature imposed on all of the solid boundaries acts as a source for the evolution of the scalar; thus, the mean temperature will not be constant in time but will evolve asymptotically from T * m = 0 to T * m = 1 (the whole fluid is at the wall temperature). This behavior is shown in Figure 6 for stirring configuration 3 and three wall velocity modulation conditions: nonmodulated, continuous and non-continuous (see Figure 2). In comparison, if we had considered the evolution of a material of dye for the same mixing process, the mean value of the dye concentration field would have been constant. From the temporal evolution of the standard deviation and mean temperature shown in Figure 6, it is evident that a non-continuous modulation of the velocity field enhances the mixing efficiency. After a transition time, the standard deviation of the temperature field decays exponentially for each modulation case; however, the decay rate is much higher for the non-continuous modulation condition (about a decade of difference on σ at 1500 s). While the temperature contrast decays exponentially, the mean temperature increases asymptotically towards the wall temperature T * = 1. This behavior is easily understandable when we observe the temperature fields shown in Figure 7 for the three cases of angular velocity wall modulation. After several periods (t = 120 s), we find that only the non-continuous modulation condition avoids the persistence of large, non-heated zones. The temperature fields of non-modulated and continuously modulated cases are quite similar and follow very well the closed streamlines observed in Figure 5 (f). Moreover, these closed streamlines define a barrier for the evolution of the scalar towards the interior of the cylinder tank. Indeed, if the temperature difference is confined in a relatively large boundary along each wall, both the temperature gradient normal to the wall and the heat flux will decrease. Thus, the efficiency of the whole heat exchange will be reduced. In the case of non-continuous modulation, four parabolic points appear along the wall when the tank is stopped (see Figure 5(b)). It is from these parabolic points that strips of hot fluid will be pulled out from the wall. In Figure 7(c), it is evident that two of these points are diametrically opposed on the wall of the tank. At this time (t = 120 s), the tank has been at rest since 15 s. This phenomenon of heat extraction from the walls for non-continuous stirring protocols is also clearly evident from Figure 8, which present simultaneous portraits of streamlines asso- ciated with the corresponding temperature field. The hot heat fluid strips extracted from the tank (at 120 and 132 s) and rod walls (at 140 s) towards the interior of the mixer are clearly identifiable. Streamline portraits taken at different instants of the same period of modulation show crossings when superimposed. These conditions are generic for chaotic flow. As a consequence, the temperature field exhibits varied complex patterns. Recently Hertzsch et al. [START_REF] Hertzsch | DNA Microarrays: design principles for maximizing ergodic, chaotic mixing[END_REF] have shown for a pulse source-sink chaotic system that the need to have a great transversality in streamline crossing is not always an optimal condition to improve global mixing. Hence, more theoretical studies are required to clearly understand the complex mechanism of the scalar gradient orientation with the strain axes in response to a disturbance of the velocity field imposed by the walls. This hovewer, is beyond the scope of the present work. In the following sections, only the condition of non-continuous modulation is investigated.

Effect of wall rotation direction and period length

In Figure 9, the three stirring configurations are compared over time for a modulation period of 15 s. We observe that the standard deviation of temperature decreases exponentially in each case. The relatively high standard deviation for stirring configuration 3 indicates the formation of large KAM islands in the mixer; the mixing is nevertheless good outside of these KAM regions.

In order to quantify the role of the modulation period length on mixing for the three stirring configurations, we plotted the composite indicator A T σ defined in section 4 for four modulation periods 15, 30, 60 and 100 s. We can observe in Figure 10 that this indicator reveals differences between stirring configurations mainly for low periods of modulation. An optimum seems to be a modulation period of 30 s for stirring configurations 1 and 2 and 60 s for stirring configuration 3. The latter combination is best even though configuration 3 gives, by far, the worst result at τ = 15 s. Thus, no clear preferential configuration can be deduced, but we can conclude that:

• short periods of modulation (about 15 s) should be avoided to achieve thermal mixing; 2. The modulation period is τ = 15 s.

• configuration 1 does not give the best performance regardless of the period length;

• for a period length τ = 100 s, none of the configurations give optimum results.

Some explanations of these observations can be proposed by the examination of the evolution of the temperature field as the onepresented in Figure 8. For configuration 3 with a low period, the extracted hot fluid strips from the tank boundary during the rotation of the rods are not sufficiently long; therefore, these strips are brought back to the tank when the tank begins to rotate in the opposite direction. On the other hand, a period that is generally too long reduces the frequency of flow pattern alternation as well as that of stretching and folding operations.

Temperature Probability Distribution Functions (PDFs)

We focus here on the probability distribution functions of the dimensionless temperature T * for the whole section filled by the fluid. Examples of such T * fields are shown in Figures 7 and8.

Figure 11 displays the PDFs at two different times for the three stirring configurations given in Table 2. For stirring configurations 1 and 2, which exhibit good mixing performances, a peak located is observed near the lowest temperatures. This peak moves in time to the right while its height increases. The shift of the PDF towards the right illustrates that the cold temperature disappears during the mixing process (and does so completely for good stirring configurations (configurations 1 and 2)). For stirring configuration 3, the PDF of the dimensionless temperature field remains quite flat even if the range of temperature is reduced. This result confirms the one observe in Figure 7 for the same stirring configuration. An important feature of this thermal chaotic flow can be shown in Figure 12 by the PDFs of the rescaled dimensionless temperature:

X = T * -T * m σ (17)
The PDFs, which are plotted for different times during the mixing process but at the same phase of the period, are superimposed; this is the signature of a strange eigenmode, which is characterized by the production of persistent patterns in the flow. These patterns arise from a combination of stretching, folding and thermal diffusion. The mechanism of their apparition and persistence will be detailed in the next section.

Thermal strange eigenmodes

The mechanism we investigate here is the spreading of temperature strips emerging from the walls in the bulk fluid. The flow has been chosen to produce chaotic mixing without the appearance of unmixed island (stirring protocol 2 and non-continuous modulation). The width s of a hot strip is imposed by an equilibrium between the stretching and folding and thermal diffusion. This width defines a length scale which is given by [START_REF] Gouillart | Slow decay of concentration variance due to no-slip walls in chaotic mixing[END_REF]:

s = α λ 1/2 ( 18
)
where λ is a typical rate at which the equilibrium operates and represents a stretching rate. In our flow the stretching is not constant over the whole fluid domain, thus the widths of the hot strips evolve following local stretching rates. The effect of thermal diffusion can be understood by taking into account the average characteristic size s of a computational cell (near 0.7 mm) and the thermal diffusivity coefficient (see section 2. scalar strip that falls below this size will cease to be visible after a time 1 λ = s 2 α ≈ 3 s that corresponds, for example, to a tenth of a period if τ = 30 s. Thus, the diffusive boundaries of the tracer strip interpenetrate rapidly to avoid the creation of very thin striations. The spatial structure of the temperature field is then smooth, which is due to the relatively high value of the thermal diffusivity and is also the reason why only relatively large temperature striations are observed. By comparison, the spatial structures of the concentration patterns generally present more lamellar structures due to a lowest value of mass diffusivity. This also is responsible for the nature of the pattern of the thermal strange eigenmode shown in the snapshots of Figure 13. The pattern is the same at each periodic time but the amplitude of the dimensionless temperature tends towards 1; this is also confirmed on the vertical cut-lines shown in Figure 14, where the temperature gradients correspondingly decrease. These curves were obtained by plotting the values of T * along the vertical symmetry axis of the mixer. In Figure 15, the temperature values for three of the snapshots of Figure 13 are rescaled between the T * min and T * max of each snapshot. With this scaling we observe that the patterns of the thermal eigenmode are exactly the same, which was not clearly obvious in Figure 13. Normalized vertical cut-lines of dimensionless temperature X (Eq. 17) are shown in Figure 16 for the same cases as those of

Effect of eccentricity

Figure 17 shows the difference of efficiency when changing the eccentricity. The largest eccentricity configuration (ε = 30 mm) gives the highest most probable value for the mean temperature but at the same time presents a broad tail at the left of the peak that extends towards the lower temperatures. This left tail is due to the persistence of non-heated zones, as clearly seen in Figure 18 (c).

The smallest eccentricity (ε = 17.5 mm) condition does not show a left tail (see Figure 17) but presents the lowest most probable value for the mean temperature. The intermediate rod eccentricity condition (ε = 25 mm) provides the best heating and temperature homogenization. This choice is a good compromise amongst the three cases.

Temperature gradient evolutions

In order to study and compare the evolution of the temperature gradients for different stirring protocols, we followed the temperature scalar dissipation indicator, which was defined as:

χ g = 1 S tot Stot ∇T * 2 dS = 1 c A c c A c ∇T * 2 c ( 19
)
where S tot is the total surface of the fluid in the mixer and A c is the surface of a mesh element.

In Figure 19 convergence of the temperature field to a recurrent spatial pattern characterized by a strange eigenmode (the patterns described in section 5.3 above). For the two other stirring protocols (continuous and the low modulation period of stirring configuration 3), the decay rates of the temperature gradients are slower than exponential. In the case of continuous rotation of the walls, the temperature gradients near the walls are imposed by the existence of the closed streamlines parallel to these boundaries. The temperature scalar transport evolves in these wall regions only by diffusion across streamlines, and the temperature gradient runs perpendicular to the streamlines. This mechanism has as a consequence the inhibition of thermal chaotic mixing; it is similar, in a certain way, to that observed by Gouillart et al. [START_REF] Gouillart | Walls Inhibit Chaotic Mixing[END_REF][START_REF] Gouillart | Slow decay of concentration variance due to no-slip walls in chaotic mixing[END_REF] for the decay of the concentration of a scalar blob submitted to chaotic flow in a closed vessel with a no-slip boundary condition at the outer wall. However, in our case, continuous rotation of the wall is responsible for the non-asymptotic regime of thermal scalar gradient decay.

Contrary to the case of the dissipation of a concentration scalar (a blob of dye for example), which shows two consecutive phases, a phase of a creation of the gradients and a phase of their dissipation [START_REF] Rothstein | Persistent patterns in transient chaotic fluid mixing[END_REF], only a phase for scalar dissipation exists for temperature scalar gradients. The main factor responsible for this difference is that thermal diffusion is much higher by comparison to molecular diffusion. Thus, thermal equilibrium is rapidly reached and temperature scalar striations are not conserved, contrary to concentration striations that remains in fluid and are thinner and numerous as well as stretching and folding operates. Nevertheless, the global structure of this field has a self-similar pattern. While the temperature field presents striations with the same spatial scale, even if their scalar magnitude is decreasing, the concentration field exhibits striations of decreasing spatial scale before the diffusion makes them disappear. Thus, there is a first phase of scalar gradient creation in the case of concentration mixing that is not present in the case of thermal mixing.

Conclusion

We have numerically investigated the coupled mixing and heating performances induced by chaotic advection in a 2D two-rod mixer. Specific geometrical (eccentricity) and dynamical (stirring protocols, period length) parameters have been identified that promote a high degree of thermal mixing. We have shown the importance of the near wall flow topology for the development of efficient scalar temperature transport inside the mixer. For better understanding of the enhancement mechanism, different statistical tools and mixing and energy indicators have been used. A new composite indicator has been defined as the ratio of the mean temperature of the fluid (which represents the energy extracted from the walls) to the standard deviation of the fluid temperature (which corresponds to the level of homogenization of the scalar). We have also found that chaotic thermal mixing gives rise to singular recurrent spatial structures of the normalized non-dimensional temperature, while the mixing contrast tends to decay slowly with time. The main conclusion of our work is that maximizing heat transfer from the wall boundaries requires that the walls (i.e., tank or rods here) move alternately. In this way, one avoids the development of closed streamlines near the wall, which reduce the temperature gradients at the wall-fluid interfaces, and prevents the build-up of a wall-confined hot fluid zone. Continuous modulation of the wall is not sufficient; to produce effective chaotic mixing, the existence of stagnation points on the walls is necessary to create separatrices from which heteroclinic tangling can develop. Thus, the main advice to retain for industrial applications is the imperative to alternate the rotation of the stirring rods and the cylinder tank with non-continuous modulation (the rods are stopped when the cylinder tank rotates and vice-versa). Further studies will explore the constant wall heat flux boundary condition to study the influence of the thermodependence of the fluid physical properties on the thermal chaotic mixing. We will also address the problem of scalar temperature transport for rheologically complex fluids inside this thermal chaotic mixer.

Figure 1 :

 1 Figure 1: Sketch of the two-rod mixer.

Figure 2 :

 2 Figure 2: Temporal modulation of the angular velocities of the rods and tank for continuous and non-continuous stirring cases. Here the modulation period, τ , is 30 s.

Figure 3 :

 3 Figure 3: A computational cell and its neighbours.

Figure 4 :

 4 Figure 4: Computational mesh of the mixer.

Figure 5 :

 5 Figure 5: Generic steady-state streamline flow topologies. Upper half: without rotation of the outer cylinder: (a) counter-rotation and (b) co-rotation of the rods and (c) rotation of the outer cylinder tank without rod rotation. Lower half: with rotation of the outer cylinder for the stirring configurations 1 (d), 2 (e) and 3 (f) given in Table 2 (from left to right).

Figure 6 :

 6 Figure6: Temporal evolution of the standard deviation and mean value of the temperature for stirring configuration 3 and three wall velocity modulation conditions: non-modulated, continuous and non-continuous (see Figure2).

Figure 7 :

 7 Figure 7: Dimensionless temperature field at t = 120 s for stirring configuration 3 and different flow modulations: (a) non-modulated, (b) continuous modulation and (c) non-continuous modulation (period of modulation τ = 30 s).

Figure 8 : 3 Figure 9 :

 839 Figure 8: Evolution of flow patterns for stirring configuration 1 and non-continuous modulation of the wall velocities at different times during one period of the mixing process. Upper half streamlines. Lower half: dimensionless temperature fields. (Period of modulation τ = 30 s).

Figure 10 :

 10 Figure 10: Composite mixing indicator AT σ as a function of the modulation period for a total mixing time of 1600 s.

Figure 11 :

 11 Figure 11: Probability distribution functions (PDFs) of T * at two different times for τ = 15 s.

Figure 12 :

 12 Figure 12: PDFs of the rescaled dimensionless temperature X for stirring configuration 2 at different times for τ = 15 s.

Figure 13 :

 13 Figure 13: Snapshots showing the evolution of the temperature field towards a strange eigenmode for stirring configuration 2 and modulation period τ = 30 s.

Figure 14 .

 14 The convergence to a thermal strange eigenmode is confirmed. The asymptotic transport properties of the time periodic flow are governed by this eigenmode.

Figure 14 :

 14 Figure 14: Vertical cut-lines of the dimensionless temperature at periodic times: 40, 70, 100, 130, 430 and 730 s (stirring configuration 2, τ = 30 s).

Figure 15 : 20 -Figure 16 :

 152016 Figure 15: Snapshots showing the pattern of the strange thermal eigenmode for stirring configuration 2 and a modulation period τ = 30 s. In these plots, the temperature scalars are rescaled between T * min and T * max .

  , we present the time evolution of this indicator for different stirring configurations. For two stirring protocols (stirring configuration 2 and stirring configuration 3 for τ = 60 s) we observe perfectly exponential decreases. These cases correspond to the rapid

Figure 17 :

 17 Figure 17: PDFs of non-dimensional temperature for stirring configuration 1 and three eccentricities: ε = 17.5 mm, ε = 25 mm and ε = 30 mm (τ = 30 s).

Figure 18 :

 18 Figure 18: Dimensionless temperature fields for increasing mixer eccentricity at t = 220 s for stirring configuration 1. (a) ε = 17.5 mm (b) ε = 25 mm and (c) ε = 30 mm (τ = 30 s).

  Stir. conf. 2 τ=60s, Stir. conf. 3 τ=15s, Stir. conf. 3 τ=30s, Stir. conf. 3 Continuous

Figure 19 :

 19 Figure 19: Evolution of the scalar temperature dissipation indicator for different stirring configurations.

Table 1 :

 1 . Fluid properties.

	Dynamic viscosity (µ)	1.5 P a.s
	Density (ρ)	990 kg.m -3
	Thermal conductivity (k) 0.15 W.m -1 K -1
	Specific heat (c p )	1000 W.kg -1 K -1
	Prandtl number (Pr)	10 4

Table 2 :

 2 The three possible stirring configurations related to the sign of angular velocity.


Weijun Zhou
email: zhouwj@zucc.edu.cn

Member, IEEE Ning Liu
email: ning.liu@femto-st.fr

Member, IEEE Yongxin Wu
email: yongxin.wu@femto-st.fr

Ramirez Estay

Member, IEEE Yann Le Gorrec
email: yann.le.gorrec@femto-st.fr

Member, IEEE Hector Ramirez
email: hector.ramireze@usm.cl

Energy-Based Modeling and Hamiltonian LQG Control of a Flexible Beam Actuated by IPMC Actuators

Keywords: Port-Hamiltonian system, IPMC, flexible beam, passivity-based LQG control VOLUME, 2016

I. INTRODUCTION

The port-Hamiltonian system (PHS) formulation and passivity-based control have been widely used and demonstrated to be effective for the modeling, analysis and control of nonlinear systems [START_REF] Maschke | Port controlled Hamiltonian systems: modeling origins and system theoretic properties[END_REF], [START_REF] Maschke | On the Hamiltonian Formulation of Nonholonomic Mechanical Systems[END_REF]. The classical port Hamiltonian system was first introduced for the finite-dimensional system described by ordinary differential equations [START_REF] Duindam | Modeling and Control of Complex Physical Systems -The Port-Hamiltonian Approach[END_REF] with application to mechanical, electrical and chemical systems using port energy variables. Later, the authors in [START_REF] Van Der Schaft | Hamiltonian Formulation of Distributed Parameter Systems with Boundary Energy Flow[END_REF] discussed the infinite-dimensional port-Hamiltonian system with a new geometric structure called Stokes-Dirac structures. The wellposedness of the infinite-dimensional PHS was investigated in [START_REF] Le Gorrec | Dirac Structures and Boundary Control Systems Associated with Skew-symmetric Differential Operators[END_REF], which shows that the passive propriety is very helpful to prove the well-posedness of the PHS.

It is of great interest to use the PHS formalism for complex multi-physics systems. Under the PHS framework, physical natures are considered: energy conservation, dissipation and energy exchanges between system components. The PHS enables the direct and clear interconnection of different components of the complex system via energy ports. Furthermore, when associated with structure-preserving model reduction and passivity-based control techniques, the port-Hamiltonian framework with physical interpretation is of great interest. Many studies have been dedicated to passivity-based control for multi-physics physical systems [START_REF] Macchelli | Energy Shaping of Distributed Parameter Port-Hamiltonian Systems Based on Finite Element Approximation[END_REF]- [START_REF] Liu | Energy based modeling of ionic polymer metal composite actuators dedicated to the control of flexible structures[END_REF].

Ionic polymer-metal composites (IPMCs) are widely used because of their core advantages: ease of fabrication, fast response, large strain, and low actuation voltage [START_REF] Shahinpoor | Ionic polymer-metal composites: I. fundamentals[END_REF]. Thus, IPMCs have been increasingly used in various domains in recent decades. Of particular interest is the application of IPMCs in medical endoscopes. The use of smart material provides additional degrees of freedom, which contribute to the avoidance of irreversible damage and alleviation of suffering. In [START_REF] Chikhaoui | Kinematic Modeling of an EAP Actuated Continuum Robot for Active Micro-endoscopy[END_REF], a micro-endoscope actuated by IPMCs to drive the bending movement was proposed for endonasal skull base surgery.

The proposed IPMC actuated endoscope system is composed of two components: the electric dynamics of IPMCs described by an ordinary differential equation and the mechanical flexible beam described by a set of partial differential equations. Different physical properties of IPMCactuated endoscopes lead to a complex multi-physical sys-tem problem. Port-Hamiltonian modeling of an IPMC patch actuator was proposed in [START_REF] Liu | Energy based modeling of ionic polymer metal composite actuators dedicated to the control of flexible structures[END_REF]. The mechanical structure is a polyethylene polytube, modelled as a flexible beam using Timoshenko beam theory in a first approach. The complexity of the multi-physics system motivates the use of the port-Hamiltonian framework. Using the PHS formulation, all physical properties of the complex system, such as the physical nature of IPMC actuators and the mechanical flexibility of the beam, are considered. These components are considered interconnected by the PHS formalism via the energy ports.

Many studies have been dedicated to the control of IPMC actuators and flexible structures. The reader can find some classical control analysis and control design methods for flexible structures in [START_REF] He | Dynamical modeling and boundary vibration control of a rigid-flexible wing system[END_REF], [START_REF] Liu | Modeling and adaptive control for a spatial flexible spacecraft with unknown actuator failures[END_REF]. A reduced-order linearquadratic-Gaussian (LQG) controller was proposed in [START_REF] Braun | Reduced-order lqg control of a timoshenko beam model[END_REF] using an early lumping approach. However, the closed-loop stability was not investigated, and only numerical results were obtained. Of particular interest are the results in [START_REF] Wu | Reduced order LQG control design for port hamiltonian systems[END_REF], where a reduced-order LQG control design technique for high-dimensional port-Hamiltonian systems was proposed. The results are based on the LQG controller [START_REF] Jonckheere | A New Set of Invariants for Linear Systems-Application to Reduced Order Compensator Design[END_REF], [START_REF] King | Reduced Order Controllers for Distributed Parameter Systems: LQG balanced Truncation and an Adaptive Approach[END_REF]; then, an equivalent port-Hamiltonian formulation and structurepreserving model reduction were proposed. In [START_REF] Wu | Reduced order LQG control design for infinite dimensional port hamiltonian systems[END_REF], these results were generalized to the infinite-dimensional case. In [START_REF] Mattioni | Modelling and control of an IPMC actuated flexible structure: A lumped port Hamiltonian approach[END_REF], a nonlinear lumped parameter system of a flexible structure using the PHS framework of a class of 1-D IPMC actuated flexible structures was proposed. In that work, a state feedback control law with a state observer for practical implementation is proposed. However, with the increase in number of links, the proposed model becomes very complex and challenging to control design. To address these constraints, a classical passivity-based interconnection damping assignment passivity-based control (IDA-PBC) method is proposed based on a distributed parameter model in [START_REF] Zhou | Port-Hamiltonian modeling and IDA-PBC control of an IPMC-actuated flexible beam[END_REF]. However, the early lumped approach is considered in the previous work, i.e., and the control design is based on the discretized model of the distributed parameter system. The closed-loop stability has not been discussed when the designed discretized state feedback control law is implemented in the original distributed parameter system.

Unlike previous works, the proposed paper designs a passivity-based optimal LQG control of the ODE-PDE coupled system under the port Hamiltonian framework. The PDE is used to address the structure's flexible nature, and the ODEs are used to describe the main dynamics of the IPMC actuators. By choosing the weighting operators in a special manner, the designed optimal LQG controller is passive and has a Hamiltonian structure, which can guarantee the asymptotic stability of the closed-loop system even when we apply the reduced-order optimal controller on the original PDE-ODE system. This stability issue is not considered in previous works [START_REF] Braun | Reduced-order lqg control of a timoshenko beam model[END_REF], [START_REF] Wu | Modeling and control of an IPMC actuated flexible beam under the port-Hamiltonian framework[END_REF]. Furthermore, the effectiveness of the proposed control is validated by experimental results. The control law is implemented on the experimental setup to validate the effectiveness of the proposed method and achieve the desired performance of the flexible structure in both mono-actuated and multi-actuated cases. Preliminary results were reported in [START_REF] Wu | Modeling and control of an IPMC actuated flexible beam under the port-Hamiltonian framework[END_REF], where a passive LQG controller was proposed for a mono-IPMC actuated flexible structure. Only simulation results were presented.

In Section II, a PDE-ODE interconnected model in the port-Hamiltonian framework of a one-dimensional IPMC actuated flexible structure is presented. In Section III, an LQG and damping injection-based reduced-order controller is proposed to improve the dynamics of the flexible beam system. The closed-loop stability analysis is also illustrated in the same Section III. An experimental setup is applied to validate the proposed model and the proposed control law in Section IV. Some conclusions and perspectives of future works are given in SectionV.

II. IPMC ACTUATED FLEXIBLE BEAM

We modelled the IPMC actuated endoscope using the PHS framework. The considered flexible endoscope is on the 1-D spatial coordinate z from a to b, as depicted in Fig 1 . The IPMC actuator patches are attached to the endoscope and subsequently actuated to control its configuration.

Inner tube Flexible beam

A. FLEXIBLE BEAM WITH DISTRIBUTED CONTROL

Endoscopes are used for high-frequency scanning tasks in medical diagnosis. IPMC actuators are attached to the structure, which leads to a composite material structure for the overall system. When a high-frequency mode is excited and the composite material structure is considered, the shear deformation cannot be neglected [START_REF] Korayem | Kinematic and dynamic modeling of viscoelastic robotic manipulators using timoshenko beam theory: Theory and experiment[END_REF]- [START_REF] Ghariblu | Trajectory optimization of flexible mobile manipulators[END_REF]. Hence, the mechanical behaviour of the endoscope is described by a distributed Timoshenko beam model (PDE). An actual endoscope may undergo large deflection. However, in the first instance, we focus on small deformations and control design and keep the modeling as simple as possible. An experimental setup is used to validate the results, for which a Timoshenko beam model is sufficient to handle the primary dynamical behaviour. The mathematical modeling for beams undergoing large deflections can be found in [START_REF] Malatkar | Nonlinear Vibrations of Cantilever Beams and Plates[END_REF] and its PHS formulation in [START_REF] Nishida | Distributed port Hamiltonian formulation of flexible beams under large deformations[END_REF]. The extension of the linear beam model to the large deflection case using the port Hamiltonian formulation is direct.

The flexible beam is defined on the 1-D spatial domain z ∈ [a, b] and driven by the distributed bending moment, denoted by u d .This bending moment is generated by the IPMC patches when we apply the applied voltage. We consider the Timoshenko beam model [START_REF] Macchelli | Modeling and control of the Timoshenko beam. the Distributed Port Hamiltonian approach[END_REF], [START_REF] Jacob | Linear Port-Hamiltonian Systems on Infinitedimensional Spaces[END_REF] with distributed input and power conjugate output y d under the port Hamiltonian framework to be described as:

ẋ = (J -R) Lx + Bu d y d = B * Lx (1)
where

J = P 1 ∂ ∂z + P 0 is a skew symmetric differential operator on the state space X = L 2 ([a, b]; R 4 ) 1 and P1 =    0 1 0 0 1 0 0 0 0 0 0 1 0 0 1 0    , P0 =    0 0 0 -1 0 0 0 0 0 0 0 0 1 0 0 0    (2)
The state variable vector of the system x = [x 1 , x 2 , x 3 , x 4 ] T ∈ X is defined as follows: x 1 = ∂w ∂z (z, t)φ(z, t) is the shear strain, x 2 = ρ(z) ∂w ∂t (z, t) is the transverse momentum distribution, ∂φ ∂z (z, t) is the angular strain, and x 4 = I ρ ∂φ ∂t (z, t) is defined as the angular momentum distribution for z ∈ (a, b), t ≥ 0 with transverse displacement w(z, t) and rotation angle of the beam φ(z, t). We consider the total mechanical energy of the flexible beam as a Hamiltonian, which is defined by the energy variables:

H b (x) = 1 2 b a (Kx 2 1 + 1 ρ x 2 2 + EIx 2 3 + 1 Iρ x 2 4 )dz = 1 2 b a x(z) T (Lx)(z)dz. (3)
The energy of the flexible beam can also be presented in the norm form as

H b (x) = 1 2 x 2
L with the operator L, which is defined as

L = diag K 1 ρ EI 1
Iρ .

The coefficients ρ and I ρ are the mass density and mass moment of area, respectively. E, K and I are the Young's modulus, shear modulus, and moment of area, respectively.

Let us define the dissipation operator as:

R =     0 0 0 0 0 R t 0 0 0 0 0 0 0 0 0 R r    
where constants R t and R r represent the translational and angular viscous fractions, respectively. An important advantage of the port Hamiltonian formulation (1) is that one can prove the existence of a solution only by a proper parametrization of the co-energy variable Lx(z, t) at the boundary a and b [START_REF] Le Gorrec | Dirac Structures and Boundary Control Systems Associated with Skew-symmetric Differential Operators[END_REF] and a simple matrix condition. The boundary port variables can be defined by

f ∂ e ∂ = 1 √ 2 P 1 -P 1 I I Lx(b, t) Lx(a, t) (4)
1 L2 ([a, b]; R n denotes the vector space of 2 integrable functions.

From [START_REF] Le Gorrec | Dirac Structures and Boundary Control Systems Associated with Skew-symmetric Differential Operators[END_REF]Theorem 4.1], we can choose the boundary input as

u = W B f ∂ e ∂ ( 5
)
where W B is full rank and W B ΣW T B ≥ 0 with Σ = 0 I I 0 , then the operator Ax = (J -R) Lx with domain

D (A) = Lx ∈ H 1 [a, b] ; R 4 f ∂ e ∂ = KerW B
generates a contraction semigroup on X.

In the following, we consider the boundary condition according to Fig. 1. The flexible beam is clamped at one side (a) and free at the other side (b). Thus, at the fixed side, the velocity and angular velocity are zero; at the free side, the force and moment are zero. The boundary conditions are defined as the boundary input:

u b = v(a) w(a) F (b) T (b) T . ( 6
)
where the matrix W B is chosen as:

W B = 1 √ 2     -1 0 0 0 1 0 0 0 0 0 -1 0 0 0 1 0 0 1 0 0 0 1 0 0 0 0 0 1 0 0 0 1     (7)
which satisfies the condition W B ΣW T B ≥ 0. According to the boundary conditions, the power conjugate outputs at the boundary are the reaction forces F (a) and T (a) at the fixed side and the velocity and angular velocity at the free side.

y b = F (a) T (a) -v(b) -w(b) T . (8)
The distributed control of the beam is defined as follows. The input mapping over the spatial domain is denoted bt B : C i → X; u d ∈ C i is the distributed moment density applied to the beam 2 , and y d ∈ C i is the distributed angular velocity, which is the power-conjugated output of u d . Each distributed input is defined over the i-th interval on the spatial space I bi = [α i , β i ] and denoted by operator b i (z)u di (t), where b i (z) = 1 if z ∈ I bi and b i (z) = 0 elsewhere and the index i ∈ {1, 2, • • • , m} with m actuators attached to the beam. The mean angular velocity is the output over the same intervals

f di = y di = b a b i (z) 1
Iρ x 4 dz. Hence, we define the distributed input of the flexible beam as:

Bu d =     0 0 0 b(z)     u d (t) (9)
where B :

C m → X, b(z) = [b 1 (z), • • • , b m (z)] and u d (z) = [u d1 (z), • • • , u dm (z)] T .
We define its power conjugate variable as the output variable as y d = B * Lx. Then, the energy balance of the system can be computed as:

∂H b ∂t ≤ y T d u d .

B. IPMC ACTUATOR MODEL

The dynamics of IPMC actuators can be separated into three different physical phenomena: the electric part, the diffusion phenomena and the mechanical deformation dynamics. The first part comes from the electric distribution over the double electric layers between the electrodes. The diffusion phenomena are the water molecules and caution flux due to different potentials caused by the applied voltage on the double electric layers. The last dynamics of the actuators are the mechanical contributions due to their flexible structure. The main physical deformation induced by the applied voltage is mainly caused by the diffusion of cation flux between the double electrode layers (see Fig 2) [START_REF] Shahinpoor | Ionic polymer-metal composites: I. fundamentals[END_REF].

In this work, we assume that the IPMC actuator and beam are perfectly interconnected. Hence, the mechanical contribution of the IPMC patch is assumed to be a part of the Timoshenko beam model. The IPMC electric dynamics are modelled by a simplified lumped RLC control-oriented mode [START_REF] Gutta | Modeling of ionic polymer metal composite actuator dynamics using a large deflection beam model[END_REF]. The torque generated by the IPMC patch is proportional to the applied voltage.

The port-Hamiltonian formulation of the equivalent electrical model is given by:

ϕ = [ϕ 1 • • • ϕ m ] T ∈ R m , Q = [Q 1 • • • Q m ] T ∈ R m .
We define the dissipation matrices R 1 ∈ R m×m and R 2 ∈ R m×m of the IPMC actuator with resistances r 1i and r 2i of each equivalent IPMC model:

R 1 =      r 11 0 • • • 0 0 r 12 • • • 0 . . . . . . . . . . . . 0 0 • • • r 1m      , (11)
R 2 =      1 r21 0 • • • 0 0 1 r22 • • • 0 . . . . . . . . . . . . 0 0 • • • 1 r2m      . ( 12
)
The Hamiltonian of the m IPMC actuators is their electric energy defined as:

H a = 1 2 Q 2 C + 1 2 ϕ 2 L ( 13
)
with capacitance matrix C ∈ R m and inductance matrix L ∈ R m :

C =      C 1 0 • • • 0 0 C 2 • • • 0 . . . . . . . . . . . . 0 0 • • • C m      , (14)
L =      L 1 0 • • • 0 0 L 2 • • • 0 . . . . . . . . . . . . 0 0 • • • L m      . ( 15
)
The input variable u is the applied voltage on two electrodes of the actuators, and u a is the current input, which is generated from the mechanical deflection of the actuator structure. From the power conjugate viewpoint, their associated outputs are the current y = ∂Ha ∂ϕ over inductance L and voltage y a = ∂Ha ∂Q of capacitor C. According to the state variables and input-output definitions, the interconnection relation between the flexible beam dynamics and the IPMC patch dynamics is defined as:

u d u a = 0 +k -k 0 y d y a , k =      k 1 0 • • • 0 0 k 2 • • • 0 . . . . . . . . . . . . 0 0 • • • k m     
The bending moments generated by the IPMC actuators are proportional to voltage y a with constant coefficients k i . Since the interconnection is power conservative, the applied current u a on the capacitor is caused by flexible structure movement.

Considering the IPMC actuator equation [START_REF] Liu | Energy based modeling of ionic polymer metal composite actuators dedicated to the control of flexible structures[END_REF] and flexible equation (1), using the power preserving interconnections, the total system is defined by:

ẋ = J -R Bk -k T B * J -R J-R ∂H ∂x +   0 I m 0   u y = 0 I m 0 ∂H ∂x , (16) 4
VOLUME , 2016

where u, y ∈ R m , x = [x, x a ] T ,

J -R = 0 -I m I m 0 - R 1 0 0 R 2 ,
where 0 is the zero matrix with the appropriate dimension. The state space of the coupled system is

X = L 2 ([a, b]; R 4 ) × R 2m .
Due to the energy-preserving interconnections (II-B), the energy of the interconnected system can be rewritten as the sum of the energies of the two parts:

H(x, x a ) = H b (x) + H a (x a ) = 1 2 x T Qx = 1 2 x 2 L + 1 2 Q T C -1 Q + 1 2 ϕ T L -1 ϕ The energy matrix of the coupled PDE-ODE system is de- fined by Q = diag[L, L -1 , C -1 ].

III. CONTROL DESIGN BY THE LQG METHOD WITH DAMPING INJECTION

The control goal is to obtain a desired configuration of the flexible beam with the desired performance by regulating the input voltages of the IPMC patches. A control law of the endoscope system ( 16) using an LQG control design coupled with damping injection is proposed. The closed-loop control schema is shown in Fig 3.

A. LQG AND DAMPING INJECTION CONTROL DESIGN

The design is divided into two parts. First, we consider an LQG controller to eliminate the vibration of the endoscope. By definition, LQG control is a combination of linear-quadratic state estimation (LQE) and a linear quadratic regulator (LQR). Thus, the disadvantage is that it has the same order as the system itself. In our case, the LQG controller is also a coupled PDE-ODE system as the system [START_REF] Wu | Reduced order LQG control design for port hamiltonian systems[END_REF]. To implement this type of controller, its order must be reduced. However, when the lower-order (finite-dimensional) controller is applied to the PDE system, closed-loop stability cannot be guaranteed in the general case, which is known as the spillover effect. Hence, in this work, an LQG controller using the port-Hamiltonian formulation and model reduction method is applied to obtain a lower-order controller.

This controller can be easily implemented in real physical systems, while closed-loop stability is simultaneously guaranteed. Second, damping injection is applied to reduce the settling time of the system.

The LQG controller can be formulated as an observerbased dynamical feedback:

ẋ = ((J -R) Q -BK -F B * Q) x + F u c y c = K x ( 17
) with the LQG controller state x, filter gain F and optimal state feedback gain K. The two gains F and K can be computed as:

K = R-1 B * P c and F = P f QBR -1 w
with the unique solutions P c = P * c > 0 and P f = P * f > 0 of the two Riccati equations [START_REF] Curtain | Introduction to Infinite-Dimensional Systems Theory: A State-Space Approach[END_REF]:

(J -R) QP f +P f Q (J -R) * -P f QBR -1 w B * QP f +Q v = 0 (18) Q (J -R) * P c + P c (J -R) Q -P c B R-1 B * P c + Q = 0 (19) Q v and
R w denote the state and output white noise, respectively. The weighting operators Q and R are taken from the optimal control problem, which minimizes the energy cost function:

J c = +∞ 0 x 2 Q + u 2 R dt. ( 20
)
Remark 1. The two Riccati equations ( 18) and ( 19) are well posed and have unique positive definite solutions if the following holds. (i) The pair

((J -R)Q, Q 1/2 v
) is exponentially stable. (ii) The pair ( Q1/2 , (J -R)Q) is exponentially detectable [START_REF] Curtain | Introduction to Infinite-Dimensional Systems Theory: A State-Space Approach[END_REF]. By using the energy function as the Lyapunov function and considering the dissipation term R, these two conditions can easily be verified on the system [START_REF] Wu | Reduced order LQG control design for port hamiltonian systems[END_REF].

In general, LQG controllers are not passive, and we cannot preserve the Hamiltonian structure when the controller is interconnected with the plant system. However, when we take the exceptional choices of weighting operators Q and R and covariance operators Q v and R w , the closed-loop system coupled with the LQG controller has a port-Hamiltonian system, as presented in [START_REF] Wu | Reduced order LQG control design for infinite dimensional port hamiltonian systems[END_REF]:

Theorem 2 (Hamiltonian LQG method). If the control weighting operators R, Q and covariance operators R w , Q v are chosen such that R = R w . (21)
Q v z = Q -1 2QJ * P c + 2P c JQ + Q Q -1 z, (22)
the LQG controller ( 17) is passive and has a port-Hamiltonian realization. Moreover, the two unique solutions of the operator equations ( 18) and ( 19) are related by

Q -1 P c = P f Q (23)
The above theorem results in an LQG controller with a port Hamiltonian formulation. Using this Hamiltonian LQG VOLUME , 2016 control formulation, we can preserve the port Hamiltonian structure in the closed-loop system. Furthermore, this passive LQG method can provide a balanced reduction coordinate because P c P f = I. This relation implies that the state space can be separated into two parts based on different importance values of their contributions to the controller design. By taking this advantage, the design and reduction of the controller can simultaneously proceed. We define the balanced reduction coordinate for the port-Hamiltonian system as follows: Definition 3. If two Riccati equation solutions P f and P c satisfy

P f = P c = Σ = diag(σ n ) n∈N ∈ L( 2 ), (24)
where (σ n ) n∈N is a positive and non-increasing sequence with

σ 1 > σ 2 > • • • > σ n > • • • > 0.
Then, we say that the port Hamiltonian system (16) is Hamiltonian LQG balanced.

Define the transformation operator T that can diagonalize P c and P f as:

T P f T * = T + P c T + * = Σ. (25)
Then, we can denote the realization of the port-Hamiltonian system (16), i.e., Hamiltonian LQG balanced realization, as:

ẋb = (J b -R b ) Q b x b + B b u y = B * b Q b x b with x b ∈ 2 .
(26) To obtain the reduced-order system of the balanced realization [START_REF] Malatkar | Nonlinear Vibrations of Cantilever Beams and Plates[END_REF] and the reduced Hamiltonian LQG controller while preserving their passivity and Hamiltonian structure, we use the Petrov-Galerkin projection method. The details of this method are shown in [START_REF] Harkort | Stability and Passivity Preserving Petrov-Galerkin Approximation of Linear Infinite-Dimensional Systems[END_REF]. The reduced port-Hamiltonian system is written as:

ẋr = (J r -R r ) Q r x r + B r u y = B T r Q r x r with x r ∈ R n (27)
The reduced-order LQG controller is obtained using the reduced-order system [START_REF] Nishida | Distributed port Hamiltonian formulation of flexible beams under large deformations[END_REF] and Theorem 2 and formulated as a strictly passive PH control system:

ẋc = (J c -R c ) Q c x c + B c u c y c = B T c Q c x c with x c ∈ R n (28
) Then, this controller ( 28) is implemented in the original coupled PDE-ODE system [START_REF] Wu | Reduced order LQG control design for port hamiltonian systems[END_REF] to eliminate the vibration of the endoscope using the passive preserving interconnection:

u c = y u = -y c ( 29
)
To improve the response time performance, output feedback damping injection [START_REF] Van Der Schaft | L2-gain and Passivity Techniques in Nonlinear Control[END_REF] is employed:

u = -r c y, (30)
The output of system y is the current of the IPMC patch, which can be easily measured. The main purpose of using damping injection is to improve the response performance.

The dynamics of the system can be accelerated using positive damping injection, i.e., the dissipation parameter of the controller r c > 0. Furthermore, this parameter must be bounded below the natural damping coefficient r c > -r 1 to guarantee the stability of the closed-loop system. The oscillation of the flexible beam becomes more important if the response dynamics are faster. The LQG controller will be combined with the damping injection to find the desired compromise between the vibration flexible beam and the response time.

B. CLOSED-LOOP STABILITY ANALYSIS

In this part, we analyse the closed-loop stability of the original coupled PDE-ODE system ( 16) with the reduced-order controller [START_REF] Macchelli | Modeling and control of the Timoshenko beam. the Distributed Port Hamiltonian approach[END_REF]. First, the existence of a solution is considered in the following theorem. x 2 ≤ u T y and let the controller (28) be strictly passive. Then, the closed-loop system with the interconnection relation ( 29) is defined by

ẇ(t) = J cl w(t), w(0) ∈ X ( 31
)
where X = X × R n is the state space of the closed-loop system, w = x x r ∈ X and J cl : X → X is a linear operator defined by

J cl w(t) =   (J -R) Q BkQ 0 -k T B * Q (J -R) Q -B T r Q r 0 B r Q (J r -R r ) Q r   x x r (32) with D(J cl ) = H 1 ([a, b]; R 4 ) × R 2m×n . ( 33
)
Moreover, the operator A c defined by A c w = J cl w with

D(A c ) = D(J) × R n ( 34
)
generates a contraction semigroup on X.

Proof. The open-loop operator A = J -R of the system ( 16) is a generator of the contraction semigroup [START_REF] Villegas | A Port-Hamiltonian Approach to Distributed Parameter Systems[END_REF], [START_REF] Lamoline | Analysis and LQG control of infinite dimensional Stochastic port-Hamiltonian systems[END_REF].

The closed-loop operator is a consequence of the power preserving interconnection of the system ( 16) with the controller [START_REF] Macchelli | Modeling and control of the Timoshenko beam. the Distributed Port Hamiltonian approach[END_REF]. Then, from the Lumer-Phillips theorem, the closed-loop operator A c generates a contraction semigroup.

The stability of the closed-loop system [START_REF] Curtain | Introduction to Infinite-Dimensional Systems Theory: A State-Space Approach[END_REF] is considered in the following theorem.

Theorem 5. Let the controller (28) be strictly positive real. Then, the closed-loop system (31) is globally asymptotically stable, i.e., for any initial condition w(0) ∈ X, the unique solution of (31) asymptotically approaches zero, i.e., lim t→∞ w(t) X = 0.

Proof. Since the controller ( 28) is strictly positive real and the system (16) satisfies 1 2 d dt

x 2 ≤ u T y, the closed-loop system (31) has a compact resolvent [START_REF] Villegas | A Port-Hamiltonian Approach to Distributed Parameter Systems[END_REF]Thm 2.25]. We can choose the closed-loop energy function as a Lyapunov function candidate, and the asymptotic stability can be proven using the Lyapunov argument in combination with Lasalle's invariance principle. To use Lasalle's invariance principle, we need the solution trajectory of the closed-loop system to be precompact. Since the system is linear, the precompact condition for the trajectory can be reduced so that the closedloop operator has a compact resolvent. In our case, we can prove the norm of operator (λI -A c ) x X ≥ a x X with a > 0 because of the dissipativity of the closed system. Thus, using Theorem 2.25 in [START_REF] Villegas | A Port-Hamiltonian Approach to Distributed Parameter Systems[END_REF], we obtain that the closed-loop system operator has a compact resolvent. From

1 2 x r 2 Qr < -R r x 2
r , the equilibrium profile satisfies u r = y r = 0, and from Remark 1, it reduces to zero.

IV. EXPERIMENTAL METHOD

In this section, we will use the experimental setup to validate the proposed model and show the effectiveness of the proposed control design method. In Subsection IV.A, the flexible beam and the IPMC actuator dynamics are identified and validated from the experimental data. Subsection IV.B uses the positive damping injection method to track the periodic signal reference. Subsections IV.C and IV.D present experimental results of the passive LQG controller implementation to the flexible beam for both mono-and multi IPMC patch actuated cases.

The experimental setup is depicted in Fig. 4. A dSPACE MicroLabBox compact prototyping unit and a computer with MATLAB Simulink are used to acquire the measurements and implement the designed control law in the system by generating a voltage across the IPMC patches. The displacement measurement is acquired by a KEYENCE (LK-G152) laser sensor to identify unknown parameters.

The dimension of the flexible structure L is 160 mm, the width W is 5 mm, and the beam thickness T is approximately 0.20 mm. The mass density ρ is 936kg/m 3 . The inertia moment of area I and angular moment of inertia I ρ can be calculated by physical parameters L and W since the Timoshenko beam model is used, and their values are 4.7 × 10 -15 m 4 and 4.34 × 10 -12 kg/m, respectively. Remark 6. The controlled beam in the experiment is actuated by the IPMC patches, which leads to a composite structure. The actuation of the beam is driven by the torques generated from the IPMC actuators. Hence, the Timoshenko beam assumptions fit quite well with the proposed experimental setup even if the considered beam is quite thin.

The Timoshenko beam model is discretized for the identification procedure in MATLAB ® . A structure-preserving discretization method (mixed finite elements method [START_REF] Golo | Hamiltonian Discretization of Boundary Control Systems[END_REF]) is used to preserve the port-Hamiltonian structure of the model.

Sequential quadratic programming (SQP) and trust-regionreflective algorithms ('fmincon') are used to identify the unknown parameters. This optimal nonlinear model identification algorithm can be found in the identification toolbox, which has been implemented in the MATLAB software ® . 2.

The blocking force of the IPMC has been measured by different researchers and shown in [START_REF] Shahinpoor | Ionic polymer-metal composites: I. fundamentals[END_REF], [START_REF] Gutta | Modeling of ionic polymer metal composite actuator dynamics using a large deflection beam model[END_REF]. It has an almost

F U = 3.75 × 10 -4 N/V ( 35
)
where F is the measured blocking force and U is the applied voltage. The bending moment of the IPMC actuator can be computed by the blocking force as follows [START_REF] Gutta | Modeling of ionic polymer metal composite actuator dynamics using a large deflection beam model[END_REF]:

M = 2L a 3 F (36)
From ( 35) and ( 36), we can find the coupling constant

k i = M U = 7.5 × 10 -6 N.m/V (37)
In the setup in Fig. 4, the flexible structure is actuated by one IPMC patch on the clamped side. Fig. 6 illustrates the identification result of the IPMC actuator with u = 1.5V .

B. TRAJECTORY TRACKING OF THE FLEXIBLE BEAM

Trajectory tracking of the flexible beam tip is an important control problem. In this part, this problem will be investigated by using the experimental setup. We employ the damping injection [START_REF] Gutta | Modeling of ionic polymer metal composite actuator dynamics using a large deflection beam model[END_REF], where output y is the current. The following figure shows the experimental validation of the proposed control law for reference tracking. In Fig. 7, the tracked reference is a periodic signal. The open-loop response cannot reach the reference since the response time of the open-loop system is too slow with respect to the reference. The control law ( 30) is used to track the reference. The solid red line shows the experimental result. We observe that the closed-loop system tracks the reference signal with good precision.

C. CONTROL OF THE IPMC ACTUATED FLEXIBLE BEAM

The settling time for the open-loop system is approximately 11s, as shown in Fig. 6. We illustrate the closed-loop performances with the designed state feedback controls in Section III. To apply the Hamiltonian LQG method defined in Theorem 2, the control weighting operator is

Q = Q * b B b B * b Q
, where the operator and matrix are defined in the balanced realization of the system (26) and R = I. The cost function is defined as shown in equation [START_REF] Mattioni | Modelling and control of an IPMC actuated flexible structure: A lumped port Hamiltonian approach[END_REF], which is the input output energy function in this case. To guarantee the closed loop passivity and stability, we compute the covariance operators R w and Q v using ( 21) and [START_REF] Wu | Modeling and control of an IPMC actuated flexible beam under the port-Hamiltonian framework[END_REF], which are defined in Theorem 2. For implementation purposes, the designed full-order LQG controller is reduced in a structure-preserving manner to the lower-order controller, as shown in [START_REF] Macchelli | Modeling and control of the Timoshenko beam. the Distributed Port Hamiltonian approach[END_REF].

Fig. 8 illustrates the comparison between 2 different control laws for position assignment and open-loop response.

We observe that with simple positive damping injection, the raising time is approximately 6 s instead of 11 s for the open-loop system. However, we observe a significant oscillation of the displacement due to the vibration of the flexible beam, which may gradually irreversibly destroy the flexible structure. By using LQG control plus positive damping injection control, we observe that the response time is reduced to 6 seconds. Furthermore, compared to only using the positive damping injection, the vibration has been significantly reduced.

D. MULTI-ACTUATION OF A FLEXIBLE BEAM

To work in the human body environment, the endoscope should be capable of being driven to the desired shape. In this part, a multi-actuation case with two patches is considered, as shown in Fig. 9. The first IPMC patch is attached at the fixed side of the beam to achieve the largest deflection of the beam's free tip. The second patch cannot be placed too close to the free end because the beam would be deformed by the actuator and its wires. Hence, the patch is attached at the middle of the structure. Then, the distributed input map B of the flexible beam model (1) can be defined as: 3 and b 1 (z) = b 2 (z) = 0 elsewhere. The control weighting operators are defined in the same manner as the mono-actuation case, as shown in IV-C.The control objective in this case is to drive the beam to a desired shape (position) as fast as possible while reducing the beam vibration. The desired shape is presented in Fig. 9. We first use the clamped side IPMC actuator to move the flexible beam such that the beam tip is displaced 5mm. Then, the second IPMC actuator is used to change the beam configuration in the middle. In this case, we use the second actuator to drive the beam tip to 10 mm.

B =     0 0 0 0 0 0 b 1 (z) b 2 (z)     (38) where b 1 (z) = 1 if z ∈ [0, 0.2L] and b 2 (z) = 1 if z ∈ [0.4L, 0.6L]
In Fig. 10, we show the beam tip displacements (red dashed line in the top figure) and middle displacement (solid blue line). The clamped side actuator is activated at 3 seconds, and the second actuator is actuated at 15 seconds. We see that the tip displacement is changed by the two actuators, while the middle displacement depends only on the clamped side actuator. The second actuator has no impact on the middle displacement except for a slight vibration at approximately 15 seconds. Finally, we implement the proposed control law in the same experimental setup. The experimental results illustrated the mono-actuation and multi-actuation cases to show the effectiveness of the presented control methods.

A perspective of this work is to consider the uncertainty of the parameters due to environmental reasons and the disturbance from external perturbations. During the experiment, the actuation of the IPMC actuator is sensitive to the humidity of the working conditions. The robustness of the control law should be considered in future work. Moreover, in [START_REF] Liu | Energy based modeling of ionic polymer metal composite actuators dedicated to the control of flexible structures[END_REF], a complete IPMC actuator model, where the diffusion phenomena were considered, was proposed using the port-Hamiltonian approach. The control design based on this complete actuator model should be considered in the future.

FIGURE 1 .

 1 FIGURE 1. Flexible structure actuated by IPMC

FIGURE 2 .

 2 FIGURE 2. Physical principle of the IPMC actuator and its simplified electrical model

FIGURE 3 .

 3 FIGURE 3. Closed-loop schema

Theorem 4 .

 4 Let the state of the open loop system of[START_REF] Wu | Reduced order LQG control design for port hamiltonian systems[END_REF]

Fig 5

 5 Fig 5 illustrates the result obtained by the identification. With a fitting percentage of 89.69%, the curve fitting of the proposed model and the experimental data are validated. The identified parameters of Young's modulus, shear modulus and two dissipation constants are shown in Table1.

  Fig 5 illustrates the result obtained by the identification. With a fitting percentage of 89.69%, the curve fitting of the proposed model and the experimental data are validated. The identified parameters of Young's modulus, shear modulus and two dissipation constants are shown in Table1.

FIGURE 5 .

 5 FIGURE 5. Identification curve fit of the flexible beam

FIGURE 6 .

 6 FIGURE 6. IPMC actuated compliant endoscope

FIGURE 7 .

 7 FIGURE 7. Reference tracking

FIGURE 8 .

 8 FIGURE 8. Positive damping injection control and LQG+positive damping injection control

3 L

 3 is the length of the flexible beam

FIGURE 9 .FIGURE 10 .

 910 FIGURE 9. Multiple actuations of the flexible beam and the desired shape of the flexible beam

TABLE 1 .

 1 Identified parameters of the flexible beam E Young's modulus 4.14 × 10 9 P a K shear modulus 1.418 × 10 9 P a Rt Traversal viscous fraction 2 × 10 -5 kg.m 3 /s Rr Angular viscous fraction 1 × 10 -5 kg.m/s

TABLE 2 .

 2 Parameters of the IPMC actuator

	C	Capacitance	5.8 × 10 -2 F
	r 1	Resistance r 1	29.75 Ω
	r 2	Resistance r 2	700 Ω
	linear relation with respect to the applied voltage:

VOLUME , 2016

C i is an i-dimensional complex space.VOLUME , 2016

VOLUME , 2016

VOLUME , 2016

VOLUME , 2016

This research was supported by Zhejiang Provincial Natural Science Foundation of China under grant no. LQ21F030008, the EIPHI Graduate School (contract: ANR-17-EURE-0002), ANR IPMACTS project (contract: ANR-21-CE48-0018) and Bourgogne Franche Comté region of France (contract: ANER 2018Y-06145).