Mathieu Bajodek

Frédéric Gouaisbaut

Alexandre Seuret

Design of a dynamical controller as a quasi-predictive control law for linear time-delay systems

come

Design of a dynamical controller as a quasi-predictive control law for linear time-delay systems

I. INTRODUCTION

Time-delays often occur during modeling of systems which process numerical or bio-physical information. Controlling such systems without taking into consideration the delay may induce degradation of performances and instability of the closed-loop. To bypass the problem, the so-called Smith predictor or dead time compensator was originally introduced in 1959. It consists in predicting the future of the state (see [START_REF] Smith | A controller to overcome dead time[END_REF]) to design after a lag time a control law as if there is no delay. Rewritten as a backstepping method on the transport equation in [START_REF] Krstic | Delay compensation for nonlinear, adaptative and PDE systems[END_REF], this nominal predictive control law is considered to be an appropriate solution (see [START_REF] Karafyllis | Predictor feedback for delay systems : Implementations and approximations[END_REF]) to stabilize the system for any given delay under the condition of controllability of the system without delay.

Nevertheless, as presented in [START_REF] Assche | Some problems arising in the implementation of distributed-delay control laws[END_REF], when the infinitedimensional part of this control law, modeled as an integral term, is realized by a dynamical system, the implementation can be unsafe. Indeed, as explained in [START_REF] Palmor | Time delay compensation : Smith predictor and its modifications, The control handbook[END_REF], if the state matrix of the plant is not Hurwitz, it leads to an unstable polezero compensation in an internal loop. For these reasons, modifications are often added to the original Smith predictor to extend the stability properties to a larger range of systems (see [START_REF] Padhan | Modified Smith predictor based cascade control of unstable time delay processes[END_REF][START_REF] Normey-Rico | Dead-time compensators : A survey[END_REF]) or to reinforce the robustness with respect to uncertainties (see [START_REF] Masui | Modified state predictive control of continuous-time systems with input delay[END_REF]). Then, as no guarantee of stability of the closed-loop system are proposed, the design of an implementable predictive control law is still an open problem.

Regarding the literature, many researchers have circumvented the problem by developing numerical approximations 1 Mathieu Bajodek is with LAAS-CNRS, Univ de Toulouse, 7 avenue du Colonel Roche, 31077 Toulouse, France mbajodek@laas.fr 2 Frédéric Gouaisbaut is with LAAS-CNRS, Univ de Toulouse, 7 avenue du Colonel Roche, 31077 Toulouse, France fgouaisbaut@laas.fr 3 Alexandre Seuret is with LAAS-CNRS, 7 avenue du Colonel Roche, 31077 Toulouse, France aseuret@laas.fr like rectangular or trapeze methods. Nevertheless, the closedloop reveals itself to be a neutral time-delay system which contains an infinite number of characteristic roots in a vertical line (see [START_REF] Manitius | Finite spectrum Assignment problem for systems with delays[END_REF][START_REF] Michiels | The effect of approximating distributed delay control laws on stability[END_REF]). Focusing on the case in which the line has a positive real part, the paper [START_REF] Mondié | Approximation of control laws with distributed delays: a necessary condition for stability[END_REF] pointed out a sufficient condition of instability of the system looped with the approximated control law. It means that, for a given stabilizing feedback gain, there exists a maximal allowable delay from which the implementation becomes unstable. These results have raised the question of the design of implementable controllers from the choice of the state feedback gain (e.g. feedback synthesis [START_REF] Barreau | Static state and output feedback synthesis for time-delay systems[END_REF]) to the use of additional states (e.g. low-pass filter [START_REF] Mondié | Finite Spectrum assignment of unstable timedelay systems with a safe implementation[END_REF]).

In this paper, we proposed a new dynamical controller based on the first coefficients on Legendre polynomials of the integral term contained in the predictive control law. Then, the controller as a finite-dimensional model can be computed numerically and the closed-loop system seen as a time-delay system. Analyzing it, stability properties and robustness with respect to model uncertainties can be estimated with techniques developed before (mainly [START_REF] Seuret | Hierarchy of LMI conditions for the stability analysis of time-delay systems[END_REF][START_REF] Seuret | Complete quadratic Lyapunov functionals using Bessel-Legendre inequality[END_REF][START_REF] Bajodek | Insight into stability analysis of time-delay systems using Legendre polynomials[END_REF] dealing with Bessel-Legendre inequality). Two examples are tested at the end to illustrate that this new dynamical controller can be stable for a larger range of delays and more robust than other implementation methods.

Notations : Throughout the paper, N denotes the set of integers, R the set of real numbers, R n the n dimensional Euclidian space and R n⇥m the set of n ⇥ m real matrices. For any square matrix A 2 R n⇥n , A T represents its transpose, tril(A) its lower triangular part and H(A) = A + A T . Moreover, diag(d 0 , . . . , d N ) is the diagonal matrix defined by its diagonal coefficients (d 0 , . . . , d N ) and the symbol ⌦ traduces a Kronecker product. Lastly, the set of square-integrable functions from ( h, 0) to R n is noted L 2 ( h, 0; R n ) and the notation u t (✓) stands for u(t + ✓), for all t 0 and all ✓ 2 ( h, 0).

II. DESIGN OF THE DYNAMICAL CONTROLLER

In this section, the methodology to design the finitedimensional dynamical controller is developed. The construction is based on the nominal predictive control law given by for instance Krstic in [START_REF] Krstic | Delay compensation for nonlinear, adaptative and PDE systems[END_REF] which is in a second step approximated using a polynomial based approach. It consists of the truncation of the Fourier-Legendre serie of the input and the computation of the required first quasi-Legendre coefficients through a finite dynamical system. Then, our new dynamical controller can be seen as a quasi-predictive control law and its state representation which is easy to compute is obtained.

A. Nominal predictive control law

Let us consider the following linear system (1) subject to an input delay,

ẋ(t) = Ax(t) + Bu(t h), (1)
where x(t) 2 R n , u(t) 2 R m , dynamical state matrix A 2 R n⇥n and input matrix B 2 R n⇥m . The delay h > 0 is assumed to be constant but possibly uncertain. The nominal control law (2) is a Volterra integral equation of the second kind.

u(t) = K 0 B B B @ e hA x(t) + Z 0 h e ✓A Bu t (✓)d✓ | {z } ◆(t) 1 C C C A . (2)
As presented in [START_REF] Krstic | Delay compensation for nonlinear, adaptative and PDE systems[END_REF], it ensures asymptotic stability of the origin x = 0 under the assumption A + BK Hurwitz with K 2 R m⇥n the control gain. However, as an infinitedimensional law, the integral part ◆(t) cannot be implemented directly.

In order to understand the induced problem, have a look at the closed-loop system in Laplace domain in Fig. 1. One inspects the transfer function G(s) in between U (s) and I(s) (respectively Laplace transforms of u(t) and ◆(t)) to obtain

G(s) = Z 0 h e ✓(sIn A) Bd✓. (3)
Considering its realization, one can use the most used Smith predictor form, G(s) = (sI A) 1 (I n e h(sIn A) )B, which turns out to be tricky to compute for unstable A because it introduces a pole-zero cancellation on each eigenvalue of A.

Focusing on numerical approximations [START_REF] Manitius | Finite spectrum Assignment problem for systems with delays[END_REF][START_REF] Michiels | The effect of approximating distributed delay control laws on stability[END_REF], taking

G(s) = 1 N P N 1 k=0 e k N h(sIn A) B
, other implementation issues appear. The controller needs a large amount of data and becomes unstable from a low maximal delay. But, in the following, A might be unstable and these limitations are exceeded.

Inspired by what has been proposed in the literature and motivated by the promising results of Legendre polynomials to give sufficient stability properties for time-delay systems as [START_REF] Seuret | Hierarchy of LMI conditions for the stability analysis of time-delay systems[END_REF] suggests, one decided to approximate the term ◆(t) by considering a truncated Fourier-Legendre decomposition of the signal u t (i.e. working on its first Legendre polynomial coefficients).

B. Approximated control law on Legendre polynomials

Define the following approximated finite-dimensional control law

u(t) = K e hA x(t) + ◆ N (t) , (4)
where the approximated integral term ◆ N (t), defined by Equation ( 5), approximates ◆(t) on Legendre polynomials as ruled on Theorem 1,

◆ N (t) = C N U N (t). (5)
From one side, the vector

U N (t) = ⇥ u T 0 (t) . . . u T N (t) ⇤ T collocates each Legendre coefficients from 0 to N . For each k 2 N, these coefficients u k (t) 2 R m , projection of the transported signal u t 2 L 2 ( h, 0; R m ) on the k-th Legendre polynomial L k are recalled below u k (t) = Z 0 h L k (✓)u t (✓)d✓,
where,

L k : 8 > < > : [ h, 0] ! R ✓ ! ( 1) k k X l=0 ( 1) l ( k l )( k+l l )( ✓ + h h ) l ,
taking the notation ( k l ) for the binomial coefficient. From the other side, matrix C N gathers each associated weighting coefficients and is defined by

8 > > > > < > > > > : k = Z 0 h e ✓A L k (✓)d✓, 8k 2 N, I N = diag(1, . . . , 2N + 1), C N = ⇥ 0 . . . N ⇤ ( 1 h I N ⌦ B).
Remark 1: For computational issues, in the case A is nonsingular, the matrices sequence ( k ) k2N can be calculated recursively using Legendre polynomials properties of [START_REF] Gautschi | Orthogonal polynomials, quadrature, and approximation: Computational methods and software[END_REF].

8 < :

0 = A 1 (I n e hA ), 1 = A 1 (I n + e hA ) 2 h 0 , k+1 = k 1 2(2k+1) h A 1 k . Indeed, a recursion relation on Legendre polynomials gives for all k 2, L 0 k+1 (✓) L 0 k 1 (✓) = 2(2k+1) h L k (✓).
Remark 2: One can also note that if A is nilpotent at order N ⇤ , then for each k > N ⇤ , k = 0 n and the approximation ◆ N at order N = N ⇤ gives the exact value of ◆(t). It means that N ⇤ +1 Legendre coefficients are enough to compute the integral term.

The following theorem ensures that the simple convergence from ◆ N (t) to ◆(t) holds.

Theorem 1: The approximated integral ◆ N (t) converges to the expected integral ◆(t).

Proof: From the definition of C N , the error ◆N (t) between ◆(t) and ◆ N (t) is equal to

◆N (t) = Z 0 h e ✓A B u t (✓) N X k=0 2k + 1 h L k (✓)u k (t)d✓ ! .
Hence, using Cauchy-Schwartz inequality, we have

|◆ N (t)| 2 = | Z 0 h e ✓A Be N (t, ✓)d✓| 2 ,  h Z 0 h e T N (t, ✓)B T e ✓(A T +A) Be N (t, ✓)d✓,
where e N (t, .) is the Fourier-Legendre error of the signal u t on the orthogonal basis

(L k ) k2N of L 2 ( h, 0; R m ) defined as e N (t, ✓) = u t (✓) N X k=0 2k + 1 h L k (✓)u k (t).
Furthermore, there exists ↵ > 0 such that B T e ✓(A T +A) B  ↵I m , for all ✓ 2 ( h, 0). Then,

|◆ N (t)| 2  h↵ Z 0 h e T N (t, ✓)e N (t, ✓)d✓. Assuming u t 2 L 2 ( h, 0; R m ), the L 2 -norm convergence of
Fourier-Legendre error e N (t, .) concludes the proof. This theorem justifies that ◆ N (t) is a well-chosen term to approximate ◆(t). Therefore, u(t) defined by Equation ( 4) remains close to the nominal one described by Equation ( 2). Otherwise, this result just carries the problem of implementation over to realize U N . Indeed, besides having an accurate approximant of the integral term, the own dynamics of Legendre coefficients U N (t) described by Equation ( 6) cannot be implemented safely. Indeed, by derivation, it represents N integrators in cascade and makes appear u(t h).

( UN (t) = h 1 A N U N (t) + B N u(t) B ⇤ N e N (t, h), e N (t, h) = u(t h) h 1 C ⇤ N U N (t), (6)
with 8 > > > > < > > > > :

1 N = ⇥ 1 . . . 1 ⇤ T , 1 ⇤ N = ⇥ ( 1) 0 . . . ( 1) N ⇤ T , L N = tril(1 N 1 T N 1 ⇤ N 1 ⇤T N ), A N = (L N + 1 ⇤ N 1 ⇤T N )I N ⌦ I m , B N = 1 N ⌦I m , B ⇤ N = 1 ⇤ N ⌦I m , C ⇤ N = (1 ⇤T N I N )⌦I m .
The proof of this dynamical equation can be found in [START_REF] Bajodek | Insight into stability analysis of time-delay systems using Legendre polynomials[END_REF] and relies on integration by parts and Legendre properties at the boundaries.

Then, a last modification is needed to construct a safe implementable dynamical controller.

C. Quasi-predictive control law and its state representation

This last step is made up of the construction of a stable model close to system [START_REF] Padhan | Modified Smith predictor based cascade control of unstable time delay processes[END_REF]. In previous work [START_REF] Bajodek | Insight into stability analysis of time-delay systems using Legendre polynomials[END_REF], putting aside the error e N (t, h) done at order N on u(t h), it is proven that the finite-dimensional model in between u(t) and u(t h) obtained is a Pade approximant of the single delay h (i.e. the function e hs in Laplace domain). Moreover, having a look to numerical approximation of time-delay systems on Legendre polynomials we found that such techniques are also used to compute characteristic roots as accurate as desired (see [START_REF] Vyasarayani | Spectral approximations for characteristic roots of delay differential equations[END_REF]). Following this idea, removing this model error, the dynamical controller of state X c 2 R m(N +1) described by Equation ( 7) is finally proposed.

(

Ẋc (t) = h 1 A N X c (t) + B N u(t), u(t) = K e hA x(t) + C N X c (t) , (7)
Instead of realizing the transfer function G(s) or approximating it by delayed samples [START_REF] Manitius | Finite spectrum Assignment problem for systems with delays[END_REF][START_REF] Michiels | The effect of approximating distributed delay control laws on stability[END_REF][START_REF] Mondié | Approximation of control laws with distributed delays: a necessary condition for stability[END_REF], we decided to work on a transfer function G N (s) associated to the state representation

✓ h 1 A N B N C N 0 ◆ .
Remark 3: One can remark that, if h ! 0, then u(t) ! Kx(t). Indeed, the state X c (t) converges quickly to zero and, by Taylor expansions, we have

C N ! h!0 ⇥ B 0 . . . 0 ⇤ .
To sum up, the input u(t) is given by the nominal control law [START_REF] Krstic | Delay compensation for nonlinear, adaptative and PDE systems[END_REF] where the integral term ◆(t) is approximated by Equation ( 5) and where Legendre coefficients U N (t) are replaced by quasi-Legendre coefficients X c (t) which follows a dynamics chosen arbitrary close to the one of U N (t). Thus, this new finite-dimensional dynamical controller, constructed smartly, is then implementable and can be seen as a quasipredictive control law. It is a finite-dimensional model with an input x(t) and an output u(t) and of state representation

✓ h 1 A N + B N KC N B N Ke hA KC N Ke hA ◆ = ✓ A c B c C c D c ◆ .
In the next paragraph, the stability and robustness of the closed-loop system with our finite-dimensional dynamical controller is investigated.

III. STABILITY AND ROBUSTNESS OF THE CLOSED-LOOP SYSTEM

In this section, one focuses on stability and robustness properties of our closed-loop system and compares it with other time-delay compensators.

A. Closed-loop system

The design of a finite-dimensional dynamical controller leads to a closed-loop system which can be modeled as a time-delay system.

Indeed, the closed-loop system verifies Equation [START_REF] Masui | Modified state predictive control of continuous-time systems with input delay[END_REF].

⇠(t) = A⇠(t) + B d C d ⇠(t h), (8)
where the state

⇠(t) =  x(t) X c (t)
belongs to R with = n + m(N + 1) and matrices A, B d , C d defined for given matrices A, B, K and constant delay h by

A =  A 0 B c A c , B d =  B 0 , C d = ⇥ D c C c ⇤ .
The characteristic roots of sI A + A d e hs can directly be estimated through an algorithm based on Tchebychev discretization developed and explained in [START_REF] Wu | Reliably computing all characteristic roots of delay differential equations in a given right half plane using a spectral method[END_REF]. A spectral analysis is possible giving an idea of the performances and the stability with an error margin made by discretization. To go further into details, thanks to a Lyapunov-Krasovskii functional, a fine enough sufficient condition of stability and robustness with respect to the delay can be given.

B. Sufficient condition of stability

With similar calculation as presented in [START_REF] Seuret | Hierarchy of LMI conditions for the stability analysis of time-delay systems[END_REF][START_REF] Seuret | Complete quadratic Lyapunov functionals using Bessel-Legendre inequality[END_REF], a rewriting of the hierarchical sufficient condition of stability can be given at an order M 2 N as it was already done in [START_REF] Bajodek | Insight into stability analysis of time-delay systems using Legendre polynomials[END_REF].

Theorem 2: If there exist symmetric positive definite matrices P 2 S +m(M +1) , S 2 S m , R 2 S m and a delay h

such that  H(E(h)P F ) + S + R E(h)P G G T P T E(h) T S < 0, ( 9
)
where M is the order of the linear matrix inequality and

8 > > > > > > > < > > > > > > > : E(h) = diag(I , hI m(M +1) ), F =  A B d C ⇤ M B M C d A M , G =  B d C ⇤ M , H = ⇥ C d 1 T M I M ⌦I m ⇤ , S = H T SH, R =  hC T d RC d 0 0 hI M ⌦R ,
then system ( 8) is asymptotically stable, for the delay h.

Proof: Following the proof given in [START_REF] Bajodek | Insight into stability analysis of time-delay systems using Legendre polynomials[END_REF], consider the Lyapunov-Krasovskii functional candidate V given by [START_REF] Michiels | The effect of approximating distributed delay control laws on stability[END_REF].

V (t, ⇠) =  ⇠(t) ⇠ M (t) T P  ⇠(t) ⇠ M (t) + V S (t, ⇠) + V R (t, ⇠), ( 10
) with 8 > > > > > > > > > > > > > > < > > > > > > > > > > > > > > : ⇠(t) =  X(t) X c (t) , ⇠ M (t) = 2 6 4 R 0 h C d ⇠(t + ⌧ )L 0 (⌧ )d⌧ . . . R 0 h C d ⇠(t + ⌧ )L M (⌧ )d⌧ 3 7 5, V S (t, ⇠) = Z 0 h C d ⇠(t + ⌧ ) T S C d ⇠(t + ⌧ ) d⌧ ⇠ T M (t)( 1 h I M ⌦ S)⇠ M (t), V R (t, ⇠) = Z 0 h (h + ⌧ ) C d ⇠(t + ⌧ ) T R C d ⇠(t + ⌧ ) d⌧,
The positivity of V is ensured by the positive definiteness of P , S and R and application of Bessel lemma on V S . The derivative of V along the trajectories of system ( 8) is bounded from above by 2 4

⇠(t) 1 h ⇠ M (t) e M (t) 3 5 T  H(E(h)P F ) + S + R E(h)P G G T P T E(h) T S 2 4 ⇠(t) 1 h ⇠ M (t) e M (t) 3 5
,

with the error e M (t) = C d ⇠(t h) C ⇤ M ⇠ M (t).
Calculation details relied on application of Bessel lemma can be found in [START_REF] Bajodek | Insight into stability analysis of time-delay systems using Legendre polynomials[END_REF]. Therefore, if the linear matrix inequality (9) is satisfied, the origin of system [START_REF] Masui | Modified state predictive control of continuous-time systems with input delay[END_REF] with that delay h is asymptotically stable by application of the Lyapunov-Krasovskii theorem.

Remark 4: Notice that following [START_REF] Seuret | Hierarchy of LMI conditions for the stability analysis of time-delay systems[END_REF], the principle of hierarchy applied and if inequality ( 9) is satisfied for a delay h at order M ⇤ , then this inequality is also verified at order M > M ⇤ .

Remark 5: According to Remark 3, since A+BK is Hurwitz and because the linear matrix inequality is continuous in h, we expect to find an interval of delay which assert stability from 0 to a maximal bound.

Lastly, assuming that the closed-loop system is stable for a nominal delay h, one proposes a theorem which asserts the stability of the system (8) for some delays around.

C. Robustness with respect to delay uncertainties

Consider an unknown delay h belonging to the interval

[h h, h + h].
A criterion of robusteness is required to ensure stability of the whole system with a controller implemented for a delay h for all constant delays in the pocket [h h, h + h]. Thanks to the previous Lyapunov-Krasovskii approach, stability results can be extended to robustness. One defines the matrix (h) which appears in the linear matrix inequality (9) for a given h.

(h) =  H(E(h)P F ) + S + R E(h)P G G T P T E(h) T S It is clear that (h) is affine in h.
From this statement, it is possible to rewrite the theorem proposed in [START_REF] Seuret | Complete quadratic Lyapunov functionals using Bessel-Legendre inequality[END_REF]. The following revisited theorem gives a linear matrix inequality which ensures stability for all constant delays in the pocket [h , h + ].

Theorem 3: If there exist symmetric positive definite matrices P 2 S +m(M +1) , S 2 S m and R 2 S m and two positive scalars h and h + such that

(h ) < 0, (h + ) < 0, (11)
then system (8) is asymptotically stable, for all delays in the interval [h , h + ].

Proof: Consider the Lyapunov-Krasovskii functional candidate V given by [START_REF] Michiels | The effect of approximating distributed delay control laws on stability[END_REF]. The positivity of V is again verified by the positive definiteness of P , S and R. Then, from the linearity of (h) in h, the derivative of V for system (8) leads to

V (t, ⇠)  (h) = h h h + h (h + ) + h + h h + h (h )
Thus, the negativity of (h ) and (h + ) guarantees that the origin of system ( 8) is asymptotically stable for any constant delay h 2 [h , h + ].

Note that similar results could be found for robustness with respect to a time-varying delay which is an important issue (see [START_REF] Brech-Pietri | Commande robuste de systèmes à retard variable[END_REF]) or to model uncertainties on state matrix A or input matrix B.

IV. EXAMPLES

In this last section, the stability and robustness of our dynamical controller on two classical examples are evaluated. These simulations are conducted on purpose on unstable systems because if A was stable, it would be better to set up usual tools [START_REF] Normey-Rico | Dead-time compensators : A survey[END_REF]. From the examples, to get closer to the entire set of stability of our closed-loop system, a criterion of stability based on Bessel-Legendre as Theorem 2 is required and more precisely at an order M N because the design of our dynamical control is related to the N + 1 first Legendre coefficients. On Fig. 2, using Theorem 2 at order M = 5, the pointwise maximum allowable delay is given for each example and for a controller at order N 2 {0, . . . , 5}. It is compared with the delay from which classical rectangular or trapeze approximation methods fail to compute the nominal predictive control law (see the horizontal line). It is the one given by the stability of the nominal controller itself (i.e. det(I m KG(s)) has at least one zero with positive real part) which is a necessary condition of stability of the closed-loop according to [START_REF] Mondié | Approximation of control laws with distributed delays: a necessary condition for stability[END_REF].

In the sense that for stability, the impact of K turns out to be minor, we decided to fix K for desired properties.

A. First example

Consider a 1-dimensional system [START_REF] Barreau | Static state and output feedback synthesis for time-delay systems[END_REF], widely studied and known to be practically unstable by usual approximations of the nominal control law with K = 2 since h > 0.405 and to have no stabilizing feedback gain K since h > 0.693 as mentioned in [START_REF] Mirkin | On the approximation of distributed-delay control laws[END_REF].

ẋ(t) = x(t) + u(t h). (12)
First, taking K = 2, an estimation of the spectrum assignment of the closed-loop system with the dynamical controller at order N = 1 is computed for h = 0.4, h = 1 and h = 3 (see Fig. 3). The closed-loop system is stable for a larger bound of delays in comparison with approximation techniques. Indeed, thanks to Theorem 2 at order M = 5 and increasing pointwise the delay h, it appears that the system with our finite-dimensional dynamical controller is asymptotically stable for 0 < h < 1.3 from order N = 1 (see Fig. 2a).

Furthermore, for a delay equal to 0.4, 1 and 3 the eigenvalues have respectively positive real parts lower than 1, 0 and higher than 0. It means that there is a first interval of delay in which the performance given by A+BK are conserved and a second interval in which the closedloop stays stable but the performances are degraded. This is confirmed by simulation in time of our closed-loop system with the dynamical controller (4) with K = 2 for the corresponding delays. This is highlighted by Fig. 4 where the output x(t) is drawn in continuous line. The system with initial value x(0) = 1 is converging to the origin x = 0 for h 2 {0.4, 1} even if for h = 1 the exponential decay is reduced. Notice that, by rectangular approximation method of the integral, the output would diverge for h = 1.

Finally, Theorem 3 until order M = 5 gives a criterion of robustness with respect to delay uncertainties. For a controller at order N = 1 designed respectively for h = 0.4 and h = 1, the closed-loop system is proven to be stable on the pocket [h h, h + h] with h = 0.3 and h = 0.1. On Fig. 4, dotted lines represent the output of the controlled system with an error done on the delay randomly chosen on [ h, + h] . One can see that our closed-loop system at order N = 1 is robust with respect to delay uncertainties. Moreover, once the closed-loop system is stable for h, the performances are deteriorating while the error is increasing. We also noticed that the bound h is independent of the order N . Same statements could be done with respect to uncertainties on the state or input matrix.

B. Second example Consider now system (1) with

A = h 2 0 1 1 2 2 0 1 1 i , B = h 0 0 1
i , and with a gain K designed using a LQR (linear quadratic regulator) method as in [START_REF] Krstic | Delay compensation for nonlinear, adaptative and PDE systems[END_REF]. Notice that A is unstable.

The stability margin with respect to a constant delay h is computed with Theorem 2 for controllers from order N = 0 to N = 5 on Fig. 2b. The system remains stable until Fig. 5: Second example : Simulation in time with controller at order N = 1. h = 0.64 for N 1 which means a bound much larger than the one usually computed (h = 0.3). It is significant to note that N = 1 is enough to reach the maximum allowable delay because it also means that a smaller finite-dimensional set than classical approximation techniques (where hundreds of samples are taken into consideration) are required to design the controller. But, to balance these promising results, even if this dynamical controller is easy to compute and enlarges the interval of stability compared to the one reached by rectangular approximations, the stability of the closed-loop is still not guaranteed for all delays.

Finally, as for the previous example, Theorem 5 ensures some robustness with respect to uncertainties on the delay h. On Fig. 5, for a dynamical controller designed for h = 0.3 at order N = 1, the behavior of systems having a delay equal to 0.3 (continuous line) and 30 randomly chosen delays in 0.3 ± 0.1 (dotted lines) are drawn. The simulation in time confirms that the outputs x(t) is converging from ⇥ 1 0 0 ⇤ T to ⇥ 0 0 0 ⇤ T and that one more time we get the worst performances for the upper bound h + h.

V. CONCLUSIONS

A finite-dimensional dynamical controller has been designed through several steps based on an interpolation of the nominal control law by first Legendre polynomials. Techniques based on Legendre polynomials has been then proposed to ensure its stability with respect to the delay and its robustness as well for uncertain delays. As in many papers in the literature, one obtains also a closed loop system which is asymptotically stable up to a prescribed upper bound. Keeping the same feedback state representation, future work will focus on the choice of the control gain to optimize the closed-loop performances and robustness properties. In contrast to this late lumping approach, an early lumping approach could also be investigated.

Fig. 1 :

 1 Fig. 1: Block diagram of the closed-loop system with the nominal controller.

Fig. 2 :

 2 Fig. 2: Maximal allowable delay given by Theorem 2.

Fig. 3 :

 3 Fig. 3: Characteristic roots of the controlled system (N = 1).

Fig. 4 :

 4 Fig. 4: First example : Simulation in time with controller at order N = 1.


Léo C C Mesquita

Aymeric Vié
email: aymeric.vie@centralesupelec.fr

Laurent Zimmer

Sébastien Ducruix

Numerical analysis of flame shape bifurcation in a two-stage swirled liquid burner using Large Eddy Simulation

Keywords: Large Eddy Simulation, Staged combustion, Bifurcation, Combustion instability, Vortices

A flame shape bifurcation in the liquid-fueled two-stage swirled BIMER combustor is studied using Large Eddy Simulations. This combustor, developed at the EM2C Laboratory to study Lean Premixed Prevaporized (LPP) burners, is composed of a two-stage injection system: a central swirled pilot stage fueled with a pressure-swirl atomizer, to sustain a piloting flame, and an outer swirled stage fed with a multi-point injection, to generate the LPP regime. After ignition in the pilot-only operating condition, a V flame is

stabilized near the Inner Shear Layer (ISL). When switching to multipointonly injection, a flame shape transition is observed, and the flame bifurcates into a M-shape. In this work, we identify the mechanisms that lead to this bifurcation, and we show that the transition is driven by a complex coupling between the flame, the chamber acoustics and the ISL vortices. By switching to a multipoint-only injection, the fuel is essentially given to the ISL flame, which is mainly premixed. Because of the increased heat release and

Introduction

Lean Premixed Prevaporized (LPP) burners are prone to extinction, flashback or thermo-acoustic instabilities [START_REF] Ducruix | Combustion dynamics and instabilities: Elementary coupling and driving mechanisms[END_REF][START_REF]Combustion Instabilities in Gas Turbine Engines: Operational Experience, Fundamental Mechanisms, and Modeling[END_REF]. A candidate for the control of LPP burners is staged injection, which splits the fuel supply in several swirler stages. The two-stage BIMER combustor of EM2C laboratory has been developed to expand the knowledge on LPP staged burners.It is composed of an outer stage fed with a multipoint injection system, which goal is to generate the LPP regime thanks to fast atomization and evaporation; and of a central stage fueled with a pressure-swirl atomizer that sustains a rich diffusionlike pilot flame close to the injection system, aiming to stabilize the burner operation [START_REF] Providakis | Characterization of the coherent structures in swirling flames stabilized in a two-stage multi-injection burner: influence of the staging factor[END_REF][START_REF] Renaud | Flame shape transition in a swirl stabilised liquid fueled burner[END_REF].

Thanks to the versatility of this staged injection system, thermo-acoustic instabilities and flame shape bifurcations can be triggered and studied depending on the injection regime and a hysteretic behaviour can lead to multistable operating points. Indeed, experimental and numerical studies have exhibited three flame archetypes and three bifurcations between these flames when varying the fuel distribution between the stages [START_REF] Providakis | Characterization of the coherent structures in swirling flames stabilized in a two-stage multi-injection burner: influence of the staging factor[END_REF][START_REF] Renaud | Flame shape transition in a swirl stabilised liquid fueled burner[END_REF][START_REF] Renaud | High-speed diagnostics for the study of flame stabilization and transient behaviour in a swirled burner with variable liquid-fuel distribution[END_REF][START_REF] Renaud | Bistable behaviour and thermoacoustic instability triggering in a gas turbine model combustor[END_REF][START_REF] Cheneau | Characterization of the hysteresis cycle in a two-stage liquid-fueled swirled burner through numerical simulation[END_REF][START_REF] Cheneau | Numerical study of flame shapes and structures in a two-stage two-injection aeronautical burner with variable fuel staging using eulerian large eddy simulations[END_REF]: a V-flame, an M-flame and a Tulip-flame. In [START_REF] Renaud | Flame shape transition in a swirl stabilised liquid fueled burner[END_REF][START_REF] Renaud | Bistable behaviour and thermoacoustic instability triggering in a gas turbine model combustor[END_REF], the authors studied two of these flame shape bifurcations. Starting from a V-flame with dual injection and ceasing to feed the pilot injector with fuel leads to an M lifted flame. Reinjecting fuel from the pilot burner recovers the V-flame by flashback. In [START_REF] Cheneau | Characterization of the hysteresis cycle in a two-stage liquid-fueled swirled burner through numerical simulation[END_REF], this transition has been numerically evidenced for the operating conditions of [START_REF] Providakis | Characterization of the coherent structures in swirling flames stabilized in a two-stage multi-injection burner: influence of the staging factor[END_REF] and some mechanisms have been suggested to explain it. Here, the objective of this study is to investigate further the mechanisms driving this type of transition using an updated LES numerical setup with a finer mesh and a Lagrangian description for the liquid fuel to better reproduce the phenomena observed experimentally in [START_REF] Renaud | High-speed diagnostics for the study of flame stabilization and transient behaviour in a swirled burner with variable liquid-fuel distribution[END_REF]. Compared to [START_REF] Cheneau | Characterization of the hysteresis cycle in a two-stage liquid-fueled swirled burner through numerical simulation[END_REF], we focus on the operating conditions of [START_REF] Renaud | Flame shape transition in a swirl stabilised liquid fueled burner[END_REF][START_REF] Renaud | High-speed diagnostics for the study of flame stabilization and transient behaviour in a swirled burner with variable liquid-fuel distribution[END_REF][START_REF] Renaud | Bistable behaviour and thermoacoustic instability triggering in a gas turbine model combustor[END_REF], for which a similar transition is observed and a larger set of validation data is available. This paper is organized as follows. First, the combustor and the operating conditions are presented. The numerical setup is showed and a brief description of the flame ignition and initialization is given. Finally, the transition results are presented, and the main events and mechanisms controlling the transition are discussed.

Description of the setup

Experimental setup and operating conditions

BIMER is a labscale combustor designed to investigate the operation of a swirling two-stage injector. An inlet feeding line leads the air to a cylindrical plenum, which passes through the injection system and then reaches a rectangular combustion chamber (500 × 150 × 150mm 3 ). The chamber contains three water-cooled walls (entrance, top and bottom) and two large lateral silica windows for optical access. The swirling injector is composed of two stages (Fig. 1): the central stage is called the pilot stage. Including a pressurized nozzle it creates a hollow cone of droplets with a spray half angle of 30 • . It naturally admits approximately 15 % of the total air mass flow rate, and a geometric swirl number of 0.6 is imposed. The main external stage is called the multipoint stage, as fuel is injected through 10 holes of 0.3 mm diameter, placed at the second swirler vanes exits to perform a jet in cross-flow injection and improve mixing. This stage admits 85 % of the air mass flow rate, imposing a geometric swirl number of 1.

The operating conditions match the experimental data of Renaud [START_REF] Renaud | High-speed diagnostics for the study of flame stabilization and transient behaviour in a swirled burner with variable liquid-fuel distribution[END_REF].

The burner operates at atmospheric pressure, with an air mass flow rate of 43 g.s -1 preheated at 433K, and a total fuel mass flow rate of 1.64 g.s -1 , leading to a global equivalence ratio of 0.6, and a thermal power close to 73 kW. These parameters are kept constant, while the fuel distribution between the stages varies. This variation is measured by the staging factor α, which is the ratio between the fuel mass flow rate injected through the pilot stage and the total fuel mass flow rate. The staging factor is varied from 15% to 0%.

Numerical setup

AVBP, an unstructured code co-developed by CERFACS and IFPEN [START_REF] Gicquel | Large eddy simulations of gaseous flames in gas turbine combustion chambers[END_REF] that solves the 3D compressible reacting multi-species Navier-Stokes equations, is used. The LES sub-grid scales are modeled using the WALE model [START_REF] Nicoud | Subgrid-scale stress modelling based on the square of the velocity gradient tensor[END_REF]. A 3 rd order in time and space Two-step Taylor-Galerkin (TTGC) scheme is used. Artificial diffusion is implemented to treat numerical perturbations and stiff gradients [START_REF] Colin | Development of high-order taylor-galerkin schemes for les[END_REF]. Gaseous inlet and outlet are treated with the Navier-Stokes Characteristic Boundary Conditions (NSCBC) [START_REF] Poinsot | Boundary conditions for direct simulations of compressible viscous flows[END_REF]. Walls are considered as non-slipping and adiabatic, except for the two water-cooled walls which are treated as isothermal.

A Lagrangian modeling is used for the liquid phase, using one numerical parcel per physical droplet. Droplet diameter distributions have been extracted from experimental studies [START_REF] Renaud | High-speed diagnostics for the study of flame stabilization and transient behaviour in a swirled burner with variable liquid-fuel distribution[END_REF] for both stages. Evaporation is mod-eled with the Abramzon-Sirignano model [START_REF] Abramzon | Droplet vaporization model for spray combustion calculations[END_REF]. Details of the implementation can be found in [START_REF] Hannebique | Large eddy simulation of reactive two-phase flow in an aeronautical multipoint burner[END_REF]. Interaction of droplets with the walls is set to slip. This aims at reproducing in a very simplistic way a film-like behavior at the divergent exit observed in the experiments [START_REF] Renaud | High-speed diagnostics for the study of flame stabilization and transient behaviour in a swirled burner with variable liquid-fuel distribution[END_REF]. This setup and the injection strategy for the pilot spray were validated in [START_REF] Mesquita | Proceedings of the ASME Turbo Expo 2018[END_REF]. The chemistry is described by the BFER reduced scheme for kerosene [START_REF] Franzelli | A two-step chemical scheme for kerosene-air premixed flames[END_REF], under the unity Lewis number assumption, and the flame-turbulence interaction is treated using the Thickened Flame model (TFLES) [START_REF] Poinsot | Theoretical and Numerical Combustion, elearning@CERFACS[END_REF], where flame wrinkling is calculated using Charlette's formalism [START_REF] Charlette | A power-law flame wrinkling model for LES of premixed turbulent combustion part II: dynamic formulation[END_REF]. An unstructured mesh composed of 132 million tetrahedra is used (elements size can be seen in Fig. 1), keeping the thickening factor between 3 and 7.

Initialization procedure

To reproduce the experimental procedure, the flame is ignited at a pilotonly injection condition (α = 100%), which generates a V-shaped flame.

From this state, the fuel staging is slowly modified, changing α at the rate of 1% per millisecond of simulated physical time, until α reaches the value of 15%. The resulting V-shaped flame at α = 15% is shown in Fig. 2a and presents two main features: (1) a rich diffusion flame downstream the pilot injector mostly burning in a diffusion regime, characterized by a strong interaction between spray and flame and associated with high temperatures (around 2300K); and (2) further downstream, the remaining of the flame burning mostly premixed close to the global equivalence ratio, leading to lower temperature, near 1600K, confirming the observations made in [START_REF] Cheneau | Characterization of the hysteresis cycle in a two-stage liquid-fueled swirled burner through numerical simulation[END_REF].

This premixed flame is stabilized on the Inner Shear Layer (ISL) between the fresh and burned gases and the V-shape of the flame is linked with the flow topology created by the Conical Vortex Breakdown, present in this condition. Thus, the Central Recirculation Zone (CRZ) is large and open, starting upstream the pilot injection and reaching the center of the chamber.

Analysis of the flame shape transition

From the stable V state, α is reduced from 15% to 0% at a rate of 1% per ms, meaning that the pilot injection is progressively decreased as the multipoint injection is increased. In the following, the transition is carefully analyzed and the phenomena at stake are identified.

Flame shape evolution

In Fig. 2, the flame shape evolution is presented from the initial V-state to the final M-state. Between Fig. 2a and 2b, the pilot injection has been stopped (at t = 15ms) and the portion of the flame anchored by the central spray extincts. In Fig. 2b, a large toroidal vortex is visible (identified by arrows), increasing in size in the middle of the flame, in the portion stabilized at the ISL. These large vortices, not observed during the initial V state, occur at the axial position of the CRZ bubble contraction (Fig. 2b). They are directed inwards, leading to a strong oscillation of the inner flame and CRZ diameters. This first vortex is convected downstream, restoring the flame to its original shape and showing a periodic oscillatory behavior. In the following period (Fig. 2c), a new vortex, subsequent to that in Fig. 2b, starts forming. This vortex pushes the inner flame and CRZ inwards even further (Fig. 2d). However, this time, they do not recover their original diameter after the vortex is convected downstream. In the following period, a third vortex arrives and becomes large enough to further push the inner flame inside the CRZ and contract its diameter down to zero, leading to the flame merging in the CRZ (Fig. 2e). The flame is then separated into two parts by the vortex and each one evolves independently: an internal flame, once stabilized over the ISL is blown towards the pilot injector by the CRZ negative velocity (Fig. 2e) and eventually extinguishes (Fig. 2f); and an external flame, anchored at the OSL, keeps this anchoring and crosses the CRZ (as a consequence of the transition), and stabilize its leading edge at the CRZ's contraction and to take the final M-shape (Fig. 2f).

The final M flame is shown in Fig. 2f. It is completely stabilized inside the combustion chamber, burning mostly in lean premixed regime, with a lower temperature near 1600K. After the transition, the flow topology returns to the original non-reactive one: instead of the large open recirculation zone, the CRZ has a round bubble-like shape, characteristic of the Bubble Vortex Breakdown. The CRZ is filled with fuel vapor and air coming from the multipoint stage, which mixes with air coming from the pilot stage to produce a mixture close to the global equivalence ratio. A pressure drop at the CRZ is also observed, like in non-reacting conditions. The flame is stabilized by the OSL and the leading point is localized at the bottle-neck region of the CRZ, where the fresh gases reenter the CRZ and quenches the flame, preventing it to propagate and burn the mixture inside the CRZ.

Combustion instability and inner flame evolution

The acoustic activity inside the combustion chamber (Fig. 3.a) shows that, as soon as the pilot injection starts being reduced (t = 0 ms), the pressure fluctuations increase, reaching the maximum amplitude at the flame transition and, then, decreases. This pressure oscillation is in phase with the global heat release rate, meaning that the Rayleigh criterion is fulfilled and that a thermo-acoustic instability occurs, locked on the longitudinal quarterwave mode of the chamber (with a characteristic frequency of around 290 Hz [START_REF] Renaud | High-speed diagnostics for the study of flame stabilization and transient behaviour in a swirled burner with variable liquid-fuel distribution[END_REF]). This increase in thermo-acoustic activity is due to the change in fuel distribution and the fact that the flame now features LPP characteristics. These points will be addressed in more detail in the next section.

Simultaneously with the increase in acoustic and heat release rate oscillations, a larger flame surface oscillation is visible as the pilot injection is stopped (Fig. 3.b). The maxima of the flame surface evolution coincide with the presence of the toroidal vortices shown in Fig. 2, linking these vortices, which are at the origin of the contraction of the inner flame and of the transition, with the heat release rate oscillations. Therefore, to better study the flame shape transition, the evolution in time of the inner flame, defined as the portion of the flame that is stabilized over the ISL (as seen in Fig. 2), is closely observed. To analyze the evolution of its diameter, the inner flame is approximated by a circle, which has been verified as a reasonable assumption because of its axisymmetry. The beginning of the transition is clearly identified by the fast contraction of the ellipse at t = 15 ms, synchronized with the CRZ diameter contraction, attesting that both evolve very similarly.

Inner flame and ISL vortices dynamics

It is now clear that this transition is driven by the inner flame evolution and the way it interacts with the acoustic mode and the inner shear layer.

However, the phenomena involved here seem different from the classical driving mechanisms usually generating heat release rate fluctuations in premixed swirling flames subjected to incident acoustic perturbations. Here, the typ-ical flame tip roll-up [START_REF] Candel | Dynamics of swirling flames[END_REF] with a rotation direction from the center to the outwards of the chamber is not observed. Alternatively, as shown in Fig. 2b, a vortex rotating inwards is observed at the center of the flame, and not at its tip. It is rather the interaction between the acoustic mode and the ISL helical vortices that drives the dynamics in heat release rate oscillations by wrinkling the flame. Such mechanisms were also observed and analyzed by [START_REF] Caux-Brisebois | Thermoacoustic velocity coupling in a swirl stabilized gas turbine model combustor[END_REF][START_REF] Steinberg | Flow-flame interactions causing acoustically coupled heat release fluctuations in a thermo-acoustically unstable gas turbine model combustor[END_REF][START_REF] Steinberg | Effects of flow structure dynamics on thermoacoustic instabilities in swirlstabilized combustion[END_REF] in premixed lifted methane flames.

In the BIMER combustor, when the staging factor is decreased from α=15% to 0%, 15% of the fuel mass flow rate is redirected towards the multipoint injection. This changes the spatial distribution of fuel in the chamber as visible in Fig. 2, transferring fuel from the central nozzle spray region to the ISL region. As the flame must adapt to this new distribution, the flame front moves from a region burning in a stable rich diffusion regime to another where it is lean, premixed and strongly interacting with the ISL. Figure 3.a shows the evolution of the heat release rate as a function of time during the transition. One clearly observes a change in amplitude of oscillations of the heat release rate and of the flame surface at t = 6 ms, after the change in staging is initiated. This is then followed by an increase in acoustic pressure oscillations. This increase in the (longitudinal) acoustic pressure oscillations leads to axial velocity oscillations that excite the helical vortices in the ISL, probably associated with Kelvin-Helmholtz hydrodynamic instabilities. This can be particularly seen in Fig. 4, where the evolution of positive axial velocity in the jet entering the chamber and of negative velocity inside the CRZ are plotted simultaneously with the CRZ diameter. Indeed, before t = 15 ms, both velocities and the CRZ diameter are oscillating in phase. Thus, when the axial velocity increases due the passage of the acoustic wave, the CRZ diameter also increases and the jet (and flame), then, opens, corresponding to the classical flame opening oscillations found in swirled flames subject to acoustic oscillations [START_REF] Candel | Dynamics of swirling flames[END_REF]. Nonetheless, after t = 15 ms, when the staging factor reaches zero and the toroidal vortices at the ISL begin to increase in size (as shown in Fig. 2), the velocities and the CRZ diameter begin to evolve in phase opposition. Therefore, as the acoustic pressure wave passes, it increases both positive axial velocity in the fresh gases jet and the negative velocity in the CRZ, but this now leads to the increase in the size of the toroidal vortex at the ISL, that pushes the CRZ inwards, reducing its diameter. The toroidal vortex is then the link between the axial velocity oscillations and the opposite variation of the CRZ diameter.

To understand how the acceleration of the flow by the acoustic pressure wave leads to the creation of the toroidal vortex at the ISL, the ISL dynamics in the initial stable case (α=15%) is analyzed. Studying the stable V flame at α=15% one can better observe how the vortices are formed and how they interact between themselves, with the flame, and with the pressure and flow oscillations due to the acoustics. Figure 5 zooms in the divergent region, focusing on the ISL helical vortices. As these vortices are convected downstream by the flow and an acoustic pressure wave passes (accelerating the flow accordingly), the helical vortices tend to merge until they create a larger vortex, locked on the acoustic frequency, as shown in Figure 5. See [START_REF] Lieuwen | Unsteady Combustor Physics[END_REF] and references therein for details on this mechanism. With the triggering of the combustion instability due to the change on the staging factor and continuous increase in acoustic pressure amplitude, this merging of helical vortices is intensified, thus creating the toroidal vortices shown in Fig. 2 that will lead to the flame shape transition.

Regarding the interaction between the ISL vortices and the flame, in Fig. 5, it is possible to see how the ISL helical vortices wrinkles the flame, especially on the top part of the figure. Therefore, changing the staging factor to α=0% will increase the quantity of fuel over the ISL and, thus, the heat release rate of the premixed flame over the ISL. This will naturally result in the increase of the heat release rate oscillations, due to the natural wrinkling the flame caused by the ISL helical vortices. And, as the combustion instability leads to the merging of the helical ISL vortices in a large toroidal vortex, the later increases even more the wrinkling, re-increasing the acoustic pressure oscillations and closing the coupling mechanism of the instability. The occurrence of the toroidal vortices matches the flame surface variations (as shown in Fig 3), showing their role in altering the flame surface and modifying the heat release rate. This seems to be a case of a classical vortex-driven instability [START_REF] Ducruix | Combustion dynamics and instabilities: Elementary coupling and driving mechanisms[END_REF][START_REF] Poinsot | Vortexdriven acoustically coupled combustion instabilities[END_REF], but in a swirled ISL stabilized flame, where the heat release rate fluctuates due to vortex driven flame wrinkling [START_REF] Caux-Brisebois | Thermoacoustic velocity coupling in a swirl stabilized gas turbine model combustor[END_REF][START_REF] Steinberg | Flow-flame interactions causing acoustically coupled heat release fluctuations in a thermo-acoustically unstable gas turbine model combustor[END_REF][START_REF] Steinberg | Effects of flow structure dynamics on thermoacoustic instabilities in swirlstabilized combustion[END_REF].

Changes in flow topology

As the vortices grow, they push the CRZ limits inwards. Figure 6.a presents the time evolution of the CRZ radius at a distance x = 0.015 m from the chamber entrance. As already discussed, before the transition, the CRZ radius slightly oscillates, because of the ISL vortices and the flow angle oscillations due to the interactions between the acoustic fluctuations and the swirler [START_REF] Candel | Dynamics of swirling flames[END_REF]. By this mechanism, the flame is expected just to oscillate around its stable opening, however, after the pilot injection is stopped, the radius of the CRZ progressively decreases and the transition occurs.

One key question is why do this progressive reduction occurs? The increasingly large toroidal vortices do not explain the CRZ being progressively reduced, they should only act as an oscillating mechanism, pushing the CRZ inwards, that would afterwards return to its original diameter. However, after reaching a certain value, the CRZ radius inevitably decreases and leads to the change in flow topology. The flow is forced to move from a stable topology to another stable one, by gathering sufficient energy from the vortices responsible for the transition. Below is proposed a discussion on a possible mechanism causing the flow topology bifurcation.

The toroidal vortices push the CRZ, consequently reducing its radius and, thus, the incoming fresh gases jet radius, as one can see on the image sequence of Fig. 2. Considering angular momentum conservation, the reduction of the radius of rotating annular mass implies an increase in angular velocity. Figure 6.b confirms this behavior, as the angular velocity increases in magnitude following the radius reduction. To understand how the swirling flow dynamics affect the radial velocity u r , the radial momentum conservation is analyzed. When dealing with vortex breakdown, the centrifugal force and radial depression terms are commonly considered as the most relevant ones for the radial velocity budget ( [START_REF] Lucca-Negro | Vortex breakdown: a review[END_REF][START_REF] Krause | A contribution to the problem of vortex breakdown[END_REF]). This same assumption is used here.

Back to the transition study, CRZ radius reduction and the swirling velocity increase imply an augmentation of the centrifugal force. As a consequence, due to radial momentum conservation, when reaching a critical radius value, this increase in centrifugal force causes an increase in radial pressure gradient. This creates a new equilibrium state for the system, stabilizing the CRZ at the current reduced diameter and not returning the CRZ to its original diameter. The difference between centrifugal force and radial pressure gradient is overall positive until t = 20ms (Fig. 6c), when the toroidal vortices are large enough to start pushing the CRZ inwards. From this point on, approximately 3ms before the transition instant (when the swirling velocity reaches the maximum absolute value before the transition), the difference between centrifugal force and radial pressure gradient becomes negative (Fig. 6c), creating a depression inside the CRZ and progressively reducing the radial velocity component (Fig. 6d). These phenomena still evolve until the CRZ diameter reaches the non-reacting CRZ diameter, then the pressure drop pull combined with the final vortex push are responsible for deviating part of the fresh gas jet inside the CRZ and recreating the fresh gas bubble seen in the non-reacting case. This pushes the flame inwards and changes the flow topology from a conical vortex breakdown to a bubble vortex breakdown and the flame shape from a V-ISL-stabilized-shape to an M-shaped flame stabilized by the OSL and the vortex breakdown bubble.

Conclusion

The LES of a two-stage swirl stabilized liquid-fueled burner have been performed to analyse a singular flame shape transition. This transition is triggered by a change in the injection regime to make the burner operate in LPP conditions. A strong coupling between flame, acoustics and flow dynamics is observed, creating an intricate mechanism. Changing the spatial distribution of fuel triggers a combustion instability, due to the wrinkling of the flame caused by the ISL helical vortices. The increase in acoustic oscillations will favor the natural merging of the ISL helical vortices to create a large toroidal vortex. This vortex then wrinkles more the flame, increasing the amplitude of the instability and closing the loop. This phenomenon selfamplifies at each iteration until a sufficiently large toroidal vortex completely contracts the inner flame, pushes the CRZ inwards and deviates the fresh gases jet inside the CRZ. This causes the flow topology to change from a conical Vortex Breakdown to a Bubble Vortex Breakdown and, thus, the flame to transition from a V-shaped to M-shaped flame. It is, hence, the existence of two possible flow topologies in this burner that will ultimately be the reason for the flame shape transition. The instability triggers the mechanism responsible to transfer enough energy to the flow, enabling it to change its topology. The flame shape changes, then, just to adapt to this new flow topology. The non-existence of multiple possible flow topologies should, therefore, be verified as a primary condition to avoid this type of transition.

Figure 1 :

 1 Figure 1: XZ plane view of the injector at the entrance of the combustion chamber colored by the mesh size. Cyan crosses and arrow indicate, respectively, the pilot swirlers and injector, while the red crosses and red arrows indicate respectively multipoint swirlers and two multipoint injection holes (out of scale). The pressure probe is identified by a red rectangle. Flow from left to right.

Figure 3 :Figure 4 :

 34 Figure 3: (a) time evolution of the integrated heat release rate (red) and the pressure (blue), measured by a probe at the experimental microphone position at (x,y,z)=[0.25, 0, 0.075]m. (b) time evolution of heat release rate iso-surface at a value of 10 7 W level. The letters correspond to the flames shown in Fig. 2. The vertical black line indicates the transition instant. (c) Time evolution of the equivalent radius of the interpolated circle over the inner flame at position X=0.03m (blue). The data stops when the circle diameter reaches zero, moment defined as flame shape transition. The CRZ radius (dotted red line) indicates the correlation between the inner flame and the CRZ. The dot-dashed horizontal line indicates the average CRZ radius for the M flame at this position.

Figure 5 :

 5 Figure 5: (a,b,c) XZ plane instantaneous snapshots of the pseudo-streamlines colored by the pressure showing the vortex-merging in the ISL of the stable V flame. Red lines represent the iso-contour of heat release rate at 10 7 W/m 3 . (d) Normalized PSD of the Q-criterion measured at the position of the vortex 1+2 in (c). Acquisition frequency and time of 10 5 Hz and 40 ms; PSD performed with the Welsh method with 2 equal length Hanning windows.

22

 22

Figure 6 :

 6 Figure 6: Dynamics of a point over the CRZ at x=0.015 m: (a) CRZ radius, (b) velocity swirl component, (c) balance between centrifugal force and pressure radial gradient, (d) velocity radial component. Signals are filtered with a 4th order Butterworth low-pass filter and a frequency cut at 600Hz. The dashed vertical line identifies the moment at which the inner flame has merged in the CRZ.

Acknowledgments

We thank CERFACS for sharing the AVBP solver. This work was granted access to the HPC resources of CINES under the allocation A0012B00164 made available by GENCI, and of the mesocentre computing center of Cen-traleSupélec and Ecole Normale Supérieure Paris-Saclay supported by CNRS and Region Ile-de-France (http://mesocentre.centralesupelec.fr/).


Henry W Orton

Jan Stanek

Tobias Schubeis

Dylan Foucaudeau

Claire Ollier

Adrian W Draney

Tanguy Le Marchand

Diane Cala De Paepe

Isabella C Felli

Diane Cala-De Paepe

Roberta Pierattelli

Sebastian Hiller

Wolfgang Bermel

Guido Pintacuda

Protein NMR resonance assignment without spectral analysis: 5D SOlid-state Automated Projection SpectroscopY (SO-APSY)

Keywords: solid-state NMR, proton detection, NMR resonance assignment, projection spectroscopy, automation

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

One of the major bottlenecks for the widespread use of NMR spectroscopy in structural biology is the difficulty to automate data acquisition and analysis of complex protein spectra. Specific information on individual nuclei is typically achieved by the combination of several 3D spectra which correlate shifts of sequential or proximal 1 H, 15 N and 13 C nuclei. [1] Protein spectra with higher (>3) dimensionality may yield multiple nuclear correlations within fewer experiments, increasing resolution, simplifying resonance assignment and thus pushing forward the size limits of the proteins amenable to site-specific NMR studies. [2] The acquisition of such spectra however requires alternative approaches to sample the indirect time space, since regular sampling leads to impractical experimental durations and data matrices. Randomised non-uniform sampling (NUS) [3] is a sparse sampling technique proposed to solve the first problem. NUS aims to reconstruct a high-dimensional spectrum from data acquired at randomly chosen time points, but it requires particular expertise in processing and analysis. Yet the issue of the size of the data matrix remains. Automated projection spectroscopy (APSY), [4] is an alternative tool which allows direct inference of a high-dimensional peak list from a number of lower order projection spectra (typically 2D). APSY represents an elegant concept which allows to bypass spectral analysis and delivers an output that directly contains the positions of all resonances. As such it is in principle amenable to widespread access even by inexperienced spectroscopists.

APSY in solution NMR spectroscopy has suffered from a major drawback: the sensitivity required to implement APSY can only by provided by globular proteins with molecular sizes smaller than about 20 kDa (fast tumbling), [5] or by intrinsically disordered proteins. [4b, 6] Figure 1 illustrates the rapid loss of sensitivity with increasing molecular size for a 5D experiment in solution of a protonated and highly-deuterated molecule respectively. The deadlock is removed in proton-detected solidstate NMR at fast magic-angle spinning (MAS). Narrow proton signals, high sensitivity, and efficient coherence transfers can be obtained at high magnetic fields, independent of molecular mass, at 60 kHz MAS in deuterated proteins, reprotonated at the exchangeable sites, [7] and at 100 kHz MAS and above in fullyprotonated samples. [8] As illustrated in Figure 1, the efficiency of hyperdimensional experiments in solids outperforms the solution case above 12 and 33 kDa respectively for a protonated and highly-deuterated protein. As demonstrated below, this makes five dimensional correlations by SOlid-state APSY (SO-APSY) feasible on microcrystalline and fibrillar proteins, which in turn allows their expeditious automated backbone resonance assignments.

We have developed a 5D experiment (H)NCOCANH (Figure 2a) which transfers coherences from each amide proton to the amide proton of the previous residue. This unidirectional coherence pathway is ensured by two cross-polarization (CP) steps between 1 H and 15 N nuclei, two selective CP steps between bonded 15 N and 13 C / 13 C' nuclei, and one 13 C -13 C' scalar transfer. In contrast to previous implementations, [9] this experiment evolves the chemical shift of each of the nuclei ( 13 C , 1 H spectra. The curves show the calculated efficiency of a 5D amide proton linking experiment on a protonated molecule in solution (pink dotted curve), of the TROSY based variant on a highly deuterated molecule in solution (brown dashed curve) and of the analogous experiment in the solid state at fast MAS (cyan solid curve). Calculations assumed scalar transfers in solution, with transverse relaxation modelled using Redfield theory, and a mix of scalar and dipolar transfers in solids (see SI for details).

[a] Dr. J. [4a]

13 C', 15 N) traversed along the protein backbone. The experiment generates a 5D dataset where each peak correlates the chemical shift of five consecutive backbone nuclei (Figure 2b). It can be acquired as a SO-APSY experiment by collecting a number of 2D projections, each characterized by projection angles , , and within the 4D indirect time space (Figure 2c-e).

The method was first established on a sample of uniformly 15 N-and 13 C-labelled and fully protonated microcrystalline GB1 protein (56 aa), packed in a 0.7 mm rotor and spun at 111 kHz on a 23.5 T NMR spectrometer ( 1 H frequency of 1 GHz). Under these conditions, the efficiency the 5-step transfer is ensured by long R1 and coherence lifetimes, and the sensitivity amounts to ~20% with respect to a basic experiment correlating 1 H- 15 N nuclei. The SO-APSY series was acquired using a newly developed TopSpin module for about 1.5 days at which point convergence was achieved with 24 projections (no additional 5D peaks were found by including the last projections). The first four projections, where the chemical shift of a single heteronucleus is evolved in the indirect dimension, are shown in Figure 3. Contrary to fully dimensional spectral data obtained with NUS, identifying resonances in APSY is straightforward and does not require prior knowledge. [10] The Geometrical Analysis of PROjection (GAPRO) algorithm [4a] identifies robust peaks, as long as they are supported (resolved) in at least Smin projections (here 7), and the on-the-fly analysis yielded a list with all 54 expected 5D peaks (Table S8).

The experiment was then acquired on two, more challenging targets, a much larger microcrystalline protein (the 32 kDa Cu(I)-Zn loaded dimeric superoxide dismutase, SOD) [START_REF] Knight | Proc. Natl. Acad. Sci[END_REF] and a fibrillar aggregate (D76N 2 microglobulin, 2mD76N). [START_REF] Marchand | [END_REF] Both samples were uniformly 15 N, 13 C and 2 H labelled with 100% backexchange of labile amide protons, and packed in a 1.3 mm rotor spun at 60 kHz MAS on a 1 GHz spectrometer.

For SOD, which contains 153 residues, but several dynamic regions and proline residues, a maximum of 115 peaks can be expected in the 5D experiment. The SO-APSY series was accomplished in 32 projections in 3 days. Despite the significant degree of peak overlap in the first (orthogonal) projections (Figure 3, middle row), the tilted projections contributed a total of 95 peaks. For 2mD76N, spectral overlap does not come from protein size (only 62 residues S20-R81 are observed in solids), but from limited dispersion of chemical shifts and increased linewidths, characteristic for fibrillar samples. The SO-APSY series was acquired with 32 projections in about 2.5 days. GAPRO identified 46 correlations, proving that consistent 5D peak coordinates can be determined even within regions of severe overlap in orthogonal 2D projections (Figure 3, bottom).

The 5D experiment has the particular advantage over conventional 3D triple-resonance spectra that a single peak encodes amide shifts from two neighboring residues. This alone allows sequential linking of amide resonances along the protein backbone by matching shared 15 N frequencies. In cases where 15 N resolution is limited or an amide proton is unobservable, assignment of sequential residues can be accomplished with the acquisition of a second 5D experiment, (H)NCACONH. It performs a similar coherence transfer along the protein backbone, but in a reverse (i.e. forward) direction (Figure 4a), with comparable efficiency (Table S1). Such experiment was indeed acquired as a SO-APSY series for the three samples and yielded an equivalent number of signals (Table S9). Peaks from both experiments overlay in the 15 N(i)-15 N(i+1) plane, and effectively connect consecutive 1 H N (i) and 1 H N (i+1) amide protons (two-way 15 N linking). [9] A pair of forward and backward 5D data sets also correlate consecutive amide (i and i+1) protons to pairs of 13 C (i)- 13 C'(i) frequencies, effectively building quartets of shifts ( 15 N(i), 13 C (i), 13 C'(i), 15 N(i+1)) that unambiguously link amide protons at every second step of the assignment protocol, as shown in Figure 4c.

Five-dimensional peak lists are inconvenient for a conventional manual analysis, and their potential for sequential assignment is optimally exploited with a powerful automated assignment algorithm. We employed the FLYA software [13] which is particularly flexible with respect to customized correlation schemes (see the SI for FLYA definitions, parameters and input data). In the case of GB1, a complete assignment of all backbone nuclei was obtained when providing two 5D peak lists from a pair of forward and backward SO-APSY experiments. The FLYA output schematic which compares the assigned peaks to a reference assignment [8a] is shown in Figure 5a.

In the case of SOD, the automated analysis by FLYA reveals 8 segments of robust ("strong") assignments (Figure S3a), with 280 chemical shifts consistent with published data. [START_REF] Knight | Proc. Natl. Acad. Sci[END_REF] Assignments were not found in regions 23-27 and 129-140, N pairs (odd steps) or quartets of 15 N(i), 13 C (i), 13 C' and 15 N(i+1) chemical shifts (even steps).

known to be dynamic, [START_REF] Knight | Proc. Natl. Acad. Sci[END_REF] and signals were sparse in the regions 50-70 and 100-110 due to the presence of five proline residues. An almost complete assignment can be recovered for these regions (Figure 5b), if the 5D peak lists are supplemented with information on residue type, as for example that contained in 3D C -N-H correlations (in total 346 correct assignments). About 4% of chemical shifts were incorrectly assigned when compared to a reference assignment obtained by conventional proton-detected 3D spectra. This reflects an inherent uncertainty that arises during automated assignment, which must be considered when inferring firm conclusions on protein site-specific dynamics or interactions. FLYA analysis of 5D peak lists obtained for 2mD76N leads to three regions of reliable assignments at the fibril core (27-31, 34-45, 59-65), while many regions with reasonably complete observed correlations are assigned only tentatively (Figure S3b). We attribute this result to the significantly smaller distribution of chemical shifts in this sample and their decreased precision (increased line widths). Similar to SOD, this issue can be addressed by providing supplementary input for the automatic analysis. Notably, an additional unassigned 3D C -N-H peak list allows to extend the assignment to most of the fibril core, including regions 46-49, 56-58, 67-71 and 73-75 (Figure 5c) giving rise to 166 correctly assigned chemical shifts and 5 misassignments (3%).

Although a small number of incorrect assignments were identified within this study when compared to manually assigned 3D spectra, they do not form clusters spanning more than 1 or 2 residues, and are frequently at the extremities of long tracts of assigned residues. In effect, these misassigned shifts do not affect the prediction of dihedral angles , , which is the primary structural information encoded in 1 H N , 15 N, 13 C , 13 C and C' shifts. Indeed, when using the widely applied TALOS+ algorithm, [14] the dihedral angles obtained from chemical shifts are virtually identical, with only 1 outlier ( for residue T137 of SOD) among 214 angles reliably predicted in both reference and FLYA-derived data sets for SOD and 2mD76N (Figure 5d,e). Additionally, it has been shown that reliable structural modelling can be obtained with up to 10% missing or incorrect proton chemical shifts provided that robust algorithms (such as CANDID [15] ) are employed. [16] Therefore we would expect 5D APSY of yet unassigned proteins in the solid state to offer appropriate data for subsequent modelling of structure and dynamics.

15 N i+1 ppm 13 CO i ppm 13 CA i ppm 15 N i ppm In summary, we have presented the first example of ultrahigh dimensional (>3) automated projection spectroscopy in solid-state NMR (SO-APSY) and demonstrated its power with respect to automatic resonance assignment of samples of different molecular size and aggregation state. The approach was illustrated with the design of two new 5D pulse sequences that correlate backbone nuclei with amide protons. A module to acquire 2D projections and to reconstruct 5D peak lists was designed within the widespread software TopSpin, and is now readily accessible to the NMR community. SO-APSY tolerates higher peak overlap and yields data prone to automatic analysis, thus accelerating the identification of sequential correlations. We anticipate that in a context where even faster MAS rates and longer coherence lifetimes are available, SO-APSY techniques will become a key for rapid and simple spectroscopic analysis, broadening the range of targets amenable to NMR.

1 H i / ppm 1 H i / ppm 1 H i / ppm 1 H i / ppm D O S 1 B G β2m D76N 50 nm 100 µm Cα C N O H α H C β N (a) (b)

Figure 1 .

 1 Figure 1. Sensitivity of hyperdimensional NMR plotted over molecular weight in solution and in the solid state with respect to 1D1 H spectra. The curves show the calculated efficiency of a 5D amide proton linking experiment on a protonated molecule in solution (pink dotted curve), of the TROSY based variant on a highly deuterated molecule in solution (brown dashed curve) and of the analogous experiment in the solid state at fast MAS (cyan solid curve). Calculations assumed scalar transfers in solution, with transverse relaxation modelled using Redfield theory, and a mix of scalar and dipolar transfers in solids (see SI for details).

Figure 2 .

 2 Figure 2. 5D (H)NCOCANH experiment and schematic for SO-APSY acquisition and analysis. a) Time diagram of RF irradiation. High-power 90° and 180° are represented by black and white rectangles, respectively, and selective pulses by bell shapes. Cross-polarization transfers are represented as open shapes on two parallel channels, and grey rectangles indicate heteronuclear decoupling. The five chemical shift evolution periods are highlighted as colored fields. Additional details can be found in the SI. b) Coherence transfer pathway. c) Time-domain sampling scheme for the SO-APSY experiment in the indirect 4D space, where each projection is characterized by three polar angles , , and . d) Frequency-domain orthogonal projections. e) Reconstruction of the peak positions in the 5D space by Geometrical Analysis of PROjection (GAPRO algorithm).

Figure 3 .

 3 Figure 3. Orthogonal SO-APSY projections from 5D (H)NCOCANH experiment of microcrystalline GB1 and SOD, and fibrils of 2mD76N. These projections exclusively evolve a single nucleus in the indirect dimension: from left to right, intraresidue (i)15 N,13 C ,13 C', and sequential (i+1)15 N with respect to the directly detected amide proton (i) shifts. Crosses indicate experimental peaks identified by GAPRO algorithm, and peak labels indicate residue assignments obtained from automated analysis (see below).

Figure 4 .

 4 Figure 4. Sequential resonance assignment from a pair of 5D backward and forward experiments. (a) Coherence transfer pathway of a forward (H)NCACONH experiment. (b) Time diagram of RF irradiation, with the same conventions of Figure 2a. (c) Schematic of resonance assignment, obtained by matching either 15 N, 1 HN pairs (odd steps) or quartets of 15 N(i), 13 C (i),13 C' and 15 N(i+1) chemical shifts (even steps).

Figure 5 .

 5 Figure 5. Graphical summary of automated assignments of microcrystalline GB1 (a), SOD (b) and fibrils of 2mD76N (c) by FLYA, superimposed to the amino acid sequence and to a sketch of secondary structure elements. Assignments in agreement or incompatible with the reference are colored in green and red, respectively. Assignments with low confidence levels in FLYA are colored in light green or in grey, depending on whether they can or cannot be confirmed by manual inspection. Additional assignments found by FLYA are shown in dark blue. Black arrows represent all correlations observed in the two 5D experiments. Dotted lines connecting grey arrows ( -strands) and cylinders (helices) denote flexible regions where resonances are missing or sparse. (d,e) Correlation of backbone dihedral angles , predicted by TALOS+ based on reference and automatically obtained 1 H N , 15 N, 13 C , 13 C and 13 C' chemical shifts for SOD (d) and 2mD76N (e).

Acknowledgements

The work was funded by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (ERC-2015-CoG GA 648974), by the CNRS (IR-RMN FR3050), and by the EU-project iNext (GA 653706). HO was supported by the Westpac Foundation with a Future Leaders Scholarship, and JS by the EC's REA with a MSCA fellowship (GA 661799). The project was co-financed by the Polish National Agency for Academic Exchange (contract No PPN/PPO/2018/1/00098).