Ahmed Khalifa

Olivier Kermorgant

Salvador Dominguez

Philippe Martinet

Vehicles Platooning in Urban Environments: Integrated Consensus-based Longitudinal Control with Gap Closure Maneuvering and Collision Avoidance Capabilities

This paper proposes a distributed longitudinal controller for car-like vehicles platooning that travel in an urban environment. The presented control strategy combines the platoon maintaining, gap closure, and collision avoidance functionality into a unified control law. A consensus-based controller designed in the path coordinates is the basis of the proposed control strategy and its role is to achieve position and velocity consensus among the platoon members taking into consideration the nature of the motion in an urban environment. For platoon creation, gap closure scenario is highly recommended for achieving a fast convergence of the platoon. For that, an algorithm is proposed to adjust the controller parameters online. A longitudinal collision between followers can occur due to several circumstances. Therefore, the proposed control strategy considers the assurance of collision avoidance by the guarantee of a minimum safe inter-vehicle distance. Convergence of the proposed algorithm is proved in the different modes of operations. Finally, studies are conducted to demonstrate and validate the efficiency of the proposed control strategy under different driving conditions. To better emulate a realistic setup, the controller is tested by an implementation of the car-like vehicles platoon in a vehicular mobility simulator called ICARS, which considers the real vehicle dynamics and other platooning staff in urban environments.

I. INTRODUCTION

The car-like vehicles platooning in urban environments is one of the current trends to solve the problem of the cars redistribution in the car sharing system. This is achieved by utilizing vehicle platoon guided by professional drivers who come to pick up and drop off vehicles over the stations [START_REF] Khalifa | Vehicles platooning in urban environment: Consensus-based longitudinal control with limited communications capabilities[END_REF].

In the literature, variety of platoon models are presented including linear and nonlinear models. Based on different assumptions, three linear models are used including single integrator model [START_REF] Lin | Optimal control of vehicular formations with nearest neighbor interactions[END_REF], double integrator model [START_REF] Santini | A consensus-based approach for platooning with intervehicular communications and its validation in realistic scenarios[END_REF], and third order model [START_REF] Chehardoli | Third-order leaderfollowing consensus protocol of traffic flow formed by cooperative vehicular platoons by considering time delay: constant spacing strategy[END_REF]. Few studies have addressed the problem of vehicles platoon in an urban environment, e.g., [START_REF] Li | Minimum fuel control strategy in automated car-following scenarios[END_REF]. In such studies, the authors uses a first order model of the platoon which is widely different from the real vehicle dynamics. In our previous research [START_REF] Khalifa | Vehicles platooning in urban environment: Consensus-based longitudinal control with limited communications capabilities[END_REF], a second order dynamic model is utilized and transformation to its corresponding one in the path coordinates is given so that we can design a controller to satisfy control objectives defined in the operational coordinates.

Adaptive Cruise Control (ACC) is presented as the simplest platoon controller, and it is role is to maintain a safe distance with the preceding vehicle. Recently, platoon control is designed based on a distributed controller known as Cooperative Adaptive Cruise Control (CACC) [START_REF] Vanderwerf | Modeling effects of driver control assistance systems on traffic[END_REF] which needs information from the neighbor vehicles to achieve a global consensus. Six main approaches are presented to design a platoon controller including Linear [START_REF] Chehardoli | Third-order leaderfollowing consensus protocol of traffic flow formed by cooperative vehicular platoons by considering time delay: constant spacing strategy[END_REF], Optimal [START_REF] Fardad | Sparsity-promoting optimal control for a class of distributed systems[END_REF], H ∞ -based [START_REF] Gao | Robust control of heterogeneous vehicular platoon with uncertain dynamics and communication delay[END_REF], Sliding Mode [START_REF] Wu | Distributed sliding mode control for multi-vehicle systems with positive definite topologies[END_REF], Model predictive [START_REF] Dunbar | Distributed receding horizon control of vehicle platoons: Stability and string stability[END_REF], and Consensus-based controller [START_REF] Santini | A consensus-based approach for platooning with intervehicular communications and its validation in realistic scenarios[END_REF]. However, none of these works tackles the control in the Curvilinear coordinates. Furthermore, they assume that the leader navigates with constant velocity which is not feasible to the motion in urban environments in which it travels with variable velocities. Therefore, in our previous research [START_REF] Khalifa | Vehicles platooning in urban environment: Consensus-based longitudinal control with limited communications capabilities[END_REF], the authors propose a controller that can achieve asymptotic stability of the tracking errors and the platoon string stability considering variable velocity of leader.

For platoon creation, gap closure scenario, between a merging vehicle and an in-platoon vehicle or at the end of a splitting scenario, is required. Hence, a vehicle should be able to cooperatively close or open a gap. The gap closure requirement is to increase the response speed in order to decrease the gap closing time. In other words, at large inter-vehicular distance, a reduced damping can be desirable to allow significant velocity differences for gap closing performance. Therefore, in this research, platooning and gap closing functionality is combined into a single control design.

Recently, control problems with collision avoidance have got large progress. Many approaches have been developed in this field. Firstly, Rule-based approach (e.g., Artificial Potential Field [START_REF] Zavlanos | Potential fields for maintaining connectivity of mobile networks[END_REF]). Secondly, optimization-based approach (e.g., Model Predictive Control [START_REF] Dunbar | Distributed receding horizon control of vehicle platoons: Stability and string stability[END_REF]). Thirdly, with the aid of time-optimal path-constrained trajectory planning techniques, for instance, the authors in [START_REF] Shen | Platoon of autonomous vehicles with rear-end collision avoidance through time-optimal pathconstrained trajectory planning[END_REF] present a trajectory planning algorithm to generate time optimal velocity profiles, which assures the rear-end collision avoidance and considers both velocity and acceleration constraints for the leader and follower vehicles. In this approach, the problem is handled at the path generation level, i.e., the system has a control on the path of the leader and followers. However, this is not the case for our application in which the followers follow the path generated from the leader that has a free motion. Fourthly, in [START_REF] Ferrara | Collision avoidance strategies and coordinated control of passenger vehicles[END_REF], the vehicles are equipped with a collision detection and avoidance module, like pedestrians crossing the road or other vehicles, at the low level controller. Consensus-based controller with collision avoidance capability applied to quad-rotors formation control is introduced in [START_REF] Kuriki | Consensus-based cooperative formation control with collision avoidance for a multi-uav system[END_REF] in which the authors propose collision avoidance strategy based on the potential field method. However, they use a potential function that has a switching effect (i.e., the transition from the non-collision case to collision case and vice versa is discontinuous). In our study, the platoon longitudinal controller is supported by a collision avoidance capability so that a safe inter-vehicle distance, in case of abnormal situation occurs such as sudden braking of a vehicle inside the platoon, can be guaranteed. Therefore, the main high-level longitudinal control law is modified in order to make the platoon able to guarantee a safety region to keep a minimum inter-vehicle safe distance and to avoid the collision. Moreover, we use a potential function that has smooth and local effect during the transitions.

The proposed controller is implemented in a home made fast prototyping environment called ICARS that is used for simulation. It is developed under Robot Operating System (ROS) and it has features such as real vehicle dynamics and models (engine, transmission, braking systems, etc.), both low and high level measurements units with noises, low level communication module between different parts of the vehicle and the on-board computer, real maps and urban environments with different road conditions (curvature, bumps, etc.), vehicle communication modules, and fully compatible with embedded computer architecture.

Summarizing, the main contributions of this research are,

• design of a gap closure algorithm;

• design of a longitudinal anti-collision strategy which guarantees a safety zone between the vehicles; • combining platooning, gap closure, collision avoidance functionalities into a single control law; • proof of the system convergence.

II. PROBLEM STATEMENT

A. Kinematic Model of Car-like Vehicle

The kinematic model of a car-like vehicle can be simplified to that of a bicycle model [START_REF] Khalifa | Vehicles platooning in urban environment: Consensus-based longitudinal control with limited communications capabilities[END_REF], see Fig. 1, in which the left and right wheels are considered as a single wheel at the center of front and rear axles. 1) Cartesian Coordinates: Let {V }, O v -x v y v , represents the vehicle body-fixed reference frame with its x-axis is in the vehicle forward direction and its origin is at the center of its rear axle, O v , see Fig. 1. The vehicle configuration, with respect to the world-fixed inertial reference frame, {G}, O -xy, is represented by the generalized coordinates p = [x y θ] ∈ R 3 , where [x y] T is the position and θ represents its orientation. The vehicle steering angle is given by δ. Assuming that the wheels cannot slip sideways, the vehicle velocity in the body x-direction is v, and it is zero in the body y-direction. The front wheel is mounted at a distance L from the rear wheel. The kinematic model is given by

ẋ = v cos(θ), (1a) ẏ = v sin(θ), (1b) θ = v L tan(δ). (1c)
2) Curvilinear Coordinates: The vehicle kinematics can be represented in the Curvilinear coordinates (s r ψ) T , as illustrated in Fig. 1. The path to be tracked is C. The path length s (Curvilinear abscissa) is defined between the closest point M to O v . The angular deviation of the vehicle with respect to the path is resented by ψ = θ -θ c , where θ c is the angle between path tangent at M and the x-axis. The lateral deviation, r, is defined as the signed orthogonal distance from the rear axle center, O v , to point M .

Given the geometry in Fig. 1, the kinematic model in the path coordinates is given as

ṡ = v cos(ψ) 1 -rκ(s) , (2a)
ṙ = v sin(ψ), (2b)
ψ = v tan(δ) L - κ(s) cos(ψ) 1 -rκ(s) , (2c)
where κ(s) is the path curvature at point M .

B. Longitudinal Dynamics

The vehicle longitudinal dynamics are affected by many factors such as the engine, drive line, brake system, aerodynamics drag, tire friction, rolling distance, and gravitational force, etc [START_REF] Wu | Distributed sliding mode control for multi-vehicle systems with positive definite topologies[END_REF]. For the vehicle dynamics, let us assume: Assumption 1. The body is rigid and symmetric.

Assumption 2. The throttle and brake inputs are combined as one control input.

Assumption 3. The effect of pitch and yaw motions is neglected.

Therefore, the longitudinal dynamic equation can be represented by

M v a + C v v + G v = F v , (3)
where a = v denotes the vehicle acceleration in the body frame, and M v , C v , G v , and F v ∈ R are terms representing the vehicle Inertia, Coriolis and Gravity effects, and input force, respectively. For the purpose of longitudinal HLC, The longitudinal LLC is considered as an inverse model-based control law which is given by

F v = M v µ + C v v + G v (4)
where µ is the new input signal to be designed. Therefore, the vehicle linear longitudinal dynamic model is given as

η = µ, (5)
Since we aim at platooning in an urban environment, relation between the vehicle dynamics in the Curvilinear and Cartesian coordinates will be given as follows.

Let us rewrite (2a) as

ṡ = Jv, (6)
Thus, the vehicle acceleration in the Curvilinear coordinates, η, is given as

η = Ja + Jv, (7)
where η = q = s and q = ṡ are the vehicle acceleration and velocity in the path coordinates respectively. Let us denote the control signal in the path coordinates as u. If the relation between u and µ is given by

µ = 1 J (u -Jv), (8)
then, the path coordinates longitudinal dynamics are given by ṡ = q, (9a) q = η = u.

(9b)

C. Platoon Longitudinal Model

Consider N + 1 platoon vehicles as illustrated in Fig. 2, which consists of a leading vehicle (called as the leader with index 0) and N following vehicles (called followers with index i; i = 1, 2, . . . , N ). The platoon runs on urban and flat road, and the vehicles can share information by either sensor-based or communication-based way.

2: Car-like vehicle platoon representation.

Since the vehicles considered in the car sharing system are small passenger cars, one can assume that they have equal dynamics. The in-platoon vehicle dynamic model in the Curvilinear coordinates is given by ṡi = q i , (10a) qi = u i .

(10b)

For vehicles string, let us define the inter-vehicle distance as d i = s j -s i that is the actual relative Curvilinear distance between vehicle i and its predecessor j = i -1. In addition, the desired inter-vehicle distance between vehicles i and j is represented as d r,i . Many spacing policies are proposed in the literature [START_REF] Li | An overview of vehicular platoon control under the four-component framework[END_REF]. For higher traffic flow, the constant spacing policy method is utilized in our research.

D. Control Objectives

The platoon can navigate in an urban environment/Curvilinear path and track a leading vehicle that is either autonomous/human driven, and each follower vehicle, i, can receive leader position, velocity, and acceleration represented by s 0 , q 0 , η 0 respectively and predecessor position, s i-1 . Under these circumstances, the objective is to combine platooning, gap closing, and collision avoidance functionality into a single control design. Therefore, we aim at achieving the following goals:

• convergence of the tracking errors as well as string stability; • gap closure requirement: to increase the response speed in order to decrease the gap closing time for large intervehicle distances; • collision avoidance: to guarantee a safety region between the vehicles (keep a minimum safe intervehicle distance in case of an abnormal situation).

III. PLATOON CONTROL DESIGN

Let G = (V, E, M) be a digraph describing the information flow topology among the in-platoon vehicles. The graph has set of nodes V = 1, . . . , N and set of edge E ⊆ V × V. Each graph node represents a vehicle whose dynamics are defined by [START_REF] Dunbar | Distributed receding horizon control of vehicle platoons: Stability and string stability[END_REF]. Let M = [m ij ] N ×N be adjacency matrix associated with G. Its elements are defined as follows: if vehicle i can obtain information from vehicle j, then m ij = 1; otherwise m ij = 0. Therefore, it represents graph topology. An augmented digraph Ḡ is also defined in order to involve the leader in the network.

We propose a control law given by

u i = u p i + u c i , (11)
where u p i is the part of the control signal in charge of the platoon creation and maintaining, while u c i is responsible for the collision avoidance.

A. Design of u

p i u p i = η 0 + b(q 0 -q i ) + N j=0 k ij m ij (s j -s i -d r,ij ), (12)
where, b and k ij are tuning parameters, and d r,ij is the desired inter-vehicle distance between vehicles i and j.

Since we assume a PLF topology, ( 12) can be reformulated as u p i = η 0 + be q,i0 + k 0 e s,i0 + k 1 e s,i ,

where e s,i0 = s 0 -s i -d r,i0 , e q,i0 = q 0 -q i , e s,i = s js i -d r,i , and k 0 and k 1 are controller parameters. 1) Internal Stability Analysis: If the control law ( 13) is applied to the platoon system [START_REF] Dunbar | Distributed receding horizon control of vehicle platoons: Stability and string stability[END_REF], and defining the intervehicle distance error, e s,i , in terms of the errors with respect to leader, e s,i0 , as e s,i = s j -s i -d r,i = e s,i0 -e s,j0 , d r,i = d r,i0 -d r,j0 , j = i -1, then the error dynamics are given as ės,i0 = e q,i0 , ėq,i0 = -be q,i0 -(k 0 + k 1 ) e s,i0 + k 1 e s,j0 .

Let us define e s,0 = [e s,10 , . . . , e s,i0 , . . . , e s,N 0 ] T and e q,0 = [e q,10 , . . . , e q,i0 , . . . , e q,N 0 ] T as the vectors of position and velocity errors, respectively. Thus, the error state vector, X (t) = [e s,0 e q,0 ] T , is given as

Ẋ (t) = AX (t), ( 15
)
where A is given by A = O N I N -P -D , I n and O n denote

(n × n) identity and (n × n) null matrices, respectively, D = diag{b, . . . , b} ∈ R N ×N , and P = [p ii ] ∈ R N ×N is represented by p ii =            k 0 , j = i, i = 1, k 0 + k 1 , j = i, i > 1, -k 1 , j = i -1, i > 1, 0, otherwise. (16)
The eigenvalues of P , λ P,i , are λ P,1 = k 0 and λ P,i = k 0 + k 1 , ∀i > 1. Thus, ∀ k 0 > 0 and k 1 > 0, the matrix P is positive stable. In addition, ∀ b > 0, the matrix D is positive stable.

Exploiting the Schur's formula [START_REF] Cao | Distributed discrete-time coordinated tracking with a time-varying reference state and limited communication[END_REF], the characteristics equation of A is calculated as

det (SI 2N -A) = N i=1 S 2 + bS + λ P,i . (17)
The polynomial S 2 + bS + λ P,i is Hurwitz stable under the conditions b > 0, k 0 > 0, k 1 > 0, i.e., the matrix A is stable. Consequently, the dynamics are asymptotically stable i.e., lim t→∞ X (t) = 0.

2) String Stability Analysis: To simplify our analysis, let us define e i = s j -s i -d r,i , ėi = q j -q i , ëi = u j -u i , e i0 = s 0 -s i -d r,i0 , and ėi0 = q 0 -q i .

Form [START_REF] Ferrara | Collision avoidance strategies and coordinated control of passenger vehicles[END_REF] one can write ëi as

ëi = u j -u i = b ėj0 + k 0 e j0 + k 1 e j -b ėi0 -k 0 e i0 -k 1 e i , (18)
Writing [START_REF] Thuilot | Accurate automatic guidance of an urban electric vehicle relying on a kinematic gps sensor[END_REF] in the S-domain results in

H(S) = E i (S) E j (S) = k 1 S 2 + bS + c , (19)
where

c = k 0 + k 1 .
The L ∞ string stability can be guaranteed by satisfying h(t) 1 < 1 and h(t) > 0, where h(t) = L -1 {H(S)} is the impulse response of H(S).

The poles of H(S) are p 1,2 = -0.5b ± (0.5b) 2 -c. Therefore, to achieve h(t) > 0, the response of H(S) have to be over-damped, i.e., (0.5b) 2 > c, or to be critical-damped, i.e., (0.5b) 2 = c. However, for a higher speed response, the parameters of H(S) are designed to achieve a criticaldamped response, i.e., (0.5b) 2 = c. In such as case, we have h(t) = k 1 te -0.5bt , so ∀ k 1 > 0, we satisfy h(t) > 0.

Lemma 1. [START_REF] Darbha | String stability of interconnected systems: An application to platooning in automated highway systems[END_REF]. If G(s) is a stable proper transfer function with g(0) = δ, and if g(t) is its impulse response, then g(t) 1 = δ if and only if g(t) always has the same sign.

Exploiting Lemma 1, we have

h(t) 1 = H(0) = k 1 c = γ. (20)
Consequently, by setting k 1 = γc with γ < 1, then h(t) 1 < 1 and as a result the platoon is strong string stable.

3) Transient Response: Let us define the eigenvalues of the system matrix A as λ A,i , so from (17), λ A,i is given by

λ A,i1,2 = -0.5b ± (0.5b) 2 -λ P,i (21)
Comparing ( 21) with the standard form of the characteristic equation of second order system, then the settling time is with c = k 0 + k 1 and k 1 = γc. Therefore, during the gap closing, if the error e s,i is larger than a certain value, e u , smaller value for ζ, ζ l , is allowed to increase the rising time and a higher weight is given for e s,i than that for e s,i0 via changing the value of γ to γ u = 1. As soon as, the error approaches e l this will be an indication for the end of gap closure scenario and the value of ζ and γ converge to the values required for string stability, ζ u = 1 and γ l < 1, respectively. Consequently, a smooth commutation between the normal platooning mode and the gap closing mode can be done by the calculation of ζ and γ as follows

T s = 4τ = 8 b , damping factor ζ = b 2 √ λ . 4
ζ =          ζ u , e s,i ≤ e l , ζ l , e s,i ≥ e u , ζu-ζ l 2 1 + cos π es,i-e l eu-e l + ζ l , otherwise. (22) γ =          γ l , e s,i ≤ e l , γ u , e s,i ≥ e u ,
γu-γ l 2

1 + cos π es,i-eu eu-e l + γ l , otherwise.

(23) To evaluate the gap closure scenario, it is required to quantify the corresponding time required to reduce the gap. Thus, the absolute error is integrated over time. For such performance evaluation, let us consider a discrete time span defined by t ∈ [t 0 , t K ] with equal segments, K, of length ∆t k .

Q gc,i = K k=1 e s,i (t) ∆t k .
(24)

B. Design of u c i

In this part, the controller part u c i is designed based on the artificial potential field that guarantees the collision avoidance of vehicle i with its predecessor.

The u c i can be designed for longitudinal collision avoidance as follows.

u c i = -∇ si ϕ i = - ∂ϕ i ∂s i , ( 25
)
where ϕ i is the artificial potential field produced by vehicle i.

In the literature several potential field approaches have been studied to guarantee the collision avoidance of multiagent systems [START_REF] Zavlanos | Potential fields for maintaining connectivity of mobile networks[END_REF]. Let us define the potential field as

ϕ i = 1 Γ kc i , (26)
where

d i = s i-1 -s i , k c > 0 is a tuning parameter,
and Γ i is the potential function that is adapted as

Γ i =     1 -α |d i | 2 -d 2 s 2 1 + |d i | 2 -d 2 s 2     ρ , (27)
where, d s is the desired safe inter-vehicle distance and α = . α is chosen such that the function Γ i is differentiable at d s , so the gradient of ϕ i remains continuous everywhere. The value of Γ i lays in the interval [0, 1], becoming 0 when d i = 0, and 1 when d i > d s . Clearly, ϕ i → ∞ whenever Γ i → 0. Thus, moving in the direction of the negative gradient of ϕ i guarantees that the collisions between the cars will be avoided. The reason for choosing Γ i as in ( 27) is two fold. First, it is a smooth function that is the transition between none-collision and collision states is done smoothly. Second, it has local effect, i.e., the collision avoidance part is ineffective once the car is sufficiently far from the predecessor. The exponent ρ is chosen as ρ = 1-sign(|di|-ds) 2 to force Γ i to be constant beyond a distance d s , so the collision avoidance part will only affect the value of control signal if d i approaches d s .

The main result of this paper can be stated as follows.

Theorem 1. Consider the platoon system whose longitudinal model defined in (1), ( 2), ( 3) -(8), and (10), under assumptions 1 -3, and the proposed control algorithm described in (11), ( 13), ( 22), (23), and (25) is applied. Then, with both gap closure and collision avoidance capabilities, the platoon tracking error convergence and string stability are guaranteed as long as, ∀ γ < 1, the following conditions are satisfied: b > 0, c = 0.25b 2 , k 1 = γc, and k 0 = (1-γ)c.

Proof. Applying the control law [START_REF] Zavlanos | Potential fields for maintaining connectivity of mobile networks[END_REF] to the system [START_REF] Dunbar | Distributed receding horizon control of vehicle platoons: Stability and string stability[END_REF], then the error dynamics [START_REF] Li | An overview of vehicular platoon control under the four-component framework[END_REF] become

Ẋ (t) = AX (t) + BU c , (28)
where

B = O N I N and U c = [-u c 1 , -u c 2 , . . . , -u c N ]
T . There are two cases of operation, including normal case (i.e., no collision is detected d i > d s ) in which U c = 0, and collision case in which U c = 0. For the normal case and from the stability analysis presented in the previous sections, one can set the following conditions: ∀ γ < 1, b > 0, c = 0.25b 2 , k 1 = γc, and k 0 = (1 -γ)c, for both platoon asymptotic stability and string stability.

For the collision case (d i < d s ), U c = 0, so if one proves the boundedness of U c , one can deduce that the platoon will be stable with steady state error, i.e., lim t→∞ Ẋ (t) = 0 and lim t→∞ X (t) = 0. Thus, next, we will prove that U c is bounded.

Let us recast (26) as

ϕ i = (β i ) -ρkc , (29)
where β i is given by

β i = 1 -α |d i | 2 -d 2 s 2 1 + |d i | 2 -d 2 s 2 . (30)
Thus,

∇ si ϕ i = -ρk c (β i ) -ρkc-1 ∇ si β i (31)
From (30), (β i ) -ρkc-1 is bounded provided that |d i | = 0. in the rest of the proof will find the bound of ∇ si β i . Let us define

w i = |d i | 2 -d 2 s
, so (30) can be recast as

β i = 1 -α w 2 i 1 + w 2 i . (32)
Consequently, ∇ si β i is given by

∇ si β i = ∇ wi β i ∇ si w i = -2αw i (1 + w 2 i ) 2 2 w i + d 2 s (33)
From (33), one can deduce

|∇ si β i | < (1 + w 2 i ) 2 |∇ si β i |. (34)
The maximum value of the term (1 + w 2 i ) 2 |∇ si β i | can be easily found using the conventional method which equals to at w i = -2 3 d 2 s . Thus, the value of |∇ si β i | is bounded. As a result, the control signal, u c i is bounded and hence U c is bounded.

IV. RESULTS

A. Environment

The considered urban cars platoon with the proposed control strategy is simulated in ICARS framework. The platoon under study includes 3 follower vehicles (Model: Robotized Renault ZOE) plus a manually driven leader (Model: Renault FLUENCE electric equipped with sensors and communication module). The parameters of the vehicle is given in Table I. The vehicles are initially parked and we use the maps of Ecole Centrale de Nantes campus. The leader moves from its starting location to the destination points. During its motion its path is built online and sent via radio communication module to the followers vehicles as the reference path to follow. Radio system receives the vehicle information in the platoon one by one. Thus, an algorithm is developed in the ICARS framework to receive the information from the radio module and then synchronize them. The lateral controller proposed in [START_REF] Thuilot | Accurate automatic guidance of an urban electric vehicle relying on a kinematic gps sensor[END_REF] is responsible for producing the vehicle steering angle of enable them to follow the reference path of the leader by minimizing both the lateral and angular deviations. Results for the lateral motion are shown in Fig. 3. The longitudinal controller is

B. Normal Platooning Scenario

In this study, the leader travels with variable velocity. The vehicles are initially parked in different positions such that the initial relative distances different from desired ones. Fig. 4 illustrates the results for the leader tracking consensus. The vehicles can achieve the consensus and the Curvilinear position and speed errors converge. Results enlighten the capability of the proposed strategy to create and maintain the platoon in case of a human-driven leader (i.e., the velocity is variable). Table III illustrates the Root Mean Squared Error (RMSE) of the tracking errors which is in the acceptable ranges. Furthermore, a sudden break of the leader occurs

C. Gap Closure Scenario

To investigate the integrated controller capability in case of a gap closing scenario, follower 3 starts at a 32 m gap from its desired set point while the other followers move in normal platoon. As soon as the leader sends commands to the follower 3 to join the platoon, it will move in a gap closing mode, as e e,3 > e l , until it reaches a position where e e,3 <= e l , then it moves in normal platoon.

Fig. 5 shows the numerical results for a gap closing scenario. These results show that the proposed algorithm is able to reduce the gap closing time without violates the acceleration and velocity constraints for both passenger comfortability and safety. To quantify these results, one can use the performance index (24) which ensures that the response is improved by 11.4%.

D. Collision Avoidance Scenario

In this scenario, we test how a vehicle reacts in case of a sudden braking of a front vehicle inside the platoon. A braking maneuver performed by the follower 1 where its speed is sharply reduced at instant 45 s from 5.9 m /s to a full stop. Results in Fig. 6 show how the vehicle performs in case that collision avoidance controller is not used, i.e., u c i = 0. These results illustrate that platoon is failed and collision occurs between vehicles 2 and 1. A unified consensus-based control technique has been proposed to address different issues of platooning in urban environment, including normal platooning, gap closure, and collision avoidance. A second-order longitudinal model of the platoon in the path coordinates is presented. The controller is designed to achieve convergence of the tracking errors taking into consideration the motion in urban environments. An algorithm for adapting the controller parameters is presented to achieve gap closure scenario. Furthermore, a collision avoidance algorithm is designed and integrated with the main controller to guarantee a safety region between vehicles during their travel in the platoon. Conditions to guarantee internal stability and string stability are given. The resulting algorithm was implemented in a realistic simulator called ICARS with real vehicles models. Simulation results have validated the proposed strategy and the theoretical work.

Fig. 1 :

 1 Fig. 1: Kinematics of car-like vehicle.

  ) Design of u p i for Gap Closure Capability: As discussed in the transient analysis, the speed of the response can be designed from the relation, c = b 2ζ 2

Fig. 3 :

 3 Fig. 3: Lateral Control Results; L: Leader Actual Path, and F i : Follower i Actual Path.

Fig. 4 :

 4 Fig. 4: Platoon maintaining scenario.

1 (c) a 3 Fig. 5 :

 135 Fig. 5: Demonstration of Gap Closure Scenario for the 3 rd vehicle: comparison study. (•) g value in case of using the Gap closure algorithm

d 2 Fig. 6 :

 26 Fig. 6: Demonstration of collision avoidance capability: comparison study. (•) c denotes the value in case of using the collision avoidance capability

TABLE I :

 I Vehicles parameters for the realistic study

	Parameter/Vehicle	FLUENCE	ZOE
	Mass	1605 kg	1428 kg
	Width	1.545 m	1.945 m
	Height	1.462 m	1.562 m
	Length	4.748 m	4.084 m
	Rear-axis to Rear-end	1.140 m	0.657 m
	distance between axes	2.701 m	2.588 m
	distance between wheels	1.545 m	1.511 m
	Wheel radius	0.29 m	0.29 m
	Inertia momentum	45.0 kgm 2	28.0 kgm 2
	implemented at a rate of 100 Hz. Control parameters, given
	in Table II, are chosen such that both the internal stability
	and string stability are guaranteed.

TABLE II :

 II Controller Parameters

	Parameter	Value	Parameter	Value
	u ∈	[-6, 1] m /s 2	v ∈	[0, 8] m /s
	d r,i	10 m	γ	0.5
	e l	0.2d r,i	eu	0.8d r,i
	ζ l	0.001	b	1.6
	kc	1.5	ds	5 m

TABLE III

 III

	: RMSE of tracking errors
	RMSE / Follower i	1	2	3
	e s,i [m]	0.2103	0.0872	0.0482
	e q,i [ m /s]	0.0763	0.0297	0.0219

ACKNOWLEDGMENTS

The work presented in this paper is accomplished in the framework of the Valet project (reference: ANR-15-CE22-0013-02).


G Penelet

M Guedra

V Gusev

T Devaux

Simplified account of Rayleigh streaming for the description of nonlinear processes leading to steady state sound in thermoacoustic engines

Keywords: numbers: 43, 35, Ud,43, 25, Nm

This paper focuses on the transient regime of wave amplitude growth and stabilization occuring in a standing wave thermoacoustic engine. Experiments are performed on a simple apparatus consisting of an open ended thermoacoustic oscillator with atmospheric air as working fluid. The results show that, even in that simple device, the transient regime leading to steady state sound exhibits complicated dynamics, like the systematic overshoot of wave amplitude before its final stabilization, and the spontaneous and periodic switch on/off of the thermoacoustic instability at constant heat power supply. A simplified model is then presented which describes wave amplitude growth from the coupled equations describing thermoacoustic amplification and unsteady heat transfer. In this model, the assumption of a one-dimensional and exponential temperature profile is retained and the equations describing heat transfer through the thermoacoustic core are substantially simplified into a set of ordinary differential equations. These equations include the description of two processes saturating wave amplitude growth, i.e. thermoacoustic heat pumping and heat convection by acoustic streaming. It is notably shown that accounting for the effect of acoustic streaming allows to reproduce qualitatively the overshoot process.

I. INTRODUCTION

Thermoacoustic engines belong to a type of heat engines in which the application of a temperature gradient along an open-cell porous medium placed inside an acoustic resonator results in the onset of large amplitude self-sustained acoustic waves. These kind of engines have been studied for about three decades, but some effort is still devoted to the description of their operation. The well-established linear (or weakly nonlinear) theory of thermoacoustics [START_REF] Rott | Thermoacoustics[END_REF][START_REF] Swift | Thermoacoustic engines[END_REF] is largely used for the design of high power (typically up to a few kilowatts) thermoacoustic engines, and for the prediction of their performances with reasonable accuracy [START_REF] Ward | Interactive analysis, design and teaching for thermoacoustics using DeltaEC (A)[END_REF]. However, one can find in the litterature some experimental evidences of complicated effects during the transient regime of wave amplitude growth, which cannot be predicted by a steady-state theory. For example, the existence of a hysteretic loop in the onset and damping of the engine has been reported in both standing wave [START_REF] Chen | Experimental investigation on the onset and damping behavior of the oscillation in a thermoacoustic prime mover[END_REF] or traveling wave devices [START_REF] Jin | Characteristics study on the oscillation onset and damping of a traveling-wave thermoacoustic prime mover[END_REF]; complicated dynamics of the acoustic oscillations have also been reported in various devices, like the periodic switch on/off of thermoacoustic instability [6][START_REF] Zhou | Experimental research of thermoacoustic prime-mover[END_REF][START_REF] Yu | Investigation on the oscillation modes in a thermoacoustic Stirling prime mover: mode stability and mode transition[END_REF][START_REF] Penelet | Nontrivial influence of acoustic streaming on the efficiency of annular thermoacoustic prime movers[END_REF], the "double-threshold effect" [START_REF] Penelet | Experimental investigation of transient nonlinear phenomena in an annular thermoacoustic prime-mover: observation of a double-threshold effect[END_REF] or the "fish-bone like instability" [START_REF] Yu | Fishbone-like instability in a looped-tube thermoacoustic engine[END_REF]. All of these effects indicate that thermoacoustic engines can operate as multistable systems, and also that such complicated dynamical behaviors are due to the existence of different time scales in the process of wave amplitude saturation. That is the reason why the development of adequate simulation tools is still needed to describe the evolution of acoustic wave amplitude after the onset of thermoacoustic instability. On the one hand, direct numerical simulation [START_REF] Hamilton | Nonlinear two-dimensional model for thermoacoustic engines[END_REF][START_REF] Karpov | A nonlinear model of thermoacoustic engines[END_REF][START_REF] Yu | Study of nonlinear processes of a large experimental thermoacoustic-Stirling heat engine by using computational fluid dynamics[END_REF][START_REF] Shimizu | Numerical study of thermoacoustic Taconis oscillations[END_REF] seems to be the only way to reproduce quantitatively the effects mentioned above, but it is still limited by large computation times inherent to the complicated physics and the multiple time and space scales involved in the description of thermoacoustic engines. On the other hand, analytical models are often based on substantial approximations, but one can be motivated by the development of some phenomenological approach aiming at reproducing qualitatively the experiments in order to get deeper understanding of the operation of thermoacoustic engines.

Different approaches have been presented concerning the analytical description of the transient regime. Karpov et al. [START_REF] Karpov | Nonlinear saturation of the thermoacoustic instability[END_REF] proposed a time-domain description of the evolution of thermoacoustic instability combined with a multiple time scales method to calculate the initial wave amplitude growth and its saturation due to higher harmonics generation in the case of a standing wave engine with a fixed temperature gradient. De Waele [START_REF] De Waele | Basic treatment of onset conditions and transient effects in thermoacoustic Stirling engines[END_REF] presented a simplified model based on the lumped element description of thermoacoustic engines, and he performed calculations of the transient regime in a so-called travelling wave thermoacoustic Stirling engine [START_REF] Backhaus | A thermoacoustic Stirling heat engine[END_REF], in which the effect saturating the wave amplitude growth is the cooling effect due to acoustic oscillations in the thermoacoustic core. Penelet et al. [START_REF] Penelet | Nontrivial influence of acoustic streaming on the efficiency of annular thermoacoustic prime movers[END_REF][START_REF] Penelet | Experimental and theoretical study of processes leading to steady-state sound in annular thermoacoustic engines[END_REF] developped a theoretical model of an annular thermoacoustic engine, in which the equation characterizing the variation of acoustic pressure amplitude is combined with the equations describing unsteady heat transfer through the thermoacoustic core. They took into account the influence of various nonlinear effects such as minor losses at the edges of the stack, higher harmonics generation, heat pumping by acoustic waves, and heat convection by the so-called Gedeon streaming [START_REF] Gedeon | DC gas flows in Stirling and pulse tube cryocoolers[END_REF]. It is worth noting that only a few papers [START_REF] Penelet | Nontrivial influence of acoustic streaming on the efficiency of annular thermoacoustic prime movers[END_REF][START_REF] Karpov | A nonlinear model of thermoacoustic engines[END_REF][START_REF] Yu | Study of nonlinear processes of a large experimental thermoacoustic-Stirling heat engine by using computational fluid dynamics[END_REF][START_REF] Penelet | Experimental and theoretical study of processes leading to steady-state sound in annular thermoacoustic engines[END_REF] provide direct comparisons of calculated transient regimes with experimental data, and to our knowledge, most of the complicated effects mentioned above cannot be reproduced by the models, even for the simplest thermoacoustic oscillators which can be built, like the Sondhaus tube or the so-called thermoacoustic laser [START_REF] Garrett | Build an Acoustic Laser[END_REF].

In this paper, our objective is to provide a simplified model of a standing wave thermoacoustic oscillator, keeping in mind our objective to capture some of the complicated dynamical behaviors observed in experiments. The main approximations retained in the model concern the description of heat transfer through the thermoacoustic core, but some important effects involving the mutual influence of acoustic and temperature fields are taken into account. In particular, the heat convected by Rayleigh streaming can be included (with great simplification) in the model. In section II the experimental apparatus is described, and typical transient regimes of wave amplitude growth are presented. In section III, the theoretical model is presented, while section IV provides comparisons between experiments and theory.

II. EXPERIMENTS

The system under study is a basic thermoacoustic standing wave engine, which is quite easy to build [START_REF] Garrett | Build an Acoustic Laser[END_REF] and often used as a demonstration apparatus for graduate students. A photograph of this thermoacoustic oscillator is presented in A schematic drawing of the apparatus, comprising the definition of an appropriate system of coordinates, is presented in Fig. 2. The first step in studying this device consists in determining the heat power supply Q onset which is necessary to initiate self-sustained acoustic waves. This critical value Q onset of heat supply depends on the position x s of the stack along the waveguide, and the corresponding stability curve is presented in Fig. 3. The experimental protocol used to obtain Q onset -or more precisely the range of heat power Q comprising Q onset -in function of x s is as follows: (1) fix a position x s of the stack, without heating ; (2) apply an electrical current increment ∆I = 0.1A; (3) wait for 10 minutes (stabilization of the thermal field); (4) repeat steps (2) and (3) until the onset of thermoacoustic instability. As shown in Fig. 3, there exists an optimum position of the stack, x s ≈ 35cm, corresponding to the lower Q onset . Note that a theoretical stability curve is also presented in Fig. 3, which will be discussed in section IV.

Once the stability curve is determined, it is then possible to measure the evolution of wave amplitude growth. For each of the measurements presented in the following, the heat power supply is preliminarily set to a value Q 0 slightly below Q onset . A small ∆Q increment on heat power supply is then sufficient for the acoustic wave to be generated in the device at a frequency f ≈ c 0 /(4L) (i.e. the frequency of the first mode of the empty resonator). Typical transient regimes of wave amplitude growth are presented in Figs. 4 and5, for various increments of heat supply ∆Q, and for two positions of the stack along the waveguide. Figure 4 shows the transient and steady state waveform of acoustic pressure p(L, t), for different values of ∆Q, and when the stack is placed at position x s = 36.5 cm (i.e. next to the closed end of the resonator). The results show that there exists a systematic overshoot of wave amplitude growth before its final stabilization, and that the higher the increment ∆Q is, the faster will be reached stabilization. Figure 5 shows the measured transient regimes when the stack is placed at position x s = 26.5 cm. In that case, the evolution with time of acoustic pressure amplitude exhibits more complicated dynamics: at moderate increments ∆Q of heat power above its initial value Q 0 , the system turns on and off spontaneously and almost periodically. Also, the switch on/off period decreases when ∆Q increases, so that after some critical value ∆Q cr of the heat increment, the acoustic wave finally stabilizes to a finite value (0.16Q 0 < ∆Q cr ≤ 0.24Q 0 if x s = 26.5 cm, as depicted in Fig. 5). As mentioned before, similar switch on/off processes have already been reported in the litterature concerning both standing wave [6,[START_REF] Zhou | Experimental research of thermoacoustic prime-mover[END_REF] or travelling wave engines [START_REF] Yu | Investigation on the oscillation modes in a thermoacoustic Stirling prime mover: mode stability and mode transition[END_REF][START_REF] Penelet | Nontrivial influence of acoustic streaming on the efficiency of annular thermoacoustic prime movers[END_REF]. Moreover, some of these papers [START_REF] Zhou | Experimental research of thermoacoustic prime-mover[END_REF][START_REF] Yu | Investigation on the oscillation modes in a thermoacoustic Stirling prime mover: mode stability and mode transition[END_REF][START_REF] Penelet | Nontrivial influence of acoustic streaming on the efficiency of annular thermoacoustic prime movers[END_REF] also report that the evolution of acoustic pressure amplitude is accompanied by significant variations of the temperature field in the stack and in the thermal buffer tube (i.e. the region of the waveguide where the temperature field is heterogeneous). So, it seems to be clearly admitted that such an effect is mostly due to the reverse influence of the acoustic field on the temperature field (heat pumping by acoustic waves, acoustic streaming). We performed several measurements of the transient regimes for different positions of the stack, and we found that the closer is the stack from the open end of the resonator, the higher is the critical heat power increment ∆Q cr . Moreover, we did not observe any switch on/off process when the stack position x s is higher than x s = 31.5 cm (note however that this process might have happened if we had used lower ∆Q increments).

The observation of the switch on/off process thus depends on the stack position along the waveguide.

III. THEORY A. Thermoacoustic amplification

In the device under consideration, the frequency and amplification rate of self-sustained acoustic waves depend on the geometry of the device and on the spatial distribution of the temperature field. If the instantaneous temperature field T (x, t) is known, one can compute both the instantaneous thermoacoustic amplification rate ǫ ampl [T (x, t)] and the corresponding onset angular frequency Ω [T (x, t)], provided that the assumption of a quasisteady state, ǫ ampl << Ω, can be retained. In this paper, we derive the essential steps of the calculations without details: this problem has been treated recently [START_REF] Guedra | On the use of a complex frequency for the description of thermoacoustic engines[END_REF] in a more general situation. Assuming that harmonic plane waves are propagating along the duct, the acoustic pressure p(x, t) and acoustic volume velocity u(x, t) are written in the following form:

ξ(x, t) = ℜ ξ(x)e -jωt , (1)
where j 2 = -1, ξ may be either p or u, ξ denotes the complex amplitude of ξ, ℜ () denotes the real part of a complex number, and ω stands for the angular frequency of acoustic oscillations.

The thermoacoustic device is then separated into two parts, i.e. the thermoacoustic core comprising both the stack and the waveguide region where temperature is inhomogeneous (x ∈ [x s -d s , L], see Fig. 2), and the cold part of the waveguide (x ∈ [0, x s -d s ]) where the fluid is assumed to be at room temperature T ∞ . The propagation of acoustic waves through the thermoacoustic core can be described on the whole by an acoustical two-port relating the complex amplitudes of acoustic pressure and volume velocity at both sides:

  p(L) ũ(L)   =   ∆ pp ∆ pu ∆ up ∆ uu   ×   p(x s -d s ) ũ(x s -d s )   , (2)
where the coefficients ∆ ij depend on the geometrical and thermophysical properties of the components, on the temperature distribution T (x) and on the angular frequency ω. The determination of the coefficients ∆ ij requires to describe acoustic propagation in both the thermal buffer tube (x ∈ [x s , L]) and the stack (x ∈ [x s -d s , x s ]), in the presence of an inhomogeneous temperature field. This problem has been addressed by Rott [START_REF] Rott | Thermoacoustics[END_REF] who established the wave equation associated to this problem (e.g. Eq. (54) in ref. [START_REF] Swift | Thermoacoustic engines[END_REF]). This wave equation is here transformed into an integral Volterra equation of the second kind [START_REF] Gusev | Asymptotic theory of nonlinear acoustic waves in a thermoacoustic prime mover[END_REF],

which allows to compute the coefficients ∆ ij (ω) for a fixed temperature distribution T (x).

The details of derivation are not provided in this paper due to a need of conciseness, but the reader can report to ref. [START_REF] Penelet | Dependence of sound amplification on temperature distribution in annular thermoacoustic engines[END_REF], where both the solving process and the explicit expression of the transfer matrix (Eq. 19 in ref [START_REF] Penelet | Dependence of sound amplification on temperature distribution in annular thermoacoustic engines[END_REF]) are provided. Note also that in the remaining of the paper, the stack which actually consists of multiple square-channels of geometrical radius a/2 is assumed equivalent to a stack of multiple cylindrical channels with the same radius r s = a/2 (the cylindrical channel is formally more easy to treat than the square channel).

Then, neglecting sound radiation at the open end of the resonator (p(x = 0) = 0 ), the lossy propagation of acoustic waves in the remaining of the waveguide (x

∈ [0, x s -d s ]
) is taken into account by the reflected impedance

Z s = p(x s -d s ) ũ(x s -d s ) = jZ c tan (k (x s -d s )) , (3)
where k is the complex wavenumber defined as

k = ω c 0 1 + (γ -1) f κ 1 -f ν 1/2 , (4)
and where Z c is the characteristic impedance of the lossy duct, defined as

Z c = ρ f c 0 πr 2 i [(1 + (γ -1) f κ ) (1 -f ν )] -1/2 . (5)
In Eqs. ( 4) and ( 5), ρ f and γ stand for the density and the specific heat ratio of fluid, respectively, and the functions f ν and f κ , defined as

f ν,κ = 2r i (1 + i)δ ν,κ J 1 (1 + i) δν,κ r i J 0 (1 + i) δν,κ r i (6)
are the well-known functions [START_REF] Swift | Thermoacoustic engines[END_REF][START_REF] Arnott | General formulation of thermoacoustics for stacks having arbitrarily shaped pore cross sections[END_REF] which characterize the viscous and thermal coupling between the oscillating fluid and the waveguide walls (J n is the n th order Bessel function of the first kind; δ ν = 2µ f /(ρ f ω) and δ κ = 2κ f /ω are the viscous and thermal boundary layer thicknesses, respectively ; µ f and κ f stand for dynamic viscosity and thermal diffusivity of fluid at room temperature T ∞ , respectively). Combining Eqs. ( 2) and ( 3) with the boundary condition ũ(L) = 0 (rigid termination) leads to the characteristic equation

G [ω, T (x)] = Z s ∆ up + ∆ uu = 0, (7)
which must be satisfied. Because Eqs. ( 2) and ( 3) are derived in the frequency domain, Eq. ( 7) should only be satisfied at threshold of thermoacoustic instability (i.e. the only situation for which a steady-state solution is possible for the acoustic field, in the absence of nonlinear saturating processes), but as it is discussed in detail in ref. [START_REF] Guedra | On the use of a complex frequency for the description of thermoacoustic engines[END_REF], one may use Eq. ( 7) to solve the unsteady problem of wave amplitude growth. This is realized by allowing the angular frequency ω to have an imaginary part ǫ ampl , i.e. ω = Ω + jǫ ampl , which means that acoustic pressure is being amplified (ǫ ampl > 0) or attenuated (ǫ ampl < 0). This assumption enables to satisfy Eq. ( 7) even if the instantaneous temperature field T (x, t) does not allow to satisfy the marginal stability condition. Therefore, if the instantaneous temperature field T (x, t) is known at time t, Eq. ( 7) can be solved using conventional numerical methods in order to obtain both the instantaneous thermoacoustic amplification rate ǫ ampl [T (x, t)] and the corresponding angular frequency Ω [T (x, t)] of acoustic oscillations [START_REF] Guedra | On the use of a complex frequency for the description of thermoacoustic engines[END_REF]. Note however that the solving approach described above requires to assume that the system is under quasisteady state, which means that at the time scale of an acoustic period, the variations of the temperature field are negligible while those of the wave amplitude are not significant (ǫ ampl << Ω). Finally, once the instantaneous thermoacoustic amplification rate ǫ ampl is known, it is quite direct from Eq. ( 1) that the time variations of the peak amplitude of acoustic pressure, and in particular the peak pressure P (t) at position x = L where the microphone is flush-mounted, are described with the following ordinary differential equation:

d P dt = ǫ ampl P. (8)

B. Unsteady heat transfer

The main simplification of the model presented in this paper concerns the description of unsteady heat transfer through the thermoacoustic core. Actually, this description of heat transfer is so much simplified in the following that one cannot expect, at best, anything else than qualitative agreement between experiments and theory. It is first considered here that the glass tube is a perfect thermal reservoir: its temperature is assumed to be constant at room temperature T ∞ (even if heat is absorbed from the inside of the tube). It is also considered that the thermophysical properties (density, heat capacity, thermal conductivity) of the stack walls and of the fluid do not depend on temperature (except in section III D 1

where the velocity of acoustic streaming is estimated). Moreover, the stack, which consist of both solid walls and air, is treated as a homogeneous medium of density ρ s , thermal conductivity λ s and heat capacity C s (these average thermophysical properties are calculated from those of both ceramic and air, see Tab.I). Finally, the temperature field inside the tube is assumed to be axisymetric, and the equations describing heat transfer in the thermoacoustic core are written as follows:

x ≤ x s , ∂T (x, t) ∂t = κ s ∂ 2 T (x, t) ∂x 2 - T (x, t) -T ∞ τ s , ( 9a
) x ≥ x s , ∂T (x, t) ∂t = κ f ∂ 2 T (x, t) ∂x 2 - T (x, t) -T ∞ τ f . (9b)
In Eqs. (9), T (x, t) refers to the unsteady cross-sectional average temperature inside the tube at position x, and κ s,f stand for the thermal diffusivities of the stack (subscript "s") and of the fluid in the waveguide (subscript "f"). The last terms on the right-hand-side of Eqs. ( 9) describe heat transfer between the inside of the duct and its surrounding walls.

The phenomenological coefficients τ f,s are obtained from empirical correlations which can be found in heat transfer textbooks [START_REF] Incropera | Fundamentals of Heat and Mass Transfer[END_REF]. While obtaining the analytical expression of the parameter τ f is quite direct, the derivation of τ s is less direct, but the details are given in the appendix of ref. [START_REF] Penelet | Experimental and theoretical study of processes leading to steady-state sound in annular thermoacoustic engines[END_REF]. The numerical values of the thermophysical parameters mentioned above are summarized in Tab.I. Then, the external thermal action due to the heat power Q(t) dissipated in the Nichrome wire is taken into account through the following boundary condition:

λ s ∂T (x - s , t) ∂x -λ f ∂T (x + s , t) ∂x = Q(t) πr 2 i , (10)
describing the continuity of heat flux at the hot side of the stack, where λ s and λ f stand for the thermal conductivities of the stack and the fluid, respectively. In addition to the simplifications mentioned above, it is also assumed that the spatial distribution of the thermal field is exponential. More precisely, introducing the space variable z = x -x s , we seek a solution in the form:

z ≤ 0, T (z, t) ≈ T ∞ + ∆T (t)e z/ls(t) , (11a)
z ≥ 0, T (z, t) ≈ T ∞ + ∆T (t)e -z/l f (t) , (11b)
where ∆T (t) = T (z = 0, t) -T ∞ , and where the variables l s (t) and l f (t) are typical lengths which are representative of the actual temperature distribution and its variations with time.

Then, introducing these solutions in Eqs. ( 9) leads to the following equations:

d ∆T dt - z l 2 s ∆T d l s dt = κ s ∆T l 2 s - ∆T τ s , (12a)
d ∆T dt + z l 2 f ∆T d l f dt = κ f ∆T l 2 f - ∆T τ f . (12b)
Therefore, the assumption of a time varying exponential temperature profile in Eqs. [START_REF] Yu | Fishbone-like instability in a looped-tube thermoacoustic engine[END_REF] implies that the parameters characterizing the spatial distribution of temperature (l s , l f , ∆T ) do not depend on the space variable z, but the the problem still depends on z in Eqs. [START_REF] Hamilton | Nonlinear two-dimensional model for thermoacoustic engines[END_REF]. Carrying on in our will to give priority to simplicity at the expense of rigour, the space variable z is thus replaced in Eqs. (12a) and (12b) by its characteristic values -l s (t) and l f (t), respectively. This means that it is considered that the temperature variations at positions z = -l s and z = l f are representative of the global temperature field. Finally, accounting for this new assumption, and reporting also Eqs. (11a) and (11b) in Eq. ( 10), the description of unsteady heat transfer reduces to the following set of equations:

l s ∆T d ∆T dt + d l s dt = κ s l s - l s τ s (13a) l f ∆T d ∆T dt + d l f dt = κ f l f - l f τ f (13b) λ s l s + λ f l f ∆T = Q πr 2 i . ( 13c
)
The combination of Eqs. (13a -13c) with Eq. ( 8) allows to describe the initial start-up of self-sustained acoustic waves when the device is submitted to the external heat input Q(t). However, the equations mentioned above would describe the indefinite growth of wave amplitude (as soon as Q > Q onset ) if the nonlinear effects saturating the amplification process are not taken into account. In the following, two processes are considered to describe the stabilization of acoustic amplitude, i.e. the thermoacoustic heat flux due to acoustic oscillations along the stack, and the heat convection by acoustic streaming. Note that there exists additional mechanisms which may play a role in the saturation of wave amplitude growth, like the dissipation of acoustic energy due to geometrical singularities at the edges of the stack and at the open end of the resonator [START_REF] Peters | Acoustically induced periodic vortex shedding at sharp edged open channel ends: simple vortex models[END_REF], or like nonlinear propagation in the openended tube [START_REF] Chester | Resonant oscillations of a gas in an open-ended tube[END_REF]. Those processes are voluntarily ignored in this study, because attention is focused on the effects involving heat transport induced by high amplitude acoustic waves.

ρ f C f = 1.2 × 1003 Jm -3 K -1 λ f = 2.26 10 -2 W m -1 K -1 µ f = 1.84 10 -5 P a s Cordierite ρ c C c = 2600 × 1465 Jm -3 K -1 λ c = 3 W m -1 K -1 - Stack ρ s C s = Φρ f C f + (1 -Φ)ρ c C c λ s = Φλ f + (1 -Φ)λ c -

Heat exchange with resonator walls

Fluid/resonator † :

τ f = r 2 i ρ f C f 3.66λ f τ f = 8.29 s Stack/resonator ‡ : τ s = ρsCsr 2 i ln(2) 2 e λc(2rs+e) + 2rs 29.28λ f rs+λce τ s = 288.2 s † see Eq. (A.
3) in ref. [19] ‡ see Eq. (A.12) in ref. [START_REF] Penelet | Experimental and theoretical study of processes leading to steady-state sound in annular thermoacoustic engines[END_REF] TABLE I. Geometrical and thermophysical properties of the elements constituting the apparatus.

Note that the thermophysical properties of the fluid and the stack material are given at 300 K.

These parameters are assumed constant in the model, except in section III D 1 for ρ f , λ f and µ f when estimating the velocity of acoustic streaming (ρ f ∝ T -1 , while λ f ∝ T β and µ f ∝ T β , with β = 0.77 [START_REF] Bailliet | Acoustic streaming in closed thermoacoustic devices[END_REF]).

C. Thermoacoustic heat flux.

In any kind of thermodynamic heat engine, the production of work induces heat transfer from a hot source to a cold sink. This is also the case in thermoacoustic engines, in which the production of acoustic work is accompanied by acoustically induced heat transport [START_REF] Swift | Thermoacoustic engines[END_REF], with subsequent reduction of the externally imposed temperature gradient. Under the circumstance where a steep temperature gradient is applied along the stack, it is worth noting that this thermoacoustic heat flux is approximately proportional to the temperature gradient d x T , so that, in the frame of the short stack approximation (d s << 4L), the heat pumping by acoustic waves can be described by an effective thermal conductivity [ref. [START_REF] Yazaki | Large heat transport due to spontaneous gas oscillation induced in a tube with steep temperature gradients[END_REF], Eq. ( 9)] as follows:

λ ac = ρ f C f ω ℑ σf (s) ν -f (s) κ σ 2 -1 < u 2 (x s , t) > πr 2 i , (14)
where C f and σ stand for the isobaric heat capacity and the Prandtl number of fluid, < • • • > denotes time averaging over an acoustic period, and where ℑ(. . . ) stands for the imaginary of a complex number. The functions f

(s) ν and f (s)
κ are the functions characterizing viscous and thermal coupling between the oscillating fluid and the stack walls and they are obtained from Eq. ( 6) in which r i is replaced by the pore radius r s . The expression of the acoustically enhanced thermal conductivity λ ac can be further simplified by assuming that the acoustic field roughly corresponds to a pure standing wave oscillating at frequency f 0 = c 0 /(4L), so that the spatial distribution u(x, t) of acoustic volume velocity can be written

u(x, t) ≈ πr 2 i P (t) ρ f c 0 cos(k 0 x) sin(Ω 0 t), (15)
where k 0 = π/(2L), Ω 0 = πc 0 /(2L), and where P (t) refers to the (slowly varying) peak amplitude of acoustic pressure at the closed end of the resonator (note that the actual spatial distribution of the acoustic field is not that of a perfect standing wave, because of thermo-viscous losses and sound scattering by the thermoacoustic core). This leads to the following expression

λ ac ≈ 1 2 C f ρ f Ω 0 c 2 0 ℑ σf (s) ν -f (s) κ σ 2 -1 cos 2 (k 0 x s )P 2 = Γ λ P 2 (16)
of the acoustically enhanced thermal conductivity, which can be included in Eqs. (13c) and (13a) in the form of an additional term. Finally, the equations describing unsteady heat transfer, thermoacoustic amplification and heat pumping by acoustic waves are written as follows:

l s ∆T d ∆T dt + d l s dt = κ s + Γ κ P 2 l s - l s τ s , (17a)
l f ∆T d ∆T dt + d l f dt = κ f l f - l f τ f , (17b)
λ s + Γ λ P 2 l s + λ f l f ∆T = Q πr 2 i , (17c)
d P dt = ǫ ampl P, (17d)
where Γ κ = Γ λ /(ρ s C s ), and where ρ s C s stand for the average volumetric heat capacity of the stack (see. Tab I).

This set of equations can be transformed into a system of ordinary differential equations describing the evolutions of the parameters ∆T , l s , l f and P for appropriate initial conditions and for some external action Q(t) on the system. It is however possible to include the additional process of heat convection by acoustic streaming, as will be discussed in the following.

D. Acoustic streaming

Acoustic streaming refers to the steady mass flow induced by large amplitude acoustic oscillations [START_REF] Boluriaan | Acoustic streaming: from Rayleigh to today[END_REF]. This unavoidable nonlinear effect may be due for instance to Reynolds stresses inside viscous boundary layers, or to the difference in minor loss phenomenon in a differentially heated stack [START_REF] Mironov | Acoustic streaming related to minor loss phenomenon in differentially heated elements of thermoacoustic devices[END_REF]. Acoustic streaming is discarded in many applications involving acoustics, but it is now well-established that it must be considered in thermoacoustics, because it convects heat and tends to reduce the temperature gradient imposed along the stack. When considering acoustic streaming excitation in thermoacoustic engines, it is quite usual to dissociate two classes of streaming (though the two classes may owe their origin from the same nonlinear effects). The first class of streaming, i.e. the so-called Gedeon streaming [START_REF] Gedeon | DC gas flows in Stirling and pulse tube cryocoolers[END_REF], is driven in thermoacoustic Stirling engines containing a closed loop path for the acoustic wave propagation: this is a directional streaming accompanying wave propagation and carrying a non-zero mass flow through each cross-section of the loop. In most of the one dimensional models describing thermoacoustic engines, it is possible to account for the heat convection due to Gedeon streaming (e.g. in the freely available software Delta-EC [START_REF] Ward | Interactive analysis, design and teaching for thermoacoustics using DeltaEC (A)[END_REF] or in ref. [START_REF] Penelet | Nontrivial influence of acoustic streaming on the efficiency of annular thermoacoustic prime movers[END_REF]) and, in practical engines, this directional mass flow can be cancelled using a membrane or a jet pump. The second class of streaming, i.e. the Rayleigh streaming, is the classical boundary layer driven streaming excited at large acoustic amplitudes in resonant gas columns: unlike Gedeon streaming, it is characterized by large vortex cells which do not carry any mass flow across the resonator's cross-section (i.e. the mean mass flow rate is zero).

The spatial distribution of Rayleigh streaming can be obtained from classical theories using perturbation expansions (e.g. in ref. [START_REF] Bailliet | Acoustic streaming in closed thermoacoustic devices[END_REF]) but because it is by nature three-dimensional, it is difficult to quantify the heat it convects when the basic modeling of the thermoacoustic device is one-dimensional. Moreover, it is difficult to cancel this streaming, though one can diminish its amplitude by tapering the waveguide [START_REF] Olson | Acoustic streaming in pulse tube refrigerators: tapered pulse tubes[END_REF]. In the following, a simplified approach is presented which attempts to estimate the influence of Rayleigh streaming in the operation of a standing wave thermoacoustic engine.

Spatial distribution of acoustic streaming at threshold

Before trying to consider acoustic streaming in the equations governing heat transfer in the thermoacoustic core, it is useful to evaluate its magnitude and its spatial distribution in the present device. This can be done from the analytical model developped by Bailliet et al. [START_REF] Bailliet | Acoustic streaming in closed thermoacoustic devices[END_REF], which allows to calculate both the velocity of acoustic streaming and the second order mass flux in closed (as opposed to closed-loop) thermoacoustic devices, in the frame of a weakly nonlinear theory where successive approximations methods can be used. To do this, we use the protocol described in the following. At first, the geometrical and thermophysical parameters of the thermoacoustic device are fixed according to Tab. I. Then, the steady state temperature field T (x) is computed for a fixed heat power supply Q 0 in the absence of sound :

this is quite direct by setting P = 0 and d dt = 0 in Eqs. [START_REF] De Waele | Basic treatment of onset conditions and transient effects in thermoacoustic Stirling engines[END_REF], which leads to l s = l s 0 = √ κ s τ s ,

l f = l f 0 = √ κ f τ f and ∆T = ∆T 0 = (Q 0 l f 0 l s 0 )/ [(πr 2 i )(λ s l f 0 + λ f l s 0 )].
Then, assuming that the peak pressure amplitude P at the closed end of the resonator equals 1 Pa, the equations describing linear acoustic propagation through the device are solved to obtain the spatial distribution of the first order acoustic variables (pressure, density, temperature, axial and transverse particle velocities), which themselves are used as input parameters to compute the second order time-averaged streaming mass flow ṁ. More precisely, the spatial distribution of the axial streaming velocity in a cylindrical tube (the resonator or one stack channel) is computed from Eq. ( 16) in ref. [START_REF] Bailliet | Acoustic streaming in closed thermoacoustic devices[END_REF], and the resulting spatial distribution ṁ(x, r) of the rate at which mass flows along the axial direction is obtained using Eq. ( 11) in the same reference. Note that in the following evaluations of ṁ, the heat power supply Q 0 is fixed to the x s -dependent value Q onset (x s ) in order that the calculated thermoacoustic amplification coefficient ǫ ampl (∆T, l s , l f ) equals zero (this corresponds to the threshold of thermoacoustic instability for the considered position x s of the stack). The calculated streaming mass flow ṁ(x, r) is presented in Fig. 6, when x s is fixed to x s = 26.5 cm and Q 0 to Q 0 = Q onset (x s ) = 20.17 W : it is presented in function of both axial and transverse coordinates, in the stack region [Fig. 6(a)] and in the waveguide region [Fig. 6(b)]. The results show that the average magnitude of ṁ is of the same order of magnitude in both the stack and the waveguide, i.e. around 10 -8 kg.m -2 .s -1 when P = 1 P a. The results also show that the mass flux along the centerline of one stack pore goes rightwards while that along the centerline of the waveguide goes leftwards. Note also that in both cases the direction of mass flow reverses around

r ≈ r i,s / √ 2.
It is however worth noting that the calculated second order axial mass flow is not realistic near the stack/waveguide interface, because the stack pore or the waveguide are both treated as isolated systems. In other words, there exists complicated hydrodynamical edges effects due to the geometrical singularity at position x = x s which impact both the oscillating and the steady flows, but which are not taken into account here. Recent LDV measurements performed by Moreau et al. [START_REF] Moreau | Effect of a stack on Rayleigh streaming cells investigated by laser Doppler velocimetry for application to thermoacoustic devices[END_REF] indicate that the presence of a stack in the waveguide induces new streaming vortices in the vicinity of the stack ends. This means that heat transport by acoustic streaming in the vicinity of the stack should be very complex. In that context, the authors easily admit that the derivation presented in the following, which attempts to account for heat transport by acoustic streaming in the standing wave engine, should be considered cautiously.

Simplistic account of acoustic streaming

a. Separation of inner and outer zones. As depicted in Fig. 6, the velocity of acoustic streaming is multidimensional. Therefore, it seems arduous to account for the heat convected by acoustic streaming in the one-dimensional model described in section III B, unless putting up with a simplistic approach. In our attempt to include the effect of acoustic streaming, the actual spatial distribution of the streaming velocity in the waveguide needs to be significantly simplified, and a schematic drawing of the retained assumptions is presented in Fig. 7. In order to reproduce the vortex cell structure of acoustic streaming without considering the details of its spatial distribution, we separate the waveguide cross-section between an inner zone (r ∈ [0, r i / √ 2]) where the streaming mass flow is directed leftwards, and an outer zone (r ∈ [r i / √ 2, r i ]) where it is directed rightwards [Fig. 7(b)]. The variations with axial coordinate of the streaming velocity are also discarded: only the "characteristic" values at positions z = -l s 0 and z = l w 0 are considered (this assumption is consistent with the one used in the transition from Eqs. [START_REF] Hamilton | Nonlinear two-dimensional model for thermoacoustic engines[END_REF] to Eqs. ( 13)). Moreover, the transverse variations of ṁ are ignored: the assumption of an uniform flow is retained in both inner and outer zones and the corresponding mass flow rates are calculated from their cross-sectional averaged values

ṁ(i) f = 2×2π πr 2 i r i / √ 2 0 ṁ(l w 0 , r).rdr and ṁ(o) f = 2×2π πr 2 i r i r i / √
2 ṁ(l w 0 , r).rdr. The same assumptions are also retained inside the stack pores of radius r s : each pore is separated into an inner zone where the streaming mass flow rate ṁ(i) s is directed rightwards and an outer zone where the streaming mass flow rate ṁ(o) s is directed leftwards. Therefore, since the total mass flux across the section of the waveguide equals zero, and since the surface area of the inner zone equals that of the corresponding outer zone ( πr 2 i,s /2 = π(r 2 i,s -r 2 i,s /2) ), the counterflowing mass flow rates have the same magnitude:

ṁ(o) f = - ṁ(i) f = ṁf and ṁ(i) s = - ṁ (o)
str P 2 , ( 18b
)
where ṁ(l w 0 , r) and ṁ(l s 0 , r) are the transverse distribution of the mass flow rate at positions

x s + l w 0 and x s -l s 0 calculated from the theoretical model of Bailliet et al. [START_REF] Bailliet | Acoustic streaming in closed thermoacoustic devices[END_REF]. In Eqs. [START_REF] Backhaus | A thermoacoustic Stirling heat engine[END_REF] the quadratic dependance of ṁf,s with P is obvious because it is calculated from the time average of the products of first order (acoustic) quantities, and the parameters Γ

(s)
str and Γ (f ) str , which depend on the stack position x s and on the temperature distribution, simply correspond to the results of calculation for ṁf and ṁs when a peak amplitude of acoustic pressure P = 1 P a is assigned at position x = L. In the following, these parameters are computed for various stack positions with their corresponding temperature distribution at threshold. In particular, at position x s = 26.5cm, we get Γ (s) str ≈ 1.2 10 -8 kg 3 m -4 s -5 and Γ (f ) str ≈ 4.4 10 -9 kg 3 m -4 s -5 , and at position x s = 36.5cm (with the corresponding Q onset ≈ 14.47 W ) we get Γ (s) str ≈ 5.6 10 -9 kg 3 m -4 s -5 and Γ (f ) str ≈ 2.5 10 -9 kg 3 m -4 s -5 . Another important point concerns the time delay for streaming establishment, which may be important when investigating the dynamics of wave amplitude saturation. It is worth noting that the estimate of streaming velocity mentioned above is valid in steady state, while the stabilization of acoustic streaming after switching on the acoustic field is not instantaneous. The rough estimate θ = ρ f D 2 /(µ f π 2 ) of the characteristic time of streaming establishment was given by Amari et al. [START_REF] Amari | Transient unidirectional acoustic streaming in annular resonators[END_REF], where µ f is the dynamic viscosity of fluid and D is the diameter of the channel. Then one can evaluate the characteristic times θ s,f of streaming establishment in one stack pore and in the waveguide, respectively, which leads to θ s ≈ 5 10 -3 s and θ f ≈ 17.9 s. Clearly, the time delay for streaming establishment can be neglected in the stack, but should be taken into account in the waveguide.

From the assumptions described above, it is finally possible to include acoustic streaming in the equations describing unsteady heat transfer in the thermoacoustic core, Eqs. [START_REF] Penelet | Nontrivial influence of acoustic streaming on the efficiency of annular thermoacoustic prime movers[END_REF] , which leads to:

x ≥ x s , ∂T (i) ∂t - ṁf ρ f ∂T (i) ∂x = κ f ∂ 2 T (i) ∂x 2 - T (i) τ f (19a) x ≥ x s , ∂T (o) ∂t + ṁf ρ f ∂T (o) ∂x = κ f ∂ 2 T (o) ∂x 2 - T (o) τ f (19b) x ≤ x s , ∂T (i) ∂t + Φρ f C f ρ s C s ṁs ρ f ∂T (i) ∂x = κ s + Γ κ P 2 ∂ 2 T (i) ∂x 2 - T (i) τ s (19c) x ≤ x s , ∂T (o) ∂t - Φρ f C f ρ s C s ṁs ρ f ∂T (o) ∂x = κ s + Γ κ P 2 ∂ 2 T (o) ∂x 2 - T (o) τ s (19d) ṁs = Γ (s) str P 2 , (19e)
d ṁf dt + ṁf θ f = Γ (f ) str P 2 θ f (19f)
where T (i) (x, t) and T (o) (x, t) refer to the temperature in the inner zone and the outer zone, respectively, and where the parameter Φρ f C f /(ρ s C s ) in Eqs. (19c) and (19d) accounts for the fact that the mass flow velocity in the stack is an effective velocity which results from the only motion of the fluid in the honeycombed ceramic of porosity Φ and volumetric heat capacity ρ s C s . It should be emphasized, however, that besides their intrinsic limitations associated to the one dimensional approximation in the description of heat transfer, Eqs.

(19) are also rigorously uncorrect for two reasons: firstly, the thermal coupling between the inner and outer zones is ignored, and secondly the description of heat transfer in the inner zones should not include the term T (i) /τ s,f because there is no physical contact between the fluid and solid walls. So, Eqs. [START_REF] Penelet | Experimental and theoretical study of processes leading to steady-state sound in annular thermoacoustic engines[END_REF] are imperfect but at least one can check that setting Γ (s,f ) str = 0 yields Eqs. ( 9) which are valid in the absence of streaming.

b. Estimate of the heat taken by acoustic streaming at the stack/waveguide interface

The equations of unsteady heat transfer, Eqs. [START_REF] Penelet | Experimental and theoretical study of processes leading to steady-state sound in annular thermoacoustic engines[END_REF], must be completed by some boundary condition at the stack/waveguide interface. To do this, Eq. ( 10) is modified as follows:

λ s + Γ λ P 2 ∂ x T | x - s -λ f ∂ x T | x + s + h (f ) conv (T H -T ∞ ) + h (s) conv (T H -T ∞ ) = Q(t) πr 2 i , (20)
where

T H = T (o) (x s , t) = T (i) (x s , t) and ∂ x T | x ± s = 1 2 ∂ x T (i) | x ± s + ∂ x T (i) | x ± s .
In this equation, the heat exchange coefficients h (s) conv and h (f ) conv are representative of some heat taken away from the electrical heat resistance, by acoustic streaming, towards the stack and the waveguide, respectively. Calculating h

(s,f )
conv is actually a difficult task which implies knowing the details of both temperature and streaming velocity fields in the vicinity of the stack/waveguide interface, so that one have to estimate them very roughly. To do this, we consider that the actual problem, i.e. the presence of vortex streaming cells in both stack and waveguide, is equivalent to two independant problems where each cell is unwrapped. This is illustrated by the two schematic drawings of Figs. 7 (c) and (d) in which one focuses on the streaming cell in the waveguide region (x ≥ x s ): it is assumed that the fluid going leftwards at velocity ṁf /ρ f , upwards and downwards along the stack interface, and rightwards away from the stack at velocity ṁf /ρ f , carries the same amount of heat than the fluid which would go rightwards in a single tube of infinite length and which would cross an isothermal grid with fixed temperature T H . In other words, the heat flux ϕ conv taken by convection from the electrical heat resistance at temperature T H , Fig. 7 (c), is assumed to be equal to the heat flux taken by the fluid in the more tractable situation of Fig. 7 (d). The problem of heat transport in that last situation is indeed quite easy to solve: assuming first adiabatic condition along the tube walls, and introducing the new space coordinate y defined in Fig. 7 (d), one can write the following steady-state heat transfer equation

ṁf ρ f ∂T ∂y = κ f ∂ 2 T ∂y 2 . ( 21
)
Then, assuming that the temperature of the ingoing flow T (y = -∞) equals T ∞ , one gets

y ≤ 0, T (y) = (T H -T ∞ ) e ( ṁf y)/(ρ f κ f ) + T ∞ , (22a)
y ≥ 0, T (y) = T H . (22b)
Therefore, one can report Eqs. (22a) and (22b) in the following boundary condition

λ f ∂T (0 + ) ∂y -λ f ∂T (0 -) ∂y + ϕ conv = 0, (23)
where ϕ

(f ) conv is the heat power provided by the grid to sustain the steady-state temperature field of Eqs. [START_REF] Guedra | On the use of a complex frequency for the description of thermoacoustic engines[END_REF]. This leads to

ϕ (f ) conv = C f ṁf (T H -T ∞ ) = h (f ) conv (T H -T ∞ ) . (24)
The same approach is used to estimate the heat flux h (s) conv convected away from the electrical heat resistance to the stack, which gives:

ϕ (s) conv = φC f ṁs (T H -T ∞ ) = h (s) conv (T H -T ∞ ) . (25)
These estimates of h conv are finally reported in the boundary condition at the stack/waveguide interface, Eq. [START_REF] Gedeon | DC gas flows in Stirling and pulse tube cryocoolers[END_REF]. c. Summary. The complete problem of heat transfer in the presence of acoustic streaming is now described by Eqs. [START_REF] Penelet | Experimental and theoretical study of processes leading to steady-state sound in annular thermoacoustic engines[END_REF] and [START_REF] Gedeon | DC gas flows in Stirling and pulse tube cryocoolers[END_REF]. Therefore, following the approach of section III B, one can transform these equations into ordinary differential equations: introducing once again the space variable z = x -x s , we seek a solution to Eqs. [START_REF] Penelet | Experimental and theoretical study of processes leading to steady-state sound in annular thermoacoustic engines[END_REF] and [START_REF] Gedeon | DC gas flows in Stirling and pulse tube cryocoolers[END_REF]

in the form z ≤ 0, T (i,o) (z, t) ≈ T ∞ + ∆T (t)e z/l (i,o) s (t) , (26a)
z ≥ 0, T (i,o) (z, t) ≈ T ∞ + ∆T (t)e -z/l (i,o) f (t) . ( 26b
)
where l

(i,o)
s,f (t) are the four characteristic lengths associated to both inner and outer zones in the stack and in the waveguide. Finally, using the same approximations than those formulated in section III B, and including the equation describing thermoacoustic amplification, Eq. ( 8), leads to the following set of equations:

l (i) s ∆T d ∆T dt + d l (i) s dt = κ s + Γ κ P 2 l (i) s - l (i) s τ s - Φρ f C f ρ s C s Γ (s) str P 2 , (27a) l (o) s ∆T d ∆T dt + d l (o) s dt = κ s + Γ κ P 2 l (o) s - l (o) s τ s + Φρ f C f ρ s C s Γ (s) str P 2 , (27b)
l (i) f ∆T d ∆T dt + d l (i) f dt = κ f l (i) f - l (i) f τ f -ṁf , (27c)
l (o) f ∆T d ∆T dt + d l (o) f dt = κ f l (o) f - l (o) f τ f + ṁf , (27d)
d ṁf dt + ṁf θ f = Γ (f ) str P 2 θ f , (27e)
1 2 λ s + Γ λ P 2 1 l (i) s + 1 l (o) s + λ f 2 1 l (i) f + 1 l (o) f + C f ṁf + ΦC f Γ (s) str P 2 ∆T = Q πr 2 i , (27f)
d P dt = ǫ ampl P. (27g)
The last point concerns the calculation of the thermacoustic amplification rate ǫ ampl (together with the frequency Ω of self-sustained oscillations) which must be estimated from a monodimensional temperature field. In the following, this amplification rate ǫ ampl (t) = ǫ ampl [T (z, t)] is calculated from the average temperature field T = 1 2 T (i) + T (o) . The above set of equations [START_REF] Peters | Acoustically induced periodic vortex shedding at sharp edged open channel ends: simple vortex models[END_REF], which reduces to Eqs. [START_REF] De Waele | Basic treatment of onset conditions and transient effects in thermoacoustic Stirling engines[END_REF] in the absence of streaming (Γ (s,f ) str = 0), can be transformed into a set of ordinary differential equations and solved numerically.

IV. THEORETICAL RESULTS

Before calculating the transient regime, the theoretical modeling of section III is first used to compute the onset heat power supply Q onset in function of the stack position x s . This is done by setting P = 0 in Eqs. [START_REF] De Waele | Basic treatment of onset conditions and transient effects in thermoacoustic Stirling engines[END_REF], and by adjusting the heat supply Q(t) = Q 0 to Q onset in order that ǫ ampl = 0. The resulting theoretical stability curve Q onset (x s ) is presented in Fig.

(dashed line). The results show reasonable agreement between experiments and theory,

with unavoidable differences due to the approximations of the model.

Once the theoretical stability curve is determined, the protocol used to compute the transient regime is chosen in accordance with the experimental protocol. The initial heat supply Q(t ≤ 0) = Q 0 is chosen slightly below Q onset , while the initial peak pressure amplitude P is fixed arbitrarily to the small value P (t = 0) = 10 -8 P a. Then, at time t = 0, a power increment ∆Q is applied and the governing set of ordinary differential equations, Eqs. [START_REF] Peters | Acoustically induced periodic vortex shedding at sharp edged open channel ends: simple vortex models[END_REF], is solved numerically using a variable time-step, fourth-order Runge-Kutta method.

In Fig. 8(a), we present the results obtained for the stack position x s = 26.5 cm, where the initial heat power supply Q 0 is fixed to 18.6 , W . The same results are presented in Fig. 8(b) for the stack position x s = 36.5 cm, where the initial heat power supply Q 0 is fixed to 13.7 , W . The results depicted in Fig. 8 show that the time of occurrence of the onset of thermoacoustic instability is controlled by the power increment ∆Q, and that the amplitudes of acoustic pressure in steady state are roughly of the same order of magnitude than those observed in experiments. The models also predicts larger steady state acoustic pressures for x s = 36.5 cm than for x s = 26.5 cm, which is consistent with the experimental results of Figs 4 and5. However, when comparing in details the theoretical results with the experimental results of Figs 4 and5, it appears clearly that the model is unable to reproduce the actual dynamics of wave amplitude evolution: there exists very small overshoots of wave amplitude growth which are however much lower than those observed in experiments, and the switch on/off process is not predicted by theory.

In order to get a deeper physical insight about the engine's operation, it is useful to calculate the transient regime when each of the saturating processes is considered independantly. This is realized in Fig. 9(a) where x s is fixed to 36.5 cm and ∆Q/Q 0 = 16% (Q 0 = 13.7 W ).

For an adequate readability, the computed data are only presented from time t = 25 s to t = 100 s. The gradual evolution of the characteristic temperatures T H (t) = T (z = 0, t), T S (t) = T (-l s 0 , t), and T F (t) = T (l f 0 , t) are also presented in Fig. 9, with l s 0 ≈ 1.52 cm and l f 0 ≈ 1.36 cm. The transient regime is computed in different configurations depending on whether both thermoacoustic streaming and acoustic streaming are considered (solid lines)

or not (• and ⋄). The impact of cancelling the time delay θ f of streaming establishment is also examined ( ). From the analysis of the evolution of acoustic pressure P (t), it firstly appears that the major contribution to wave saturation is due to the thermoacoustic heat pumping (dashed lines, •) which, however, cannot be considered as the only contributor to wave saturation (otherwise dashed lines would coincide with solid lines). The only contribution of acoustic streaming (dash-dotted lines, ⋄) leads to a steady state acoustic pressure which is not significantly larger than the one due to the only thermoacoustic pumping: this means that, from our simplified model, acoustic streaming contributes to wave saturation.

Another interesting point concerning the only contribution of acoustic streaming to wave saturation is that there exists a significant overshoot of wave amplitude growth, clearly visible around t ≈ 35 s. Moreover, the results show that if the characteristic time θ f of streaming establishment is discarded (θ f = 0, ), the overshoot disappears. From the analysis of the evolution of the temperatures T H , T S and T F , it appears that the wave amplitude growth is accompanied by a diminution of T H leading to a diminution of the thermoacoustic amplification rate. However, it is difficult to distinguish the different curves from the only analysis of T H while significant differences clearly appear from the gradual evolutions of both T S and T F . This clearly means that thermoacoustic amplification is not only controlled by the temperature difference across the stack but also by the shape of the temperature field in the entire thermoacoustic core.

V. CONCLUSION

We presented an experimental and theoretical study dealing with the onset of selfsustained acoustic waves in a standing wave thermoacoustic engine. The results clearly show that even in the simple thermoacoustic device considered in this study, a rather simple modeling is unable to reproduce the complicated dynamics of wave amplitude evolution observed in experiments. The model indicates that the major contribution to wave saturation is due to thermoacoustic heat pumping by acoustic waves, but also that forced convection due to acoustic streaming should be worth considering. Though the complicated dynamics of wave amplitude evolution are not completely reproduced by the model, it seems reasonable to suggest that the effects observed in experiments could be due to the gradual evolution of the temperature field along the thermoacoustic core. The results also indicate that acoustic streaming impacts the temperature field within a time scale which is significantly larger than the time scale ǫ -1 ampl of thermoacoustic amplification: this could explain the overshoot process and, perhaps, the switch on/off process. It is thus challenging to improve the theoretical description of the device in order to get a better agreement with experiments. In our opinion there are two major drawbacks in the simplified model presented in this paper. The first one concerns our rough estimate of the heat exchange coefficients h (s,f ) conv associated to the heat taken by the fluid near the stack/waveguide interface. The second drawback concerns the rough simplification of the equations describing unsteady heat transfer : by transforming these partial differential equations into ordinary differential equations, one may lose important information linked to the details of the temperature distribution and its evolution with time. It seems interesting for future works to give up the simplification mentioned above, and to solve the unsteady heat transfer equations using a finite difference numerical scheme.

Finally, it would be interesting to include in the analysis the additional processes of wave saturation (nonlinear propagation in the open-ended tube and minor losses at the edges of the stack).

Fig. 1 .

 1 photograph of this thermoacoustic oscillator is presented in Fig.1. It consists of a glass tube (length L=49 cm, inner radius r i =26 mm, outer radius r e =30 mm) open to free space

||

  ṁ(x s + l w 0 , r)| r dr = Γ ṁ(x s -l s 0 , r)| r dr = Γ (s)

FIG. 1 .FIG. 2 .FIG. 3 .FIG. 4 .FIG. 5 .

 12345 FIG. 1. (a) Photograph of the experimental apparatus. (b) photograph of the hot side of the stack.

FIG. 7 .

 7 FIG. 7. Schematic representation of the simplistic approach used to account for heat convection by acoustic streaming. (a): schematic representation of the streamlines associated to acoustic streaming ; (b): schematic representation of the streaming velocity field inside both stack and waveguide ; (c) and (d): the heat convected from the stack/waveguide interface is estimated by unwrapping the vortex cells and by considering that the interface is equivalent to an isothermal grid at constant temperature T H crossed by a steady flow.

ACKNOWLEDGMENTS

This research is funded from the French National Agency for Research under grant agreement ANR 2010 BLAN 0302 01. The authors are also indebted to Denis Bruneau from the I2M laboratory (UMR CNRS 5295) for helpful discussions about the modeling of heat transfer.