Franco Fusco
email: franco.fusco@ls2n.fr

Olivier Kermorgant
email: olivier.kermorgant@ls2n.fr

Philippe Martinet
email: philippe.martinet@inria.fr

Constrained Path Planning using Quadratic Programming

Sampling-based planning algorithms have been extensively exploited to solve a wide variety of problems. In recent years, many efforts have been dedicated to extend these tools to solve problems involving constraints, such as geometric loop-closure, which lead the valid Configuration Space (CS) to collapse to a lower-dimensional manifold.

One proposed solution considers an approximation of the constrained Configuration Space that is obtained by relaxing constraints up to a desired tolerance. The resulting set has then non-zero measure, allowing to exploit classical planning algorithms to search for a path connecting two given states. When constraints involve kinematic loops in the system, relaxation generally bears to undesired contact forces, which need to be compensated during execution by a proper control action.

We propose a new tool that exploits relaxation to plan in presence of constraints. Local motions inside the approximated manifold are found as the result of an iterative scheme that uses Quadratic Optimization to proceed towards a new sample without falling outside the relaxed region. By properly guiding the exploration, paths are found with smaller relaxation factors and the need of a dedicated controller to compensate errors is reduced. We complete the analysis by showing the feasibility of the approach with experiments on a real platform.

I. INTRODUCTION

Planning a motion via exact techniques for a robot with many degrees of freedom (dof) is hard and often impractical. Sampling-based planners overcome computational issues by exploring only a finite subset of the Configuration Space (CS) and checking for collisions at the considered states, rather than attempting to obtain a full representation of the C-obstacle regions. In particular, algorithms like Probabilistic Road Maps (PRM) [START_REF] Kavraki | Probabilistic roadmaps for path planning in high-dimensional configuration spaces[END_REF] and Rapidly-exploring Random Trees (RRT) [START_REF] Lavalle | Rapidly-exploring random trees: A new tool for path planing[END_REF] approximate the valid CS of a system by means of a graph whose nodes correspond to configurations and edges to motions between pairs of samples. To expand the graph, random configurations can be obtained by drawing each component independently from a given random distribution, while local motions are often created using linear interpolation.

These techniques can be applied to a number of contexts, but when dealing with tasks like sliding an object on a surface or carrying a load with two arms their performances become unsatisfactory. The reason is that in presence of such category of constraints the CS degenerates to a manifold with zero-measure. The probability of randomly sampling a valid configuration belonging to an n-dimensional subset using uniform distributions is directly related to the measure of the set itself, being null in case of lower-dimensional manifolds. Furthermore, linear interpolation fails to connect valid samples, since it produces straight paths that likely fall outside the valid CS.

Different strategies can be found in the literature to tackle the problem and extend randomized techniques to the constrained setup, to generate both feasible samples [START_REF] Cortes | Sampling-based motion planning under kinematic loop-closure constraints[END_REF] and valid motions. A projection approach is used by the Constrained Bi-Directional Rapildy-exploring Random Tree (CBiRRT) [START_REF] Berenson | Manipulation planning on constraint manifolds[END_REF] which exploits the gradient descent algorithm to push invalid configurations onto the manifold. Two samples are connected by alternately performing a linear step toward the goal and a projection of the resulting (infeasible) configuration back to the valid set. Other recent works [START_REF] Jaillet | Path planning under kinematic constraints by rapidly exploring manifolds[END_REF], [START_REF] Voss | Atlas+ x: Sampling-based planners on constraint manifolds[END_REF] rely on the construction of an Atlas whose charts consist in local linearizations of the constraint manifold. A planner is then required to explore the charts rather than the manifold, using linear interpolation while moving inside a single chart.

The concept of relaxation [START_REF] Bonilla | Sample-based motion planning for soft robot manipulators under task constraints[END_REF], [START_REF] Bonilla | Noninteracting constrained motion planning and control for robot manipulators[END_REF] consists in allowing a small constraint violation during planning. As a result, the valid CS is approximated by a full-dimensional set, which can then be explored by standard sampling-based planners. The technique was exploited to plan motions for compliant closed-loop systems, for whom infeasible samples lead to undesired internal forces. These efforts are compensated by a proper control action which steers the configuration back to the constraint manifold.

Since relaxation highly relies on system's compliance and on the control action, it might be infeasible in many scenarios involving rigid robots. These systems would need paths with lower constraint violation directly from the planning step, which can be achieved only by reducing the considered relaxation factor. However, by doing so the topology of the valid Configuration Space changes to a set of extremely narrow passages and the planning time increases significantly.

We propose a new approach based on relaxation which allows to find, in a shorter amount of time, paths characterized by lower constraint violation. The objective is to reduce the necessity of a control action during execution. We support the feasibility of the approach by including real experiments on a dual-arm system that carries an object while avoiding obstacles.

The rest of the paper is organized as follows. In Section II constraints are formally introduced and the relaxed Configuration Space is defined, while the planning algorithm is detailed in Section III. Real experiments involving geometric loop-closure constraints were performed to test the planner, and are discussed in Section IV. We report our conclusions and future work in Section V.

II. CONSTRAINTS AND RELAXATION

In this Section constraints are introduced, distinguishing between two classes, and a formal definition of the relaxed manifold is presented. We then detail the formulation of loopclosure constraints exploited in real experiments.

The Configuration Space C is assumed to be an ndimensional subset of R n , i.e., C ⊂ R n , n corresponding to the dof of the system. The configuration vector, denoted as q = q 1 • • • q n T , is assumed to have bounded

coordinates q i ∈ [q i,min , q i,max ] (i = 1, • • • , n).
Constraints are distinguished in two classes: (1) differentiable and (2) non-differentiable. The first category is represented by a set of n 1 functions in the form C i : C → R, for whom the gradient vector is assumed to be well defined and available in analytic form. They define the constraint manifold as the set M = {q ∈ C : C i (q) = 0, ∀i}. Examples from this category are kinematic loop-closure or orientation constraints on the end-effector.

The second class of constraints is instead described by n 2 functions D i : C → {0, 1}. A configuration q satisfies a constraint in this form if D i (q) = 1. This formulation can be used, e.g., to describe the collision state of a robot, and in general those constraints that do not cause the CS to collapse to lower-dimensional sets.

Relaxation is introduced by considering n 1 constants ε i > 0, each one quantifying the maximum allowed violation of the corresponding differentiable constraint. This bears to a natural definition of the relaxed Configuration Space as:

C R = {q ∈ C : |C i (q)| ≤ ε i ∩ D i (q) = 1 ∀i, i } (1)
As functions C i provide a numerical value for constraints violation, we adopt in the rest of this work the terminology error to refer to the value of constraints at a given configuration q. We therefore introduce the vector e = C 1 (q) • • • C n1 (q)

T to group all constraints violations and introduce the Jacobian matrix J e = ∂e ∂q , whose i-th row is given by the gradient (as a row vector) of C i . In order to explore C R , the Jacobian is used to locally approximate M, and motions are performed with the objective of minimizing the norm of e.

0 T l 0 T l T r F l F r F 0
Fig. 2. Reference frames and transformations in a dual-arm system.

We conclude this section by showing a possible choice of constraint functions to describe loop-closure for a dual-arm system rigidly grasping an object with both end-effectors. Fig. 2 shows a sketch of a dual-arm system, wherein three frames are considered: F 0 (world reference frame), F l (left end-effector frame) and F r (right end-effector frame). We assume the full configuration vector to be partitioned as q = q T l q T r T .

Rigid grasp of an object constrains relative translation and rotation of frames F l and F r , that is, the transformation matrix l T r needs to be constant. Let t be the desired translation and l R r the desired rotation matrix in l T r . It is then possible to define the constraint error as:

e . = t l + 0 R l t -t r θu (2)
t l and t r being the position of the end-effectors in world frame, and θu the angle-axis decomposition of the rotation matrix 0 R T r 0 R l l R r . Finally, the Jacobian writes:

J e = MJ l -J (v) r B 0 R T r J (ω) l -B 0 R T r J (ω) r (3)
where: J l is the Jacobian matrix of the left end-effector,

J (v) r
the right arm Jacobian of linear velocity, J

l , J

(ω) r are the left and right Jacobians of angular velocity. Matrices M and B are given by:

M = I 3 -0 R l t × (4)
B = I 3 - θ 2 [u] × + 1 - sinc θ sinc 2 θ 2 [u] 2 × (5)
I 3 representing the 3-by-3 identity matrix and [•] × the skewmatrix operator. The derivation of the Jacobian is detailed in the Appendix.

III. PLANNING WITH CONSTRAINTS

In this section the relaxation-based planner is presented. The core of the algorithm, presented in Section III-A, is a local router that tries to reach a given configuration starting from an initial sample. It produces a discrete set of configurations, each one being obtained after solving a Quadratic Optimization. In Section III-B we show how the connection algorithm is integrated in a Bi-directional RRT to explore the manifold and connect the start configuration to the goal one. We finally propose some post-processing techniques to refine a path in Section III-C.

A. Local Router

The idea behind the local routing algorithm is that the motion toward a desired configuration can be generated by optimizing both the current distance to the local goal and the error vector e. The two objectives will be in general in contrast, since to approach the goal a system will often be required to leave the constraint manifold. Quadratic Programming (QP) is then exploited to find a compromise that allows to shorten the distance towards the goal without excessively violating the constraints.

To enforce also non-differentiable constraints such as collision checking, we require the optimization to produce a sample that is not far away from the current configuration. If at the new point q new the condition D i (q new ) = 1 is satisfied for every non-differentiable constraint, a new step can be performed using QP. Thus, the local motion is generated and validated iteratively by repeating the two phases.

In the sequel, we firstly show how a single optimization is formulated and solved and then we detail the full iterative scheme.

1) Quadratic Programming: For the purposes of this work, a QP optimization problem is modeled as follows:

q (j+1) = arg min q (j+1) ∈C Q (j) q (j+1) -v (j) 2 (6)
subject to:

A (j) q (j+1) ≤ b (j) (7)
The optimization variable, q (j+1) , corresponds to the new configuration that will be appended to the local motion. Since the router will produce a number of intermediate samples, we use the index j to denote the last generated configuration. The objective is modeled via Q (j) ∈ R m×n and v (j) ∈ R m . The superscript (j) is used to underline that these quantities are constant during the optimization, fully determined by the value of q (j) . The optimization is subject to the set of linear inequalities given by (7) (A (j) and b (j) being of dimension p × n and p respectively).

In a first instance, the minimization of q (j+1) -q g 2 + α e (j+1) lin 2 could be considered, q g being the local goal configuration, α a diagonal weighting matrix and e (j+1) lin the linearized error corresponding to q (j+1) , given by: e (j+1) lin = e (j) + J (j) e q (j+1) -q (j) [START_REF] Bonilla | Noninteracting constrained motion planning and control for robot manipulators[END_REF] This formulation however presents a drawback: in early steps the distance to q g will be higher, and it will reduce after each iteration. The relative weight between the two terms in the objective would then be non-uniform during the iterative procedure, resulting in generally higher violations at the beginning of a local motion. A better behavior can be obtained by minimizing in the objective the distance to an intermediate point q (j) g , obtained using linear interpolation from q (j) to q g , i.e., q (j) g = q (j) +δ q g -q (j) , for a given scalar δ ∈ (0, 1] selected to ensure homogeneous relative weight of the two objectives.

The objective is then defined, in a single iteration, by the following values:

Q (j) = I n α J (j) e v (j) = q (j) g α J (j) e q (j) -e (j) (9)
To complete the definition of a QP instance, the set of inequalities ( 7) is now detailed. As we assume each coordinate of the configuration vector to be bounded, as mentioned in Section II, a first set of inequalities is exploited to enforce these limits:

q min ≤ q (j+1) ≤ q max (10)
A further set of inequalities imposes a maximum local displacement on each coordinate during a single iteration:

q (j) -β k ∆q ≤ q (j+1) ≤ q (j) + β k ∆q (11)
the entries of vector ∆q being positive constants and introducing β ∈ (0, 1) and k ∈ N. These inequalities play a double role in the optimization: on one hand, they prevent the distance between subsequent configurations from being too large. This is fundamental in order to consistently check non-differentiable constraints. Considering collision as an example, small obstacles could be wrongly overtaken if the path is too sparse, as depicted in Fig. 3(a). On the other hand, the inequalities can be used to reduce the impact of linearization introduced in (8). Depending on the manifold, a first-order approximation may not well reflect the topology of the valid set, as exeplified in Fig. 3(b). Thus, a step guided by the error Jacobian might lead to a sample that is indeed outside the relaxed CS. Whenever this happens, it is possible to shrink the allowed displacement by increasing k and to repeat the optimization step.

(a) (b) Fig. 3. Pictorial justification of inequalities [START_REF] Sucan | The Open Motion Planning Library[END_REF]. In (a) an obstacle (in gray) is not detected due to a large step, and two configurations are wrongly connected. In (b) a manifold is represented as the bold black line, and it is linearized around two points. Due to the high curvature at the red point, a new configuration could easily fall outside the relaxed region (in light blue) even with a small step, while around the green point the linear approximation is valid in a larger area.

2) Iterative Scheme: The motion connecting two configurations q s and q g is found by iteratively applying the procedure illustrated above to produce new intermediate samples. Algorithm 1 contains the pseudo-code of the local router. At the beginning of an iteration, the QP instance is initialized by computing Q (j) , v (j) and A (j) . An inner cycle tries to generate a new sample q (j+1) for increasing values of the parameter k. If the new configuration falls outside the relaxed CS, the optimization is repeated using a smaller maximum step size.

Algorithm 1 QP-based Motion Validator (QPMove) 1: QPMove(qs, qg) : 2: q (0) ← qs 3: σ ← { } 4: d (0) ← qaq b 5: for j=0 to jmax do 6: b (j) ← GET QP B VECTOR(q (j) , k) 13:

k ← 0 7: Q (j) , v (j) , A (j) ← INIT QP(q (j) , q b ) 8:
q (j+1) ← SOLVE QP(Q (j) , v (j) , A (j) , b (j) )
14: e (j+1) ← EVALUATE ERROR(q (j+1) ) 15:

k ← k + 1

16:

while not(-ε ≤ e (j+1) ≤ ε)

17:

if ∃i : Di q (j+1) = 0 then 18:

return "stop", σ 19:

end if 20: σ ← σ ∪ q (j+1)
21:

d (j+1) ← q (j+1) -q b 22: if d (j+1) ≤ dmin then

23:

return "success", σ

24:

end if

25:

if d (j+1) ≤ d (j) and d (j) -d (j+1) ≤ ∆d (-) then

26:

return "stop", σ

27:

end if

28:

if d (j+1) ≥ d (j) and d (j+1) -d (j) ≥ ∆d (+) then

29:

return "stop", σ 30: end if 31: end for 32: return "stop", σ

As soon as a new configuration lying in the relaxed region is found, non-differentiable constraints are evaluated at the sample. If any of them is violated the routine must be interrupted and the local path is returned, otherwise q (j+1) can be appended to the local path σ.

Four criteria are finally used to check if the local router should stop:

1) If the current distance d (j+1) to the goal configuration is smaller than a given threshold, the connection attempt is considered valid. 2) If an iteration produced a very small improvement in terms of distance d (j) to the target, the search is stopped, as a stationary point might have been reached. 3) As not only the distance to the goal is optimized, but also the linearized error, a new configuration might be located farther from the goal than the previous sample. Such configurations are rejected if the increment in distance exceeds the limit. 4) A limit of iterations is set as well, to prevent spending too much time on the connection of two samples. Start configurations are represented as blue circles, while the goals as blue squares. The constrained Configuration Space is a sinusoidal wave. Linear interpolation cannot produce any valid sample, while our strategy is able to proceed toward the goal. The projection approach moves slowly due to the fact that the goal is found almost orthogonally to the constraint.

Fig. 4 contains an example of motion computed using the proposed router. Starting from the given sample, linear interpolation cannot produce any feasible motion toward the goal configuration due to the small relaxation. Instead, using the QPMove function it is possible to move inside C R by keeping e sufficiently small (Fig. 4(a)). A comparison can be made with respect to classical projection approaches, e.g., CBiRRT. These planners usually perform motions by firstly stepping toward the goal and then re-projecting the sample on the manifold. In the case of Fig. 4(b) the goal is found in a direction that is almost perpendicular to the constraint curve at the start configuration. This is an unlucky condition for a projection approach, since a linearly interpolated sample is reprojected almost at the same initial position of the previous configuration. Despite being more expensive in terms of computation time, our approach usually explores C R faster than a projection technique, not blocking in these harder configurations.

B. Sampling-based Planner

The connection routine presented above can be easily inserted in classical Sampling-based planning algorithms. We decided during tests to use a bidirectional RRT [START_REF] Kuffner | RRT-connect: An efficient approach to single-query path planning[END_REF], replacing the extend function with our custom router. In its standard version, the bidirectional RRT firstly picks a random configuration and then selects the nearest neighbor from the current tree to attempt a connection. The motion is however limited to a maximum length, which may prevent from reaching the random sample. In this study, it was considered a more aggressive approach, firstly proposed in [START_REF] Lavalle | Rapidly-exploring random trees: Progress and prospects[END_REF], that does not limit the length of an extension step.

Regarding state sampling, we exploit the simplest uniform sampler that can extract any point from the Configuration Space of the system with equal probability. This is justified by the fact that the connection algorithm can produce new motions even when heading toward infeasible samples. In these cases, it will simply stop when the current branch cannot be further extended. However, since the motions are always constrained on the manifold, the new branches that are added to the trees could cover the relaxed CS non uniformly. In the future, it will be worth investigating the use of better sampling techniques.

C. Post Processing

It is well-known that randomized planners usually find rather irregular paths which would result in jerky trajectories during execution. To smooth and shorten the solution path, a number of shortcuts is therefore attempted. The operation works by selecting two nodes at random from the final path and trying to connect them. We again exploit QPMove for this purpose, and the new local motion, if successfully found, replaces the existing one if it reduces the total length of the overall path.

Finally, as we handle constraints using a relaxation approach their violation will not be completely nullified along local paths. A further post-processing phase is thus performed to refine a solution by better enforcing the constraints at each intermediate sample. This operation is carried out by solving a number of QP instances formulated in a similar way as the ones considered before. The main difference is that regarding the objective the local goal is set to be the initial path sample. This allows to optimize the error without drastically change the path. It was noted, however, that in some cases this operation might invalidate one or more samples. As an example, if a path point was initially located really near to an obstacle, after some iterations the path could be moved in such a way that collision is not satisfied anymore. To deal with such problem, we try to re-optimize the initial path point by slightly changing the position of the local goal in the opposite direction in which collision or other non-differentiable constraints were invalid.

IV. EXPERIMENTAL RESULTS

We present in this Section some experimental results obtained using the proposed planner. The algorithm was tested to find motions for a dual-arm system featuring two Kuka LWR4 arms (Fig. 1), having 14 dof in total. The robot is required to move an object grasped with both hands from an initial configuration to a given one while avoiding obstacles. The bi-manual grasp adds six constraints to the system, due to the presence of a closed kinematic loop in the structure, which is described via (2) and (3).

The Open Motion Planning Library [START_REF] Sucan | The Open Motion Planning Library[END_REF] was used to implement planning-related capabilities, while ROS [START_REF] Quigley | ROS: an open-source robot operating system[END_REF] and MoveIt! [START_REF] Sucan | Moveit![END_REF] are used to send and execute planning requests. Obstacles are assumed to be static and with known pose, as perception of the environment is out of the scope of this study.

Different simulations were run to experimentally tune the parameters of QPMove. Regarding the choice of β, values ranging between 0.7 and 0.85 gave the best results. It was noted that high values of constraint weights in the objective increase the planning time, as less importance is given to the exploration of C R . Motions get often stuck due to very small improvements or by reaching the maximum number of iterations. Concerning the choice of the relaxation factor, rather tolerant violations were allowed during the planning phase (around 1 cm for translation and 1 • for rotation). Although these values would result inappropriate during execution, by properly tuning the weights for error minimization it is possible to achieve better results in practice. In addition, while attempting to smooth and simplify the path by shortcut operations higher weights can be used. As a result violations can be reduced even without the final error refinement operation to be less than 1 mm and 1 mrad in average.

In Fig. 5 some frames taken from a real experiment are shown. The robot is required to lift a metallic bar while avoiding some obstacles, two of them forming a narrow passage. Along the planned path the average violation is of 0.27 mm and 0.05 • without running the final refinement process (only shortcuts were performed). By repeatedly attempt to solve the query, planning times were recorded. The average runtime is 4.168 s (excluding the post-processing phase), with a minimum of 0.931 s. Standard relaxation techniques were tested for a comparison, and no feasible plan was found within 60 s, using relaxation factors of 100 mm and 17.2 • . An example of error evolution during actual execution is shown in Fig. 6. It can be noted that violations increase significantly with respect to the values obtained after the planning phase, raising to a maximum of 5 mm in translation and 0.6 • in rotation. This results from two factors: the absence of a control loop that regulates the error during the actual motion and a rough estimation of the dynamic parameters of the end-effectors. Depending on the final application these values may not be acceptable, and thus a regulating control action would be required to compensate them.

V. CONCLUSIONS

In this work we proposed a new approach that exploits Quadratic Programming to reduce constraints violations in a relaxation-based context. Exploiting the Jacobian matrix of constraints functions allows to formulate the problem of moving toward a random sample as a local optimization of the distance to the configuration and the error associated with differentiable constraints.

As a result, paths of higher quality can be found in a shorter time, reducing the need of control actions at the execution level. The feasibility of our algorithm is confirmed by successful experiments on a real platform.

A drawback of this approach is the presence of a number of parameters which need to be tuned experimentally, and whose value can affect significantly the planning time.

A naive sampling technique was exploited to generate random configurations in the full Configuration Space of the robot. A better sampling approach able to uniformly cover the constraint manifold might lead to improvements in planning time and quality of the resulting path, and should therefore be investigated and integrated in the proposed approach.

APPENDIX

We show in the sequel the derivation of (3) starting from error's definition. The translation component is easily identified by considering that as t is constant, it holds:

ṫl + ˙ 0 R l t -ṫr = v l + ω l × 0 R l t -v r (12)
By definition, the linear/angular velocities v l , v r and ω l are linked with joint vector's derivatives by robot's Jacobian matrices J l and J r , thus leading to the first part of the error Jacobian. Regarding rotation, one needs to evaluate ˙ θu. The derivative of the angle-axis representation of a rotation matrix a R b is in general given by B a ω b/a , with a ω b/a representing the angular velocity vector of the generic frame b with respect to a (expressed in frame a). The form of matrix B was already given in [START_REF] Jaillet | Path planning under kinematic constraints by rapidly exploring manifolds[END_REF]. In the considered case the rotation matrix for which the decomposition is evaluated is 0 R T r 0 R l l R r . It is possible to consider a virtual frame F r rigidly attached to the origin of the left end-effector, characterized by the rotation l R r = l R r . In this sense 0 R T r 0 R l l R r represents the rotation of the virtual frame F r with respect to F r , and therefore:

˙ θu = B r ω r /r = B 0 R T r (ω r -ω r ) (13)
In addition, as F r and F l are rigidly linked to each other, ω r = ω l . Substituting ω l = J

(ω) l ql and ω r = J

(ω) r qr bears to the final relation.

Fig. 1 .

 1 Fig. 1. The experimental platform at LS2N: a dual-arm robot featuring two Kuka LWR manipulators with a Schunk SDH (left hand) and Schunk SVH (right hand).

Fig. 4 .

 4 Fig.[START_REF] Berenson | Manipulation planning on constraint manifolds[END_REF]. Motions generated in C R after 10 iterations by QPMove (a) and a projection approach (b). Start configurations are represented as blue circles, while the goals as blue squares. The constrained Configuration Space is a sinusoidal wave. Linear interpolation cannot produce any valid sample, while our strategy is able to proceed toward the goal. The projection approach moves slowly due to the fact that the goal is found almost orthogonally to the constraint.

Fig. 5 .

 5 Fig. 5. Execution of a planned path, from initial configuration (left) to the final one (right). The full execution is included in the attached video.

Fig. 6 .

 6 Fig. 6. Evolution of constraint errors during real execution. Translation and rotation components are reported respectively in the top and bottom pictures.

ACKNOWLEDGMENT

This work was carried out in the framework of the PROMPT project, a project funded by RFI Atlanstic 2020. Parts of the equipment used here were funded by the project ROBOTEX, reference ANR-10-EQPX-44-01.


Petra Bragée

Florence Mazier

Peter Rosén

Anna Brostrom

Wilhelm Granéli

Dan Hammarlund

Anna Broström

Recent brownification of South Swedish lake waters -an effect of climate change or land use?

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

A. Introduction

Dissolved organic carbon (DOC) content, especially terrestrial humic matter (DHM), and water colour of many lakes and streams have increased substantially during the last decades (Evans et al. 2004, Roulet and[START_REF] Roulet | Browning the waters[END_REF]. Brownification not only reduces the quality of these waters as drinking water supplies and for recreational purposes, but lake biodiversity is also expected to be effected. The aim of this project is to evaluate the potential underlying causes of the recent brownification of South Swedish lakes. A multiproxy approach is applied (fig 3 .) for reconstruct on of the biogeochemical and biological changes within selected lakes, and the vegetation/land use changes within their catchments for the past c. 1000 years. Here we present the study area, main methods, preliminary results and links to Fredh et al. project.

The study area is situated in southern Sweden (fig 1 .). The selection of lakes was based on; 1. recent water colour trends (fig 2.), 2. size and 3. catchment properties. Lake Åbodasjön is 50 ha and the catchment is influenced by anthropogenic activities. Lake Lindhultsgöl is 7,1 ha and the surroundings consist of forests and wetlands. Field work was carried out in April 2008 and sediments were obtained with a Kajak corer and a Russian peat corer.

Lake Åbodasjön

TOC concentration

Past changes in lake-water total organic carbon (TOC) concentration are based on near-infrared spectroscopy (NIRS) of the lake sediments. Recent studies have shown that NIRS of lake sediments exhibit good correlations with TOC concentrations of lake-water [START_REF] Rosén | Total organic carbon (TOC) of lake water during the Holocene inferred from lake sediments and near-infrared spectroscopy (NIRS) in eight lakes from northern Sweden[END_REF]. The method utilizes the fact that the organic fraction of lake sediments has a distinctive NIR signal that can be summarized using multivariate statistical tools. New models for South Swedish lakes with TOC content variations up to 30 mg/l are in progress.

Catchment vegetation

A novel modelling approach, the Landscape Reconstruction Algorithm (LRA), will be used to quantify catchmentscale vegetation changes based on pollen data. Past land cover in terms of absolute areas of different landscape units (forest, open land, pastures, cultivated areas, etc.) will be estimated (see also poster by Daniel Fredh et al.).

E. Next steps

• Establish an accurate chronology.

• Apply the same approach and methods for the sediment sequence from Lake Lindhultsgöl.

• Quantitative vegetation reconstructions using the LRA modelling.

• Sample the lakes with a freezecorer to enable continuous and more highly resolved subsampling.

• The flagellate Gonyostomum semen (fig. 5.) is an invasive algae which is favoured by increased lake water DHM content [START_REF] Cronberg | Mass development of the flagellate Gonyostomum semen (Raphidophyta) in Swedish forest lakes -an effect of acidification?[END_REF]. If its resting cysts can be found in the lake sediments, information about its occurence and its connection to different changes in the lake, TOC concentration for example, might be evaluated. Contact information

Figure 1 .

 1 Figure 1. Location of the study sites including Lake Fiolen, used for the LRA modelling (see Daniel Fredh´s poster).

Figure 3 .

 3 Figure3. Methods applied to the sediment sequence from Lake Åboda. Table of the methods for the sediment of the Lake Åboda. Green =finished, Yellow=in progress and Red=planned.

Figure 4 .

 4 Figure4. Radiocarbon dates, NIRS-inferred lake-water TOC concentration, magnetic susceptibility and preliminary pollen data from Lake Åboda plotted against depth below the water surface.

Figure 2 .

 2 Figure 2. The two lakes show distinct increases in lake-water TOC concentration in the upper c. 4 cm of the sediments.

Chronology

Accurate chronologies will be based on radioisotope modelling ( 210 Pb and 137 Cs) in combination with radiocarbon dating. Core correlation is based on highly resolved proxy data.

B. Study sites


Dorian Dupuy

Adrien Toutant
email: adrien.toutant@univ-perp.fr

Françoise Bataille

Equations of energy exchanges in variable density turbulent flows

come

Introduction

This paper addresses the energy exchanges in turbulent flows with highly variable fluid properties. This covers flows with a high Mach number (high speed flows), such as the flows around a high-speed aircraft, or through a high speed jet or a nozzle [START_REF] Gatski | Compressibility, turbulence and high speed flow[END_REF], and low Mach number flows submitted to a strong temperature gradient, found for instance in heat exchangers, propulsion systems or nuclear or concentrated solar power plants [START_REF] Bataille | Eddy viscosity and diffusivity modeling[END_REF][START_REF] Mathelin | Theoretical investigation of some thermal effects in turbulence modeling[END_REF][START_REF] Daguenet-Frick | Numerical investigation of a ceramic high-temperature pressurized-air solar receiver[END_REF][START_REF] Serra | Thermal large eddy simulation in a very simplified geometry of a solar receiver[END_REF][START_REF] Serra | High-temperature gradient effect on a turbulent channel flow using thermal large-eddy simulation in physical and spectral spaces[END_REF][START_REF] Serra | Turbulent kinetic energy spectrum in very anisothermal flows[END_REF][START_REF] Toutant | Turbulence statistics in a fully developed channel flow submitted to a high temperature gradient[END_REF]. The study of the energy exchanges between the different parts of total energy is a useful tool for both turbulence modelling and the fundamental understanding of turbulence. More detailed information is obtained through the study of the energy exchanges in the spectral domain [START_REF] Lumley | Spectral energy budget in wall turbulence[END_REF][START_REF] Domaradzki | Local energy transfer and nonlocal interactions in homogeneous, isotropic turbulence[END_REF][START_REF] Domaradzki | Energy transfer in numerically simulated wall-bounded turbulent flows[END_REF][START_REF] Zhou | Degrees of locality of energy transfer in the inertial range[END_REF][START_REF] Zhou | Interacting scales and energy transfer in isotropic turbulence[END_REF][START_REF] Domaradzki | Energy transfer in numerically simulated wall-bounded turbulent flows[END_REF][START_REF] Marati | Energy cascade and spatial fluxes in wall turbulence[END_REF][START_REF] Dunn | Analysis of the energy budget in turbulent channel flow using orthogonal wavelets[END_REF][START_REF] Bolotnov | Spectral analysis of turbulence based on the DNS of a channel flow[END_REF][START_REF] Cimarelli | Paths of energy in turbulent channel flows[END_REF][START_REF] Cimarelli | Spectral enstrophy budget in a shear-less flow with turbulent/non-turbulent interface[END_REF][START_REF] Cimarelli | Sources and fluxes of scale energy in the overlap layer of wall turbulence[END_REF][START_REF] Cimarelli | Cascades and wallnormal fluxes in turbulent channel flows[END_REF][START_REF] Lee | Spectral analysis on reynolds stress transport equation in high firewall-bounded turbulence[END_REF][START_REF] Mizuno | Spectra of energy transport in turbulent channel flows for moderate reynolds numbers[END_REF]. However, while kinetic energy is fundamental property of any flow, it is not the case of its decomposition into turbulence kinetic energy and mean kinetic energy.

In incompressible flows with constant fluid properties, such decomposition is unique. The averaged kinetic energy is decomposed clearly, unambiguously and straightforwardly into the sum of two contributions: the kinetic energy of the mean motion associated with the mean velocity and the kinetic energy of the turbulent motion associated with the velocity fluctuation [see e.g. [START_REF] Wilcox | Turbulence modeling for CFD[END_REF][START_REF] Pope | Turbulent Flows[END_REF]. In compressible flows with highly variable density, this analysis is hindered by additional density-velocity correlations. The decomposition of kinetic energy becomes more complex and arbitrary. It is even more difficult in the spectral domain. The choice ultimately depends on the physical role given to the additional density-velocity correlations with respect to what constitutes the mean motion and the turbulent motion [START_REF] Chassaing | Variable density fluid turbulence[END_REF]. The most popular and successful decomposition extends the incompressible decomposition to the compressible case through the introduction of a density weighted averaging. This decomposition was widely developed by Favre [START_REF] Favre | The equations of compressible turbulent gases[END_REF][START_REF] Favre | Statistical equations of turbulent gases[END_REF][START_REF] Favre | Turbulence: Space-time statistical properties and behavior in supersonic flows[END_REF]. Since, it has been used extensively by various authors [START_REF] Lele | Compressibility effects on turbulence[END_REF][START_REF] Huang | Compressible turbulent channel flows: DNS results and modelling[END_REF][START_REF] Guarini | Direct numerical simulation of a supersonic turbulent boundary layer at Mach 2.5[END_REF][START_REF] Maeder | Direct simulation of turbulent supersonic boundary layers by an extended temporal approach[END_REF][START_REF] Pirozzoli | Direct numerical simulation and analysis of a spatially evolving supersonic turbulent boundary layer at M=2.25[END_REF]. Another approach, the mixed weighted decomposition, mixes density weighted averaging and Reynolds averaging. It was first introduced by Bauer et al. [START_REF] Bauer | A numerical method and an extension of the korst jet mixing theory for multispecie turbulent jet mixing[END_REF] and further studied by Ha Minh et al. [START_REF] Ha Minh | A new approach to the analysis of turbulent mixing in variable density flows[END_REF][START_REF] Ha Minh | The turbulence modelling of variable density flows-a mixed-weighted decomposition[END_REF]. In this formulation, kinetic energy is seen as the product of the velocity and the density weighted velocity. In a third method, kinetic energy is decomposed using a change of variable based on the density square root weighted velocity. This decomposition was first proposed by Yih [START_REF] Yih | Exact solutions for steady two-dimensional flow of a stratified fluid[END_REF] then adopted by various authors [START_REF] Kida | Energy and spectral dynamics in decaying compressible turbulence[END_REF][START_REF] Rodriguez | Compressible rotational flows generated by the substitution principle[END_REF][START_REF] Cook | Energy transfer in Rayleigh-Taylor instability[END_REF][START_REF] Aulery | Energy transfer process of anisothermal wall-bounded flows[END_REF][START_REF] Aulery | Spectral analysis of turbulence in anisothermal channel flows[END_REF]. This change of variable allows the study of kinetic energy to be extended easily to the spectral domain. Finally, Chassaing [START_REF] Chassaing | Une alternative à la formulation des équations du mouvement turbulent d'un fluide à masse volumique variable[END_REF] [see also [START_REF] Cousteix | Turbulence models for compressible flows[END_REF][START_REF] Aupoix | Introduction to turbulence modelling for compressible flows[END_REF][START_REF] Chassaing | The modeling of variable density turbulent flows. a review of first-order closure schemes[END_REF][START_REF] Chassaing | Variable density fluid turbulence[END_REF] suggested the decomposition of kinetic energy using the Reynolds averaging. From a modelling perspective, the use of the unweighted averaging may be advantageous in low Mach number flows, in which the energy conservation acts as a constraint on the divergence of the velocity [START_REF] Nicoud | Conservative high-order finite-difference schemes for low-Mach number flows[END_REF]. The square of the fluctuating velocity (without the density) is also encountered for instance in the modelling of two-phase flows [START_REF] Toutant | DNS of the interaction between a deformable buoyant bubble and a spatially decaying turbulence: a priori tests for LES two-phase flow modelling[END_REF] or in variable density flows, provided the momentum equation is divided by the density before averaging [START_REF] Sidharth | Filtered velocity based LES of mixing in high speed recirculating shear flow[END_REF]. In a variable density setting, the use of the Reynolds averaging necessarily leads to the decomposition of kinetic energy into three parts, called ternary decomposition. The kinetic energy is thus split into turbulence kinetic energy, mean kinetic energy and a mixed kinetic energy, related to both the mean and turbulent motions. However, we believe the underlying idea behind the ternary decomposition has not been taken to its logical conclusion as no interaction between the mixed kinetic energy and another part of total energy was identified.

This paper aims to establish a new formulation of the energy exchanges between the different parts of total energy in a ternary decomposition that gives to the mixed kinetic energy a full role. The formulation is compared to the formulation of Chassaing [START_REF] Chassaing | Une alternative à la formulation des équations du mouvement turbulent d'un fluide à masse volumique variable[END_REF] and the differences between the two formulations with regard to the physical interpretation of the terms are discussed. We then take the decomposition further and split the density into a mean and fluctuating part. This leads to the definition of the mean density part of total energy and the fluctuating density part of total energy. The mean density turbulence kinetic energy, product of the mean density and the half-trace of the velocity fluctuation correlation tensor, appears in the mean density part of the decomposition as exchanging energy with the other parts of total energy. This quantity is approximately equal to the turbulence kinetic energy in flows satisfying Morkovin's hypothesis [START_REF] Morkovin | The Mechanics of Turbulence[END_REF].

Once the new formulation of the energy exchanges established, we focus more specifically on the mean density turbulence kinetic energy. We establish its evolution equation in spectral domain, recognizing that the mean density turbulence kinetic energy has with the Reynolds averaging a clear spectral equivalent. The spectral equation extends the spatial equation to the spectral domain, associating to each spatial term a spectral equivalent. To the knowledge of the authors, this has not been achieved in the literature for variable density flows. A purely spectral term that redistributes the energy between scales is identified, as in the work of Lee and Moser [START_REF] Lee | Spectral analysis on reynolds stress transport equation in high firewall-bounded turbulence[END_REF] and Mizuno [START_REF] Mizuno | Spectra of energy transport in turbulent channel flows for moderate reynolds numbers[END_REF]. In order to carry out the Fourier transform, we consider a flow with two homogeneous and periodic directions. This does not lead to a loss of generality as the equations given may easily be adapted to a flow with one or three homogeneous directions.

The complete representation of the energy exchanges between the different parts of total energy is presented in section 2 and the equation of the mean density turbulence kinetic energy in the spectral domain in section 3.

2 Energy exchanges between the different parts of total energy in the ternary decomposition

General considerations

In this section, we describe a new formulation of the energy exchanges between the different parts of total energy in a ternary decomposition. We will establish the formulation obtained from the decomposition of velocity, but not density, with the Reynolds averaging, referred to as the one-stage formulation in this paper, and from the decomposition of both the velocity and density, referred to as the twostage formulation in this paper. The two-stage formulation is required to write the spectral equation of the mean density turbulence kinetic energy. We first define here a few useful quantities and give some general remarks on the derivation of the formulation.

The total energy per unit volume ρ(E + I) is a conservative quantity. Its components however are not as they exchange energy among themselves. In the following, the evolution equation of each part of total energy in the ternary representation will be given and we will identify the energy exchanges between these quantities. Many consistent formulations of the energy exchanges can be proposed. The formulation was devised according to the following criteria:

• Each term of the formulation must be either interpreted as a conservative energy transfer or an interaction with exactly one of the other parts of total energy.

• If a term is to be interpreted as a conservative energy transfer, it must be written in a conservative form, that is as a divergence; otherwise, it must be written in a non-conservative form.

• The formulation must be symmetrical, in particular with respect to the manner in which it deals with fluctuations and statistically averaged quantities.

• The formulation must correctly behave when considering a limit case such as laminar, homogeneous or incompressible flows. In particular, a quantity that becomes instantaneously equal to zero must not be associated with any energy exchange.

We consider a non-relativistic compressible flow with highly variable fluid properties under the continuity hypothesis. Without loss of generality, no body forces are taken into account which means gravity is neglected and there is no heat source. The flow is governed by the Navier-Stokes equations under the following form [START_REF] Garnier | Large eddy simulation for compressible flows[END_REF]:

• mass conservation ∂ρ ∂t + ∂ρU j ∂x j = 0, (1)
• momentum conservation

∂ρU i ∂t + ∂ρU j U i ∂x j = ∂Υ ij ∂x j , (2)
• energy conservation

∂ρI ∂t + ∂ρU j I ∂x j = ∂ ∂x j λ ∂T ∂x j + Υ ij ∂U i ∂x j , (3)
with ρ the density, T the temperature, I the internal energy per unit mass, t the time, U i the i-th component of the velocity, Υ ij the component of the total stress tensor with the i and j indices and x i the Cartesian coordinate in the i-th direction. Einstein summation convention is used. The total stress tensor Υ ij is given by the contributions of the viscous shear stress tensor and of the pressure stress. We will keep the total stress tensor undissociated throughout this paper because the pressure and viscous contributions are formally similar.

One-stage formulation

The instantaneous total energy per unit volume ρ(E +I) is the sum of the instantaneous kinetic energy per unit volume ρE and the internal energy per unit volume ρI. In the ternary decomposition, the kinetic energy is decomposed into three parts by splitting the velocity into a mean and fluctuating part [following 43], namely U i = U i + u i , where the overline ( ) denotes the statistical average and the prime symbol ( ) the fluctuating part. We use a lowercase u for the velocity fluctuation for a better visual differentiation but there is no further underlying differences. We obtain [START_REF] Chassaing | Une alternative à la formulation des équations du mouvement turbulent d'un fluide à masse volumique variable[END_REF]

ρE = 1 2 ρU i U i = ρE + ρe + ρe, (4)
with ρE = 1 2 ρU i U i the mean kinetic energy, associated with the mean motion, ρe = 1 2 ρu i u i the turbulence kinetic energy, associated with the turbulent motion, and ρe = ρu i U i the mixed kinetic energy, associated with both the mean and turbulent motion. This results in a fourfold decomposition of total energy. The evolution equation of total energy ρ(E + I) is given by

∂ρ(E + I) ∂t = Φ c + Φ Υ + Φ λ , (5)
with: Φ c the convection, Φ Υ the transfer by the total stress and Φ λ the transfer by conduction, given by:

Φ c = - ∂ρU j (E + I) ∂x j , (6)
Φ Υ = ∂Υ ij U i ∂x j , (7)
Φ λ = ∂ ∂x j λ ∂T ∂x j . ( 8
)
These three terms are conservative terms that represent the transfer of total energy by three different physical phenomena. With the fourfold decomposition of total energy, this equation becomes four equations associated with each part of total energy. This has two effects. First, the various conservative energy transfer terms are distributed among the four parts of kinetic energy. Secondly, additional nonconservative terms emerge. From the decomposition of total energy into kinetic energy and internal energy appears a new term E that represents the interaction between these two quantities. From the decomposition of kinetic energy into three terms appears two new terms P and P that represent an interaction between the different parts of kinetic energy.

The energy exchanges between the four parts of total kinetic energy may be written as:

∂ρE ∂t = Φ c + Φ Υ + P + E, (9)
∂ρe ∂t = ϕ c + ϕ Υ + P + , (10)
∂ρe ∂t = ϕ c -P -P, (11)
∂ρI ∂t = Φ T,c + Φ λ -E -, (12)
where we identify the following terms:

• the convection Φ c , decomposed into: Φ c associated with the mean kinetic energy, ϕ c associated with the turbulence kinetic energy, ϕ c associated with the mixed kinetic energy and Φ T,c associated with internal energy,

Φ c = - ∂ρU j E ∂x j , (13)
ϕ c = - ∂ρU j e ∂x j , (14)
ϕ c = - ∂ρU j e ∂x j , (15)
Φ T,c = - ∂ρU j I ∂x j , (16)
• the transfer by the total stress Φ Υ , decomposed into: Φ Υ associated with the mean kinetic energy and ϕ Υ associated with the turbulence kinetic energy,

Φ Υ = ∂Υ ij U i ∂x j , (17)
ϕ Υ = ∂Υ ij u i ∂x j , (18)
• the interaction between the turbulence kinetic energy and the mixed kinetic energy P,

P = -ρu i U j ∂U i ∂x j + ρu i U j ∂U i ∂x j -ρu i 1 ρ ∂Υ ij ∂x j , (19)
• the interaction between the mean kinetic energy and the mixed kinetic energy P,

P = -ρU i U j ∂u i ∂x j + ρU i U j ∂U i ∂x j -ρU i 1 ρ ∂Υ ij ∂x j , (20)
• the interaction between kinetic energy and internal energy E, decomposed into: E associated with the mean kinetic energy and associated with the turbulence kinetic energy,

E = -Υ ij ∂U i ∂x j , (21)
= -Υ ij ∂u i ∂x j . ( 22
)
This set of equations is represented in a schematic form in figure 1.

In the limit case U = 0, the mean kinetic energy E and the mixed kinetic energy e vanish as do all energy exchange terms associated with these two quantities. The formulation reduces to the sole exchanges between turbulence kinetic energy and internal energy and describes the well-known rate of decay of turbulence kinetic energy and the paired gain of internal energy [START_REF] Batchelor | The theory of homogeneous turbulence[END_REF]. Similarly, the formulation reduces to the exchanges between mean kinetic energy and internal energy in the limit case u = 0.

We compare the energy exchanges in the present formulation to the ternary decomposition of Chassaing [START_REF] Chassaing | Une alternative à la formulation des équations du mouvement turbulent d'un fluide à masse volumique variable[END_REF]. While Chassaing [START_REF] Chassaing | Une alternative à la formulation des équations du mouvement turbulent d'un fluide à masse volumique variable[END_REF] did not explicitly give the energy exchanges between the different parts of total energy, he gave enough information to identify the energy exchanges without ambiguity. The formulation of Chassaing [START_REF] Chassaing | Une alternative à la formulation des équations du mouvement turbulent d'un fluide à masse volumique variable[END_REF] may be written using the notations of this paper as:

∂ρE ∂t = Φ c + Φ Υ + Y + P + X -Y -X + E, (23)
∂ρe ∂t = ϕ c + ϕ Υ + P -X + X + , (24)
∂ρe ∂t = ϕ c -Y -P -X -P + X -Y , (25)
∂ρI ∂t = Φ T,c + Φ λ -E -, (26)
with:

X = -ρu i u j ∂U i ∂x j , (27)
Y = - ∂ρu j e ∂x j . (28)
It is represented in a schematic form in figure 2. The formulation presented in this paper and the formulation of Chassaing [START_REF] Chassaing | Une alternative à la formulation des équations du mouvement turbulent d'un fluide à masse volumique variable[END_REF] are mathematically equivalent. However, the energy exchanges identified and the interpretation given to the terms are different. The differences can be attributed to two main changes. First, the present formulation associates to the four parts of total energy a convective term related to both the mean and turbulent motion, that is of the transport by advection and diffusion. However, the formulation of Chassaing [START_REF] Chassaing | Une alternative à la formulation des équations du mouvement turbulent d'un fluide à masse volumique variable[END_REF] only associates an advective term to the mixed kinetic energy ρe, but no diffusive term. Due to this, the term Y , diffusion of mixed kinetic energy in the present formulation, instead appears in the evolution equation of the mean kinetic energy ρE and is interpreted as the power of the Reynolds stress through the mean motion. In addition, there is in the formulation of Chassaing [START_REF] Chassaing | Une alternative à la formulation des équations du mouvement turbulent d'un fluide à masse volumique variable[END_REF] a direct energy exchange X between turbulence kinetic energy and mean kinetic energy whereas any interaction between turbulence kinetic energy and mean kinetic energy occurs through the mixed kinetic energy in the present formulation. Both of these changes modify the energy exchange associated with mixed kinetic energy.

The formulation of Chassaing [START_REF] Chassaing | Une alternative à la formulation des équations du mouvement turbulent d'un fluide à masse volumique variable[END_REF] is due to these two differences more similar to the classical incompressible representation of the energy exchanges, in which the term Y appears in the evolution equation of the mean kinetic energy and where there is a direct energy exchange between turbulence kinetic energy and mean kinetic energy. However, some elements suggest that the present formulation is more physical. Indeed, it associates a full convective term to each part of kinetic energy. We consider this as necessary as it is part of the material derivative. Besides, the formulation is more symmetrical with regard to the manner it deals with fluctuations and statistically averaged quantities. The formulation is not modified if the statistical average operator and the fluctuating part operator are substituted in equations ( 13) to [START_REF] Favre | Statistical equations of turbulent gases[END_REF] and in the definitions of ρE, ρe and ρe. Finally, the formulation may be used to consider the instantaneous energy exchanges as it does not rely on simplifications only valid in the statistically averaged case. This ensures the consistency of the formulation, in the sense that the energy exchanges in both the instantaneous and statistically averaged cases are well-defined and are not conflicting. This consistency is important to give a physical interpretation to the energy exchange, as this lets us consider the statistically averaged energy exchanges as the statistical average of the associated instantaneous energy exchanges.

Two-stage formulation

The ternary decomposition of kinetic energy is taken further with the decomposition of density into a mean and fluctuating part ρ = ρ + ρ . Namely, total energy ρ(E + I) is decomposed into a mean density part ρ(E + I) and a fluctuating density part ρ (E + I). Similarly, the mean kinetic energy ρE, the turbulence kinetic energy ρe, the mixed kinetic energy ρe and the internal energy per unit volume ρI are decomposed into a mean density part, ρE, ρe, ρe and ρI respectively, and a fluctuating density part, ρ E, ρ e, ρ e and ρ I respectively. This results in a eightfold decomposition of total energy. In this eightfold decomposition of total energy, the terms of the ternary formulation are decomposed further in a mean and fluctuating density part. Any term α identified in the one-stage formulation is split into two terms in the two-stage formulation: α 0 , associated with the mean density part of total energy and α 1 , associated with the fluctuating density part of total energy. Moreover, additional terms appear that represent an interaction between the mean and fluctuating density part of total energy.

The energy exchanges between the eight parts of total energy may be written as:

∂ρE ∂t = Φ c 0 + Φ Υ 0 + P 0 + Z c + Z Υ + E 0 , (29) ∂ρ
E ∂t = Φ c 1 + Φ Υ 1 + P 1 -Z c -Z Υ + E 1 , (30)
∂ρe ∂t = ϕ c 0 + ϕ Υ 0 + P 0 + ζ c + ζ Υ + 0 , (31)
∂ρ e ∂t = ϕ c 1 + ϕ Υ 1 + P 1 -ζ c -ζ Υ + 1 , (32)
∂ρe ∂t = ϕ c 0 -P 0 -P 0 + ζ c , (33)
∂ρ e ∂t = ϕ c -P 1 -P 1 -ζ c , (34)
∂ρI ∂t = Φ T,c 0 + Φ λ 0 + Z T,c + Z λ -E 0 -0 , (35) ∂ρ
I ∂t = Φ T,c 1 + Φ λ 1 -Z T,c -Z λ -E 1 -1 , (36)
with

Φ c 0 = - ∂ρU j E ∂x j , Φ c 1 = - ∂ρ U j E ∂x j , ϕ c 0 = - ∂ρU j e ∂x j , ϕ c 1 = - ∂ρ U j e ∂x j , ϕ c 0 = - ∂ρU j e ∂x j , ϕ c 1 = - ∂ρ U j e ∂x j , Φ T,c 0 = - ∂ρU j I ∂x j , Φ T,c 1 = - ∂ρ U j I ∂x j , Φ Υ 0 = ∂ ∂x j Υ ij ρ ρU i , Φ Υ 1 = ∂ ∂x j Υ ij ρ ρ U i , ϕ Υ 0 = ∂ ∂x j Υ ij ρ ρu i , ϕ Υ 1 = ∂ ∂x j Υ ij ρ ρ u i , Φ λ 0 = ∂ ∂x j ρ ρ λ ∂T ∂x j , Φ λ 1 = ∂ ∂x j ρ ρ λ ∂T ∂x j , P 0 = -ρu i U j ∂U i ∂x j + ρu i U j ∂u i ∂x j -ρu i 1 ρ ∂Υ ij ∂x j , P 1 = -ρ u i U j ∂U i ∂x j + ρ u i U j ∂u i ∂x j -ρ u 1 ρ ∂Υ ij ∂x j , P 0 = -ρU i U j ∂u i ∂x j + ρU i U j ∂U i ∂x j -ρU i 1 ρ ∂Υ ij ∂x j , P 1 = -ρ U i U j ∂u i ∂x j + ρ U i U j ∂U i ∂x j -ρ U i 1 ρ ∂Υ ij ∂x j , E 0 = - ρΥ ij ρ ∂U i ∂x j , E 1 = - ρ Υ ij ρ ∂U i ∂x j , 0 = - ρΥ ij ρ ∂u i ∂x j , 1 = - ρ Υ ij ρ ∂u i ∂x j ,
and where we identify the following terms:

• the energy dilatation correlation Z c ,

Z c = (E + I) ∂ρU j -ρU j ∂x j , (37)
decomposed into: Z c associated with the mean kinetic energy, ζ c associated with the turbulence kinetic energy, ζ c associated with the mixed kinetic energy and Z T,c associated with the internal energy,

Z c = E ∂ρU j -ρU j ∂x j , (38)
ζ c = e ∂ρU j -ρU j ∂x j , (39)
ζ c = e ∂ρU j -ρU j ∂x j , (40)
Z T,c = I ∂ρU j -ρU j ∂x j , (41)
• the interaction between the mean and fluctuating density part of total energy by the total stress Z Υ ,

Z Υ = -Υ ij U i ∂ ∂x j ρ ρ , (42)
decomposed into: Z Υ associated with the mean kinetic energy and ζ Υ associated with the turbulence kinetic energy,

Z Υ = -Υ ij U i ∂ ∂x j ρ ρ , (43)
ζ Υ = -Υ ij u i ∂ ∂x j ρ ρ , (44)
• the interaction between the mean and fluctuating density part of total energy by conduction Z λ .

Z λ = -λ ∂T ∂x j ∂ ∂x j ρ ρ . ( 45
)
This set of equations is represented in a schematic form in figure 3.

The two-stage formulation includes the mean density turbulence kinetic energy. In the remaining part of the paper, we will establish its evolution equation in the spectral domain.

Spectral equation of the mean density turbulence kinetic energy

The investigation of the mean density turbulence kinetic energy evolution equation provides information on the energy exchanges associated with this quantity in the spatial direction. The study may be extended to the spectral domain by establishing its spectral evolution equation. This investigation permits to give the effect of the energy exchanges with regard to the size of the turbulent structures.

In order to write the spectral equation, we need to consider a flow with at least one direction of homogeneity. Without loss of generality, we consider a turbulent flow with two homogeneous and periodic directions, x and z. The inhomogeneous direction is denoted y. The dimensions of the domain in the x, y and z directions are denoted L x , L y and L z respectively. Since the flow is periodic in the x and z directions, we perform the Fourier transform in the x and z directions only and leave the y direction untransformed. Any physical quantity g(x, y, z) can be expressed as a Fourier series,

g(x, y, t) = ∞ p,q=-∞ g(k p,q , y, t)e ikp,q•x , (46)
where p and q are positive or negative integers, x = (x, z) is the position vector in the xOz plane and

k p,q = k = (k x , k z ) = ( 2πp Lx , 2πq
Lz ) is the position vector in the k x Ok z plane. The Fourier coefficients of the Fourier series expansion of g are denoted with the hat operator ( ) and are given by [START_REF] Lesieur | Turbulence in Fluids. Fluid Mechanics and Its Applications[END_REF]:

g(k, y, t) = 1 L x L z Lx 0 Lz 0 g(x, y, t)e -ik•x dx. ( 47
)
The half-trace of the velocity fluctuation correlation tensor e is equal to half the correlation C(r, y, t) = u i (r, y, t)u i (x + r, y, t)

in the particular case r = 0, that is e(y, t) = 1 2 C(0, y, t). We can express C as a Fourier series from [START_REF] Serra | High-temperature gradient effect on a turbulent channel flow using thermal large-eddy simulation in physical and spectral spaces[END_REF]. The Fourier coefficients C can be written as [START_REF] Lesieur | Turbulence in Fluids. Fluid Mechanics and Its Applications[END_REF][START_REF] Pope | Turbulent Flows[END_REF]

C(k, y, t) = u i * (k, y, t) u i (k, y, t), (49)
where the asterisk ( * ) denotes the complex conjugate. Henceforth, we denote ě and call spectral equivalent of e the half of the spectral correlation C,

ě = 1 2 u i * u i . (50)
From similar arguments, we shall associate to each term of the evolution equation of the half-trace of the velocity fluctuation correlation tensor a spectral equivalent.

Spectral equations of the turbulence kinetic energy were given by Domaradzki et al. [START_REF] Domaradzki | Energy transfer in numerically simulated wall-bounded turbulent flows[END_REF], Marati et al. [START_REF] Marati | Energy cascade and spatial fluxes in wall turbulence[END_REF], Dunn and Morrison [START_REF] Dunn | Analysis of the energy budget in turbulent channel flow using orthogonal wavelets[END_REF], Bolotnov et al. [START_REF] Bolotnov | Spectral analysis of turbulence based on the DNS of a channel flow[END_REF], Lee and Moser [START_REF] Lee | Spectral analysis on reynolds stress transport equation in high firewall-bounded turbulence[END_REF], Mizuno [START_REF] Mizuno | Spectra of energy transport in turbulent channel flows for moderate reynolds numbers[END_REF] in the incompressible case and Aulery et al. [START_REF] Aulery | Energy transfer process of anisothermal wall-bounded flows[END_REF] in the variable density case. The present decomposition gives a clear one-to-one correspondence between the terms of the spectral and spatial decompositions. To this end, a purely spectral term with no spatial contribution has to be considered. The decomposition is similar to the decomposition of Mizuno [START_REF] Mizuno | Spectra of energy transport in turbulent channel flows for moderate reynolds numbers[END_REF] at the incompressible limit.

The spectral evolution equation of the mean density turbulence kinetic energy may be written as:

∂ě ∂t = φc 0 + φΥ 0 + P0 + ζc + ζΥ + ˇ 0 , (51)
where we identify the following terms:

• the convection φc 0 , associated with the spatial convection ϕ c 0 ,

φc 0 = -Re   - 1 2 ∂ρ u i * u i U j ∂x j   , (52)
• the transfer by the total stress φΥ 0 , associated with the spatial transfer by the total stress ϕ Υ 0 ,

φΥ 0 = Re ∂ ∂x j ρ u i * Υ ij ρ , (53)
• the interaction with mixed kinetic energy P0 , associated with the spatial interaction with mixed kinetic energy P 0 ,

P0 = Re -ρ u i * u j ∂U i ∂x j , (54)
• the kinetic energy dilation correlation ζc , associated with the spatial kinetic energy dilation correlation ζ c ,

ζc = Re   1 2 u i * u i ∂ρU j -ρU j ∂x j   , (55)
• the interaction with fluctuating density kinetic energy by the total stress ζΥ , associated with the spatial interaction with fluctuating density kinetic energy by the total stress ζ Υ ,

ζΥ = Re   u i * Υ ij ∂ ∂x j ρ ρ   , (56)
• the interaction between kinetic energy and internal energy ˇ 0 , associated with the spatial interaction with internal energy 0 ,

ˇ 0 = Re   -ρ ∂u i ∂x j * Υ ij ρ   , (57)
• the purely spectral term Ξ, with no contribution in the spatial domain,

Ξ = Re   1 2 ∂ρu i ∂x j * u i u j - 1 2 u i * u j ∂ρu i ∂x j   , (58)
where Re denotes the real part operator.

The terms are decomposed in order to have a one-to-one correspondence with the terms of the spatial decomposition. The inverted hat operator (ˇ) is used to indicate the spectral equivalent of a spatial term. The spectral and spatial terms are tied closely. Given a spectral term written in the form Re( u i ' * a), for any a, the associated spatial term is u i a. The spectral term comes from the Fourier coefficients of the spatial two-point correlation between u i and a [START_REF] Lesieur | Turbulence in Fluids. Fluid Mechanics and Its Applications[END_REF].

The purely spectral term Ξ has no associated spatial term. The summation over the whole wavenumber space of its spectra is zero. In other words, this term has no effect on the spatial balance of kinetic energy but contributes to the interscale redistribution of kinetic energy.

Conclusion

The ternary decomposition of kinetic energy gives another angle of approach to the study of energy exchanges in turbulent flows. Based on the classical Reynolds averaging, the decomposition leads to the definition of a turbulence kinetic energy, a mean kinetic energy and a mixed kinetic energy. This new term, specific to the formulation, is related to the interaction between the turbulent motion and the mean motion. In the formulation, any energy exchange between the turbulence kinetic energy and the mean kinetic energy goes through the mixed term, which adds its contribution to the exchange. The formulation is decomposed further in order to include the mean density turbulence kinetic energy, product of the mean density and the half-trace of the velocity fluctuation correlation tensor. This is done by splitting the density in a mean part and a fluctuating part. Contrary to the turbulence kinetic energy, the equation of the mean density turbulence kinetic energy can easily be extended to the spectral domain. This associates a spectral equivalent to each spatial term and adds a purely spectral term to the decomposition.

Figure 1 -

 1 Figure 1 -Schematic representation of the energy exchanges between the four parts of total kinetic energy (in the one-stage formulation). An arrow represents an interaction between two quantities. The orientation of the arrow gives the sign of the term according to classical thermodynamic convention. (The term is positive in the right-hand side of the evolution equation of the quantity towards which the arrow points and negative in the other). The conservative energy transfers are represented within square brackets.

Figure 2 -

 2 Figure 2 -Schematic representation of the energy exchanges between the four parts of total kinetic energy according to the formulation of Chassaing et al. [10]. Refer to the caption of figure 1 for some indications on how to read this diagram.

Figure 3 -

 3 Figure 3 -Schematic representation of the energy exchanges between the eight parts of total kinetic energy (in the two-stage formulation). Refer to the caption of figure 1 for some indications on how to read this diagram.