Johan Peralez
email: johan.peralez@gmail.com

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

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

Paolino Tona
email: paolino.tona@ifp.fr

Antonio Sciarretta
email: antonio.sciarretta@ifp.fr

Editor S Di

Cairano J Peralez

ON CONTROL SYSTEMS TECHNOLOGY 1 Organic Rankine Cycle for Vehicles: Control Design and Experimental Results

Keywords: Control design, feedforward, heavy-duty vehicles, observer, Rankine cycle, waste heat recovery (WHR)

The system considered here is an organic Rankine cycle for recovering waste heat from a heavy-duty diesel engine. Because of the highly transient conditions these systems are subject to, control plays a fundamental role to enable the viability and efficiency of those systems. In this context, this paper investigates the problem of control design for superheating (SH) and pressure at evaporator outlet. Based on a moving boundary heat exchanger model, a first controller, which consists of a dynamic feedforward combined to a gain-scheduled PID, is implemented on the pump speed to maintain the SH close to the set-point value. Experimental results illustrate the enhanced performance in terms of disturbance rejection. Then, a second controller based on nonlinear state estimation is proposed. This is a nonlinear feedback law, which allows to adjust the evaporating pressure to time-varying demand with a good accuracy.

I. INTRODUCTION

I N RECENT years, growing awareness of environmental and energy issues has steadily driven research and innovation in the automotive sector. Internal combustion engines (ICEs) are now much more energy efficient, thanks, in particular, to the introduction of new architectures (turbocharging and gasoline direct injection) and to the development of more and more efficient engine control systems. However, despite these significant advances in technology, efficiency of ICE remains limited. Thus, even a modern engine cannot convert more than 50% of fuel energy into useful work. The remaining portion is lost through coolants and exhaust gases [START_REF] Stanton | Systematic development of highly efficient and clean engines to meet future commercial vehicle greenhouse gas regulations[END_REF].

A potential solution to recover this waste heat exists, especially from the heat sources with the highest temperature levels: namely, the tailpipe exhaust gases and the exhaust gas recirculation circuit, when present. Among the waste heat recovery (WHR) technologies, compounding the engine with a Rankine bottoming cycle has been widely considered as the most promising solution [START_REF] Stanton | Systematic development of highly efficient and clean engines to meet future commercial vehicle greenhouse gas regulations[END_REF].

Rankine systems for WHR in vehicles are based on the same principle also used for electricity production in most power plants. A pump circulates a working fluid in a closed loop where an external source supplies heat, via a heat exchanger (or a series of exchangers). Vaporized fluid expands in a turbine or an expander to produce mechanical power. Vapor is then cooled by a condenser, which transfers heat to an external cold sink.

Typically, Rankine cycle-based WHR systems for transport applications can be represented as in Fig. 1. The major differences with stationary applications lie with the strongly transient behavior of the hot source and with the limited cooling potential on board. Both these differences depend on the driving conditions of the vehicle.

To be viable, these systems must be lightweight and compact. Thus, heating, vaporization, and superheating (SH) of the fluid usually take place in a single heat exchanger: the evaporator. In most cases, Rankine systems for automotive applications are designed to produce electricity via a generator connected to the auxiliary network and/or an energy storage system, although produced mechanical power can be used directly via a mechanical connection to the transmission, as in [START_REF] Howell | Development of an ORC system to improve HD truck fuel efficiency[END_REF].

A. Motivation

The WHR systems based on the Rankine thermodynamic cycle have been the focus of intensive research for road vehicles over recent years. For example, BMW [START_REF] Horst | Dynamic heat exchanger model for performance prediction and control system design of automotive waste heat recovery systems[END_REF], [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],

Honda [START_REF] Endo | Study on maximizing exergy in automotive engines[END_REF], and Ford [START_REF] Hussain | Organic Rankine cycle for light duty passenger vehicles[END_REF] worked on this topic for cars, and Cummins [START_REF] Nelson | Exhaust energy recovery[END_REF], Caterpillar [START_REF] Kruiswyk | An engine system approach to exhaust waste heat recovery[END_REF], Daimler Trucks [START_REF] Aneja | Exhaust heat driven Rankine cycle for a heavy duty diesel engine[END_REF], and Volvo [START_REF] Espinosa | Transient organic Rankine cycle modelling for waste heat recovery on a truck[END_REF] for trucks. The interest of manufacturers is justified by announced reductions in fuel consumption ranging from 5 to 10% depending on the system and the driving cycle [START_REF] Nelson | Exhaust energy recovery[END_REF].

As reported in [START_REF] Sprouse | Review of organic Rankine cycles for internal combustion engine exhaust waste heat recovery[END_REF], which presents a state of the art on organic Rankine cycle (ORC) for ICE, most of the research papers are dedicated to design issues (structure, components, and working fluid) and to potential assessment (good examples are [START_REF] Teng | A Rankine cycle system for recovering waste heat from HD diesel engines-WHR system development[END_REF]- [START_REF] Grelet | Transient performance evaluation of waste heat recovery Rankine cycle based system for heavy duty trucks[END_REF]). In addition, although control plays a fundamental role to enable viability and energy efficiency of Rankine cycle WHR systems for automotive applications, papers on control issues are surprisingly a few in numbers.

Despite the development of several Rankine system prototypes for road transport, with a few of them even making their way to demo vehicles, it is still unclear how viable and efficient these systems can be in the real-world applications. In practice, these systems often struggle to attain satisfactory performance over a broad range of operating (transient) conditions. Eventually, if the heat recovery system has to be shut down too often and/or it takes too long to attain power production conditions, expected gains will vanish.

In this context, controlling the working fluid state (namely, SH and pressure) at the evaporator outlet is a key issue when controlling an ORC system, as cycle efficiency and system safety depend on it. This is the focus of this paper.

B. Related Work

In the (scarce) publications for such applications, including experimental results, the authors systematically put forward the difficulty of controlling the working fluid state (SH and pressure) at the evaporator outlet. In addition, as shown in experimental studies, such as [START_REF] Horst | Dynamic heat exchanger model for performance prediction and control system design of automotive waste heat recovery systems[END_REF], [START_REF] Endo | Study on maximizing exergy in automotive engines[END_REF], and [START_REF] Park | A Rankine cycle system for recovering waste heat from HD diesel engines-Experimental results[END_REF], the control of Rankine systems for mobile applications is far from trivial.

Endo et al. [START_REF] Endo | Study on maximizing exergy in automotive engines[END_REF] present a Rankine system for WHR from a spark ignition engine using water as working fluid and a volumetric expander to produce electric power. Experimental results are provided for a decentralized control system composed of two proportional controllers with feedforward.

More recently, the study [START_REF] Horst | Dynamic heat exchanger model for performance prediction and control system design of automotive waste heat recovery systems[END_REF] also presents experimental results based on a steam Rankine cycle. Here, the main originality is the use of a first-principle model for computation of a (static) feedforward on the pump speed.

Global energy management approaches for vehicles' integrating Rankine systems have also been presented, namely, in [START_REF] Merz | On the optimal thermal management of hybrid-electric vehicles with heat recovery systems[END_REF] (for a passenger car) and [START_REF] Willems | Integrated energy & emission management for heavy-duty diesel engines with waste heat recovery system[END_REF] (for a heavyduty vehicle), based on the Pontryagin minimum principle. More recently, [START_REF] Peralez | Optimal control for an organic Rankine cycle on board a diesel-electric railcar[END_REF] addresses the problem of maximizing the power produced by an ORC on board a diesel-electric rail car, first offline via dynamic programming, then online via dynamic real-time optimization. But the experimental assessment of such (high-level) approaches still requires an efficient low-level control, allowing to tightly control the working fluid state at the evaporator outlet.

On the more general topic of ORC for WHR operating with variable heat sources (not necessarily for transport applications), Hou et al. [START_REF] Hou | Supervisory predictive control of evaporator in organic Rankine cycle (ORC) system for waste heat recovery[END_REF] and Zhang et al. [START_REF] Zhang | Dynamic modeling and multivariable control of organic Rankine cycles in waste heat utilizing processes[END_REF], [START_REF] Zhang | Constrained predictive control based on state space model of organic Rankine cycle system for waste heat recovery[END_REF] apply control strategies based on linear models (Linear Quadratic Regulator or Model Predictive Control control), validated around one operating point. Reference [START_REF] Quoilin | Dynamic modeling and optimal control strategy of waste heat recovery organic Rankine cycles[END_REF] deals with the modeling and control of an ORC on larger operational ranges. However, the hot source variations used for simulation are much slower than those observed at the exhaust of an automotive engine, especially in terms of mass flow rate.

A somewhat richer literature exists on dynamic modeling and control of vapor-compression cycles, the reverse of Rankine cycle [START_REF] He | Modeling of vapor compression cycles for multivariable feedback control of HVAC systems[END_REF]- [START_REF] Alleyne | Advances in energy systems modeling and control[END_REF]. More particularly, in [START_REF] Rasmussen | Gain scheduled control of an air conditioning system using the Youla parameterization[END_REF], a solid approach is presented for multivariable H ∞ synthesis.

C. Paper Contribution and Organization

This paper addresses the control design of an ORC with a kinetic turbine. For this system, a hierarchical and modular control structure has been designed, implemented, and validated experimentally on an engine test-bed cell on a heavyduty road cycle. First, a baseline strategy for SH regulation (with the pump mass flow rate as the manipulated variable) is improved by combining a dynamic feedforward term to a commonly used gain-scheduled PID approach. The robustness of the feedforward term computation (based on model inversion) is demonstrated, allowing an embedded implementation. Compared with the preliminary version of this paper presented in [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 [START_REF] Peralez | Control design for an automotive turbine Rankine cycle system based on nonlinear state estimation[END_REF], experimental results, representative of a long-haul truck mission, illustrate here the enhanced performance in terms of disturbance rejection.

This approach is then extended to the multivariable case, with the use of the evaporator bypass as an additional actuator, allowing the tracking of a pressure set point while maintaining the previous SH regulation performances. The proposed scheme combines an additional nonlinear controller with an implicit extended Kalman filter (EKF) for wall temperature estimation. An original contribution lies with the use of the exhaust bypass to recover the degree of freedom lost, because the turbine is kinetic and not of the positive-displacement type.

This paper is organized as follows. Section II describes the main components of the system (actuators and available sensors) and details the main control requirements. Section III briefly describes the reference model from which is derived a nonlinear reduced model for control purpose. The SH controller is then designed in Section IV and validated with experimental results. In Section V, a new multivariable control scheme based on an observer is introduced, and the results are validated in simulation. Finally, the conclusions are drawn in Section VI.

II. SYSTEM DESCRIPTION

The system considered here is shown in Fig. 2. It is an ORC system for a heavy-duty diesel engine using a turbine for the expansion of the working fluid. A description of the main system variables follows, using the nomenclature given in Table I.

A. Inputs-Outputs

Four actuators are available, as shown in Fig. 3. 1) The pump, whose speed N pump allows to control the working fluid mass flow entering the evaporator. It will be used to control the SH (defined as the difference between the fluid temperature and its evaporation temperature) at evaporator outlet.

2) The evaporator bypass Vo evap controlling the fraction of exhaust gas entering the evaporator. It will be used to control the pressure in the evaporator. The choice of such a decentralized control for SH and pressure control is based on two-time-scale considerations (working fluid dynamics are much faster than wall thermal dynamics). Hence, as detailed in Section V, SH will be assumed perfectly regulated by pump mass flow, allowing to encompass the model complexity.

3) The turbine speed N turb entering the load control system.

In the present case of a kinetic turbine, its influence on working fluid condition at evaporator outlet can be neglected while it provides a degree of freedom to optimize the power production. This optimization is not within the scope of this paper (the interested reader can refer to [START_REF] Bao | A review of working fluid and expander selections for organic Rankine cycle[END_REF] and references therein). 4) The turbine bypass Vo turb intended to protect the turbine from fluid condensation. It will be kept open as long as fluid condition at evaporator outlet is not safe for the turbine while it should be kept closed during production mode. Hence, it will be used as an ON-OFF actuator. As shown in Fig. 3, measurements of the pressuretemperature pairs ( p, T ) between each component are available (corresponding to the corner points 1, . . . , 4 of the thermodynamic cycle in Fig. 4). The available measurements can be used to estimate key output variables, such as SH, subcooling, or enthalpy, from fluid thermodynamic properties.

Respecting exhaust gas conditions, the temperature T exh is measured while an estimation of the mass flow ṁexh is provided by the engine control unit.

By contrast, no sensors are available inside the heat exchangers. In particular, the temperature of the wall is not measured.

B. Control Requirements

The main objective of the supervision and control system is to maximize the production of electric energy during vehicle usage, in the presence of variations in engine exhaust gas conditions.

This objective can be addressed by a two-level closed-loop control strategy [START_REF] Peralez | Optimal control for an organic Rankine cycle on board a diesel-electric railcar[END_REF], [START_REF] Quoilin | Dynamic modeling and optimal control strategy of waste heat recovery organic Rankine cycles[END_REF], where a real-time optimization provides (optimal) set points to the low-level controller, which computes the control values. This paper focuses on the evaporator outlet pressure and SH low-level controllers required by the supervisory control proposed in [START_REF] Peralez | Optimal control for an organic Rankine cycle on board a diesel-electric railcar[END_REF] for a diesel-electric railcar, which could also be applied to the present case study.

Even without the optimization layer (in which case the set point can be chosen constant), the proposed low-level controller remains relevant as it ensures system safety. Indeed, several output constraints must be satisfied, in particular, at the evaporator outlet, to protect the turbine from fluid condensation (requiring the SH to be strictly positive), and to avoid an unsafe pressure level.

III. MODELING

This section describes the modeling of the ORC components needed to predict the fluid state at evaporator outlet (characterized by the pressure p and SH values) during nominal operation, that is, when superheated vapor feeds the turbine. A reduced model is then presented, for control design.

A. Evaporator Model

In systems, implementing the Rankine thermodynamic cycle, the working fluid enters the evaporator in a liquid state and exits in a superheated vapor state. A classic 1-D representation of heat exchangers (along the fluid displacement direction z) is given by mass and energy balances for the fluid and the wall [START_REF] Wei | Dynamic modeling and simulation of an organic Rankine cycle (ORC) system for waste heat recovery[END_REF]. The evolution of fluid and wall variablesnamely fluid mass flow ṁ, specific enthalpy h, pressure p, and wall temperature T w -is then described by

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ V f L ∂ρ ∂t + ∂ ṁ ∂z = 0 V f L ∂(ρ h) ∂t - V f L ∂ p ∂t + L ∂ ṁ ∂z = S f L α f (T w -T f ) c w ρ w V w ∂ T w ∂t = S f α f (T w -T f ) + S w α exh (T exh -T w ) (1)
where V f , V w , S f , and S w are (constant) design parameters. Density ρ and temperature T f of the fluid are nonlinear maps of p and h. T exh , representing the exhaust gas temperature, is considered as a time-varying input for the evaporator model. α exh (resp. α f ) is the heat exchange coefficient between the wall and the exhaust gas (resp. the working fluid). Notice that the above model is a system of three coupled nonlinear partial differential equations. Such a representation is difficult to deal with for control and estimation purposes. Therefore, a discretization of system ( 1) is used. Among the methods for discretization of system (1), the use of moving boundaries (MBs) results in a relatively low order model. As shown in Fig. 5, the MB model monitors the length of each fluid phase along the evaporator: the normalized zone lengths L 1 , L 2 , and L 3 track the liquid, two phase and vapor zones, respectively. Fluid and wall dynamics are then captured by the seven state variables (T w,1 , T w,2 , T w,3 , L 1 , L 2 , h 3 , and p). A more detailed representation of this model may be found in [START_REF] Wei | Dynamic modeling and simulation of an organic Rankine cycle (ORC) system for waste heat recovery[END_REF]. In the following, it is referred as the reference model that will be used for simulation purposes in Section V-D.

In an ORC (that is with a dry working fluid [START_REF] Kang | Design and experimental study of ORC (organic Rankine cycle) and radial turbine using R245fa working fluid[END_REF]), the fluid enters the condenser in a vapor state and exits in a liquid state or in a two-phase state depending on the cooling conditions and on the pressure imposed by the tank, acting as a separator. Applying the MB approach to the low-pressure (LP) part of the circuit generally yields a high-order hybrid model (the number of states changes depending on fluid conditions at the condenser outlet). However, as it will be justified in Section III-C, condenser modeling is not necessary for our control purposes.

B. Pump and Turbine Models

Since pump and turbine dynamics are very fast compared with exchanger dynamics, they are modeled by algebraic equations. The (positive displacement) pump produces a mass flow rate proportional to its rotational speed

ṁpump (t) = ρ pump (t)η pump V pump N pump (t) (2)
where its volumetric efficiency η pump can be considered constant in nominal conditions. At the evaporator outlet, the fluid in a vapor state expands through the turbine nozzle.

For nominal mass flow rate values, the fluid then reaches supersonic speeds that allows to neglect the influence of the outlet pressure [START_REF] Feru | Experimental validation of a dynamic waste heat recovery system model for control purposes[END_REF] ṁevap,out

(t) = C d S turb 2ρ( p evap,out (t), SH(t)) p evap,out (t) (3)
where S turb is the turbine orifice equivalent area and C d is the discharge coefficient.

C. Nonlinear Control-Oriented Model

Designing a controller based on a full model of the Rankine cycle proves a complex and challenging problem. In the literature, model-based control generally copes with the complexity of a reference model, such as the one described earlier, by linearization around an operating point [START_REF] Luong | Linear quadratic integral control of an organic Rankine cycle for waste heat recovery in heavy-duty diesel powertrain[END_REF] or a set of operating points [START_REF] Feru | Model predictive control of a waste heat recovery system for automotive diesel engines[END_REF]. However, [START_REF] Horst | Dynamic heat exchanger model for performance prediction and control system design of automotive waste heat recovery systems[END_REF] and [START_REF] Peralez | Towards model-based control of a steam Rankine process for engine waste heat recovery[END_REF] underline the need of considering the nonlinear behavior of the Rankine system for control design, as its static gains and response times strongly vary with the operating conditions. Reference [START_REF] Rasmussen | Dynamic modeling and advanced control of air conditioning and refrigeration systems[END_REF] attempts to cast the evaporator control problem for a similar system (a vapor-compression cycle for refrigeration) in a linear parameter varying framework; this approach is made very difficult by the large number of scheduling variables involved.

Here, a two-time-scale dynamic behavior is assumed where wall temperatures capture the slow dynamics. Moreover, only the high-pressure (HP) parts of the ORC will be used (which includes the pump, the evaporator, and the turbine).

Indeed, expression (3) outlines that the effect of fluid conditions at the turbine outlet can be neglected. Furthermore, the thermal inertia of the tank slows down the variations of pump inlet temperature. Such considerations allow to decouple the HP part of the cycle from the LP part (which includes the condenser and the tank). Evaporator inlet working fluid temperature T f,in is then considered as a (measured) disturbance for the HP part.

The model can be further reduced assuming the fluid to be at thermodynamic equilibrium. For a steam Rankine cycle, [START_REF] Peralez | Towards model-based control of a steam Rankine process for engine waste heat recovery[END_REF] shows that slow evaporator dynamics are due to wall temperatures and that a third-order model can be effectively used for disturbance rejection.

The reduced model of the evaporator is derived from the physical equations given in the following. To simplify notations, all the subscripts of variables relating to the evaporator are dropped.

1) Working Fluid Mass Balance: Assuming working fluid at a steady state leads to consider a homogeneous mass flow (denoted by ṁ in the following) along the evaporator equal to the pump mass flow

ṁin,i (t) = ṁout,i (t) = ṁpump (t), i = 1, . . . , 3.
2) Working Fluid Energy Balance: Energy balance equilibrium for working fluid is written for each of the three zones (indexed by i )

0 = ṁ(t) (h in,i (t) -h out,i (t)) + Q f,i (t) L i (t) (4)
where

Q f,i (t) = S f α i (T w,i (t) -T f,i (t))
represents the heat transfer from wall and S f is the wall-fluid exchange area. T f,2 (the two-phase zone temperature) is the evaporation temperature that only depends on pressure. The liquid zone temperature can be computed as the average of evaporator inlet temperature and boiling temperature

T f,1 ( p(t), t) = 0.5 T f,in (t) + 0.5 T f,2 ( p(t)).
The vapor zone temperature is expressed according to the SH value

T f,3 ( p(t), SH(t)) = T f,2 ( p(t)) + 0.5 SH(t).
Notice that a (constant) coefficient α i is chosen for each zone, reflecting significant differences in heat transfer efficiency induced by the working fluid conditions.

3) Wall Energy Balance: Wall energy balance yields the following (dynamic) equation for each zone:

m w c w dT w,i dt = Qexh,i (t) -Q f,i (t) (5)
where

Qexh,i (t) = ṁexh (t) c exh 1 -exp - α exh S exh ṁexh (t) c exh × [T exh (t) -T w,i (t)]
represents the heat transfer from exhaust gas and S exh is the wall-exhaust gas exchange area.

4) Differential Algebraic Equation System:

The abovementioned equations are completed by interface equations. Fluid enthalpies and mass flows at the inlet of zones 2 and 3 correspond to those at the outlet of zones 1 and 2, that is

h in,i+1 (t) = h in,i (t), i = 1, 2.
For the two-phase zone, inlet and outlet enthalpies correspond to saturation values

h in,2 ( p(t)) = h l ( p(t)), h out,2 ( p(t)) = h v ( p(t))
where h l and h v are, respectively, the saturated liquid enthalpy and the saturated vapor enthalpy that, therefore, only depend on pressure. Eventually, balance and interface equations form a differential algebraic equation (DAE) system with three dynamic states (T w,1 , T w,2 , and T w,3 ).

Remark 1: Notice that the above model is only relevant for nominal conditions as it does not handle the problem when one or two zones do not exist (during start-stop operations). However, the model-based controller detailed in Section IV will be completed to manage start-stop procedures. Details on this overall strategy will be given in Section IV-B.

In the case of a positive-displacement expander, the additional degree of freedom yields an explicit system, and a model inversion can be computed analytically, allowing effective tracking of pressure at evaporator outlet [START_REF] Peralez | Towards model-based control of a steam Rankine process for engine waste heat recovery[END_REF]. In Section IV, this inversion-based approach will be adapted to the case of a kinetic turbine and will be used to analyze the robustness of its computation scheme.

IV. SUPERHEATING CONTROL

In this section, the objective is to design a controller of the SH at the evaporator outlet. Effective SH control is a key issue when controlling an ORC system, as cycle efficiency and turbine safety depend on it. As pointed out in [START_REF] Yamamoto | Design and testing of the organic Rankine cycle[END_REF], SH must be kept as low as possible to ensure good ORC efficiency when using high molecular weight organic fluids. But, SH must always remain positive to prevent the formation of droplets that could damage the turbine. In a nominal mode, when power is produced, the expander bypass is closed (for the sake of efficiency) and only one actuator is fast enough to tightly control SH at evaporator outlet: the pump speed N pump or equivalently, from relation (2), the pump mass flow rate ṁpump .

The ORC system considered has a second (slower) actuator, the exhaust bypass, which can be used to reduce the hot gas enthalpy flow rate so as to limit the pressure in the circuit. This slower control loop, detailed in Section V, acts as a disturbance to the faster SH control loop.

A. Nonlinear Inversion-Based Control

The SH control scheme is described in Fig. 6. The focus is on improving the baseline control strategy using a model-based approach. The improvements come from an extensive system identification campaign, which allows model-based tuning of PID controllers and, more particularly, from a dynamic feedforward term computed from a nonlinear reduced model of the HP part of the system.

1) Feedback Controller: Input-output dynamics are identified at different operating points via a set of linear models. The good fitting obtained with damped second-order transfer functions (see Fig. 7) justifies the use of a gain-scheduled PID controller (and not just a PI) in the feedback path. PID controllers are tuned with the well-known internal model control (IMC) method (more details on IMC method can be found in [START_REF] Ang | PID control system analysis, design, and technology[END_REF]). Following the gain scheduling idea presented in [START_REF] Rasmussen | Gain scheduled control of an air conditioning system using the Youla parameterization[END_REF] in the context of an air conditioning system, the feedback control value is finally obtained by a linear interpolation of these local controllers, using the pressure value as a scheduling variable.

2) Inverse Reduced Model: The evaporator model derived in Section III-C is merged with the nozzle mass flow equation ( 3) and inverted to be used in the feedforward path for disturbance rejection. Evaporator mass flow computed from this reduced model then constitutes the feedforward part u d of the control input (see Fig. 6), as detailed thereafter. 1Notice that, among the disturbances appearing in the model and which can be gathered into the vector d = [T exh ṁexh T f,in ] T , T exh and T f,in are measured, while an estimation of ṁexh is provided by the engine control unit.

Providing that pressure p is known, the dynamic part of the resulting system is described by the following (explicit) system:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ dT w,1 dt = -β 1 (T w,1 -T f,1 ( p, d)) -η exh (T w,1 -T exh ) dT w,2 dt = -β 2 (T w,2 -T f,2 ( p)) -η exh (T w,2 -T exh ) dT w,3 dt = -β 3 (T w,3 -T f,3 ( p, SH)) -η exh (T w,3 -T exh ) (6)
where

β i = S f α i m w c w , i = 1, . . . , 3 η exh = ṁexh c exh m w c w 1 -exp - α exh S exh ṁexh c exh while p can be computed from the (implicit) algebraic part ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ṁ = C d S turb √ 2ρ( p, SH) p L 1 = ṁ h l ( p) -h in (d)) S f α 1 (T w,1 -T f,1 ( p, t)) L 2 = ṁ h v ( p) -h l ( p) S f α 2 (T w,2 -T f,2 ( p)) L 3 = ṁ h out ( p, SH) -h v ( p) S f α 3 (T w,3 -T f,3 ( p, SH)) L 1 + L 2 + L 3 = 1. (7)
In order to ensure that a robust implementation of the controller is possible, the implicit system (7) needs to be analyzed. This is done in Proposition 1 that applies for the physical domain defined in Assumption 1. For the sake of simplicity, SH SP will be assumed constant and equal to 30 K in the following, but the results remain valid for any other strictly positive values of the SH.

Assumption 1: The following physical constraints are fulfilled:

p < p max (8a) T min f,in < T f,in < T max f,in (8b) T f,i < T w,i , i = 1, 2, 3 ( 8 c )
where p max , T min f,in , and T max f,in are positive constants, characterizing the nominal conditions for the working fluid. 2Remark 2: The physical constraints (8a) and (8b) imposed by Assumption 1 are satisfied for the considered system and working fluid in normal conditions (larger pressures or inlet temperatures do not match safety conditions and/or efficiency considerations and would lead to system shutdown by the control system supervisor).

On the other hand, the physical constraint (8c) means that the fluid is heated by the wall. This might not be verified during start-stop operations. However, under nominal conditions (for which the control law is designed), none of these conditions is restrictive.

Proposition 1: Consider the algebraic system (7) applied to the working fluid R245fa, and assume that Assumption 1 holds for the following operation conditions:

p max = 25 bar, T min f,in = 10°C, T max f,in = 40°C.
Then, the residual function ϕ defined by

ϕ = 1 -L 1 -L 2 -L 3
is strictly decreasing with respect to p, i.e., (∂ϕ/∂ p) < 0.

The proof of Proposition 1 is postponed to the Appendix. Finally, since ϕ decreases monotonically with p, one can numerically solve ( 7), e.g., by a bisection method, ensuring a robust (inline) implementation of the controller. In other words, given a disturbance vector d and a state vector T w computed by the integration of (stable) equation system ( 6), one can compute the corresponding pressure value p. The feedforward part u d is then deduced from the nozzle mass flow equation

u d = C d S turb 2 ρ( p, SH) p.
Eventually, the pump mass flow u d is transposed into a pump speed control signal assuming the linear relation [START_REF] Howell | Development of an ORC system to improve HD truck fuel efficiency[END_REF].

B. Experimental Results

In this section, experimental results demonstrate that the dynamic feedforward allows significant gains in performance when coupled to the commonly used PID approach for SH regulation. Experimental validation was performed by means of the Rankine system plugged on a heavy-duty diesel engine shown in Fig. 8.

1) SH Regulation Performance: In Fig. 9, we compared the two control strategies (with or without feedforward), for the same variations of exhaust gas conditions. These transient profiles [Fig. 9(a)] have been obtained by varying engine torque and speed in order to compare the two strategies on a large panel of operating conditions.

The feedforward action demonstrates significant performance improvements [Fig. 9(c)]. Indeed, the maximum error (SH-SH SP |) is only about 1.9 K in this case compared with about 10 K without feedforward. Notice that this would allow to reduce considerably the SH set point, thus increasing cycle efficiency, while keeping the same safety margins during power production [START_REF] Yamamoto | Design and testing of the organic Rankine cycle[END_REF].

2) SH Regulation on a Representative Heavy-Duty Road Cycle:

The SH controller detailed in Section IV-A (designed for nominal conditions) must be complemented to manage start-stop procedures. The considered supervision structure for the Rankine pilot process is similar to that described in [START_REF] Tona | Supervision and control prototyping for an engine exhaust gas heat recovery system based on a steam Rankine cycle[END_REF]. For the sake of completeness, safety and management modes used during the experiment are detailed in the following.

1) Stop: This mode is maintained as long as the exhaust gas is not hot enough. Evaporator bypasses the evaporator and pump speed is brought to zero. 2) Cold Idle: Exhaust gas enters the evaporator while pump circulates fluid at a constant mass flow (open loop). 3) Hot Idle: This mode is enabled when superheated vapor appears at the evaporator outlet (SH>0). The SH control scheme described in Fig. 6 is then used (closed loop). 4) Nominal Mode: This mode is enabled when SH is steady and large enough to ensure turbine safety. Power production is then possible.

control strategy is finally validated on a realistic exhaust conditions profile, corresponding to a heavy-duty road cycle (Fig. 10). SH is tightly controlled along the entire duration of the cycle (lasting over 50 min), demonstrating the robustness of the proposed controller. Moreover, thanks to an early stabilization of the SH, the nominal mode is enabled during most of the experiment, thus creating the conditions for continuous power production [Fig. 10(c)].

Notice that to ensure a safe pressure level in the circuit, the evaporator bypass was used to limit the hot gas enthalpy flow rate [Fig. 10(a)].

This heuristic use of the bypass, which does not consider the system dynamics, presents two main drawbacks. First, the hot gas flow may be restricted abruptly in response to increases in engine load. Second, it does not allow the tracking of a varying pressure set point, which may be required by an overall Rankine system supervisor (for instance, in the presence of limited cooling capacities, see [START_REF] Peralez | Optimal control for an organic Rankine cycle on board a diesel-electric railcar[END_REF]).

This trajectory-tracking problem is addressed in Section V.

V. PRESSURE CONTROL

A. Problem Formulation

In Section IV, it was experimentally shown that a nonlinear model-based monovariable control approach is more efficient than the commonly used PID approach for SH regulation. In this section, this approach is extended to the multivariable case, with the use of the evaporator bypass as a second actuator, which also allows the tracking of a pressure set point while ensuring a safe level for SH. Interestingly, despite the importance of this bypass, which is mandatory for Fig. 10. Experimental result. Validation on a realistic profile of exhaust gas conditions. (a) Measured disturbances for the SH control. Working fluid temperature T f,in at evaporator inlet and exhaust gas conditions (temperature T exh and mass flow ṁexh,evap entering the evaporator). (b) SH is tightly controlled (using the pump speed N pump ) while a safe pressure level is ensured. (c) Automatic startup system achieves nominal operating conditions in about 1 min. safety reasons, only heuristic control approaches have been reported [START_REF] Horst | Dynamic heat exchanger model for performance prediction and control system design of automotive waste heat recovery systems[END_REF], [START_REF] Feru | Experimental validation of a dynamic waste heat recovery system model for control purposes[END_REF]. Now, let us consider that the output of the system is y = [y 1 y 2 ] T = [SH p] T (corresponding to the working fluid state at evaporator outlet) and the dynamic state is T w = [T w,1 T w,2 T w,3 ] T in accordance with the reduced model derived in Section III-C. In what follows, the fraction of exhaust gas entering the evaporator is assumed to be proportional to the valve position Vo evap . Consequently, the chosen control vector is u = [u 1 u 2 ] T = [ ṁ η exh ] T where the new second manipulated variable is Rather than choosing the valve position as a control variable, the selection of η exh allows to obtain a system affine in the control as it is shown in the following. This two-input two-output control problem is treated in two steps, assuming a two-time-scale dynamic behavior. First, the control design proposed in Section IV-A is used to regulate the SH y 1 with pump mass flow u 1 . Indeed, experimental results in Section IV-B have shown that a control based on a gainscheduled PID with a dynamic feedforward allows a tight SH regulation even in the presence of fast transient conditions on exhaust gas. Then, a second nonlinear controller for the pressure y 2 will be designed that assumes that the (fast) SH dynamics are perfectly regulated.

η exh = Vo evap ṁexh c exh m w c w 1 -exp - α exh S exh Vo evap ṁexh c exh .
As it can be noticed, the proposed control requires the knowledge of all system states (i.e., T w i measurements). To overcome this problem, in Section IV, an open-loop model was used for the estimation of needed states. In the following, instead of the open-loop model, an observer is proposed for the state vector estimation, which will improve the robustness of the closed loop against model uncertainties. The control scheme described in Fig. 6 is taken up and completed, as shown in Fig. 11.

Assuming that the pump mass flow is fixed by the value of SH, system ( 6) and ( 7) can be rewritten the following implicit form:

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ dT w dt = f (T w , p, d) + u 2 g(T w , d) ϕ(T w , p, d) = 0 y = p (9)
where

f (T w , p, d) = ⎛ ⎝ β 1 (T f,1 ( p, d) -T w,1 ) β 2 (T f,2 ( p) -T w,2 ) β 3 (T f,3 ( p) -T w,3 ) ⎞ ⎠ g(T w , t) = ⎛ ⎝ T exh -T w,1 T exh -T w,2 T exh -T w,3 ⎞ ⎠ (10)
with

β i = S f α i m w c w , i = 1, . . . , 3 and
ϕ(T w , p, d) = 1 -(L 1 + L 2 + L 3 )
with

L 1 = u 1 h l ( p) -h in (d) S f α 1 (T w,1 -T f,1 ( p, d)) L 2 = u 1 h v ( p) -h l ( p) S f α 2 (T w,2 -T f,2 ( p)) L 3 = u 1 h out ( p) -h v ( p) S f α 3 (T w,3 -T f,3 ( p))
.

In what follows, the main focus will be on the design of the above-mentioned observer and pressure controller, based on system (9).

B. Control Design

From Proposition 1, (∂ϕ/∂ p) < 0 for all (T w , p, d) in the physical domain defined by Assumption 1. Hence, (∂ϕ/∂ p) -1 is well defined on this domain, and the implicit model given by system (9) can be transformed [START_REF] Hammouri | An observer design for a class of implicit systems[END_REF] into an ordinary differential equation system

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ dT w dt = f (T w , p, d) + u 2 g(T w , d) d p dt = - ∂ϕ ∂ p (T w , p, d) -1 ∂ϕ ∂ T w (T w , p, d) × [ f (T w , p, d) + u 2 g(T w , d)] y = p (11)
where

∂ϕ ∂ T w (T w , p, d) = S turb C d ρ( p) p × ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ h l ( p) -h in (d) S f α 1 (T w,1 -T f,1 ( p, d)) 2 h v ( p) -h l ( p) S f α 2 (T w,2 -T f,2 ( p)) 2 h out ( p) -h v ( p) S f α 3 (T w,3 -T f,3 ( p)) 2 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠
. [START_REF] Sprouse | Review of organic Rankine cycles for internal combustion engine exhaust waste heat recovery[END_REF] Notice that the above system is affine in the control u 2 . The feedback controller introduced in Proposition 2 takes advantage of this structure.

Proposition 2: Let k be a positive real number, then the control3

u 2 = -ϕ T T w f (T w , p, d) + k( p SP -p) ϕ T T w g(T w , d) (13)
asymptotically stabilizes p to the set point p SP for all (T w , p, d) in the physical domain defined by Assumption 1. Moreover, the input u 2 and the states of system (9) remain bounded.

In order to prove Proposition 2, Lemma 2 is needed. Lemma 1: For every T f,i (i = 1, . . . , 3), T exh bounded and for the input u 2 , which takes values in some bounded set U ⊂ R + , the solution T w of system

dT w dt = f (T w , p, d) + u 2 g(T w , d) ∀ T wi (t 0 ) ∈ R + (14) IEEE TRANSACTIONS ON CONTROL SYSTEMS TECHNOLOGY
is uniformly bounded with respect to t in R + , i.e., system [START_REF] Nelson | Exhaust energy recovery[END_REF] is bounded-input bounded-output (BIBO) stable.

In the following, we first start with the proof of Lemma 1; then the proof of Proposition 2 will be given.

Proof (Lemma 1): Considering the positive definite quadratic function V (T wi ) = T 2 wi /2, the derivative of V along the trajectories of the dynamical system ( 14) is given by

dV dt = T w,i (β i (T f,i -T w,i ) + u 2 (T exh -T w,i )) = -T 2 w,i (β i + u 2 ) + T w,i (β i T f,i + u 2 T exh ) ≤ -2aV + b √ 2V ≤ - √ V (2a √ V -b √ 2) ( 15
)
where a = inf(β i + u 2 ) and b = sup(β i T f,i + u 2 T exh ). Notice that a and b are strictly positive and bounded. From [START_REF] Grelet | Transient performance evaluation of waste heat recovery Rankine cycle based system for heavy duty trucks[END_REF], it can be deduced that (dV /dt) is strictly negative for T w , such that {|T w,i | > b/a, i = 1, . . . , 3}. Hence, solutions starting from the set {V (T w,i ) ≤ (b/a) 2 /2} at an initial time t 0 remain there for t > t 0 , and system ( 11) is BIBO stable. Proof (Proposition 2): Let u 2 be defined as in [START_REF] Teng | A Rankine cycle system for recovering waste heat from HD diesel engines-WHR system development[END_REF]. Since Assumption 1 holds, u 2 is bounded. Then, the second equation of system [START_REF] Nelson | Exhaust energy recovery[END_REF] becomes

d p dt = -kϕ -1 p ( p SP -p). (16)
Moreover, from Proposition 1, it is deduced that ϕ p < 0. So, obviously, pressure dynamics in ( 16) are asymptotically attracted by p SP . Notice that the proof of the asymptotic stability assumes that the model is perfect. In practice, parameter k in the control design [START_REF] Teng | A Rankine cycle system for recovering waste heat from HD diesel engines-WHR system development[END_REF] will be chosen sufficiently large to deal with model uncertainties.

C. Observer Design

The problem of wall temperature estimation was addressed by Luong [START_REF] Luong | Modeling, estimation, and control of waste heat recovery systems[END_REF], where he proposed to apply an EKF on the MB reference model described in Section III-A. Unfortunately, the limits of computational resources, particularly, stringent in an embedded context, prevent the use of such a complex model. Therefore, in what follows, an observer based on the reduced model is designed.

Let us consider the class of implicit system of the following form, which includes system (9):

⎧ ⎪ ⎨ ⎪ ⎩ ẋ = F(x, z, u) (x, z) = 0 y = h(x, z) (17)
where

y ∈ R p , u ∈ U ⊂ R m , (x, z) ∈ R n × R d , F
, h, and = ( 1 , . . . , d ) T are assumed to be sufficiently smooth with respect to their arguments, and

∂ ∂z (x, z) is of full rank ∀(x, z) ∈ M ( 18
)
where M is the set of zeros of

M = {(x, z) ∈ R n × R d , s.t. (x, z) = 0}. ( 19
)
For the implicit system of index 1 [START_REF] Merz | On the optimal thermal management of hybrid-electric vehicles with heat recovery systems[END_REF], Nikoukhah [START_REF] Nikoukhah | A new methodology for observer design and implementation[END_REF] has designed an implicit observer described by a dynamic system coupled with an optimization problem permitting to solve the algebraic constraint at each time instant.

Definition 1: An implicit observer for system [START_REF] Merz | On the optimal thermal management of hybrid-electric vehicles with heat recovery systems[END_REF], which converges for every input u ∈ U, is a dynamical system of the form

⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ d x dt = F( x, ẑ, l, y, u) dl dt = L(l, u), l ∈ O is an open subset of R k ( x, ẑ) = 0 (20)
such that for every u ∈ U, the following holds. 1) ∀l(0) ∈ O, the trajectory (l(t)) t ≥0 lies into a compact subset of O.

2) ∀r > 0; ∀(x(0), z(0)) ∈ M; ∃r > 0; ∀( x(0), ẑ(0))) ∈ M, ( x(0) -x(0), ẑ(0) -z(0)) < r implies ( x(t) - x(t), ẑ(t) -z(t)) < r , for all t ≥ 0. 3) ∀(x(0), z(0)) ∈ M; ∀( x(0), ẑ(0))) ∈ M, lim t →+∞ ( x(t) -x(t), ẑ(t) -z(t)) = 0.
Here, conditions 1) and 2) ensure the stability of the observer, and condition 3) guarantees its convergence.

Based on this definition and using the notations z = p and x = T w with (T w , p) ∈ R 3 × R and y = p, the gain l(t) is given by a Riccati equation, and an implicit EKF for system (9) is given by [START_REF] Nikoukhah | A new methodology for observer design and implementation[END_REF]

-[45] ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ d Tw dt = f ( Tw , p) + u 2 g( Tw , p) -SC T R -1 ( p -p) ϕ( Tw , p) = 0 d S dt = AS + S A T -SC T R -1 C T S + Q S(0) = S(0) T > 0 ( 21
)
where Q is a constant symmetric positive definite matrix, R is a real positive constant, and

A = ∂( f + u 2 g) ∂ T w ( Tw , p,d) , C = ϕ -1 p ϕ T T w .

D. Simulation Results

To demonstrate the effectiveness of the control structure described in Fig. 11, this section presents a set of simulation results 4 with realistic operating conditions. After the observer validation, the controller is tested on the seven-state DAE reference model described in Sections III-A and III-B, showing remarkable robustness in the presence of model uncertainties.

1) Observer Validation: Observer (21) performance in the presence of measurement noise is shown in Fig. 12. System simulation is conducted under transient conditions. Observer ( 21) is initialized with initial conditions different from the reference model. The states of the observer are found to converge faster to their target values than the states of the open-loop model, if they both start with the same initial conditions. Moreover, the observer provides an efficient Simulation results under varying inlets. Reference model (solid lines), observer (dashed lines), and reference model with observer initial conditions (dotted lines). filtering of pressure measurement (Fig. 13) that will be useful to the control feedback.

2) Closed-Loop Control Assessment: Experimental data of exhaust gas conditions collected from the heavy-duty road cycle [Fig. 10(a)] were used to evaluate the closed-loop performances under realistic conditions. Moreover, in order to assess the robustness of the control, errors were introduced on two main model parameters: an error of 10% on the heat exchange coefficient between exhaust gas and evaporator wall, and of 20% on the wall volume. equivalently, power production demand-cannot be satisfied due to insufficient exhaust gas heat flow.

Notice that, in this simulation, an arbitrary pressure setpoint trajectory was used to illustrate the tracking capability of the low-level control, without attempting any direct WHR optimization. However, the proposed scheme opens up the possibility to implement a two-level closed-loop control strategy [START_REF] Peralez | Optimal control for an organic Rankine cycle on board a diesel-electric railcar[END_REF], [START_REF] Quoilin | Dynamic modeling and optimal control strategy of waste heat recovery organic Rankine cycles[END_REF], where an energy management system would provide optimal set points to the above-described low-level controllers. Indeed, in the context of transport applications, due to limited cooling capacities [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], [START_REF] Espinosa | Transient organic Rankine cycle modelling for waste heat recovery on a truck[END_REF], it cannot be assumed that the optimal strategy is to minimize the exhaust bypass.

Finally, to emphasize the benefits of the observer, a second simulation is conducted where state estimation is performed via an open-loop model. Fig. 17 compares the resulting SH regulation performance: the maximum error is more than 8 K in this case, while it is about 5 K when the observer is used. Thus, observer was found to improve the robustness of the model-based control strategy. This would allow to reduce the SH set point and thus, once again, increase cycle efficiency while keeping the same safety margins.

VI. CONCLUSION

Effective SH control is a key issue when controlling an ORC system for engine WHR. In this paper, the control design for an ORC with a kinetic turbine has been presented and assessed.

The baseline strategy for SH regulation was improved by combining a dynamic feedforward term to a commonly used gain-scheduled PID approach. The robustness of the feedforward term computation (based on model inversion) has been demonstrated, allowing a real-time implementation. Experimental results, representative of a long-haul truck mission, illustrate the enhanced performance in terms of disturbance rejection.

This approach has been extended to the multivariable case, with the use of the evaporator bypass as an additional actuator, providing the possibility to track a second variable (a pressure set point) while maintaining the previous SH regulation performance. The proposed scheme combines an additional nonlinear controller to an implicit EKF for wall temperature estimation. This opens up the possibility to implement a two-level closed-loop control strategy [START_REF] Peralez | Optimal control for an organic Rankine cycle on board a diesel-electric railcar[END_REF], [START_REF] Quoilin | Dynamic modeling and optimal control strategy of waste heat recovery organic Rankine cycles[END_REF], where an energy management system would provide set points to a low-level controller. Future work will be directed toward the experimental assessment of the overall strategy, combining the low-level controller proposed herein with the above-mentioned energy management system.

APPENDIX PROOF OF PROPOSITION 1

Since ϕ = 1 -L 1 -L 2 -L 3 , to prove that (∂ϕ/∂ p) < 0, it suffices to prove the following conditions:

1) (∂ L 3 /∂ p) > 0, 2) (∂ L 2 /∂ p) > 0, 3) (∂ L 1 /∂ p) > 0.
1) From the algebraic system (7), the following equation is deduced:

L 3 = N 3 D 3 with

Fig. 1 .

 1 Fig. 1. Rankine WHR system layout for automotive application.

Fig. 3 .

 3 Fig. 3. ORC system inputs-outputs.

Fig. 4 .

 4 Fig. 4. Rankine cycle for a dry fluid.

Fig. 5 .

 5 Fig. 5. MBs layout for the evaporator.

Fig. 6 .

 6 Fig. 6. Controller with inverse of the nonlinear control-oriented model in the feedforward path.

Fig

  Fig. SH-N pump pair identification on two different operating points, using under damped second-order transfer functions.

Fig. 8 .

 8 Fig. 8. Picture of the Rankine system plugged on a heavy-duty diesel engine.

Fig. 9 .

 9 Fig. 9. Experimental result. Dynamic inverse model improve SH control. (a) Disturbances (exhaust gas conditions). (b) Control performance: PID (alone). (c) Control performance: PID with dynamic inverse model in feedforward path.

Fig. 11 .

 11 Fig. 11. Proposed closed-loop structure for control of SH and p. Observer feeds state estimation to model-based controller.

Fig. 12 .

 12 Fig. 12.Simulation results under varying inlets. Reference model (solid lines), observer (dashed lines), and reference model with observer initial conditions (dotted lines).

Fig. 13 .

 13 Fig. 13. Zoomed-in pressure plot of Fig. 12.

Fig. 14

 14 Fig. 14. Estimation (dashed lines) of reference model unmeasurable states (solid lines).

Fig. 15 .Fig. 16 .

 1516 Fig. 15. Pressure set point p SP tracking and SH regulation around SH SP .

Fig. 17

 17 Fig. 17. SH regulation error (SH-SH SP ) with observer (solid blue line) and with state estimation performed by an open-loop model (dashed green line).

  Fig. 17. SH regulation error (SH-SH SP ) with observer (solid blue line) and with state estimation performed by an open-loop model (dashed green line).

N 3 =

 3 C d S turb (2ρ( p, SH) p) 1/2 [h out ( p, SH)h v ( p)] D 3 = S f α 3 [T w,3 -T f,3 ( p, SH)].

Fig. 18 .D 2 with N 2 = 1 with N 1 =

 182211 Fig.[START_REF] Willems | Integrated energy & emission management for heavy-duty diesel engines with waste heat recovery system[END_REF].Functions N 1 , N 2 , and N 3 increase with pressure p, for the considered working fluid R245fa and SH = 30 K .

TABLE I NOMENCLATURE

 I Fig. 2. Schematic of the ORC system in the experimental setup.

In the following, the implicit time dependencies are dropped for compactness.

The numerical values considered for this application (namely for the organic fluid R245fa) are specified in Proposition 1.

To simplify the presentation, the following notations are used : ϕ p = (∂ϕ/ p)(T w , p, d) and ϕ Tw = (∂ϕ/T w )(T w , p, d).

Notice that the experimental setup was not available when this part of the study was conducted.


Didier Bodin

Gilles Pijaudier-Cabot

Chantal De

La Roche

Jean-Michel Piau

Armelle Chabot

Continuum Damage Approach to Asphalt Concrete Fatigue Modeling

Keywords: Asphaltic concrete, Fatigue, Damage, Cracking, Pavements, Constitutive relations

A nonlocal damage model is proposed to predict the behavior of pavement fatigue cracking. This constitutive relation has been implemented in a finite-element code, along with a selfadaptive jump-in-cycle procedure for high cycle fatigue computations. Strain localization analysis shows that during uniaxial fatigue tests, bifurcation due to strain softening occurs much later than in monotonic tests. The incorporation of an internal length into the constitutive model is advocated since the model should encompass loading histories with very different amplitudes of cycles, in which localization may still occur. The influence of the internal length on the fatigue life of bending beams is also investigated. Calibration of the damage model is performed after thermal effects have been evaluated and accounted for in a simplified way, uncoupled to damage. Parameter identification is performed in bending and uniaxial tests. The resulting calibrated constitutive relation is found to yield a good description of several different uniaxial tests.

Introduction

Repeated stress and strain induced by traffic loads can be very damaging to asphalt pavements. Cyclic mechanical loading leads to fatigue damage in the bottom layers. Accumulation of these small degradations can initiate macrocrack propagation through the pavement body from the bottom to the top. Design methods have to be developed to avoid the initiation and propagation of cracks, at least over the life of service of the structures. For asphalt pavements, design methods often rely on standardized material testing. In the French design method, the principle is to compare the strain level calculated at the bottom of the asphalt layer to admissible strain, which is strain leading to failure for 1 million cycles during a standardized fatigue test, modified by a shift factor that takes into account in situ conditions that are more complex than in laboratory test conditions ͑healing, temperature, frequency, etc.͒.

Most laboratory fatigue tests consist of applying constant amplitude sinusoidal displacement or force at the boundary of a bituminous mixture sample ͑see, e.g., SHRP 1994;[START_REF] Di Benedetto | Fatigue of bituminous mixtures: Different approaches and RILEM interlaboratory tests[END_REF][START_REF] Partl | Current research projects of the next RILEM TC PEB performance testing and evaluation of bituminous materials[END_REF]. During the tests, the variation in global stiffness is monitored. It is defined as the ratio of the force amplitude applied to the displacement amplitude. Because it is difficult to relate laboratory test results to field data, such experiments provide, in most cases, groundwork for comparisons between several different types of materials only ͑mix, granulometry, etc.͒. Rest periods can also be integrated into continuous fatigue tests to investigate more realistic fatigue performance ͑Kim and Little 1990; [START_REF] Breysse | Fatigue tests on bituminous composites: Balance between damage and recovering[END_REF] that include healing phenomena which are not dealt with here.

Computational prediction of fatigue life is quite a challenge because permanent strain and variations in stiffness are due to several related physical processes. There is dispersed microcracking in the binder ͑Lefeuvre 2001͒, and viscoelastic behavior induces an increase of temperature due to viscous dissipation during cyclic tests associated with important thermosensitivity of the binder ͑de La Roche and Marsac 1996͒. During a typical test, the global stiffness, calculated as the ratio of the force to the displacement amplitude, decreases following three regimes of evolution ͑Fig. 1͒: During phase I, a fast decrease in stiffness is observed. It is followed by phase II that corresponds to a quasilinear decay in stiffness. Fracture occurs during phase III, due to damage acceleration and to ultimately to the propagation of macrocracks.

There are several approaches by which to model asphalt concrete damage. It can be treated directly with viscoelasticity coupled to damage ͑Lefeuvre et al. 2000͒. Kim and Little ͑1990͒, Lee and Kim ͑1998͒, and Lee et al. ͑2000͒ have proposed models in which damage and viscoelasticity are uncoupled, with the help of a corresponding principle given by Shapery ͑1984͒ that transforms viscoelastic analysis into elastic analysis. For sinusoidal loading histories, Alimami ͑1987͒, Piau ͑1989͒, and Ullidtz et al. ͑1997͒ provided an interpretation of the variation in material stiffness due to damage as variation of the complex modulus. It is this variation of the complex modulus due to damage that is considered here. Therefore the viscoelasticity of the material will not be considered in this paper.

This work focuses on the description of dispersed microcracks in material during fatigue load histories. The constitutive model is The influence of the internal length, which enters in the nonlocal approach, is investigated on a cantilever beam subjected to controlled displacement cycles. Full thermomechanical coupling is not considered, but the mechanical response of specimens is corrected in order to account for the variation in structural stiffness due to a rise of temperature in a simplified way. The principle behind the correction is outlined. Model calibration and comparisons with test data are also provided.

Mechanical Damage Modeling

The mathematical model used to describe mechanical damage is an elasticity based damage model for fatigue. It is inspired by the constitutive relations of Mazars ͑1984͒ and Pijaudier-Cabot and Bazant ͑1987͒. Compared to existing proposals, e.g., by Lee et al. ͑2000͒, the constitutive relations are simpler since they rely on a scalar, isotropic, damage model which is widely used for concrete and other quasibrittle materials. It should be mentioned that for concrete at least, isotropy of damage may not be a major restriction in terms of modeling capabilities. Although microcracks are always pretty much oriented with respect to the loading system, the isotropic damage approach provides reasonable predictions when damage is due to tensile strain or stress especially ͑Fichant et al. 1999͒. The influence of microcracking is introduced via the damage variable d that ranges from 0 to 1. The stress strain relation reads

i j ϭ 1ϩ E͑1Ϫd ͒ i j Ϫ E͑1Ϫd ͒ ͓ kk ␦ i j ͔ or i j ϭC i jkl kl ϭ͑1Ϫd ͒C i jkl 0 kl (1)
E and ϭYoung's modulus and Poisson ratio of the undamaged material, respectively; i j and i j ϭstrain and stress components; C i jkl and C i jkl 0 ϭdamaged and initial ͑elastic͒ secant stiffness of the material; and ␦ i j ϭKronecker symbol. The evolution of damage depends on the amount of tensile stress that the material experiences during mechanical loading. It is expected that repeated tensile stress is the main source of local decohesion in asphalt concrete. The damage growth criterion is based on a modified Rankine criterion ͑Patza ´k and Jira ´sek 2001͒ with zero threshold damage growth, in which is the average ¯of equivalent strain ˜defined as

˜ϭͱ ͚ iϭ1 3 ͩ ͗ i ͘ ϩ E͑1Ϫd ͒ ͪ 2 (2)
where ͗ ͘ ϩ ϭMacauley bracket; and i ϭprincipal stress. The av- erage equivalent strain ¯is ¯͑x ͒ϭ 1

V r ͑ x ͒ ͵ ⍀ ͑xϪs ͒ ˜͑s ͒ds, with V r ͑ x ͒ϭ ͵ ⍀ ͑xϪs ͒ds (3)
where ⍀ϭvolume of the structure; V r (x)ϭrepresentative volume at point x; and (xϪs)ϭweight function:

͑xϪs ͒ϭexp ͩ Ϫ 4ʈxϪsʈ 2 l c 2 ͪ (4)
and l c ϭinternal length of the nonlocal continuum. The rate of damage growth is defined as a function of the equivalent strain rate ͑Paas et al. 1993͒:

d ˙ϭ f ͑ d ͒ ¯␤͗ ថ ͘ ϩ (5)
where f (d)ϭfunction of damage, which will be defined later on, and ␤ϭmodel parameter.

When the state of strain is homogeneous over the specimen, e.g., in uniaxial tension or compression, the nonlocal equivalent strain, Eq. ͑3͒, is, by definition, equal to the local equivalent strain, Eq. ͑2͒. It is possible to integrate Eq. ͑5͒ over the history of loading, which is restricted in this case to displacement controlled cycles of constant amplitude:

dd f ͑ d ͒ ϭ ¯␤͗d ¯͘ϩ ⇒F͑d ͒ϭ ͵ 0 d dd f ͑ d ͒ ϭ ͚ cycle 1 N ¯a ␤ϩ1 ␤ϩ1 (6)
where ¯aϭamplitude of the equivalent strain over one cycle; and F(d)ϭscalar function of damage. Let us define fatigue life as the number of load cycles needed prior to reaching a critical value of damage d crit . Eq. ͑6͒ provides the corresponding allowable number of cycles N crit of constant amplitude ¯a :

N crit ϭ F͑d crit ͒͑ ␤ϩ1 ͒ ¯a ␤ϩ1 with F͑d crit ͒ϭconstant (7)
The value of d crit is understood as a parameter that is chosen by the user. It is related to the allowable relative decrease in stiffness of the material when subjected to uniaxial fatigue. In the French design method, this maximum decrease in stiffness could serve as an end of life criterion, with the help of corrective terms that take into account in situ conditions. From this expression, it can be also noticed that Ϫ(␤ϩ1) corresponds to the slope of the Wo ¨hler curve in a log-log diagram.

Fatigue tests at constant displacement amplitude have been carried out on a single bituminous mix ͑0/6 mm dense asphalt concrete containing 6.4% of 50/70 pure bitumen͒, part of a current interlaboratory test program within the framework of the RILEM TC 180 PEB ͑Partl et al. 2000͒. Tension compression tests have been performed on cylindrical specimens (height ϭ120 mm, diameterϭ80 mm) with sinusoidal displacement cycles corresponding to the following constant strain amplitudes: ͑80, 90, 100, 140, 150, 180 m/m͒. Bending tests on cantilever trapezoidal samples (heightϭ250 mm, large clamped base ϭ56 mm, small loaded baseϭ25 mm and thicknessϭ25 mm) have been also performed. These bending tests were displacement controlled. The displacement cycles were such that the maximum ͑absolute value͒ strain in the specimen, evaluated from beam theory, takes the following values: 140, 180, and 220 m/m. Both types of tests have been performed at temperature of 10°C and frequency of 10 Hz. For each test, the decrease in structural stiffness of the specimen was recorded as a function of the number of cycles. Fig. 2 shows a typical Wo ¨hler curve obtained for asphalt subjected to tension compression cycles of constant strain amplitude. In Fig. 2, the fatigue life is defined for a critical value, d crit ϭ0.3. The slope of the Wo ¨hler curve is 5 approximately. The fact that a straight line is recovered in a log-log diagram is a justification of the mathematical form of the damage evolution law in Eq. ͑5͒. The function F(d) needs, however, to be determined so that the evolution of stiffness with the number of cycles is recovered. This evolution of damage should be able to describe the three regimes of damage growth shown in Fig. 1. It is also quite important to devise a form of this function in which the duration of the three regimes of damage growth can be adjusted. The function used by Paas ͑1990͒ and Bodin et al. ͑2002͒,

f ͑ d ͒ϭCd ␣ (8)
does not permit such control. Fig. 3 shows a fit of the same test data as in Fig. 2 with this equation. The first two regimes are captured very well. In bending tests, the occurrence of the third regime, with significant acceleration of damage growth, cannot be described properly ͑see [START_REF] Bodin | A continuum damage approach of asphalt concrete fatigue tests[END_REF]. In fact, we observed that it was not possible to get a good description of the three phases of damage growth at the same time with the evolution of damage given by Eq. ͑8͒. In order to correct this weakness, a new function F(d) is defined as

F͑d ͒ϭ␣ 1 ͭ 1Ϫexp ͫͩ Ϫd ␣ 2 ͪ ␣ 3 ͬͮ (9)
where ␣ 1 , ␣ 2 , and ␣ 3 ϭthree model parameters. Fig. 4 shows the variation in material stiffness in a uniaxial fatigue test ͑sinusoidal displacement cycles with constant amplitude͒ for different values of these model parameters. These curves will be quite helpful later on for parameter calibration. The major influence of ␣ 1 is on fatigue life ͑here the number of cycles up to complete failure with zero stiffness͒. Parameter ␣ 2 essentially controls the duration of the three regimes, and keeps the fatigue life constant. The larger it is, the larger the growth of damage over a given number of cycles.

The third parameter, ␣ 3 , acts on the shape and steepness of the evolution of damage, in the second regime, where the rate of damage growth per cycle is almost linear. The larger ␣ 3 , the lower the slope and the greater the damage growth per cycle. The present constitutive relation has three parameters involved in evolution law for damage, in addition to the elastic constants and internal length that are in the nonlocal equivalent strain. Before describing calibration of the model, let us first discuss finiteelement implementation of this model in the context of high cycle fatigue computations.

Finite-Element Implementation

The above constitutive equations have been implemented into the finite-element code Castem 2000. The implementation does not differ from that of a standard, nonlocal isotropic damage model ͑see, e.g., Pijaudier-Cabot 1995͒. A secant stiffness algorithm is implemented, which avoids computation of the tangent material stiffness which has a large bandwith and is nonsymmetric ͑Jira ´sek and Patza ´k 2002͒. For high cycle fatigue computations, the fa- tigue life can be greater than 1 million cycles. Therefore, it does not seem to be realistic to perform a finite-element computation over the entire number of cycles. A minimum of 10 increments per cycle is needed to integrate the evolution of damage correctly over a single cycle ͑in the case of uniaxial tension compression͒ and the total number of increments for a typical complete fatigue computation would amount at least to 10 million, which is prohibitive for structural analysis. The so-called ''jump-in-cycle'' procedure devised by Lemaitre ͑1996͒ and improved by Peerlings et al. ͑2000͒ has been implemented.

Jump-in-Cycle Procedure

The jump-in-cycle method is based on projection of the increment of damage computed within one cycle, onto a finite number of cycles. The number of cycles, extended to a continuous variable, is taken as a second time scale, in addition to the one that is used to perform computation within a single cycle incrementally. From an initial state defined at a number of cycles NϪ1 and for a state of damage d defined at each material point at this time, computation of the increment damage over one cycle is performed. This increment of damage is then transformed into a rate of damage per cycle:

ץd ץN Ϸ ͵ cycle N d ˙dtϭG͑d͑N ͒͒ (10)
Integration of damage growth between cycle N and cycle N ϩ⌬N is approximated with the help of a trapezoidal rule:

d͑Nϩ⌬N ͒ϭd͑ N ͒ϩ 1 2 ͓G͑ d͑N ͒͒ϩG͑ d͑Nϩ⌬N ͔͒⌬N (11)
In the present form, integration is implicit since it requires the computation of G͓d(Nϩ⌬N)͔. This term is evaluated following Euler forward prediction,

G͑Nϩ⌬N ͒ϷG͑ N ͒ϩ ץG͑N ͒ ץN ⌬N (12)
When this expression is substituted in Eq. ͑11͒, which is linearized to first order in ⌬N, an explicit integration scheme is recovered along with evaluation of error e which is propagated from one increment of cycles ⌬N to the next one ͑see Peerlings et al.

2000͒. The increment of cycles ⌬N is adapted so as to limit this error. The increment of cycles is chosen so that the error is kept constant:

⌬Nϭ ץd ץG (13)
where ϭfixed parameter in the computation (ϭ5ϫ10 Ϫ2 ).

The increment of cycles is chosen so that Eq. ͑13͒ is satisfied at the Gauss point of the finite-element mesh where growth of damage over a cycle is largest. Near failure, ⌬N decreases greatly as the damage growth tends to become localized in a small region of the structure and eventually form a crack.

Strain and Damage Localization

The constitutive relations described earlier incorporate an internal length. This additional parameter is classically introduced in strain softening models in order to control the occurrence of strain and damage localization and to avoid strain that localizes in a Here we perform a simple, one-dimensional analysis of localization in order to show the influence of the internal length in fatigue problems.

Consider a one-dimensional specimen, in a state of homogeneous damage and strain, subjected to load cycles of amplitude a . This initial state is defined by the value of damage d 0 and strain 0 at an arbitrary point of a specific load cycle. The stress strain relation, in the rate format, about this initial state follows from Eq. ͑1͒:

˙ϭ͑1Ϫd 0 ͒E ˙Ϫd ˙E 0 ϭ͑1Ϫd 0 ͒E ϪE 0 f ͑ d 0 ͒ ¯0 ␤ 1 V r ͑ x ͒ ͵ ⍀ ͑xϪs ͒ 8 ͑ s ͒ds (14)
We consider an infinite specimen so that the normalized weight function is independent of the coordinate x. Furthermore, damage is assumed to grow during the small perturbation ͑the strain and strain rates are positive͒. Equilibrium requires that the incremental stress is homogeneous so one arrives at the following integrodifferential equation:

͑ 1Ϫd 0 ͒E d ˙͑x ͒ dx ϪE 0 f ͑ d 0 ͒ ¯0 ␤ 1 V r ͑ x ͒ ͵ ⍀ ͑xϪs ͒ d ˙͑s ͒ dx dsϭ0 (15)
where the local equivalent strain rate has been replaced by the local strain rate inside the integral because they are equal in the one-dimensional case. Strain localization occurs when bifurcation is observed in the rate problem. This means that Eq. ͑15͒ has several, nontrivial solutions.

One may look for harmonic solutions to this equation in terms of harmonic perturbations: u ˙(x)ϭA exp(Ϫix) where Aϭunknown constant; and ϭangular frequency. Substitution of this expression of the displacement rate into Eq. ͑15͒ yields the following characteristic equation:

A 2 exp͑Ϫix ͓͒͑1Ϫd 0 ͒EϪE 0 f ͑ d 0 ͒ ¯0 ␤ ¯͑l c , ͔͒ϭ0 (16)
in which ¯(l c ,)ϭFourier transform of the normalized weight function. A similar result was obtained by Pijaudier-Cabot and Benallal ͑1993͒ in a more general ͓three-dimensional ͑3-D͔͒ context. Eq. ͑16͒ shows that bifurcation occurs if

¯͑l c , ͒ϭ ͑ 1Ϫd 0 ͒E E 0 f ͑ d 0 ͒ ¯0 ␤ (17)
For a fixed value of internal length, this equation provides the angular frequency, or wavelength, of the localized mode. Note that if the internal length vanishes, which is the case in a standard local damage model, the wavelength of the localized mode is arbitrary when Eq. ͑16͒ is satisfied. It follows that the boundary value problem has an infinite number of solutions and is ill posed. Same as in monotonic loading, this is a sufficient reason for introducing an internal length in the model, since one cannot be sure that the localization condition will ever be satisfied for any arbitrary fatigue load history.

In the course of monotonic loading, bifurcation occurs at the onset of softening, at peak stress. The major difference between fatigue and monotonic loading is that the amplitude of the strain remains small during fatigue cycles. The denominator of the right hand-side term in Eq. ͑17͒ is always very small. The result is plotted in Fig. 5, which shows when localization may occur in the cases of fatigue and monotonic load histories. The values of the parameters in the evolution law for damage are provided in Table 1. We have considered in Fig. 5 the local damage case (l c ϭ0), which corresponds to the situation where localization occurs soonest in the course of the loading history ͑see Pijaudier-Cabot and Benallal 1993͒. As far as the bifurcation criterion is concerned, taking a nonzero value of the internal length would not change the comparison between monotonic and fatigue load histories substantially. The localized deformation modes would, however, be very different. Localization occurs when the load history ͑strain versus damage curve͒ intersects the bifurcation cri-terion. For a monotonic increase in strain, it occurs for 0 Ϸ8ϫ10 Ϫ3 and for a value of initial damage of d 0 Ϸ0.2. For fatigue cycles with amplitude of a ϭ140ϫ10 Ϫ6 , bifurcation occurs when damage is close to 1. If the amplitude increases, bifurcation occurs for smaller and smaller values of initial damage in the fatigue problem. Strain localization, at least in uniaxial tests, is almost avoided. Then, one might ask about the requirement of introducing an internal length in the model whose aim is to restore well posedness of the boundary value problem upon localization. Again, an internal length is still needed since it cannot be proved that for any general loading localization will never occur. The model should encompass fatigue loading histories with very different amplitudes, including large ones for which bifurcation due to softening will occur. In this last case, without an internal length, failure without energy dissipation might be encountered upon localization, the same as in usual monotonic loading situations. This feature, along with the ill posedness of boundary value problems at the onset of localization would render finite-element computations mathematically and physically meaningless.

Influence of Internal Length

Besides control of the localized mode at the onset of localization, which results from Eq. ͑17͒, the internal length still has an influence on the growth of damage in boundary value problems where the state of strain is not homogeneous. In order to exhibit this influence, we consider the case of a cantilever beam subjected to sinusoidal displacement of fixed amplitude at its free end. Fig. 6 shows the problem and the finite-element ͑FE͒ mesh that has been used for the plane stress computations. The model parameters are reported in Table 1. Fig. 7 shows the evolution of the overall stiffness of the beam as a function of the number of cycles for several values of internal length: l c ϭ0, 9, 18, and 36 mm. Whenever the state of strain is not homogeneous, there is a well known sensitivity of the model predictions to the internal length ͑see, e.g., [START_REF] Pijaudier-Cabot | Nonlocal damage theory[END_REF][START_REF] De Borst | Fundamental issues in finite element analysis of localisation of deformation[END_REF]. In a

␣ 1 ␣ 2 ␣ 3 ␤ 5ϫ10 Ϫ19 0.5 3.0 5
monotonic loading problem, the peak load increases and the softening slope decreases when the internal length increases. In a fatigue problem, Fig. 7 shows that the larger the internal length, the larger the fatigue life time because the strain, which is not homogeneous, is averaged over a larger area as this length is increased ͑consequently the maximum nonlocal equivalent strain decreases͒.

It should be stressed that such sensitivity cannot be used for calibrating the internal length from test data on a single specimen size. In fact, it has been shown for monotonic loads that there is an infinite set of model parameters, with arbitrary values of inter-nal length that can fit a single structural response ͑three-point bend tests͒ with the same very good accuracy ͑Le Belle ´go et al. 2003͒. To arrive at proper calibration of the model parameters, i.e., a unique set of model parameters, including the internal length, tests on several specimens subjected to different boundary conditions, with different geometries, or size effect tests must be available ͑Carmeliet 1999͒. Another possibility is to use an approximate value of the internal length. It is this second possibility that will be used next.

Calibration and Comparisons with Test Data

Bituminous materials exhibit viscoelastic behavior, which leads to energy dissipation during each loading cycle of the fatigue tests. This dissipated energy is transformed into heat, leading to an increase of temperature inside the specimen. Dissipated energy, thermal exchange at specimen boundaries, and thermal properties of the material introduce a coupled thermomechanical effect. In uniaxial tension, compression fatigue tests, for instance, a 1°C rise may cause 5% loss of stiffness for usual conditions of temperature and frequency. We will evaluate the influence of thermal effects in a simplified uncoupled way. Although the viscoelastic response of the material is not explicitly modeled in this work ͑since we look at variation of the complex modulus͒, its consequences in terms of energy dissipation and thermomechanical effects will be taken into account for a proper comparison with experiments.

Evaluation of the Temperature Increase during Fatigue Tests

Let us consider that the structural stiffness of the specimen K(N) measured at cycle N can be written as

K͑N ͒ K͑0 ͒ ϭ͑1ϪD th ͒͑ 1ϪD mech ͒ (18)
where the variations in specimen stiffness due to thermal effects and mechanical damage are denoted D th and D mech , respectively. Note that these variations result from a decrease of material stiffness due to the increase of temperature and from the distribution of damage d in the specimen. Evaluation of the influence of thermal effects D th will now be performed independently from the calculation of mechanical damage. This coupled thermomechanical problem without any damage has been computed with the finite-element program Ce ´sar-LCPC, and the results have been validated by infrared thermography measurements for bituminous mixtures in the past ͑de La Roche et al. 1998͒. The viscoelastic model used to describe bituminous mixture behavior is from work by Huet ͑1963͒ and Sayegh ͑1967͒. The model uses the ''equivalence principle'' between the frequency f and the temperature :

E͑ f , ͒ϭE͓ f ,͑ ͔͒ϭE real ϩiE imag ϭE 0 ϩ E ϱ ϪE 0 1ϩ␦͑i f t͒ Ϫk ϩ͑i f t͒ Ϫh (19)
E 0 and E ϱ ϭlimit of the complex modulus when f tends towards 0 and infinity, respectively. k and h exponents ͑such that 1Ͼh ϾkϾ0) are related to the behavior of the ratio E imag /E real when f tends towards infinity ͑resp. f tends towards 0͒. ␦ϭa dimensionless constant. Function ϭ()ϭthe characteristic time, a function of the temperature, which Huet and Sayegh suggested approximating with Arrhenius or Eyring type law. log() is usually defined as the following parabolic function:

For the mix studied, the model parameters are the following: E ϱ ϭ29,700 MPa, E 0 ϭ25 MPa, hϭ0.675, kϭ0.21, ␦ϭ2.16, A 0 ϭ1.71, A 1 ϭϪ0.393 and A 2 ϭ0.00226. For the fatigue test conditions ͑10°C, 10 Hz͒, the complex modulus of the mix is 11,600 MPa and its phase angle is 15.1°. The Poisson ratio is considered to remain constant and ϭ0.35. The thermal constants are as follows: heat conductivity 0.8 10 6 W/m/°C and volumetric heat capacity 1.6 J/m 3 /°C. The exchange coefficients that give the heat flow at the boundaries of the specimen are determined in order to reach good agreement between measured and calculated temperatures at the surface of the specimen. The results of the simulation on tension compression tests are compared to thermal gauge measurements at the surface of the samples and presented in Fig. 8 ͑for more details also see [START_REF] Bodin | Prediction of the Intrinsic damage during bituminous mixes fatigue tests[END_REF]. Table 2 summarizes the relative decrease of structural stiffness computed for uniaxial tension compression tests and for bending tests on cantilever beams at several amplitudes. It can be observed that the variation in the structural stiffness of the specimen due to the thermal effect is almost negligible in bending tests, and slightly higher in tension compression tests. It has to be pointed out that the thermal calculations did not integrate the decrease in stiffness due to damage that occurred during fatigue tests. This assumption of noncoupled effects between damage and temperature is less and less true as damage grows and energy dissipation due to damage is converted into heat.

Calibration of Fatigue Test Data

Calibration of the evolution law for damage, Eqs. ͑5͒, and ͑9͒, can be performed directly in tension compression tests and bending tests with small amplitudes ( a ϭ80ϫ10 Ϫ6 and 140ϫ10 Ϫ6 m/m, respectively͒. Strain localization analysis shows that it is possible to use the tension-compression tests under the assumption that the strain distribution is homogeneous throughout the specimen, since bifurcation occurs for values of damage that are close to 1. Unfortunately, the set of available experimental data in uniaxial tension does not exhibit the third regime of damage growth. The test control and the displacement transducers did not allow data acquisition during this very fast regime so experimental measurements are lacking. Hence, this regime of damage growth, which is important for calibration of parameter ␣ 1 , is not accessible for the calibration from uniaxial tests and other experiments; here bending beams subjected to fatigue are required.

With the help of the Matlab standard least square fitting function, the fit of the uniaxial experiments, for the smallest amplitude after thermal correction is performed, yields an initial set of material parameters. Uniaxial tests are used because integration of the evolution of damage can be performed analytically ͓Eqs. ͑7͒ and ͑9͔͒ and standard fitting techniques, which do not need any finite-element calculation, can be implemented without any difficulty.

This set of model parameters serves as initial input for finiteelement simulation of bending tests and is corrected in order to capture the third regime of damage growth properly. To do this, the sensitivity analyses performed earlier ͑Fig. 4͒ are helpful. Parameter ␣ 3 is kept equal to its value obtained from tension compression tests. Since parameter ␣ 1 has a direct influence on fatigue life ͑Fig. 4͒, it is modified so that the computations provide the correct number of cycles up to complete failure. In order to balance the influence of the variation of ␣ 1 , ␣ 2 is adjusted to keep the correct evolution of damage during phases I and II of damage growth ͑Bodin 2002͒. Fig. 9 shows the result of this calibration in the bending and uniaxial tests. The set of model parameters obtained is given in Table 3. For the asphalt mix studied, the internal length l c was chosen as equal to three times the maximum aggregate size, the same as for concrete ͑Bazant and Pijaudier-Cabot 1989͒. Since thermal effects are very small in the bending test ͑see Table 2͒, they are neglected.

In addition to this fit, the gray zone in Fig. 9 delineates lower and upper bound predictions that are a result of a fit of the extreme values in uniaxial tension compression ͓Fig. 9͑b͔͒. Compared to the previous fit, which provides average specimen responses ͑Table 3͒, parameters ␣ 1 and ␣ 2 were changed significantly by keeping the third model parameter almost constant. These lower and upper bounds show how sensitive model responses in the bending problem are, depending on inherent dispersion of uniaxial tests. It appears that bending tests are more appropriate for model calibration. Dispersion of the test data is smaller and data up to complete failure are more easily accessible. Such calibration relies on inverse analysis and finite-element computations. It is therefore more difficult to perform, and it still requires an initial set of model parameters to start an inverse analysis which reaches an optimal fit that can be obtained from uniaxial tests.

The simulation of uniaxial tests at amplitude 100ϫ10 Ϫ6 m/m is shown in Fig. 10, including the responses obtained for the upper and lower bound fits discussed above. The correction due to thermal effects was reported in order to exhibit the influence of the temperature which is rather limited in view of scattering of the experimental data. Overall, this simulation provides good agreement with test data. Strain localization analysis in a one-dimensional setting shows the difference between the failure processes in monotonic loading and fatigue loading. In monotonic loading, bifurcation occurs at peak stress, while it occurs very close to complete failure ͑damage close to 1͒ in fatigue tests. Therefore, according to the proposed constitutive relations, the third regime of damage, with an increasingly growing rate of damage near failure, is not solely due to structural effects inherent to localization of damage.

An internal length is still introduced to ensure the consistency of the constitutive relations upon localization of damage, which may still occur when the amplitude of loading is very high as, for instance, in the case of notched specimens. The internal length controls the size of the damage localization zone, and the larger this parameter, the longer the fatigue life of bending beams.

The model parameters were calibrated on bending beams and uniaxial test data. Prior to calibration, the thermal part of the decrease in stiffness due to thermomechanical effects was evaluated by a simplified approach. For the material studied and laboratory fatigue conditions tested ͑10°C, 10 Hz͒, the variation of stiffness due to thermal effects remains under 4% of the structural stiffness of the specimens when the maximum strain value is below 100ϫ10 Ϫ6 for tension compression tests and 220ϫ10 Ϫ6 for cantilever tests. The model provides reasonable simulations of tension-compression tests for several amplitudes of loading.

It should be pointed out that, in the calibration process, the internal length was chosen as a function of the maximum aggregate size, following experimental data obtained for concrete. This estimate of the internal length needs to be further justified, since the microcracking process in bituminous mixes may be quite different from that in concrete. Size effect tests on geometrically similar cantilever beams should be able to provide a better estimate. Since the internal length controls the size of the damage zone, the response of geometrically similar beams should not be geometrically similar. Inverse analysis of such size effect tests should provide a more accurate determination of this model parameter.

1

  Laboratoire Central des Ponts et Chaussées, Route de Bouaye, B.P. 4129, F-44341 Bouguenais Cedex, France 2 R&DO-GeM, Institut de Recherche en Génie Civil et Mécanique, Ecole Centrale de Nantes, B.P. 92101, F-44321 Nantes Cedex 3, France based on an elastic isotropic continuum damage model, which expresses the decrease of elasticity parameters with loading. A simple model has been chosen and implemented in a finiteelement code in order to highlight performance trends. The complete model is described next in this paper. Then the problem of high cycle fatigue computations in the context of finite-element analyses is given. The method of Peerlings et al. ͑2000͒ is implemented in the finite-element code Castem 2000. In order to capture ultimate failure and strain localization consistently, a nonlocal approach to the damage model is advocated.

Fig. 1 .

 1 Fig. 1. Evolution of stiffness during typical fatigue test. Notice the three regimes: in the first damage evolution decelerates and the second is almost linear. In the third, damage growth accelerates.

Fig. 2 .

 2 Fig. 2. Wo ¨hler curve for typical bituminous mix loaded in uniaxial tension compression

  region of zero volume with a vanishing energy dissipation at failure ͑Pijaudier-Cabot and Bazant 1987; de Borst et al. 1993͒. In addition, nonlocal constitutive relations prevent spurious mesh dependence and exhibit a structural size effect, which in the case of monotonic loads provides a good indicator for obtaining the internal length from experiments ͑see Le Belle ´go et al. 2003͒.

Fig. 4 .

 4 Fig. 4. Sensitivity of parameters in damage model to uniaxial fatigue (␣ 1 top, ␣ 2 middle, ␣ 3 bottom͒

Fig. 5 .

 5 Fig. 5. Evolution of localization criterion during monotonic and fatigue uniaxial tests

Fig. 6 .Fig. 7 .

 67 Fig. 6. Cantilever beam subjected to high cycling fatigue. Geometry of the specimen ͑a͒ and finite-element mesh ͑b͒. The specimen is clamped at the bottom and the upper end is subjected to sinusoidal horizontal displacement.

Fig. 8 .

 8 Fig. 8. Increase in sample surface temperature as a function of the number of cycles: experimental data and FE simulations for tension-compression tests ͑after Bodin et al. 2003͒

Conclusion

  An elasticity based nonlocal damage model for fatigue failure of asphalt concrete was presented. Compared to the existing model of Paas et al. ͑1993͒, a different form of evolution law for damage was implemented in order to capture more accurately the third regime of damage growth, just prior to complete failure. This constitutive relation was implemented in a finite-element code, along with a jump-in-cycle method for finite-element computations of high cycle fatigue tests.

Fig. 9 .Fig. 10 .

 910 Fig. 9. Calibration of model parameters from cantilever bending tests ͑top͒ and upper and lower bounds from uniaxial test for amplitude of 80ϫ10 Ϫ6 m ͑bottom͒; gray zone shows lower and upper bound fits

Table 1 .

 1 Model Parameters for Computation of Cantilever Beam for Several Internal Lengths

Table 2 .

 2 Maximum Correction Factors due to Thermal Effects

	Variable	Tension compression	Two-point bending
	max (10 Ϫ6 )	80	100	150	180	140	180	220
	D th (10 Ϫ2 )	2.05	3.2	6.4	10.7	1.7	2.4	3.9

Table 3 .

 3 Set of Model Parameters Obtained from Calibration

	Limit	␤	␣ 1 (10 Ϫ19 )	␣ 2 (10 Ϫ1 )	␣ 3
	Upper	5	4.73	5.04	3.07
	Average fit		6.84	5.0	3.00
	Bottom		15.8	5.43	3.10

Acknowledgments

The experimental results presented were performed within the framework of a scientific collaboration between the Laboratoire Central des Ponts et Chaussee ´s ͑LCPC͒ and ENTPE, France. They are part of an international interlaboratory program for fatigue tests of bituminous mixtures conducted under the auspices of RILEM, TC 180 PEB. This work was part of LCPC research project Fatigue Damage of Bituminous Mixtures.