J Peralez
email: johan.peralez@ifpen.fr

P Tona

M Nadri
email: nadri@lagep.univ-lyon1.fr

P Dufour
email: dufour@lagep.univ-lyon1.fr

A Sciarretta
email: antonio.sciarretta@ifpen.fr

Optimal Control for an Organic Rankine Cycle on board a Diesel-Electric Railcar

Keywords: Dynamic programming, backward reachability, real-time optimization, waste heat recovery, Rankine cycle, vehicle energy management

This paper addresses the problem of maximizing the power produced by an organic rankine cycle (ORC) waste heat recovery system on board a diesel-electric railcar. A simplified model of the system allows the formulation of an optimal control problem that can be solved via dynamic programming (DP). To increase the smoothness and the accuracy of the solution obtained offline using the implementation of DP known as level-set DP, an improved version is developed, making use of adaptive grids for discretization. The analysis of the resulting optimal trajectory for the ORC control problem provides useful insight for both control design and system design. Based on these results, the optimal control problem is reformulated allowing online implementation via dynamic real-time optimization. The proposed approach is validated on a realistic simulator, showing significant benefits in the amount of recovered energy when compared with the classical, quasi-static approach found in ORC literature.

Introduction

With new technologies continuously being introduced to further reduce vehicle fuel consumption, new control problems arise that must be efficiently solved to make the most of the innovation potential. This is for instance the case of engine waste heat recovery (WHR) systems based on the Rankine thermodynamic cycle. Applications of this technology exist, at different levels of maturity, for cars, trucks, trains and ships (see Sprouse and Depcik [START_REF] Sprouse | Review of organic rankine cycles for internal combustion engine exhaust waste heat recovery[END_REF] for a recent overview). For some applications, such as WHR via organic rankine cycle (ORC) from Diesel engines of long-haul trucks, very interesting fuel consumption reductions have been reported (up to 6%, according to Stanton [START_REF] Stanton | Systematic development of highly efficient and clean engines to meet future commercial vehicle greenhouse gas regulations[END_REF]). As shown for instance in Peralez et al. [START_REF] Peralez | Towards model-based control of a steam rankine process for engine waste heat recovery[END_REF], Horst et al. [START_REF] Horst | Dynamic heat exchanger model for performance prediction and control system design of automotive waste heat recovery systems[END_REF] and Xie and Yang [START_REF] Xie | Dynamic behavior of rankine cycle system for waste heat recovery of heavy duty diesel engines under driving cycle[END_REF], control of Rankine-based WHR systems for vehicular applications is far from being trivial.

In these applications, contrary to conventional (stationary) Rankine cycles, hot source dynamics depends on driving conditions and acts as a (highly-transient) time-varying disturbance. In this context, an effective control system is essential to attain satisfactory performance over a broad range of operating conditions, where "satisfactory" means safe, first of all, but also energetically efficient. In the still scarce literature on control design for Rankine WHR (see for instance the comprehensive overview presented in [START_REF] Sprouse | Review of organic rankine cycles for internal combustion engine exhaust waste heat recovery[END_REF]), experimentally-validated strategies proposed so far for transport-related applications (such as Endo et al. [START_REF] Endo | Study on maximizing exergy in automotive engines[END_REF], Peralez et al. [START_REF] Peralez | Improving the control performance of an organic rankine cycle system for waste heat recovery from a heavy-duty diesel engine using a model-based approach[END_REF]) focus on low-level control to ensure safe operation.

Few papers address supervisory control of Rankine WHR in the framework of integrated powertrain energy and emission management. The solution based on the Pontryagin minimum principle (PMP) presented in Merz et al. [START_REF] Merz | On the optimal thermal management of hybrid-electric vehicles with heat recovery systems[END_REF] includes a very simple dynamic Rankine model but cannot be implemented on line as such, while the heuristic strategy, also inspired by PMP, proposed in Willems, Frank et al. [START_REF] Willems | Integrated energy and emission management for diesel engines with waste heat recovery using dynamic models[END_REF] does not use any explicit WHR model in the controller. A few other studies in ORC literature address the optimization of an economic objective, mostly using quasisteady-state models. Among them, Quoilin et al. [START_REF] Quoilin | Dynamic modeling and optimal control strategy of waste heat recovery Organic Rankine Cycles[END_REF][START_REF] Quoilin | Thermo-economic optimization of waste heat recovery organic rankine cycles[END_REF] focused on the optimal evaporating temperature in the case of an ORC with a volumetric expander, while Manente et al. [START_REF] Manente | An organic rankine cycle off-design model for the search of the optimal control strategy[END_REF]) considered an ORC with air-cooled condenser. In the field of ORC stabilization, Hou et al. [START_REF] Hou | Supervisory predictive control of evaporator in organic rankine cycle (ORC) system for waste heat recovery[END_REF], Zhang et al. [START_REF] Zhang | Constrained predictive control based on state space model of organic rankine cycle system for waste heat recovery[END_REF] used a dynamic (linear) model to design a model predictive controller: a cost function associated to stabilization performance was then optimized.

On the other hand, several problems of transient control of energy production (see e.g Guzzella and Sciarretta [START_REF] Guzzella | Vehicle propulsion systems[END_REF]) on board vehicles in the presence of constraints have been addressed and solved via optimal control methods and in particular via dynamic programming (DP) (Bellman [3]).

However, to the authors' knowledge, no solution for online optimization of power production taking into account Rankine dynamics has been reported so far.

This paper addresses the problem of maximizing the power produced by an ORC WHR system on board a diesel-electric railcar. 1 An experimentally-validated, reduced-order dynamic model allows the formulation of an optimal control problem that can be solved via DP. To increase the smoothness and the accuracy of the solution obtained using the state-of-the-art implementation of DP known as Level-Set DP (Elbert et al. [START_REF] Elbert | Implementation of dynamic programming for n -dimensional optimal control problems with final state constraints[END_REF]), a novel version of that algorithm, called Adaptive-Grid DP, is developed. The analysis of the resulting optimal trajectory for the ORC control problem provides useful insight for the design of an online control strategy, based on dynamic real-time optimization (DRTO, see e.g Würth et al. [START_REF] Würth | Neighboring-extremal updates for nonlinear model-predictive control and dynamic real-time optimization[END_REF]).

This paper is organized as follows. Section 2 introduces the application and its specificities. A simplified model of the system and the formulation of the optimal control problem are proposed. Section 3 details the novel Adaptive-Grid DP algorithm -illustrated with a simple energy management example in Appendix A -providing an offline solution to the optimal control problem. Based on these results, the optimal control problem is reformulated in Section 4, allowing online implementation via DRTO. This approach is validated on a realistic simulator, showing significant increase in the amount of the energy recovered during a typical mission, when compared with the classical approach found in ORC literature, in which a quasi-static model is used for optimization.

Modeling and problem formulation

Problem statement

The system under consideration is represented in Figure 1. It is a prototype ORC system for WHR from a diesel engine generating set ("Power Pack") manufactured by Alstom Transport for installation in diesel-electric railcars. A pump pressurizes an organic fluid (R365mfc) in liquid state and circulates it through a closed circuit at flow rate * m. Via a by-pass valve, a fraction V o,e of engine exhaust gas is fed to a heat exchanger (evaporator) which transfers gas heat to the working fluid. The vaporized fluid then enters a turbine, which converts some of the kinetic energy into work, driving a generator connected to the Power Pack electrical network. Lastly, the fluid is recondensed before returning to the pump at low pressure. Cooling is provided by a fan with controllable air flow * m A . Notice that, although cooling is rarely taken into account in ORC optimization problems, it represents a significant cost in energy balances for ORC applications on board vehicles (Horst et al. [START_REF] Horst | Prediction of dynamic rankine cycle waste heat recovery performance and fuel saving potential in passenger car applications considering interactions with vehicles' energy management[END_REF]). As the ratio between the power consumed by the auxiliaries and the power used for traction is relatively large on board trains (which benefit from a low aerodynamic friction and from a low rolling resistance coefficient), it can be safely assumed that all the power produced by the WHR system is fed to the electrical network. Consequently, the objective here is to find the control strategy that maximizes the energy J recovered by the WHR system along a given mission of duration t f :

J = t f 0 P net (t) dt, (1)
where P net is given by the power produced by the turbine, P turb , minus the powers consumed by the pump and the cooling system, P pump and P A ,

P net (t) = P turb (t) -P pump (t) -P A (t) (2)
The mission profile considered herein, a portion of a typical intercity train trip, results in a set of external conditions for the ORC system, namely for the engine exhaust gas temperature T exh (t) and mass-flow * m exh (t). A large variation of ambient temperature T A (t) is also considered as shown in Figure 2. To solve this problem, a model is needed to predict the evolution of the variables in the expressions of powers [START_REF] Bao | A review of working fluid and expander selections for organic rankine cycle[END_REF]. Model dimension should be small so as not to be confronted with the well-known "curse of dimensionality" (Powell [START_REF] Powell | Approximate Dynamic Programming: Solving the curses of dimensionality[END_REF]). The so-called moving boundaries method (see e.g. McKinley and Alleyne [START_REF] Mckinley | An advanced nonlinear switched heat exchanger model for vapor compression cycles using the moving-boundary method[END_REF]) allows a realistic dynamic representation of heat exchangers with a limited number of state variables. Cheng et al. [START_REF] Cheng | Nonlinear observer design for two-phase flow heat ex-changers of air conditioning systems[END_REF] proposes a model that considers an average wall temperature, while Peralez et al. [START_REF] Peralez | Improving the control performance of an organic rankine cycle system for waste heat recovery from a heavy-duty diesel engine using a model-based approach[END_REF] shows experimentally that, for an ORC, the wall temperatures have the slowest dynamics. Thus, it is possible to neglect the (much faster) fluid dynamics, thereby reducing the number of state variables.

In the following subsection, the thermodynamic cycle is first explained and described by mathematical expressions, then these dependencies are reduced by assumptions and simplifications, resulting in an implicit dynamic model.

Two-state ORC model

The thermodynamic cycle of the working fluid is represented in Figure 3. In the following, subscripts 1 to 4 will refer to the working fluid at each component outlet, namely the evaporator, the turbine, the condenser and the pump.

Recoverable mechanical power from the turbine is usually expressed (Bao and Zhao [START_REF] Bao | A review of working fluid and expander selections for organic rankine cycle[END_REF]) in terms of a constant mechanical efficiency η turb , of the working fluid mass-flow * m and of the enthalpy gradient h 2 -h 1 by

P turb = * m(h 2 -h 1 ) η turb . (3)
Similarly, the efficiency η pump can be considered constant and the power consumed by the pump can be written (Borsukiewicz-Gozdur [START_REF] Borsukiewicz-Gozdur | Pumping work in the organic rankine cycle[END_REF]) in terms of the pressure gradient p 1 -p 3 and the fluid density ρ 3 as following

P pump = * m p 1 -p 3 ρ 3 η pump . (4)
The power consumed by the condenser fan can be considered proportional to the delivered air flow * m A (Manente et al. [START_REF] Manente | An organic rankine cycle off-design model for the search of the optimal control strategy[END_REF]) and is given by

P A = * m A k A , (5)
where k A is the specific power consumption of the condenser fan.

Considering the thermodynamic cycle shown in Figure 3, fluid superheating at evaporator outlet (point 1) is assumed to be perfectly regulated at a constant value by the pump mass-flow. The fluid at the pump inlet (point 3) is assumed to be in the saturated liquid state. The slow dynamics of the system, i.e. the thermal dynamics of exchangers walls, are taken into account while the fluid is considered to be in equilibrium. [START_REF] Dufour | On nonlinear distributed parameter model predictive control strategy: on-line calculation time reduction and application to an experimental drying process[END_REF] where M w is an equivalent thermal capacity, S f the wall-fluid exchange area, S A the wall-air exchange area and

η A = 1 -exp - α A S A * m A c p,A
reflects the decrease in heat exchange efficiency for large air mass flows (McKinley and Alleyne [START_REF] Mckinley | An advanced nonlinear switched heat exchanger model for vapor compression cycles using the moving-boundary method[END_REF]), whereas η exh can be considered constant because of the lower mass flow involved.

In the following, Tf,e and Tf,c will be assumed equal to the evaporation temperatures respectively of the high-pressure and of the low-pressure sides of the ORC. Tf,e (resp. Tf,c ) is then an unimodal function of p 1 (resp. of p 3 ). Likewise, h 1 and ρ 1 depend only on p 1 , while h 3 and ρ 3 depend on p 3 .

Fluid energy equilibrium in the exchangers yields:

       0 = h 1 -S f,e ᾱf,e Tw,e -Tf,e * m -h 3 0 = h 3 -S f,c ᾱf,c Tw,c -Tf,c * m -h 2 , (7)
where the mass flow * m is considered homogeneous throughout the circuit. Assuming supersonic flow condition, the pressure at turbine outlet is neglected (Feru et al. [START_REF] Feru | Experimental validation of a dynamic waste heat recovery system model for control purposes[END_REF]) and the mass flow can be expressed as follows:

* m = c d S 2 ρ 1 p 1 , (8)
where c d is a (constant) discharge coefficient, and S the equivalent section of the turbine nozzles. Assuming a constant isentropic efficiency η is of the turbine

h 2 = h 1 -η is (h 1 -h is,2 ), ( 9
)
where h is,2 corresponds to the enthalpy whose entropy is equal to that of point 1 and pressure is equal to that of point 3. h is,2 then depends on p 1 and p 3 .

Using the previous expressions for * m and h 2 , system (7) can be written in the following implicit form:

0 = ϕ( Tw,e , Tw,c , p 1 , p 3 ). ( 10
)
At each instant, from the knowledge of state variables Tw,e and Tw,c , the implicit equation ( 10) can be solved numerically by iterating over p 1 and p 3 , as illustrated in Figure 4. Through the computation of fluid properties the expression of wall thermal dynamics can be found [START_REF] Dufour | On nonlinear distributed parameter model predictive control strategy: on-line calculation time reduction and application to an experimental drying process[END_REF], as well as other key variables (such as * m in Figure 4). The next subsection discusses model validation for transient operating conditions.

Experimental validation of the reduced model

In order to validate the single-state model of the heat-exchanger, the evaporator model response to transient external conditions is compared to experimental data. Transient conditions of exhaust gas entering the evaporator are shown in Figure 5(a).

The numerical value for the average heat exchange coefficient ᾱf,e has been fitted against the pump mass flow * m experimental data. Figure 5(b) compares model predictions of key working fluid variables at evaporator outlet, namely the mass flow * m, the pressure p 1 and the temperature T 1 , with the corresponding measurements. Relative average deviations are found to be about 3.3% for working fluid mass flow, 2% for pressure and less than 1% for temperature. This accuracy, demonstrated in transient conditions, motivates the use of the simplified model for a dynamic optimization purpose.

Optimal control formulation for ORC-based WHR on board diesel-electric railcars

Considering the mathematical model of an ORC described above, the formulation of the problem stated in Section 2.1 for WHR on board diesel-electric railcars will be developed in the following.

The aim is to find a control strategy that maximizes the energy recovered J, defined by equations ( 2)-( 5), along a given mission of duration t f . The cost functional J to be maximized depends on time, control inputs u = {V o,e , * m A } and state vector x = { Tw,e , Tw,c } whose dynamics are defined by the two-state ORC model equations ( 6)- [START_REF] Guzzella | Vehicle propulsion systems[END_REF]. This maximization problem must be solved under the following safety constraint on pressure p 1 whereas bounds on the control inputs are set by technical considerations: the fraction of exhaust gas (u 1 ) is bounded by [0, 1], whereas the air mass flow provided by the fan (u 2 ) is between 0 and 4 kg s -1 .

Then, the optimal control problem of the ORC can be formulated under the following minimization form:

min u(t) t f 0 -P turb (x) + P pump (x) + P A (u) dt, (11)
m exh c p,exh (T exh -x 1 )/M w,e S f,c ᾱf,c ( Tf,c -x 2 )/M w,c + 1 -exp - α A S A u 2 c p,A u 2 c p,A (T A -Tw,c )/M w,c   , (12)
x(0) = x 0 , (13)
       0 = h 1 -S f,e ᾱf,e
Tw,e -Tf,e * m -h 3

0 = h 3 -S f,c ᾱf,c Tw,c -Tf,c * m -h 2 , (14)
p 1 (x(t)) ≤ 25 bar, ∀t ∈ [0, t f ], (15)
u(t) ∈ [0, 1] × [0, 4], ∀t ∈ [0, t f ]. ( 16
)
Remark 1. Pump mass flow * m is not included as a decision variable in the above formulation as it is used to maintain superheating at a constant value. Indeed, it is usually accepted that, in ORCs, superheating should be maintained as low as possible in order to maximize cycle efficiency (Sprouse and Depcik [START_REF] Sprouse | Review of organic rankine cycles for internal combustion engine exhaust waste heat recovery[END_REF], Bao and Zhao [START_REF] Bao | A review of working fluid and expander selections for organic rankine cycle[END_REF]). Based on two-time-scale considerations (working fluid dynamics is much faster than wall thermal dynamics), superheating is assumed perfectly regulated by pump mass flow (Peralez et al. [START_REF] Peralez | Improving the control performance of an organic rankine cycle system for waste heat recovery from a heavy-duty diesel engine using a model-based approach[END_REF] showed experimentally that such assumption is realistic, even in transient driving conditions).

Offline dynamic optimization

The main advantage of DP above other optimization methods is that it gives sufficient conditions for global optimality. Its main drawback is that the computational load grows exponentially with the number of state variables and control inputs (Powell [START_REF] Powell | Approximate Dynamic Programming: Solving the curses of dimensionality[END_REF], Lee and Wong [START_REF] Lee | Approximate dynamic programming approach for process control[END_REF]). DP is usually applied writing the system dynamics in a discrete form and the accuracy of the solution increases with the number of points considered. With the basic implementation, the necessary computational effort rapidly becomes untenable, especially in the presence of active state constraints as in the ORC problem.

In this section, the dynamic programming principle and the state-of-the-art implementation of DP known as Level-Set DP are first recalled. Then, in order to increase the smoothness and the accuracy of the solution while keeping an acceptable computational load, a novel version of that algorithm, called Adaptive-Grid DP, is developed. The solution obtained for the ORC optimal control problem is then analyzed, in order to get insight for the online strategy developed in Section 4.

Level-set Dynamic Programming principle 3.1.1. Dynamic Programming principle

The generic optimal control problem to solve for a dynamic system with n state variables and m control inputs can be written as min

u(t) J (u (t)) (17) s.t. ẋ(t) = f (x(t), u(t), t) (18) x(0) = x 0 (19) x(t) = X (t) ⊆ R n (20) u(t) = U(t) ⊆ R m (21)
where

J (u (t)) = g (x(t f )) + t f 0 h (x(t), u(t), t) dt, (22)
is the cost functional. Since dynamic programming is a numerical algorithm, discretization of time, state space, and control space is required. Let the discrete-time model be

x k+1 = f k (x k , u k ), k = 0, 1, ..., N -1 (23)
where x k ∈ X k is the state vector and u k ∈ U k the control signal vector.

Based on the principle of optimality introduced by Bellman [3], DP proceeds backwards in time to evaluate the optimal cost-to-go function J k (x i ) at every node x i of the discretized time-state space Sundström et al. [START_REF] Sundström | On implementation of dynamic programming for optimal control problems with final state constraints[END_REF]:

Step 1: Initialization

J N (x i ) = g N (x i ), for x i ∈ X N ∞ else. ( 24
)
Following steps: Backward iteration for k = N -1 to 0

J k (x i ) = min u k ∈U k {g k (x i , u k ) + J k+1 (f k (x i , u k ))}. ( 25
)
The optimal control is the argument that minimizes the right-hand side of equation ( 25) for each x i at time index k, which yields the optimal control policy π = {u 0 (x), u 1 (x), ..., u N -1 (x)}. This map is used to find the optimal control signal during a forward simulation of model [START_REF] Petit | Optimal drive of electric vehicles using an inversion-based trajectory generation approach[END_REF] starting from a given initial state x 0 . Since the resulting state trajectory does not generally coincide with the nodes x i , interpolation is used.

Level-set algorithm

A major issue to consider when implementing a DP algorithm, is the definition of the costto-go function for an infeasible state vector. Grid points that are not backward-reachable should have infinite cost as in [START_REF] Powell | Approximate Dynamic Programming: Solving the curses of dimensionality[END_REF]. However, this causes numerical problems as illustrated in Figure 6. When running step k of backward iteration, the cost-to-go J k+1 (x) is known for all grid points x i . In order to evaluate the cost-to-go J k (x p ), all possible control candidates u ∈ U k are applied. In Figure 6, the point f k (x p , u 1 ) is in the backward-reachable space. But since the value of J k+1 (x i+1 ) is infinite, interpolation will lead to consider that the cost-to-go of x i+2 k is infinite too. The classical method to deal with this problem consists in using a large but finite value for an infeasible state vector. The value of this penalty represents a critical parameter that is difficult to calibrate. Nevertheless, the penalty method results in a steep gradient of the cost-to-go function near the boundary and cannot completely solve the aforementioned problem.

Another approach is based on the concept of backward-reachable space, as proposed in [START_REF] Back | Predictive powertrain control for hybrid electric vehicles[END_REF][START_REF] Sundström | On implementation of dynamic programming for optimal control problems with final state constraints[END_REF] for single-input single-state systems. For larger systems, the level-set DP algorithm, introduced by Elbert et al. [START_REF] Elbert | Implementation of dynamic programming for n -dimensional optimal control problems with final state constraints[END_REF], evaluates the backward-reachable space in parallel with performance index evaluation. A level-set function I is considered that associates a real value to the state x: where I(x) ≤ 0 if x is in reachable space and I(x) > 0 else. Such a function can represent nonconvex regions, moreover, the number of state variables and control inputs is not limited. The level-set function is evaluated for grid points x i at each time index k yielding a grid evaluation I k x . At the next backward iteration step, interpolation is used to determine whether a point x is backward reachable or not. It is assumed that the final state constraint is given as a target set T , which is defined by a level-set function h(x) h :

I : X ⊆ R n → R, (26)
X N → R, where X N ⊆ R n (27) T = {x ∈ X N |h(x) ≤ 0}. ( 28
)
Then the DP principle is applied, where the cost-to-go function to minimize is the level-set function I:

Step 1: Initialization

I N (x i ) = h(x i ). ( 29
)
Following steps: Backward iteration for k = N -1 to 0

I k (x i ) = min u k ∈U k {I k+1 (f k (x i , u k ))}. ( 30
)

A new adaptive-grid DP algorithm

At the expense of an acceptable increase of computational load, the level-set DP algorithm introduced in Elbert et al. [START_REF] Elbert | Implementation of dynamic programming for n -dimensional optimal control problems with final state constraints[END_REF] can be modified to provide a smoother and more accurate solution. Differences with original algorithm consist in the second discretization step of state space in the backward algorithm, defined by ( 33)-(37), and in the finer discretization of control space in the forward algorithm, defined by (44).

Adaptive-grid DP backward algorithm

To improve the accuracy of the solution, the backward algorithm introduced in Elbert et al. [START_REF] Elbert | Implementation of dynamic programming for n -dimensional optimal control problems with final state constraints[END_REF] is modified as follows.

Step 1: Initialize k=N and the level-set and cost-to-go functions as:

I N (x i ) = h(x i ) (31) J N (x i ) = g N (x i ). ( 32
)
Considering that x i ∈ R n let x i = x i (1) ... x i (n) . A second discretization of state space XN is introduced:

XN = {x 1 N , x2 N , ..., xq N }, (33)
where each j th component of vector xi N is bounded below by x N (j) and above by xN (j) as follows:

x N (j) = min i∈{1,...,q}

{x i N (j) | 0 < I N (x i N )} (34) xN (j) = max i∈{1,...,q} {x i N (j) | 0 < I N (x i N )} (35)
As illustrated in Figure 7 this second discretization may improve the accuracy of level-set and cost-to-go estimations:

ĨN (x i ) = h(x i ) (36) JN (x i ) = g N (x i ). ( 37
)
Step 2: Reduce k by 1 and update the level-set function by

I k (x i ) = min u k ∈U k { Ĩk+1 (f k (x i , u k ))}. (38)
A second discretization Xk is computed in the same way as for the first step. Then the level-set function is updated by

Ĩk (x i ) = min u k ∈U k { Ĩk+1 (f k (x i , u k ))}. ( 39
)
Step 3: For each grid point xi , find the set of control signals for which the system trajectory ends up inside the backward-reachable space at the next time step

U F k (x i ) = {u k ∈ U k | Ĩk+1 (f k (x i , u k )) ≤ 0}, (40)
and the one control candidate that minimizes the level-set function * u k (x i ) = arg min

u k ∈U k { Ĩk+1 (f k (x i , u k ))}. ( 41
)
Step 4: Update the optimal cost-to-go by the following rule: if at least one valid control candidate is found, i.e., U F k (x i ) = ∅, then calculate the cost-to-go based upon the optimal candidate Jk (x i ) = min

u k ∈U F k (x i ) {g k (x i , u k ) + Jk+1 (f k (x i , u k ))}. ( 42
)
If, however, the grid point is not backward-reachable, then calculate the cost-to-go based on the control input *

u k (x i ) Jk (x i ) = g k (x i , * u k ) + Jk+1 (f k (x i , * u k )) ( 43
)
and repeat steps 2 -4 until k = 0.

Adaptive-grid DP forward algorithm

To increase smoothness, the forward algorithm introduced in Elbert et al. [START_REF] Elbert | Implementation of dynamic programming for n -dimensional optimal control problems with final state constraints[END_REF] is used with a finer discretization of the control space:

Ũk = {u 1 k , ..., u r k U k , ..., u r k } (44)
Step 1: Initialize k = 0 and x 0 at initial condition.

Step 2: Increase k by 1 and find the feasible control candidates

Ũ F k (x k ) = {u k ∈ Ũk | Ĩk+1 (f k (x k , u k )) ≤ 0}, (45)
Step 3: Find the optimal control input

u 0 k (x k ) = arg min u k ∈ Ũ F k {g k (x k , u k ) + Jk+1 (f k (x k , u k ))} (46)
Step 4: Simulate the system using the optimal control input

x k+1 = f k (x k , u o k ) (47)
and repeat steps 2 -4 until k = N .

Remark 2. Note that the computational effort involved in the finer discretization Ũk (44) is limited. Since steps 2 -4 of the algorithm are applied to one value of x k , the number of additional calls to the model is N (r -r).

Before application to the ORC case study, the new algorithm has been validated on a simple energy management example, where an analytic solution exists. In this case, as shown in Appendix A, the results obtained with the adaptive grid DP are close to the analytic solution even with a coarse discretization.

DP solution for ORC-based WHR on board diesel-electric railcars

The ORC optimal control problem exposed in Section 2.4 is solved using the new adaptive-grid DP algorithm with numbers of grid points set to N = N x1 = N x2 = N u1 = N u2 = 21 and with a sampling time set to 1 second. Figure 8 shows the system trajectory obtained. High and low pressures (respectively in the evaporator and in the condenser) periodically rise in response to the increase in the exhaust gas mass flow from the mission profile (Figure 2). In response, evaporator by-pass is open up to limit the high pressure at the safety constraint [START_REF] Lee | Approximate dynamic programming approach for process control[END_REF]. On the other hand, the second control input exhibits a noncausal behavior: for instance, it can be observed that a few seconds before the external condition * m exh rises, * m A decreases at first. Such a behavior can be explained considering that the efficiency of heat exchange with ambient air improves when condenser wall temperature (and consequently pressure p 3 ) is higher.

In order to evaluate the benefits of the new algorithm, the ORC optimal control problem was also solved with the original level-set DP algorithm. It was found that the number of grid points N = N x1 = N x2 = N u1 = N u2 = 34 was required to reach an equivalent accuracy. It resulted in a computational time about 4.5 times smaller for the adaptive grid DP. Despite the reduction in computational load obtained with the adaptive-grid algorithm, the execution times of the DP are still too long for online implementation. However, the DP solution remains very useful as an analysis tool to understand (optimal) ORC behavior. From the results above, one can infer that tracking control of condenser pressure can be obtained using the air mass flow * m A to regulate condenser cooling, which would in turn ensure the pressure ratio required to run the turbine efficiently. On the other hand, the evaporator by-pass can be used to enforce safety limits in terms of pressure and temperature.

The DP solution also highlights the influence of ambient temperature on the overall efficiency of the ORC system. In fact, the power recovered by the turbine is substantially reduced for high ambient temperatures (during the end of the simulation), whereas the power consumed by the pump is scarcely affected. These considerations will be used as design guidelines for the online supervisory control which is the object of section Section 4.

Other uses of the DP approach

Solving the optimal control problem introduced above can also prove useful at system design stage, especially considering that transient behavior is accounted for in system performance evaluation. The problem can be formulated for variants of the system under consideration, obtained with different sets of model parameters. Thanks to the good accuracy achievable by the adaptivegrid DP algorithm at low levels of discretization, it becomes feasible to study the impact of these (numerous) parameters within reasonable simulation times.

To illustrate this, the influence of a key design parameter -namely the equivalent turbine nozzle surface S -on the ORC efficiency is investigated. The optimal control problem is then solved for different values of S and the corresponding efficiencies are compared in term of mean net recovered power Pnet (Figure 9). Such simulations help us to conclude that an optimal value for the surface S may be found. Working fluid choice or exchangers sizing -taken into account in the two-state model described in Section 2.2 -are examples of other design parameters whose influence could be studied with this approach.

Since the proposed improvement does not rely on a specific structure, application of the new algorithm should also be beneficial to other optimal control problems. Yet, it has not been proved that the novel "adaptive-grid" algorithm consistently outperforms the original level-set algorithm in different contexts than those presented in this article.

Online supervisory control

One way to design an online supervisory control for the ORC is to directly solve the problem exposed in Section 2.4. Unfortunately, the limits of computational resources, particularly stringent in an embedded context, prevent achieving accurate real time solutions. Computational issues are mainly due to: i) the non-linearities of the model, ii) the implicit structure of the model, iii) the multi-input multi-output nature of the optimal control problem. Also, all disturbances need to be known over the whole cycle.

The aim of the following section is to investigate how to reformulate the problem in order to overcome these issues and find a real-time implementable (sub) optimal control strategy, regardless of the cycle considered. This is done in three steps. As a first step, through a new choice of dynamic state variables, it is shown that the ORC model can be written in an explicit form. This results in saving computation time. Moreover, the new state variables (the pressures p 1 and p 3 ) are directly measurable, avoiding observability issues. Then, insight given by the DP solution (Section 3.3) is used to simplify the optimal control problem formulation. In what follows, it is inferred that the evaporator by-pass can enforce the pressure (safety) limit and that the air mass flow should be used to optimize cycle efficiency. Finally, the main result is stated giving a new dynamic real-time optimization strategy.

Explicit two-state ORC model

To transform the implicit ORC model given by the system ( 6)-( 10) into an ordinary differential equation (ODE) system of the following form

       Ṫw =f (T w , p, t) ṗ = - ∂ϕ ∂p (Tw,p) -1 ∂ϕ ∂T w (Tw,p) f (T w , p, t), (48)
where T w = Tw,e Tw,c

T , p = p 1 p 3 T , f = S f,e ᾱf,e ( Tf,e -Tw,e )/M w,e + V o,e η exh * m exh c p,exh (T exh -Tw,e )/M w,e S f,c ᾱf,c ( Tf,c -Tw,c )/M w,c + η A * m A c p,A (T A -Tw,c )/M w,c , , (49)
∂ϕ ∂p =      -(h 1 -h 3 ) ∂ * m ∂p 1 - * m ∂h 1 ∂p 1 -S f,e ᾱf,e ∂ Tf,e ∂p 1 * m ∂h 3 ∂p 3 * m ∂h 2 ∂p 1 - ∂ * m ∂p 1 (h 3 -h 2 ) * m ∂h 2 ∂p 3 - ∂h 3 ∂p 3 -S f,c ᾱf,c ∂ Tf,c ∂p 3      , (50)
and

∂ϕ ∂T w = S f,e ᾱf,e 0 0 S f,c ᾱf,c , (51)
one has to check that ∂ϕ ∂p is of full rank (Hammouri and Nadri [START_REF] Hammouri | An observer design for a class of implicit systems[END_REF]).

Proposition 1. Let Ω be the set which represents the physical domain:

Ω = {(T w , p) ∈ R 2 × R 2 , p 1 ≤ p 3 , p 1 < p crit , p 3 > 0},
where p crit refers to the critical pressure for the considered working fluid. For every (T w , p) ∈ Ω, ∂ϕ ∂p is of full rank.

Proof 1. Consider the following notation

∂ϕ ∂p T a b c d , where                        a = -(h 1 -h 3 ) ∂ * m ∂p 1 - * m ∂h 1 ∂p 1 -S f,e ᾱf,e ∂ Tf,e ∂p 1 b = * m ∂h 2 ∂p 1 - ∂ * m ∂p 1 (h 3 -h 2 ) c = * m ∂h 3 ∂p 3 d = * m ∂h 2 ∂p 3 - ∂h 3 ∂p 3 -S f,c ᾱf,c ∂ Tf,c ∂p 3 .
To show that ∂ϕ ∂p is well defined, it is sufficient to show that its transpose is strictly diagonally dominant, namely that |a| > |b| and |d| > |c|.

Considering the thermodynamic cycle of Figure 3 and assuming that (T w , p) ∈ Ω, it can be deduced that

|a| = (h 1 -h 3 ) ∂ * m ∂p 1 + * m ∂h 1 ∂p 1 + S f,e ᾱf,e ∂ Tf,e ∂p 1 . (52)
Indeed, h 1 and Tf,e (the saturation temperature) are increasing functions of p 1 . Moreover, h 1 > h 3 and the product S f,e ᾱf,e is positive. Finally, since the mass flow can be expressed as

* m = c d S 2 ρ 1 p 1 ,
it is deduced that * m is an increasing function of p 1 .

Using the same steps as above,

|b| = * m ∂h 2 ∂p 1 - ∂ * m ∂p 1 (h 3 -h 2 ) . (53)
From (53) and (52), it yields

|a| -|b| = (h 1 -h 2 ) ∂ * m ∂p 1 + * m ∂h 1 ∂p 1 - ∂h 2 ∂p 1 + S f,e ᾱf,e ∂ Tf,e ∂p 1 .
From Figure 3, the resultant work (h 1 -h 2 ) > 0 which may be recovered by the turbine is an increasing function of the high pressure p 1 , i.e

∂h 1 ∂p 1 - ∂h 2 ∂p 1 > 0. Then |a| > |b|.
Using the same approach, it comes

|c| -|d| = - * m ∂h 2 ∂p 3 + S f,c ᾱf,c ∂ Tf,c ∂p 3 .
Using the expression of * m, it comes that the following condition on system sizing is required

S f,c c d S < * m ∂h 2 ∂p 3 ᾱf,c ∂ Tf,c ∂p 3 -1 . ( 54
)
In practice, this condition is realistic because the right-hand part of (54) evaluated in the physical domain is very large.

Consequently, the condition |d| > |c| is verified. This achieves the verification of Assumption1. Now, considering system (7), T w can be deduced from p 1 and p 3

T w (p 1 , p 3 ) =     * m h 1 -h 3 S f,e ᾱf,e + Tf,e * m h 3 -h 2 S f,c ᾱf,c + Tf,c     . (55)
Then, T w can be removed from the state vector and the ORC dynamics can be described by the following ODE equation:

ṗ = - ∂ϕ ∂p (Tw,p) -1 ∂ϕ ∂T w (Tw,p) f (T w , p, t). (56)

Evaporator pressure control

The original optimal problem with the two decision variables {V o,e , * m A } will be split into two simpler sub-problems. First, a closed-loop controller for pressure p 1 acting on V o,e is designed. Then a dynamic real-time optimization strategy with * m A as a decision variable is described. In order to take into account safety limits, a closed-loop is built to regulate the high-pressure p 1 . For this, the specific structure of dynamic model (56) will be exploited, exhibiting that it is affine in V o,e , the evaporator bypass control input.

u d 1 = -A 11 f 1 (p, Tw,c ) -A 12 f 2 ( Tw,c , u 2 , t) + ṗ1 d A 11 g 2 ( Tw,e , t) (59)
From Section 4.1 it can be easily shown that A 11 > 0 on the physical domain. Moreover g 2 (.) does not vanish for * m exh > 0 and T exh > Tw,e , that is when the diesel engine is on. In consequence u d 1 is well defined. Assuming from DP results analysis, that p 1 should be as high as possible (under the safety limit p max ), the desired high-pressure dynamics can be stated as:

ṗ1 d = k c (p max 1 -p d 1 ) with k c > 0. ( 60
)
Taking into account the physical constraint on by-pass opening (i.e u 1 ∈ [0, 1]), then it comes the following closed-loop model:

           ṗ1 d = k c (p max 1 -p d 1 ) with k c > 0 ṗ3 = A 21 f 1 (p, Tw,c ) + A 22 f 2 ( Tw,c , u 2 , t) + A 21 g 2 ( Tw,e , t)u 1 u 1 =    1, for u d 1 > 1 0, for u d 1 < 0 u d 1 else. ( 61
)
Remark 3. k c , in the desired high-pressure dynamics (60), is a design parameter of the supervisory control. A large value of k c should be used to maintain p as close to p max as possible. However, a lower value of k c reduces the difficulties in solving the DRTO based on the above closed-loop model.

Dynamic real-time optimization

The noncausal (optimal) behavior of the second control input * m A shown by the DP solution detailed in Section 3.3, illustrates the potential benefit of considering the future exhaust gas conditions. As in the context of an ORC WHR system on board a diesel-electric railcar, future driving conditions can be known in advance over a short horizon, a prediction-based strategy -namely dynamic real-time optimization (DRTO) -is explored here. The formulation of the DRTO problem is similar to the moving horizon formulation used in nonlinear model-predictive control, although an economic objective is chosen to provide economically optimal operation at all times (see e.g Würth et al. [START_REF] Würth | Neighboring-extremal updates for nonlinear model-predictive control and dynamic real-time optimization[END_REF]). Here, the economic objective is the net energy production along a moving horizon and the problem can be defined as follows: max

u j 2 (t) t j f t j 0 P net (p, u), dt (62) s
.t. p(t j 0 ) = pj (63)
where the (closed-loop) dynamics of p is defined by system (61) and P net by equations ( 2)-( 5).

The optimization problem is solved on the time horizon t j 0 , t j f (on which driving conditions are assumed to be perfectly predicted) at each sampling instant t j 0 considering the (measured) value of the process output pj .

Notice that the state constraint ( 15) is directly managed by the closed-loop model (61). On the other hand, input constraints ( 16) are handled as described in Dufour et al. [START_REF] Dufour | On nonlinear distributed parameter model predictive control strategy: on-line calculation time reduction and application to an experimental drying process[END_REF], that is by defining a new unconstrained decision variable d j using the invertible transformation:

u j 2 = g(d j ) = g moy + g amp tanh d j -g moy g amp , (64)
where g moy = (u max 2 -u min 2 ) /2 and g amp = (u max 2 + u min 2 ) /2. The pressures trajectory, resulting from DRTO solution, provides the set-point of a low-level controller as illustrated in Figure 10. This trajectory planning is implemented on the current sampling interval, and the optimization horizon is then shifted by the sampling interval ∆t. Notice that a low-level superheating controller -that is not presented here for the sake of clarity -must also be implemented. Details on the pressure and superheating low-level controllers used in this study can be found in Peralez et al. [START_REF] Peralez | Improving the control performance of an organic rankine cycle system for waste heat recovery from a heavy-duty diesel engine using a model-based approach[END_REF] and Peralez et al. [START_REF] Peralez | Control design for an automotive turbine rankine cycle system based on nonlinear state estimation[END_REF]. The relevance of the DRTO method is shown in Section 4.4. The benefits of using a dynamic model for supervisory control are evaluated by comparing DRTO with a simpler control strategy based on static optimization, akin to those presented in the studies mentioned in the beginning of the paper (Quoilin et al. [START_REF] Quoilin | Thermo-economic optimization of waste heat recovery organic rankine cycles[END_REF][START_REF] Quoilin | Dynamic modeling and optimal control strategy of waste heat recovery Organic Rankine Cycles[END_REF], Manente et al. [START_REF] Manente | An organic rankine cycle off-design model for the search of the optimal control strategy[END_REF]). In static optimization, a quasi-steady-state model derived from the dynamic model (Section 2.2), where time derivatives are set to zero, is used.

Application of dynamic real-time optimization

The supervisory control structure of DRTO illustrated in Figure 10 is applied to a realistic Rankine system simulator (based on the moving-boundary principle, see e.g Peralez et al. [START_REF] Peralez | Improving the control performance of an organic rankine cycle system for waste heat recovery from a heavy-duty diesel engine using a model-based approach[END_REF]). The method described in Section 4.3 is implemented with a moving horizon length of 60s and sampling time values of ∆t = 1s and ∆t = 0.02s. Notice that the chosen value for ∆t remains small when compared to the response time of the simplified model used for the DRTO.

The benefits of using a dynamic model for the supervisory control are evaluated by comparing DRTO with the simpler control strategy based on static optimization. It is shown (Figure 11) that the increase in the amount of recovered energy J is of about 7%.

In order to demonstrate the relevance of the splitting approach introduced in section 4.2, the DRTO problem with original 2 × 2 MIMO model formulation was also solved. It resulted in a similar level of accuracy (with an increase of about 0.1% for J), whereas the computational times were dramatically increased. Thus, for the splitting approach the average and maximum runtime per step was found to be of 0.18 s and 0.32 s whereas it was of 0.43 s and 0.90 s for the MIMO formulation. This computing performance, guarantees the real-time use of the proposed algorithm.

Notice that state variables and control inputs trajectories resulting from the -closed loop -DRTO (Figure 10) are close to the -ideal -DP solution (Figure 8). Because the control model is not perfect (i.e it is different from the simulation model) unavoidable differences in the results remain.

Conclusion

In this study, an original optimal control problem is formulated in order to find the maximum power produced by an organic Rankine cycle waste heat recovery system on board a diesel-electric railcar, along a given mission profile. The development of an ad-hoc reduced-order dynamic model, validated on experimental data, makes solving the optimization problem computationally possible using a novel version of the state-of-the-art level-set DP algorithm, based on an adaptive grid.

The analysis of DP results leads to a dynamic real-time optimization problem with a single decision variable, which can be solved online and implemented as a supervisory control for the ORC system. Although providing theoretical proof of its stability is a difficult task, a future experimental implementation could demonstrate the relevance of the proposed strategy. However, the results obtained when it is applied to a realistic Rankine system simulator show a significant increase in the amount of recovered energy with respect to the benchmark provided by a simpler (but state-of-the-art) supervisory control based on static optimization.

Finally, the offline and online control strategies presented in this article are, in the author's opinion, of general interest, though they have been designed for an ORC with some distinctive -more or less favorable -features, namely in terms of cooling system, power management and mission predictability. Their application to other vehicular ORCs, on board road vehicles for instance, appears perfectly possible, though some adaptations may be needed.

Figure 1 :

 1 Figure 1: ORC system under investigation. Three control inputs are available: pump mass flow * m, air mass flow

Figure 2 :

 2 Figure 2: Input disturbances for ORC problem.

Figure 3 :

 3 Figure 3: Rankine cycle for a dry fluid

Figure 4 :

 4 Figure 4: p 1 , p 3 and * m maps, function of Tw,e and Tw,c

( a )

 a External conditions. (b) Model prediction and measurement.

Figure 5 :

 5 Figure 5: Single-state evaporator model: experimental validation

Figure 6 :

 6 Figure 6: Numerical issues near boundary surface.

Figure 7 :

 7 Figure 7: A first estimation I k of level-set function is computed on the whole state space (left plot). A second discretization (right plot) allows to improve the accuracy of both level-set and cost-to-go functions.

Figure 8 :

 8 Figure 8: Adaptive-grid DP solution: time evolution of the ORC system characteristics

Figure 9 :

 9 Figure 9: Influence of the equivalent turbine nozzle surface S in term of mean net power recovered Pnet .

Figure 10 :

 10 Figure 10: Two-level closed-loop control strategy. DRTO provides the ORC pressures set-point p ref (with a sampling time ∆t) to the low-level controller which computes the control values -corresponding to evaporator by-pass and air fan actuators -at a (faster) sampling time ∆t.

Figure 11 :

 11 Figure 11: Result of dynamic real-time optimization (solid blue lines) compared with static optimization approach (dash green lines).

Figure A. 13 :

 13 Figure A.13: Optimal trajectory (x o , v o ) and evolution of backward-reachable space (evaluated by adaptive-grid DP).

Figure A. 14 :

 14 Figure A.14: Algorithm accuracies vs number of model calls for the simple energy management example.

Part of this study, namely the offline approach, has been partially presented in the conference paper Peralez et al.[START_REF] Peralez | Optimal control of a vehicular organic rankine cycle via dynamic programming with adaptive discretization grid[END_REF]

Acknowledgments

The authors wish to thank ANR, the French National Research Agency, for partially funding this work within the framework of the project TRENERGY ("Train Energy Efficiency via Rankinecycle Exhaust Gas Heat Recovery", # ANR-12-VPTT-09), and ALSTOM Transport, for providing application data.

Appendix A. A simple energy management example

To illustrate the benefits of the new algorithm, a first optimal control problem with an analytic solution is studied. It concerns energy management for an electric vehicle.

Appendix A.1. System description

Neglecting the aerodynamic friction and considering the road slope null, the simplified continuoustime dynamics of the vehicle in the longitudinal direction x can be written as in Petit and Sciarretta [START_REF] Petit | Optimal drive of electric vehicles using an inversion-based trajectory generation approach[END_REF]:

where the control variable u is a percent torque demand, h 0 and h 1 are constant parameters depending on the rolling resistance coefficient, the motor maximum torque, vehicle mass, transmission ratio and wheel radius.

The on-board electric power consumption can be written as

. Optimal control problem

The aim is to find a control strategy that minimizes the power consumption under the constraints that the vehicle must reach a destination point at a distance D in a given time t f starting from a given point, at rest. In order to formulate this problem for DP solving, the continuoustime model is discretized using an Euler forward approximation with a time step T s . The optimal control problem can be formulated as follows:

where

Ts . The final constraints are formulated in equations (A.8). In the following, a particular solution is sought for a final distance D = 200 and a final time t f = 60. The time discretization is chosen to be T s = 0.2s (then N = 300), with x = 205, v = 8, v = 0.3.

Appendix A.3. DP solution

The problem stated above is solved using the basic, the original level-set and the level-set with two-step discretization DP algorithms.

State variable spaces are first discretized respectively with N v = 81 and N x = 106, while the control space is discretized with N u = 81. Figure A.12 shows the resulting system trajectory and the corresponding control inputs for each of the three algorithms. They are compared to the analytic solution (a parabola) presented by Petit and Sciarretta [START_REF] Petit | Optimal drive of electric vehicles using an inversion-based trajectory generation approach[END_REF]. The results obtained using the new algorithm are closer to the analytic solution, even with a relatively coarse discretization. In term of cost, the relative mean errors of the three algorithms are found to be of 1.4%, 0.78% and 0.19% respectively.

The evolution of backward-reachable space boundaries, where level-set function vanishes, helps to analyze those results. Figure A.13 shows that in the first backward steps, the optimal trajectory (x o , v o ) is close to Ĩk estimations. Here, for accurate estimation of boundaries it is important not to consider reachable points close to the optimal solution as infeasible. After those first backward steps, the optimal trajectory moves away from the boundaries. However, the reachable state space remain restricted. Thus, the second discretization step of the proposed algorithm can improve the accuracy of cost-to-go function.

Two-step discretization involves a larger computational effort for a given set of parameters (N v , N x , N u ). This computation cost can be expressed in terms of number of model calls for evaluation. For example, for the level of discretization used above, basic and original level-set DP both require about 2.085


Mathieu Grandemange

Olivier Cadot

A Courbois

Vincent Herbert

Denis Ricot

Tony Ruiz

Rémi Vigneron

A study of wake effects on the drag of the Ahmed's squareback model at the industrial scale

Keywords: Aerodynamics, separation, bi-stability, control, optimisation, induced drag

d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Car aerodynamics, and more generally ground vehicle aerodynamics is of major interest for future energy saving developments. Due to severe constraints, such as functional shapes, security, and comfort, vehicles have a high aerodynamics drag [START_REF] Hucho | Aerodynamics of road vehicles[END_REF][START_REF] Choi | Aerodynamics of heavy vehicles[END_REF]. They belong to the so called bluff bodies associated with strong flow separation and unsteady wake. There are therefore very good reasons to address this problem fundamentally. The work of Ahmed et al. [START_REF] Ahmed | Some salient features of the time-averaged ground vehicle wake[END_REF] describes the effect of the afterbody shape of a simplified road vehicle on the flow topology and on the drag. In addition to the large recirculating region of low pressure due to the full flow separation at the base, their work proves the critical influence of the slant angle of the rear window. A drag reduction of 8% compared to the squareback configuration (0 • case) is obtained for a slant inclination of 12.5 • . Beyond this optimal value the drag increases continuously and the worst case (slant angle close to 30 • ) is associated with a 50% increase in drag compared to the squareback geometry. Over a critical angle close to 30 • , the drag decreases abruptly to recover values close to the 0 • configuration. In terms of flow topologies, the blunt afterbody is responsible for a massive separation on the base. For slant angles below 30 • , the flow remains attached on the rear window (at least partially) and a pair of counter-rotating vortices develops from the sides of the slanted face [START_REF] Lienhart | Flow and turbulence structure in the wake of a simplified car model[END_REF][START_REF] Beaudoin | Cavitation as a complementary tool for automotive aerodynamics[END_REF][START_REF] Krajnović | Flow around a simplified car, part 1: large eddy simulation[END_REF][START_REF] Krajnović | Flow around a simplified car. Part 2. Understanding the flow[END_REF]: the intensity of these streamwise vortices gradually grows as the slant angle increases and consequently the drag rises. Above 30 • , the flow separates upstream of the rear window and a massive recirculation region without energetic streamwise vortices, similar to the 0 • case, is established. The flow properties of the squareback geometry are investigated both numerically and experimentally in [START_REF] Duell | Experimental study of a ground vehicle body unsteady near wake[END_REF][START_REF] Bayraktar | Experimental and computational investigation of Ahmed body for ground vehicle aerodynamics[END_REF][START_REF] Krajnović | Numerical study of the flow around a bus-shaped body[END_REF][START_REF] Grandemange | Reflectional symmetry breaking of the separated flow over threedimensional bluff bodies[END_REF][START_REF] Grandemange | Turbulent wake past a three-dimensional blunt body. Part 1. Global modes and bi-stability[END_REF]. The geometry used in the experiments of Ahmed et al. [START_REF] Ahmed | Some salient features of the time-averaged ground vehicle wake[END_REF] is often used to explore drag reduction strategies strived for industry. Thus, literature reports the use of various different control devices. Much work has been devoted to the 25 • slant angle configuration; significant drag reductions are achieved through passive or active control strategies such as splitter plates [START_REF] Gilliéron | Aerodynamic drag reduction by vertical splitter plates[END_REF], flaps [START_REF] Beaudoin | Drag and lift reduction of a three-dimensional bluff body using flaps[END_REF][START_REF] Fourrié | Bluff body drag reduction using a deflector[END_REF], boundary layer streaks [START_REF] Pujals | Drag reduction of a three-dimensional bluff body using coherent streamwise streaks[END_REF] or even pulsed jets [START_REF] Bruneau | Active procedures to control the flow past the ahmed body with a 25 • rear window[END_REF]. Diverse strategies equally provide interesting base pressure recovery in the squareback case: splitter plates [START_REF] Khalighi | Experimental and computational study of unsteady wake flow behind a bluff body with a drag reduction device[END_REF], porous devices [START_REF] Bruneau | Coupling active and passive techniques to control the flow past the square back ahmed body[END_REF] or active control [START_REF] Littlewood | Aerodynamic drag reduction of a simplified squareback vehicle using steady blowing[END_REF][START_REF] Barros | Control of three-dimensionnal blunt body wake using low and high frequency pulsed jets[END_REF][START_REF] Rouméas | Analysis and control of the near-wake flow over a squareback geometry[END_REF].

The experiments of Littlewood & Passmore [START_REF] Littlewood | The optimization of roof trailing edge geometry of a simple squareback[END_REF] depict the effect of a small chamfer at the upper trailing edge of a squareback geometry. The geometry of the after body is similar to the squareback Ahmed configuration. Despite the small size of the chamfer in the streamwise direction (only 15% of the base height), an optimal drag reduction of 4.4% is obtained for a chamfer angle of 12 • . The parabolic-like dependence between the chamfer angle and the drag is reminiscent of the results of Ahmed et al. [START_REF] Ahmed | Some salient features of the time-averaged ground vehicle wake[END_REF] for moderate slant angles where an optimal 8% drag reduction is obtained for an angle of 12.5 • . Thus, the significant drag reduction obtained by Littlewood & Passmore [START_REF] Littlewood | The optimization of roof trailing edge geometry of a simple squareback[END_REF] with a small chamfer highlights the high impact on force of the flow orientation at the trailing edge. Similarly, Grandemange et al. [START_REF] Grandemange | Effect on drag of the flow orientation at the base separation of a simplified blunt road vehicle[END_REF] studied the effect of the angle of small flaps (or spoilers) located at the top and the bottom of the squareback model trailing edge. Their parametric study led to an optimized configuration, with an angle Φ T = 9.2 • for the top flap and Φ B = -7.4 • for the bottom flap. The drag reduction mechanism was explained by a simple two dimensional effect of wake thinning, hence reducing the bluffness of the body [START_REF] Roshko | Perspectives on bluff body aerodynamics[END_REF]. More importantly, the presence of the optimum was related to a drag penalization because of the presence of longitudinal vortices associated with 3D separation at the salient lateral edges of the flaps. This approach using chamfered shapes is particularly interesting since it corresponds to the method applied by the car manufacturers on real vehicles to limit the drag. Indeed, the angles of the spoiler and of the diffuser (when one is present) are set empirically to reach the optimum drag point given the vehicle parameters (global shape, ground clearance...), the characteristic length of these devices remaining small compared to the vehicle size.

Recently, Grandemange et al. [START_REF] Grandemange | Reflectional symmetry breaking of the separated flow over threedimensional bluff bodies[END_REF][START_REF] Grandemange | Turbulent wake past a three-dimensional blunt body. Part 1. Global modes and bi-stability[END_REF] found that the Ahmed body presents bi-stable wake positions related to a reflectional symmetry breaking in the flow solution. This bi-stability creates large cross flow force fluctuations and probably additional drag. The existence of bi-stability is not only a fact of the Ahmed body, it is related to the aspect ratio of the rectangular blunt trailing edge of the 3D body together with the ground clearance [START_REF] Grandemange | Bi-stability in the turbulent wake past parallelepiped bodies with various aspect ratios and wall effects[END_REF]. For instance, it has has been observed for a double backward facing step geometry [START_REF] Herry | Flow bi-stability downstream of three-dimensional double backward facing steps at zero-degree slideslip[END_REF]. In car aerodynamics, Lawson et al. [START_REF] Lawson | An investigation of the flow characteristics in the bootdeck region of a scale model notchback saloon vehicle[END_REF] observed such behaviour for a realistic car model at 1/4 scale. For the Ahmed body, the bi-stability was found particularly sensitive to steady disturbances when placed in the middle of the recirculation bubble. A vertical cylinder, having the body height and only 5.6% the body width in diameter is able to suppress the bi-stability together with producing a drag reduction estimated to be between 5% and 8% [START_REF] Grandemange | Turbulent wake past a three-dimensional blunt body. part 2. experimental sensitivity analysis[END_REF].

The goal of this article is to investigate these new flow understandings and their corresponding drag improvements at the industrial scale while considering two main questions: (i) can the substantial drag contribution of the longitudinal vortices near the lateral edges of the chamfers be evidenced ? 2

(ii) can the substantial drag contribution due to the wake bi-stability be evidenced ?

The simplified Ahmed body is still employed, but in real flow condition with road effect and rotating wheels. The article is organized as follows. The experimental set-up is first described in section 2. The parametric study of small chamfers at the top and bottom of the trailing edge in section 3.1 leads to an optimal drag configuration. This optimal drag model is then used as a reference to test drag reduction methods in the next two sections. In section 3.2, the bi-stable behavior of the wake is evidenced through sideslips measurements and then controlled to successfully reduce the drag. In section 3.3, the longitudinal vortices are suppressed by introducing chamfers on the lateral sides and proves their important drag contribution to the optimal configuration with top and bottom chamfers. Finally, the concluding remarks are presented in Sec. 4.

Experimental setup

Figure 1 illustrates the actual set-up for the experiments. It is similar to the squareback Ahmed gemetry used in [START_REF] Ahmed | Some salient features of the time-averaged ground vehicle wake[END_REF], but is 4 times larger, has a chamfered appendix (shown at (ii) in the diagram) and has a top rear edge chamfered at an angle Φ T > 0 (shown at (i) in the diagram). Also, the original cylindrical supports are replaced by proper wheels, and the distance to the floor is larger. In [START_REF] Ahmed | Some salient features of the time-averaged ground vehicle wake[END_REF], the separating distance is 0.186 the body height while it is 0.165 in the present study. The body is tested in the full scale wind-tunnel of GIE S2A at Montigny-Le-Bretonneux especially designed for automotive aerodynamics. A full description of the facility may be found in [START_REF] Waudby-Smith | The GIE S2A full-scale aeroacoustic wind tunnel[END_REF]. Briefly, the test section is a 3/4 open jet with a cross-section of 24 m 2 . The body is placed on a large turntable to reproduce sideslip situations with a yaw angle β as defined in figure 1(b). A wide belt is entrained at the wind velocity between the rotating wheels to introduce a real road effect on the underbody flow.

The measurements systems are provided by the GIE S2A. The model is fixed by four thin supports (two downstream the front wheels, and two upstream the rear wheels) on a 6 components balance which provides both the mean value and the standard deviation of each component of the force, say F i and Std(F i ) with i ∈ {x, y, z} denotes respectively, the drag, the side and the lift force component. These quantities are obtained from recordings at a sampling frequency of 10 Hz during 40 s when a statistically steady state is reached; part of the fluctuations are ascribed to measurement uncertainty but the values of Std(F i ) still give qualitative information on the steadiness of the forces. Parietal static pressure are obtained on the geometry at 47 locations: 35 taps are placed on the base as presented in figure 1(c), the 12 others are distributed in the middle of the top, bottom and lateral faces at x * = -0.22, -0.43 and -1.30. When additional information is needed, other pressure sensors are set at some precise locations. Standard pressure measurements are obtained from signal sampled at 2.5 Hz, usually averaged over 40 s.

Measurements in the wake use three 18 holes spherical probes mounted on a three-axis traversing system [START_REF] Waudby-Smith | The GIE S2A full-scale aeroacoustic wind tunnel[END_REF]. The probes are "omniprobe" models from Aeroprobe Corporation having diameters of 9.53 mm. The measurement is performed continuously at an acquisition rate of 2.5 Hz while the traversing system is moving at a constant speed of 75 mm.s -1 . Hence, one measurement point is obtained each 30 mm that sets the spatial resolution. It gives the mean velocities U x , U y and U z and its gradients after processing as well as the static pressure. The precision of the velocity is given by the manufacturer to be less than 1% of the measured velocity modulus, and the angular precision to be less than 0.4 • . The automatic displacement allows to get these data in planes iso-x, iso-y or iso-z. O sets the origin of the coordinate system. The blue dots locate the pressure taps on the base. For the reference squareback Ahmed geometry (SB) (see table 1), there is no appendix (ii) and

Φ T = 0 • in (i).
For all the tests, the inlet velocity is set at U 0 = 33.3 m s -1 . The Reynolds number based on the height H of the geometry (see figure 1a) is Re = U 0 H ν = 2.5 10 6 . In the following, a * denotes the non-dimensional value of any quantity a(x, y, z, t) made dimensionless by a combination of the height H and the inlet velocity U 0 . The force and the pressure coefficients are defined as :

C i = F i 1 2 ρS U 2 0 ; C p = p -p 0 1 2 ρU 2 0 ,
with i referring to the force component, ρ is the air density, S = H × W is the frontal area of the Ahmed body and p 0 is the free stream pressure. The coefficients are given to be accurate in the absolute range ±0.002.

Results

The force coefficients of the squareback Ahmed model are displayed in table 1. The drag coefficient is 0.309 which is larger than the value of 0.250 obtained in [START_REF] Ahmed | Some salient features of the time-averaged ground vehicle wake[END_REF] despite the fact that the Reynolds numbers of our study is only twice the Reynolds number in [START_REF] Ahmed | Some salient features of the time-averaged ground vehicle wake[END_REF]. The larger drag is then ascribed to the differences in the ground clearances, the rotating wheels and the road effect. The lift coefficient C z is slightly negative. The fluctuation of the side force coefficient S td(C y ) is very large compared to the values of the fluctuations in lift and drag. It is related to the wake bi-stability [START_REF] Grandemange | Reflectional symmetry breaking of the separated flow over threedimensional bluff bodies[END_REF][START_REF] Grandemange | Turbulent wake past a three-dimensional blunt body. Part 1. Global modes and bi-stability[END_REF] as it will be shown and discussed in section 3.2.

Drag optimisation by flow orientation at the top and the bottom of the base

This section studies at industrial scales the drag dependence on the flow orientation at the base presented in [START_REF] Grandemange | Effect on drag of the flow orientation at the base separation of a simplified blunt road vehicle[END_REF]. Our objective here is to obtain an optimized Ahmed model with substantial drag reduction compared to that of the squareback geometry by chamfering the trailing edge. The top edge of the squareback model is chamfered with positive or negative angles as indicated in the dashed box (i) in figure 1. For practical reasons, we had to introduce an appendix at the bottom edge of the squareback model in order to orientate the flow as indicated in the dashed box (ii) in figure 1. The geometry of the chamfers is a circle arc with an horizontal tangent at the upstream edge. The angles Φ T and Φ B are defined by the line passing through the upstream and the downstream edge of the chamfers as depicted in figure 1(a). As a first step, the bottom angle is set to ϕ B = 0 • and the effects of the top angle are considered. The forces as functions of ϕ T are presented in figure 2. They confirm the trends and amplitudes obtained in [START_REF] Grandemange | Effect on drag of the flow orientation at the base separation of a simplified blunt road vehicle[END_REF] as the drag is measured as a 2 nd order polynomial fit while the lift presents an affine relationship. Besides, the rear lift1 , denoted C zr , is equally an affine function of the top slant angle corroborating the fact that the lift evolutions are mostly ascribed to pressure modification on the afterbody. However, one can see a slight difference of slope between C z and C zr in figure 2; this indicates that there are still some small variations of the lift on the front axle.

The pressure distribution on the afterbody in the plane y * = 0 are presented in figure 3 for different top angles ϕ T . As in [START_REF] Grandemange | Effect on drag of the flow orientation at the base separation of a simplified blunt road vehicle[END_REF], the effect is clear on the pressure distribution on the top face. The curvature of the geometry is transferred to the potential flow. This results in increases or decreases in pressure depending on the sign of ϕ T , just upstream of the slanted surface for x * > -1.0. The effect is also clear on the base pressure, especially on the gradient in the z direction. Finally, one can observe some modifications of the pressure on the lower face evidencing the coupling effects of the top and bottom angles.

The quadratic dependence between the drag and ϕ T combined with the affine dependence between the lift and ϕ T leads to a relationship between the drag and the lift that can be expressed by

C x -C x0 = α(C z -C z0 ) 2 .
(

) 1
These relationships are presented in figure 4(a) for different bottom angles ϕ B . As in the experiments in [START_REF] Grandemange | Effect on drag of the flow orientation at the base separation of a simplified blunt road vehicle[END_REF], all the parabole are self-similar with a parameter α = 1.27 independent of ϕ B at first order (see figure 4b). It is remarkable to see that the value of α is identical to the one obtained in [START_REF] Grandemange | Effect on drag of the flow orientation at the base separation of a simplified blunt road vehicle[END_REF] which proves that the force evolutions are of similar amplitude. Moreover, this is consistent with the geometrical interpretation of α in the theory of induced drag.

These variations of force intensities are associated with important modifications of the flow topology. The flows obtained with the two extreme lift configurations are compared to the flow obtained for ϕ T = ϕ B = 0 • in figures 5, 6 and 7. In the plane y * = 0, the results plotted in figure 5 show that the top -bottom equilibrium of the flow is strongly modified. The recirculation bubble is almost symmetric for ϕ T = ϕ B = 0 • while it is significantly asymmetric in the extreme cases. So, the different couples of angles considered in these experiments should explore all the wake orientations between figures 5(a) and 5(c).

The three-dimensional effects play a first order role in the drag. Figure 6 presents the streamwise vorticity in the plane x * = 0.43 for the three different configurations. These maps confirm the formation of streamwise vortices downstream of the corners of the base. The dashed lines locate the contour U x * = 0.5 to discriminate the energetic vortices in the high velocity flow from the chaotic vorticity obtained in the recirculation bubble, mostly due to measurement uncertainty. Besides, for the case ϕ T = ϕ B = 0 • in figure 6(b), there are clear vortices downstream of the underbody in spite of the absence of bottom inclination; they may be ascribed to the presence of rotating wheels.

Further downstream, at x * = 2.60, figures 7(a) and 7(c) evidence that the concentrated vortices merge from each side to form a single pair of large scale counter-rotating vortices for the configurations (ϕ T = -5 • , ϕ B = -10 • ) and (ϕ T = 10 • , ϕ B = 12 • ). The presence of these vortices and their directions of rotation are consistent with the lift measured at C z = -0.22, C z = -0.05 and C z = 0.20 for the configurations (ϕ

T = -5 • , ϕ B = -10 • ), (ϕ T = 0 • , ϕ B = 0 • ) and (ϕ T = 10 • , ϕ B = 12 • ) respectively.
As a consequence, the general trends and topologies ascribed to the different angles reported in [START_REF] Grandemange | Effect on drag of the flow orientation at the base separation of a simplified blunt road vehicle[END_REF] remain valid at higher Reynolds numbers. The optimal drag configuration obtained for ϕ T = 7.5 • and ϕ B = -5.0 • (O model in table 1) with a 5.8% drag reduction is now considered. The corresponding flow is presented in figure 8. The velocities in the plane y * = 0 show that the flow is slightly asymmetric with a reduced recirculation region in comparison to the case ϕ T = ϕ B = 0 • (see figure 8a). The streamwise vorticity in the very near wake is visible in figure 8(b); it evidences the presence of counter-rotating vortices downstream of the corners of the flaps. The vortices emitted from the same side of the geometry are very likely to attenuate each other which results in the absence of large-scale vortices further downstream: only residual traces of vortices are visible in the streamlines in figure 8(c).

As in the model of drag sources [START_REF] Grandemange | Effect on drag of the flow orientation at the base separation of a simplified blunt road vehicle[END_REF] deduced from the drag optimization using orientated flaps at the top and bottom trailing edge, the optimized geometry results in a compromise between the reduced bluffness at the rear and the formation of streamwise vorticity at both ends of the chamfers. This can be seen by comparing the natural case in figure 6(b) to the optimized geometry in figure 8(b). The height of the separated region depicted with the dashed line is smaller for the optimized case, but the streamwise vortices are more intense. The idea in [START_REF] Grandemange | Effect on drag of the flow orientation at the base separation of a simplified blunt road vehicle[END_REF] is to ascribe a source of drag decrease by the wake thinning, which takes into account two dimensional effects only, and a source of drag increase resulting from the drag induced by the streamwise vortices which takes into account three-dimensional effects. The latter depends on the square of the lift force exerted on the chamfers as for the finite size effect of wings in aeronautics.

To gain more understanding of the flow, control experiments explore both the sensitivity of this optimized configuration, especially regarding a possible bi-stable behavior in section 3.2 and some perspectives to limit the formation of the streamwise vortices due to the chamfered trailing edge in section 3.3.

Sideslip effect, wake bi-stability and its control

One of the capabilities of this wind tunnel is the possibility to set a slight and accurate yaw angle; this provides a quantitative parameter of the setup asymmetry. The force measurements allow a detailed description of the effects of a small yaw angle. The drag and the fluctuations of forces as functions of β are plotted in figure 9. The drag evolution can be decomposed into two regions. There is an affine drag increase as the modulus of β increases for | β |> 0.5 • whereas the drag saturates (given the accuracy of the measurements) on a plateau in the range | β |< 0.5 • . Note that the linear trend is equally reported in the experiments of [START_REF] Gohlke | Experimental analysis of flow structures and forces on a threedimensional bluff body in constant cross-wind[END_REF] and is likely to extend up to |β| ≈ 15 • . The fluctuations of drag and lift are independent of the yaw angle but the side force presents strong fluctuations for | β |< 0.5 • in correspondence to the region of the saturated drag value. These fluctuations are governed by the coexistence of two mirror states (bi-stability) in the wake as evidenced in [START_REF] Grandemange | Turbulent wake past a three-dimensional blunt body. Part 1. Global modes and bi-stability[END_REF]. The pressure distributions on the base for β = -0.7 • , 0 • and 0.8 • are plotted in figure 10. Both extreme angles present a predominance of one of the two asymmetric states, using the same terminology as in [START_REF] Grandemange | Turbulent wake past a three-dimensional blunt body. Part 1. Global modes and bi-stability[END_REF], the state #N is predominant for the negative yaw angle in figure 10(a) and state #P is predominant for the positive yaw angle in figure 10(c). For the nil yaw angle, the residual asymmetry in figure 10(b) indicates that both states are not equiprobable, but with a small predominance of state #P. Hence, the asymmetry introduced by a slight yaw angle β modifies the states proportion #P and #N. For yaw angles such that | β |> 1 • , the wake is locked in one state which explains the low and constant side force fluctuation. This is similar to the wake manipulation in [START_REF] Grandemange | Turbulent wake past a three-dimensional blunt body. part 2. experimental sensitivity analysis[END_REF] using a small control cylinder that is able to change the proportion of each state #P and #N keeping a constant drag of the Ahmed body. In the present case, the drag is also insensitive to the states proportion because of the presence of the plateau in the range of small yaw angles. Now, as in the laboratory experiments of [START_REF] Grandemange | Turbulent wake past a three-dimensional blunt body. part 2. experimental sensitivity analysis[END_REF], the drag associated with the bi-stability may be quantified using a vertical disturbance. They showed that a vertical cylinder placed in the middle of the recirculation region suppresses the bi-stability which enables a drag reduction to be achieved. To confirm the sensitivity of the asymmetric states to such a disturbance, a similar experiment is performed here at the industrial scale. A vertical control cylinder is placed at x C * = 0.87 and y C * = 0; it has a square cross section with a length matching the height of the body and a side of 50 mm, i.e. 4.3% of the body height. Practically, it is supported by a rod of 20 mm in diameter, fixed on the wind tunnel floor. As in [START_REF] Grandemange | Turbulent wake past a three-dimensional blunt body. part 2. experimental sensitivity analysis[END_REF], the drag of the control cylinder is not measured. The control cylinder position is at the end of the sensitive region obtained in [START_REF] Grandemange | Turbulent wake past a three-dimensional blunt body. part 2. experimental sensitivity analysis[END_REF]. The associated force measurements are given in table 1 and plotted in figure 9 with the square symbols for β = 0 • . The drag is slightly decreased, by approximately 1.7% compared to the optimized model, but a significant effect is obtained on the fluctuations of the side force. The levels of Std(C y ) are reduced from 0.015 to 0.007, the latter value being close to the expected value for a stable configuration. A smaller but significant effect is reported using a smaller control cylinder of diameter 20 mm, i.e. 1.7% of the body height, even if this disturbance is much less efficient in stabilizing the wake and the drag reduction hardly measurable (see figure 9).

One can note that the use of a disturbance leads to the minimum drag point presumed from the linear dependencies on the yaw angle for | β |> 1.0 • . However, the drag reductions obtained in these experiments after suppression of the bi-stability are clearly smaller than the 5 to 8% expected from the results of [START_REF] Grandemange | Turbulent wake past a three-dimensional blunt body. part 2. experimental sensitivity analysis[END_REF]. This could be ascribed to a non-optimal position of the disturbance 2 . However, we do not rule out that the disturbance may also produce other effects, not yet identified, independent from the bi-stability suppression and sensitive to the Reynolds number. Further analysis is needed for definite conclusions.

Flow orientation at the four sides of the base with cavity

The squareback trailing edge is now slightly extended downstream (by a length of 100 mm) as depicted in figure 11. Four configurations are studied, they are listed in table 2.

For the first configuration denoted C, the extension is straight on the four sides. A drag reduction of 3.6% is obtained with C x = 0.298 for the C model against C x = 0.309 for the squareback Ahmed body (table 1). It is known that deep cavities at the base improve the drag [START_REF] Kruiswyk | Effect of a base cavity on subsonic near-wake flow[END_REF], which is particularly employed as a major drag reducer for trucks [START_REF] Browand | Fuel savings by means of flaps attached to the base of a trailer: field test results[END_REF]. The second configuration, called C-TB, uses the flow orientation at the top and bottom trailing edge. Both the top and the bottom extensions are thus chamfered as in the previous part, but with larger angles, of Φ T = -Φ B = 10 • , than the optimized Ahmed body. For these larger chamfers, the C-TB geometry has no drag improvement (see table 2) since it is not the optimum. Particularly, we can see the presence of longitudinal vortices in figure 12, that are associated with three dimensional separations produced by global pressure differences between the body faces. The pressure distribution on the sides of the geometry at x * = -0.22 presented in figure 13 supplies a measurement of these global pressure levels on each face. We can see in figure 13(a) that the pressure levels on the C-TB model are only affected at the top and bottom faces compared to the C model. Because of the curvature introduced by the flow orientation at the edges, the pressure level on the top and bottom faces is lower than that on the left and right faces. The pressure differences are responsible for the three dimensional separation and vortex formation around each corner at the base. 2 for the different models configuration and force coefficient measurements. The arrows locate the pressure taps; the pressure coefficient is always negative and quantified by the distance from the surface (the length for C p = 1 is indicated) and the dotted gray lines locate the contours C p = -0.2.

We turn now to the C-LR model for which only both lateral edges are chamfered with Φ L = -Φ R = 10 • (see table 2). In that case, only the pressure levels on the left and right faces are significantly decreased compared to the C model (figure 13b). It results in three dimensional separation producing longitudinal vortices at each corners but with opposite vorticity to that of the C-TB model (see figure 12b). A drag reduction of 1.7% is obtained for the C-LR model.

The very interesting result is for the flow orientation at all the edges Φ T = -Φ B = 10 • and Φ L = -Φ R = 10 • , better known as the "boat tail" geometry. The C-TBLR model now combines the two previous geometries and, as a result also combines the pressure levels as can be noticed in figure 13. The pressure distributions on the top and bottom faces are identical to these of C-TB model, and the left and right pressure distributions are identical to these of the C-LR model. In opposition, the drag reduction is not the sum of the DR of both the C-TB and C-LR models, but is 7.4% that is much larger ! The explanation is given by the equilibrium of the pressure on all the faces due to the boat tail effect, that suppresses the three dimensional separations, and hence the longitudinal vortices as shown in figure 13(c). This result is consistent with the model proposed in [START_REF] Grandemange | Effect on drag of the flow orientation at the base separation of a simplified blunt road vehicle[END_REF] because it proves the role of the induced drag in the presence of the drag optimum when the wake flow is oriented on two opposite sides only. Eventually, the drag reduction obtained with the C-TBLR model is 10.7% compared to the reference squareback Ahmed model.

Concluding remarks

The control of the wake past the squareback Ahmed geometry at Re H = 2.5 10 6 extends the results of the previous laboratory experiments to industrial scales. First, the flow presents a dependence on the angles of the trailing edges similar to the one depicted in [START_REF] Grandemange | Effect on drag of the flow orientation at the base separation of a simplified blunt road vehicle[END_REF] both in terms of trends and amplitude. The drag as a function of the top and bottom angles can be analyzed as a two dimensional contribution related to the apparent rear bluffness of the geometry and three-dimensional contributions from both the local vortices emitted at the corners of the base and the large-scale vortices associated with the lift. An optimized geometry is found presenting a 5.8% drag reduction compared to the squareback Ahmed model. The bi-stable behavior of the wake of the reference geometry [START_REF] Grandemange | Turbulent wake past a three-dimensional blunt body. Part 1. Global modes and bi-stability[END_REF] is still observed in the present optimized configuration. It is characterized by a plateau of drag for very small yaw angles of sideslip, typically for |β| < 1.0 • , with important fluctuations of the side force. Based on the sensitivity analysis at the laboratory scale detailed in [START_REF] Grandemange | Turbulent wake past a three-dimensional blunt body. part 2. experimental sensitivity analysis[END_REF], the use of a vertical disturbance at the appropriated location, i.e. near the center of the recirculation region, leads to a suppression of the bi-stable behavior and some additional drag reduction at the industrial scale. For this controlled case the drag reduction is 7.4% compared to the squareback Ahmed model. Hence 1, 6% in DR is ascribed to the presence of the control cylinder.

In parallel, an efficient method to limit the formation of the local longitudinal vortices of the optimized geometry is to orient simultaneously the flow from the sides of the geometry in order to retrieve a relative axisymmetry in the afterbody flow. In that case a drag reduction of 10.7% is obtained compared to the squareback Ahmed model. This drag reduction involves, a cavity effect bringing 3.6% drag reduction, and the boat tail effect that together make the cumulative DR of 10.7%. The cavity effect appears to be an efficient method, in our case it is coupled to a strong reduction of the side force fluctuation indicating an inhibition of the wake bi-stability. However, further experiments are needed to clarify the role of the wake bi-stability in the ability of a deep cavity to reduce the drag of the squareback Ahmed model.

Figure 1 :

 1 Figure 1: Layout of the squareback Ahmed geometry used for the investigations. Side view in shown at (a), top view at (b) and back view at (c) .O sets the origin of the coordinate system. The blue dots locate the pressure taps on the base. For the reference squareback Ahmed geometry (SB) (see table1), there is no appendix (ii) and Φ T = 0 • in (i).

Figure 2 :

 2 Figure 2: Drag (a) and lift (b) for different angles ϕ T with ϕ B = 0 • ; the dashed lines are 1 st and 2 nd order polynomial fits.

Figure 3 :

 3 Figure 3: Pressure coefficient levels on the afterbody in the plane of symmetry for different angles ϕ T with ϕ B = 0 • . The blue dots locate the pressure taps. The pressure coefficient is always negative and quantified by the distance from the surface (the length for C p = 1 is indicated).

Figure 4 :

 4 Figure 4: Drag vs. lift (a) and data centered on their respective extremum (b) as a parametric function of the top slant angle ϕ T for different ϕ B : +, ϕ B = 12 • ; ⃝, ϕ B = 0 • ; ⋄, ϕ B = -5 • ; △, ϕ B = -10 • ; ---, 2 nd order polynomial fits.

Figure 5 :

 5 Figure 5: Velocity in the plane y * = 0 for ϕ T = -5 • and ϕ B = -10 • (a), ϕ T = 0 • and ϕ B = 0 • (b), ϕ T = 10 • and ϕ B = 12 • (c).

Figure 6 :

 6 Figure 6: Streamwise vorticity in the plane x * = 0.43 for ϕ T = -5 • and ϕ B = -10 • (a), ϕ T = 0 • and ϕ B = 0 • (b), ϕ T = 10 • and ϕ B = 12 • (c). The dashed line is the contour U x * = 0.5.

Figure 7 :

 7 Figure 7: Velocity in the plane x * = 2.60 for ϕ T = -5 • and ϕ B = -10 • (a), ϕ T = 0 • and ϕ B = 0 • (b), ϕ T = 10 • and ϕ B = 12 • (c).

Figure 8 :

 8 Figure 8: Flow around the optimal configuration ϕ T = 7.5 • and ϕ B = -5 • (O model in table 1): velocity in the plane y * = 0 (a), streamwise vorticity in the plane x * = 0.43 (b) and velocity in the plane x * = 2.60 The dashed line in (b) is the contour U x * = 0.5.

Figure 9 :

 9 Figure 9: Drag (a) and fluctuations of forces (b) as a function of a slight yaw angle β: ×, reference SB model; +, O model;, OD-1 model (disturbing square cylinder of 5 cm side) ; ⃝, OD-1 model (disturbing circular cylinder of 2 cm diameter). The disturbing cylinders are placed at x C * = 0.87 and y C * = 0. See table1for the different models configuration and force coefficient measurements.

Figure 10 :

 10 Figure 10: Base pressure distribution for different sideslip yaw angles of the O model: β = -0.7 • (a), β = 0 • (b) and β = 0.8 • (c).

Table 2 :

 2 Mean and standard deviation of the force coefficients for the squareback model extended with a straight cavity (C), and 3 configurations of flow orientation combinations at the trailing edge. The C-TB configuration has top and bottom chamfers (Φ T = 10 • ,Φ B = -10 • ), C-LR has left and right chamfers (Φ L = 10 • ,Φ R = -10 • ) and C-TBLR is combining the two previous configurations with chamfers at all edges.

Figure 11 :

 11 Figure 11: Scheme of the C model geometry with the cavity at the trailing edge; the blue dots locate the pressure taps added for these experiments.

Figure 12 :

 12 Figure 12: Velocity (top) and streamwise vorticity (bottom) in the plane x * = 0.35 downstream of the upper corner of the base: C-TB model (ϕ T = -ϕ B = 10 • and ϕ L= ϕ R = 0 • ), (a); C-LR model (ϕ T = ϕ B = 0 • and ϕ L = -ϕ R = 10 • ), (b); C-TBLR model (ϕ T = -ϕ B = ϕ L = -ϕ R = 10 • , (c)).See table2for the different models configuration and force coefficient measurements.

Figure 13 :

 13 Figure 13: Levels of pressure coefficients around the geometry at x * = -0.22 for the different combinations of flow orientations at the trailing edge. Comparisons between C model and C-TB model (a); C-LR model (b); C-TBLR model (c). See table2for the different models configuration and force coefficient measurements. The arrows locate the pressure taps; the pressure coefficient is always negative and quantified by the distance from the surface (the length for C p = 1 is indicated) and the dotted gray lines locate the contours C p = -0.2.

Table 1 :

 1 Mean and standard deviation of the force coefficients for the squareback model (SB), and 3 configurations of wake control: optimized model (O) with top and bottom chamfers (Φ T = 7.5 • ,Φ B = -5 • ), the optimized model disturbed with vertical cylinder of different diameter sizes placed in the centre of the recirculating region (O-D1 and O-D2). The SB model is used as the reference to compute the drag reductions (DR).

	Configuration	C x	DR	C y	C z	S td(C x ) S td(C y ) S td(C z )
	SB : squareback (ref.)	0,309	-	0,001 -0,047	0.004	0.015	0.006
	O : optimized	0.291 -5.8% -0.003 -0.014	0.003	0.012	0.005
	O-D1 : disturbance 5cm 0.286 -7.4% 0.002 -0.009	0.004	0.007	0.005
	O-D2 : disturbance 2cm 0.291 -5.8% 0.004 -0.014	0.003	0.012	0.005

The rear lift is the part of the lift measured on the rear axle of the geometry.

The cylinder is placed manually so an exploration of a domain (x C * , y C * ) is long and tedious.

Acknowledgments

The authors wish to thank the CNRT R2A, National Centre for Technologic Research for Aerodynamics and Aeroacoustics of ground vehicles for providing the full scale wind tunnel of GIE S2A as well as the financial support for the postdoc position of A. Courbois.