Numerical Study of the Magnetic Damping Effect on the Sloshing of Liquid Oxygen in a Propellant Tank

: Nowadays, the use of bafﬂe plates is anticipated to be one of potential devices used to dampen the sloshing of propellant in rocket tanks. However, some of previous studies reported that the use of a bafﬂe plate may cause larger pressure ﬂuctuations in the tank. In this study, aiming at damping the sloshing without a bafﬂe plate, we paid attention to the characteristic that liquid oxygen is paramagnetic and numerically investigated damping effect of a magnetic ﬁeld when liquid oxygen sloshing occurs. An incompressible gas–liquid two-phase ﬂow of gaseous oxygen and liquid oxygen was assumed in a spherical spacecraft tank with a diameter of 1 m in a non-gravitational ﬁeld, and a triangular impact force was assumed to be imposed as the excitation force. In addition, an electric circular coil was placed outside the spherical tank to generate a static magnetic ﬁeld. For the sake of simplicity, the effect of heat was not taken into consideration. As a result of computation, the sloshing was damped to a certain extent when the magnetic ﬂux density at the coil center was 1.0 T, and a sufﬁcient damping effect was obtained by setting it to 3.0 T. In fact, it is anticipated that less than 3.0 T is sufﬁcient if the coil is placed on the tank surface. This may contribute to damping of the movement of the center of gravity of a spacecraft and prevention of mixing of ullage gas into the piping.


Introduction
In a dynamic acceleration environment caused by changes in the thrust and attitude of a spacecraft, the liquid with a free surface in the propellant tank oscillates violently. This phenomenon is called sloshing. Besides, the pressure inside the tank and the center of gravity of the spacecraft fluctuate. Because the propellant occupies most of the weight of the spacecraft, the fluctuation of the center of gravity due to the sloshing adversely affects the attitude control, and there is a risk of misfiring or malfunction [1,2]. Furthermore, when very low temperature liquid and relatively high temperature gas coexist inside of the tank, heat exchange and phase change between gas and liquid are promoted by the sloshing, and the tank pressure drops quickly, resulting in a unstable supply of propellant and in a cavitation at an entry of a turbopump [3][4][5]. There is another risk that the thrust control becomes difficult. One of the potential tools to prevent the sloshing is to use a device called baffle plate, which is installed inside the tank. Research on baffle plates has been carried out since the 1960s [6], and the optimal shape of the baffle plate is still being studied vigorously [7][8][9]. However, according to several previous studies, the baffle plate can suppress the kinematic fluctuation of the liquid surface during the sloshing in a tank, but it cannot necessarily suppress the thermal fluctuation occurring between the gas and the liquid phases. The pressure fluctuation in the tank with a baffle plate may be greater than that without a baffle plate [10]. The other previous study reports that the reason for this is that the liquid phase becomes turbulent due to the inner edge of the baffle plate and the temperature stratification is disrupted. It was found that the temperature exchange between the gas and the liquid phases was promoted by cold liquid exposure to the interface [11]. This may have a serious impact on flight feasibility for a spacecraft currently equipped with baffle plates. If the sloshing can be suppressed without a baffle plate, the thermal stratification will not be destroyed and the dynamic and thermal behaviors with sloshing may be well controlled simultaneously. However, at present there is little accumulation of knowledge about the propellant control without a baffle plate.
In this study, we paid attention to the fact that liquid oxygen, which is one of the propellants widely used in rockets, is a paramagnetic substance. We numerically investigated the damping effect of the magnetic field when liquid oxygen sloshing in a tank occurred.

Computational Model
A computational model in this study is shown in Figure 1. We assume an incompressible gas-liquid two-phase flow of gaseous oxygen and liquid oxygen in a spherical spacecraft tank with a diameter of 1 m under zero gravity. A triangle acceleration force whose maximum value is 0.40 G (0.40 × 9.8 m/s 2 ) is imposed in the x direction, as shown in Figure 2. In this computation, the effects of heat and wettability are not considered. The computation area is a cube containing a spherical tank. It is assumed that the magnetic field is a static magnetic field generated by a circular coil placed outside the spherical tank. Each computational case is performed by varying the position of the coil and the magnetic field strength gauged at the center of the coil. Computational conditions of the coil are shown in Table 1.

Initial Condition
Considering the most dangerous situation on the attitude of the spacecraft, we focused on the case that the initial liquid phase was located in the lower half of the tank under the impact force excited from the lateral direction. The liquid surface is assumed to be flat to make it easier to capture the tendencies of liquid level fluctuations. In fact, it should be noted that in zero gravity the liquid surface would not be flat due to surface tension and wettability. At the initial condition, the excitation force is zero at t = 0 ( Figure 2) and the static magnetic field generated by the coil has been already imposed.

Computational Grid
For the discretization of the computational domain, a regular staggered grid is used in the Cartesian coordinate system. As shown in Figure 3, scalar quantities such as pressure p, density ρ, Level set function φ and VOF function C are defined at the center of the cell (i, j), and vector quantities such as velocity u, v, w and magnetic flux density B x , B y , B z are defined at the interface of the cells (i + 1/2, j), (i, j + 1/2).

Model of Two-Phase Flow
Since the gas-liquid interface is inherently discontinuous, such computation, including the discontinuous interface, tends to be unstable, especially when the density ratio of gas and liquid is very large. Therefore, it is necessary to provide a transitional region with a certain thickness for the discontinuous gas-liquid interface in physical properties, such as density, viscosity and susceptibility. For the estimation of the surface normal force, the continuum surface force (CSF) model [12] was employed in this study. The following shows the equations expressing the physical properties (density ρ, viscosity µ, magnetic susceptibility χ) of gas-liquid two-phase flow with a transitional region.
Here, the subscript L represents the liquid phase, G represents the gas phase and φ is the level set function measured from the interface. H α (φ) is the smoothed Heaviside function expressing the transitional region, and α is the interface thickness. In this computation, α = 1.75∆x, where ∆x is the grid interval.
However, in this study, the density difference between gas and liquid is large, so the calculation of the surface tension was stabilized by using the density-scaled balanced CSF model proposed by Yokoi [13]. H scaled α (φ) represents the density-scaled Heaviside function. Figure 4 shows the Heaviside function and the density-scaled Heaviside function. The density-scaled Heaviside function is deviated to the liquid phase compared to the simple smoothed Heaviside function. By reducing the surface normal force acting on the gas phase, gas velocities near the interface can be suppressed even for the case of high density ratios.

Method to Capture the Dynamic Surface
The coupled level set and volume of fluid (CLSVOF) method possesses an advantage of surface force estimation and volume preservation, since it is a method coupled between the level set method and the Volume Of Fluid (VOF) method [14]. In the present CLSVOF method, the level set function is used for estimation of the surface force, whereas the VOF method is used for calculation of the advection equation. In this study, the phenomena are assumed to be occurred in a zero-gravity environment where the interface shape is greatly deformed and the influence of surface tension is significant. Therefore, the CLSVOF method was used to capture the surface deformation accurately in this study, and the Tangent of Hyperbola for INterface Capturing/Weighed Line Interface Calculation (THINC/WLIC) method [15] was used for advection of the VOF function. The coupling formula between the level set function and the VOF function was defined as shown in Equation (6), where α is the width of the transitional interface of the CSF (continuum surface force) model and φ is the level set function. It should be noted that the VOF function, especially in the case of severe breaking of surface waves, significantly increases or decreases the liquid phase volume due to numerical diffusion and the expansion of the interfacial transitional region. To circumvent this problem, we employed a volume correction method proposed by Sakuraba et al. [16].

Governing Equation
The following shows the dimensional governing equations used in this study.

Continuity equation
Navier-Stokes equation Advection equation (density conservation) Gauss's law for magnetic field ∇ · b = 0 (13) The notation e exc in Equation (8) is the unit vector in the direction of excitation force, and the notation a(t) is acceleration and varies over time according to Figure 2. Equation (9) is the surface normal force in the density-scaled balanced CSF model, and Equation (10) is the magnetizing force [17].
The above dimensional governing equation was made dimensionless using the dimensionless variables and dimensionless numbers defined below, where subscript φ represents physical property depending on the level set function.
• Dimensionless variables • Dimensionless numbers The dimensionless numbers are the Galilei number Ga representing the strength of gravity; the Laplace number Γ represents the strength of surface tension; and the dimensionless number M represents the strength of the magnetic field. Note that i c was calculated from the magnetic flux density at the center of the single-turn coil b c = µ m i c /2r. The following shows the dimensionless governing equations.
Definition of physical properties Navier-Stokes equation Advection equation Biot-Savart law Gauss's law for magnetic field ∇ · B = 0 (23) Figure 5 shows the acceleration as a function of dimensionless time for the model of lateral impact force. At τ = 0.079 × 10 −6 , the acceleration takes its peak, and at τ = 0.158 × 10 −6 , the excitation force reaches zero. The entire computational time is τ = 5.0 × 10 −6 (stated later), and the acceleration can be regarded as an impact force. For discretization, the Euler explicit method was used for the time derivative term, the third-order upwind difference (UTOPIA) [18] was used for the advection term in the Navier-Stokes equation, and the second-order center difference was used for the others.
As a method of obtaining an approximate solution of the line integral in Equation (12), the composite Simpson's rule was used. In fact, magnetic flux density obtained from the Biot-Savart law B Biot (X) does not necessarily satisfy the Gauss's law for magnetic field (Equation (23)), so it is necessary to correct B Biot (X). In order to satisfy Equation (23), a scalar potential φ p was introduced as shown in Equation (24).
Taking the divergence of both sides, Because of the Gauss's law for magnetic field (∇ · B = 0), This is the Poisson equation about φ p , and iterative correction was repeated until this numerical error became below a certain value. The HSMAC method was used to repeat iterative correction in this study.

IB (Immersed Boundary) Method
The tank shape in the square grid in the Cartesian coordinate system was represented using a immersed boundary (IB) method with linear interpolations of velocity, as shown in Figure 6 and the following equations from (27) to (30).
U IP is calculated by linear interpolations with nearest eight IPs (image points). The way to find U IP in the two-dimensional case is described below. Figure 7 shows a schematic of the grid around the IP point. Defining the velocity components in the X direction above and below the IP point are U up and U low , the linear interpolations of Equations (28) and (29) can be used to calculate U up and U low .
with U up , U low from the above equations and a linear interpolation below, U IP can be calculated.
V IP can be calculated likewise.

Boundary Condition
On the surface of the spherical tank, the no-slip condition and the no inflow/outflow condition were applied. In other words, all the velocity components are zero on the boundary, and the pressure gradient is also zero. In addition, the effect of wettability is not taken into account in this computation because the velocity vector near the tank wall is circumvented by the IB method to represent the wall surface. In fact, since the effect of wettability should be significant in zero gravity [19,20], the method for expressing wettability together with the IB method is required as one of the future issues.

Pressure Correction
In general CFD computation, iterative calculation of pressure correction occupies most of the whole computational time. Therefore, reduction of this iterative corrections greatly contributes to reduce the consumed time. Especially in the multiphase flow problem, since a fluid with a large density ratio is handled, the coefficient matrix of the pressure Poisson equation becomes bad condition, and enormous time is required for convergence. As a countermeasure, the pressure correction was performed by using the Projection method compatible with the IB method. In addition, the multigrid method was applied to reduce the number of corrections together with using C++AMP as a programming language and parallel computation on GPU.

Parameters for Computation
Of the parameters for the present computation, the physical property values are shown in Table 2, and the dimensionless parameters are shown in Table 3. Physical properties refer to the values of gaseous oxygen and liquid oxygen at 90 K under atmospheric pressure [21,22].

Validation of the Numerical Method
In order to verify the validity of this computational code, we computed the oscillation of a droplet in zero gravity. Since the validity of the magnetic field calculation method and the gas-liquid two-phase flow calculation method has been verified by previous studies [14,23], in this section we confirm whether the program could be implemented correctly and the number of grids be sufficient.
We assumed that the case where a cubic droplet at the initial state was gently ejected into a weightless field, and the subsequent change was obtained by numerical calculation. The calculation parameters are shown in Table 4. The results are shown in Figure 8. Due to the surface normal force acting at the surface of the droplet, the part with a large curvature entered inside and the interface deformed from the cubic shape to the octahedral shape. The advection of the interfacial shape was successfully performed while keeping its symmetry with respect to the X, Y and Z directions. No results contradicted the physics, and we were convinced that the implementation was completed correctly and the number of grids was sufficient.

Results and Discussion
The computational results of Case 0, 1, 2, 3, 4 are shown in Figures 9-13 respectively. In the visualizations of Cases 1, 2, 3 and 4, the visualization of Case 0 is additionally shown in translucent red for comparison. The translucent sphere represents a spherical tank, and the liquid phase is visualized for the value of the VOF function C ≥ 0.5. The vector represents the magnetic flux density vector.
In Case 0 (no coil), the liquid phase moved along the wall due to the lateral excitation force, and the interface was greatly deformed. Additionally, a small amount of liquid phase remained on the left wall of the tank due to the influence of the IB method. In Case 1 (D/L = 1.0, H = −0.2, b c = 1.0 T), the amount of liquid phase moving along the tank was smaller and the moving speed was slower than that in Case 0, but the interface was not settled. This must be probably because the position of the coil was far from the tank location, and the magnetizing force intensity in the tank was not sufficient. In Case 2 (D/L = 2.0, H = 0.0, b c = 0.1 T), the results were similar to those of Case 0, and no significant difference was observed. It is considered that the magnetizing force exerted in the tank was weak and the effect of damping of the sloshing was insufficient when the magnetic flux density at the center of the coil was 0.1 T. On the other hand, in Case 3 (D/L = 2.0, H = 0.0, b c = 1.0 T), the liquid level fluctuated and droplets appeared from the liquid surface near the center of the tank, but the most liquid phase remained in the lower part of the tank. Under the computational conditions in this study, if the magnetic flux density at the center of the coil was about 1.0 T, the sloshing was suppressed to a certain extent, which might contribute to the damping of the movement of the center of gravity of the spacecraft and the prevention of mixing of ullage gas into piping. Additionally, comparing between Case 1 and Case 3, a tendency was suggested that when the magnetic flux density at the center of the coil was the same, the effect of suppressing the sloshing was higher when the coil location was closer to the tank than when the coil diameter was equal to the tank. In Case 4 (D/L = 2.0, H = 0.0, b c = 3.0 T), the movement of the liquid phase was mostly suppressed, and the liquid surface near the center of the tank only remained wavy. In this computation, considering that the coil was not placed near the tank, at least 3.0 T at the coil center was required for sufficient damping of the sloshing for the tank with a diameter of about 1 m. But it is anticipated that the magnetic field less than 3.0 T is even effective, if the coil actually can be attached on the tank surface. In these computational results, the tendency of damping effect of the sloshing due to differences in the coil position and the magnetic field strength at the coil center was shown, but these computations were performed when the coil was set outside the fluid domain and the coil was not near the tank. This is because if the coil is placed in the fluid domain, the denominator (|R − X| 3 ) diverges to infinity when calculating the magnetic field according to the Biot-Savart law, and the calculation does not converge. Furthermore, the unwanted behavior that a small amount of liquid phase remained on the wall was observed. Since the accuracy of the IB method is restricted by the grid width, the undesirable behavior can be reduced by increasing the number of grids, but the computational time becomes longer. In order to construct a model that can predict more real phenomena, there is room for improvement in the speed of computation and the representation of arbitrary wall shape by the IB method. In addition, we did not introduce wettability in this study, so we need to construct a useful method that can apply wettability to arbitrary wall shape.
On the subject of feasibility of 3 T, since most of superconducting magnets can produce up to 10 T, it is easy to produce 3 T. However, considering the current resistance of the coils, it is highly possible that the coil needs to be superconducting in order to produce 3 T, and therefore it is neither efficient nor economical to install a superconducting system on a spacecraft. Since the sloshing in the present case can be dampened to some extent at 1 T, a superconducting system is not necessarily required, but the details of the tank size, coil shape design, and electric current value should be designed carefully as the future issues.
As other future tasks, it is necessary to demonstrate the magnetic damping effect when the coil comes into contact with the tank, and quantitatively to evaluate not only the movement and deformation of the liquid surface but also the movement of the center of gravity in the tank. In addition, the pressure field and the thermal behavior in the tank have to be predicted by introducing a beneficial model of phase change.

Conclusions
In this paper, damping effect by the magnetizing force for the sloshing of liquid oxygen was discussed by using numerical computations. The following conclusions were obtained.

•
When the magnetic flux density at the center of the coil was the same, the sloshing damping effect tended to be higher when the coil surface was closer to the tank than when the coil diameter was equal to the tank.

•
When the magnetic flux density at the center of the coil was 1.0 T, the sloshing was dampened to a certain extent, and by setting the magnetic flux density at the center of the coil to 3.0 T, a sufficient sloshing damping effect was obtained. Actually, it is anticipated that even if less than 3.0 T, the damping effect is sufficient by placing the coil closer to the tank. This has the potential to contribute to the damping of the center of gravity movement of spacecrafts and the prevention of mixing of ullage gas into the piping. • Future tasks include speeding up the computation, improving the IB method, and introducing wettability to an arbitrary shaped wall. Besides, it is necessary to evaluate quantitatively the shift of the center of gravity in the tank and the thermal behavior when the sloshing is dampened by a magnetic field.