Romain Postoyan
email: romain.postoyan@univ-lorraine.fr.

Lucian Busoniu

Dragan Nesic
email: dnesic@unimelb.edu.au.

Jamal Daafouz
email: jamal.daafouz@univ-lorraine.fr.

L Bus ¸oniu

D Nešić

Stability of infinite-horizon optimal control with discounted cost

come   L'archive ouverte pluridisciplinaire

I. INTRODUCTION

Optimal control selects control inputs so as to minimize a cost incurred during the system operation ( [START_REF] Kirk | Optimal Control Theory: An Introduction[END_REF]). In this paper, we focus on optimal control in discrete time over an infinite horizon, with general nonlinear system dynamics as well as general stage costs. In this general setting, optimal control is a very powerful framework ( [START_REF] Bertsekas | Stochastic Optimal Control: The Discrete Time Case[END_REF]), able to address decision-making problems not only in control engineering, but also in artificial intelligence, operations research, economy, medicine, etc. In discounted optimal control, stage costs are weighted by an exponentially decreasing term γ k , where γ ∈ (0, 1) is the discount factor and k is the time step. The discounted setting is popular in many fields, such as reinforcement learning [START_REF] Sutton | Reinforcement Learning: An Introduction[END_REF], [START_REF] Szepesvári | Algorithms for Reinforcement Learning[END_REF], [START_REF] Bus ¸oniu | Reinforcement Learning and Dynamic Programming Using Function Approximators[END_REF] and dynamic programming [START_REF] Lewis | Reinforcement Learning and Adaptive Dynamic Programming for Feedback Control[END_REF]. For example, a novel class of predictive optimistic algorithms [START_REF] Munos | The optimistic principle applied to games, optimization and planning: towards foundations of Monte-Carlo tree search[END_REF], which provides a very useful relationship between computation invested and nearoptimality, only makes sense in the discounted case.

A core practical question is whether the discounted optimal control law stabilizes the system. Some results exist in the economy literature, showing local stability in the continuous-time case e.g., [START_REF] Rodriguez | On the local stability of the solution to optimal control problems[END_REF], [START_REF] Sorger | Local stability of stationary states in discounted optimal control systems[END_REF], while, in the dynamic programming field, the analysis is usually tailored to the specific cost function considered or the specific algorithm used, see for instance [START_REF] Balakrishnan | Issues on stability of ADP feedback controllers for dynamical systems[END_REF], [START_REF] Boussios | Semiglobal nonlinear stabilization via approximate policy iteration[END_REF], [START_REF] Sokolov | Complete stability analysis of a heuristic ADP control design[END_REF]. To the best of our knowledge, the stability properties in the general discounted case are not yet understood. The authors of [START_REF] Kiumarsi | Reinforcement Q-learning for optimal tracking control of linear discrete-time systems with unknown dynamics[END_REF] highlight the difficulties remaining even in the linear case. The main difficulty comes from the impact of the value of γ on the system stability. We will show that even for a simple linear example, γ needs to be sufficiently close to 1 to ensure stability.

Motivated by this insight, we develop a stability analysis for discounted optimal control with general nonlinear dynamics and stage costs, building on the undiscounted problem (γ = 1) as a nominal case. Contrary to most aforementioned references, we investigate global asymptotic stability properties. Furthermore, we define stability using a generic measure as in [START_REF] Grimm | Model predictive control: for want of a local control Lyapunov function, all is not lost[END_REF], which allows one to address the classical equilibrium point stability as a particular case, but also set stability for instance. We first provide sufficient conditions to ensure a global asymptotic stability property for the undiscounted problem, like in [START_REF] Kreisselmeier | Numerical nonlinear regulator design[END_REF]. Then, our main result guarantees that this stability property is preserved semiglobally and practically where the adjustable parameter is γ. In other words, for any basin of attraction, the system solutions initialized in this basin will converge to any neighborhood of the target set provided γ is sufficiently close to 1. We only impose standard continuity and compactness conditions on the problem ( [START_REF] Bertsekas | Stochastic Optimal Control: The Discrete Time Case[END_REF], [START_REF] Hernández-Lerma | Further topics on discretetime Markov control processes[END_REF]), together with a controllability property of the system as in [START_REF] Grimm | Model predictive control: for want of a local control Lyapunov function, all is not lost[END_REF]. Our analysis is inspired by [START_REF] Grimm | Model predictive control: for want of a local control Lyapunov function, all is not lost[END_REF], where the undiscounted finite-horizon case was considered.

Often, the stage cost is taken bounded, see for instance [START_REF] Munos | The optimistic principle applied to games, optimization and planning: towards foundations of Monte-Carlo tree search[END_REF], [START_REF] Sutton | Reinforcement Learning: An Introduction[END_REF], [START_REF] Szepesvári | Algorithms for Reinforcement Learning[END_REF], [START_REF] Bus ¸oniu | Reinforcement Learning and Dynamic Programming Using Function Approximators[END_REF] and Chapter 4 of [START_REF] Bertsekas | Stochastic Optimal Control: The Discrete Time Case[END_REF]. This has important practical advantages, such as ensuring bounded value functions, as well as a contraction property that guarantees the exponential convergence of important algorithms such as value iteration and policy evaluation. By strengthening the controllability requirement, we explain how our framework applies to a type of bounded stage cost, obtained by saturating an originally unbounded function. Finally, for general continuous bounded stage costs, we show that the Lyapunov (value) function is continuous. This property is essential to endow the system with some nominal robustness, see [START_REF] Grimm | Examples when nonlinear model predictive control is nonrobust[END_REF], [START_REF] Kellett | On the robustness of KL-stability for difference inclusions: smooth discrete-time Lyapunov functions[END_REF]. Note that when stage costs can grow unbounded, the value function is usually only shown to be lower semicontinuous [START_REF] Hernández-Lerma | Further topics on discretetime Markov control processes[END_REF], [START_REF] Cavazos-Cadena | The discounted method and equivalence of average criteria for risk-sensitive Markov decision processes on Borel spaces[END_REF], while e.g., [START_REF] Blume | Characterization of optimal plans for stochastic dynamic programs[END_REF] shows continuity for the specific case of concave stage costs.

Next, after introducing some preliminaries in Section II and stating the problem in Section III, we provide the nominal, undiscounted analysis in Section IV. Stability with discount is studied in Section V, while Section VI specializes the analysis to bounded stage costs, and Section VII shows the continuity of the Lyapunov function. Conclusions are given in Section VIII.

II. PRELIMINARIES Let R = (-∞, ∞), R := R ∪ {-∞, ∞}, R ≥0 = [0, ∞), R >0 = (0, ∞), Z ≥0 = {0, 1, 2, . . .}, Z >0 = {1, 2, . . .}. For (x, y) ∈ R n+m , (x, y) stands for [x T , y T ] T . A function χ : R ≥0 → R ≥0 is of class K if it
is continuous, zero at zero and strictly increasing, and it is of class The definitions below can be found in [START_REF] Rockafellar | Variational analysis[END_REF]. A function

K ∞ if in addition it is unbounded. A continuous function χ : R 2 ≥0 -→ R ≥0 is of class KL if for each t ∈ R ≥0 , χ(•, t) is of class K,
f : R n → R is lower semicontinuous at x ∈ R n if lim inf x→x f (x) := lim δց0 [ inf x∈B(x,δ) f (x)] = f (x), where B(x, δ)
is the closed ball of R n , centered at x of radius δ ≥ 0. We say that f is lower semicontinuous on X ⊆ R n , when it is lower semicontinuous at any x ∈ X. Note that when f is continuous, it is also lower semicontinuous.

Let S : R n ⇉ R m denote a set-valued mapping. The domain of S is defined by dom S := {x ∈ R n : S(x) = ∅}. We recall the notions of outer and inner limits of S at x ∈ R n : lim sup

S(x).

The set-valued mapping S is continuous at x ∈ R n when S(x) → S(x) as x → x, and it is continuous on X ⊆ R n when it is continuous at any x ∈ X. Define the image of a set V under the mapping S by S(V) = x∈V S(x). The mapping S : R n ⇉ R m is locally bounded when for any x ∈ R n , for some neighborhood V of x, the set S(V) ⊂ R m is bounded.

III. PROBLEM STATEMENT

Consider the system

x(k + 1) = f (x(k), u(k)) (1)
with state x ∈ R n and input u ∈ U(x) ⊆ R m , where U(x) is a non-empty set of admissible inputs associated to state x (as in [START_REF] Hernández-Lerma | Further topics on discretetime Markov control processes[END_REF], [START_REF] Keerthi | An existence theorem for discretetime infinite-horizon optimal control problems[END_REF] for instance), and n, m ∈ Z >0 . Let W := {(x, u) : x ∈ R n and u ∈ U(x)}. Define an infinite-length sequence of control inputs u = (u 0 , u 1 , . . .), in which the control input at time k ∈ Z ≥0 is given by u(k) = u k . We denote the solution to (1) at the k th -step starting at state x and with the input sequence u as φ(k, x, u| k ), where u| k := (u 0 , . . . , u k-1 ) is the truncation of u to the first k ∈ Z >0 steps, and we use the convention φ(0, x, u| 0 ) = x. We focus on the scenario where the sequence of control inputs is given by the solution to an optimization problem. In particular, we consider the cost function

J(x, u) := ∞ k=0 γ k ℓ(φ(k, x, u| k ), u k ), (2)
where ℓ : W → R is the stage cost and γ ∈ (0, 1) is the discount factor. The objective is to minimize the cost J.

When inf u J(x, u) is well-defined (which is assumed to be the case in this study), we call it the optimal value function and we denote it

V (x) := inf u J(x, u). (3)
Further, when the infimum in ( 3) is achieved by an infinitelength input sequence, we denote this sequence u * (x) and call it the optimal solution at state x. Thus,

V (x) = J(x, u * (x)).
The objective of this study is to infer about the stability property of system (1) when controlled by the optimal sequence u * (x) (which may be non-unique), when the latter exists. Our approach consists in considering, as a nominal case, the corresponding undiscounted problem, i.e. when the cost function is given by ( 2) with γ = 1

J(x, u) := ∞ k=0 ℓ(φ(k, x, u| k ), u k ), (4)
and the optimal value function is V (x) := inf u J(x, u). In this undiscounted setting, we can guarantee the stability of the system with an optimal sequence of inputs u * (x) (which will be ensured to exist), defined as a sequence so that V (x) = J(x, u * (x)). We then study under which conditions it is possible to (approximately) maintain this stability property in the discounted case, when γ ∈ (0, 1).

Our approach is motivated by the study of a simple example which shows that γ indeed needs to be appropriately selected to preserve stability.

Example 1: Consider the scalar system x(k 2 , where γ ∈ (0, 1) and u is an infinite-length sequence of inputs. The optimal solution is given by the feedback law

+ 1) = 2x(k) + u(k) and the discounted quadratic cost J(x, u) = ∞ k=0 γ k x(k) 2 + u(k)
u * = Kx with K = -2 1 + 2 5γ -1 + (5γ -1) 2 + 4γ -1 -1 (see Sec- tion 4.2 in [2]
). The origin of the closed-loop system is globally exponentially stable if and only if 2 + K ∈ (-1, 1), which is equivalent to γ ∈ (γ ⋆ , 1) where γ ⋆ = 1 3 . Hence, γ needs to be large enough, otherwise the optimal feedback law does not stabilize the system.

IV. PRELIMINARIES ON THE UNDISCOUNTED PROBLEM

We make the following assumption on the undiscounted problem, which guarantees the existence of an optimal solution according to [START_REF] Keerthi | An existence theorem for discretetime infinite-horizon optimal control problems[END_REF].

Assumption 1: The following holds. (i) f is continuous on W and ℓ is lower semicontinuous and nonnegative on W. (ii) For any compact set C, the set {(x, u) ∈ W : x ∈ C} is compact. (iii) For any x ∈ R n , there exists an admissible sequence of inputs u such that J(x, u) is bounded. Assumption 1 corresponds to the conditions of Theorem 1 in [START_REF] Keerthi | An existence theorem for discretetime infinite-horizon optimal control problems[END_REF]. Item (ii) of Assumption 1 is satisfied e.g., when U(x) = U for any x ∈ R n with U compact; more general conditions are given in Theorem 2 in [START_REF] Keerthi | An existence theorem for discretetime infinite-horizon optimal control problems[END_REF]. Item (iii) of Assumption 1 is verified for instance when ℓ satisfies a general exponential controllability property introduced in Definition 3 in [START_REF] Grimm | Model predictive control: for want of a local control Lyapunov function, all is not lost[END_REF], as will be formalized in the lemma below. This property is stated with respect to a generic continuous, positive definite function σ : R → R ≥0 , which will serve later as a state measure (as in [START_REF] Grimm | Model predictive control: for want of a local control Lyapunov function, all is not lost[END_REF]). The function σ can be defined e.g., as | • |, | • | 2 or by the distance to a set.

Lemma 1: Consider system (1) and suppose that ℓ is globally exponentially controllable to zero with respect to σ : R n → R ≥0 , i.e. there exist M, λ > 0 such that for any x ∈ R n there exists an admissible infinite-length control input sequence u such that ℓ(φ(k, x, u| k ), u k ) ≤ M σ(x)e -λk for any k ∈ Z >0 . Then item (iii) of Assumption 1 holds. Proof. Let x ∈ R n and take the sequence u that satisfies the condition of Lemma 1. It then holds that, for any

N ∈ Z >0 , N k=0 ℓ(φ(k, x, u| k ), u k ) ≤ N k=0 M σ(x)e -λk ≤ M σ(x) 1 -e -λ . (5)
Since the inequalities above hold for any N ∈ Z >0 , they also hold when N → ∞, and

J(x, u) = ∞ k=0 ℓ(φ(k, x, u| k ), u k ) ≤ M σ(x)
1e -λ , which implies the satisfaction of item (iii) of Assumption 1.

The conditions below come from [START_REF] Grimm | Model predictive control: for want of a local control Lyapunov function, all is not lost[END_REF] where the stability of model predictive control problems with finite-horizon undiscounted cost is investigated. These same conditions will allow us to guarantee a global asymptotic stability property for system (1) when the cost is infinite-horizon and undiscounted, similarly to [START_REF] Kreisselmeier | Numerical nonlinear regulator design[END_REF].

Assumption 2: Let σ : R n → R ≥0 be continuous and positive definite. There exist ρ, α ∈ K ∞ such that the following holds for any (x, u) ∈ W

V (x) ≤ α(σ(x)) and ℓ(x, u) ≥ ρ(σ(x)). (6)
The first inequality in ( 6) is verified e.g., when ℓ is globally exponentially controllable to zero with respect to σ : R n → R ≥0 . Indeed, in that case, we have that

V (x) ≤ α(σ(x)) with α(σ(x)) = Mσ(x)
1-e -λ in view of the proof of Lemma 1. The second inequality in ( 6) is a detectability property of σ from ℓ (see Definition 1 in [START_REF] Grimm | Model predictive control: for want of a local control Lyapunov function, all is not lost[END_REF]), which is satisfied for instance when ℓ is quadratic in its first argument and σ(•) = | • |. This detectability property can be relaxed as in SA3 in [START_REF] Grimm | Model predictive control: for want of a local control Lyapunov function, all is not lost[END_REF]; we will address this extension in future work.

The theorem below formalizes the existence of an achievable sequence which minimizes the cost function in (4) (item (a)) and states a global asymptotic stability property for system (1) when this sequence is applied (item (b)).

Theorem 1: Consider system (1), the cost function J in (4), and suppose Assumptions 1-2 hold. Then the following holds.

(a) For any x ∈ R n there exists an optimal sequence of inputs u * (x). (b) There exists β ∈ KL such that for any x ∈ R n and

k ∈ Z ≥0 , σ(φ(k, x, u * (x)| k )) ≤ β(σ(x), k).
Proof. Assumption 1 allows us to directly apply Theorem 1 in [START_REF] Keerthi | An existence theorem for discretetime infinite-horizon optimal control problems[END_REF] which ensures the satisfaction of item (a) of Theorem 1. Let x ∈ R n and consider u * (x). According to the Bellman equation,

V (x) = ℓ(x, u * 0 ) + V (φ(1, x, u * (x)| 1 )). (7)
Hence, in view of Assumption 2,

V (φ(1, x, u * (x)| 1 )) -V (x) = -ℓ(x, u * 0 ) ≤ -ρ(σ(x)) ≤ -ρ • α -1 (V (x)).
(8) We obtain by induction, for any k ∈ Z ≥0 ,

V (φ(k + 1, x, u * (x)| k+1 )) -V (φ(k, x, u * (x)| k )) ≤ -ρ • α -1 (V (φ(k, x, u * (x)| k ))). (9)
As a consequence, according to [START_REF] Nešić | A framework for stabilization of nonlinear sampled-data systems based on their approximate discrete-time[END_REF], there exists β ∈ KL such that for any

k ∈ Z >0 , V (φ(k, x, u * (x)| k )) ≤ β(V (x), k). Using the fact that ρ(σ(x)) ≤ ℓ(x, u * 0 ) ≤ V (x) ≤ α(σ(x)) in view of Assumption 2, we deduce that item (b) of Theorem 1 holds with β(s, k) = ρ -1 • (β(α(s), k)) for any s, k ∈ R ≥0 .
Item (b) of Theorem 3 is equivalent to the usual (ε, δ)definition of global asymptotic stability and is standard in the nonlinear systems literature. It means that the state remains within a region that depends on the magnitude of the initial state, but shrinks with time.

V. ANALYSIS OF THE DISCOUNTED PROBLEM

Even in the discounted case, Assumption 1 ensures the existence of an optimal sequence of inputs for any initial condition, according to the proposition below.

Proposition 1: Consider system (1), the cost function J in (2), and suppose Assumption 1 holds. For any x ∈ R n there exists an optimal sequence of inputs u * (x). Proof. Let γ ∈ (0, 1). As before, conditions a)-d) of Theorem 1 in [START_REF] Keerthi | An existence theorem for discretetime infinite-horizon optimal control problems[END_REF] are verified in view of items (i)-(ii) of Assumption 1. Let x ∈ R n . According to item (iii) of Assumption 1, there exists u such that J (x, u) is finite. Hence, since γ ∈ (0, 1) and ℓ is nonnegative in view of item (i) of Assumption 1, for any

N ∈ Z >0 0 ≤ N k=0 γ k ℓ(φ(k, x, u| k ), u k ) ≤ N k=0 ℓ(φ(k, x, u| k ), u k ) ≤ J(x, u). (10)
Since these inequalities hold for any N , they also hold as N → ∞. We deduce that J(x, u) ∈ [0, J(x, u)]. Hence, J(x, u) is bounded and condition e) of Theorem 1 in [START_REF] Keerthi | An existence theorem for discretetime infinite-horizon optimal control problems[END_REF] holds. We apply Theorem 1 in [START_REF] Keerthi | An existence theorem for discretetime infinite-horizon optimal control problems[END_REF] to conclude that the infimum in (3) is well-defined and achievable for any x ∈ R n .

The theorem below gives Lyapunov-based properties, from which we will later deduce a stability property for system (1) when applying an optimal sequence of inputs.

Theorem 2: Consider system (1) and the cost function J in (2). The following holds under Assumptions 1-2.

(a) For any x ∈ R n and any γ ∈ (0, 1),

ρ(σ(x)) ≤ V (x) ≤ α(σ(x)). (11)
(b) For any γ m ∈ (0, 1), there exists Υ γm ∈ KL such that for any x ∈ R n , γ ∈ (γ m , 1), and

k ∈ Z ≥0 V (x(k + 1)) -V (x(k)) ≤ -ρ(σ(x(k))) +Υ γm (σ(x(k)), γ 1-γ ), (12)
where x(k) = φ(k, x, u * (x)| k ) with x ∈ R n and u * (x) the associated optimal sequence. Proof. Let γ ∈ (0, 1), x ∈ R n , V (x) is well-defined and there exists an optimal sequence u * (x) = (u * 0 , u * 1 , . . .) according to Proposition 1. Since ℓ is nonnegative according to item (i) of Assumption 1 and (6) holds,

V (x) ≥ ℓ(x, u * 0 ) ≥ ρ(σ(x)). (13)
On the other hand, since γ ∈ (0, 1),

V (x) = inf u ∞ k=0 γ k ℓ(φ(k, x, u| k ), u k ) ≤ inf u ∞ k=0 ℓ(φ(k, x, u| k ), u k ) = V (x) (14)
and according to Assumption 2,

V (x) ≤ V (x) ≤ α(σ(x)). (15)
From ( 13) and ( 15), we have that item (a) of Theorem 2 holds. We now prove that item (b) of Theorem 2 holds. Let γ m ∈ (0, 1), γ ∈ (γ m , 1), x ∈ R n and u * 0 be defined as above. We denote φ(1, x, u * (x)| 1 ) = f (x, u * 0 (x)) as x + for the sake of convenience. According to the Bellman equation,

V (x) = ℓ(x, u * 0 ) + γV (x + ), (16)
therefore

V (x + ) -V (x) = -ℓ(x, u * 0 ) + (1 -γ)V (x + ). (17) Let x + := f (x, u * 0 (x))
where u * (x) = (u * 0 (x), u * 1 (x), . . .) is an optimal sequence of inputs for the undiscounted cost. We note that V

(x) = inf u J(x, u) ≤ J(x, (u * 0 (x), u * (x + )) = ℓ(x, u * 0 ) + γV (x + ). From (16), γV (x + ) ≤ ℓ(x, u * 0 ) + γV (x + ) ≤ ℓ(x, u * 0 ) + γV (x + ). Hence γV (x + ) ≤ α(σ(x)) + γV (x + ) from (15) and since ℓ(x, u * 0 ) ≤ V (x) ≤ α(σ(x)). From (9), it follows directly that V (x + ) ≤ (I -ρ • α -1 )(V (x)). Using V (x) ≤ α(σ(x))
from Assumption 2, we deduce that V (x + ) ≤ (α-ρ)(σ(x)). Consequently, γV (x + ) ≤ α(σ(x)) + γ(αρ)(σ(x)). Then, by using this expression in [START_REF] Lewis | Reinforcement Learning and Adaptive Dynamic Programming for Feedback Control[END_REF], we have

V (x + ) -V (x) ≤ -ρ(σ(x)) + (1 -γ) γ -1 α(σ(x)) + (α -ρ)(σ(x)) ≤ -ρ(σ(x)) + (γ -1 -1) 2α(σ(x)) -ρ(σ(x))
, where γ -1 ≥ 1 was used to obtain the last inequality. We assume without loss of generality that 2αρ ∈ K (we can always upper-bound 2αρ by a class-K function if this is not the case by following Lemma 4.3 in [START_REF] Khalil | Nonlinear systems[END_REF] as this function is zero at zero, nonnegative and continuous). We also note that [START_REF] Balakrishnan | Issues on stability of ADP feedback controllers for dynamical systems[END_REF]. By induction on the time, we deduce that item (b) of Theorem 2 holds with Υ γm (s, t) = min{t -1 , γ -1 m -1} 2α(s)ρ(s) for any s, t ≥ 0.

0 ≤ γ -1 -1 = min{γ -1 -1, γ -1 m -1} since γ ∈ (γ m ,
The following theorem means that, if the optimal solution for the undiscounted cost globally asymptotically stabilizes the system (in the sense that item (b) of Theorem 1 holds), then the discounted optimal solution will also stabilize the system but semiglobally and practically, where the adjustable parameter is γ.

Theorem 3: Consider system (1), the cost function J in (2), and suppose Assumptions 1-2 hold. There exists β ∈ KL such that for any δ, ∆ > 0, there exists γ ⋆ ∈ (0, 1) such that for any γ ∈ (γ ⋆ , 1) and x ∈ {z ∈ R n : σ(z) ≤ ∆}, the solution φ to system (1) with initial condition x and the optimal sequence of inputs u * (x) satisfies

σ(φ(k, x, u * (x)| k )) ≤ max{β(σ(x), k), δ} ∀k ∈ Z ≥0 . ( 18
)
Proof. We want to use similar arguments as in the proof of Corollary 1 in [START_REF] Grimm | Model predictive control: for want of a local control Lyapunov function, all is not lost[END_REF] to obtain the desired result. For that purpose, we need to consider an autonomous system. For any x ∈ R n , take the optimal sequence u * (x). We write system (1) with input sequence u * (x) as x(k + 1) = f (x(k), k).

Adding the time k as an extra variable, we obtain the autonomous system x = F (x) where x = (x, k) and F (x) = ( f (x, k), k + 1). We define Ṽ (x) := V (x) and σ(x) := σ(x) for any x = (x, k). From Theorem 2, it holds that, for any

γ m ∈ (0, 1), γ ∈ (γ m , 1), x ∈ R n+1 , ρ(σ(x)) ≤ Ṽ (x) ≤ α(σ(x)). (19)
and, for any k ∈ Z ≥0 ,

Ṽ ( φ(k + 1, x)) -V ( φ(k, x)) ≤ -ρ(σ( φ(k, x))) +Υ γm (σ( φ(k, x)), γ 1-γ ), (20)
where φ(k, x) is the solution to system x(k + 1) = F (x) at time k, with initial condition x. Like in the proof of Corollary 1 in [START_REF] Grimm | Model predictive control: for want of a local control Lyapunov function, all is not lost[END_REF], we then derive that for any δ, ∆ > 0, there exists γ ⋆ ∈ (γ m , 1) such that for any γ ∈ (γ ⋆ , 1) and x ∈ {z ∈ R n+1 : σ(z) ≤ ∆}, the solution φ to system x(k + 1) = F (x) with initial condition x satisfies

σ( φ(k, x)) ≤ max{β(σ(x), k), δ} ∀k ∈ Z ≥0 . (21)
Since σ(x) = σ(x) for any x = (x, k), we deduce from ( 21) that ( 18) holds. This theorem is one of the main results of the paper. It is interesting to note the analogy with the results in [START_REF] Grimm | Model predictive control: for want of a local control Lyapunov function, all is not lost[END_REF] where the finite horizon N takes the place of γ 1-γ in the infinitehorizon discounted problem. Informally, quantity γ 1-γ can be thought of as an 'effective horizon' of the discounted problem. Thus, while [START_REF] Grimm | Model predictive control: for want of a local control Lyapunov function, all is not lost[END_REF] shows stability for horizons N greater than some lower bound N ⋆ , we show it for γ > γ ⋆ , and thus for effective horizons larger than γ ⋆ 1-γ ⋆ . The value function V is used in this section as a Lyapunov function to analyse stability. In practice, it is essential to work with a Lyapunov function which is continuous in order to guarantee some nominal robustness for the system (see [START_REF] Kellett | On the robustness of KL-stability for difference inclusions: smooth discrete-time Lyapunov functions[END_REF]). The assumptions made so far do not a priori allow us to assert that V is continuous. In Section VII, we will provide conditions which will allow us to prove the continuity of V . These conditions include the boundedness of the stage cost, so we first explain, in Section VI below, how our stability guarantees can be applied to the case of bounded stage costs, which is also of practical importance in its own right.

VI. BOUNDED STAGE COST

In many cases, the stage cost is taken bounded, see for instance Chapter 4 in [START_REF] Bertsekas | Stochastic Optimal Control: The Discrete Time Case[END_REF] as well as [START_REF] Bus ¸oniu | Reinforcement Learning and Dynamic Programming Using Function Approximators[END_REF], [START_REF] Munos | The optimistic principle applied to games, optimization and planning: towards foundations of Monte-Carlo tree search[END_REF] to mention a few. We explain how the results of the previous sections can be applied when the stage cost ℓ is obtained by saturating a stage cost ℓ ′ which verifies the assumptions of Section IV.

We thus consider an a priori unbounded stage cost ℓ ′ : W → R and a state measure σ ′ : R n → R ≥0 . We denote the undiscounted cost function with stage cost ℓ ′ as J ′ . We focus on the scenario where Assumption 1-2 hold for system (1) with cost function J ′ . In addition, we assume that the properties below are verified.

Assumption 3:

(i) There exist M, λ > 0 such that for any x ∈ R n there exists an infinite-length control input sequence u ′ such that ℓ ′ (φ(k, x, u ′ | k ), u ′ k ) ≤ M σ ′ (x)e -λk for any k ∈ Z >0 . (ii) For any µ > 0 there exists k µ ∈ Z ≥0 such that for any

x ∈ R n there exists an infinite-length control input sequence u µ such that ℓ ′ (φ(k, x,

u µ | k ), u µ k ) ≤ µ for any k ∈ Z ≥0 with k ≥ k µ .
Item (i) of Assumption 3 corresponds to the exponential controllability property considered in Lemma 1. Item (b) of Proposition 2 is an uniform finite-time controllability property of ℓ ′ . It is likely possible to relax this condition to hold semiglobally, i.e. to make k µ depend on the ball of initial conditions; this extension is left for future work.

We now saturate ℓ ′ by an arbitrary constant ν > 0. In that way, the stage cost becomes ℓ(x, u) := min{ℓ ′ (x, u), ν} ∀(x, u) ∈ W. [START_REF] Sokolov | Complete stability analysis of a heuristic ADP control design[END_REF] In this case, we no longer consider the state measure σ ′ but σ which is defined as

σ(x) := min{ρ ′ (σ ′ (x)), ν} ∀x ∈ R n , (23)
where ρ ′ ∈ K ∞ is such that ρ ′ (σ ′ (x)) ≤ ℓ ′ (x, u) for any (x, u) ∈ W (such a function exists since Assumption 2 is assumed to hold when the cost function is given by J ′ ). We denote the undiscounted cost function with stage cost ℓ as J. The proposition below states that Assumptions 1-2 hold for system (1) with cost function J. Hence, the results of Section IV apply in this case with the state measure σ in [START_REF] Sorger | Local stability of stationary states in discounted optimal control systems[END_REF]. The proof is omitted due to space limitations. Proposition 2: Suppose the following holds.

(a) Assumptions 1-2 are verified for system (1) with cost function J ′ .

(b) Assumption 3 is verified. Then Assumptions 1-2 are guaranteed for system (1) with the cost function J and the conclusions of Theorems 1-3 hold with σ defined in [START_REF] Sorger | Local stability of stationary states in discounted optimal control systems[END_REF].

In view of Proposition 2, the state measure σ in ( 23) is ensured to semiglobally practically converge to the origin in the discounted case (since the conclusions of Theorem 3 apply). As a consequence, so does ρ ′ (σ ′ ) since σ becomes equal to ρ ′ (σ ′ ) as it approaches the origin, when δ is selected smaller than ν in Theorem 3. We therefore see that the semiglobal practical convergence of σ ′ , equivalently of ρ ′ (σ ′ ) as ρ ′ ∈ K ∞ , is preserved when saturating the stage cost function ℓ ′ . That means for instance that, when σ ′ (•) = | • |, the system solution converges to a neighborhood of the origin, which can be rendered as small as desired by decreasing δ (i.e. by increasing γ), when the stage cost is saturated as in [START_REF] Sokolov | Complete stability analysis of a heuristic ADP control design[END_REF]. The main difference appears for the stability. Indeed, to saturate ℓ ′ does not allow us to bound σ ′ (equivalently ρ ′ (σ ′ )) for all time by a class-K ∞ function of its initial value. However, we know that σ ′ is bounded as it converges to a neighborhood of the origin as time grows, which may be sufficient in many problems.

VII. CONTINUITY OF THE VALUE FUNCTION V

The purpose of this section is to provide conditions under which V is ensured to be continuous. In addition, we will guarantee the existence of a stationary, state-feedback optimal control law.

Assumption 4: The stage cost function ℓ is continuous on W and bounded, i.e. there exists ν > 0 such that ℓ(x, u) ≤ ν for any (x, u) ∈ W. Furthermore, the set-valued mapping x → U(x) is continuous and locally bounded. We are ready to state the main result of this section.

Proposition 3: Consider system (1), the cost function J in (2) with any γ ∈ (0, 1), and suppose Assumptions 1 and 4 are satisfied. Then the following holds.

(a) The optimal value function V is continuous on R n . (b) There exists an optimal solution which is given by a stationary, state-feedback law, i.e. there exists κ : R n → R m such that any x ∈ R n , κ(x) ∈ U(x) and the optimal sequence u * (x) = (u * 0 , u * 1 , . . . ) can be constructed as follows:

u * k = κ(φ(k, x, u * (x)| k )), for any k ∈ Z ≥0 .
Proof. We begin by observing that the boundedness of the stage cost from Assumption 4, its nonnegativity from item (i) of Assumption 1, and the existence of an optimal solution from Proposition 1, ensure the satisfaction of Assumption C of [START_REF] Bertsekas | Stochastic Optimal Control: The Discrete Time Case[END_REF] (the contraction assumption) by Proposition 4.10 of [START_REF] Bertsekas | Stochastic Optimal Control: The Discrete Time Case[END_REF]. This assumption in turn allows us to apply Propositions 4.1, 4.2, and 4.3 of [START_REF] Bertsekas | Stochastic Optimal Control: The Discrete Time Case[END_REF], which are essential to this proof.

We define the Bellman mapping T : B → B, where [T (v)](x) := min u∈U (x) [ℓ(x, u) + γv(f (x, u))]. Item (c) of Proposition 4.1 in [START_REF] Bertsekas | Stochastic Optimal Control: The Discrete Time Case[END_REF] shows that T is a contraction of modulus γ, i.e. that T (v) -T (v ′ ) ≤ γ vv ′ for any v, v ′ ∈ B, recalling also that • denotes the infinity norm. Due to item (a) of Proposition 4.2 in [START_REF] Bertsekas | Stochastic Optimal Control: The Discrete Time Case[END_REF], we obtain that the optimal value function V is the unique fixed point of this mapping, i.e. V = T (V ). By item (c) of Proposition 4.2 in [START_REF] Bertsekas | Stochastic Optimal Control: The Discrete Time Case[END_REF], the repeated application of T converges to V , that is, the sequence v 0 = 0, v n+1 = T (v n ) satisfies lim n→∞ Vv n = 0, with n ∈ Z ≥0 .

Using these facts, the following standard result also holds, by an immediate induction Vv n ≤ γ n Vv 0 = γ n V , so that the convergence of v n , n ∈ Z ≥0 , to V is uniform. We now prove by induction that v n is continuous for any n ∈ Z ≥0 . This assertion is true for n = 0 since v 0 = 0, suppose then that it holds for some n ∈ Z ≥0 . We have v n+1 : x → min u∈U (x)

[ℓ(x, u) + γv n (f (x, u))]. Since U is continuous, locally bounded (in view of Assumption 4), dom U = R n (by assumption), and ℓ + γv n is continuous, we apply Lemma 2 given in the Appendix to conclude that v n+1 is continuous. Consequently, v n is continuous for any n ∈ Z ≥0 . Hence, we have a sequence of continuous functions that converges uniformly to V , which allows us to reach our final conclusion that V is continuous, according to the uniform limit theorem.

VIII. CONCLUSIONS

We have shown that the solutions of a discounted optimal control problem preserves in a certain sense the stability properties of the undiscounted solution, by adjusting the discount factor. Further, the stability guarantees remain valid for a type of bounded stage costs, and the Lyapunov value function can be shown to be continuous when the stage costs are bounded, under mild assumptions.

An immediate line of future work is relaxing the detectability condition on the stage cost in (6) along the lines in [START_REF] Grimm | Model predictive control: for want of a local control Lyapunov function, all is not lost[END_REF], and the finite-time controllability condition as explained after Assumption 3. A second direction is motivated by the fact that practical algorithms rarely achieve exactly the optimal solution, ensuring instead only some bounded nearoptimality. Further, methods originating in artificial intelligence often require the inputs to be discretized, see [START_REF] Sutton | Reinforcement Learning: An Introduction[END_REF], [START_REF] Szepesvári | Algorithms for Reinforcement Learning[END_REF], [START_REF] Bus ¸oniu | Reinforcement Learning and Dynamic Programming Using Function Approximators[END_REF]. It would therefore be interesting to take into account these sources of errors, aiming to ensure stability for a wide class of practical algorithms.

  and, for each s ∈ R >0 , χ(s, •) is decreasing to zero. We denote by B the Banach space of bounded real functions on R n with the infinity norm v = sup x∈R n |v(x)|.

is the set of strictly increasing functions from Z ≥0 to Z ≥0 . When lim sup x→x S(x) = lim inf x→x S(x), we write the limit as lim x→x

  : ∃{x n } n∈Z ≥0 ∃{u n } n∈Z ≥0 s.t. x n →

			S(x) :=
			x→x
	with u n ∈ S(x n )}	and	n→∞ lim inf x→x S(x) x, u n → n→∞ u :=
	u : ∀{x n } n∈Z ≥0 s.t. x n → n→∞	x, ∃φ ∈ N ∃{u φ(n) } n∈Z ≥0
	s.t. u φ(n) →

u n→∞ u with u n ∈ S(x n ) , where N

Their work was partially supported by the ANR under the grant COMPACS (ANR-13-BS03-0004-02) and the European 7th Framework Network of Excellence 'Highly-complex and networked control systems' (HYCON2) (grant agreement no. 257462).

L. Bus ¸oniu is with the Department of Automation, Technical University of Cluj-Napoca, Romania, lucian.busoniuutcluj.ro. His work was supported by a grant of the Romanian National Authority for Scientific Research, CNCS-UEFISCDI, project number PNII-RU-TE-2012-3-0040, and by an internal grant of the Technical University of Cluj-Napoca. D. Nešić is with the Department of Electrical and His work was supported by the Australian Research Council under the Discovery Projects and Future Fellowship schemes.

APPENDIX

The following lemma is needed to prove Proposition 3. Its proof is omitted due to space constraints.

Lemma 2: Let ϕ : R n+m → R ≥0 be continuous and the set-valued mapping U : dom U ⇉ R m be continuous and locally bounded and dom U = R n . Then W : x → inf u∈U (x) ϕ(x, u) is continuous.


David Coulette
email: david.coulette@ipcms.unistra.fr

Giovanni Manfredi
email: giovanni.manfredi@ipcms.unistra.fr

An Eulerian Vlasov code for plasma-wall interactions

In this work we present a new Eulerian kinetic 1D3V code devoted to the study of plasma-wall interactions and sheath structure. The code solves the Vlasov-Poisson system for two or more kinetic species on a one-dimensional spatial grid between two limiting plates, in the presence of a uniform magnetic field tilted with respect to the normal to the plates. A source-sink term in the Vlasov equation leads to a stationary solution comprising a magnetic and a collisional presheath, as well as the usual Debye sheath in front of the wall. Thanks to a nonuniform spatial grid the sheath structure can be resolved accurately. Several advection schemes are implemented and the code is parallelized. Here, the code is used to illustrate the problems of a stationary sheath-presheath structure in front of a negatively biased wall, with and without a tilted magnetic field.

Introduction

Despite the constant improvement of power and accessibility of computing resources, the accurate numerical simulation of plasma-wall interactions remains a challenging task. The intermixing of various time and space scales, dynamical regimes (from collisional to ballistic) and mutually-interacting multiple species, along with their nontrivial interactions with materials, lead to highly complex and costly modeling. As of today, the most comprehensive simulations are massively parallel PIC/Monte Carlo codes [1]. Such simulations have the advantage that atomic processes can be modeled in a straightforward way, but require heavy computational resources, since large numbers of macro-particles must be used to obtain reasonable accuracy. Accuracy is particularly hard to maintain in low-density regions such as the Debye sheaths. In contrast, Eulerian (Vlasov, Fokker-Planck) codes have the potential to fare better in those situations, as the distribution function is sampled with the same accuracy in all phase-space regions. However, they are computationally more demanding, and their treatment of collisions is more costly and less intuitive compared to PIC codes.

In this paper, we report on the latest advancements in the development of an Eulerian code for the study of the plasma-wall transition. We adopt a simplified geometry restricted to one dimension in space, with a uniform and constant magnetic field, while retaining full 3D velocity space description. The current implementation is a parallelized extension and generalization of the codes used in [2,3,4] to arbitrary number of species. As an illustration, we will show a few numerical results on the stationary sheath-presheath structure in front of a negatively biased wall, with and without a tilted magnetic field.

Model and code structure 2.1. Physical model

We consider a plasma bounded by two material plates (Fig. 1). The system is assumed to be invariant along the two coordinates parallel to the plates. The x direction, normal to the plates, may either represent the distance along the principal axis of a small linear machine, or the curvilinear coordinate along an open magnetic field line between divertor plates in a tokamak. The system is magnetized by a constant and uniform magnetic field B 0 = B 0 (sin αe x + cos αe y ). For each species, the dynamics is governed by the collisional Vlasov equation

∂ t f s + v x ∂ x f s + q s m s (E + v × B 0 ) • ∇ v f = C s + S s , (1)
where f s (x, v) is the s th species distribution function, C s stands for the sum of all contributions from collisional processes, S s is the sum of all independent sources, and E = -∂ x φ is the electrostatic field obtained self-consistently from the Poisson equation The distribution functions solution of Eq. ( 1) are assumed to be zero at all boundaries in velocity space (±v max ). In the x direction, it is assumed that the incoming flux at the plates is zero (perfectly absorbing plate), but more complex boundary conditions can also be implemented to simulate particle-wall interactions such as physical sputtering or secondary electron emission. In some cases, only a half-infinite plasma will be considered, between one plate and a bulk position where the plasma is supposed to be at equilibrium.

-∂ xx φ = s q s 0 f s d 3 v. (2)
Boundary conditions for the electrostatic potential at the plates may either be mixed Neumann/Dirichlet or Dirichlet/Dirichlet, which allows for the treatment of both floating and externally biased walls. When a full kinetic description is not necessary, the dynamics of some species (most often electrons) can follow a simplified model such as a Boltzmann law

n s (x) = n s0 exp [-(q s (φ -φ ref )/k B T s ],

Numerical methods and implementation

Equation ( 1) is solved using a the following time-splitting scheme:

1 :

∂ t f s + v x ∂ x f s = 0 x advec. (∆t) 2 :
-∆φ -1 0 ρ = 0 Poisson solver 3 :

∂ t f s + qs ms (E x -B y v z )∂ vx f s = 0 v x advec(∆t) 4 : ∂ t f s + qs ms (B x v z )∂ vy f s = 0 v y advec(∆t) 5 : ∂ t f s + qs ms (B x v y -B y v x )∂ vz f s = 0 v z advec.(∆t) 6 : ∂ t f s + ν(f s -F 0 ) = 0 BGK relaxation (∆t) (3)
where the advections in all four phase-space directions (steps 1, 3-5) and the source and collision terms (step 6) are treated sequentially. Collisional processes are modeled using a single BGK relaxation operator acting as a restoring force that drives the distribution function towards a fixed Maxwellian

F 0 = n 0s exp(-m s v 2 /2T 0s ). Thus C s = -ν(f s -F 0 )
, where ν is the relaxation rate.

All of the conservative advection equations take the simple linear form

∂ t f + u ξ ∂ ξ f = 0, where ξ ∈ {x, v x , v y , v z }
and u ξ does not depend on ξ. They are integrated in time using a finite volume scheme [START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF] with reconstruction by primitive. Earlier versions of the code [2,3,4] used the flux-balanced scheme [START_REF] Fijalkow | [END_REF]. Several higher-order schemes have now been implemented and tested: (i) a fixed stencil polynomial interpolation with a choice of classical limiter functions, (ii) an ENO [7] scheme of arbitrary order, (iii) and a version of the PFC3 scheme [8] whose limiters have been modified to be used on a nonuniform grid. The latter is the default scheme in all results presented hereafter. It was selected as a good compromise between accuracy and computing effort.

The parallelization strategy relies on domain decomposition in the x dimension, implemented using MPI, and parallelization over species on each spatial sub-domain using OpenMP. This design does not fully exploit the parallelism that could be extracted by domain-decomposing the velocity dimensions. The rationale for this choice is to allow for the possibility of implementing more realistic collision operators. In that context, the full non-local coupling between velocity spaces of interacting species could lead to detrimental communication costs. Parallelization over velocity dimensions, either by domain decomposition or fine-grained OpenMP parallelism, will thus be considered only once the collision operators framework is implemented.

Scalability tests were performed to evaluate the MPI parallel performance of the current design. The first (Fig. 2a) measures the so-called "weak" scaling: starting from a given spatial sub-domain size n x = 64 with one MPI process N p = 1, global domain size is constrained to grow linearly with the number of processes (n x = 64N p ). The efficiency e = t(1)/t(N p ) is measured against the ideal case, for which the runtime would be constant. The results show an acceptable scaling: the communication overhead is nearly constant in the tested N p range. From a user point of view, it means that global domain size can be extended with no penalty.

The second test measures the "strong" scaling of efficiency. The global problem size is kept constant, while the number of sub-domains (and MPI processes) grows. The efficiency e = t(1)/[N p × t(N p )] is measured against the ideal case where runtime would scale as 1/N p . The results (Fig. 2b) show that the efficiency degrades and saturates with growing number of processes. The cost the Poisson solver [step 2 in Eq. ( 3)], which is parallelized using a block-decomposition scheme [9], stays negligible. Scalability is thus mostly limited by the x advection (step 1), for which the overhead brought forth by buffer exchange between subdomains, redundant computation on overlapping regions, and all-to-all communications required by the PFC3 scheme become non-negligible with respect to the ideally parallel steps 3-6. From a user's perspective, and for "small" spatial domain sizes (typically a few hundred points), one can expect a maximal speed-up of about 10 to 20 without wasting computing resources.

Numerical results

We tested our code on two well-documented steady-state problems: the plasma-wall transition with and without a tilted magnetic field. In both cases, we considered a half-infinite plasma between an absorbing plate (located at x = 0) and a bulk plasma at equilibrium (x = L). The domain length is chosen to encompass several mean-free paths relatively to the BGK relaxation operator, so as to provide a smooth transition in the pre-sheath. At t = 0, the ion distribution is set to a spatially homogeneous Maxwellian with bulk parameters. During the evolution, governed by the Vlasov-Poisson equation, the ions are lost to the wall and the sheath forms rapidly. The system then evolves slowly until the balance between particle transport and the BGK relaxation is established in the presheath. The results shown in the forthcoming sections are all taken in this asymptotic steady-state regime.

Static sheath in a linear machine

First, we consider the steady state of a small linear machine filled with a non-magnetized argon plasma, leading to a simple 1D1V model. An analytical kinetic treatment of similar systems in a collisionless regime can be found in Refs. [11,12]. The bulk plasma parameters (n e = n i = 5.5 × 10 15 m -3 , k B T i = 0.027eV , k B T e = 2.5eV ≈ 100k B T i ) were chosen to match those of Ref. [10]. In the actual machine, the wall potential is not driven by the bulk electrons, but by the highly energetic primary electrons emitted by the heated filament and then accelerated by a fixed discharge potential (-50V to -100V). The wall bias then reaches values of the order of the discharge potential. As the density of primary electrons is three orders of magnitude lower than that of bulk electrons, their direct contribution to the space charge density is neglected. In order to mimic their impact on wall bias, we force the wall potential to values of the order of the discharge potential. This assumption makes it possible to reproduce the ion velocity profile observed in the experiment [10].

The bulk electrons are treated using a Boltzmann model. Additional simulations with kinetic electrons showed good agreement with the results obtained using the Boltzmann law. All collisional processes affecting the ions are modeled using a single BGK operator that forces the ion velocity distribution function (IVDF) to relax towards a Maxwellian with bulk plasma parameters, with a relaxation rate ν ranging from ν = 10 -4 ω pi to ν = 10 -2 ω pi . Although the model is simple, the actual parameters impose quite stringent requirements on the phase space resolution if one wants to describe accurately the evolution of the ion distribution function up to the wall. The required velocity range is governed by the wall potential, which is rather large compared to the ion bulk temperature. In the sheath, the ion distribution becomes very narrow, which requires even higher resolution, as both cell averaging and numerical diffusion tend to smear out the distribution profile. For this set of test cases, a (n x , n vx ) = (1024, 1024) grid has been used. The cell width in x ranges from ∆x = 0.01λ D near the wall at x = 0 to ∆x = 100λ D at the right boundary (x ≈ 2 × 10 5 λ D ) where the plasma is at equilibrium

The density and temperature profiles are shown in Fig. 3. In the quasi-neutral pre-sheath, the ion behavior is dominated by the balance between plasma flow towards the wall and relaxation driven by the BGK term. The electric field tends to accelerate the ions towards the wall, whereas the BGK term tries to drag the ion distribution towards the equilibrium Maxwellian. The net result is that the ion velocity distribution widens in the presheath, leading to an increased kinetic temperature (Fig. 3b), defined as the velocity spread of the distribution function. The widening of the ion distribution in the presheath is clearly visible in Fig. 4a at position x = 385λ D .

In the Debye sheath, characterized by strong charge separation (Fig. 3a) and strong electric fields (Fig. 4b), the effect of the BGK term is negligible, so that the ion trajectories are essentially ballistic. The ion distribution reduces to a cold beam that is strongly accelerated towards the wall, as can be seen in Fig 4a between x = 18λ D and x = 0. The thinning of the ion distribution is reflected in a drop of the ion kinetic temperature (Fig. 3b,inset).

In order to check that the chosen grid and interpolation scheme do not introduce any artificial widening of the distribution function, either by cell averaging or numerical diffusion, the simulation results are compared against an analytical ballistic model in the sheath (Fig. 5). Starting from a reference distribution function at the entrance of the sheath (x = 34.3λ D ), the ion distribution functions closer to the wall are obtained by conservation of the distribution along the characteristics. Although a slight erosion of the IVDF maximum can be seen (mostly due to the action of the limiters), the IVDF geometry, and most importantly its asymmetry, are well preserved by the PFC3 scheme, and spurious oscillations are controlled.

v x c -1 s 0.0 0.2 0.4 0.6 0.8 1.0 1.2 IVDF (a.u) 0λ D 5λ D 10λ D 18λ D 385λ D 45563λ D (a)
0 10 20 30 40

xλ -1 D 2.0 1.5 1.0 0.5 0.0 E x (λ D e)(k B T e ) -1 (b)

Magnetized deuterium plasma with tilted magnetic field

In order to test the operation of the code with a full 3D description in velocity space, we consider a deuterium plasma with kinetic ions and Boltzmann electrons. The setup is similar to the cases studied in Refs. [3,4,13,14], with the geometry depicted in Figure 1. The discharge parameters are T e = T i = 5 eV, ω Ci = 0.05 ω pi , α = 15 • , ν = 0.002 ω pi . The phase space is discretized on a (256, 96, 96, 96) grid. The cell width ranges from ∆x = 0.5λ D to ∆x = 500λ D , and the ion velocity range is [-6v thi ; 6v thi ] in all directions. Albeit rather coarse, particularly in the sheath, this resolution is sufficient to obtain the main features of the collisional pre-sheath/magnetic pre-sheath/sheath transitions described in Ref. [13].

In a non-magnetized (or with normal incidence, α = 90 • ) and weakly collisional plasma, the transition from the equilibrium bulk plasma to the quasi-ballistic regime in the sheath occurs in the collisional pre-sheath, where collisional/ionization processes compete with the ion flow towards the wall. In the magnetized plasma considered here, the magnetic field prevents the ions to flow directly towards the wall, particularly for grazing incidence. The overall transition from bulk to wall can be seen by looking at the evolution of the steadystate ion distribution function in the (v x , v y ) plane (Fig. 6). In the collisional pre-sheath, as the ion mean velocity is slowly accelerated up to the sound speed, the magnetic field maintains the flow aligned with its own direction (Fig. 6a-c). In the magnetic pre-sheath, while the plasma is still quasi-neutral (Fig. 7b), the magnitude of the electric field is sufficient to break the alignment of the flow with the magnetic field (Fig. 6d). Eventually, the electric field dominates in the Debye sheath (Fig. 6e-f), where the ions are strongly accelerated in the direction normal to the wall.

The transition from the B-dominated pre-sheath to the E-dominated sheath is clearly visible in Fig. 7a, where we compare the magnitude of the various forces acting in the x direction. The Debye sheath entrance corresponds not only to the beginning of the nonneutral region (Fig. 7b), but also to the point where the electric force starts dominating over the Lorentz force v × B.

Conclusion and perspectives

We presented the first developments of an Eulerian 1D3V code devoted to the study of plasmawall interactions. Building up from earlier versions [3,4], significant modifications were made in order to adapt the code to its planned applications. Particle fluxes reaching the wall are the source of complex dynamics, including non-trivial feedback mechanisms, so that an accurate description of the distribution function must be preserved up to the wall. This necessity triggered the switch from the flux balance method to the PFC3 advection scheme. In addition, the present version can also address dynamical transient, and not only steady-state solutions, as was the case before [3,4]. Finally, the parallelization over spatial sub-domains and different species opens the way to more demanding applications with acceptable run times. Some standard cases, both in 1D1V and 1D3V, were used to check the validity and parallel performance of the code. The next steps will concentrate on the fully kinetic regime, in particular the comparison of Boltzmann and kinetic electron models for near-grazing incidence (α = 1-5 • ) of the magnetic field. It is also planned to incorporate more accurate collision models than the simple BGK relaxation used so far. This is a major step, as it will introduce nontrivia integral terms in velocity space. Finally, more general plasma-wall effects should also be modeled, such as secondary-electron emission and sputtering from the wall [4].

F ext x × eλ D k B T e -q i E x q i (v ×B) x -ν m i v x (

Figure 1 .

 1 Figure 1. Model geometry (a) and non-uniform grid near the wall (b)

Figure 3 .

 3 Figure 3. Argon discharge steady-state with 50V wall bias and collision rate ν = 10 -4 ω pi : particles densities (a); ion temperature profile (b). The insets show zooms of the respective quantities.

Figure 4 .

 4 Figure 4. Argon discharge with -50V wall bias and collision rate ν = 10 -4 ω pi . Evolution of the IVDF (normalized to its peak value) from the bulk plasma to the wall (a); distances from the wall are indicated above each peak. Electric field profile near the wall (b).

Figure 5 .

 5 Figure 5. Argon discharge with a -50V wall bias and ν = 10 -4 ω pi . Evolution of the IVDF in the sheath: numerical results (continuous lines); IVDF computed from the reference profile at x = 34.3λ D using a ballistic model (crosses). The vertical lines mark the ballistic model peak velocity. All distributions are normalized to the peak value of the reference distribution at x = 34.3λ D .

Figure 6 .

 6 Figure 6. Contour plots of the v z -averaged ion distributions normalized to their peak value (indicated in the upper left corner); dashed lines indicate the B 0 and e z × B 0 directions. Bulk plasma (a); collisional pre-sheath (b-c); magnetic pre-sheath (d); Debye sheath (e); wall (f). The orange curve in (f) retraces the evolution of the peak of the distribution from the bulk plasma to the wall.

Figure 7 .

 7 Figure 7. Deuterium plasma with tilted magnetic field, α = 15 • . Profiles of the various forces in the momentum balance equation (a). Note that the friction force -νm i v x caused by the BGK term is negligible over the whole domain. Density profiles for electrons and ions (b).

  in which cases the Poisson equation becomes nonlinear and is solved iteratively. Parallel efficiency of x domain decomposition (MPI): N p is the number of MPI processes, t(N p ) is the run duration. Efficiency t(1)/t(N p ) scaling with growing global domain size n x = 64 × N p and (n vx , n vy , n vz ) = (64, 64, 64) (a) and efficiency t(1)/[N p × t(N p )] scaling with a constant-size global domain n x = 1024, (n vx , n vy , n vz ) = (96, 96, 96) (b).

		1.0					1.0
		0.8					0.8
	t(1)=t(N p )	0.4 0.6					0.4 0.6 £t(N p )) t(1)=(N p
		0.2					0.2
		0.0	16	32	48	64	1 0.0	8	16	32
				N p					N p
			(a)					(b)
	Figure 2.

Acknowledgements

The authors acknowledge the support of the French Agence Nationale de la Recherche (ANR) under reference ANR-12-BS09-028-01.


Charles Maragna
email: c.maragna@brgm.fr

Xavier Rachez

Innovative Methodology to Compute the Temperature Evolution of Pile Heat Exchangers

Keywords:

Energy geostructures such as heat exchanger piles couple the structural role of geostructures with heat and cold supply via shallow geothermal energy. This combination makes it possible to cut down the investment costs of ground heat exchangers (GHE). Thermal dynamic simulations require numerical models of pile heat exchangers to run over a reasonable amount of time. In this perspective semi-analytical models seem interesting. The paper presents an approach to semi-analytical modeling of pile heat exchangers. This approach relies on three elements: First, correlations are established to describe the evolution of wall pipe temperature under constant heat load in presence of underground water flow. These correlations take into account the underground water flow in the vicinity of the pile, allowing the computation of the temperature evolution over both the short and long terms. Second, a resistive-capacitive (RC) circuit is developed to account for the thermal inertia of the pile concrete without having to mesh its geometry. Third, the RC circuit is combined to a heat balance over the heat carrier fluid and to correlations to compute temperature evolution in the fluid. Predictions of this semi-analytical model are compared with those of a fully-discretized finite elements model. A good agreement between both models is reached. Acquisition of in-situ data is foreseen to validate both models.

INTRODUCTION

Geothermal energy is becoming widely used to supply buildings with renewable heating and cooling. This development has been taking place in a context of policies fostering the reduction of green-houses gas emission. For instance the total number of geothermal heat pumps installed in Europe has been recently estimated to reach 1,140,000 and to represent an equivalent thermal power of 14,000 MWth (Observ"ER, 2013). In this context, a diversification of the heat sources and sinks of ground-sourced heat pumps has been researched. Energy geostructures such as piles heat exchanger (Figure 1) appear to be an attractive technology as costs for both geostructures and ground heat exchanger (GHE) realization are mutualized which allows cutting down the costs of GHE achievement. The dual role of the energy geostructures (mechanical and energetics) makes their design more complex and challenging [START_REF] Laloui | Energy Geostructures -Innovation in Underground Engineering[END_REF]. The thermo-mechanical behaviour of energy geostructures has been and is still extensively studied [START_REF] Suryatriyastuti | Understanding the temperature-induced mechanical behaviour of energy pile foundations[END_REF][START_REF] Di Donna | Thermo-mechanical aspects of energy piles[END_REF], and helps developing thermo-mechanical design tools of energy piles. The thermal performance of energy piles has been less studied in details and still needs improvement. For instance, the energy pile thermal design may be inadequate if the same design approaches as for borehole heat exchangers is used for piles [START_REF] Bourne-Webb | Observed Responses of Energy Geostructures[END_REF].

A thermal dynamic simulation is often necessary to optimize the energetic system made of ground heat exchangers, heat pump and building. To do so it appears mandatory to rely on numerical models of pile heat exchangers (PHE) that run over a reasonable amount of time. The resolution of the heat and mass balance equations with numerical techniques such as the finite element (FE) often leads to computation times that are not compatible with engineering practices. Analytical solutions are an alternative.

However the range of these solutions is often limited due to the fact that no known analytical solution describes both the cylindrical GHE geometry and the surrounding underground water flow correctly [START_REF] Wagner | Analytical approach to groundwater-influenced thermal response tests of grouted borehole heat exchangers[END_REF].

The paper presents a new approach to semi-analytical thermal modeling of pile heat exchangers. This approach relies on three innovative elements. After a short presentation of assumption (section 2.1.), the paper first presents correlations that have been established over the temperature evolution of the wall pipe under constant heat load, referred to as step response (section 2.2). These step responses have been designed to correctly describe the temperature evolution in the vicinity of the pile. Second, a resistive-capacitive (RC) circuit is developed to account for the thermal inertia of the pile concrete (section 2.3). Third, the RC circuit is combined to a heat balance over the heat carrier fluid and to correlations for step response to obtain a semi-analytical model. A numerical scheme is used to solve temperatures in the fluid, pile and surrounding ground: Temperatures at the nodes of the RC circuit and in fluid are solved with an implicit Euler method, while the evolution of pile wall temperature is computed by discretizing the convolution of time-variable heat flux with the step response (section 2.4). Finally, predictions of the semianalytical model model are compared with those of a fully-discretized finite elements model (section 3).

DESCRIPTION OF THE INNOVATIVE METHODOLOGY 2.1. Assumptions

Throughout the whole paper, the following assumptions are made:

(i) The ground is regarded as a homogenous porous media saturated with water where Darcy"s law applies. (ii) At the scale of observation, the concrete material properties are supposed to be homogeneous. As the hydraulic conductivity of concrete for piles would typically be of the order of 10 -12 m/s, the concrete is considered as an impermeable material. (iii) Strains originating from thermal stresses are not taken into account in this paper. Both pile concrete and soil matrix are therefore considered to be non-deformable. (iv) Physical properties of the materials (underground water, soil matrix, PHE heat carrier fluid) do not depend upon temperature. (v) The initial, non-disturbed temperature T 0 is constant in the whole domain; far away from the pile heat exchanger temperature remains constant. In what follows we simply set T 0 =0 K. (vi) Underground water flows and heat transfer occurs in a horizontal plan. Underground water flow is determined in stationary regime.

Moving infinite cylindrical source model. New approach

This section aims at developing a "step response" function that describes the temperature evolution at the pile wall under a constant linear heat flux (i.e. per borehole length unit, in W.m -1 ) in presence of underground water flow (cf. Figure 2): (1) denotes the thermal conductivity of the porous media (W.K -1 .m -1 ).

accounts for moving infinite cylindrical source model. Water flow in the soil matrix is driven by hydraulic head gradient, according to Darcy"s law:

(2)

Where denotes the Darcy velocity (m.s -1 ), the hydraulic conductivity (m.s -1 ), the hydraulic head (m). Under the assumption of incompressible water, the equation for water mass conservation reads:

( ) (3)
Heat in a homogenous porous medium with ground water flow is transported by conduction and advection. The partial derivative equation for energy conservation reads:

( ) ( ) (4)
Where the subscripts and respectively denote the ground and the underground water. The volume-specific heat capacity ( ) (J.K.m -3 ) of the porous media is related to the volume-specific heat capacities of the matrix ( ) and of the water ( ) by:

( ) ( ) ( ) (5)
refers to the total porosity of the porous media.

As the set of partial derivative equations ( 3) and ( 4) are linear, they can be normalized. Two dimensionless numbers, Fourier Fo and Peclet Pe, are introduced. We took as reference length the radius of the pile . Hence a normalized length of is:

(6)

Fo is the Fourier number characterizing the ratio of diffused heat to stored heat:

(7)
is a normalized time. This involves that in dynamic simulation Fo replaces the time t.

is the thermal diffusivity of the ground (m².s -1 ) defined by:

( ) (8)
The Peclet number is a parameter accounting for the ratio of the heat transferred by convection to the heat transferred by conduction:

( ) (9)
Higher Peclet numbers lead to non-negligible convection compared to conduction. denotes the "undisturbed" Darcy flow, i.e. far away from the pile, where the impermeable pile doesn"t deflect underground water flow.

The normalized equation for heat transport reads:

(10) Normalizing the heat equation has a huge advantage: It makes it possible to greatly decrease the number of parameters involved in the problem. The temperature T at a point of coordinates (x, y) depends on five parameters (( ) , , ( ) , , ) and three variables (x, y, t) while the normalized temperature T* depends on a single parameter (Pe) and three variables (x*, y*, Fo).

The step response

is the normalized temperature averaged on the pile wall perimeter . Therefore, depends on one parameter (Pe) and one variable (Fo) :

∫ (11)
For practical engineering problems, the temperature evolution at the pile wall can be computed based if the step response is known thanks to Equation (1).

Carslaw and Jaeger [START_REF] Carslaw | Conduction of Heat in Solids[END_REF] analytically determined a step response in case heat is only transferred by conduction (no convection, i.e. Pe=0) where the cylindrical geometry of the borehole is accounted for. This step response G ICS is referred to as the infinite cylindrical source (ICS) model:

∫ ( ) (12)
Here refers to the classical function evaluated at the pile wall, i.e. for the normalized distance from the pile axis . The same authors provided another step response where conduction and convection are taken into account, but where the pile cross-section degenerates to a point. This step response is here denoted G MILS and is known as the moving infinite line source (MILS) model. It can be expressed as a function of Fo and Pe as:

( ) ( ( ) ) (13)
The MILS model replaces the cylindrical geometry of the pile by a simple line. Therefore, it does not correctly account for the cylindrical geometry of the pile. Consequently the temperature evolution on short duration is not described well. Moreover, Wagner et al. [START_REF] Wagner | Analytical approach to groundwater-influenced thermal response tests of grouted borehole heat exchangers[END_REF] reported that using the MILS model at high Peclet numbers may lead to underestimate the wall pile temperature evolution in the long term. They attributed this observation to the fact that the MILS model does not account for the fact that the pile is an obstacle to underground water flow, since the MILS model consider the pile cross-section as a point. In the MILS model, underground water flows "straight ahead"; in other words, the MILS model does not take into account the deflection of the underground water flow by the impermeable pile wall. Therefore, the MILS model leads to an overestimation of the contribution to heat transfer by convection.

It appeared extremely difficult to find an analytical solution to the coupled PDE ( 3) and ( 4) to determine the normalized temperature T* in the vicinity of the ground and in fine the step response . To get around the problem, the COMSOL-Multiphysics software was used to solve the PDEs for Fourier number ranging from Fo=0 to Fo=10 4 . COMSOL-Multiphysics is a versatile tool for the resolution of coupled partial derivative equations based on the finite elements method. The domain mesh and an example of computed normalized temperature field T* and streamlines in the borehole surrounding can be visualized in Figure 3. A parametric study for Pe values ranging from 0.05 to 1.5 was carried out. A maximum time t max corresponding to the final Fourier number Fo max =10 4 and the Darcy velocity can respectively be computed thanks to Equation ( 7) and ( 9):

(14) ( ) (15)
A maximum time t max corresponding to the final Fourier number Fo max =10 4 and the Darcy velocity corresponding to the higher Peclet number (Pe=1.5) are reported in Table 1 for typical values of pile heat exchangers radius and underground characteristics. The ground is considered made of sand saturated with water. Its thermal properties were taken from the recommended values given in the Swiss guideline for BHE sizing SIA 384/6 (SIA, 2010). Based upon these parameters, the computed MICS step responses are valid for Darcy velocities in the order of 100 m.y -1 and for durations up to several years. The range of investigated Fo and Pe seems to be adequate for most practical engineering applications.

Input data Output data

Pile radius rb (m) Volumespecific heat capacity of the underground water ( ) Figure 4 shows computed step responses for the three hereby mentioned models, for Pe=0.1 (conduction-dominated regime), Pe=0.5 (mixed regime) and Pe=1.5 (convection-dominated regime). It can be noticed that:

(MJ.K -1 .m -3 ) Volume- specific heat capacity of the ground ( ) (MJ.K -1 .m -3 ) Thermal conductivity of the ground (W.K -1 .m -
(i) MICS step responses fits ICS step responses well at the beginning of the solicitation (i.e. for low Fo number). This behavior was expected since both models account for the cylindrical geometry of the pile. Therefore, both ICS and MICS correctly model the temperature evolution in the vicinity of the pile, which is solicited at the beginning of the simulation.

(ii)

In the presence of underground water flow (i.e. Pe>0) the step response converges to a steady-state value. (iii)

At low Pe number (Pe=0.1), steady-state values of MILS and MICS step responses are close. (iv)

At higher Pe number (Pe=0.5 and Pe=1.5), the MICS model gives higher estimation of the steady-state step responses than the MILS model. It seems realistic as the MICS takes into account the fact that the impermeable pile is an obstacle to the underground water flow which deflects it. However the MILS model underground considers that water flows "straight ahead" and is not deflected by the pile geometry, since its geometry is degenerated to a point. Qualitatively, this agrees with the results provided by Wagner et al. [START_REF] Wagner | Analytical approach to groundwater-influenced thermal response tests of grouted borehole heat exchangers[END_REF]. Based on our results, the MILS model compared to the MICS model underestimates the steady-state step responses by 7 % for Pe=0.5 and by 30 % for Pe=1.5.

Development of a new resistive-capacitive circuit

A thermal dynamic simulation of a building, heat pump and ground heat exchangers requires a short time step (typically one hour or less) to capture the evolution in heat and cold demand. Therefore, an adequate description of the transient thermal behavior of the pile is needed. Basically, two approaches are possible: either the whole cross-section of the pipe is discretized and the heat equation is solved with numerical techniques such as finite elements (FE), or the temperature evolution in the is computed by analytical models. Fully discretized models often lead to extensive computation times, while most analytic models do not take into account the thermal inertia (capacity) of the concrete pile [START_REF] Bauer | Diersch: Thermal resistance and capacity models for borehole heat exchangers[END_REF].

Development of analytic resistive-capacitive (RC) models that could deal with the thermal inertia of the concrete pile have been reported (Bauer et al., 2011, Shirazi and[START_REF] Shirazi | Thermal capacity effects in borehole ground heat exchangers[END_REF]. In RC circuits, the temperature evolution in the concrete is accounted for by nodes. Nodes are given a capacity (J.K -1 .m -1 ) and are connected by resistances (K.m.W -1 ). Such existing models seem to be limited to piles equipped with two or four pipes. However, pile heat exchangers can be equipped with six or eight pipes, depending upon the pile radius and engineering constraints [START_REF] Laloui | Energy Geostructures -Innovation in Underground Engineering[END_REF]. Besides, it shall be noticed that RC models have been developed for borehole heat exchangers characterized by a typical borehole radius r b in the range 6-10 cm, while the pile heat exchangers investigated in our study are characterized by r b in the typical range 15-30 cm. As the thermal capacity of a pile section is proportional to its surface, it seems necessary to correctly describe the temperature evolution within the pile. Therefore, disposing on accurate RC models appears to be a problem of first importance to achieve accurate dynamic simulations of pile heat exchangers.

To overcome these obstacles we developed a RC circuit where the concrete is divided in as many zones as pipe in the pile. For the sake of clarity, we limit the paper to the case where the pile is equipped with four pipes (cf. Figure 5). The model can be easily extended to the case there are more pipes. Every outer face of a pipe, denoted (i=1,2,3,4) is connected to the opposite portion of the borehole wall by a branch made of three resistances , , . Two capacities and are inserted at the subsequent nodes. This portion of the circuit aims at describing the temperature evolution in the outer part of the pile, i.e. between the pipe and the borehole wall. The choice of the number of resistances and capacities results from a trade-off between an accurate description and a high complexity of the model. Every pipe is connected to the central part of the pile by a branch with two resistances and in the middle of which lies a capacity

. Finally, the central part of the pile is represented by a capacity while interactions between two close pipes (such as and ) are accounted for by a capacity lying between two resistances .

To compute the eleven parameters of the RC model we followed a procedure where the underlying idea is to solve the equation in a fully-discretized model and to "tune" the parameters of the RC model to fit the evolution in energy computed by the fullydiscretized model. First, the heat equation inside the pile is normalized, so that to decrease the number of relevant parameters, in a similar way to what we did in section 2.2. Second, a set of dynamic simulations characterized by specific boundary conditions are determined. These simulations are run for a FE model developed in COMSOL-Multiphysics; evolution of energy and fluxes are computed for each simulation. In parallel, a RC circuit computing the temperature evolution at the nodes has been developed in the MATLAB software. Finally, parameters of the RC model are determined by fitting the energy computed by the RC model to the energy computed by the FE model.

The Table 2 gives an example of computed values of the RC parameters for a PHE characterized by a radius of 0.3 m. Thermal properties of the pile concrete are similar to the data used by Suryatriyastuti [START_REF] Suryatriyastuti | Understanding the temperature-induced mechanical behaviour of energy pile foundations[END_REF] and Gashti [START_REF] Gashti | Numerical modelling of thermal regimes in steel energy pile foundations: A case study[END_REF].

Physical parameters of the pile (input data)

Thermal conductivity of the concrete λc

Volume-specific heat capacity of the concrete ( )

Borehole pile radius rb

Outer radius of the pipe rp Distance between two opposite pipes s

1.8 W.K -1 .m -1 2.11 MJ.K -1 .m -3 0.3 m 0.01 m 0.3 m
Computed RC parameters (output data)

Resistances (K.m.W -1 ) Capacities (kJ.K -1 .m -1 ) R 1 R 2,1 R 2,2 R 3,1 R 3,2 R 3,3
2.158 0.135 1.343 0.220 0.119 0.021 6.9 82.0 47.6 32.6 3.9

Table 2: Computed parameters of the RC circuit

Development of a semi-analytical model to compute the evolution of fluid temperature (SA model)

A numerical semi-analytical model (SA model) has been developed to compute the evolution of temperatures in the fluid, pile and surrounding ground. "Semi-analytical" refers to the fact that the temperature at the pile wall is computed based on a step response G. Therefore, no numerical resolution of the heat equation in the surrounding ground is carried out.

The RC circuit is coupled to a heat balance over the fluid (cf. Figure 6-a). The power exchanged by meter of pile (W.m -1 ) between the fluid and the surrounding ground at each time step is an input data of the model. A heat balance over the heat carrier fluid gives:

̇ (16)
denotes the depth of the PHE, ̇ the flow rate in the pipe (kg.s -1 ), the specific heat capacity of the heat carrier fluid (J.K -1 .kg), and the PHE inlet and outlet temperatures respectively. We introduce , a mean temperature of the fluid inside the pipes:

(17)

In the paper we limit to the assumption that the fluid temperature is the same in every pipe. In other words, no thermal "shortcircuit" between cold and hot heat-carrier fluid is taken into account. Due to symmetry only one quarter of the RC circuit is modeled. Therefore, at the node L (i.e. at the heat-carrier fluid) the applied power is defined by:

(18)

Where R p denotes the thermal resistance of the pipe. The temperature evolution can be expressed by a heat balance as:

[ ] { } [ ]{ } { } (19)
Where { } is a vector containing the temperature at the eight nodes:

{ } { } ( 20
)
[ ] is the capacity matrix (J.K -1 .m -1 ), [ ] the conductance matrix (W.K -1 .m -1 ), { } has the dimension of a power vector (W.m -1 ).

[ ], [ ], { } are filled by the heat balance at every node. The temperatures are computed at discretized times t n = n Δt, where Δt is the time step of the modification. The node B 3 deserves a particular attention as it is connected to the ground. To compute its temperature, we used the superposition principle [START_REF] Ingersoll | Theory of the ground pipe heat source for the heat pump[END_REF]:

( ∑ ( ) ) (21)
Where denotes the step response G evaluated at t k = k Δt. In what follows we used the MICS model, i.e.

. The Peclet number is set to Pe=0. In other words underground water doesn"t flow ( . is the total heat flux exchanged at the time step k between the pile and the surrounding ground. An algorithm was developed in MATLAB to solve the system of equations ( 19) with a basic implicit Euler method [START_REF] Mohammadi | Pratique de la simulation numérique[END_REF]. For both models, the sum of concrete power and ground power is constant and equal to the applied power (50 W.m -1 ). This constitutes a verification that both models conserve power. (ii)

Both models give really close evolution of the internal powers. This is an element of validation of the SA model. (iii)

At the beginning of the solicitation, is equal to , while is equal to 0, since transient thermal transfer occurs in the pile concrete and no power is transferred to the surrounding ground. As the time goes, pile concrete warms up and heat starts to be transferred to the surrounding ground through the pile wall; as a consequence increase and decreases. (iv) It takes 7.5 hours to reach =25.0 W.m -1 , i.e. for the concrete pile to be half loaded. After 20 hours =37 W.m -1 , this means that the concrete is loaded at 77 % (=38/50*100). It takes 51 hours for the concrete to be loaded at 95 % (out of the graph). These results suggest that heat transfer in pile heat exchanger is a transient phenomenon that lasts over dozens of hours. It highlights the fact that, depending on the level of precision foreseen, heat transfer in a pile should be treated with relevant models taking into account concrete thermal inertia. A resistive-capacitive (RC) circuit has been developed to deal with the thermal inertia of the concrete pile. The parameters of the RC circuit are "tuned" to fit the energy evolution computed by a fully discretized model. The RC circuit has been integrated in a semi-analytical model (SA model) and combined with the G MICS functions to compute the evolution of fluid temperature. The system of equation is solved with an implicit Euler scheme.

A comparison of the SA model, a fully discretized finite-elements (FE) model and a model that overlooks the thermal capacity of the pile (1-R model) has been carried out. Fluid temperatures computed by the SA and FE models are really close, which is an element of validation of the numerical implementation of the developed SA model. The 1-R model overestimates the fluid temperature compared to the FE and SA models, by 2.6 K after 1 hour, by 1.8 K after 10 hours. Besides, it takes 7.5 hours to half load the pile concrete and it takes 51 hours for the concrete to be loaded at 95 %.

The semi-analytical model seems to be an interesting tool for dynamic simulations of heat pumps connected to pile heat exchangers. In-situ tests achieved in the GECKO project will deliver experimental data to validate this model. Further developments include the computation of RC parameters for a representative number of pipe geometry and equipment, and the integration of the semi-analytical model into an environment of thermal dynamic simulation such as TRNSYS.

ACKNOWLEDGEMENT

The work described in this paper is a part of a research project ""GECKO"" (geo-structures and hybrid solar panel coupling for optimized energy storage) which is supported by a grant from the French National Research Agency (ANR). It is an industrial project, involving an engineering office, three public companies and two research laboratories in civil and energy engineering sector: ECOME, BRGM, IFSTTAR, CETE Nord Picardie, LGCgE-Polytech"Lille and LEMTA-INPL. The authors would like to express their gratitude to their partners for their support in this ongoing project.

NOMENCLATURE

T temperature (°C) Subscripts power per meter of pile (W.m -1 ) b borehole wall c concrete thermal conductivity (W.K -1 .m -1 ) ground Darcy velocity (m.s -1 ) w underground water density (kg.m -3 ) fl heat-carrier fluid specific heat capacity (J.K -1 .kg -1 ) in inlet volume-specific heat capacity (J.K -1 .m -3 ) out outlet thermal diffusivity (m².s -1 ) 0 undisturbed conditions Fourier number at/r b ² Peclet number thermal resistance (K.m.W -1 ) Superscripts capacity (J.K -1 .m

Figure 1 :

 1 Figure 1: Example of a pile heat exchanger (PHE) before installation in ground (with the courtesy of Auxerre city -France)

Figure 2 :

 2 Figure 2: Scheme of the modeled domain and boundary condition

Figure 3 :

 3 Figure 3: (a) Domain mesh. (b) Computed normalized temperature T* and streamlines in the borehole surrounding at normalized time (Fourier number) Fo = 1000. In this simulation the Peclet number Pe was set to 0.5. The water flows rightwards.

Figure 4 :

 4 Figure 4: Computed step responses for Pe = 0.1, 0.5, 1.5 as a function of the Fourier number Fo (normalized time). ICS: infinite cylindrical source; MILS: moving infinite line source; MICS: moving infinite cylindrical source (new approach).

Figure 5 :

 5 Figure 5: (a) cross-sectional view of the PHE geometry. (b) The developed RC model. The RC circuit that we developed has 11 parameters: (i) Six resistances: , , , , , . (ii) Five capacities: , , , , .

Figure 6 :

 6 Figure 6: (a) New semi-analytical model, referred to as "SA" model. (b) One-resistance model, referred to as "1-R" model.

Figure 7 :

 7 Figure 7: Evolution of mean temperature fluid T fl computed by the three models in normal scale (a) and in semi-logarithmic scale (b).

Figure 8

 8 Figure8respectively give the evolution of the internal power (W.m -1 ) of the concrete and of the ground computed by SA and FE models. It can be emphasized that:(i)For both models, the sum of concrete power and ground power is constant and equal to the applied power (50 W.m -1 ). This constitutes a verification that both models conserve power. (ii)Both models give really close evolution of the internal powers. This is an element of validation of the SA model. (iii)At the beginning of the solicitation, is equal to , while is equal to 0, since transient thermal transfer occurs in the pile concrete and no power is transferred to the surrounding ground. As the time goes, pile concrete warms up and heat starts to be transferred to the surrounding ground through the pile wall; as a consequence increase and decreases. (iv) It takes 7.5 hours to reach =25.0 W.m -1 , i.e. for the concrete pile to be half loaded. After 20 hours =37 W.m -1 , this means that the concrete is loaded at 77 % (=38/50*100). It takes 51 hours for the concrete to be loaded at 95 % (out of the graph). These results suggest that heat transfer in pile heat exchanger is a transient phenomenon that lasts over dozens of hours. It highlights the fact that, depending on the level of precision foreseen, heat transfer in a pile should be treated with relevant models taking into account concrete thermal inertia.

Figure 8 :

 8 Figure 8: Ground internal power and concrete internal power computed by the SA and FE models for the 40 first hours.4. CONCLUSIONA new MICS (moving infinite cylindrical source) model has been developed to compute the temperature evolution at the pile wall under steady underground water flow. These step responses G MICS functions are parameterized by the Peclet number Pe, the ratio of the heat transferred by convection to the heat transferred by conduction. As it accounts for the cylindrical geometry of the pile, the MICS model correctly describes the temperature evolution over the short term. It takes into account the fact the pile is an obstacle to the underground flow and therefore gives higher estimation of the long-term temperature that the MILS (moving infinite line source) model. Our results suggest that the MILS model underestimates the steady-state step responses by 7 % for Pe=0.5 and by 30 % for Pe=1.5 compared to the MICS model.

Table 1 : Typical values of pile heat exchangers radius and underground characteristics ("Input data") and corresponding maximum investigated time t max and Darcy velocity ("Output data") for and Fourier Fo=10 4 and Pe=1.5.

 1

	Fourier	Peclet number	Maximum time	Darcy velocity
	number Fo	Pe	tmax (h)	(m.y -1 )
	1 )

COMPARISON OF THE SEMI-ANALYTICAL MODEL WITH A FINITE-ELEMENT MODEL

A comparison of three models was carried out for parameters given in Table 3: (i)

A numerical model ("FE" model) where the PHE cross-section is fully discretized. The assumption of a homogenous mean fluid temperature remains. (ii)

The semi-analytical model ("SA" model) presented in section 2.4. (iii)

A one-resistance model, where the thermal capacity of the pile is neglected ("1-R" model, cf. Figure 6-b). This model takes into account the cylindrical geometry of the borehole. It lies upon a borehole resistivity R b given by: (

In each model the initial temperature is T 0 = 0 K. The power exchanged by meter of pile is set to 50 W.m -1 : The fluid cools down while the concrete and ground warm up. 50 W.m -1 is a typical value for thermal response tests (TRT).

PHE characteristics Ground characteristics

Thermal resistance of the pipe Rp

Pile depth H Thermal conductivity of the ground λm Volume-specific heat capacity of the ground ( ) Figure 7 shows the evolutions of the mean temperature fluid T fl for the three models. It can be emphasized that:

Heat-carrier fluid characteristics PHE solicitation

(i) All models converge to the same value at t=100 days. It was expected since once the transient thermal transfer within the pile is over the thermal transfer is driven by conduction in the ground. Each model takes into account the same conductivity λm=2.3 W.K -1 .m -1 and volume-specific heat capacity of the ground ( ) = 2.4 MJ.K -1 .m -3 .

(ii) FE and SA models give close results. After one hour, the mean fluid temperature inside the pipes T fl is respectively estimated at3.7 K and 3.4 K by the FE and SA models, i.e. a temperature difference of 0.3 K. After 10 hours, the difference reduces to 0.2 K since the models respectively returns 5.6 K and 5.8 K. The good agreement between the two models is an element of validation of the RC circuit, the computation procedure of the RC parameters as well as the development and implementation of the semi-analytical model. (iii) 1-R model greatly overestimate T fl compared to FE elements. After one hour, the 1-R model gives 6.3 K against 3.7 K for the FE model, which means an overestimation of almost 2.6 K. After 10 hours, the 1-R model returns 7.6 K against 5.8 K; the overestimation is 1.8 K. It should be kept in mind that these quantitative results are relative to the hypothesis for thermal properties of both concrete and ground (respectively λc=1.8 W.K -1 .m -1 , ( ) =2.11 MJ.K -1 .m -3 and λm=2.4 W.K -1 .m -1 , ( ) =2.3 MJ.K -1 .m -3 ). However, it seems the 1-R model is not relevant for computation of the temperature evolution over short duration, typically of the order of some hours.