A Dalmon

M Lepilliez

A Tanguy

Pedrono

B Busset

H Bavestrello

J Mignot

Direct numerical simulation of a bubble motion in a spherical tank under external forces and microgravity conditions

Keywords: bubble dynamics, computational methods, drops and bubbles

We present, in this paper, numerical simulations of bubble sloshing in a spherical tank, resulting from a tank rotation around a fixed axis in microgravity conditions. This configuration is of great interest in space applications where sloshing can have harmful effects on the stability of satellites. Depending on the dimensionless numbers characterising this phenomenon, our study is focused on the motion and the deformation of a bubble, initially at rest, which is set in motion when the manoeuvre is starting until it reaches a constant rotation speed around the axis. It is shown in this article that, during the first stage of the manoeuvre, the motion of the bubble is essentially driven by the inertial force that depends on the angular acceleration. Next, when the angular velocity is increasing, the centrifugal force being dominant, the trajectory of the bubble is pushed towards the direction between the centre of the tank and the axis of rotation. Finally, when the angular velocity becomes constant, the bubble, reaching a quasi-steady position, is deformed and pressed against the solid boundary of the tank. A quantified description of these phenomena is proposed through a parametric study varying the essential dimensionless numbers, i.e. the Bond number based on the angular velocity, and another Bond number based on the angular acceleration. As the temporal evolution of the forces acting on the satellite wall is of utmost importance for designing satellites and manoeuvres, we also present an analysis characterising the latter. We also detail the first comparisons between the numerical simulations and the Fluidics experiment performed in the International Space Station (ISS) in microgravity conditions. Thanks to these comparisons, we can validate the simulations in configurations of interest.

Introduction

The sloshing of a liquid phase contained in a tank under the effect of an external force is a topic of high interest in many applications. Several disruptions can affect the stability of a satellite. They may originate from the space environment, as solar winds, magnetic fields or atmospheric residues. Furthermore, they may come from the satellite itself such as inertial reactions due to propulsion or internal mechanisms, construction defects, sloshing generated by flexible components (solar panel or antenna) or fluids inside the tank. It is well known [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] that the magnitude of fluid sloshing disruption is one of the most important of all those cited. Indeed, during its lifetime in space, a satellite will execute a wide range of operational manoeuvres, such as reorientation of an Earth observation device or solar panels. During these manoeuvres, inertial forces acting on propellant tanks will induce sloshing of a helium bubble formed inside these tanks after the latter gas was injected to maintain a sufficient pressure. By generating forces and torques on the structure, this bubble sloshing can deteriorate the quality of satellite imaging. Substantial efforts are undertaken to predict this unexpected phenomenon in order to propose technical solutions to decrease the amplitude of bubble sloshing in microgravity conditions. Several previous studies have been conducted to predict slosh behaviour inside a satellite tank, as in [START_REF] Abramson | The dynamic behaviour of liquids in moving containers, with applications to space vehicle technology[END_REF] and [START_REF] Dodge | The new dynamic behaviour of liquids in moving containers[END_REF], where analytical models were proposed, respectively based on a one-dimensional mass-spring system or a pendulum system. The motion of the centre of mass was deduced from a second-order linear equation which gives the sloshing frequency, the damping ratio and allows calculation of the sloshing forces. These models provide a correct order of magnitude of the sloshing parameters for high values of Bond number (Bo > 10) and are currently used in the space industry, in particular for high thrust manoeuvres. However, they become inefficient for applications involving lower values of the Bond number where surface tension effects cannot be neglected, leading to nonlinear equations. Several theoretical studies have also been proposed recently [START_REF] Faltinsen | A multimodal method for liquid sloshing in a twodimensional circular tank[END_REF][START_REF] Faltinsen | Steady-state liquid sloshing in a rectangular tank with a slat-type screen in the middle: quasilinear modal analysis and experiments[END_REF][START_REF] Ikeda | Nonlinear liquid sloshing in a square tank subjected to obliquely horizontal excitation[END_REF][START_REF] Faltinsen | Multimodal analysis of weakly nonlinear sloshing in a spherical tank[END_REF] in configurations involving circular, spherical or rectangular tanks. Considering the complexity of these phenomena where a strongly perturbed interface can interact with a complex geometry in nonlinear regimes, numerical simulations, as the ones presented in [START_REF] Veldman | The numerical simulation of liquid sloshing on board spacecraft[END_REF], [START_REF] Liu | A numerical study of three-dimensional liquid sloshing in tanks[END_REF], are complementary to theoretical studies in order to push back the limits of our knowledge in the field of sloshing. However, none of the previous studies have considered the bubble motion in a spherical tank subjected to an external force, as can occur during the manoeuvring of a satellite in microgravity conditions. Little experimental data exist on bubble sloshing in microgravity conditions, and significant advances are expected by using the direct numerical simulation to investigate the low Bond number regimes which are of great interest for satellite station keeping, rather than high thrust manoeuvres. Consequently, the present study is focused on configurations involving external forces which have the same order of magnitude or a smaller magnitude than the capillary forces, i.e. for a low Bond number varying between 0.05 and 5.

The objectives of this paper are as follows; (i) first, by providing accurate comparisons between well-resolved three-dimensional numerical simulations and microgravity conditions experiments, a strong validation of the overall numerical methodology is given in the considered configurations, (ii) secondly, depending on the dimensionless numbers, a detailed numerical study of the bubble motion is presented both in the first unsteady phase and in the steady phase of the manoeuvring, finally (iii) a numerical investigation on the steady bubble deformation depending on the centrifugal Bond number is detailed. As this parameter is of utmost importance for the design of satellites, our study also includes an analysis on the forces acting on the tank wall during the manoeuvre.

In § 2, the tank configuration in microgravity conditions is presented. The dimensionless analysis of the governing equations highlights the dominating effects coming from the inertial forces which are expressed with two Bond numbers. Some of the physical hypotheses we make are also justified in this section, such as the incompressibility of both phases and the perfectly wetting contact angle model. Section 3 details the formalism and the numerical methods we implement to model the fluids' behaviour inside the tank. Then the validation of our numerical code with experimental data collected in the International Space Station is presented in § 4. The evolution of the centre of mass of the fluids inside the tank is investigated in § 5. After presenting a space convergence study and an example of bubble behaviour during a manoeuvre, the trajectory of the centre of mass of the fluid medium is analysed depending on the dimensionless parameters. A focus on the oscillation regime and on the equilibrium position of the centre of mass is presented. The forces acting on the solid frontier at the steady state are also investigated in this section. Finally, § 6 depicts the bubble shape evolution through time. The capillary wave and the first bounces of the bubble against the tank wall are studied and the equilibrium shape of the bubble is described.

Formalism and description of the tank modelling

We propose in this section, a description of the satellite tank configuration and the manoeuvres we enforce. The analysis of the dimensionless Navier-Stokes equations emphasises the dominating effects and physical assumptions are justified.

Dimensionless Navier-Stokes equations

Consider a domain Ω = Ω f ∪ Ω s with boundary ∂Ω. The regions Ω f and Ω s represent the fluid and solid regions, respectively. The boundary between the fluid and the solid, denoted by Γ s , is a sphere of diameter D t which symbolises the tank wall. The incompressible Navier-Stokes equations for Newtonian single-phase flows can be written

∇ • u = 0 in Ω f , (2.1) ρ ∂u ∂t + (u • ∇)u = -∇p + ∇ • (2µD) + F ext in Ω f , (2.2)
where t is the time, ρ the fluid density, u = (u, v, w) the velocity field, µ the fluid viscosity, p the pressure, F ext the volume forces induced by the motion of the satellite and D is the rate-of-deformation tensor defined as

D = ∇u + ∇u T 2 . (2.3)
We now divide the fluid domain Ω f into two parts, a liquid and a gas region denoted by Ω l and Ω g , respectively. (We consider a liquid region and a gas region for the sake of presentation, but this work is equally applicable to two liquids.) The interface between the two fluids is called Γ and n is its outward normal. Likewise, n s corresponds to the outward normal of the solid's boundary Γ s . Initially, the liquid-gas interface is a sphere, the liquid is in contact with the tank wall and the gas forms a spherical bubble of diameter D 0 in the centre (see figure 1). The density and viscosity fields are discontinuous at the interface and the effect of the surface tension σ must be considered.

In order to understand the prominent physical effects involved in this configuration, we present here the expressions of the external forces acting in the Navier-Stokes equations (2.1) and (2.2). The equations are solved in the frame of the tank which is a non-Galilean referential when the satellite is rotating. Figure 1 shows that satellite tanks are usually not aligned with the rotation axis. As a consequence, external forces, induced by the rotation of the tank around the rotation axis, such as inertial and Coriolis effects (although the latter forces are weak in the configurations discussed in this study), are taken into account. The Coriolis acceleration a 1 is expressed as follows,

a 1 = 2Ω × u, (2.4)
where Ω = Ωe z is the angular velocity vector of the satellite. The inertial acceleration a ce is the sum of the Euler acceleration a 2 , the centrifugal acceleration a 3 and the acceleration a R between the centre of mass of the satellite and the tank.

a ce = a 2 + a 3 + a R with a 2 = Ω × OM and a 3 = Ω × (Ω × OM), (2.5a,b)
where Ω = Ωe z the angular acceleration vector, OM the position vector of the fluid particle, defined as OM = xe x + (L + y)e y + ze z with L the lever arm (distance between the centre of the tank and the axis of rotation). As the tank is fixed to the satellite, a R is null. After defining the forces induced by the change of referential, we give now the expression of the Navier-Stokes equations in the frame of a rotating tank,

∇ • u = 0, (2.6) ρ ∂u ∂t + (u • ∇)u = -∇p + ∇ • (2µD) + ρ(a 1 + a ce ).
(2.7)

To define dimensionless Navier-Stokes equations in our configuration, the latter two expressions can be specified with the dimensionless variables,

ρ * = ρ ρ l , µ * = µ µ l , p * = pD 0 σ , u * = u ρ l D 0 σ , t * = t σ ρ l D 3 0 , (2.8a-e) ∇ * = D 0 ∇, a * 1 = a 1 γ , a * 2 = a 2 γ a , a * 3 = a 3 γ , (2.9a-d)
with γ and γ a characteristic accelerations of the flow that are defined below. Then the Navier-Stokes equations become

∇ * • u * = 0, (2.10) ρ * ∂u * ∂t * + (u * • ∇ * )u * = -∇ * p * + Oh∇ * • (2µ * D * ) + Bo i ρ * a * 2 + Bo c ρ * (a * 1 + a * 3 ),
(2.11)

with Oh the Ohnesorge number and two Bond numbers, defined as

Oh = µ 2 l ρ l σ D 0 1/2 , Bo c = ρ l γ D 2 0 σ , Bo i = ρ l γ a D 2 0 σ . (2.12a-c)
The Bond numbers compare the inertial effects induced by the satellite motion to surface tension effects, and the Ohnesorge number compares viscosity to surface tension effects. The two different accelerations, γ and γ a , are respectively based on the angular velocity Ω, such as γ = Ω 2 L, and the angular acceleration γ a = ΩL. As a result, the following two Bond numbers can be defined to characterise our physical system,

Bo c = ρ l Ω 2 LD 2 0 σ , Bo i = ρ l ΩLD 2 0 σ , (2.13a,b)
where Bo c is the Bond number related to the angular velocity responsible of centrifugal forces, and Bo i the Bond number depending on the angular acceleration.

Description of the tank and the manoeuvres

The spherical tank has a diameter of D t = 0.585 m and is not aligned with the rotation axis of the satellite. The lever arm, for which the value is L = 0.4 m, is the distance between the centre of the tank and the axis of rotation. The properties of the fluids included in the tank (close to the actual properties of fluids used in space applications) are ρ l = 1004 kg m -3 , µ l = 9.13 × 10 -4 kg m -1 s -1 for the liquid phase, and ρ g = 2.41 kg m -3 , µ g = 1.99 × 10 -5 kg m -1 s -1 for the gas phase. The surface tension is equal to σ = 0.03325 N m -1 . The bubble is initially located at the centre of the tank. As the wall of the tank is perfectly wet, the static contact angle is set to θ s = 0 • , preventing the formation of a contact line on the solid frontier. This assumption is justified in § 2.3.

The manoeuvres consist in a rotation of the satellite around its vertical axis e z . They can be split in two steps, first a constant acceleration phase, Ω = Cste = 0 and Ω(t) = Ωt until the final angular velocity Ω is reached. In the second step, the angular acceleration is stopped Ω = 0, and the final angular velocity is maintained to its final value Ω. These angular acceleration and final velocity are respectively related to the impulsional and centrifugal Bond numbers defined in (2.13). Figure 2 illustrates the evolution of Bo c and Bo i throughout the manoeuvre. The figure on the left presents two configurations with the same Bo i to reach different values of the final Bo c , whereas the figure on the right shows two different Bo i to obtain the same Bo c in the final step of the manoeuvre.

The range of the Bond numbers are described in tables 1 and 2 in the case of a half-filled tank. The corresponding angular accelerations and velocities are low values, which correspond to low inertial manoeuvres.

Assumptions of the physical model

In this section, we justify two important assumptions that have been made in our physical model, i.e. the fluids' incompressibility and the zero contact angle condition. First, the effect of gas compressibility is investigated through the Mach number. We define two characteristic velocities v 1 and v 2 which depend on the angular velocity Ω and the angular acceleration Ω imposed during the manoeuvre of the tank,

v 1 = ΩL and v 2 = ΩL. (2.14a,b)
Considering helium gas, the velocity of sound c he is computed using the ideal gas law c he = √ γ RT/M with γ = 1.4 the adiabatic index, R = 8.31 J mol -1 K -1 the molar gas constant, M = 4 g mol -1 the molar mass of the helium and T ≈ 300 K the absolute temperature. Finally, the following two Mach numbers can be considered

Ma 1 = v 1 c he = ΩL M γ RT , (2.15) Ma 2 = v 2 c he = L ΩM γ RT . (2.16)
Given the inertial parameters of tables 1 and 2, the Mach numbers always respect Ma 1 < 2.2 × 10 -5 and Ma 2 < 2.2 × 10 -5 . Therefore, we neglect the gas compressibility and consider the incompressible Navier-Stokes equations in the fluid phase.

We introduce now a discussion on the no contact angle assumption which is made in our simulations. A contact angle between a solid and a liquid-gas interface results from a force balance between the surface tension forces between the liquid-gas, liquid-solid and the solid-gas interfaces. If the surface tension characterising the solid-fluid interactions is not always well known [START_REF] Langbein | Capillary Surfaces: Shape Stability Dynamics[END_REF], indirect measurements of the difference between these two surface tension coefficients can be obtained by measuring the contact angle of a bubble or a droplet on a substrate.

As usual propellants are perfectly wetting fluids, we impose in our simulations a perfectly wetting condition, i.e. the contact angle is set equal to zero, when the liquidgas interface will be close to the tank wall. This condition prevents the formation of a contact line at the frontier between the liquid, the gas and the solid phase. As a result a thin liquid film is still entrapped between the solid wall and the liquidgas interface. Consequently, there is no hydrodynamical singularity on the viscous dissipation in this specific configuration since there is no contact line, and thus no contact line velocity. This no contact line assumption seems reasonable within the limits of low Bond numbers which is consistent with the present study. Indeed, it is clear that even if a perfectly wetting fluid is considered, in high Bond number configurations, the thin liquid film can be disrupted by strong external forces, such as gravity in ground applications. However, to the best of our knowledge, the critical Bond number for which the liquid film disruption will be observed is not known. In fact, an indirect validation of the no contact line assumption is presented in the rest of the paper by comparing our numerical simulations with experimental results from fluidics experiments both on the force acting on the tank wall and more specifically on the torque exerted by the fluid motion on the tank.

Finally, it should be noticed here that an apparent contact angle can also result from boiling phenomena, even if perfectly wetting fluids are involved, see for instance [START_REF] Stephan | Analysis of the heat transfer coefficient of grooved heat pipe evaporator walls[END_REF], [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] for more details. However as we consider here isothermal conditions, this physical effect is not relevant in the present study.

Numerical methods

Our numerical simulations have been performed by using the home-made code DIVA (dynamics of interface for vaporisation and atomisation) which is based on a wide range of numerical methods specifically dedicated to the computation of two-phase flows. In particular, the level set method [START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flow[END_REF] and the ghost fluid method [START_REF] Fedkiw | A non-oscillatory Eulerian approach to interfaces in multimaterial flows (the ghost fluid method)[END_REF]Kang, Fedkiw & Liu 2000;Liu, Fedkiw & Kang 2000;[START_REF] Sussman | A sharp interface method for incompressible two-phase flows[END_REF]) are used jointly to ensure a sharp and accurate representation of the liquid-gas interface. Computations in complex geometry can be performed by using the numerical scheme proposed in [START_REF] Ng | An efficient fluid-solid coupling algorithm for single-phase flows[END_REF] for single-phase flows and for which an extension to two-phase flows is proposed in [START_REF] Lepilliez | On two-phase flow solvers in irregular domains with contact line[END_REF]. This numerical scheme is based on a second-order numerical scheme to impose the Neumann boundary condition on the pressure field at the solid frontier while maintaining a sharp representation of the solid-fluid interface. Our numerical tool also allows computations of liquid-vapour phase change following the numerical methods described in [START_REF] Gibou | A second-order-accurate symmetric discretization of the poisson equation on irregular domains[END_REF][START_REF] Gibou | A level set based sharp interface method for the multiphase incompressible Navier-Stokes equations with phase change[END_REF], [START_REF] Tanguy | A level set method for vaporizing two-phase flows[END_REF], [START_REF] Papac | Efficient symmetric discretization for the Poisson, heat and Stefan-type problems with Robin boundary conditions[END_REF], [START_REF] Tanguy | Benchmarks and numerical methods for the simulation of boiling flows[END_REF], Rueda [START_REF] Villegas | A ghost fluid/level set method for boiling flows and liquid evaporation: application to the Leidenfrost effect[END_REF] and interaction of acoustic waves with bubbles [START_REF] Huber | A time splitting projection scheme for compressible two-phase flows. Application to the interaction of bubbles with ultrasound waves[END_REF]. It has been validated by direct comparisons with a wide range of theoretical studies [START_REF] Tanguy | Benchmarks and numerical methods for the simulation of boiling flows[END_REF][START_REF] Huber | A time splitting projection scheme for compressible two-phase flows. Application to the interaction of bubbles with ultrasound waves[END_REF]Lalanne et al. 2015b;[START_REF] Villegas | A ghost fluid/level set method for boiling flows and liquid evaporation: application to the Leidenfrost effect[END_REF] for isothermal two-phase flows and configurations involving liquid-vapour phase change. This section describe briefly the formalism and the numerical methods used to perform multiphase flows simulations by solving the Navier-Stokes equations, and imposing the suitable boundary conditions on the solid frontier when an irregular domain is involved.

Interface tracking method

In this work, 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, both to capture the interface motion between the two fluids and to impose the static boundary condition between the fluids and the solid. Signed and continuous level set functions φ and φ s are defined in the domain Ω respectively such as Ω g = {x : φ(x) < 0}, Ω l = {x : φ(x) > 0} and Γ = {x : φ(x) = 0} where Γ represents the fluid-fluid interface and Ω s = {x : φ s (x) > 0}, Ω f = {x : φ s (x) < 0} and Γ s = {x : φ s (x) = 0}, the fluid-solid interface. The fluid-solid interface does not move because the simulation is done in the referential of the tank. The motion of the liquid-gas interface Γ is captured by solving the following transport equation ∂φ ∂t

+ u • ∇φ = 0. (3.1)
The level set method benefits from attractive regularity properties in the whole domain that allow a simple, stable and efficient computation of the interface curvature. In particular, the signed distance function property can be imposed at the initial step of the computation and maintained all along the temporal integration process by applying the reinitialisation algorithm proposed in [START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flow[END_REF] at the end of every physical time step. It consists in solving the following partial differential equation over a fictitious time τ , ∂d ∂τ

+ sign(φ)(1 -∇d ) = 0. (3.2)
Next, geometrical properties such as n the outward unit normal vector to the interface and κ the local mean curvature can be computed by a simple differencing of the level set function φ. Spatial derivatives in the transport equation and the reinitialisation equation are computed with the weighted essentially non-oscillatory scheme [START_REF] Borges | An improved weighted essentially non-oscillatory scheme for hyperbolic conservation laws[END_REF], and the temporal derivatives with a third-order total variation diminishing Runge-Kutta scheme. Inside the solid domain, the following boundary condition is imposed on the level set function,

n s • ∇φ| Γ = cos θ s , (3.3)
where θ s is the contact angle formed on the triple line between the liquid, the gas and the solid. Let us note here that θ s = 0 • and thus, the right-hand side is equal to 1. Based on the previous works of [START_REF] Fedkiw | A non-oscillatory Eulerian approach to interfaces in multimaterial flows (the ghost fluid method)[END_REF], this boundary condition can be imposed by extending the level set function φ inside the solid domain, following the approach proposed in [START_REF] Lepilliez | On two-phase flow solvers in irregular domains with contact line[END_REF].

The motion of the liquid-gas interface occurs at the beginning of each time step and allows us to update the viscosity and density fields {µ n+1

; ρ n+1 }. They can be expressed with the following expressions

ρ(φ) = ρ l ρ g ρ l + 1 - ρ g ρ l H(φ) , (3.4) µ(φ) = µ l µ g µ l + 1 - µ g µ l H(φ) , (3.5)
where H(φ) is the Heaviside function, defined as H(φ) = 0 in the gas and H(φ) = 1 in the liquid. For cells crossed by the interface, a sharp interpolation of the density and the viscosity is used following the general guidelines proposed in [START_REF] Liu | A boundary condition capturing method for poissons equation on irregular domain[END_REF] and [START_REF] Sussman | A sharp interface method for incompressible two-phase flows[END_REF].

Projection method

The incompressible Navier-Stokes equations for Newtonian two-phase flows can be described in different ways, whether the jump conditions between the two fluids are imposed as smoothed singular functions, or as sharp jump conditions locally applied to the density, viscosity and pressure field. In what follows, we briefly present one of the formulations for two-phase flows, adapted from the projection method of [START_REF] Chorin | A numerical method for solving incompressible viscous flow problems[END_REF] and first introduced by [START_REF] Sussman | A sharp interface method for incompressible two-phase flows[END_REF]. Some other formulations are thoroughly detailed and compared by [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] and Lalanne et al. (2015b). The authors demonstrate that the following formulation is best suited for our study. The implicit temporal discretisation, referred to as the ghost fluid conservative method with an implicit scheme (GFCMI) in [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] is used in this study. First, given a velocity u n at time t n = n t, an auxiliary field u * can be computed for a time step t without considering the pressure component

ρ n+1 u * -t∇ • (2µ n+1 D * ) = ρ n+1 u n -t (u n • ∇)u n - F ext ρ n+1 , (3.6)
which leads to a large linear system where the three velocity components are coupled.

As it is a diagonally dominant linear system, it can be simply and quickly solved (approximately 10 iterations or less in most configurations) by using a Gauss-Seidel algorithm.

Then the pressure field p n+1 serves as the scalar potential function of the Hodge decomposition, satisfying a Poisson equation with homogeneous Neumann boundary conditions on ∂Ω and non-homogeneous Neumann boundary condition on

Γ s ∇ • ∇p n+1 ρ n+1 = ∇ • u * t + ∇ • σ κnδ Γ ρ n+1 , (3.7)
where κ and δ Γ are respectively the local curvature of the interface and a multidimensional Dirac distribution located at the interface. It is well known that the resolution of the Poisson equation is the most expensive step in numerical simulations of two-phase flows. Consequently, to speed up our simulations, we have implemented the black box multigrid method from [START_REF] Dendy | Black box multigrid[END_REF] and MacLachlan, Tang & Vuik (2008) which allows a fast resolution of pseudo-Poisson equations that contain irregular terms either in the matrix coefficients (density jump) or on the right-hand side (pressure jump) of the linear system. These singular source terms are approximated using the general schemes developed in [START_REF] Liu | A boundary condition capturing method for poissons equation on irregular domain[END_REF] for maintaining a sharp representation of the liquid-gas interface. Finally, the fluid velocity u n+1 is defined at the new time step t n+1 as the projection of u * onto the divergence-free space

u n+1 = u * - t ρ n+1 (∇p n+1
σ κnδ Γ ).

(3.8) 3.3. Numerical methods for the fluid-solid boundary As our simulations imply a solid irregular frontier immersed in the computational domain, the velocity field must respect the impermeability and no-slip conditions

n s • ∇p ρ Γ s = n s • (u * -u s + σ κnδ Γ )| Γ s , (3.9)
with u s the solid velocity.

Comparison with the Fluidics experiment data

Very few experimental data exist for the study of low Bond number manoeuvres in a microgravity environment because of the long time period required that cannot be achieved with drop tower or parabolic flight experiments. In 2017, the Fluidics experiment was sent by the CNES, the French national space agency, to the International Space Station (ISS) to quantify the sloshing phenomenon inside a spherical tank [START_REF] Mignot | Fluid dynamic in space experiment[END_REF]). To validate the numerical code, we implement the case of the Fluidics experiment and compare our numerical results with the experimental data.

Description of the Fluidics experiment

The experiment, represented in figure 3, consists of a spherical transparent tank containing a safe substitute of propellant and helium, connected by an arm to a motor in the centre. The global system is integrated into a metallic housing and fixed on the seat rack rails of the ISS. The tank is subjected to angular velocity profiles similar to figure 2 followed by a deceleration and a zero velocity phase. Forces and torques generated on the tank are measured by sensors and bubble deformations are observed thanks to cameras. Two tanks have been sent to the ISS with a filling ratio of respectively 50 % and 75 %. The same angular acceleration Ω and velocity Ω are enforced on the tanks. Equation (2.13) returns the Bond numbers of each tank as Bo i = 23 and Bo c = 76 for the half-filled tank and Bo i = 15 and Bo c = 48 for the tank with a filling ratio of 75 %. The values of the Bond numbers are higher than those of § 2.2 because of the noise of the experimental data. Indeed, the data from the force and torque sensor contain an important level of noise which prevents the data recovery from low Bond numbers manoeuvres.

Space convergence

The simulations of the Fluidics case are computed with a mesh of 128 cells in the 3 directions. To prove the space convergence, we consider the half-filled tank because the Bond numbers, Bo i = 23 and Bo c = 76, are higher than those of the tank with τ fill = 75 %. The evolution of the force in the y-direction and the torque in the z-direction through time are shown in figure 4, for three different meshes, 64, 128 and 256 nodes in each direction. We observe that the global evolution of the force through time is well predicted by the three meshes. The average force during the constant angular velocity phase remains identical but the magnitude and the frequency of the oscillations depend weakly on the mesh size. A finer grid leads to slightly stronger oscillations around the mean value. The same behaviour appears for the evolution of the z-direction torque. Finally, we can conclude that the two finer grids give results that are close enough to allow us to compare the experimental data with the results from the mesh with 128 cells in each direction. The force in the y-direction, denoted by F y , increases quickly at the beginning because of the angular acceleration. When the constant angular velocity phase is reached, the force oscillates around a constant value. This value corresponds to the centrifugal force exerted by the liquid on the tank wall. The average value of the dimensionless force is higher for the tank with τ fill = 75 % because more liquid is subjected to the centrifugal acceleration. The force evolution with the filling ratio is discussed in § 5.6. Moreover, the damping of the oscillations of the force is more important in the τ fill = 75 % tank. Indeed, the gas bubble is smaller in this tank and more constrained by the liquid against the tank wall. Finally, the angular velocity decreases until it reaches zero leading to the cancellation of the force on the tank wall.

The global evolution of the force achieved with our simulations is similar to the force measured under microgravity conditions despite the important noise that remains in the experimental signal. Table 3 lists the average values of the dimensionless force during the constant angular velocity phase obtained with the experimental data, F exp y /σ D t , and our simulations, F num y /σ D t , for both filling ratio cases. These values are close together with an error of 1.8 % for the half-filled tank and 3.8 % for τ fill = 75 % which consolidates the validity of the numerical simulations. The oscillations are difficult to compare because of the noise of the experimental data. Nevertheless, the first oscillations of the half-filled tank case, which magnitudes are higher than the level noise, are well predicted by the numerical simulations.

Study of the torque exerted on the tank wall

The torque data are also compared to our numerical results. The sensor is located at the interface between the rotating arm and the tank. The torque data allow us to validate the force exerted by the fluids on the tank wall and the position of the centre of mass. The torque T is integrated in our simulations as

T = Ω f SM × (∇p -σ κnδ Γ ) dV. (4.2)
with SM the vector between the sensor position and the fluid particle.

Figure 6 compares the torque in the z-direction from the Fluidics data to our numerical results for the two tanks cases. During the acceleration phase, we observe a positive peak of the z-direction torque which is due to the component of the acceleration in the x-direction depending on Bo i in (2.7). We observe the same peak with the opposite sign during the deceleration phase, at t * = 11.5 for the half-filled tank and at t * = 16 for the tank with τ fill = 75 %. During the constant velocity phase, the z-direction torque oscillates around zero. It represents the oscillation of the bubble around its equilibrium position against the tank wall. This bubble behaviour generates forces in the x-direction which appear in the z-direction torque. The numerical simulation reproduces well the torque oscillations until the damping brings their magnitudes around the noise level. The damping of the oscillation is more important with the higher filling ratio because more liquid is subjected to the centrifugal acceleration. Finally, during the zero velocity phase at the end, the torque remains zero. The peaks and the oscillations for which the magnitude is higher than the noise level are well predicted by the numerical simulations.

As the torque depends both on the force exerted on the tank wall and on the centre of mass of the liquid, these results provide a strong validation for these variables. As a result, an indirect validation of the zero contact angle hypothesis is obtained here since the centre of mass depends mainly on the bubble volume and on the interface shape and thus on the contact angle value.

Evolution of the centre of mass

This section investigates the evolution of the centre of mass of the fluid medium under the effect of the inertial forces. First, a grid sensitivity study is presented. Then, an overview of the physical phenomena is given through an example of the fluids behaviour during a manoeuvre. The temporal evolution of the centre of mass is described for a wide range of the Bond numbers. For some of the higher values of the dimensionless parameters, an oscillatory phenomenon of the centre of mass around its equilibrium position is observed and discussed. Finally, the position of the centre of mass and the force exerted on the tank wall at the steady state are described and fitted by semi-empirical laws.

The subsections focusing on the transient regime consider the filling ratio of 50 % and the Bond numbers summarised in tables 1 and 2. The two last subsections aim the steady state, which does not depend on Bo i , and investigate the influence of the filling ratio. The filling ratios considered in the steady state are listed in table 4.

Grid sensitivity study

We present here numerical evidence about the grid sensitivity of the numerical simulations. All the simulations presented in this study have been computed with a mesh of 128 cells in each direction. In figure 7, the evolution of the coordinates of the centre of mass of the fluids throughout time has been plotted for three different meshes containing respectively 64, 128 and 256 nodes in each direction. The z-coordinate of the centre of mass does not appear in figure 7 because no external acceleration happens in the z-direction and therefore, the bubble stays centred in the z-direction. The z-coordinate of the centre of mass obtained in our numerical results is always lower than the mesh size. These results correspond to the case with Bo c = 5.80, Bo i = 5.80 and τ fill = 50 % implying the highest velocity of the centre of mass and the highest deformation of the bubble. These comparisons show a good agreement between the solutions on the latter two grids, from which we assume that the intermediate grid 128 3 is thin enough considering the spatial convergence of our computations.

An example of bubble motion

This section describes the motion of the fluids inside the tank in a specific case. It illustrates the global behaviour of the bubble gas, and consequently the centre of mass evolution, during a manoeuvre. We consider the half-filled tank case with Bo i = 1.16 and Bo c = 0.116. The first plot in figure 8 shows the trajectory of the centre of mass in the xy-plane. The other plots in figure 8 are pictures of the position of the fluids inside the tank at 6 specific times and are described below. The dark blue interface represents the gasliquid interface and the grey represents the fluid-solid interface. In addition, movie 1 (available at https://doi.org/10.1017/jfm.2018.389) displays the accelerated motion of the bubble during the manoeuvre.

Initially, the bubble is spherical and located in the centre of the tank (a). At the beginning of the manoeuvre, the bubble goes toward the axis -e x until it reaches the tank wall (b). The bubble spreads and bounces (c), moves in the direction -e y and spreads on the wall a second time (d). Finally, the bubble follows the tank wall (e) until it reaches its final position around the x-axis ( f ). The motion of the centre of mass is the opposite of the global motion of the bubble. The first bouncing of the bubble leads to the loop in the trajectory of the centre of mass. In all the pictures in figure 8, a thin film of liquid is always trapped between the gas bubble and the tank wall because of the perfectly wetting condition.

The detailed analysis of the trajectories for different Bond numbers is the aim of the next subsection.

Trajectory of the centre of mass

The trajectory of the centre of mass of the fluids depends on the temporal evolution of the inertial acceleration a i = a 1 + a 2 + a 3 , expressed here as a function of the two main dimensionless numbers Bo i and Bo c characterising our physical system,

a i = σ ρ l D 2 0 L [Bo i (-(L + y)e x + xe y ) -Bo c (xe x + (L + y)e y )].
(5.1)

The trajectories are presented in figure 9. The evolution of each coordinate of the centre of mass throughout time is given in the Appendix.

At early flow times, Bo i dominates the flow, regardless of the final value of the Bo c . As predicted by (5.1), the centre of mass moves along the opposite direction to the acceleration as L x,

a i - σ ρ l D 2 0 L Bo i (L + y)e x .
(5.2)

In the case where the rotational acceleration and the rotational velocity are small, the acceleration is stopped before the bubble reaches the tank wall. The rest of the flow entirely relies on Bo c . This can be observed in figure 9(a,b), where for a given Bo c , there is no distinction between the different values of Bo i . After the end of the acceleration phase, the rest of the trajectory can be predicted by

a i = - σ ρ l D 2 0 L Bo c (xe x + (L + y)e y ).
(5.

3)

The deformation of the bubble during the first contact is described in § 6.2. Then the bubble stays close to the tank wall due to the xe x term, to finally reach its final equilibrium value (where x cm (t eq ) = 0, and y cm (t eq ) is discussed in § 5.5).

In figure 9(c), the trajectories are similar, except for the lowest value of Bo i . In this case, the final value of Bo c is not reached before the bubble touches the tank wall. After the bounce, the terms corresponding to the angular acceleration are still not equal to zero and the centrifugal acceleration is more important in the direction -e x . The bubble tends to move faster to the tank wall in this direction. The parameter γ = Bo i /Bo c is defined to determine if the angular acceleration phase is ended before the bubble first spreads on the wall. We report that the bubble is not subjected to any Bo i contribution after the spreading of the bubble against the tank wall for γ > γ lim = 2L/D t . This limiting value γ lim is obtained geometrically considering the simplified motion of the centre of mass at early flow time. This result predicts that all trajectories are similar if γ > γ lim . Actually, if the acceleration phase ends before the bubble reaches the tank wall, the angular velocity during the spreading is constant and identical for a given Bo c . This means that the resulting trajectory of the bubble remains similar. Yet, another phenomenon can disrupt the accuracy of this predictive parameter: when the acceleration is stopped, a capillary wave on the bubble interface is observed. This is discussed in § 6.1 on the early motion of the bubble.

In figure 9(d-g), it is shown, for a given Bo c , that the trajectories evolve depending on the Bo i . For higher values of Bo c , more potential energy is transmitted through the centrifugal acceleration to the system. As a result, for Bo c > 0.58, the centre of [0.058; 0.116; 0.290; 0.580; 1.16; 2.90; 5.80] from (a-g).

mass of the fluids goes beyond the equilibrium point in the x-direction. Indeed, the formulation of the inertial acceleration at this time of the flow is given by (5.3). If the centre of mass goes beyond x = 0, the inertial acceleration tends to correct it by forcing the centre of mass in the opposite direction (if x cm < 0, a i • e x > 0, and if x cm > 0, a i • e x < 0). This correction of the inertial acceleration is the source of an oscillatory phenomenon, discussed in § 5.4. 5.4. Oscillations of the centre of mass of the fluids upon reaching equilibrium state Figure 9 exhibits that for high Bo c and for all the values of Bo i , an oscillatory motion appears around the x-axis before the system stabilises. Upon reaching its equilibrium state, when the bubble goes beyond x = 0, the inertial acceleration corrects the trajectory: if x cm < 0, a i • e x > 0 and if x cm > 0, a i • e x < 0. This results in an oscillatory motion of the centre of mass, as observed in figure 10. Features of the bubble inside the tank at three different times are shown in figure 11.

The oscillatory phenomenon can be approached by the analytical model of a onedimensional mass-spring system with a damping term described by (5.4) with ω 0 the natural pulsation of the mass-spring system, and ξ a damping coefficient. The solution of this equation is

∂ 2 x ∂t 2 + 2ξ ω 0 ∂x ∂t + ω 2 0 x = 0,
x(t) = A 0 exp(-ξ ω 0 t) sin(ωt + Φ),
(5.5)

with ω = 1 -ξ 2 ω 0 , A 0 the initial value of the amplitude of the oscillation, and Φ the phase. Only the results for high Bo c are presented because periods are very long and are not entirely captured in the simulations with Bo c < 1. Moreover, the study is limited to cases for Bo i > 1 because the observed oscillations are not related to the angular acceleration phase. The parameters ξ and ω 0 are deduced from the fitting of our numerical results with (5.5) using a nonlinear least squares method and are reported in table 5. An example of the oscillation phase is circled in red in figure 10. The oscillations produced during manoeuvres with other Bond numbers are shown in figure 19(e-g) of the Appendix. It corresponds to the trajectory of the centre of mass of the fluids at the final stage of the manoeuvre when it oscillates around an equilibrium state. The two first periods of oscillations are not considered in the comparison between the theoretical model and the numerical results, because it corresponds to the early time of the manoeuvre when the bubble motion, the bouncing and the capillary waves participate in the motion.

Table 5 shows that when Bo c increases, the natural pulsation and the damping ratio also increase. Indeed, when the centrifugal acceleration is important, the bubble spreads on the wall with a high magnitude and is maintained against the wall. We observe that Bo i does not affect ω 0 . In the same way, the damping ratio stays relatively close to the same value regardless of Bo i . This value depends essentially on the Bo c .

Equilibrium position of the centre of mass

As the equilibrium state is reached, the resulting position of the centre of mass of the fluids corresponds to x cm = 0 and differs in the y-direction for different Bo c and filling ratios, as illustrated by the markers in figure 12 and plots of the Appendix. The impulsional Bond number Bo i become irrelevant in the steady state.

The y-position of the centre of mass increases with the angular velocity until it approaches an asymptotic value because of the tank wall. The evolution of y-position with the Bond number has been fitted with several functions: polynomial, exponential and logarithmic functions. We show that the best candidate is the logarithmic law. Indeed, at low Bond number, the evolution of the centre of mass is highly impacted by the Bond number because the inertial effects are lower than the surface tension effects. At higher Bond numbers, we are close to a saturation phenomenon due to the tank geometry and this leads to a weaker effect of the Bond number on the bubble shape. It is more complicated to spread an already deformed bubble. This saturation phenomenon is well depicted by the logarithmic function. On the contrary, the linear dependency on the filling ratio is clear in figure 12. It reveals that for each Bond number, the amount of liquid affects proportionally the position of the centre of mass. The position of the centre of mass in the y-direction is denoted by y cm and estimated by

y cm D t (t eq ) = C 1 (1 -τ fill ) ln(C 2 Bo c + 1), (5.6)
where τ fill = V liq /V tank is the filling ratio of the tank, C 1 and C 2 correspond to constants obtained with a least squares method: C 1 = 2.62 × 10 -2 and C 2 = 3.35 × 10 4 in our cases. The law is plotted with dashed lines in figure 12 and gives a good approximation of the numerical results.

5.6. Force generated on the tank wall at the equilibrium state The theoretical force F th (t eq ) corresponds to the reaction to the centrifugal acceleration at the equilibrium state, F th (t eq ) = -V t ρa i (t eq ) dV ρ l Ω 2 (L + y cm (t eq ))V l e y ,

(5.7) with a i the sum of all the external accelerations, V l the volume of liquid inside the tank by neglecting ρ g against ρ l . The volume of liquid is deduced from the filling ratio and the tank diameter as V l = τ fill V t = τ fill (π/6)D 3 t . We can write the dimensionless force with the parameters Bo c and τ fill ,

F th σ D t (t eq ) π 6 Bo c 1 + y cm (t eq ) L τ fill (1 -τ fill ) 2/3 e y , (5.8)
with D 0 = D t (1 -τ fill ) 1/3 . Section 5.5 exhibits that the final position of the centre of mass of the fluids y cm (t eq ) depends on Bo c and τ fill . We must consider the data in figure 12 or the approximation of (5.6) to obtain the force generated on the tank wall. In our numerical simulations, the integration of the pressure inside the tank (see (4.1)), allows deducing the value of the force in the y-direction at the equilibrium state. This computed force and the theoretical force of (5.8) are plotted, respectively with symbols and dashed lines, in figure 13.

We observe that the force evolves linearly with the Bond number and as a power law with the filling ratio. Indeed, the centrifugal force depends directly on the square of the angular velocity and the mass of the fluid inside the tank. Finally, the approximations given by (5.6) and (5.8) give really close results with the computed force. FIGURE 13. Evolution of the y-direction generated force, for different Bo c and filling ratios.

Evolution of the bubble shape

This section is focused on the bubble shape, from the initial sphere located at the centre of the tank until its equilibrium shape spread against the tank wall. First, the capillary wave, observed for some configurations, is described through the variation of the bubble surface. Then, the bubble behaviour is analysed during the first contact with the tank wall for several pairs of Bond number. Finally, the description of the equilibrium shape of the bubble, as its spreading diameter and its thickness, is done for different values of Bo c and τ fill .

Similarly than § 5, the two first subsections consider a filling ratio of 50 % and the Bond numbers summarised in tables 1 and 2. The last subsection does not deal with Bo i and investigates the influence of the filling ratios listed in table 4.

Capillary wave

During the transient regime, different phenomena can be observed, such as a capillary wave on the liquid-gas interface when the angular acceleration Ω stops, or when the bubble reaches the tank wall and spreads on it. The deformation of the bubble is analysed by considering the temporal evolution of the dimensionless surface energy variation E s (t) defined by

E s (t) = E s (t) -E s (t 0 ) E s (t 0 ) . (6.1)
For low values of Bo c and Bo i , a first variation of the surface energy is observed before the bubble reaches the tank wall, as depicted in figure 14(a), whereas it does not appear for higher Bo c values, as in figure 14(b). This variation is generated by Bo i evolving from its original value to zero as the final centrifugal velocity is reached. This time is denoted by t cw thereafter. The frequency of this capillary wave is equal to the Rayleigh frequency f R , defined by

f R = σ ρ l D 3 0 . (6.2)
Table 6 contains the amplitude of the variation of the surface energy due to capillary waves for different Bo c and Bo i , when measurable. From these data, we can conclude

Bo i = 0.116 - - - - - - - Bo i = 0.58 5.08 × 10 -5 1.75 × 10 -4 - - - - - Bo i = 1.16 6.46 × 10 -5 2.52 × 10 -4 1.81 × 10 -3 - - - - Bo i = 2.32
8.51 × 10 -5 2.88 × 10 -4 1.85 × 10 -3 ----Bo i = 5.80 9.12 × 10 -5 3.05 × 10 -4 1.88 × 10 -3 ----TABLE 6. Variation of the dimensionless surface energy at the capillary wave E s (t cw ) when measurable, and -otherwise.

that Bo c has a great influence on the amplitude of the capillary wave. Moreover, increasing Bo i tends to approach E s (t cw ) to an asymptotic value depending on Bo c . For the lower value of Bo i , the amplitude is too low to be measurable. This phenomenon also occurs for higher values of Bo c , probably with a stronger amplitude, but as the bubble has already reached the tank wall it cannot be accurately measured when γ < γ lim .

6.2. First bounce of the bubble When the bubble reaches the tank wall, its deformation leads to a sharp increase of the variation of the dimensionless surface energy of the bubble. The highest peaks in figure 14 correspond to several spreadings of the bubble on the tank wall. The first spreading on the wall has the greatest magnitude and, for all cases, the maximum of the surface energy variation at that time are presented in figure 15.

The dimensionless surface energy variation follows a linear evolution in cases for which the acceleration phase ends before the bubble reaches the tank wall for the first time. These configurations correspond to γ > γ lim and it can be approximated by Bo i = 0.116 7. Summary of early flow time phenomena: • for configurations when the capillary wave is clearly observed, + for γ < γ lim .

+ + + + + + Bo i = 0.58 • • + + + + Bo i = 1.16 • • • + + + Bo i = 2.32 • • • + + Bo i = 5.80 • • • + TABLE
In the same way as for the trajectories, Bo i does not seem to influence the spreading of the bubble on the tank wall when γ > γ lim . The constant velocity phase remains similar for a given Bo c which leads to a similar deformation on the wall. In the case where γ γ lim , figure 15 shows that the maximum of the dimensionless surface energy variation reaches a limiting value for each Bo i . Indeed, for a given centrifugal acceleration, the same angular velocity is reached during the spreading of the bubble and thus the same amount of momentum is introduced into the system. This explains that Bo c does not influence the spreading as long as the acceleration phase is in progress. A summary of the phenomena observed at early flow time is proposed in table 7.

Shape of the bubble at the equilibrium state

In this section, the final shape of the bubble is described and analysed, with the spreading diameter D b and the thickness of the bubble y b . Depending on the centrifugal Bond number Bo c and the filling ratio τ fill , the shape of the bubble can change, as observed in figure 16. We note that the third parameter Bo i has no influence in the equilibrium state. An increase of Bo c leads to an important spreading of the bubble on the tank wall and a shrinking of its thickness. Because of the variation of the filling ratio in our numerical simulations, these parameters are divided by the initial diameter of the bubble to become dimensionless D 0 = D t (1 -τ fill ) 1/3 .

The dimensionless thickness of the bubble is represented with symbols in figure 17 for different values of the filling ratio and the Bond number. First, for Bo c < 0.58, the thickness decreases with the Bond number linearly and is not modified by the filling ratio. This corresponds to the first regime when the surface tension is dominant and the deformation of the bubble is small. For Bo c > 0.58, the thickness of the bubble keeps decreasing with the Bond number but with a slower magnitude. This case, where the inertial effects must be considered, corresponds to the second regime.

As already seen in § 5.5, the thickness evolution has been approached with several functions and the logarithmic function has given the best to fit the data. The same saturation effect, like the position of the centre of mass, limits the influence of the Bond number on the bubble shape. The thickness tends toward an asymptotic value for high Bo c because of the size of the tank. This limit depends on the filling ratio as we can see on the left plot of figure 17. Indeed, when τ fill < 0.5, the bubble equilibrium thickness approaches its initial value. For τ fill > 0.5, the bubble is small compared to the tank and spreads with a higher magnitude on the wall. For each Bond number, the thickness of the bubble evolves linearly with the filling ratio. The dimensionless bubble thickness is approached by This equation allows us to predict the behaviour of the bubble thickness with the Bond number and filling ratio for the two regimes. Equation (6.4) is plotted with dashed lines in figure 17.

y b D 0 (t eq ) = 1 + (C 3 + C 4 τ fill ) ln(C 5 Bo c + 1), ( 6
Figure 18 presents the evolution of the dimensionless spreading diameter of the bubble against the tank, for different filling ratios and different Bo c . As for the bubble thickness, the variation of the spreading diameter D b /D 0 can be divided into two regimes depending on the Bond number. When Bo c < 1.16, the spreading diameter evolves linearly with the Bond number and remains independent from the filling ratio. For higher Bond numbers, the spreading diameter increases with a slower magnitude and we assume that it will reach an asymptotic value depending on the tank diameter and the filling ratio. This dependency on the Bond number, like that of the bubble thickness, is well predicted by the logarithmic function. The same low impact of the filling ratio is noticeable for Bo c = 2.902 but a local increase of the spreading diameter is observed for τ fill < 0.5. As for the bubble thickness, the limiting value is closer to the initial spreading diameter for a low filling ratio. By considering the regime independent of the filling ratio, the spreading diameter can be approached by

D b D 0 (t eq ) = 1 + C 6 ln(C 7 Bo c + 1), (6.5)
with C 6 = 8.25 × 10 -2 and C 7 = 4.34 × 10 -1 in our simulations and is depicted with dashed lines in figure 18.

Summary

We have discussed the motion and deformation of a bubble in a spherical tank rotating around a fixed axis in microgravity conditions. Based on threedimensional direct numerical simulations, our analysis is focused on low Bond number configurations for which theoretical models are not available due to the prominence of nonlinear surface tension effects. Very few experiments exist to predict the fluids behaviour because of the long time required in microgravity conditions. Direct numerical simulation is one of the promising ways to understand the fluid behaviour during low Bond number manoeuvres. Because of the low inertial accelerations we enforce on the tank, the gas compressibility has been neglected. We considered the perfectly wetting condition in our model because of the properties of the typical propellants used in the space industry. This assumption is reasonable for low inertial manoeuvres for which the Bond numbers stay lower than a critical value. Further experiments are required to determine this value. Nevertheless, the comparison with experimental data showed that considering the no contact line assumption predicts well the fluids' behaviour for higher values of the Bond numbers. Our numerical model have been validated through a detailed spatial convergence study. The data from the Fluidics experiment, carried out in the International Space Station, have been presented in a high Bond number case and compared to our results. The experiment is a low rate centrifuge which applies a rotation around the motor axis of a spherical tank containing a safe substitute of propellant and helium gas. Sensors measure the force and torque exerted by the fluids on the tank wall in microgravity conditions. We showed that the space convergence is still valid in the Fluidics experiment cases.

The comparison of the force and torque showed that we predict well the different phases of the manoeuvre. The centrifugal force exerted during the constant angular velocity phase is obtained with an error lower than 4 % for two different tanks. The torque peaks corresponding to the acceleration and deceleration phases are similar in both cases. The oscillations of the bubble around its equilibrium position have been investigated through the torque evolution. The magnitude and the frequency of the torque oscillations measured in the ISS comply with our numerical simulations and the same damping of the oscillations throughout time is observed. This shows that the simulations predict well the phenomena observed in space and validates the assumptions made. The evolution of the centre of mass has been investigated through a parametric study on two dimensionless Bond numbers respectively based on the angular velocity and the angular acceleration. We showed that the trajectory of the centre of mass is quite different depending on the final angular velocity. For low values of Bo c , the bubble goes toward the tank wall, bounces and follows the tank wall until it reaches its final position. We observed in this case that the trajectory does not depend on the angular acceleration. We defined a criterion γ = Bo i /Bo c > γ lim = 2L/D t which quantifies that there is no influence of Bo i on the trajectory of the centre of mass. For higher values of Bo c , the spreading and bouncing of the bubble against the tank wall is more important. Several bounces against the tank wall may happen and we observed a clear influence of the second Bond number Bo i . In this case, the bubble oscillates around its equilibrium position before reaching it. We approached the oscillations phase with the solution of a one-dimensional mass-spring system with a damping term. We obtained the natural pulsation of the oscillation and the damping coefficient for different values of the Bond numbers. We showed that the natural pulsation and the damping of the oscillations increase with Bo c . The equilibrium position of the centre of mass was investigated with Bo c and the filling ratio τ fill . Its position corresponds to x = 0 and the value y depends linearly on the filling ratio and follows a logarithmic law with the Bond number. Indeed, the amount of liquid affects linearly the position of the centre of mass and increasing the angular velocity tends to spread the bubble against the tank wall with a higher magnitude. As a consequence, spreading an already deformed bubble against a tank wall is more complicated. This saturation phenomenon is well predicted by a logarithmic law. We deduced the force generated on the tank wall from the position of the centre of mass. It evolves linearly with the Bond number and follows a power law with the filling ratio. In a last part, we were interested in the shape of the bubble throughout time. We observed at early flow time capillary waves for specific cases. This happens when the angular acceleration stops because the final angular velocity is reached. We observed a variation of the surface energy before the bubble reaches the tank wall for the first time in the case of low values of Bo c and Bo i . The amplitude of the variation increases greatly with Bo c and raising Bo i leads the amplitude to reach an asymptotic value. We also focused on the first spreading of the bubble because it corresponds to the higher deformation of the bubble throughout time. The surface variation increases with Bo c and we observed that this evolution follows a linear evolution with Bo c in the case γ > γ lim . For the other cases, the surface energy variation reaches a limiting value for each Bo i . The equilibrium shape of the bubble was investigated through its thickness and spreading diameter. We divided the range of Bo c in two regimes, the first one where the surface tension dominates and the deformation of the bubble is quite small and the second one where the inertial effects become more important. In the first regime, the bubble thickness decreases linearly with the Bond number and does not vary with the filling ratio. We observed the opposite effect on the spreading diameter which increases linearly with the Bond number and stays constant with the filling ratio. During the second regime, the bubble thickness decreases with a lower magnitude because the spreading becomes more important and the bubble comes closer to its asymptotic shape. This evolution is well predicted by a logarithmic law for both the influence of the spreading diameter and the bubble thickness with the Bond number.

  FIGURE 1. (Colour online) Typical positioning of tanks in a satellite and inside view of a tank.

  FIGURE 2. (Colour online) Examples of different velocity profiles, with the evolution of Ω throughout the manoeuvre.

FIGURE 3 .

 3 FIGURE 3. (Colour online) The Fluidics experiment.

FIGURE 4 .

 4 FIGURE 4. Evolution of the y-direction force and z-direction torque generated by the fluids in the half-filled tank for Bo c = 76 and Bo i = 23 and three different meshes: 64 3 , 128 3 and 256 3 .

4. 3 .

 3 Figure5confronts the dimensionless force in the y-direction measured by the sensor with the one obtained by computation during the manoeuvre on the two tanks.

FIGURE 5 .

 5 FIGURE 5. Comparison of the y-direction force from the experimental data and our numerical simulations with (a) τ fill = 50 % and (b) τ fill = 75 %.

FIGURE 6 .

 6 FIGURE 6. Comparison of the z-direction torque from the experimental data and our numerical simulations with (a) τ fill = 50 % and (b) τ fill = 75 %.

FIGURE 7 .

 7 FIGURE 7. Evolution of the coordinates of the centre of mass of the fluids for Bo c = 5.80, Bo i = 5.80, τ fill = 50 % and three different meshes: 64 3 , 128 3 and 256 3 .

τ

  FIGURE 8. (Colour online) Trajectory of the centre of mass of the tank in the xy-plane and bubble positions at 6 different times with Bo i = 1.16, Bo c = 0.116 and τ fill = 50 %.

FIGURE 9 .

 9 FIGURE 9. Trajectories of the centre of mass of the fluids in the xy-plane with τ fill = 50 %, for different Bo i and Bo c =[0.058; 0.116; 0.290; 0.580; 1.16; 2.90; 5.80] from (a-g).

  FIGURE 10. (Colour online) Position of the centre of mass of the fluids in the x-direction for Bo c = 5.8 and Bo i = 5.8. The oscillation phase is circled in red.

FIGURE 12 .

 12 FIGURE 12. Position of the centre of mass of the fluids in the y-direction for different Bo c and filling ratios.

  FIGURE 14. (Colour online) Evolution of the surface energy of the interface of the bubble for Bo i = 1.16, (a) Bo c = 0.116 with the capillary deformation circled in red and (b) Bo c = 0.58.

  E s (t b ) = 0.0145Bo c , (6.3) with t b is the time of the bouncing defined as max(E s (t)) = E s (t b ).

FIGURE 15 .

 15 FIGURE 15. Maximum variation of the surface energy of the bubble interface at the first bounce. Filled markers correspond to γ > γ lim and empty markers to γ γ lim .

FIGURE 17 .

 17 FIGURE 16. (Colour online) Final shape of the bubble for (a) Bo c = 0.116 and (b) Bo c = 5.8, with τ fill = 50 %.

  .4) with C 3 = -2.73 × 10 -2 , C 4 = -5.02 × 10 -2 and C 5 = 2.95 in our simulations.

FIGURE 18 .

 18 FIGURE 18. Evolution of the thickness of the bubble, for different Bo c and filling ratios.

FIGURE 19 .

 19 FIGURE 19. The x-coordinate of the centre of mass of the fluids throughout time with τ fill = 50 %, for different Bo i and Bo c =[0.058; 0.116; 0.290; 0.580; 1.16; 2.90; 5.80] from (a-g).

FIGURE 20 .

 20 FIGURE 20. The x-coordinate of the centre of mass of the fluids throughout time with τ fill = 50 %, for different Bo i and Bo c =[0.058; 0.116; 0.290; 0.580; 1.16; 2.90; 5.80] from (a-g).

TABLE 1 .

 1 Inertial parameters used for this study for τ fill = 50 %: Bo i and Ω.

	Bo c	0.058 0.116 0.29	0.58	1.16	2.90	5.80
	Ω (10 -2 rad s -1 ) 0.472 0.668 1.055 1.493 2.111 3.338 4.721

TABLE 2 .

 2 Inertial parameters used for this study for τ fill = 50 %: Bo c and Ω.

TABLE 3 .

 3 Comparison of the mean y-direction forces in both cases of the Fluidics experiment.

TABLE 5 .

 5 Characteristics of the oscillatory phenomenon.

		Bo c = 1.16	Bo c = 2.90	Bo c = 5.80
	Bo i	ω 0 (rad s -1 )	ξ	ω 0 (rad s -1 )	ξ	ω 0 (rad s -1 )	ξ
	1.16	0.0277	0.1161	0.0402	0.1286	0.0555	0.1312
	2.32	0.0272	0.1241	0.0404	0.1268	0.0553	0.1257
	5.80	0.0290	0.1062	0.0403	0.1265	0.0555	0.1302

Acknowledgements

The authors thank Airbus Defence and Space and CNES (French national space agency) for their funding. The authors gratefully acknowledge the French National Agency (ANR) for financial support of the postdoctoral study of M.L. for the development of a massively parallel version of the Black-Box multigrid method to solve linear systems in the frame of the COALA project ANR-15-CE06-0013. The authors gratefully acknowledge the CNRS (Centre National pour la Recherche Scientifique) for funding compensation of teaching hours to S.T. during the year 2017-2018. Some of the computational time was provided by the scientific group CALMIP (project no. P17035), the contribution of which is greatly appreciated.

Supplementary movie

Supplementary movie is available at https://doi.org/10.1017/jfm.2018.389.

Appendix

This appendix contains the evolution of the x-coordinate (figure 19) and y-coordinate (figure 20) of the centre of mass of the fluids throughout time, for different values of Bo i and Bo c , and for a filling ratio of 50 %.

In figures 19 and 20, we observe distinctly that for the lower values of Bo c and Bo i , the coordinates of the centre of mass and their evolution through time are identical (see § 5.3). By increasing the centrifugal Bond number, the differences between the trajectories become more important. In figure 19(e-g), the oscillations of the centre of mass around the x-axis are clearly visible for the higher values of Bo i . As in § 5.4, we observe that the time period of one oscillation decreases with Bo c and does not seem to vary with Bo i for its three higher values. Moreover, the oscillations disappear sooner in time with a higher value of Bo c . This is shown in § 5.4 with the computation of the damping ratio. Finally, we notice that the final position of the centre of mass does not depend on Bo i . In the plots of figure 19, the final x-coordinate approaches