Towards New Controller Design Oriented Models of

Propellant Sloshing in Observation Spacecraft

Anthony Bourdelle, Laurent Burlion, Jean-Marc Biannic, Hélène Evain,
Sabine Moreno, Christelle Pittet, Alexis Dalmon, Sébastien Tanguy, Tarek
Ahmed-Ali

To cite this version:

Anthony Bourdelle, Laurent Burlion, Jean-Marc Biannic, Hélène Evain, Sabine Moreno, et al.. To-
wards New Controller Design Oriented Models of Propellant Sloshing in Observation Spacecraft. AIAA
Scitech 2019 Forum, Nov 2019, SAN DIEGO, United States. ฀10.2514/6.2019-0115฀. ฀hal-02502690฀

HAL Id: hal-02502690

https://hal.archives-ouvertes.fr/hal-02502690
Submitted on 9 Mar 2020

HAL is a multi-disciplinary open access L’archive ouverte pluridisciplinaire HAL, est

archive for the deposit and dissemination of sci- destinée au dépôt et à la diffusion de documents
entific research documents, whether they are pub- scientifiques de niveau recherche, publiés ou non,
lished or not. The documents may come from émanant des établissements d’enseignement et de
teaching and research institutions in France or recherche français ou étrangers, des laboratoires
abroad, or from public or private research centers. publics ou privés.
 Towards New Controller Design Oriented Models of
Propellant Sloshing in Observation Spacecraft

A. Bourdelle ∗ , L. Burlion † , and J.-M. Biannic ‡

ONERA - the French Aerospace Lab, Toulouse, France

H. Evain § , S. Moreno ¶ and C. Pittet ‖

CNES - the French Space Agency, Toulouse, France

A. Dalmon ∗∗ and S. Tanguy ††

Institut de Mécanique des Fluides de Toulouse, Toulouse, France

T. Ahmed-Ali ‡‡
University of Normandie and ENSICAEN, Caen, France

The design of challenging and ambitious space missions entails a tightening of spacecraft
pointing constraints. Among the many perturbations that have to be addressed, the sloshing
of the on-board propellant is a complex issue. Recent developments in Computational Fluid
Dynamics, supported by in-situ experiments like Fluidics or Sloshsat-FLEVO, open up ways
for the characterization of liquids behavior inside tanks in microgravity. This knowledge can
be applied in the context of spacecraft modeling and control. For this purpose, we present here
a new approach to model the disruptive sloshing dynamics affecting spacecraft during attitude
maneuvers. With the aim of mitigating slosh effects, this model will be used to design robust
attitude controllers.

I. Introduction
luid dynamics defines sloshing as the movement of the free surface of a fluid inside tanks or containers. This
F complex nonlinear dynamics leads to disruptive effects which have to be addressed when a mass of fluid can move
aboard vehicles. Sloshing then appears in tanker ships or trucks, but also in spacecraft. Indeed, orbital maneuvers
(station-keeping, relocation and de-orbiting) require thrusters that consume liquid propellant. Since it is a strong
constraint for the spacecraft lifespan, the mass of the on-board liquid propellant is a non-negligible part of the total
spacecraft mass, making sloshing effects even more significant. Furthermore, sloshing is a low frequency phenomenon,
it can compromise the system stability and it makes controller design more complex.

As a result of the coupled fluid-structure dynamics, the spacecraft undergoes disturbing efforts that alter both its
attitude and orbit. In rare cases, propellant slosh dynamics can lead to serious consequences, from mission delay (NEAR,
1998 [1]) to mission loss (ATS-5, 1969 or Falcon 1, 2007 [2]). Several methods exist to deal with sloshing disruptive
effects, passive and active ones. Hence, spacecraft tanks generally have baffles and membranes for slosh mitigation
[3, 4]. In addition, time margins can be applied to let propellant come back at rest after aggressive maneuvers, or the
desired angular velocity profiles can be smoothed to lessen the fluid excitation. Among active methods, frequency-based
controllers are used to filter slosh dynamics [5].

∗ PhD Student, ONERA/DTIS, [email protected]

† Research Engineer, ONERA/DTIS
‡ Research Director, ONERA/DTIS
§ Research Engineer, CNES
¶ Research Engineer, CNES
‖ Research Engineer, CNES
∗∗ PhD student, IMFT/INTERFACE
†† Lecturer-researcher, IMFT/INTERFACE
‡‡ University Professor

1
 However these methods have some drawbacks :

• more complex tanks are also heavier

• time margins and smoothed profiles reduce mission availability
• frequency-based controllers can be rather hard to tune and sensitive to modeling errors [5]

The constraints on spacecraft pointing stability and pointing accuracy are becoming more stringent. Hence, it is
necessary to develop a model of propellant sloshing that is representative enough to allow the design of a suitable
attitude controller, which efficiently deal with sloshing and lighten the previously mentioned drawbacks.

Surface tension arises from the fluid molecular interactions, it creates and maintains the interface between two
media. In microgravity, slosh dynamics becomes more complex as surface tension effects cannot be neglected in front of
gravity [6]. Therefore many researches have been conducted in order to understand sloshing, particularly in space [7, 8].
It is difficult to reproduce microgravity conditions in laboratories (typically with zero-g flights or drop towers) and
analytical descriptions of sloshing effects in such an environment are very complex. In-situ experiments dedicated to the
study of sloshing were led, such as Sloshsat-FLEVO (ESA), Spheres (NASA) and Fluidics (ESA). The flight data have
been used to adjust and validate CFD models. Due to their high computational cost, these CFD models are unsuitable
to be used in a closed loop with the Attitude Control System (ACS), especially for the design of controllers. Thus,
simplified Linear Time Invariant (LTI) models were developed [7, 8], but they do not take into account the dependency
of the fluid response to the spacecraft maneuver.

In this paper we first state the problem, going through various kinds of models, some suitable for validation and
others for design purposes. We then present our new modeling approach for propellant slosh dynamics in spacecraft, in
which the fluid behavior is modeled by a nonlinear second order system with varying frequency and damping ratio.
Such a model, exploiting similarities with any standard representation of LTI flexible modes [9], is expected to ease the
design of attitude control systems with sloshing. The last part presents our identification procedure, applied on data
computed with the CFD code DIVA (IMFT), and the subsequent results.

II. Problem statement

Attitude Control Systems have to ensure a pointing requirement after an attitude maneuver despite perturbations
such as sloshing. To do so, controller-design and validation models are used during spacecraft ACS development.

The most representative and accurate models are those used in CFD. These models describe the coupled fluid-
structure dynamics in microgravity with the Navier-Stokes equations extended by terms that depend on the propellant
surface tension and the spacecraft inertial forces (linked to the angular velocity and acceleration). Recent developments
led to the accurate computation of these equations. For instance, the CFD codes DIVA (Institut de Mécanique
des Fluides de Toulouse) [10, 11] and COMFLO (University of Groningen) [12] are able to solve fluid-structure
problems in microgravity. In particular, they can compute the forces and torques applied by the fluid on the tank. An
interesting point is that both DIVA and COMFLO solvers were improved and validated, using respectively Fluidics
and Sloshat-FLEVO experiments. However CFD solvers have a high computational cost, which makes difficult the
closed-loop implementation with ACS (cf. figure 1, in which ΓF is the fluid torque, ΓC is the actuators torque and ΓP is
the torque from other perturbations, such as solar radiation pressure).

Fig. 1 System block-diagram with a closed-loop CFD model

2
 A first solution is to use CFD in open loop as shown in figure 2. As there are no longer any interactions between the
satellite dynamics and the solver, the implementation and the use of this model are easier. Yet it is only valid when the
spacecraft angular velocity Ω is close to the desired one Ωd . It requires the execution of many simulations to cover er-
rors, still at high computational cost. Therefore, this model can be suitable for controllers validation but not for their design.

Fig. 2 System block-diagram with an open-loop CFD model

In order to design controllers, simplified mechanical-inspired models were developed. In these models, the fluid
behavior is described by a mechanical system such as spring-mass, pendulum [13], free-mass or mass constrained on a
surface [14, 15]. Equivalent Mechanical Models can be used to describe a complete satellite system including solar
panels [16]. They are either based on a linearized fluid theory or on (semi-) empirical laws [17, 18]. These models,
along with the passives methods mentioned in section I, have been successfully used for decades, for launchers and
satellites [19]. Such models can be described with the standard representation of LTI flexible modes (cf. eq. 1 and 2).
ΓF = LT ηÜ (1)
ηÜ + CS ηÛ + KS η = −L Ω Û (2)
In which η ∈ RnS are the nS flexible modes, L ∈ R3×nS is the matrix of the modal participation of the flexible
modes, CS = diag(2ξSi ωSi ) and KS = diag(ωS2 i ), where ωSi and ξSi are respectively the natural frequency and the
damping ratio of the i-th mode.

Equivalent Mechanical Models can be used in close-loop (cf. figure 3), providing a true interaction with the
rigid spacecraft. Yet in these models the mechanical parameters (generally constant or possibly time-varying while
considering propellant consumption [20]) do not evolve along with the spacecraft state.

Fig. 3 System block-diagram with a closed-loop LTI model

By using CFD solvers it is possible to characterize sloshing disruptive efforts depending on the spacecraft current
maneuver. Such a study has been carried out with DIVA at IMFT [10], with a focus on the dependency of the liquid
response to the centrifugal (eq. 3) and the impulse (eq. 4) Bond numbers.

ρΩ2 LD02 Û
ρΩLD 2
cent imp 0
BO = (3) BO = (4)
σ σ
Where Ω is the spacecraft angular speed around the considered axis of rotation, L is the distance between the tank
and the rotation axis (lever arm), ρ is the fluid density, σ is the fluid surface tension and D0 is the pressurization gas
bubble initial diameter (linked to the filling ratio).

3
 (a) Tank enforced angular velocity

(b) Torque ΓZ along the Z-axis

Fig. 4 Example of DIVA results for a 4.72 × 10−2 rad/s steady-state velocity

The figure 4 shows the time evolution of the sloshing torque along the rotation axis, for different values of the
enforced angular acceleration. The torque discontinuities that appear at the start of each response are due to the square
shape of the acceleration . As the system reaches a constant angular velocity, the torque goes to zero. The only remaining
effort applied by the fluid on the tank comes from the centrifugal effect that pushes the fluid against the tank wall,
which does not create a torque along the Z-axis. The data allow us to suppose that the torque could be approximated by
a damped quasiperiodic function, with a frequency and a damping ratio that are functions of the spacecraft angular
velocity and acceleration.

4
 Usually, in Attitude Control Systems [21–23], the guidance loop provides an angular velocity profile in triangular
shape or truncated triangular shape (cf. figure 5) for each satellite frame axis, similar to the profiles considered in the
IMFT study. These profiles are computed to realize the attitude maneuver in near-minimal time given the spacecraft
maximum angular velocity and acceleration constraints. The desired attitude profile is then obtained by integration of
the speed guidance profile. By derivation, a bang-bang or bang-stop-bang torque input to be realized by the actuators is
readily deduced. Consequently we will focus on these profiles.

(a) Bang-bang profile (b) Bang-stop-bang profile

Fig. 5 Angular velocity (dashed) and command torque (solid) profiles

III. Modeling of propellant Slosh Dynamics as a Nonlinear Second Order System

By studying the Navier-Stokes equations in microgravity, we can identify the parameters affecting sloshing :

• tank filling ratio

• gravity vector orientation w.r.t. the spacecraft, linked to the attitude θ
• liquid properties, e.g. density, viscosity, surface tension
• tank geometry and position inside the spacecraft
• angular speed Ω and acceleration Ω Û

The system on which we apply our modeling method is the one that have been used in the IMFT study. We consider
that thrusters are not used during attitude maneuvers, so the filling ratio is constant. This hypothesis is relevant since
propellant is saved for orbit control, then reaction wheels and control moment gyros are preferred for both attitude
control and attitude guidance maneuver. We also suppose that gravity effects can be neglected (microgravity) and that
the propellant properties do not change. The tank position inside the spacecraft is fixed, and we consider a rigid tank.
Hence the remaining variables are the angular speed and acceleration.

Let split the liquid-filled satellite in two subsystems: the dry satellite and the liquid propellant. By doing so, in order
to study the satellite dynamics, we have to consider only the perturbations exerted by the liquid on the spacecraft, instead
of the propellant behavior. From the preceding discussion, we propose to model the sloshing torque as the output of a
second order system with a frequency and a damping ratio as functions of the satellite angular velocity and acceleration,
thus we obtain the system shown in figure 6. Doing so, equations 1 and 2 become :

ΓF = f (η, t) or f (η,
Ü t) (5)
Û t)ηÛ + KS (Ω, Ω,
ηÜ + CS (Ω, Ω, Û t)η = −g(Ω, Ω, Û t) (6)

Once the shape (e.g. linear or polynomial) of the functions f , g, CS and KS is chosen, the necessary parameters
can be identified using CFD results. Note that we introduced a dependency to the time t, which has appeared to be
necessary during the identification procedure development.

The fluid block can be easily added to an already existing satellite model as a feedback, making it plug and play.
Even though the first mode energy contribution is the highest, the representation we use allows to add as many sloshing
modes as necessary. Thus the model complexity, e.g. the state-vector size, can be adjusted at will.

5
 Fig. 6 System block-diagram of the proposed closed-loop model

This approach can be considered as a generalization or an abstraction of Equivalent Mechanical Models, but it
does not have to deal with the choice of equivalent mechanical models parameters (e.g. masses, stiffness coefficient,
pendulum length and moment of inertia) which are generally taken constant and based on empirical or semi-analytical
expression often only valid for axisymmetric tanks.

We will now describe our identification method and illustrate the precision of our identified model with respect to
CFD solver results on a given example.

IV. Identification and Numerical Results

The CFD tool considered in this study to model propellant sloshing is the DIVA code (Dynamics of Interface for
Vaporization and Atomization) developed in the Fluid Mechanics Institute of Toulouse. It is based on the level-set
method of Osher and Sethian [24] to describe propellant and gas distributions within the computational domain.
A transport equation is solved at each time step to capture the motion of the fluids throughout time. A second
level-set function describes the static interface between the solid tank wall and the fluids in the tank referential. The
incompressible Navier-Stokes equations for Newtonian two-phase flows are solved using the Ghost Fluid Conserva-
tive viscous Method with an Implicit scheme (GFCMI) of Lepilliez et al. [25]. This method uses the projection
method of Chorin [26] adapted to two-phase flows by Sussman et al. [27]. The boundary conditions at the tank
wall are considered within the computational domain and consist in a subcell Dirichlet boundary condition on the
velocity and a subcell Neumann boundary condition on the pressure. The static contact angle is enforced at the
tank wall by extending the level-set function within the solid domain [25]. The numerical solver has been used to
compute propellant sloshing in spherical tanks in microgravity conditions for rotational manoeuvres at low Bond num-
bers [11] and exhibits a good agreement with data from the Fluidics experiment led in the International Space Station [28].

The data sets used in this study [11] correspond to a spherical tank with a diameter of 0.585 m, the lever arm
is 0.4 m. The fluids that have been considered have properties close to the actual properties of fluids used in space
applications. The propellant density is 1004 kg m−3 and its viscosity is 9.13 × 10−4 kg m−1 s−1 . The pressuriza-
tion gas density is 2.41 kg m−3 and its viscosity is 1.99×10−5 kg m−1 s−1 . The surface tension is equal to 0.03325 N m−1 .

The initial conditions are a tank at rest with the gas as a bubble located at the center of the tank. The static contact
angle is set to zero. The maneuver is a rotation along the Z-axis following a bang-stop profile. We consider only the
Û along this axis. We have 42 data sets, each one corresponds to a bang-stop profile
torque ΓZ and the rotation (θ, Ω, Ω)
with an angular velocity in Ωmax (rad/s) and an acceleration in Ω Û max (rad/s2 ) (cf. eq. 7 and 8).

Ωmax = 4.720 × 10−3, 6.680 × 10−3, 1.055 × 10−2, 1.493 × 10−2, 2.111 × 10−2, 3.338 × 10−2, 4.721 × 10−2 (7)
 

Û max = 4.500 × 10−5, 9.000 × 10−5, 2.230 × 10−4, 4.450 × 10−4, 8.920 × 10−4, 2.229 × 10−3
 
Ω (8)

We are in a black-box case, without any precise idea of the functions f , g, Cs and Ks (cf. eq. 5 and 6). Even if
we describe the system with a nonlinear second order system, it is obvious that stronger nonlinearities occur. Yet, we
can capture some of these nonlinearities by making the identification on several angular velocity and time intervals.
We then get an evolution of the parameters. As we are not able to describe their behavior, we chose to average the
parameters (and indirectly the nonlinearities) to find the best trade-off.

6
 Fig. 7 System in equilibrium configuration (the colored part is the liquid

While examining the simulations results, we noted that the discontinuities in the sloshing torque, arising from the
discontinuities in the acceleration profiles (cf. section II), were not well enough taken into account. Indeed it is a really
strong nonlinear effect. To override this problem, we split the model on both sides of the discontinuity. By doing so we
get two submodels, one before the discontinuity (bang part) and one after the discontinuity (stop part). This kind of
problem can be then addressed with switching control techniques.

In order to chose the functions f , g, Cs and Ks , we initially suppose that g and f are linear. As numerical integration
is very sensitive to initial condition, we prefer to use numerical derivation. Thus the sloshing torque ΓZ is directly the
flexible state η. The resulting Linear Parameter Varying system to identify is then :

ΓZ = η (9)
ηÜ + CS ηÛ + KS η Û
= −AS Ω − BS Ω (10)

We consider ranges of interest for :

• the angular acceleration, Ω Û r ng = [−ΩÛ N , . . . , −Ω

Û 1, 0, Ω
Û 1, . . . , Ω
Û N ], N ∈ N
• the angular velocity, Ωr ng = [−Ω M , . . . , −Ω1, 0, Ω1, . . . , Ω M ], M ∈ N
• the time, Tr ng = [0, T1, . . . , TP ], P ∈ N

By using the data from each profiles, we then construct tensors C, K, A and B. For the model before the discontinuity,
the (i, j)-th element of each tensor corresponds to the value of the associated function evaluated at (Ω Û i, Ω j ). For the
model after the discontinuity, the (i, j, p)-th element of each tensor corresponds to the value of the associated function
evaluated at ( Ω Û 0, j , Ω j , Tp ). Ω
Û 0, j ∈ Ω
Û r ng is the acceleration that brought the tank to the velocity Ω j . It is not possible
to have Ω j > 0 and Ω Û 0, j < 0, and Ω j < 0 and Ω Û 0, j > 0, hence the use of the absolute value of the acceleration. The
procedure is illustrated by figures 8 and 9.

(a) Before the discontinuity (time invariant) (b) After the discontinuity (time varying) at time Tp ∈ Tr ng

Fig. 8 Illustration of matrices construction and the used subscripts

7
 Fig. 9 Splitting of bang-stop profiles

A. Identification in the case Ω > 0 and Ω Û > 0 → Model before the discontinuity
This case allows us to fill the (L, R) blocks of each matrix (cf. figure 8) of the model before the discontinuity.
It corresponds to the part of the profiles with a constant acceleration and a velocity ramp. We consider the profiles
with the same angular acceleration Ω Ûj ∈Ω Û max . The number of velocity intervals N is chosen small enough to catch
nonlinearities and to avoid over-averaging, and large enough to still have enough data to perform the identification. For
this case, it is sufficient to use time invariant parameters.

On each velocity interval [Ωi, Ωi+1 ], i ∈ [1, . . . , M], we solve a linear identification problem Y = PX. The
measurement vector is given by Y = ΓÜ Z (t), which is computed using central finite differences. The observation matrix
Û
is X = ΓÛ Z (t), ΓZ (t), Ω(t), Ω(t)

. The parameter vector is P = [C, K, A, B]. In order to ensure that the values of C and
K make sense from a physics point of view, we bound them using maximal and minimal values of sloshing frequency
and damping ratio. The sloshing frequency of each data set is computed using spectral analysis with a Fast Fourier
Transform. We also add constraints to enforce both initial and final conditions on each interval. The resulting constrained
optimization problem is then solved with a linear least square method by using the MATLAB lsqlin function. As it is a
convex problem, we get a global result.

Each interval is common between different profiles with the same maximum angular acceleration. Hence the result
for a given velocity interval is taken as the average of the values obtained for the same interval on different profiles.
Û j . To cover
Finally we get the value of the parameters for the average velocity on the interval Ωi/2 and the acceleration Ω
the whole ranges Ωr ng and Ω Û r ng , we interpolate using the MATLAB function interp2 with a cubic-spline method. Note
that even though the model before the discontinuity is time invariant, the parameters are not constant and change along
the angular velocity ramp.

B. Identification in the case Ω > 0 and Ω Û = 0 → Model after the discontinuity

This case allows us to fill the (L, R) blocks of each tensor (cf. figure 8) of the model after the discontinuity. It
corresponds to the part of the profile that with a zero acceleration and a constant velocity. We consider profiles with the
same steady-state velocity Ωi ∈ Ωmax and various accelerations Ω Û 0, j ∈ k Ω
Û r ng k. The methodology is almost the same
as before. We define time intervals of a given length to catch nonlinearities without over-averaging, then we solve the
constrained linear least square problem. We get the value of the parameters for the velocity Ωi and the acceleration Ω Û 0, j
Û
at any given time in Tr ng . We then interpolate to cover the ranges Ωr ng and Ωr ng .

C. Identification in the case Ω = 0 and Ω Û > 0 → Model before the discontinuity

This case allows us to fill the (L, M) blocks of each matrix (cf. figure 8) of the model before the discontinuity. Good
results are obtained by extrapolating or duplicating values from (L, R) for the coefficients C(L, M) , K(L, M) and A(L, M) .
B(L, M) is then computed by using the initial conditions, which are the same for each profile with the same angular
acceleration. For this case it is also sufficient to use time invariant parameters.

8
D. Identification in other cases
Using the property that the system is symmetrical, we have (for instance for C):

Û i, Ω j , Tp ) = C(−Ω
• C(Ω Û i, −Ω j , Tp )
Û
• C(Ωi, 0) = C(−Ωi, 0)Û
• C(0, Ωi ) = C(0, −Ωi )

We can then fill the blocks (U, L), (U, M) and (L, L). If we suppose that the system is reversible, we flip the data
along the time axis and use the identification procedure to fill the blocks (U, R) and (L, L). As we are able to fill the
matrices blocks for which we do not have data, by using the system properties, our method requires less CFD simulations
and thus reduce the duration of identification and validation steps.

In order to go from the matrices or tensors to the functions we can apply curve-fitting techniques, for instance a least
square method applied on a polynomial function of the angular speed, acceleration and time. For validation purposes
we can skip this step and use the SIMULINK Look Up Table block.

The whole identification code execution time is about two minute of time. It is fast enough to allow us to tune
quickly the identification parameters (intervals lengths) to find the best fitting.

E. Results
For each data set, we compute the error ǫ between the torque given by our model ΓZm and the data ΓZ . In order to
quantify the quality of our estimation, we get the relative error ||ǫ ||2 /||ΓZ ||2 using L 2 norm. The tables 1 and 2 show
the results for both parts of the model, before and after the discontinuities.

Ωmax,1 Ωmax,2 Ωmax,3 Ωmax,4 Ωmax,5 Ωmax,6 Ωmax,7

Û max,1
Ω 0.5384 7.1728 5.0418 3.3715 3.1191 3.0924 3.9071
Û max,2
Ω 0.3856 0.3395 1.8239 1.6297 1.0071 0.8505 1.2408
Û max,3
Ω 0.0398 0.2123 0.0783 0.3534 0.1537 0.2346 0.4516
Û max,4
Ω 0.2333 0.5601 0.7921 0.7847 0.6367 0.5613 2.0376
Û max,5
Ω 0.1871 0.9303 1.1248 1.3472 1.5432 1.0736 1.2309
Û max,6
Ω 0.1138 0.5654 0.1559 0.3000 0.2963 0.4635 1.7179
Table 1 Relative error (%) for each profile before the discontinuity

Ωmax,1 Ωmax,2 Ωmax,3 Ωmax,4 Ωmax,5 Ωmax,6 Ωmax,7

Û max,1
Ω 8.7077 3.9964 1.3782 0.8366 7.8829 52.3987 8.4349
Û max,2
Ω 11.8112 3.7038 2.0413 0.7709 2.5303 10.7775 5.7483
Û max,3
Ω 6.0635 2.8097 3.3712 2.1829 3.7826 3.9710 11.3312
Û max,4
Ω 11.7110 2.9681 7.4316 3.3607 4.4995 2.4865 1.9818
Û max,5
Ω 3.7777 4.2716 3.5500 8.1793 2.6146 2.8900 2.3253
Û max,6
Ω 5.5645 8.4743 4.0912 1.8475 4.0507 2.6784 4.3383
Table 2 Relative error (%) for each profile after the discontinuity

Results may change along with SIMULINK model configuration parameters, in particular the relative tolerance and
the maximum time step. There is a trade-off between precision and computation speed.

9
 For the part before the discontinuity the results are excellent, the mean relative error is 1.23%, with a maximum at
7.17% and a minimum at 0.04%. For the part after the discontinuity, the results are almost all under 10%. The mean
relative error is 5.89%, with a maximum at 52.4% and a minimum at 0.77%. The second highest relative error is 11.8%.
The worst case is presented in the figure 11a, we can see that it is still usable while considering robust control. Still, we
are currently working on some improvements, mostly about tuning our algorithm parameters. Several other results are
presented in figures 10 and 11 .

(a) Case (1, 2) - 7.17% (b) Case (1, 3) - 5.04%

(c) Case (2, 7) - 1.24% (d) Case (3, 1) - 0.04%

Fig. 10 Example of the performances of the model "before the discontinuity"

This modeling approach for propellant sloshing in spacecraft is rather simple, yet it represents well the complex
physics of the problem. As we introduce more physical knowledge into the model, the controller will be more adapted
and will respond specifically to the sloshing dynamics. The main drawback of this method is the calculation time
required to execute a CFD solver, but the computed data can also be used for validation purposes (cf. figure 2). However,
as we split the system, the identification procedure has to be done one time for each kind of spacecraft tank. It can then
be conducted in a spacecraft development life cycle. This method is particularly appropriate for satellite bus product
lines, the identification and the controller design steps having to be done only once for several almost identical satellites.
The model, along with the identification method, can be applied to various kind of tanks, even those for which we cannot
use a symmetry hypothesis. In this case the model simply requires more CFD simulations. Multi-axis attitude maneuver
cases are under study.

V. Conclusion and perspectives

In this paper, a new control design oriented modeling approach of sloshing dynamics in spacecraft has been
presented. The fluid dynamics is efficiently approximated by a generalized second order system with varying frequency
and damping ratio which are written as functions of the spacecraft angular speed and acceleration. Using results from
the CFD solver DIVA of IMFT, we have been able to successfully identify our model parameters. By addressing the
dependency of the fluid response on the spacecraft inertial forces, the proposed model includes more knowledge, which
makes it more accurate. By doing so the highly nonlinear nature of sloshing can be reproduced by our model.

10
 (a) Case (1, 6) - 52.39% (b) Case (2, 1) - 11.81%

(c) Case (2, 4) - 0.77% (d) Case (3, 1) - 6.06%

Fig. 11 Example of the performances of the model "after the discontinuity"

Yet, our approach is simple enough to design robust attitude controllers. The developed procedure also allows us to
characterize identification uncertainties. These uncertainties, along with the ones on both the measurement and the
modeling of the dry satellite, will be addressed using robust control techniques. The resulting controller will be able to
deal with the errors and to ensure the spacecraft stability in spite of it.

VI. Acknowledgement
This study was co-funded by ONERA and CNES. We are grateful to Alexis Dalmon and Sebastien Tanguy (IMFT)
for their collaboration. We also thank Alexandre Janot (ONERA/DTIS) and Jean-Sébastien Schotté (ONERA/DADS)
for their advice, respectively in identification and in fluid dynamics.

References
[1] Hoffman, E. J., Ebert, W., Femiano, M., Freeman, H., Gay, C., Jones, C., Luers, P., and Palmer, J., “The near rendezvous burn
anomaly of december 1998,” Applied Physics Laboratory, Johns Hopkins University, Tech. Rep, 1999.

[2] SpaceX, “Falcon Demo Flight 2. Flight Review Update,” Tech. rep., July 2007.

[3] Dodge, F. T., “Engineering study of flexible baffles for slosh suppression (NASA CR-1880),” Tech. rep., 1971.

[4] Tam, W., Dommer, K., Wiley, S., Mosher, L., and Persons, D., “Design and manufacture of the MESSENGER propellant tank
assembly,” 38th AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit, 2002, p. 4139.

[5] Charbonnel, C., “H∞ and LMI attitude control design: towards performances and robustness enhancement,” Acta Astronautica,
Vol. 54, No. 5, 2004, pp. 307–314.

[6] Ibrahim, R. A., Liquid sloshing dynamics: theory and applications, Cambridge University Press, 2005.

11
 [7] Abramson, H. N., Bauer, H. F., Brooks, G. W., and Chu, W.-H., “The Dynamic Behavior of Liquids in Moving Containers,
With Applications to Space Vehicle Technology (NASA-SP-106),” Tech. rep., 1966.

[8] Dodge, F. T., et al., The new "Dynamic Behavior of Liquids in Moving Containers", Southwest Research Inst. San Antonio, TX,
2000.

[9] Mazzini, L., Flexible Spacecraft Dynamics, Control and Guidance, Springer, 2015.

[10] Lepilliez, M., “Simulation numérique des ballotements d’ergols dans les réservoirs de satellites en microgravité et à faible
nombre de Bond,” Ph.D. thesis, Université Paul Sabatier-Toulouse III, 2015.

[11] Dalmon, A., Lepilliez, M., Tanguy, S., Pedrono, A., Busset, B., Bavestrello, H., and Mignot, J., “Direct numerical simulation of
a bubble motion in a spherical tank under external forces and microgravity conditions,” Journal of Fluid Mechanics, Vol. 849,
2018, pp. 467–497.

[12] Veldman, A. E., Gerrits, J., Luppes, R., Helder, J. A., and Vreeburg, J., “The numerical simulation of liquid sloshing on board
spacecraft,” Journal of Computational Physics, Vol. 224, No. 1, 2007, pp. 82–99.

[13] Reyhanoglu, M., and Hervas, J. R., “Nonlinear control of a spacecraft with multiple fuel slosh modes,” Decision and Control
and European Control Conference (CDC-ECC), 2011 50th IEEE Conference on, IEEE, 2011, pp. 6192–6197.

[14] Berry, R. L., and Tegart, J. R., “Experimental study of transient liquid motion in orbiting spacecraft (NASA-CR-144003),” Tech.
rep., 1975.

[15] Vreeburg, J., and Chato, D., “Models for liquid impact onboard Sloshsat FLEVO,” Space 2000 Conference and Exposition,
2000, p. 5152.

[16] Gasbarri, P., Sabatini, M., and Pisculli, A., “Dynamic modelling and stability parametric analysis of a flexible spacecraft with
fuel slosh,” Acta Astronautica, Vol. 127, 2016, pp. 141–159.

[17] El-Kamali, M., “Ballottement des Liquides avec Tension Superficielle : Etudes Statique et Dynamique,” Ph.D. thesis,
Conservatoire National des Arts et Métiers, 2010.

[18] Utsumi, M., “A mechanical model for low-gravity sloshing in an axisymmetric tank,” Journal of applied mechanics, Vol. 71,
No. 5, 2004, pp. 724–730.

[19] Enright, P. J., and Wong, E. C. ., “Propellant Slosh Models for the Cassini Spacecraft,” Tech. rep., Jet Propulsion Laboratory,
California Institute of Technology, 1994.

[20] Hervas, J. R., and Reyhanoglu, M., “Control of a spacecraft with time-varying propellant slosh parameters,” Control, Automation
and Systems (ICCAS), 2012 12th International Conference on, IEEE, 2012, pp. 1621–1626.

[21] Somov, Y., Butyrin, S., Somov, S., and Hajiyev, C., “Attitude Guidance, Navigation and Control of Land-survey Mini-satellites,”
IFAC-PapersOnLine, Vol. 48, No. 9, 2015, pp. 222–227.

[22] Liu, X., Xin, X., Li, Z., Chen, Z., and Sheng, Y., “Near minimum-time feedback attitude control with multiple saturation
constraints for agile satellites,” Chinese Journal of Aeronautics, Vol. 29, No. 3, 2016, pp. 722–737.

[23] Sidi, M. J., Spacecraft dynamics and control: a practical engineering approach, Cambridge University Press, 1997.

[24] Osher, S., and Sethian, J. A., “Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi
formulations,” Journal of computational physics, Vol. 79, No. 1, 1988, pp. 12–49.

[25] Lepilliez, M., Popescu, E. R., Gibou, F., and Tanguy, S., “On two-phase flow solvers in irregular domains with contact line,”
Journal of Computational Physics, Vol. 321, 2016, pp. 1217–1251.

[26] Chorin, A. J., “A numerical method for solving incompressible viscous flow problems,” Journal of computational physics,
Vol. 2, No. 1, 1967, pp. 12–26.

[27] Sussman, M., Smith, K. M., Hussaini, M. Y., Ohta, M., and Zhi-Wei, R., “A sharp interface method for incompressible
two-phase flows,” Journal of computational physics, Vol. 221, No. 2, 2007, pp. 469–505.

[28] Mignot, J., et al., “Fluid dynamics in space experiment,” 68th International Astronautical Congress (IAC), IAC, 2017.

12