Alexis Dalmon

Mathieu Lepilliez

Sébastien Tanguy

Romain Alis

Elena Roxana Popescu

Rémi Roumiguié

Thomas Miquel

Barbara Busset

Henri Bavestrello

Jean Mignot

• Mathieu Lepilliez

S´ebastien Tanguy

R´emi Roumigui´e

Sébastien Tanguy

Comparison between the FLUIDICS experiment and direct numerical simulations of fluid sloshing in spherical tanks under microgravity conditions

Keywords: Fluid sloshing, Microgravity experiment, Direct Numerical Simulation

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

Introduction

The sloshing of fluids inside tanks is one of the most important disruptions of satellite stability in microgravity conditions [START_REF] Lepilliez | Simulation numérique des ballotements d'ergols dans les réservoirs de satellites en microgravité et à faible nombre de reynolds[END_REF]. A satellite tank is filled with Alexis Dalmon alexis.dalmon@imft.fr Extended author information available on the last page of the article. liquid propellant and gas to maintain a sufficient pressure during the lifetime of the satellite. During a manoeuvre, the motion of the satellite leads to the motion of the fluids within the tanks. The centre of mass of the fluids evolves and this generates forces and torques on the satellite structure which can deteriorate the pointing performances. For satellites requiring important accuracy on the attitude, the sloshing phenomenon have to be anticipated and technical solutions to limit its effects must be considered. Several studies have been conducted on the prediction of sloshing in microgravity conditions. Analytical models, based on a 1D mass spring system or a pendulum system have been introduced by [START_REF] Abramson | The dynamic behaviour of liquids in moving containers, with applications to space vehicle technology[END_REF] and taken over by [START_REF] Dodge | The new dynamic behaviour of liquids in moving containers[END_REF]. The parameters of the simplified systems are deduced from the tank geometry, the fluids properties and the accelerations exerted on the tank. Finally, the model allows to compute the sloshing frequencies, damping ratios and forces exerted on the tank wall for a given manoeuvre. However, these analytical models cannot be considered within the framework of low Bond numbers manoeuvres for which the non linear effects of the surface tension cannot be neglected. Very few experimental data exist because of the complexity to maintain microgravity conditions for a long period of time. Parabolic flights and drop towers enable to create several seconds of microgravity environment [START_REF] Fernandez | The cfvib experiment: control of fluid in microgravity with vibrations[END_REF][START_REF] Li | Dynamic behaviors of liquid in papartial filled tank in short-term microgravity[END_REF]) but low Bond numbers manoeuvres require several minutes for the fluids to stabilise. Finally, Direct Numerical Simulation is a promising candidate to investigate the fluids behaviour during station-keeping manoeuvres. The capillary effects and the microgravity conditions can be integrated in the numerical models without additional complexity compared to Earth gravity conditions. Numerical studies on the liquid sloshing in tanks have been performed in different configurations [START_REF] Ebrahimian | Free vibration sloshing analysis in axisymmetric baffled containers under low-gravity condition[END_REF][START_REF] Liu | A numerical study of three-dimensional liquid sloshing in tanks[END_REF][START_REF] Tang | Simulation of large-amplitude three dimensional liquid sloshing in spherical tanks[END_REF][START_REF] Tang | Improved method for implementing contact angle condition in simulation of liquid sloshing under microgravity[END_REF][START_REF] Veldman | The numerical simulation of liquid sloshing on board spacecraft[END_REF][START_REF] Yue | Large-scale amplitude liquid sloshing in container under pitching excitation[END_REF]). Nevertheless, experimental data are still needed to validate the numerical results. To provide experimental data on this topic, the French national space agency (CNES) sent an experiment to the International Space Station (ISS) to quantify the sloshing phenomenon in microgravity conditions [START_REF] Dalmon | Direct numerical simulation of a bubble motion in a spherical tank under external forces and microgravity conditions[END_REF][START_REF] Mignot | Fluid dynamic in space experiment[END_REF]. Our objective in this paper is to carry out comparisons between numerical simulations and data collected in space to validate accurately the numerical methods predicting fluid sloshing. In the "Tank Modelling" section, the description of the study case will be presented and the dominating effects will be highlighted by using the dimensionless Navier-Stokes equations. We will consider the motion of the fluids inside a tank rotating around a staggered axis in microgravity conditions. Because satellite tanks are usually not aligned with the inertial axis of the satellite, this study aims at modelling the fluids behaviour during a typical rotational manoeuvre. Then, the experimental setup will be described in the "The FLUIDICS Experiment" section with the different measuring devices. The "Description of the Numerical Tool" section details the numerical methods implemented to predict the fluids behaviour inside the spherical tank. Finally, the "Comparisons" section focuses on the comparisons between the experimental data and the numerical simulations. The temporal evolution of the force and torque exerted by the fluids on the tank wall will be investigated for different angular velocity profiles.

Tank Modelling

Dimensionless Navier-Stokes Equations

In this study, we consider the rotation of a spherical tank around the z-axis in microgravity conditions. Figure 1 presents a schematic of the study case. The tank is a sphere of diameter D t , the centre of which is labelled C. The fluid domain inside the tank is denoted by f and is displayed in blue in Fig. 1. The tank is linked to the rotation axis (located at point O) by a lever arm of length L. The direction of the lever arm defines the radial axis called e y in the tank referential. The rotation happens around the z-axis and is characterised by its angular velocity and angular acceleration ˙ .

To quantify the main effects of the manoeuvre on the fluids, we write the Navier-Stokes equations in the fluid domain

∇ • u = 0, ( 1
)
ρ ∂u ∂t + (u • ∇) u = -∇p + ∇ • 2μ D + F ext , ( 2
)
where t is the time, ρ the fluid density, μ the fluid viscosity, u = (u, v, w) the velocity field, p the pressure, F ext the is the rate-of-deformation tensor defined as

D = ∇u + ∇u T 2 . ( 3
)
The referential of the tank is non-Galilean. The inertial and Coriolis accelerations must be taken into account through the external volume force. The Coriolis acceleration corresponds to a 1 = 2 × u and the inertial acceleration is the sum of the Euler acceleration a 2 = ˙ × OM, the centrifugal acceleration a 3 = × ( × OM) and the relative acceleration a R between the tank and the rotating referential at point O. Here, the relative acceleration is zero and the position vector of the fluid particle with respect to the rotation axis is defined as

OM = OC + CM = xe x + (L + y)e y + ze z , ( 4
)
with CM = (x, y, z) the coordinates of a fluid particle in the referential of the tank. The total volume force in the tank reference F ext is expressed as follows,

F ext = -ρ i a i = ρ ⎛ ⎝ 2 v + (L + y) ˙ + x 2 -2 u -x ˙ + (L + y) 2 0 ⎞ ⎠ .
(5)

The fluid domain f is now divided in two parts: l the liquid propellant region and g the gas region. The filling ratio τ f ill represents the volume of liquid divided by the total volume of the tank,

τ f ill = V l V f . ( 6
)
The physical properties of each fluid are respectively denoted by the subscript l or g. We call the interface between the fluids regions (coloured in blue in Fig. 2) and n its outward normal. From now on, the jump of the physical properties and the effect of surface tension σ are considered at the interface. Moreover, the liquid propellants used in the space industry are perfectly wetting fluids. This forces the liquid to stay in contact with the tank wall. Therefore we consider that the gas initially takes the form of a spherical bubble of diameter D 0 centred in the tank. In Fig. 2, the interface between the fluids is represented in blue, the gas Fig. 2 Inside view of the tank at the initial time occupies the centre of the tank and the liquid is in contact with the tank wall.

Some dimensionless variables are used to extract characteristic numbers from the Navier-Stokes equations,

ρ * = ρ ρ l , μ * = μ μ l , ∇ * = D 0 ∇, t * = t σ ρ l D 3 0 , (7) u * = u ρ l D 0 σ , p * = pD 0 σ , a * i = a i γ , ( 8
)
with γ a characteristic acceleration.

Considering the Coriolis and the centrifugal accelerations, the characteristic acceleration is based on the angular velocity and is equal to 2 L. For the Euler acceleration, the characteristic acceleration is based on the angular acceleration, the value of which is ˙ L. Finally, the dimensionless Navier-Stokes equations can be written as,

∇ * • u * = 0, ( 9
)
ρ * ∂u * ∂t * + u * • ∇ * u * = -∇ * p * + Oh∇ * • 2μ * D * -Bo i ρ * a * 2 -Bo c ρ * (a * 1 +a * 3 ), ( 10
)
with Oh the Ohnesorge number and two Bond numbers defined as,

Oh = μ 2 l ρ l σ D 0 , Bo c = ρ l 2 LD 2 0 σ , Bo i = ρ l ˙ LD 2 0 σ . ( 11
)
The Ohnesorge number compares the viscosity to the surface tension effects and the Bond numbers compare the inertial effects induced by the manoeuvre to the surface tension effects. Bo c is the Bond number that focuses on the inertial effects related to the angular velocity and Bo i is the Bond number that depends on the angular acceleration.

Configuration of the Manoeuvres

The manoeuvre we consider is the typical rotation manoeuvre of satellites in orbit. As defined previously, the rotation of the tank happens around the z-axis and the angular acceleration and velocity are defined as ˙ = ˙ e z and = e z . The typical angular velocity profiles we consider are shown in Fig. 3 with the corresponding Bond numbers. They consist in a first phase with a constant angular acceleration ˙ until the angular velocity reaches its required value . From then on, the angular acceleration is set to zero and the velocity stays at its constant value. Finally, the opposite angular acceleration is enforced until the angular velocity comes back to zero and the tank remains fixed. In this study, two different filling ratios are considered: 50% and 75%. Several constant angular velocities are enforced on the tanks. They are defined with respect to the angular velocity of the first manoeuvre, denoted by 1 , and are summarised in Table 1. The angular accelerations applied during the first and third phases with opposite signs are characterised by the duration of these phases denoted by t acc with respect to the time of the first manoeuvre t 1 acc . The corresponding Bond numbers, defined in the "Dimensionless Navier-Stokes Equations" section, are computed for the 6 manoeuvres of this study and are listed in Table 1.

The FLUIDICS Experiment

The FLUIDICS (FLUId DynamICs in Space) experiment was a part of the Proxima mission of the CNES, the French space agency, lead by the astronaut of the ESA (European Space Agency) Thomas Pesquet in the International Space Station (ISS). The FLUIDICS experiment, represented in Fig. 4, is basically a slow rate centrifuge able to reproduce the fluids motion inside a satellite tank in orbit. It consists of a spherical tank containing a safe substitute of propellant and air, connected by an arm to a motor. Two spherical tanks with the same diameter D t have been used in microgravity conditions. They are made of two polycarbonate half spheres glued together. These tanks are respectively filled with 50% and 75% of 3M Novec 2704 liquid. Figure 5 is a photograph taken on earth of the tank half filled with orange Novec liquid. The liquid is an electronic card coating which has the same range of physical properties (viscosity μ l and density ρ l ) as propellants used in the space industry and complies with the ISS safety guidelines. The perfectly wetting property of these propellants is still respected with the Novec fluid. This will enable us to consider a zero contact angle condition in the simulations presented hereafter. The other fluid contained in these tanks is air, the physical properties of which are denoted by μ g and ρ g . The surface tension between the liquid and the gas is symbolised by σ .

Each tank can be plugged at the end of the rotating arm through an interface linked to a force sensor. This sensor gives the forces and torques exerted by the tank on the Fig. 5 Tank half filled with Novec liquid on earth Fig. 6 The FLUIDICS experiment fixed to the ISS structure rotating arm in the three directions. The sensor is located at a distance L se from the centre of the tank towards the rotating arm. The sensor location is depicted in Fig. 1 at point S. The rotating arm is linked to the motor located at the centre of the experiment. The lever arm L is the distance between the rotation axis and the centre of the tank. Beside the tank, the rotating arm supports two high resolution monochromatic cameras pointing at the tank location and, at the other end, the electronics needed to record data. Finally, the experiment is enclosed in a protection box to comply with safety rules and is fixed to the ISS seat tracks as shown in Fig. 6. A ventilation system cools down the inside of the protection box because of the heat released by the electronics. More details about the experiment are given in [START_REF] Mignot | Fluid dynamic in space experiment[END_REF].

The 6 manoeuvres described in Table 1 have been executed in the ISS in 2017. The manoeuvres with the halffilled tank have been done twice to ensure repeatability. The data have been post-processed to be correlated with numerical simulations. The data from the force and torque sensor have been filtered and the solid contribution of the empty tank in the measured data has been removed to quantify only the effect of the fluids on the tank wall.

Description of the Numerical Tool

Numerical Methods

In this section, the numerical solver considered in this study is presented. The home-made code DIVA (Dynamics of Interface for Vaporization and Atomization) is based on several numerical methods dedicated to the computation of two-phase flows.

The computational domain can be divided in two regions: a fluid region f corresponding to the inside of the tank and a solid region s which consists of the rest of the cubic computational domain. The interface between these regions is denoted by s and represents the tank wall. As seen in the "Dimensionless Navier-Stokes Equations" section, the fluid region is divided in two phases as well: the liquid region l and the gas region g separated by the interface denoted by .

The level-set method of [START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on hamilton-jacobi formulations[END_REF] and [START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flow[END_REF] is used to follow the interface motion of through time and to define the static position of s . Two level-set functions are defined and represent first, the signed distance function from the liquid gas interface φ and second, the signed distance from the solid fluid interface φ s . The motion of the liquid-gas interface is captured by solving the following transport equation

∂φ ∂t + u • ∇φ = 0, ( 12
)
coupled with the reinitialization algorithm of [START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flow[END_REF] to maintain the distance property of the level-set function through time. Based on the work of [START_REF] Fedkiw | A non-oscillatory eulerian approach to interfaces in multimaterial flows (the ghost fluid method)[END_REF] and [START_REF] Aslam | A partial differential equation approach to multidimensional extrapolation[END_REF], the following boundary condition on the level-set function φ must be considered within the solid domain to enforce the contact angle θ s if a triple line exists

n s • ∇φ| = -cos θ s , (13)
with n s the outward normal vector of the fluid-solid interface. As previously stated, the simulations involved in this paper will still consider a zero contact angle condition since the spherical tanks are filled with a perfectly wetting fluid. This means that there is always some liquid against the tank wall in our simulations. The zero contact angle hypothesis will be further justified by comparing the numerical results with the experimental data, especially the torque signal which depends on the interface position.

The incompressible Navier-Stokes equations for Newtonian two-phase flows are solved using the Ghost Fluid Conservative viscous Method with an Implicit scheme (GFCMI) of [START_REF] Lepilliez | On two-phase flow solvers in irregular domains with contact line[END_REF]. This method is based on the projection method of [START_REF] Chorin | A numerical method for solving incompressible viscous flow problems[END_REF] first introduced by [START_REF] Sussman | A sharp interface method for incompressible two-phase flows[END_REF]. It consists in three steps: first, considering the velocity field u n at time t n = n t, an intermediate velocity u * is computed without considering the pressure term by solving the following linear system,

ρ n+1 u * -t∇ • 2μ n+1 D * = ρ n+1 u n -t u n • ∇ u n - F ext ρ n+1 . ( 14
)
The rate-of-deformation tensor is solved implicitly to remove the time-step restriction due to the viscosity. The linear system which couples the three components of the velocity is solved using a Gauss-Seidel algorithm.

Then, the pressure field p n+1 is regarded as the scalar potential function of the Hodge decomposition and satisfies the Poisson equation with Neumann boundary conditions

∇ • ∇p n+1 ρ n+1 = ∇ • u * t + ∇ • σ κnδ ρ n+1 , ( 15
)
with κ the local curvature of the liquid-gas interface and n the outward normal vector of the interface and δ a multidimensional Dirac distribution at the interface.

The density jump and the pressure jump due to the surface tension term are taken into account with the sharp representation of the interface from [START_REF] Liu | A boundary condition capturing method for poisson's equation on irregular domain[END_REF]. The resolution of the linear system is performed using the Black Box MultiGrid method, see [START_REF] Dendy | Black box multigrid[END_REF] and MacLachlan et al. ( 2008) for more details.

Finally, the velocity at the next time step u n+1 is computed as the projection of the intermediate velocity onto the divergence-free space

u n+1 = u * - t ρ n+1 (∇p n+1 -σ κnδ ). ( 16
)
The no-slip boundary condition is considered at the tank wall within the computational domain. A subcell Dirichlet boundary condition on the velocity (zero velocity at the tank wall) and a Neumann boundary condition on the pressure (zero normal pressure gradient) are enforced using the methodology of [START_REF] Gibou | A second-orderaccurate symmetric discretization of the poisson equation on irregular domains[END_REF], [START_REF] Ng | An efficient fluid-solid coupling algorithm for single-phase flows[END_REF], and [START_REF] Lepilliez | On two-phase flow solvers in irregular domains with contact line[END_REF].

Following the general guidelines of the Ghost Fluid Method, our home-made code DIVA can also consider phase change simulations, by imposing additional jump conditions to maintain mass conservation [START_REF] Nguyen | A boundary condition capturing method for incompressible flame discontinuities[END_REF], energy conservation [START_REF] Gibou | A level set based sharp interface method for the multiphase incompressible navierstokes equations with phase change[END_REF][START_REF] Tanguy | Benchmarks and numerical methods for the simulation of boiling flows[END_REF]) and chemical species conservation [START_REF] Papac | Efficient symmetric discretization for the poisson, heat and stefan-type problems with robin boundary conditions[END_REF][START_REF] Tanguy | A level set method for vaporizing two-phase flows[END_REF][START_REF] Villegas | A ghost fluid/level set method for boiling flows and liquid evaporation: application to the leidenfrost effect[END_REF], and compressible flows [START_REF] Huber | A time splitting projection scheme for compressible two-phase flows. application to the interaction of bubbles with ultrasound waves[END_REF]. This solver has been extensively validated with other experimental data both for isothermal two-phase flows [START_REF] Lalanne | Non-linear shape oscillations of rising drops and bubbles: experiments and simulations[END_REF][START_REF] Tanguy | Application of a level set method for simulation of droplet collisions[END_REF] and for two-phase flows with phase change for impacting droplets in the Leidenfrost regime (Rueda [START_REF] Villegas | Direct numerical simulation of the impact of a droplet onto a hot surface above the leidenfrost temperature[END_REF]) and for nucleate pool boiling [START_REF] Huber | Direct numerical simulation of nucleate pool boiling at large microscopic contact angle and moderate Jakob number[END_REF][START_REF] Urbano | Direct numerical simulation of nucleate boiling in micro-layer regime[END_REF].

Mesh Convergence

To ensure the numerical simulations are well defined, a mesh convergence study on a specific manoeuvre is presented. We consider manoeuvre n • 1 (see Table 1) which corresponds to the highest values of the Bond numbers: Bo c = 76 and Bo i = 23 with the half-filled tank. The maximum angular velocity is maintained during the second phase and the filling ratio of 50% ensures an important motion of the gas bubble within the tank. The simulations have been done on three different meshes with 64, 128 and 256 cells in each direction. The force F and the torque at the sensor location T exerted by the fluids on the tank wall are computed in the simulations as

F = f (∇p -σ κnδ )dV , (17) T = f SM × (∇p -σ κnδ )dV . ( 18
)
with SM the vector between the sensor and a fluid particle as SM = xe x + (L se + y)e y + ze z . As the Reynolds number based on the velocity 2 L,

Re = ρ l 2 LD t μ l , ( 19
)
varies between 12500 and 37500 for manoeuvres n • 1 to n • 4, only the pressure component of the force and the torque exerted on the tank wall has been considered. Indeed, for these high values, the viscous component of the force can be neglected with respect to the pressure component.

Figures 7 and8 respectively show the evolution of the dimensionless force in the y-direction and the dimensionless torque in the z-direction until t * = 3 for the three meshes. The physics corresponding to the temporal evolution of the force and the torque will be described in the "Comparisons", this section only takes interest in the grid sensitivity.

The temporal evolutions of the force and the torque are similar for the three meshes. At the beginning, the increase of the force toward the same peak is predicted by all the meshes. The same oscillatory behaviour is obtained with similar frequencies and magnitudes. The mean value of the force after t * = 0.5 varies a little but stays close for the two finest grids. Similarly, the frequency and the magnitude of the torque oscillations for the three meshes are really close to each other. The meshes containing 128 and 258 cells in each directions are in good agreement and prove the spatial convergence. We will then consider the intermediate mesh 128 3 for all the simulations of this paper.

Comparisons

This section aims at comparing the data from the FLUIDICS experiment with the numerical results for the manoeuvres listed in Table 1. The comparison deals with the y-direction force and the z-direction torque. Indeed, because the system is symmetrical to the plane {z = 0}, no external acceleration happens in the z direction leading to no force nor motion of the centre of mass toward e z . Thus, the relation between the force and the torque at the sensor location becomes

T = SG × F = x g F y -(L se + y g )F x e z , ( 20
)
with SG = SC + CG and CG = (x g , y g , z g ≈ 0) the coordinates of the centre of mass in the tank referential. Because L se is higher than the maximum coordinates of the centre of mass, the torque in the z-direction has the same temporal evolution than the x-component of the force. To avoid redundant comparisons and because the sensor has a greater sensitivity to the torque than to the force, the x-direction force is not investigated at the expense of the z-direction torque.

Force in the Y-Direction

The force in the y-direction corresponds to the highest component of the forces exerted by the fluids on the tank wall. During the constant angular velocity phase, the centrifugal acceleration is the only external force that remains after the fluids stabilisation. Its asymptotic value is thus a relevant criterion to compare the experimental data with the simulations. Figure 9 shows the dimensionless force evolution measured in the ISS (in solid black) and computed by the code (in dashed blue) for manoeuvre n • 1. First, the dimensionless force increases considerably during the acceleration phase. Once the first phase is over, the angular velocity remains constant and the angular acceleration is set to zero. The volume force is then dominated by the centrifugal force and becomes

F ext = ρ ⎛ ⎝ 2 v + x 2 -2 u + (L + y) 2 0 ⎞ ⎠ . ( 21
)
Consequently, in the second phase, the liquid moves to the direction +e y because of the centrifugal force and the gas bubble spreads against the tank wall in the opposite direction. At the beginning of the second phase, the fluids are not stabilised and the gas bubble oscillates around the xaxis {x = 0}. This is due to the sign of the x-component of the centrifugal force which changes every time the centre of mass crosses the axis {x = 0}. Four pictures of an oscillation of the gas bubble predicted by the numerical code are shown in Fig. 10. These oscillations are damped with time until the bubble reaches its equilibrium position. At the steady state, the Coriolis contribution of Eq. 21 disappears and only the centrifugal force remains.

Similarly, we observe in Fig. 9 that after t acc , the force oscillates around a constant value. The magnitude of the oscillations decreases throughout time until they disappear in the simulations (in dashed blue), and until the magnitude of the oscillations reaches the noise level in the experimental data (in solid black).

Figure 11 depicts the gas bubble (in dark blue) inside the spherical tank at the equilibrium state according to the simulations. The gas bubble is spread against the tank wall Fig. 9 Comparison of the measured and computed forces in the y-direction throughout time for manoeuvre n • 1 Fig. 10 Oscillation of a gas bubble (in blue) in the spherical tank predicted with τ f ill = 50% in the direction -e y because of the centrifugal force. It is stabilised and symmetrical to the axis {x = 0}, leading to no force in the x-direction nor torque in the z-direction and the x-coordinate of the centre of mass equals zero. Because of the perfectly wetting condition, a thin film of liquid is trapped between the bubble and the tank wall, as can be observed in Fig. 11.

Finally, the angular velocity decreases during the deceleration phase, leading to the drop of the y-direction force until it reaches zero when no more external accelerations are enforced.

The two curves in Fig. 9 exhibit the same temporal evolution of the y-direction force. The mean values of the measured and computed forces during the second phase are close to each other. The measured mean value is 34.84 and the computed mean value is 34.16 leading to a difference lower than 2%. This manoeuvre has been enforced a second time in the ISS and the resulting mean value is 34.69 which gives a difference still lower than 2% between the numerical simulations and the experimental data. We observe as well that the first oscillations are well predicted by the simulations. The computed force coincides with the data for the three first oscillations. Over time, the magnitude of the oscillations decreases and the variation of the measured force reaches the level of noise. In the numerical simulations, the oscillations carry on until the end Fig. 11 Gas bubble at the steady state of the second phase of the first manoeuvre and zoom on the film of liquid at the tank wall of the second phase with an extremely low magnitude. The more detailed study of the oscillations is presented in the next subsection through the torque evolution.

Figure 12 shows the evolution of the measured and computed dimensionless force in the y-direction throughout time for manoeuvre n • 2. It consists in the same velocity profile as the first manoeuvre with the tank containing 75% of liquid. The force evolution corresponds to the same behaviour but the force value during the second phase is more important. The filling ratio is higher, which means that more liquid is inside the tank and leads to a higher centrifugal acceleration. The computation of the mean forces shows a slightly larger difference between the data and the simulations compared to the previous manoeuvre. The mean value of the measured force is 51.29 and the mean value of the computed force is 49.31 which leads to an error of 3.85%.

The damping of the oscillations is stronger than in the previous case. They disappear in the numerical simulation after t * = 4. This phenomenon is a consequence of the high filling ratio. The gas bubble is smaller and spreads with a higher magnitude on the tank wall. Because of the tank geometry, the gas bubble is highly deformed and its motion is more constrained. The low variation of the force at the beginning of the second phase prevents us from extracting more information about the oscillations.

Considering a lower angular velocity, the force evolution measured and predicted for manoeuvre n • 3 is plotted in Fig. 13. The average values of the force during the second phase stay close for the two experimental runs between 11.35 and 11.72. The comparison with the computed value of 11.37 shows a difference lower than 3% for the two runs of manoeuvre n • 3 carried out in the ISS. However, the oscillations of the bubble cannot be extracted from the experimental data because the noise level is too important with respect to the centrifugal force. With the simulation, the bubble oscillations are visible but their magnitudes do not exceed the experimental noise level. Whereas the fluid stabilisation may be questionable in Fig. 13, studying the videos helped us to conclude that a steady state has been reached at the end of the constant velocity phase. As a result, the fluctuations of the signals are only due to measuring uncertainties. One can notice that the noise amplitude remains approximately constant whatever the considered manoeuvre, even if it seems more important for lower Bond numbers.

The same velocity profile enforced on the 75%-filled tank gives the same observations. Figure 14 shows that the bubble oscillations are more damped because of the additional mass of liquid and their magnitudes are still lower than the noise level. The average value of the force increases from values around 11.5 to 16 between the two tanks. Finally, the comparison of the mean values of the force shows a relative difference of 6.08%, which is the highest error we observe in the manoeuvres.

Table 2 gathers all the average values of the measured and computed dimensionless forces for manoeuvres n • 1 to 4. Manoeuvres n • 1 and 3 have been executed twice in the ISS and are respectively denoted with the subscripts 'a' and 'b'. The angular velocities of manoeuvres n • 5 and 6 are too low to extract reliable data from the force and torque sensor and are not considered in this section.

In brief, the level of force increases with the Bond number Bo c and the filling ratio. This is explained by the centrifugal force which depends directly on the square of the angular velocity and on the volume of liquid inside the tank. For the highest Bond numbers manoeuvres, the temporal evolution of the force is well predicted by the simulations: the mean force in the second phase and the oscillations are in good agreement with the experimental data. For all the manoeuvres, the numerical simulations predict mean values of forces with a maximum difference of 6%, providing a strong validation of the proposed numerical approach.

Torque in the Z-Direction

The torque in the z direction allows us to extract other contributions of the fluids on the tank wall. In addition to the y-direction force, the z-component of the torque exhibits the effects of the angular acceleration on the fluids and emphasises the oscillations of the bubble around its equilibrium position. Considering manoeuvre n • 1, the dimensionless torque obtained by the FLUIDICS experiment in a microgravity environment and predicted by the simulations are plotted in Fig. 15 in solid black and dashed blue, respectively.

During the acceleration phase, we observe an important negative peak of the torque. This is due to the Euler acceleration which dominates the flow at the beginning of the manoeuvre. At early flow time, the angular velocity is low but the angular acceleration is already enforced on the

F ext = ρ ˙ ⎛ ⎝ L + y -x 0 ⎞ ⎠ . ( 22
)
Because the lever arm L is higher than the maximal centre of mass coordinates, the volume force is dominated by the contribution of the Euler acceleration in the xdirection. This contribution is recovered in the z-direction torque following Eq. 20 and explains the negative peak we have measured. The same phenomenon appears during the deceleration phase with a positive peak of the z-direction torque (visible at t * = 8.2 in Fig. 15) because the angular acceleration is the opposite. During the constant angular velocity phase, the oscillations due to the bubble motion are clearly visible. The torque oscillates around zero and its magnitude decreases with time until the bubble reaches its equilibrium position, centred at the axis {x = 0}. Indeed, at the steady state, the x-coordinate of the centre of mass is zero and the volume force is only directed toward the y-direction (see Eq. 21).

The comparison between the two curves of Fig. 15 shows a very good agreement on the temporal evolution of the torque measured in microgravity conditions and predicted by numerical simulations. The peaks corresponding to the acceleration and deceleration phases have the same magnitude in both cases and the sharp increase and decrease of the torque is well respected. The oscillations of the computed torque concur with the experimental data until the magnitude of the latter reaches the noise level. Both the magnitude and the frequency of the measured and predicted torque are in good agreement. The frequency of the oscillations is the aim of the next subsection. Figure 16 represents the time evolution of the dimensionless torque for manoeuvre n • 2. The same angular acceleration and velocity as for manoeuvre n • 1 are enforced on the tank with τ f ill = 75%. Similarly, two peaks of opposite signs are observed during the first and third phases. Their magnitudes are increased with the filling ratio. During the constant angular velocity phase, the observed oscillations are more damped than for the half-filled tank. However, the first predicted oscillations have the same magnitude and frequency as the measured torque but the important damping and the data noise prevent us from going further.

The temporal evolutions of the torque for manoeuvres n • 3 and n • 4 are depicted in Figs. 17 and 18, respectively. The magnitude of the torques decreases because of the lower angular velocity. The peaks corresponding to the contribution of the Euler acceleration are still visible and well predicted by the simulations but the angular acceleration being lower, the magnitudes of the peaks decrease. As the magnitude of the torque oscillations also decreases, the noise from the experimental setup becomes too large, preventing a possible comparison of the amplitude and the frequencies of the torque oscillations during the second phase. These comparisons between the numerical and experimental torque signals provide also an indirect validation of the zero contact angle hypothesis. Indeed, as the torque depends both on the force and on the interface position, the good agreement reported previously confirms the relevance of this assumption, since the interface position depends strongly on the contact angle value.

Sloshing Frequency

As previously stated, the oscillations observed at the beginning of the second phase are due to the motion of the gas bubble inside the tank before reaching its equilibrium shape and position depicted in Fig. 11. We consider the oscillations of the torque in the z-direction and apply a discrete Fourier transform to obtain the amplitude spectra of the signals. The transform is applied on the same time range t * = 7.5 for all the manoeuvres. The objective of this section is to extract the frequencies of the oscillations for each manoeuvre and to compare the predicted value to the data from FLUIDICS.

Figure 19 depicts the amplitude spectra of the z-direction torque for manoeuvre n • 1. In dashed blue, the amplitude spectrum of the torque from the simulation exhibits a main peak at f * = 2.67. No other characteristic frequency appears in the simulation. In solid black, the amplitude spectrum from the experimental data is more complex. The same peak appears for the same frequency f * = 2.67. These experimental and numerical frequencies match well within the frequency resolution f * = 0.133. For this first manoeuvre, the simulation provides an accurate prediction of the sloshing frequency. A second important peak appears for a lower frequency in the experimental data. It corresponds to the rotation frequency denoted by f = /2π . Knowing that D t /L = 0.57, the dimensionless rotation frequency can be expressed as

f * = f ρ l D 3 t σ = 1 2π D t Bo c L(1 -τ f ill ) 2/3 . ( 23
)
For manoeuvres n • 1 and 2, the dimensionless rotation frequency equals f * = 1.322 and corresponds to the lower frequency peak observed in Fig. 19. This parasitic frequency is likely due to an unbalance of the experiment and explains partly the spurious oscillations observed in the signal torque in Figs. 15 and16. Other minor peaks appear for higher frequencies and can be attributed to the other kind of perturbations such as the electronic devices, the fans, the structural vibrations... The magnitude of the peaks allows us to compare the effects of the different sources of the oscillations averaged over the constant angular velocity phase.

Figure 20 depicts the amplitude spectra of the vertical torque for manoeuvre n • 2. In dashed blue, the spectrum of the simulation shows a main frequency around f * = 2.93. In solid black, the range of frequencies [2.64; 2.94] have important magnitudes. We assume that this interval contains the sloshing frequency and some disturbing frequencies. The peak corresponding to the same rotation frequency as for the previous manoeuvre is the most important in Fig. 20. In this case, the sloshing oscillation is much more damped compared to manoeuvre n • 1 and then, the rotation frequency becomes dominant during a longer time period. This explains why the peaks corresponding to the rotation frequency exceed the sloshing oscillation frequencies. Finally, the same minor frequencies as in Fig. 19 appear around f * = 4 and f * = 5.33. As seen in Fig. 16, the torque oscillations are much more damped with the filling ratio of 75%. The interval of time for which the bubble oscillations are noticeable is reduced and the noise due to the unbalance has a greater importance in the amplitude spectrum. This explains why the only frequency for which the amplitude spectrum decreases between the first and second manoeuvres is the sloshing frequency.

Figure 21 illustrates the amplitude spectra for the third manoeuvre. The oscillation frequency is clearly visible at f * = 1.47 in the numerical results. In the experimental data, several peaks appear, represented in solid black. The lower frequency corresponds to the rotation frequency, which can be computed with Eq. 23 as f * = 0.765. The value of the predicted oscillation frequency is approached by a amplitude peak at f * = 1.6. Then, several higher frequencies appear and correspond to the noise observed in Fig. 17. Finally, considering the most filled tank, the same amplitude spectrum as in manoeuvre n • 3 is obtained with a decrease of the amplitude corresponding to the sloshing frequency. Despite the high level of noise involved in the signal in Fig. 18, the amplitude spectrum allows recovering the sloshing frequency even for this manoeuvre for which the highest perturbations are observed.

We observe that in the four amplitude spectra of Figs. 19 to 22, the same frequency f * = 5.33 appears in the experimental data. The angular velocities and the filling ratios being different for each one of these cases, this disturbance frequency is independent from the manoeuvre enforced. This may come from electronic noises, the ventilation system of the experiment or the ISS environment.

Considering the first manoeuvre, the frequencies corresponding to the rotation frequency and the constant frequency f * = 5.33 have been filtered in the frequency domain. The sources of the other frequencies observed in the amplitude spectra are not clearly identified and thus have not been filtered. After applying the inverse Fourier transform, the resulting torque signal is compared to the numerical results and presented in Fig. 23. We observe that, at early flow time, the filtered signal is quite similar to the signal in Fig. 15 with an agreement with the numerical results until t * = 3.

Table 3 summarises the oscillation frequency extracted from the amplitude spectra. The predicted and measured frequencies are in close agreement for manoeuvres n • 1 to 4. Moreover, the oscillation frequencies measured for the two runs of manoeuvres 1 and 3 are identical and listed in Table 3. This demonstrates the repeatability of the experiment and provides a validation of the simulations for For the lower Bond numbers manoeuvres (n • 5 and 6), the noise in the experimental data prevents us from extracting the force and torque temporal evolution. Nevertheless, the videos of the tank allow to observe the bubble motion throughout time. Considering the half-filled tank, the bubble oscillations are clearly visible during the constant velocity phase, even for Bo c around 6. For these ranges of Bond numbers, several oscillations (between 3 for manoeuvre n • 6 to 11 for manoeuvre n • 1) are visible during the constant velocity phase and allow to extract the characteristic time of an oscillation. For the tank with τ f ill = 75%, the oscillations are quickly damped and we cannot extract the frequency from the video.

Figure 24 represents six snapshots of an oscillation of the gas bubble inside the half-filled tank. The spherical tank is centred in the pictures and the ring corresponding to the region where the two half-sphere of the tank are glued together is recognisable. The gas bubble is located towards the lever arm direction at the bottom of the pictures such as in Fig. 11. We observe that during the oscillation, the bubble approaches the tank wall in the camera direction and so, rises in the picture until it reaches its maximum (a to d). Then, the bubble goes to the opposite direction and the interface goes down until it reaches its minimum location (e and f). Table 4 summarises the oscillation frequencies measured with the videos for all the Bond numbers exerted on the halffilled tank. They are compared with the values obtained with the Fast Fourier transform of the numerical simulations.

First, considering manoeuvres n • 1 and 3, the frequencies measured with the videos are close to these from the torque signals (see Table 3). Then, the frequencies of the simulations present a good agreement with the measured ones even for the lower Bond numbers Bo c of 6 and 13. The frequency resolution with the simulations is still equal to f = 0.133 and the accuracy when measuring the time period with the video depends on the manoeuvre. They are summarised in Table 3. Despite the margin of error, the predicted frequencies are close to the measured ones in all cases with the half-filled tank. This extends the use of the numerical code for low Bond numbers manoeuvres.

Figure 25 summarises all the oscillation frequencies obtained with the half-filled tank. As we can assume that Bo i influences weakly the oscillation phase, we take interest only in the centrifugal Bond number. For these four manoeuvres, we observe clearly that the oscillation frequency increases with Bo c .

Conclusion

This paper presented the comparison between experimental data and numerical simulations of the fluids behaviour inside a spherical tank rotating around a fixed axis in microgravity conditions. The FLUIDICS experiment in the ISS allowed to measure the force and torque exerted by the fluids on the tank wall for different manoeuvres. These data have been compared to the numerical results from our home-made code DIVA. Preliminarily, the numerical simulations have been validated through a mesh convergence study on the manoeuvre with the highest Bond numbers. Then, results on six different manoeuvres corresponding to different angular velocities and filling ratio are presented both from numerical simulations and experimental data carried out in the ISS.

The comparison of the y-direction force shows a good agreement between the measured and predicted forces. The average value of the centrifugal force is well predicted with an error lower than 6% for all the manoeuvres. This error decreases to 3% for the manoeuvres exerted on the halffilled tank. Moreover, the variations of the force at the beginning of the constant angular velocity phase are similar while their magnitudes are higher than the noise level. Considering the torque evolution, the peaks corresponding to the Euler acceleration contribution are clearly visible and their magnitudes are well predicted by the simulations. The oscillations of the gas bubble during the second phase of the manoeuvre have different magnitudes and damping ratios depending on the angular velocity and filling ratio. For the first two manoeuvres, the oscillations exceed the noise level and the measured and predicted temporal evolutions are in close agreement. By applying the discrete Fourier transform on the second phase of the torque signal, it has been possible to extract their main frequencies. The numerical simulations exhibit only one frequency corresponding to the bubble oscillation frequency. For the experimental data, further disturbing frequencies appear, such as the rotation frequency, which bring out some unbalance of the experimental setup. Higher frequencies may also come from vibrations of the structure, the electronic noises, the cooling system of the experiment or the ISS environment. However, the bubble oscillation frequencies are well predicted by the simulations within the frequency resolution of the data. In the case of the half-filled tank, the videos from the FLUIDICS experiment allowed us to extract the oscillation frequency even for the lower value of the Bond number manoeuvres for which the force and torque signals were too noisy to be used. Similarly, the comparisons with the predicted frequencies show a good agreement within a range of 4%.

The evolution of the force and torque with respect to the Bond numbers and filling ratios present a coherent set of data whether numerical or experimental results are considered. Increasing the angular velocity and the amount of liquid within the tank raises the centrifugal force exerted on the tank wall. Similarly, the damping of the bubble oscillations increases with the filling ratio, which is visible in both the force and torque evolutions. The magnitude of the torque peaks corresponding to the acceleration and deceleration phases depends on the angular acceleration and thus, on Bo i . The frequency of the bubble oscillations increases with the angular velocity.

In conclusion, the FLUIDICS experiment provided original data on the fluid sloshing in spherical tanks submitted to rotation in micro-gravity conditions. The comparison with the numerical results validates the use of two-phase flows Direct Numerical Simulation to predict the fluids behaviour at low Bond numbers manoeuvres in microgravity conditions.

Fig. 1

 1 Fig. 1 Schematic of the study case

Fig. 3

 3 Fig. 3 Angular velocity profile of the manoeuvre

Fig. 7

 7 Fig. 7Computed force in the y-direction throughout time for manoeuvre n • 1 and for three different meshes: 64 3 , 128 3 and 256 3

Fig. 8

 8 Fig.8Computed torque in the z-direction throughout time for manoeuvre n • 1 and for three different meshes: 64 3 , 128 3 and 256 3

Fig. 12 Fig. 13

 1213 Fig. 12 Comparison of the measured and computed forces in the y-direction throughout time for manoeuvre n • 2

Fig. 14

 14 Fig. 14 Comparison of the measured and computed forces in the y-direction throughout time for manoeuvre n • 4

Fig. 15

 15 Fig. 15 Comparison of the measured and computed torques in the z-direction throughout time for manoeuvre n • 1

Fig. 16

 16 Fig. 16 Comparison of the measured and computed torques in the z-direction throughout time for manoeuvre n • 2

Fig. 17 Fig. 18

 1718 Fig. 17Comparison of the measured and computed torques in the z-direction throughout time for manoeuvre n • 3

Fig. 19 Fig. 20

 1920 Fig. 19 Comparison of the measured and computed amplitude spectra of the z-direction torque for manoeuvre n • 1

Fig. 21 Fig. 22

 2122 Fig. 21 Comparison of the measured and computed amplitude spectra of the z-direction torque for manoeuvre n • 3

Fig. 23

 23 Fig. 23 Comparison of the filtered torque and the computed torque in the z-direction throughout time for manoeuvre n • 1

Fig. 24

 24 Fig. 24 Snapshots of a bubble oscillation with τ f ill = 50%

Table 1

 1 Study parameters of the manoeuvres considered View of the inside of the FLUIDICS experiment

	N •	τ f ill	/ 1	t acc /t 1 acc	Bo c	Bo i
	1	50%	1	1	76	24
	2 3 4 5 6	75% 50% 75% 50% 50%	1 1/ √ 3 1/ √ 3 1/ √ 6 1/ √ 12	1 1 1 1 0.5	48 26 1 6 13 6	15 14 9 10 14

Table 2

 2 Comparison of the mean force in the y-direction

	N •	F y /σ D t exp	F num y	/σ D t	Error
	1a	34.84	34.16		1.97%
	1b	34.69	34.16		1.55%
	2	51.29	49.31		3.85%
	3a	11.35	11.37		0.16%
	3b	11.72	11.37		2.99%
	4	15.52	16.46		6.08%

Table 3

 3 Comparison of the computed and measured frequencies of the oscillations with the amplitude spectra

	N •	f exp ρ l D 3 t /σ	f num ρ l D 3 t /σ
	1a	2.67 ± 0.133	2.67 ± 0.133
	1b	2.67 ± 0.133	2.67 ± 0.133
	2	2 . 6 7 ± 0.133	2.93 ± 0.133
	3a	1.60 ± 0.133	1.47 ± 0.133
	3b	1.60 ± 0.133	1.47 ± 0.133
	4	1 . 4 7 ± 0.133	1.60 ± 0.133

Table 4

 4 Comparison of computed and measured frequencies of the oscillations with the videos Computed and mesured frequencies with the sensor and the video for different Bond numbers and τ f ill = 50%

	N •	f exp ρ l D 3 t /σ	f num ρ l D 3 t /σ
	1	2 . 7 6 ± 0.136	2.67 ± 0.133
	3	1 . 5 2 ± 0.065	1.47 ± 0.133
	5	1 . 0 8 ± 0.057	1.07 ± 0.133
	6	0 . 8 4 ± 0.046	0.80 ± 0.133

Acknowledgements

The authors wish to thank Airbus Defence & Space and CNES (French national space agency) for their funding and especially for the financial support of the PhD study of Alexis Dalmon. The successful collaboration between Airbus Defence & Space and CNES has provided original data by sending in the ISS the FLUIDICS experiment which achievement has represented a great challenge.

Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


domain: spi phys

text:
Julien Queyrel
email: julien.queyrel@onera.fr

Xavier Boulanger
email: xavier.boulanger@onera.fr

Laurent Castanet
email: laurent.castanet@onera.fr

James Nessel
email: james.a.nessel@nasa.gov

Michael Zemba
email: michael.j.zemba@nasa.gov

Torgeir Prytz
email: torgeir@ksat.

Preliminary Results of the THOR7 Propagation Experiment in the North Pole Region

A propagation experiment funded by ESA with support from CNES has been deployed in Svalbard (Northern Norway, 80° latitude North) in order to study propagation impairments at Ka-band caused by the troposphere at high latitudes and low elevation angles in representative conditions for GEO satcoms systems in the Artic region. In this context, a NASA Glenn Research Centre beacon receiver, operated by KSAT in Norway, measures the signal from the THOR7 satellite at a 2.65° elevation angle since March 2016. ONERA processes and analyses the data since the beginning of the experiment. At this stage 2 years of propagation data from the satellite have been gathered. An Atmospheric Numerical Simulator relying on the high resolution meteorological model (WRF-ARW) is used to generate 3D states of the atmosphere during a concurrent period to the experiment at high spatial & temporal resolutions. Propagation calculations based upon those outputs are then conducted to derive simulated attenuation time series between the beacon receiver and THOR7. Results show a fairly good statistical agreement between the measured and simulated attenuation time series whereas the simulated statistics seem to underestimate the overall attenuation for all percentages of time.

Introduction

The North Pole area is of increasing importance due to natural resources such as fishing industry and oil/gas exploration, and shipping routes are nowadays opening in those regions as a result of global warming. Together with them comes a growing need for communications at Ka-band and above, that could be fulfilled which GEO satellites or possible LEO satcom constellations (OneWeb, LEOSAT, Starlink). Those regions are also of great interest for increasing the amount of data transferred by Earth Observation data downlink systems because visibility periods of the satellite from the Earth station happen more often within a day for a polar orbit. For decades, various Earth-space propagation experiments have been carried out above 10 GHz (SIRIO [START_REF] Carassa | The Sirio programme and its propagation and communication experiment[END_REF], ETS-II [START_REF] Ishida | Propagation experiment with Japanese satellite, ETS-II (KIKU-2)[END_REF], ACTS [START_REF] Crane | ACTS propagation experiment: experiment design, calibration, and data preparation and archival[END_REF], INTELSAT [START_REF] Allnutt | INTELSAT propagation experiments: the focus and results of recent campaigns[END_REF], COMSTAR [START_REF] Cox | Results from the 19and 28-GHz COMSTAR satellite propagation experiments at Crawford Hill[END_REF], OLYMPUS [START_REF] Poires Baptista | OPEX Reference Book on Attenuation Measurement and Prediction[END_REF], ITALSAT [START_REF] Giannone | The Italsat propagation experiment[END_REF], ALPHASAT [START_REF] Koudelka | Q/V-band communications and propagation experiments using ALPHASAT[END_REF]) and allowed significant progress in the understanding and the modelling of the propagation channel at mid-latitudes. The communication link with a GSO satellite in the Artic region however suffers from specific propagation conditions that cannot be easily investigated at those latitudes: very low elevation (down to 5°) and specific impairments (large and small scale refraction, wet snow attenuation…) are not (or poorly) taken into account in conventional propagation models. In order to investigate the Ka propagation channel in those regions, ESA funded a propagation experiment with the THOR-7 GEO satellite (𝜑 = 1°W), from the NASA Glenn Research Centre beacon receiver (20.198 GHz) installed at the Svalbard Satellite teleport (SvalSat) operated by KSAT in the Sptizberg island of the Svalbard archipelago in Norway [ref]. With the help of NASA, ONERA is in charge of the data analysis of the beacon signal level and surface meteorological parameters measurements. Extra funding from the French space agency (CNES) allows ONERA to investigate further the low elevation large scale refractive effects by conducting concurrent atmospheric simulation based on the high-resolution numerical meteorological model WRF-ARW. Extensive propagation calculations allow the generated 3D state of the atmosphere to be converted into total attenuation time series, as seen by the beacon receiver at Svalsat. This paper describes the propagation experiment conducted with the THOR7 satellite at Svalsat, (section 2) and the methodology used to estimate the attenuation time series and statistics acquired for 2 years between April 2016 and March 2018 (section 3). Section 4 gives a detailed description of the meteorological numerical model used to generate the atmosphere and the computations allowing for the retrieval of the total attenuation. Finally, section 5 gathers the results of the comparison between the generated and measured attenuations series.

Ground station & data acquisition

The experiment ground station consists in a beacon receiver designed by the NASA Glenn Research Centre [START_REF] Houts | Design of a Ka-band propagation terminal for atmospheric measurements in polar regions[END_REF] that operates in the Ka-band at the 20.198 GHz beacon frequency of the THOR-7 satellite which is seen at a 2.65° elevation angle from SvalSat. A Ka-band radiometer (Radiometrics PR-2230) that can observe multiple channels between 22 and 30 GHz is also available on site. Other meteorological instruments such as an optical rain-gauge (Optical Scientific Instruments, ORG-815), a tipping bucket rain-gauge and a weather station equipped with an IR cloud sensor complement the equipment deployed for this experiment. The beacon signal (circular polarization) is reflected by a 1.2 meter Cassegrain antenna (𝐺 = 45 dBi), then both co-and cross-polar signals are amplified through a LNA (𝐺 = 44 dB, 𝑁𝐹 = 1.4 dB) and down-converted into the final intermediate frequency (𝐼𝐹 = 70 MHz) to be digitalized and processed by a SDR FFT-software implementing a modified version of a Quinn-Fernandez frequency estimation [START_REF] Zemba | Frequency estimator performance for a softwarebased beacon receiver[END_REF].

3 Data processing

Noise power inversion

Due to the inherent difficulties of the propagation campaign architecture in Svalbard (extremely low elevation angle of 2.65°) and issues with accumulated wet snow on the radiometer that induce a wrong inversion of radiometric data, a novel approach [START_REF] Nessel | Results from three years of Ka-band propagation characterization at Svalbard, Norway[END_REF] for measuring the sky brightness temperature along the line-of-sight path from the beacon receiver to the THOR 7 satellite was investigated during the experimental campaign. Using an adaptive notch filter to remove the signal from the power spectrum of the beacon receiver, the total noise power can be integrated to provide an estimate of the impact of the sky noise on the total system noise.

The following steps are preformed to calibrate the noise power and convert it into a sky brightness temperature: (i) first a 1-hour block timescale average of the 10 Hz noise power and beacon receiver time series is performed, (ii) days with stable clear sky attenuation are identified (rain gauge tipping bucket, minimal variation of signal), (iii) ECMWF ERA5 data are used over the entire month on those days to define a monthly mean clear sky attenuation and mean radiating temperature value, (iv) a clear sky level offset is calculated from the receiver power level for each month in order to calibrated coarsely the beacon-derived attenuation, (v) the RTE is used to convert beacon-derived attenuation to equivalent brightness temperature, (vi) scaling coefficients (𝑎, 𝑏) are obtained by fitting linearly the 1hour beacon-derived brightness temperature to the 1-hour noise power levels: 𝑇 𝑠𝑘𝑦 = 𝑎 × 𝑁 𝑝𝑤𝑟 + 𝑏. Results of the procedure for two months are illustrated in Figure 1. Comparisons are made between the direct line-of-sight noise power inversion approach and the scaled (10° to 2.65°) radiometric measurement and show good agreement in most cases.

Total attenuation time series

In order to derive the total attenuation time series, a visual inspection of the data is required on a daily basis. The inspection allows to spot spikes and invalid data acquisitions, but ultimately aims at identifying precipitation events (𝑡 𝑠𝑡𝑎𝑟𝑡 and 𝑡 𝑒𝑛𝑑 ) when scattering effects are not negligible anymore. Several sources of information are used for this identification: the beacon data itself, the tipping bucket and optical rain gauge measurements, the radiometric measurements associated with the noise power inversion results, and the integrated liquid water contents. Outside of those events the beacon data can be scaled to total attenuation using the absorbing attenuation derived from the noise power inversion, whereas inside those events a linear interpolation of the absorbing attenuation levels at 𝑡 𝑠𝑡𝑎𝑟𝑡 and 𝑡 𝑒𝑛𝑑 is performed.

Experimental results

From April 2016 to March 2018 valid in-excess (resp. total) attenuation measurements have been collected with an availability of 92.03% (resp. 91.98 %) of the time: 93.27% (93.17 %) the first year and 90.79% (90.79 %) the second year. Surface parameters, including tipping bucket rain gauge measurements are available for 91.71 % of the time (91.67 % the first year, 91.75 % the second year), and the data collected from an optical rain gauge are available for 71.76 % of the time (79.63 % the first year, 63.88 the second year). In Figure 2 examples of attenuation (in-excess and total) are shown for two distinct months. Various CCDF of the rainfall rate are plotted in the left panel of Figure 3: the optical rain gauge CCDF is computed using the total rainfall (dotted curve), using only the data when 𝑇 𝑠𝑢𝑟𝑓 > 0°C (cyan curve), or using only the data except when the condition code of the instrument is 'snow' (cyan dotted curve). The differences between the curves show the difficulty in estimating the actual rainfall rates: the presence of snow and the impact of freezing temperature on the instruments operations are supposed to be responsible for those behaviours.

Statistical analyses have also been conducted for the in-excess and total attenuation measured time series. After removing the scintillation using a raised cosine filter with a 0.01 Hz cut-off frequency [START_REF] Rytir | Clear-air scintillation and multipath for low-elevation high-latitude satellite communication links[END_REF], the CCDF plotted in the right panel of Figure 3 are derived. Monthly statistics of the attenuation have also been computed: they show an increase in attenuation exceeded for any percentages during the summer period. 4 Atmospheric Numerical Simulator (ANS)

Generation of the atmosphere

In order to compute the propagation conditions along the THOR7 link towards the beacon receiver, an atmospheric numerical simulator has been used along with an electromagnetic module [START_REF] García | Use of WRF model to characterize propagation effects in the troposphere[END_REF], [START_REF] Jeannin | Atmospheric channel simulator for the simulation of propagation impairments for Ka band data downlink[END_REF]. The high-resolution numerical weather prediction model WRF-ARW (v3.4) [START_REF] Skamarock | A Description of the Advanced Research WRF Version 3[END_REF] is used to generate a 4D description of the atmosphere around the beacon receiver at a high temporal (5 and spatial resolution (2 km). The datacubes thus produced are approximately 160 km x 160 km horizontal squares and 36 vertical model levels span ~ 20 km of atmosphere. The meteorological model is initialized with the ECWMF ERA-Interim re-analysis database which provides surface and pressure levels parameters every 6 hours with a 0.75° resolution. These same re-analysis data are also used by WRF-ARW as boundary conditions every 6 hours. In order to mimic realistic passed states of the atmosphere, simulations are limited to one day with a previous 12 hours spin-up period (the results of which are deleted latter-on) giving the opportunity to the model to provide stable outputs. In order to reach the 2km resolution, 3 nested domains are used: a coarser domain with a 30 km resolution on a polar stereographic projection and static fields initialized directly by way of an interpolation from the ERA-Interim database, and two successive nests at 6 km and 2 km resolution (see Figure 5).

Attenuation retrieval

The atmospheric parameters of interest (𝑃, 𝑇, 𝜌 𝑣𝑎𝑝 , 𝜌 𝑐𝑙𝑜𝑢𝑑 , 𝜌 𝑟𝑎𝑖𝑛 ) are then converted by an electromagnetic module written in Python into specific attenuations (dB/km) of the various constituents of the troposphere known to produce attenuation effects at Ka-band frequencies: water vapour, oxygen, clouds and rain (see integrated columnar contents in Figure 6). At this point in the study, attenuation due to snow is not modelled by the electromagnetic module. The physical models used by the module are: the line-by-line spectral model of Recommendation ITU-R P.676-11 (annex 1, §1) for oxygen (𝛾 𝑜𝑥 ) and water vapour (𝛾 𝑣𝑎𝑝 ) and the model based on Rayleigh scattering of Rec. ITU-R P.840-7 (annex 1, §2) for clouds droplets (𝛾 𝑐𝑙𝑜𝑢𝑑 = 𝐾 𝑙 (𝑓, 𝑇) × 𝜌 𝑐𝑙𝑜𝑢𝑑 ). Rain specific attenuation (𝛾 𝑟𝑎𝑖𝑛 ) is derived from Mie's theory (extinction cross-section 𝜎 𝑒𝑥𝑡 (𝐷)) assuming spherical droplets and a Marshall-Palmer raindrop size distribution (𝑁(𝐷) = 𝑁 0 𝑒 -Λ𝐷 ) which slope parameter Λ is fixed by the precipitating water content of each pixel cube (𝜌 𝑟𝑎𝑖𝑛 ) [START_REF] Hulst | Light scattering by small particles[END_REF]:

𝛾 𝑟𝑎𝑖𝑛 = 4.343 × 10³ ∫ 𝜎 𝑒𝑥𝑡 (𝐷)𝑁(𝐷)𝑑𝐷 𝐷 𝑚𝑎𝑥 0
The Mie's coefficients used to infer 𝜎 𝑒𝑥𝑡 (𝐷) being computationally expensive to derive, a look-up table linking 𝑇 and 𝜌 𝑟𝑎𝑖𝑛 to 𝛾 𝑟𝑎𝑖𝑛 is used. In order to compute the total attenuation suffered by the signal when reaching the beacon receiver at low elevation angle (𝜃 ∼ 2.65°) large scale geometric refractive effects are taken into account. First, for each time step a vertical slice of atmosphere, directed toward the satellite (azimuth ψ), is extracted from the WRF datacube (see Figure 7). Then, a dichotomy algorithm is used to estimate the actual elevation angle and geometric path between the beacon receiver and the satellite by propagating a geometric ray through the various layers of atmosphere within the slice. To that aim, the refractive index of all points within the slice are computed from 𝑃, 𝑇 and 𝜌 𝑣𝑎𝑝 (cf. Rec. ITU-R P.453-13) and Snell/Descartes' geometric laws are used as suggested in the spherically layered atmosphere model of ITU-R P.676-11 (annex 1, §2.2) to propagate the signal for one layer to the next, starting on Svalsat. The process is repeated for various elevation angles until the elevation angle allowing to join the beacon receiver and the satellite is eventually found within an acceptable error margin. The specific attenuation is then integrated (iii) along the newly found geometric path to recover the total attenuation (dB) at the current time step. This way, 2 years of simulated attenuation time series (Δ𝑡 = 5 min) concurrent with the beacon measurements are generated.

Figure 7 -Vertical slice of atmosphere in the direction of THOR7 @ azimuth ψ (illustration)

Results and conclusion

Comparisons between simulated and measured parameters have been conducted both for meteorological surface parameters and tropospheric attenuation.

For year 1 and year 2 of the experiment, measured and simulated surface pressure and temperature at Svalsat and are in very good agreement (see time series example in Figure 8). The total attenuation time series (measured, scintillation not removed, and simulated) show a fairly good agreement (see example for 2 months in Figure 9) Finally, the CCDFs of simulated total attenuation have been computed for both years and the comparison with the measured CCDFs is shown in Figure 10. The results for the two years in Figure 10 show a clear underestimation of the attenuation exceeded by the Atmospheric Numerical Simulator, for nearly all percentages of time (≤ 60-70%). The effect of snow -not considered within the simulations generated by the ANSis probably the main reason of the under prediction of the simulated CCDF with respect to the measured one. Looking at the CCDFs of solid and liquid precipitations in the left panel of Figure 3, it is clear that snow should be a significant contributor to precipitation rates, and wet snow should be therefore a strong contributor to total attenuation at these frequencies [START_REF] Linlor | Permittivity and attenuation of wet snow between 4 and 12 GHz[END_REF]. Moreover, the removal of scintillation by filtering the measured attenuation with a fixed cutt-off frequency could also be responsible for an overestimation of the CCDFs' exceeded attenuation, as it may fail to remove all scintillation from the data.

In order to improve the ANS simulations, various efforts are being made: (i) a snow attenuation physical model is currently developed to be integrated into the electromagnetic module, and (ii) the use of a more resolved database from the ECMWF (eg. ERA5: 0.25°x0.25°, every hour) as an input to WRF-ARW is also being investigated.

Figure 1 -

 1 Figure 1 -Correlation between the 𝑇 𝑠𝑘𝑦 measured from the beacon attenuation and the noise power integrated on the 1 MHz beacon reciever bandwidth for October & November 2016.

Figure 2 -

 2 Figure 2 -In-excess and total attenuation time series measured by the THOR7 beacon receiver after data processing (August 2016 & January 2017)

Figure 3 -

 3 Figure 3 -Experimental CCDFs for rainfall rates (left) and in-excess and total attenuation (right) during the 2 years (April 2016 -March 2018) of the experiment

Figure 4 -

 4 Figure 4 -CCDF of scintillation (fade, enhancement and standard deviation) (left) and monthly standard deviation (right)

Figure 5 -

 5 Figure 5 -The 3 nested geographical domains of the WRF-ARW simulation (left), and position of the high-resolution domains d03 with direction of the THOR7 satellite (right)

Figure 6 -

 6 Figure 6 -Integrated water vapour (left), cloud (middle) and rain (right) content on the last domain centered around Svalsat at a specific time step of the simulation (2017/06/06T08:30)

Figure 8 -

 8 Figure 8 -Surface pressure (up) and temperature (bottom) measured by the meteorological station at Svalsat (black) and generated by the atmospheric model (red).

Figure 9 -

 9 Figure 9 -Total attenuation time series measured by the beacon receiver (grey curve, scintillation not removed), derived from ERA-Interim data every 6h (green cross), derived from ERA5 data every hour (blue curve) and computed from the simulated atmosphere generated by WRF-ARW (red curve)

Figure 10 -

 10 Figure 10 -Total attenuation CCDF measured (black) and simulated (red) for both years


Julie Sliwak

Miguel F Anjos
email: anjos@stanfordalumni.org

Lucas Létocart
email: letocart@lipn.univ-paris13.fr

Jean Maeght

Emiliano Traversi
email: traversi@lipn.univ-paris13.fr

Improving Clique Decompositions of Semidefinite Relaxations for Optimal Power Flow Problems

Keywords: Clique decomposition, Optimal Power Flow, Semidefinite Programming

Semidefinite Programming (SDP) provides tight lower bounds for Optimal Power Flow problems. However, solving large-scale SDP problems requires exploiting sparsity. In this paper, we experiment several clique decomposition algorithms that lead to different reformulations and we show that the resolution is highly sensitive to the clique decomposition procedure. Our main contribution is to demonstrate that minimizing the number of additional edges in the chordal extension is not always appropriate to get a good clique decomposition.

Introduction

The Alternating Current Optimal Power Flow (ACOPF) problem is a nonconvex optimization problem whose purpose is to provide an operating point of the power network minimizing generation costs. Solving this problem is essential for RTE, the French transmission system operator, because ACOPF problems are subproblems of decision problems involving millions of euros, e.g. operational planning and grid development. Yet there is still no efficient method to solve ACOPF problems to global optimality.

Many convex relaxations have been proposed to evaluate the quality of feasible solutions computed with local methods [START_REF] Molzahn | A Survey of Relaxations and Approximations of the Power Flow Equations[END_REF]. SDP relaxations often provide tight lower bounds [START_REF] Lavaei | Zero duality gap in optimal power flow problem[END_REF] that are useful to prove global optimality. Some promising global optimization methods are based on these relaxations : Godard et al. propose an adaptation of the Mixed-Integer Quadratic Convex Reformulation method to ACOPF problems in [START_REF] Godard | Novel approach towards global optimality of optimal power flow using quadratic convex optimization[END_REF], Gopalakrishnan et al. present a branch-and-bound approach using SDP relaxations in [START_REF] Gopalakrishnan | Global optimization of optimal power flow using a branch & bound algorithm[END_REF] and Josz et al. apply the Lasserre hierarchy in [START_REF] Josz | Application of the moment-sos approach to global optimization of the opf problem[END_REF] to achieve global optimality. All these methods depend on large-scale SDP problems being solved efficiently.

Many algorithms are proposed in the literature to solve SDP problems. In this paper, we focus on interior-point methods that are the most reliable and accurate algorithms as far as we know. Interior-point methods solve efficiently small-to-medium-sized SDP problems but do not scale well for large-scale problems because of the Hessian equation that implies forming and factorizing a fully-dense matrix at each step. Yet, sparse large-scale problems can be tackled exploiting sparsity. In this paper, we focus on clique decomposition techniques [START_REF] Fukuda | Exploiting sparsity in semidefinite programming via matrix completion i : General framework[END_REF][START_REF] Nakata | Exploiting sparsity in semidefinite programming via matrix completion ii : Implementation and numerical results[END_REF][START_REF] Vandenberghe | Chordal graphs and semidefinite optimization[END_REF] that perform well for ACOPF problems. Jabr [START_REF] Rabih | Exploiting sparsity in sdp relaxations of the opf problem[END_REF] and Molzahn et al. [START_REF] Daniel K Molzahn | Implementation of a large-scale optimal power flow solver based on semidefinite programming[END_REF] have each proposed a way to use clique decompositions for ACOPF problems and both propose to work on problems formulated in complex variables. However, there are many other possible clique decomposition procedures that provide different reformulations of a SDP problem. Most of the clique decomposition algorithms seek to minimize the number of added edges in the chordal extension. It results in decompositions with many small cliques for ACOPF problems and many linking constraints are required to handle overlaps between cliques. As linking constraints can slow down the resolution, Nakata et al. [START_REF] Nakata | Exploiting sparsity in semidefinite programming via matrix completion ii : Implementation and numerical results[END_REF] and Molzahn et al. [START_REF] Daniel K Molzahn | Implementation of a large-scale optimal power flow solver based on semidefinite programming[END_REF] propose heuristics for merging cliques in order to reduce their number. Yet, merging cliques means adding more edges in the chordal extension. Therefore, classical clique decomposition approaches may not always be the most appropriate since they focus on the number of added edges, regardless of important criteria such as the size of the largest clique or the number of linking constraints.

In this paper, we show that different clique decompositions are not equivalent in terms of resolution by comparing different chordal extension algorithms on RTE [START_REF] Josz | Ac power flow data in matpower and qcqp format : itesla, rte snapshots, and pegase[END_REF] and other MAT-POWER [START_REF] Daniel Zimmerman | Matpower : Steady-state operations, planning, and analysis tools for power systems research and education[END_REF] datasets. Our main contribution is to demonstrate in two ways that computing a chordal extension with a minimum number of added edges is not a relevant choice for ACOPF problems. The first test consists in a comparison between reformulations coming from the SDP problem formulated in complex variables and reformulations coming from the same problem but formulated in real variables. The conclusions of this comparison are that less edges are necessary in the real case but the reformulations are more performant in the complex case. The second test relies on a new clique combination algorithm that improves a given decomposition by adding more edges in the chordal extension. The edges are added in such a way as to minimize the number of linking constraints while keeping small cliques. Computational tests show that this algorithm significantly speeds up the resolution of the largest OPF instances. Both tests highlight the fact that the number of added edges is not always the best criterion to minimize, which proves that the computation of the chordal extension is a question that deserves to be investigated. This paper is organized as follows. Section 2 describes briefly the OPF problem and the classical rank relaxation. Common clique decomposition techniques are presented in section 3 along with computational comparisons. Section 4 compares reformulations coming from the complex SDP problem and reformulations coming from the real SDP problem. Section 5 details our clique combination algorithm and resulting enhancements. Section 6 concludes the paper.

ACOPF problem and rank relaxation

In this section, we present briefly the ACOPF formulation and its classical semidefinite rank relaxation.

ACOPF formulation

The ACOPF problem is defined on a power transmission network that can be modelled as an oriented graph T (N, B) where N represents the electrical buses and B the branches (transmission lines, transformers). Let us denote G ⊂ N the subset of generator buses. The ACOPF problem is defined as follows :

min g∈G c g Re(S gen g ) + k g s.t. S gen n = S load n + b∈B -(n) S dest b (v) + b∈B + (n) S orig b (v) ∀n ∈ N (v min n ) 2 ≤ |v n | 2 ≤ (v max n ) 2 ∀n ∈ N P min g ≤ Re(S gen g ) ≤ P max g ∀g ∈ G Q min g ≤ Im(S gen g ) ≤ Q max g ∀g ∈ G |i orig b (v)| 2 ≤ (i max b ) 2 ∀b ∈ B |i dest b (v)| 2 ≤ (i max b ) 2 ∀b ∈ B S gen n = 0 ∀n ∈ N \G v n ∈ C, S gen n ∈ C ∀n ∈ N (ACOPF)
where Re(z) stands for the real part of the complex number z and Im(z) for the imaginary part. All constant parameters are in bold.

S orig b (v) = v o(b) i orig b (v) S dest b (v) = v d(b) i dest b (v) (1)
Note that we use thermal limits modelled with current and linear generation costs to have a formulation closer to what is used in practice at RTE.

Rank relaxation

Since costs are linear, the ACOPF problem can be written in terms of voltage variables only. The resulting problem is a Quadratically Constrained Quadratic Program (QCQP) with complex variables that can be expressed in a compact way :

   min v H Q 0 v s.t. v H Q p v ≤ a p ∀p = 1..m v ∈ C n (2)
The classical rank relaxation is constructed by introducing the Hermitian matrix W = vv H . This constraint is equivalent to W 0 and rank(W ) = 1. The rank relaxation is the SDP problem obtained by suppressing the rank constraint :

   min Q 0 • W s.t. Q p • W ≤ a p ∀p = 1..m W 0 (3)
This SDP problem with complex numbers can be converted into a SDP problem with real numbers using the rectangular representation and a symmetric matrix of size 2n × 2n.

Clique decomposition analysis

In this section, we review the clique decomposition techniques based on matrix completion introduced by Fukuda et al. in [START_REF] Fukuda | Exploiting sparsity in semidefinite programming via matrix completion i : General framework[END_REF][START_REF] Nakata | Exploiting sparsity in semidefinite programming via matrix completion ii : Implementation and numerical results[END_REF]. Then we show that the resolution with interior-point methods is highly sensitive to the clique decomposition.

Clique decomposition framework

The clique decomposition relies on the matrix completion theorem [START_REF] Vandenberghe | Chordal graphs and semidefinite optimization[END_REF]. This allows to replace the Positive Semidefinite (PSD) constraint on a big matrix by several PSD constraints on smaller submatrices at the price of adding linking constraints between these submatrices.

Let us define A as the aggregate sparsity pattern of a SDP problem, i.e., the matrix of nonzeros entries in the data matrices (objective and constraint matrices). Let us denote G A as the graph associated to A. The matrix completion theorem can be applied if and only if G A is chordal. A chordal graph is a graph with no induced cycle of four vertices or more. A clique is a subset of vertices that are all connected together. A maximal clique is a clique that is not included in another clique.

The general framework for clique decomposition contains three steps. The first step consists in computing a chordal extension H of the aggregate sparsity pattern G A . The second step consists in determining the list of maximal cliques in H, denoted by L = {C 1 , ..., C r }. These maximal cliques define the submatrices that must be positive semidefinite in the reformulated SDP problem. From the list of maximal cliques L, a clique graph W can be defined with nodes corresponding to maximal cliques and weighted edges between each pair of vertices defined by the number of shared nodes between the two corresponding cliques. The third step consists in computing a clique tree U by computing a maximum-weight spanning tree for W. This clique tree allows to specify linking constraints between submatrices.

It is clear that the clique decomposition depends on the chordal extension H but it is not completely clear how to precisely define a "good" decomposition from a practical point of view. One reasonable choice is to use decompositions that minimize the number of additional edges. However, finding the minimal chordal extension is NP-complete [START_REF] Yannakakis | Computing the minimum fill-in is np-complete[END_REF]. Bergman et al. have recently proposed an exact model with exponentially many constraints to find a minimal extension [START_REF] Bergman | On the minimum chordal completion polytope[END_REF] but it is only usable for small instances. For this reason, several efficient heuristics for adding few edges are available in the literature [START_REF] Heggernes | Minimal triangulations of graphs : A survey[END_REF]. The most commonly used is based on a Cholesky factorization of G A whose rows and columns have been permuted according to a minimum or approximate minimum degree ordering [START_REF] Patrick R Amestoy | An approximate minimum degree ordering algorithm[END_REF]. The ordering has a significant impact on the fill-in, that is, on the number of edges added to G A . Once H is computed, the list of maximal cliques can be computed in linear time [START_REF] Robert E Tarjan | Simple linear-time algorithms to test chordality of graphs, test acyclicity of hypergraphs, and selectively reduce acyclic hypergraphs[END_REF]. Finally, the most widely used exact algorithm to compute the clique tree is Prim's algorithm [START_REF] Clay | Shortest connection networks and some generalizations[END_REF] but there are other algorithms and they can lead to different optimal clique trees. However, the decompositions obtained usually do not differ much : the number of linking constraints do not differ significantly from one tree to another because the trees have the same overall weight.

We show in the following subsection that different chordal extension algorithms give different decompositions and that these different decompositions are not equivalent as regards resolution time.

Numerical comparison of two chordal extension heuristics for SDP relaxations of ACOPF formulated in complex variables

There are several algorithms to compute chordal extensions and the clique decompositions can vary greatly from one algorithm to another, which has an impact on the resolution. In this subsection, we compare clique decompositions coming from two different algorithms : a Cholesky factorization with an AMD ordering and the Minimum Degree (MD) heuristic [START_REF] Berry | The minimum degree heuristic and the minimal triangulation process[END_REF] that is based on a dynamical computation of a minimum degree ordering of the vertices. We have picked these two algorithms because they are among the best heuristics known. Both algorithms are applied on the SDP relaxations of ACOPF formulated in complex variables, therefore providing maximal cliques with complex variables. In order to solve the SDP problems, these cliques are finally converted into cliques with real variables by doubling them.

Table 1 presents some results of this comparison on MATPOWER instances with more than 1000 buses. The tests were carried out on a Processor Intel R Core TM i7-6820HQ CPU @2.70GHz using the module MathProgComplex.jl [START_REF] Sliwak | A julia module for polynomial optimization with complex variables applied to optimal power flow[END_REF] with Julia 1.0.3. and the SDP solver MOSEK 8.1.0.72 was used for the resolution. We also used the packages JuMP.jl [START_REF] Dunning | Jump : A modeling language for mathematical optimization[END_REF], Light-Graphs.jl [START_REF] Fairbanks | Juliagraphs/lightgraphs.jl : an optimized graphs package for the julia programming language[END_REF] and MetaGraphs.jl.

Let nc M D , nlc M D and t M D be respectively the number of maximal cliques, the number of linking constraints and MOSEK resolution time obtained with the MD heuristic. Similarly, let nc AM D , nlc AM D and t AM D be respectively the number of maximal cliques, the number of linking constraints and MOSEK resolution time obtained with Cholesky and AMD. Table 3.2 presents some ratios between these three quantities. All results are presented in percentages. This table is not exhaustive but some conclusions can be drawn from it. First, clique decompositions are different from one algorithm to another. In particular, there are large differences in resolution time even if the decompositions have about the same number of cliques. Theses differences can be explained by the differences in the number of linking constraints, especially for the biggest instances. Another factor can be the size of the biggest clique in the decomposition. It is also possible to generate decompositions which differ both in the number of cliques and in the number of linking constraints so differences in resolution time can be much more impressive. In particular, in preliminary tests we experienced that there exists "bad" clique decompositions by computing several "bad" orderings for Cholesky (e.g. random ordering, maximum degree ordering). In this case, a "bad" clique decomposition means a decomposition with many linking constraints, which can lead to memory issues with interior-point methods. This confirms the idea that the computation of the chordal extension is a point worth studying.

Instance nc M D nc AM D -1 nlc M D nlc AM D -1 t M D t AM
The results shown in this section are a good motivation to continue exploring the impact of the clique reformulation on the resolution. The next section focus on the impact of computing the clique decomposition on SDP relaxations of ACOPF problems formulated in complex variables or on SDP relaxations of ACOPF problems formulated in real variables.

Comparison of clique decompositions computed in the complex and in the real framework

SDP relaxations of ACOPF problems are naturally expressed with complex numbers but in practice, in order to be solved by the available solvers, they are rewritted in terms of only real variables. The process of rewriting the model with real variables leads to a different (and bigger) aggregate sparsity pattern A. Such a A could be itself decomposed in a way that would not be possible with the original formulation with the complex variables. Therefore, there are two possibilities for the clique decomposition : either applying the procedure on the complex SDP formulation and converting the complex cliques to real cliques or directly applying the clique decomposition procedure on the real SDP formulation. However, computing the chordal extension on the complex or the real problem is not equivalent. We first show it theoretically on a small MATPOWER [START_REF] Daniel Zimmerman | Matpower : Steady-state operations, planning, and analysis tools for power systems research and education[END_REF] example, LMBM3.

Let us define G c (respectively G r ) as the graph associated to the aggregate sparsity pattern A of the complex (respectively real) SDP relaxation. For any ACOPF instance as defined in section 2, G c is the network graph T (N, B). Let us denote H c (respectively H r ) the chordal extension computed from G c (respectively G r ) with a given algorithm. The chordal extension H c can be converted to real numbers in the same way as the complex graph G c can be converted to the real graph G r . Let us denote H c real this conversion of H c . H c real is also a chordal extension of G r . The objective of this section is to compare H c real and H r . For LMBM3, G c is a clique of size 3, denoted by K 3 and G r is the graph with the 6 following nodes V r = {1 Re , 1 Im , 2 Re , 2 Im , 3 Re , 3 Im } and the 12 following edges :

E r = {(1 Re , 2 Re ), (1 Re , 2 Im ), (1 Im , 2 Re ), (1 Im , 2 Im ), (1 Re , 3 Re ), (1 Re , 3 Im ), (1 Im , 3 Re ), (1 Im , 3 Im ), (2 Re , 3 Re ), (2 Re , 3 Im ), (2 Im , 3 Re ), (2 Im , 3 Im )}.
Since G c is complete and thus chordal, H c = G c = K 3 regardless of the chordal extension algorithm used. Clearly, there is one maximal clique with 3 nodes and no linking constraints. When converted to real numbers, we get H c real = K 6 , the clique of size 6, which accounts to adding 3 edges to G r . This chordal extension gives one maximal clique with 6 real nodes and no linking constraints. On the other hand, G r is not chordal because there are several induced subgraphs of size 4, e.g. 1 Re -2 Im -1 Im -2 Re . It suffices to add the two edges (2 Re , 2 Im ) and (3 Re , 3 Im ) to G r to get a minimal chordal extension H r1 . H r is chordal because 1 Re -1 Im -2 Re -2 Im -3 Re -3 Im is a perfect elimination ordering. This chordal extension gives two maximal cliques : {2 Re , 2 Im , 3 Re , 3 Im , 1 Re } and {2 Re , 2 Im , 3 Re , 3 Im , 1 Im } and 10 linking constraints.

This small example proves that it is not equivalent to compute the clique decomposition on the complex problem or on the real problem from the theoretical point of view : more edges are added in the complex case to get a chordal extension. This idea is confirmed by computational tests that show that there are on average 50% more edges added when the chordal extension is computed in the complex case rather than in the real case for MATPOWER instances with more than 1000 buses. However, adding less edges in the chordal extension is not necessarily advantageous from the numerical resolution point of view. Indeed, adding less edges in the chordal extension can mean having more cliques and more linking constraints as shown in Table 2. This table presents comparisons between clique decomposition coming from the complex case (computed on G c ) and clique decomposition coming from the real case (computed on G r ). More precisely, both clique decompositions are based on a chordal extension obtained by a Cholesly factorization with an Approximate Minimum Degree (AMD) [START_REF] Patrick R Amestoy | An approximate minimum degree ordering algorithm[END_REF] ordering. Let L c (respectively L r ) be the list of maximal cliques obtained from the chordal extension H c (respectively H r ). The second column represents the ratio for the number of cliques defined as |L r | |L c | -1. Similarly, the third column represents the ratio for the number of linking constraints, i.e., nlc r nlc c -1 with nlc c (respectively nlc r ) the number of linking constraints when the clique decomposition is done for the complex (respectively real) problem. Finally, the fourth column represents the ratio for MOSEK resolution time, i.e., t r t c -1 with t c (respectively t r ) the resolution time with the clique decomposition coming from the complex (respectively real) problem. All results are presented in percentages and the tests were carried out in the same manner as in subsection 3.2.

Instance |L r | |L c | -1 nlc r nlc c -1 t r t c -1 case1888rte 51% 24% 17% case3012wp 39% 15% 4% case6468rte 40% 19% 58% case9241pegase 48% 18% 28% case13659pegase 59% 20% 23% TAB. 2 -Comparison of clique decompositions on G c or G r
This table shows that the clique decomposition procedure on G r gives approximatly 50% more cliques and 20% more linking constraints than on G c , which slows down the resolution (the computing time increases by at least 4% more and up to 58% more). Therefore, even if less edges are added to the graph G r when working directly in real numbers, it seems better to work in complex numbers. This comparison demonstrates that computing a chordal extension with a minimum number of added edges is not necessarily profitable. The next section presents a clique combination algorithm that leads to the same conclusion.

A new clique combination approach

Several heuristics do already exist to compute chordal extensions but to the best of our knowledge, they all try to minimize the number of added edges. However, nothing proves that a chordal extension with few added edges gives a decomposition easier to solve. In contrast, minimizing the size of the biggest clique or the number of linking constraints could be more interesting because these quantities have an impact on the resolution. To support this claim, we propose a new clique combination algorithm that allows to reduce the number of linking constraints in a given decomposition while keeping small matrices. This algorithm consists in merging some cliques that are adjacent in the clique tree by solving the problem (IP) below. There is one binay variable per edge in the clique tree. If this variable equals 1, it means that the cliques at the ends of the edge will be merged. We allow at most one combination per clique and combinations that result in a clique of size greater than a certain size are forbidden. The objective is to maximize the number of linking constraints that will be eliminated thanks to the combinations. More precisely, the model is the following : This algorithm has been tested with the solver Xpress for S max = 50 and the results are presented in Table 3. This table represents the evolution of the total resolution time depending on the number of times the combination algorithme (IP) is used. This table shows that applying this algorithm at least once improves significatively MOSEK resolution time and the best time is obtained when the algorithm is applied 4 times. Therefore, it is worth combining cliques, i.e., adding more edges in the chordal extension, to improve the performance of a given decomposition for ACOPF problems.

Conclusion and future work

Sections 4 and 5 show that, for ACOPF problems, there is no real reason to build chordal extensions with as few added edges as possible since adding more edges can lead to better performance. These results encourage us to change the approach for computing chordal extensions and in the future we will focus on the formulation of an adapted heuristic that computes performant chordal extensions for ACOPF problems.

  c g and k g represent respectively the linear cost and the constant cost of the generator bus g ∈ G. S load n represents the load. v min n and v max n are bounds on the voltage magnitude at bus n. on the active (respectively reactive) power at generator bus g ∈ G. i max b represents the current limit for branch b ∈ B. The variables are the voltage v n and the power S gen n at each bus n ∈ N with S gen n = 0 for all buses n ∈ N \G. Both are complex variables. B -(n) is the set of entering branches at bus n and B + (n) the set of exiting branches. o(b) stands for the origin of branch b, d(b) for the destination of branch b. The currents i orig b (v) and i dest b (v) of a branch b ∈ B are linear functions of v o(d) and v d(b) depending on physical characteristics. The power on the lines are defined as follows :

max 0. 5

 5 e∈E nc(C e )(nc(C e ) + 1)x e s.t. size(C e )x e ≤ S max ∀e ∈ E e∈E(i) x e ≤ 1 ∀i ∈ L x e ∈ {0, 1} ∀e ∈ E (IP) with -U = (L, E) the clique tree computed for a given decomposition ; nc(C e ) = |C i C j | the number of common nodes for e = (C i , C j ) ∈ E ; size(C e ) = |C i | + |C j | -|C i C j | the size of the merged clique for e = (C i , C j ) ∈ E ; -S max a parameter that defines the maximal size of merged cliques ; -E(i) the set of edges that have i as extremity. Number of times (IP) time (seconds) 3397 2316 1807 1691 1533 1710 1826 TAB. 3 -Total MOSEK resolution time for MATPOWER instances

There are other possibilities for H r , e.g., if the edges (1Re, 21m) and (3Re, 3Im) are added, but the results will be the same.