Lamberto Dell'elce

Jean-Baptiste Caillau

Jean-Baptiste Pomet

On the convergence of time-optimal maneuvers of fast-oscillating control systems

For a control system with one fast periodic variable, with a small parameter measuring the ratio between time derivatives of fast and slow variables, we consider the Hamiltonian equation resulting from applying Pontryagin maximum principle for the minimum time problem with fixed initial and final slow variables and free fast variable. One may perform averaging at least under normalization of the adjoint vectors and define a "limit" average system. The paper is devoted to the convergence properties of this problem as the small parameter tends to 0. We show that using the right transformations between boundary conditions of the "real" and average systems leads to a reconstruction of the fast variable on interval of times of order 1/ε where ε is the small parameter. This is only evidenced numerically in this paper. Relying on this, we propose a procedure to efficiently reconstruct the solution of the two point boundary problem for nonzero ε using only the solution of the average optimal control problem.

I. INTRODUCTION

Periodic averaging is extremely common. For a differential equation in a certain number of slow variables and one periodic fast variable, with a small parameter ε > 0 measuring the ratio, it states that the distance between solutions of the original and average systems is of order ε on a time-interval of order 1/ε. See for instance [START_REF] Sanders | Averaging methods in nonlinear dynamical systems[END_REF], [START_REF] Arnold | Mathematical methods of classical mechanics[END_REF].

Here, motivated by low-thrust orbital transfer in space mechanics (see also [START_REF] Edelbaum | Optimum power-limited orbit transfer in strong gravity fields[END_REF], [START_REF] Geffroy | Généralisation des techniques de moyennation en contrôle optimal -Application aux problèmes de rendez-vous orbitaux en poussée faible[END_REF], [START_REF] Dargent | Averaging technique in T_3D an integrated tool for continuous thrust optimal control in orbit transfers[END_REF], [START_REF] Bonnard | Riemannian metric of the averaged energy minimization problem in orbital transfer with low thrust[END_REF], [START_REF] Bombrun | Asymptotic behavior of the time optimal orbital transfer for low thrust 2-body control system[END_REF]), we consider minimum-time optimal control for a control system [START_REF] Sanders | Averaging methods in nonlinear dynamical systems[END_REF], and use averaging on the Hamiltonian equation resulting from Pontryagin Maximum Principle. Averaging for optimal control was also studied in [START_REF] Chaplais | Averaging and deterministic optimal control[END_REF].

As already noted in [START_REF] Rouot | Méthodes géométriques et numériques en contrôle optimal et applications au transfert orbital à poussée faible et à la nage à faible nombre de Reynolds[END_REF], one has to be careful in sorting fast and slow variables among the adjoint variables in this equation, and averaging only holds with properly normalized adjoint vectors. Then, convergence is obviously in a weak sense because, contrary to Cauchy problems, two point boundary value problems do not have unique solutions (and in fact, in the problem with fixed slow variables and free fast variables, more and more solutions appear as ε goes to zero).

We study this convergence and give a procedure to recover the solutions of the optimal control problem from solving the average equations only. We implement transformations between initial and final conditions of the slow variables of the average and original system; these are the same as in second order averaging mentioned above (in the context of optimal control they were already discussed in [START_REF] Dargent | Averaging technique in T_3D an integrated tool for continuous thrust optimal control in orbit transfers[END_REF]). We include a reconstruction of the fast variable in the average system.

Although there is no proof per se, we do witness very clearly, in our numerical experiments, a reconstruction of the fast variable, with an error that is possibly not of order 1 with respect to ε, but that tends to zero. Although a second order term in the average dynamics would be needed to achieve this in the classical Cauchy problem, simulations are extremely clear. We then describe an original procedure to recover the initial conditions of the adjoint vectors as well as of the fast variables from solving the average system. This procedure relies on the assumption that reconstruction of the fast variables is effective. One of its interest is that it explains very well the fact that there are more and more solutions as ε → 0. Another and more practical interest is that it provides very accurate initial guesses for the original two point boundary value problem.

II. AVERAGING FAST-OSCILLATING OPTIMAL CONTROL

SYSTEMS

The time-optimal problem of slow-fast dynamical systems, which is referred to as "original problem" hereafter, is introduced first. Then, its average counterpart is defined.

A. Original problem

We consider the following slow-fast controlled dynamics where ε is to be understood as a small parameter, scaling the slow part:

İ = ε F 0 (I, ϕ, ε) + m i=1 u i F i (I, ϕ, ε) , φ = ω(I) + ε G 0 (I, ϕ, ε) + m i=1 u i G i (I, ϕ, ε) , |u| = u 2 1 + • • • + u 2 m ≤ 1, (1)
The slow state variables, I, live on some manifold M of dimension n, and there are fixed endpoint conditions I 0 , I f . Conversely, ϕ is the fast variable, and all the data (that is the vector fields F i and G i ) are 2π-periodic with respect to it. There are no endpoint conditions neither on ϕ ∈ R, nor on its angle class in S 1 R/2πZ. We assume that the pulsation ω(I) is uniformly bounded on M by some positive constant. The control u is prescribed to the unit Euclidean ball of R m . We are interested in the minimum time control of this system. Pontryagin maximum principle implies that time minimizing trajectories are projections on M of integral curves of the Hamiltonian system defined on T * (M → R) by H(I, ϕ, p I , p ϕ , ε) := p ϕ ω(I) + εK(I, ϕ, p I , p ϕ , ε),

where

K := H 0 + m i=1 H 2 i ,
and

H i (I, ϕ, p I , p ϕ , ε) := p I F i (I, ϕ, ε) + p ϕ G i (I, ϕ, ε)
for i = 0, . . . , m. Namely,

İ = ε ∂K ∂p I , ṗI = -p ϕ ∂ω ∂I -ε ∂K ∂I , φ = ω(I) + ε ∂K ∂p ϕ , ṗϕ = -ε ∂K ∂ϕ . (2)
The initial and final phases ϕ(0) and ϕ(t f ) are free, so that we look for integral curves of H such that p ϕ (0) = p ϕ (t f ) = 0 (transversality condition of Pontryagin maximum principle). Hence, the solution of the two-point boundary problem is achieved by finding a zero (t f , p I0 , ϕ 0 ) of the shooting function

S ε (t f , p I0 , ϕ 0 , I 0 , I f ) =   I (t f , I 0 , ϕ 0 , p I0 , 0, ε) -I f p ϕ (t f , I 0 , ϕ 0 , p I0 , 0, ε) |p I0 | -1  
(3) parameterized by I 0 , I f and ε. Here, I (t, I 0 , ϕ 0 , p I0 , 0, ε) and p ϕ (t, I 0 , ϕ 0 , p I0 , 0, ε) denote the unique solution at time t of the Hamiltonian system, Eq. ( 2), with initial conditions I 0 , ϕ 0 , p I0 , p ϕ (0) = 0. The constraint |p I0 | -1 is associated to the unknown t f , i.e., the time needed to reach I f and p ϕ (t f ) = 0 starting from the aforementioned initial conditions. Such normalization is always possible by homogeneity. The compact notation S ε (t f , p I0 , ϕ 0 ) is used in the reminder whenever boundary conditions I 0 and I f are understood.

B. Average problem

The term -p ϕ ∂ω ∂I in the dynamics of the adjoint state p I is not a priori small, but rather of order zero with respect to ε. Its presence therefore prevents us from using standard averaging techniques to approximate this Hamiltonian system when ε goes to zero. Instead, let us introduce the re-scaling s := εt, ψ := εϕ, and define

H(J, ψ, p J , p ψ ) := p ψ ω(J) + K(J, p J ) where K = H 0 + K 0 with K 0 (J, p J ) := 1 2π 2π 0 m i=1 H 2 i (J, ϕ, p J , p ϕ = 0, ε = 0) dϕ, and
H 0 (J, p J ) := p J , 1 2π 2π 0 F 0 (J, ϕ, ε = 0) dϕ .
The "average" Hamiltonian H is smooth on Ω, the open complement of

Σ := (Σ), : T * M × R → T * M, Σ := {(I, p I , ϕ) ∈ T * M × R | (∀i = 1, m) : p I , F i (I, ϕ, ε = 0) = 0}.
So is the original Hamiltonian H, for small enough ε, on the complementary of a thickening of Σ.

The variable ψ does not appear in H, so for any value of the first integral p ψ , one can consider the two point boundary value problem in the rescaled time s (we denote = d/ds),

J = ∂K ∂p J , p J = -p ψ ∂ω ∂J - ∂K ∂J ,
for fixed endpoints J 0 , J f in M . Solving this problem is equivalent of finding a zero (s f , p J0 ) of the average shooting function

S(s f , p J0 , J 0 , J f , p ψ ) = J(s f , J 0 , p J0 , p ψ ) -J f |p J0 | -1
parameterized by J 0 , J f and p ψ . In analogy with the notation used in Eq. ( 3), J(s, J 0 , p J0 , p ψ ) denotes the unique solution (assumed to be defined up to this point) at time s of the average Hamiltonian system with initial conditions J 0 , p J0 and p ψ . Given a solution (s f , p J0 ), a simple integration allows to compute

∆ψ(J 0 , J f , p ψ ) = ψ(s f ) -ψ(0)
because ψ is cyclic. Note that, in general, the shooting problem has several solutions that need to be compared (e.g., by using additional optimality conditions), so that we should explicit the dependence of ∆ψ with respect to the particular solution retained. For the sake of simplicity, we assume that for a given p ψ , there is only one solution, hence the dependence on p ψ of (s f , p J0 ), and so of ∆ψ.

III. CONVERGENCE OF THE TRAJECTORIES

After providing a practical motivation of this study emphasizing difficulties in the analysis of the convergence of a given solution of Problem (3) for ε → 0, we discuss how the small parameter, ε, the average adjoint to the fast variable, p ψ , and the boundary phases, ϕ(0) and ϕ(t f ), relate each others. Then, we outline an algorithm that is capable of approximating all the different solutions of the original problem for a given ε by leveraging on solutions of its average counterpart only. This approximation can serve as a reliable initial guess to solve the original problem.

A. Motivation of the study

The numerical example discussed in Section IV is now used to disclose the fundamental questions that motivated our study (we are only interested to qualitative behaviors in this paragraph; all details of the case study are provided in Section IV): First, we note that the realization of this experiment entails numerical difficulties because t * f = O ε -1 and small ε requires a stringent tolerance of the integrator. Hence, it is not possible to approach ε = 0 too closely. However, Fig. 1 suggests that the convergence to a single solution of the average problem does not occur, and that the continued branch seems to approach to a subset of solutions of the average problem (not a single one) parameterized by p ψ . This is because ϕ * 0 , circulates in S 1 . This observation lends itself to the following question.

Let t * f , p
What p * ψ should be used, and how does it relate to ε and ϕ * 0 ? Inspection of Fig. 1 reveals that average solutions can be parameterized by p * ψ . In addition, the closest average solution to t * f /ε, p * I0 , i.e., the one that should serve as initial guess to the original problem, depends on the value of ε. Modeling the relationship between ε, p * ψ , and ϕ * 0 is the key to predict the evolution of the continued trajectory of Fig. 1.

Should we consider the same boundary conditions for both the original and the average problems?

The relation between ε, p * ψ , and ϕ 0 (discussed in the reminder of this section) requires the refinement of the boundary conditions of the average problem. In fact, using (J 0 , J f ) = (I 0 , I f ) does not yield an accurate enough approximation of the original trajectory, as depicted in Fig. 2.

Here, the change of the initial phase is sufficient to obtain remarkably different trajectories of the original problem. The near-identity transformation discussed in Section III-B is mandatory to generate average trajectories that remain sufficiently close to the original ones.

B. Considerations

Using the same endpoint conditions for the original and average problems, i.e., (J 0 , J f ) = (I 0 , I f ), does not provide with satisfactory matching between original and average trajectories as emphasized above. A crucial step in our approach is to refine the approximation by making an O(ε) correction on these endpoints. To this end, let us define near identity transformations

Î = J + εν I (J, p J , ϕ), pI = p J + εν p I (J, p J , ϕ), pϕ = εp ψ + εν pϕ (J, p J , ϕ),
using the usual ansatz: for each variable X (X equal to I, p I or p ϕ ) of the original system, ν X is searched so that the reconstruction X from the average system is an approximation of order two in ε of its dynamics:

Ẋ(t) = Ẋ(t) + O(ε 2 ).
This leads to the following choices, first for ν pϕ :

∂ν pϕ ∂ϕ = - 1 ω(J) ∂K ∂ϕ (J, ϕ, p J , p ϕ = 0, ε = 0).
Note that K is evaluated at p ϕ = 0, which is compatible with the O(ε 2 ) discrepancy. Then ν pϕ is chosen of zero average. Similarly,

∂ν p I ∂ϕ = 1 ω(J) -ν pϕ (J, p J , ϕ) ∂ω ∂J - ∂K ∂I (J, ϕ, p J , p ϕ = 0, ε = 0) + ∂K ∂J (J, p J ) .
The crucial point in the previous expression is the addition of the term -ν pϕ (J, p J , ϕ) ∂ω ∂J that cannot be neglected. Finally,

∂ν I ∂ϕ = 1 ω(J) ∂K ∂p I (J, ϕ, p J , p ϕ = 0, ε = 0) - ∂K ∂J (J, p J ) .
Again, ν p I and ν I can be chosen of zero average. We now use these transformations to construct an approximate solution of the original problem for a given value of ε.

Transversality conditions (p ϕ (0) = p ϕ (t f ) = 0) are used to determine possible values for the endpoint phases by solving pϕ (0) = pϕ (t f ) = 0 that is, up to some order in ε, p ψ +ν pϕ (J 0 , p J0 , ϕ 0 ) = p ψ +ν pϕ (J f , p J (s f ), ϕ f ) = 0. (4)

These two equations usually have several solutions, all parameterized by p ψ . Let (ϕ 0 , ϕ f ) be one of them. We eventually want to improve our current estimation which is based on a solution of the average system for endpoints J 0 = I 0 , and J f = I f . Actually, we should have solved for

J 0 = I 0 -εν I (I 0 , p J0 , ϕ 0 ), J f = I f -εν I (I f , p J (s f ), ϕ f )
where (s f , p J0 ) is our solution of the average system. Computing the sensitivities of the shooting function S(s f , p J0 , J 0 , J f , p ψ ) with respect to J 0 and J f , an update of order one in ε of its solution can be computed. In particular, ∆ψ is updated according to ∆ψ(I 0 , I f , p ψ ) + εδ∆ψ(I 0 , I f , p ψ , ϕ 0 , ϕ f ).

Note that while we kept implicit the dependence of the update δ∆ψ on the solution (s f , p J0 ) of the average boundary value problem, we made explicit its dependence on the particular choice of solution for the phases ϕ 0 and ϕ f . All in all, assuming that ∆ψ/ε provides an estimate of ϕ(t f ) -ϕ(0), we obtain the following equation, parameterized by p ψ :

∆ψ(I 0 , I f , p ψ ) ε = ϕ f (p ψ ) -ϕ 0 (p ψ )
-δ∆ψ(I 0 , I f , p ψ , ϕ 0 (p ψ ), ϕ f (p ψ )).

C. Methodology

We now propose an algorithm for the generation of reliable initial guesses for the original problem that leverage on solutions of the average system:

1) Solve the average problem

∀ p ψ ∈ p min ψ , p max ψ find s p ψ f , p p ψ J0 s.t.: S s p ψ f , p p ψ J0 , p ψ = 0 2) Find compatible endpoint phases ∀ p ψ ∈ p min ψ , p max ψ find ϕ (j) 0 , ϕ (k) f s.t.: ν p ψ J 0 , p p ψ J0 , ϕ (j) 0 = -p ψ , j = 1, . . . , n p ψ 0 ν p ψ J 0 , p J s p ψ f , ϕ (k) f = -p ψ , k = 1, . . . , n p ψ f
We omit the dependence on p ψ of ϕ

(j) 0 , ϕ (k)
f , and all other variables defined hereafter to simplify the notation.

3) Compute sensitivities

∀ p ψ ∈ p min ψ , p max ψ find δs (jk) f , δp (jk) J0 , δ∆ψ (jk) s.t.: ∂ S ∂ s f δs (jk) f + ∂ S ∂ p J0 δp (jk) J0 = - ∂ S ∂ J 0 δJ (j) 0 - ∂ S ∂ J f δJ (k) f δ∆ψ (jk) = ∂ ∆ψ ∂ s f δs (jk) f + ∂ ∆ψ ∂ p J0 δp (jk) J0 + ∂ ∆ψ ∂ J 0 δJ (j) 0 j = 1, . . . , n p ψ 0 , k = 1, . . . , n p ψ f
where δJ (j) 0

and δJ

(k) f
are evaluated by means of the near-identity transformation as

δJ (jk) 0 = -ν I I 0 , p p ψ J0 , ϕ (j) 0 , δJ (k) f = -ν I I f , p J s p ψ f , ϕ (k) f . 4) Identify solutions at ε Given ε, find p ψ s.t.: mod ∆ψ(I 0 , I f , p ψ ) ε , 2π = ϕ (k) f -ϕ (j) 0 -δ∆ψ (jk)
(5) Values of p ψ solution of Eq. ( 5) and the associated ϕ j 0 , p

p ψ J0 + εδp (jk)
J0 , and s

p ψ f + εδs (jk) f
can be used to initialize the shooting of the original problem. We stress that all steps only use the average system and the near-identity transformation. The small parameter, ε, only appears at the left-hand side of Eq. ( 5), whereas the righthand side is independent of ε.

Concerning the complexity of the methodology: Step 1 implies first the solution of one average problem for p ψ = 0 (this is the only shooting problem that has to be solved without a priori knowledge). Then, differential continuation can be used to obtain the solution in the range p min ψ , p max ψ .

Step 2 involves the zero search of algebraic equations. Step 3 requires the solution of (n + 1)-dimensional linear systems. Finally, step 4 is the identification of intersections of a planar closed curve with a discontinuous function of p ψ . Graphical insight in this last step is provided in Section IV.

IV. CASE STUDY: TIME-OPTIMAL ORBITAL TRANSFER

Orbital mechanics provides an ideal case study for this approach. We consider a Keplerian and co-planar timeoptimal low-thrust transfer. Slow variables are some set of orbital elements defining the instantaneous shape and orientation of the orbit. Slow elements used in this work are I 1 = a, I 2 = e cos(ν),

I 3 = e sin(ν),
where a, e, and ν denote the semi-major axis, eccentricity, and argument of periapsis of the orbit, respectively. The fast variable is an angle positioning the satellite on a precise point of the orbit defined by I. The mean longitude is used for this purpose. Variational equations of these variables are provided by [START_REF] Battin | An Introduction to the Mathematics and Methods of Astrodynamics[END_REF][Chap. 10.4]

İ1 =2ε I 3 1 µ (1 -I 2 2 -I 2 3 ) (I 3 sin σ -I 2 cos σ) u 1 + (1 + I 2 cos σ + I 3 sin σ) u 2 , İ2 =ε I 1 1 -I 2 2 -I 2 3 µ sin σ u 1 + I 2 + (2 + I 2 cos σ + I 3 sin σ) cos σ 1 + I 2 cos σ + I 3 sin σ u 2 ,
İ3 =ε I 1 1 -I 2 2 -I 2 3 µ -cos σ u 1 + I 3 + (2 + I 2 cos σ + I 3 sin σ) sin σ 1 + I 2 cos σ + I 3 sin σ u 2 , φ = I 3 1 µ + ε I 1 1 -I 2 2 -I 2 3 µ I 2 cos σ + I 3 sin σ 1 + 1 -I 2 2 -I 2 3 + 2 1 -I 2 2 -I 2 3 1 + I 2 cos σ + I 3 sin σ u 1 + I 2 cos σ -I 3 sin σ 1 + 1 -I 2 2 -I 2 3 2 + I 2 cos σ + I 3 sin σ 1 + I 2 cos σ + I 3 sin σ u 2 .
where the small parameter, ε, is the thrust-to-mass ratio magnitude, µ is the planetary gravitational constant, and σ is the true longitude. Conversion from true to mean longitude, ϕ = ϕ(σ, I 2 , I 3 ), is achieved by evaluating the Kepler's law for equinoctial elements. Table I lists the parameters of the TPBVP. Controls u 1 and u 2 are radial and tangential thrust, respectively. Figure 3 provides a graphical interpretation of the fourth step of the methodology. The black-dashed curve is the righthand side of the equation, which is independent of ε. Blue curves depict the left-hand term. Colored dots indicate solutions of the original problem. This shows that intersections are excellent initial guesses.

Because ε appears at the denominator of the left-hand term, the blue curve seems to circulate from the bottom to the top when ε is arguably large. Then, the number of solutions (i.e., intersections) increases when ε is very small (the slope of the curve becomes more and more vertical), as shown in Fig. 4, where a very small ε is used. We stress that each red dot of this figure denotes a converged solution of the original problem.

The quality of the estimation of the solution can be appreciated in Fig. 5. Here, solid lines are solutions of the original TPBVP obtained by numerical continuation on the parameter ε (starting from ε = ε 0 we want to reduce it). The Hampath software described in [START_REF] Caillau | Differential continuation for regular optimal control problems[END_REF] is used for this purpose. Dots denote the approximation achieved by means of our approach. The continuation parameter is not monotonically Graphical interpretation of the 4th point of the methodology (solution of Eq. ( 5)). Black-dashed lines are the right-hand side of the equation, which is independent of ε. Solid-blue curve is the left-hand side, ∆ψ/ε. Large dots are the solutions of the original problem. decreasing as shown in the close up of Fig. 5(b). We note that this zoom is centered on the value ε = 0.995 ε 0 , which is used in Fig. 3(b). Here, four intersections are present, which correspond to the four solutions obtained with the numerical continuation for that specific value of ε. Two new solutions are generated when ∆ψ/ε is tangent to the right-hand side curve of Eq. ( 5). We note that colors of the two branches are consistent with the colors of dots in Fig. depicted. Endpoints of the three trajectories coincide. The two solutions of the original problem differ in their boundary phases ϕ 0 and ϕ f . We note that the black dashed trajectory is not the one used to generate the initial guess used to solve the shooting of the original problems. Conversely, initial guesses were generated by using average trajectories with endpoints slightly modified by means of Eq. 4. This correction is mandatory to avoid the deterioration of the guess which was pointed out in Fig. 2.

V. CONCLUSION

This paper offered a methodology for the accurate approximation of time-optimal maneuvers of fast-oscillating systems. The practical interest of this contribution is the relationship that was established between the magnitude of the small parameter ε (read control) and the boundary phases, ϕ 0 and ϕ f . Specifically, the number of possible solutions and their quantitative assessment for a specific ε can be estimated using only the information given by the average solution. In addition, the proposed methodology provides insight into how new combinations of ϕ 0 and ϕ f (i.e., multiple potential local optima) are generated when ε is decreased. This qualitative behavior cannot be fully understood by inspecting solutions of the original problem owing to numerical difficulties related to the integration of trajectories for very low ε.

Fig. 1 .

 1 Fig. 1. Differential continuation of a solution of the original system for decreasing ε (colored curve). The black-dashed lines denote solutions of the average problem.

Fig. 2 .

 2 Fig. 2. Trajectory of the orbital eccentricity for ε = 10 ε 0 . Solid lines are solutions of the original problem. The black-dashed line is the solution of the average TPBVP for p ψ = 0. Initial conditions are the same for all trajectories (only ϕ(0) is different).

  Fig. 3.Graphical interpretation of the 4th point of the methodology (solution of Eq. (5)). Black-dashed lines are the right-hand side of the equation, which is independent of ε. Solid-blue curve is the left-hand side, ∆ψ/ε. Large dots are the solutions of the original problem.

Fig. 4 .

 4 Fig. 4. Solutions found for ε = 0.08ε 0 . Refer to the caption of Fig. 3 for further details.

  3. Finally, Fig 6 depicts the trajectory of the orbital eccentricity (i.e., I 2 2 + I 3 3 ) for ε = 10ε 0 . Both solutions of the original TPBVP and of the average counterpart are

Fig. 5 .

 5 Fig. 5. Convergence of the solution with respect to ε. Solid lines are solutions of the original problem. Black dots denote the estimation of the solutions obtained by means of the proposed algorithm.

Fig. 6 .

 6 Fig. 6. Trajectory of the orbital eccentricity for ε = 10 ε 0 . Solid lines are solutions of the original problem. The black-dashed line is the solution of the average TPBVP for p ψ = 0.

TABLE I PARAMETERS

 I OF THE CASE STUDY.

	Constants
	Gravitational parameter, µ Nominal thrust-to-mass ratio, ε 0	3.986 • 10 14 m 3 s 2 10 -4 N kg
	Initial conditions
	Semi-major axis, I 1	7000 • 10 3 km
	Eccentricity vector x, I 2	0.01
	Eccentricity vector y, I 3	0
	Final conditions
	Semi-major axis, I 1	7500 • 10 3 km
	Eccentricity vector x, I 2	0
	Eccentricity vector y, I 3	0

Acknowledgements. The authors thank Richard Epenoy (CNES) and Thierry Dargent (Thales Alenia Space) for many fruitful discussions.


Ziemihori Ouattara

Albert Magnin

Didier Blésès

Pascal Jay
email: pascal.jay@univ-grenoble-alpes.fr

Influence of the inclination of a plate on forces generated in flows of Newtonian and yield stress fluids

Keywords: Yield stress fluid, Newtonian fluid, Inclined plate, Drag, Lift

The aim of this article is to study experimentally and numerically the influence of the inclination of a plate placed near a wall on flows of Newtonian and yield stress fluids. The drag and lift forces were particularly studied. The inclination of the plate varies between 0° and 90°. The inertia of the fluid is considered negligible. The effects of the Oldroyd number, the shear-thinning index and the inclination of the plate on drag and lift are shown.

Experimentally, an instrumented flume was used with a Carbopol gel as a yield stress fluid model. The measured lift and drag forces are compared to the numerical results and the literature data. The elasticity and plasticity effects are discussed.

Introduction

Many industrial systems involve flows around fixed or mobile obstacles of various geometries and in particular plates. For example, flows in mixers or exchangers which are widely used in chemical engineering. These must be properly designed, and this cannot be achieved without a good understanding of the interactions. The creeping flow around a plate has been widely studied for Newtonian fluids. However, in the domain of non-Newtonian fluids and, in particular, in the case of yield stress fluids, very few studies seem to exist. To our knowledge, the determination of the forces generated on a plate in a yield stress fluid as a function of its inclination angle has never been studied both experimentally and numerically.

For this flow configuration and for a Newtonian fluid, In et al. [START_REF] In | Two-dimensional viscous flow past a flat plate[END_REF] carried out a 2D numerical study of flow around a flat plate with a variable incidence angle and Reynolds numbers varying between 1 and 30. They calculated drag and lift coefficients as a function of plate angle. In addition, several studies were carried out for the flow of a Newtonian fluid parallel [START_REF] Tamada | Low-Reynolds-number flow past a cylindrical body[END_REF][START_REF] Janssen | Flow past a flat plate at low Reynolds numbers[END_REF] and perpendicular [START_REF] Tomotika | The steady flow of a viscous fluid past an elliptic cylinder and a flat plate at small Reynolds numbers[END_REF][START_REF] Hudson | The flow of a viscous incompressible fluid past a normal flat plate at low and intermediate Reynolds numbers: the wake[END_REF] to the plate in the infinite domain.

On the other hand, for a power-law fluid or shear-thinning fluid, Wu and Thompson [START_REF] Wu | Non-Newtonian shear-thinning flows past a flat plate[END_REF] performed numerical and experimental studies of flow around an inclined plate. The effects of Reynolds number, shear-thinning index and angles of attack on drag and lift coefficients of the flat plate were presented.

For yield stress fluids, the only studies carried out to our knowledge concern the configurations for which the plate is perpendicular or parallel to the flow. For the flow of viscoplastic fluids parallel to the plate, an analytical study was performed by Piau [START_REF] Piau | Viscoplastic boundary layer[END_REF] . He proposed new approaches for determining the kinematic field and stresses for the boundary layer along the plate for creeping flows. Piau et Debiane [START_REF] Piau | The adhesive or slippery flat plate viscoplastic boundary layer for a shear-thinning power-law viscosity[END_REF] analyzed the variation in the drag coefficient and thickness of the boundary layer around a flat plate by describing the liquid region with a Herschel-Bulkley viscoplastic model and the solid region with Hooke's constitutive equation. Boujlel et al. [START_REF] Boujlel | Boundary layer in pastes-Displacement of a long object through a yield stress fluid[END_REF] studied experimentally the drag force and velocity fields around a plate parallel to the flow for a wide range of velocities. Balmforth et al. [START_REF] Balmforth | Viscoplastic boundary layers[END_REF] reexamined Oldroyd's viscoplastic boundary layer analysis as previously performed by Piau [START_REF] Piau | Viscoplastic boundary layer[END_REF] .

They proposed a general solution with an alternative theory that allows the boundary layer to be curved and of finite length. Contrary to Piau 7 , they found static and mobile rigid zones outside the boundary layer. In addition, Balmforth et al. [START_REF] Balmforth | Viscoplastic boundary layers[END_REF] indicated that it was not necessary to introduce elastic deformation below the yield stress to describe the displacement far from the plate, thus contradicting the studies of Boujlel et al. [START_REF] Boujlel | Boundary layer in pastes-Displacement of a long object through a yield stress fluid[END_REF] . For the same configuration, Ahonguio et al. [START_REF] Ahonguio | Flow of an elasto-viscoplastic fluid around a flat plate: Experimental and numerical data[END_REF] recently carried out a numerical and experimental study with an elastoviscoplastic model in which the variation in the plastic drag coefficient according to the Oldroyd numbers and the shear stress on the wall for various elastic numbers is determined. It should also be noted that the study of the effect of elasticity on the velocity profiles and the rigid zones around the plate in the flow has recently been studied numerically by Ferreira et al. [START_REF] Ferreira | External flows of elasto-viscoplastic materials over a blade[END_REF] . For flows perpendicular to plates and knife-edges, Brookes and Whitmore [START_REF] Brookes | Drag forces in Bingham plastics[END_REF][START_REF] Brookes | The Static Drag on Bodies in Bingham Plastics[END_REF] studied experimentally the static drag force of these submerged bodies in China clay suspensions.

However, the results are questionable because the rheometrical characterization of the clay suspensions is fragmentary and the slip on the walls of the obstacles is not controlled. Savreux et al. [START_REF] Savreux | Flow normal to a flat plate of a viscoplastic fluid with inertia effects[END_REF] used a Bingham model for the numerical study of the flow structure with a negligible inertia as a function of the Oldroyd number. They proposed empirical solutions for the drag coefficient for Reynolds numbers ranging from 0.01 to 30. In two recent articles, Patel and Chhabra [START_REF] Patel | Steady flow of Bingham plastic fluids past an elliptical cylinder[END_REF][START_REF] Patel | Effect of the angle of incidence on laminar forced convection from an elliptical cylinder in Bingham plastic fluids[END_REF] carried out a numerical study of the flow of a Bingham fluid around an elliptical cylinder. In a first study [START_REF] Patel | Steady flow of Bingham plastic fluids past an elliptical cylinder[END_REF] , they examined the influence of the aspect ratio. When the aspect ratio of the elliptical cylinder tends to low values, the case becomes similar to that of flow perpendicular to a plate. When the aspect ratio is high, the case becomes similar to that of a plate placed parallel to the flow. In a second work [START_REF] Patel | Effect of the angle of incidence on laminar forced convection from an elliptical cylinder in Bingham plastic fluids[END_REF] , they studied the flow around a cylinder inclined with respect to the flow for an aspect ratio equal to 0.5. These studies will be included in the comparisons analyzed below.

Studies related to the plasticity theory for applications in soil mechanics can also be found for plates perpendicular or parallel to the ultimate pull-out force. They are relatively numerous given the importance of applications. For example, Merifield et al. [START_REF] Merifield | Stability of plate anchors in undrained clay[END_REF] carried out a numerical study of the influence of the shape of the anchors on the holding capacity in clay over a range of depths and for varying widths and lengths. They calculated the values of the ultimate pullout forces at soil failure limits on the plate-shaped anchors. Bemben et al. [START_REF] Bemben | The vertical holding capacity of marine anchors in sand and clay subjected to static and cyclic loading[END_REF] carried out an experimental study of the effect of various loading conditions on the vertical holding capacity of marine anchors in sand and clay subjected to static and cyclic loading. They determined the values of the load holding capacity and the static breakout factor as a function of anchor embedment depth.

Unlike the previous works cited, this study will provide new knowledge of the inclination effects of the plate on drag and lift forces generated in the flow of a yield stress fluid and a Newtonian fluid. The plate will be completely immersed in the fluid. First, the problem and governing equations are described. The experimental approach will then be detailed with the rheological characterization of the fluids and experimental set-up. Forces acting on the plate will be measured in the Carbopol gel, which is considered to be a yield stress fluid model, and in a Newtonian glucose syrup. The forces generated on the plate will be determined as a function of the inclination angle for high Oldroyd numbers. The numerical study will then be presented. The Herschel-Bulkley model, regularized with the Papanastasiou approach, will be used. The effects of the shear-thinning index n, and the yield stress Od on the drag and lift coefficients, will be then determined. In the final part of this article, experimental and numerical results obtained will be compared with those of the literature. The effects of elasticity and plasticity will be discussed.

Theory

The two-dimensional steady-state flow of a Newtonian fluid and a viscoplastic fluid around a plate inclined at various angles is considered. Figure 1 shows the geometry of the problem. The plate of width L and thickness e is immersed in the fluid. The ratio e/L of the plate is equal to 0.1. The fluid velocity U is uniform at the inlet. Dimensions of the domain are provided in Figure 1. The flat plate is at a symmetrical distance H between the upper wall (Wall 1) and the lower wall (Wall 2). The Wall 1 corresponds to the free surface and the wall 2 to the bottom of the flume. The minimum gap equal to H = 3 L was chosen between the bottom of the flume (Wall 2) and the plate as in the experimental condition. A distance of H = 3L was taken between the plate and the upper wall 1 in the numerical simulation.

Experimentally, this distance varies between 2L and 3L with the angle of inclination. The overall distance of the domain between the upper wall 1 and the lower wall 2 is 6L.

τ = Kγ ( ) + γ if τ > γ = 0 if τ ≤ τ
where K represents the consistency coefficient, τ the yield stress, n the shear-thinning index, γ the tensor of the strain rates, γ the magnitude of the shear rate, τ the stress tensor, and τ the second invariant of the stress tensor. The yield surface is defined by the von Mises criterion. Ovarlez et al. [START_REF] Ovarlez | Three-dimensional jamming and flows of soft glassy materials[END_REF] and Martinie et al. [START_REF] Martinie | Extending yield-stress fluid paradigms[END_REF] have demonstrated that the yield stress that will be used in the experiments complies with the von Mises criterion.

The boundary conditions may be written as follows:

• At the inlet and on the wall 2: U = U and U = 0

• On the wall 1: τ = 0

• On the plate: condition of no-slip: U = 0 and U = 0

• At the outlet, an outflow condition is defined as ∀(i ,j) = 0

The following dimensionless parameters can be defined:

The shear-thinning index: n

The Oldroyd number (yield stress effect/viscous effect):

!" = # $ %(&/() ) (3)
The Reynolds number (inertia effect/viscous effect):

*+ = ,( ) & -.) % (4)
The dimensionless minimum gap:

/ = 0 ( (5)
The viscous and plastic drag and lift coefficients are defined by the following equations

1" 234 = -5 6 .( ).8 .& .) % , 1" * = 1" 234 !" (6)
(1)

1; 234 = -5 < .( ).8 .& .) % , 1; * = 1; 234 !" (7)
F D and F P represent respectively the drag and lift forces per unit length. These are the forces exerted by the fluid on the obstacle. In the Cartesian coordinate system used in Figure1, the drag force is negative, and its value will be given as an absolute value. The lift force is positive.

Experimental study

Rheological characterization

A glucose syrup produced by Nigay S.A. with a density of 1500 kg/m 3 was used as a Newtonian fluid. The glucose viscosity was measured with shear rates ranging from 0.03 to 0.1s -1 in a cone-plane geometry using a DHR3 rheometer (TA Instruments). An antievaporation cell around the sample was used. The viscosity was measured at a temperature ranging from 18 to 25°C +/-0.5°C to evaluate accurately the value of the viscosity in the flume experiment. For example, the viscosity of glucose at 20.7°C is about 2186 Pa.s. The Carbopol gel used in this study is often considered as a yield stress fluid model in fluid mechanics experiments in particular thanks to its absence of significant thixotropy (Piau [START_REF] Piau | Carbopol gels: Elastoviscoplastic and slippery glasses made of individual swollen sponges[END_REF] , Dimitriou et al. [START_REF] Dimitriou | Describing and prescribing the constitutive response of yield stress fluids using large amplitude oscillatory shear stress (LAOStress)[END_REF] ) and its transparency. Here, it was obtained by diluting a mass quantity of 0.85 wt% of Carbopol 940 manufactured by Goodrich [START_REF] Goodrich | Product Literature[END_REF] in distilled water. The translucent acidic solution obtained was then neutralized by the addition of sodium hydroxide solution to obtain a transparent gel. Rheometrical characterization was analyzed under the steady state and oscillatory conditions. In steady state, the rheological properties of the gel were determined at rotational speed, temperature and evaporation controlled by simple shear measurements using the DHR3 rheometer. Samples were taken directly from the flume at the free surface and at the bottom of the flume. The measurements were taken with a cone and plane cell for shear rates ranging from 10 -3 to 100 s -1 . The surfaces of the cell were covered with rough paper to prevent slippage, as shown in the study by Magnin and Piau [START_REF] Magnin | Cone-and-plate rheometry of yield stress fluids. Study of an aqueous gel[END_REF][START_REF] Magnin | Shear rheometry of fluids with a yield stress[END_REF] and Meeker et al. [START_REF] Meeker | Slip and Flow in Soft Particle Pastes[END_REF] . The cone used had a diameter of 50 mm with an angle of 0.035 radians.

Parameter measurement uncertainties were estimated to be 7% for the yield stress fluid and 15% for glucose syrup.

Table 1 presents the rheological properties of the Carbopol gel. The parameter values of the Herschel-Bulkley model (eq.( 1)) are obtained after fitting with rheometrical measurements. This table gives also the elastic modulus G' and the viscous modulus G'' in the linear domain below the yield stress of the material. The strain γ is the limit of the linear domain. The Carbopol gel behaves like an elasto-viscoplastic fluid as shown by Piau [START_REF] Piau | Carbopol gels: Elastoviscoplastic and slippery glasses made of individual swollen sponges[END_REF] , Dimitriou et al. [START_REF] Dimitriou | Describing and prescribing the constitutive response of yield stress fluids using large amplitude oscillatory shear stress (LAOStress)[END_REF] and Coussot et al. [START_REF] Coussot | Aging and solid or liquid behavior in pastes[END_REF] . Ahonguio et al. [START_REF] Ahonguio | Influence of surface properties on the flow of a yield stress fluid around spheres[END_REF] showed that these gels may have a non-zero first difference of normal stresses. Yarin et al. [START_REF] Yarin | Elongational behavior of gelled propellant simulants[END_REF] and Balmforth et al. [START_REF] Balmforth | Extensional dynamics of viscoplastic filaments: I. Long-wave approximation and the Rayleigh instability[END_REF][START_REF] Balmforth | Extensional dynamics of viscoplastic filaments: II. Drips and bridges[END_REF] studied the elongation behavior of these gels. The fluidization or solid-liquid transition of Carbopol gels has been the subject of recent studies by Ovarlez et al. [START_REF] Ovarlez | Three-dimensional jamming and flows of soft glassy materials[END_REF] and Martinie et al. [START_REF] Martinie | Extending yield-stress fluid paradigms[END_REF] .

τ (Pa) K (Pa.s n ) n γ(%) G'(Pa) G"(
Several authors studied the creep behavior under the yield stress (Lidon et al. [START_REF] Lidon | Power-law creep and residual stresses in a carbopol gel[END_REF] , Caton and Baravian 34 , Divoux et al. [START_REF] Divoux | Transient Shear Banding in a Simple Yield Stress Fluid[END_REF] , Divoux et al. [START_REF] Divoux | From stress-induced fluidization processes to Herschel-Bulkley behaviour in simple yield stress fluids[END_REF] ). From an experimental standpoint, Dimitriou et al. [START_REF] Dimitriou | Describing and prescribing the constitutive response of yield stress fluids using large amplitude oscillatory shear stress (LAOStress)[END_REF] introduced the isotropic kinematic hardening (IKH) model which describes the variation in yield stress during flow. The study of Fraggedakis et al. [START_REF] Fraggedakis | Yielding the yield-stress analysis: a study focused on the effects of elasticity on the settling of a single spherical particle in simple yield-stress fluids[END_REF] is interesting to read in terms of its analysis of the IKH concept and its application to numerical modeling.

Experimental set up

The experimental set-up is presented in Figure 2. The applied velocity range is varied between 0.001 and 6 mm/s with an accuracy of 4%. The flume is made of transparent Plexiglas and is 3m long, 22cm wide and 21cm high. The inclined plate is 14 cm +/-0.1mm long, 3 cm +/-0.02mm wide, and 0.3 cm +/-0.02mm thick.

To prevent slippage, the wall surface of the flume and the cylinder are coated with a sandpaper of average roughness Ra =120µm manufactured by Norton. Plate-to-wall distance is equal to G = 3. The minimum distance H between the plate and the bottom of the flume can be adjusted with an accuracy of 2%. The sensor for measuring drag force was manufactured by TE Connectivity and the sensor for measuring the lift force by Scaime. These sensors can measure forces ranging from 0 to 50N with an accuracy of about 0.05% of full scale for drag and 0.025% of full scale for lift. The parasitic drag force of the arm holding the plate was measured without the plate and for each velocity and gap value. The force due to the interaction between the plate tip and the vertical wall of the flume was also estimated for each velocity. These two parasitic forces were deducted from the measured value of the drag force with the plate in order to obtain only the drag force on the plate. It should be noted that the parasitic drag force of the plate holding arm without the plate is negligible and represents less than 2% of the overall value, and for lift, this value represents less than 3% of the overall value. The experimental uncertainties on the plastic drag and lift coefficients was 13% and 10% respectively. These uncertainties consist of 7% for yield stress, 6% for drag force and 3% for the lift force. Concerning the Oldroyd number, the uncertainty was of the order of 16%.

During the measurements from rest, a transient state where the forces increase with time is observed, followed by a steady state independent of time. The duration of the transient state depends on the imposed velocity and nature of the fluid. At a fixed velocity, the transients in our experiments can be reproduced relatively. No evidence was found of phenomena such as those observed during sedimentation of spheres in Carbopol gels [START_REF] Atapattu | Wall effect for spheres falling at small Reynolds number in a viscoplastic medium[END_REF][START_REF] Hariharaputhiran | The settling of spheres in a viscoplastic fluid[END_REF] . These authors [START_REF] Atapattu | Wall effect for spheres falling at small Reynolds number in a viscoplastic medium[END_REF][START_REF] Hariharaputhiran | The settling of spheres in a viscoplastic fluid[END_REF] pointed out that the reproducibility of the sphere terminal settling velocity is poor. Atapattu et al. [START_REF] Atapattu | Wall effect for spheres falling at small Reynolds number in a viscoplastic medium[END_REF] found that reproducible results could be obtained in Carbopol gels only after releasing 4-10 spheres. They found the settling velocity to be low for the first sphere. This velocity increases to an asymptotic value after a number of spheres were released. Hariharaputhiran et al. [START_REF] Hariharaputhiran | The settling of spheres in a viscoplastic fluid[END_REF] argue that these observations are consistent with the hypothesis of network damage caused by shear with subsequent healing. There are a number of differences between these flow conditions and ours. The present study is carried out with an imposed speed and not an imposed force as in settling. Our yield stress is higher by a factor of about 10. Indeed, our suspension is more concentrated in Carbopol microgel and a more compact structure is thus generated [START_REF] Piau | Carbopol gels: Elastoviscoplastic and slippery glasses made of individual swollen sponges[END_REF] . However, some tests have been done by changing the initial state of stress in the fluid. For example, a grid was passed through the fluid to homogenize the state of stress and the Carbopol gel structure. A notable difference in the transient of the first test was found compared to the transients obtained after the grid was passed. It is possible that this phenomenon is compatible with the observations made by the authors working on the settling of spheres in Carbopol but with less intensity. Similar phenomena of unusual transients and the effect of initial stress state were also observed on the very slow rise of bubbles in Carbopol gels [START_REF] Mougin | The significant influence of internal stresses on the dynamics of bubbles in a yield stress fluid[END_REF] . Only steady state forces have been studied in depth in the present study. It can be seen that lift and drag forces increase relatively linearly with the velocity for each angle. In Figure 3a, the drag force increases when the angle increases. In Figure 3b, the lift force is minimal at 0° then increases with the angle up to approximately 45°, then decreases to zero at 90°. The lift values are non-zero at 0° and 90° angles for high velocities. This may be due to the non-negligible thickness of the plate and to the uncertainty surrounding the angle of inclination.

Wu and Thompson 6 studied the flow of a non-Newtonian fluid around an inclined plane plate with a variable angle. They proposed analytical expressions for viscous drag Cd and lift Cl depending on the angle, defined by:

Cd = Cd (θ A B ) sin E θ + Cd (θ A ) cos E θ (8) Cl = E ICd (θ A B ) -Cd (θ A ) K sin 2θ (9)
They specified that these equations can also be applied to a shear-thinning fluid, but only in creeping flow. The drag and lift coefficients as a function of the angle can be calculated from the drag coefficient at 0° (Cd MNO(P A ) ) and 90° (Cd MNO(P A B ) ) angles. These equations will be used in section 6.

Yield stress fluid

Figure 4 shows the variation in drag force (Fig. 4a) and lift force (Fig. 4b) as a function of velocities for various inclination angles. The range of Reynolds number used is between 1.5 × 10 -9 and 1.7 × 10 -3 with Oldroyd numbers ranging from 5.5 to 179. Figure 4a shows that the variation in drag force as a function of velocity is similar, regardless of the inclination angle.

Unlike Newtonian fluids, the drag force tends towards a plateau value at low velocities corresponding to high yield stress effects. The drag force increases with velocity because the viscous forces increase with respect to the yield stress forces. The level of force increases with the inclination angle. Figure 4b shows that, just like the drag force, the lift force tends towards a plateau value when the velocity tends towards 0 and thus when the yield stress effects become high. The lift force increases with velocity when the velocity becomes higher. Figure 5 presents the measurements results for the dimensionless representation. It shows the variation in plastic drag and lift coefficient as a function of the Oldroyd number for various inclination angles. Regardless of the inclination angle, these plastic coefficients decrease and tend towards a plateau value when the Oldroyd number becomes higher. Within this limit, the coefficients Cd* and Cl* depend only on the yield stress. On the other hand, as Oldroyd numbers decrease and viscous friction becomes higher, the coefficient values increase.

In the studies carried out by Merkak et al. [START_REF] Merkak | Spheres and interactions between spheres moving at very low velocities in a yield stress fluid[END_REF] and Tokpavi et al. [START_REF] Tokpavi | Experimental study of the very slow flow of a yield stress fluid around a circular cylinder[END_REF] , the experimental data were interpolated by an empirical law defined by:

Cd * = Cd * Q + R S TU V with M = 1/(1+n ) (10)
A = Cd * Q represents the plastic contribution and B 1 the viscous part. Similarly, the equation of the plastic lift coefficient can be defined as follows:

Cl * = Cl * Q + R X TU V (11)
A E = Cl * Q is the plastic contribution of the lift and B 2 represents the viscous contribution.

The values Cd * Q and Cl * Q represent the plastic behavior of the material when Od tends towards infinity.

Figure 5 shows the variation in experimental values of the drag and lift coefficients according to the Oldroyd number and for various angles. These analytical expressions (10) and (11) are also shown in Figure 5. The values of the coefficients A 1 , A 10) and ( 11)

The values of the analytical expressions found in Table 2 are valid for an Oldroyd number between 5.5 and 179 and a negligible inertia. It was found that the A 1 values of the drag increase with the angle. The coefficient A 2 increases with the angle up to a maximum of 45° and then decreases.

On the basis of soil mechanics theory and, in particular, the calculation of anchors in a perfectly plastic soil (Raghavendra [START_REF] Raghavendra | Performance of Plate Anchors[END_REF] , Aubeny and Chi 41 ), the following expressions can also be written:

Cd * = Cd * (P A B ) sin E θ + Cd * (P A ) cos E θ ( 12
)
Cl * = E ICd * (P A B ) -Cd * (P A ) K sin 2θ (13)
Note that the maximum lift given by eq. ( 13) occurs at an angle of θ = 45°. In addition, in their studies on plate anchoring, Raghavendra [START_REF] Raghavendra | Performance of Plate Anchors[END_REF] and Anderson et al. [START_REF] Anderson | Deepwater Anchor Design PracticeVertically Loaded Drag Anchors[END_REF] proposed solutions for the drag Cd * (θ A ) and Cd * (θ A B ) which are applicable to a perfectly plastic material.

These expressions are defined by the following eq. ( 14) and ( 15):

Cd * (θ A ) = 4 + 30 \ ] (14)
Cd * (θ A B ) = 6π + 4 + 2 \ ] 2 + _ √E (15)
Figure 6 compares the experimental results with those of this theory for high yield stress effects. Equations ( 14) and (15) give in our case: Cd * (P A ) = 7 and Cd * (P A B ) = 23.81. Drag and lift were found to vary similarly with the theoretical results. In Figure 6a, the experimental values of the plastic drag and the theoretical prediction are generally close except at high angles. In Figure 6b, non-negligible differences were observed between the experimental values of lift and those of the plasticity theory of soil mechanics.

ϴ (°)

Expeirmental results

Theory Eq.( 13) (b)

Numerical modeling

A numerical study was performed to compare these predictions with the experimental results. For this modeling, the Herschel-Bulkley law (eq. ( 1)) was regularized by the wellknown Papanastasiou model [START_REF] Papanastasiou | Flows of Materials with Yield[END_REF] . This regularized model has already been used in many studies, particularly for flows around obstacles [START_REF] Savreux | Flow normal to a flat plate of a viscoplastic fluid with inertia effects[END_REF][START_REF] Tokpavi | Very slow flow of Bingham viscoplastic fluid around a circular cylinder[END_REF][START_REF] Mitsoulis | On creeping drag flow of a viscoplastic fluid past a circular cylinder: wall effects[END_REF] . It is defined by the following eq. ( 16) where m represents the regularization parameter.

# = d%e () 8) + # $ I8 +fg .h|e | K |e| j e (16)

Numerical method

Version 17.1 of the Ansys-Fluent software was used based on a finite volume approach with a double-precision solver and a second-order spatial discretization. The SIMPLEC (Semi-Implicit Method for Pressure Linked Equations-Consistent) method was adopted. An absolute convergence criterion of 10 -8 was chosen for the convergence of velocity and continuity equations. For all the numerical simulations, the Reynolds number is kept very low at 0.001.

Influence of the mesh and the shape of the plate

In order to optimize all the numerical parameters and to make the results independent of the numerical approach, a systematic study of the effect of the mesh (number of nodes), the shape of the plate (rectangular, oblong) and the type of mesh (quadrilaterals, triangles) on the plastic drag and lift coefficients was performed. For the mesh it was found that a minimal refinement of the order of 116 000 nodes is required to make the mesh independent of the drag and lift coefficient. This mesh is highly refined towards the four corners of the plate

The shape of the plate was also modified to verify whether the presence of corners, which may be singular points in the numerical modeling, had an effect on the value of the drag and lift coefficients. A plate with rounded edges was then created. The change in the geometric shape of the plate from a rectangle to that of a rounded shape does not modify the drag and the lift values to any great extent, with the maximum difference obtained being 5% for the 26.5° angle and 5.5% for the 45°angle. Singularities are present but they do not necessarily influence the values of force coefficients.

A triangular mesh was also made to verify that the type of mesh did not change the final result. Between the quadrilateral mesh and the triangular mesh a maximum difference of 3% was obtained.

The effect of varying the Reynolds number between 0.01 and 0.0001 on the drag and lift coefficients for the 45° and 26.5° angles was also studied. The values of these coefficients obtained in this range of Reynolds number are similar. Following this study, the refined quadrilateral mesh of 116 000 nodes with a plate of rectangular shape will be used to carry out the numerical computations. For the Papanastasiou 43 regularization (eq.16), the effect of the regularization parameter m on the force coefficients was also systematically studied and a value of m = 10 4 was chosen for all calculations. The boundary between yielded and unyielded zones is obtained by the condition τ = τ (1 ± ε) proposed by Burgos et al. [START_REF] Burgos | On the determination of yield surfaces in Herschel-Bulkley fluids[END_REF] . The value of ε = 0 was determined based on the study of Burgos et al [START_REF] Burgos | On the determination of yield surfaces in Herschel-Bulkley fluids[END_REF] .

Drag and lift coefficients

In this part, shear-thinning fluids and then Herschel-Bulkley yield stress fluids will be examined.

Shear-thinning fluid

Figure 7 shows the variation in viscous drag and lift coefficients as a function of the angle at negligible inertia and for various values of the shear-thinning index n at Od = 0.

In Figure 7a, the viscous drag coefficient increases with the angle from a minimum at 0° to a maximum at 90°. Figure 7b shows that, as n decreases, viscous lift Cl vis increases. It is also found that maximum lift measured at the 45° angle for a Newtonian fluid (n = 1, Od = 0) occurs at smaller angles when n decreases. From n ≤ 0.4, the maximum lift is found at 26.5°. Concerning the influence of the shear-thinning index n, the increase in n leads to an increase in the drag coefficient as shown in Figure 8a. The effect of the shear-thinning index on the drag, as shown in Figure 8b, is thus reduced when Od increases.

Yield stress fluid

The asymptotic values of theoretical plastic drag at the 0° and 90° angles given respectively by Piau and Debiane 8 and by Merifield et al. [START_REF] Merifield | Three-dimensional lower bound solutions for stability of plate anchors in clay[END_REF] are shown in Figure 8. For an infinitely thin plate in an infinite domain, and parallel to the flow, Piau and Debiane 8 proposed the following expression for plastic drag:

Cd * = 4 m1 + no E π p q / ( p ) TU r /( p ) s (17)
Figure 8 shows that the difference with our results is substantial (46%) especially for small Od values. This is due to the fact that Piau and Debiane's theory is based on an infinitely thin plate, which is not the case in this study (e/L=0.1). With a much smaller thickness (e/L=0.003), a difference of less than 7% was obtained. If the results of Patel and Chhabra [START_REF] Patel | Steady flow of Bingham plastic fluids past an elliptical cylinder[END_REF] for an elliptical cylinder of the greatest slenderness used (E = 10), and for a Bingham fluid, are applied, a value of Cd * = 6 + _.tu TU +

.tE TU .vwX is found. Formula (20) and table 11 of the Patel and Chhabra publication [START_REF] Patel | Steady flow of Bingham plastic fluids past an elliptical cylinder[END_REF] were used with our length scale. The flow for a high aspect ratio E tends towards a configuration corresponding to a plate placed parallel to the flow.

Note that the slenderness E used by Patel and Chhabra is the inverse of the ratio e/L used in our study so for E = 10 then e/L = 0.1. For Od= 5.5 a value of Cd * = 7.6 is obtained and for Od=179, Cd * = 6.11. These predictions are relatively close to the numerical solution and better than those of Piau and Debiane 8 .

For a high Od value (Od = 179), Patel and Chhabra's solution is satisfactory, among other things because the aspect ratio is close to that of the numerical simulation.

Looking at Patel and Chhabra's predictions' [START_REF] Patel | Steady flow of Bingham plastic fluids past an elliptical cylinder[END_REF] The plasticity theory that led to eqs. ( 12) -( 15) can also be applied. As can be seen in Figure 8b, a satisfactory approximation was obtained with these results.

Figures 8c and8d show the variation in plastic lift coefficient as a function of the inclination angle for Od = 5.5 and 179. The lift is zero for 0° and 90° angles and is maximum at a particular angle. In Figure 8c, the maximum lift is obtained at θ = 26.5° for n = 1. This maximum subsequently located between θ = 15° and θ = 26.5° for n = 0.5 and n = 0.4 when the Oldroyd number Od is 5.5. Thus, when the Oldroyd number Od increase to a high value (Od = 179) in Figure 8d, the maximum lift is obtained between θ = 15° and θ = 26.5 regardless of n. This inclination angle corresponding to the maximum lift is different to that obtained by the plasticity theory of eq. ( 13) (θ = 45°).

Our numerical results approximate the results of Patel and Chhabra 17 relating to the flow of a Bingham fluid around a cylinder inclined at an angle to the direction of flow. They find that the amplitude of the lift coefficient is maximun for an angle of about 30° (26.5° in our study).

There is a degree of imprecision because their results are calculated every 15°. As in our study, this maximum does not seem to be significantly influenced by the value of Od. The orders of magnitude of Cd * and Cl * cannot be compared because their aspect ratio is very different from that of our study.

When the yield stress is high (Od = 179), the curves no longer depend on n because the viscous effects are negligible. This variation as a function of n is different in the case of a shear-thinning fluid (Fig. 7b) and when the yield stress (Od = 5.5) is low (Fig. 8c). Thus, the yield stress required is inversely proportional to Cd * Q . The results from Figure 8 therefore show that when the inclination angle increases, a smaller yield stress is required to stabilize it in suspension.

Numerical and experimental comparison

Newtonian fluid

As previously seen, Wu and Thompson [START_REF] Wu | Non-Newtonian shear-thinning flows past a flat plate[END_REF] proposed laws (eqs. ( 8) and ( 9)) which can be used to calculate the force coefficients for a creeping flow. The drag coefficient values of the plate at 0° and 90° angles introduced in eqs. ( 8) and ( 9) to calculate Newtonian theory are derived from the numerical values for a gap G = 3.

Figure 9 shows a comparison of the numerical and experimental results as well as the Newtonian theory of viscous drag and lift coefficient as a function of angle for a Newtonian fluid. Figure 9a shows the increase in viscous drag when the angle increases. The maximum difference between the numerical and experimental results was 6%. By taking into account the maximum experimental uncertainties, estimated at 15%, the values are similar. Good agreement was found between Newtonian theory (eq. ( 8)) and the experimental results.

Figure 9b shows a similar change in lift between the experimental measurements, the numerical prediction and Newtonian theory of eq. ( 9). In particular, the maximum lift is obtained for θ = 45°. The maximum differences obtained are of the order of 17%. When the experimental uncertainties are taken into account, a good correlation is observed between experimental measurements, numerical predictions and Newtonian theory.

Yield stress fluid

The experimental drag and lift measurements will now be compared with the numerical simulation results using the Herschel-Bulkley model (Figure 10) for Oldroyd number values of Od = 5.5 and 179. The numerical values of the drag coefficients are generally slightly higher than experimental results, but follow the same variations. For the plastic lift (Figures 10c and10d), the numerical results are greater than the experimental results over the range 5° to 26.5°.

However, from the 45° angle, the experimental results are greater than the numerical results.

(d)

This angle of 45° corresponds to the maximum of the measured lift. Numerically it is smaller and depends on the yield stress value Od and the shear-thinning index n. As shown previously in the numerical results study, the effect of mesh type and plate shape on the drag and lift values is negligible. Any possible numerical singularity of the corners is absent.

We propose to discuss the differences between experimental results and anelastic viscoplastic modeling by referring to the recent publication of Fraggedakis et al. [START_REF] Fraggedakis | Yielding the yield-stress analysis: a study focused on the effects of elasticity on the settling of a single spherical particle in simple yield-stress fluids[END_REF] . These authors numerically modeled flows during sedimentation of a sphere in a yield stress fluid and compared their numerical predictions with experiments in Carbopol gels. The constitutive law describes the elastic effects and the hardening effects observed on Carbopol gels by Dimitriou et al. [START_REF] Dimitriou | Describing and prescribing the constitutive response of yield stress fluids using large amplitude oscillatory shear stress (LAOStress)[END_REF] . They found that the plastic drag Cd* values decreased when the elasticity of the material increased. This could explain the differences between our numerical and experimental results. However, in our study, the dimensionless Deborah number used to characterize the elasticity of a material is very small. This number can be defined by: De = ] .

• ' ′ /

. In this study, it varies between 1.274×10 -4 and 1.33×10 -4 . The elastic effects therefore seem negligible. On the other hand, this study focuses on high plastic effects with Oldroyd numbers ranging from 5.5 to 179. The phenomena of this study seem more accentuated by plasticity rather than by elasticity. This would be consistent with the analysis proposed by Fraggedakis et al. [START_REF] Fraggedakis | Yielding the yield-stress analysis: a study focused on the effects of elasticity on the settling of a single spherical particle in simple yield-stress fluids[END_REF] . In their study, they introduced the concept of isotropic kinematic hardening (IKH) proposed by Dimitriou et al. [START_REF] Dimitriou | Describing and prescribing the constitutive response of yield stress fluids using large amplitude oscillatory shear stress (LAOStress)[END_REF] . The idea of this model is that the yield stress surface of the material changes dynamically with the flow field. Fraggedakis et al. [START_REF] Fraggedakis | Yielding the yield-stress analysis: a study focused on the effects of elasticity on the settling of a single spherical particle in simple yield-stress fluids[END_REF] showed that this concept introduced in numerical modeling can be used to represent the experimental results obtained with Carbopol gels. To better understand the shift in the maximum lift, it would certainly be interesting to integrate in the numerical study models describing precisely the solid-liquid transition like the IKH model. Moreover, the influence of initial residual stresses in the gel cannot be excluded. Mougin et al. [START_REF] Mougin | The significant influence of internal stresses on the dynamics of bubbles in a yield stress fluid[END_REF] have shown the influence of these stresses in a study where they analyzed the rise of a bubble in Carbopol gels under quasi-static flow. In addition, Ahonguio et al. [START_REF] Ahonguio | Influence of surface properties on the flow of a yield stress fluid around spheres[END_REF] showed the variation in normal forces for these types of fluid as a function of the shear rate. The presence of these normal forces could in part explain the differences obtained.

Conclusion

The influence of plate inclination angle in the non-inertia flow of a yield fluid and a Newtonian fluid was studied experimentally and numerically. For yield stress fluids, the study focused on predominant plasticity effects.

From a numerical point of view, the variation in drag and lift coefficients as a function of the inclination angle, the shear-thinning index and the stress yield was established.

A good correlation for Newtonian fluids was found for the drag and lift coefficient as a function of angle between the experimental and numerical results as well as a good approximation with the prediction of the Wu and Thompson 6 model.

For yield stress fluids, the comparison between numerical prediction and experimental measurements shows that the variation in lift and drag as a function of angle is similar.

However, the maximum lift is predicted by the plasticity theory and experimental measurements at the 45° angle, whereas numerical predictions based on the Herschel-Bulkley model give a lower angle.

Analytical solutions of plastic drag and lift coefficients from the plasticity theory or boundary layer on a flat plate predict the experimental measurements relatively well. The numerical solutions give results that are relatively close to the experimental measurements except for the maximum lift angle.

The differences between experimental and numerical results do not seem to derive from elastic effects, but rather from phenomena related to plasticity. The change in yield stress or the isotropic kinematic hardening proposed by Dimitriou et al. [START_REF] Tamada | Low-Reynolds-number flow past a cylindrical body[END_REF] would certainly appear to play a role in these differences. In the future, modeling should take into account this more realistic solid-liquid transition for Carbopol gels.

Figure 1 :

 1 Figure 1: Schematic representation of the problem

Figure 2 :

 2 Figure 2: Experimental set-up

Figure 3

 3 Figure3shows the variation in the drag and lift force as a function of velocity with

Figure 3 :

 3 Figure 3: Variation in drag force (a) and lift (b) as a function of velocity for various inclination angles for a Newtonian fluid.

Figure 4 :

 4 Figure 4: Variation in drag force (a) and lift (b) as a function of velocity for various plate inclination angles for the yield stress fluid.

Figure 5 :

 5 Figure 5: Variation in experimental plastic drag (a) and lift coefficient (b) as a function of the Oldroyd number for various inclination angles. The analytical solutions (eqs. (10) and (11)) are also plotted in the figure.

Figure 6 :

 6 Figure 6: Comparative between the values of the plastic drag and lift coefficient at Od = 179 with n = 0.4 and the plasticity theory as a function of the inclination angle.

Figure 7 :

 7 Figure 7: Variation in viscous drag (a) and lift (b) as a function of the angle for a Reynolds number Re = 0.001, a gap G = 3 and for various values of n.As shear thinning increases i.e. when n decreases, the intensity of the viscous drag coefficient

Figure 8

 8 Figure8shows the influence of n on the plastic coefficients Cd* and Cl* as a function

  a value of Cd * = 20.87 + Et TU .vxX is obtained for E=0.1, with their formula (20) and their table 11. For Od= 5.5 and Od=179 values of Cd * = 31.25 and Cd * = 22.23 are obtained respectively. If these values are placed on Figure 8a, very good agreement can be seen with the computed numerical solutions of our study.

Figure 8 :

 8 Figure 8: Variation in plastic drag (a) and (b), lift coefficient (c), (d) as a function of the inclination angle for two Oldroyd numbers values Od = 5.5 and Od = 179 with n = 1 and n = 0.4 at gap G = 3.

  Additional calculations were made to determine the variation in plastic drag and plastic lift for a value of Od = 1000 and for a given angle (θ = 26.5°). Cd* values of 14.73 and 14.6 were obtained for Od=179 and Od=1000 respectively. The difference between these two values is 0.88%. The values of lift Cl* are of the order of 15.92 and 15.763 respectively for Od = 179 and Od = 1000. The difference between these two values is 0.9%. In short, the increase in yield stress effects for Od = 179 or more becomes negligible on the plastic drag and lift values since the values of Cd* and Cl* for Od = 179 and Od = 1000 are very close. This trend obtained on the values of Cd* and Cl* for these two Od values remains valid, regardless of the angle. The values of Cd* and Cl* for very high Od values correspond to the situation where the plate is motionless in the fluid. These constant values can be used to calculate the minimum yield stress to keep the plate static in the fluid. This aspect will be discussed later in this article. The influence of domain size was also studied up to G = 100 in order to check whether the wall effects are negligible with respect to the values of drag and lift at high Od values (Od = 179) and at lower Od values (Od = 5.5) for n = 1 and n = 0.4 and for θ = 26.5°.The maximum difference obtained on the values of Cd* and Cl* between the gap G = 1 and G = 100 was 3.5%. From G = 1 and Od = 5.5, the domain can be considered as infinite, which means without any interaction effect with the walls. Other complementary studies of the influence of the gap G at 90° have also confirmed that wall effects on the plastic drag are negligible from a gap G = 1. For example, a difference of 1% was found between plastic drag values for a gap G = 1 and a gap G = 3 for Od = 5.5.In addition, as previously shown, the values of Cd * Q and Cl * Q can be used to estimate the minimum value of the yield stress needed to keep the plate static in the fluid. In this case, the force created by the weight of the plate (mass (m) × constant gravity (g)) is balanced by the force exerted by the yield stress flow. The yield stress needed to maintain the motionless plate can be defined by:

Figure 9 :

 9 Figure 9: Variation in viscous drag (a) and lift (b) coefficients as a function of the inclination angle for a Newtonian fluid with negligible inertia and a gap G = 3.

Figure 10 :

 10 Figure 10: Variation in plastic drag (a) and (b), and plastic lift (c) and (d) coefficients as a function of angle for Oldroyd numbers Od = 5.5 and 179, with a gap G = 3.

Table 1 :

 1 Rheological properties of Carbopol gel

	Pa)

Table 2 :

 2 2 , B 1 and B 2 are summarized in Table 2 as a function of angle. The value of M is equal to 0.714 in this study. Values of plastic drag Cd* and lift Cl* coefficient from eqs. (

	30	(a)		90°Theory 60°Theory
				45°Theory
				26.5°Theory
				15°Theory
	20			0°Theory
	Cd *
	10
	0
	0	50	100 Od	150	200

Acknowledgments

The Laboratoire Rhéologie et Procédés is part of the LabEx Tec 21 (Investissements d'Avenir -grant agreement n°ANR-11-LABX-0030) and of the PolyNat Carnot Institute (Investissements d'Avenir -grant agreement n°ANR-11-CARN-030-01). The PhD Research project of Z. Ouattara has been made possible thanks to an Ivory Coast Scholarship.