Study of Morphology of Gas–Liquid Interfaces in Tank with Central Column in CSS under Different Gravity Conditions

Abstract

Most space shuttle fuel tanks use a center column to hold the Propellant Management Device (PMD). This paper analyzes the gas–liquid interface state in the tanks with a central column during microgravity experiments conducted in the Chinese Space Station. It launches an extended study to investigate the gas–liquid interface state under different gravity conditions. Using the perturbation method and boundary layer theory, we numerically calculated the morphology of the gas–liquid interface under varying gravity conditions based on the Young–Laplace equation. The results were then compared to those obtained from existing commercial software and were found to be consistent. Based on this, the study develops two types of calculation procedures. The first procedure generates the corresponding shape of the liquid surface by inputting the height of the liquid surface endpoints and the gravity level. The second procedure is based on the targeting method and generates the corresponding liquid surface by inputting the volume of the liquid in the storage tank and the gravity level. The procedures were used to analyze the variation of gas–liquid interface properties under different gravity conditions. This study offers theoretical support for liquid management in aerospace engineering fuel tanks.

1. Introduction

Propellant Management Devices (PMDs) have emerged as a significant area of research in aerospace engineering in recent years. Existing studies have focused on the static equilibrium configuration of liquid surfaces and the interface dynamics of propellants in a microgravity environment, providing crucial scientific references for the design of PMDs in spacecraft. With the advancement of aerospace engineering, spacecraft are tasked with increasingly complex and diverse missions, necessitating more frequent and prolonged attitude adjustments and in-orbit control. During periods of attitude and orbit control, the residual gravity within the spacecraft significantly increases, rendering the propellant in the tank no longer in a microgravity environment. This residual gravity alters the distribution of the propellant within the tank, severely impacting propulsion efficiency and engine safety. In a high-gravity environment, the effectiveness of surface tension is reduced, the gas–liquid interface contracts, and the propellant may not flow along the PMD’s guidance, potentially leading to a failure in propellant egress from the outlet. Addressing the stable operation of tanks in various residual gravity environments presents a new challenge in the design of propellant management devices. Therefore, it is necessary to reconsider the Young–Laplace equation for different gravity levels to determine the impact of gravity on the gas–liquid distribution.

Regarding the design and research of propellant tanks and PMDs in spacecraft, the current focus primarily revolves around the static equilibrium configuration of liquid surfaces, interface dynamics, and their influence on the design of devices. In the realm of static configuration of liquid surface and fluid stability, Jenson et al. conducted hand-held capillary flow contact line experiments on the International Space Station, quantifying the uncertain impact of contact line boundary conditions, which is crucial for understanding and controlling capillary phenomena in multiphase fluid systems aboard spacecraft [1]. Zimmerli et al. explored the theoretical equilibrium of liquid–gas interfaces in propellant tanks using the Surface Evolver algorithm, laying a theoretical foundation for the design of propellant tanks in NASA’s next-generation exploratory spacecraft [2]. Pylypenko et al. proposed a method for calculating the motion parameters of gas–liquid systems in space stage propellant tanks under microgravity, utilizing the latest finite element analysis tools to theoretically support the design of liquid propulsion systems [3]. Chen et al. studied fluid flow in blade-type tanks through numerical simulations and microgravity experiments, finding that the PMD in blade-type tanks effectively achieves liquid–gas interface separation and provides gas-free liquids under microgravity, which significantly impacts the design of propellant management systems whose function is fulfilled efficiently under microgravity [4]. Plaza et al. conducted numerical analyses on free surfaces with thermal capillary flows and vibrations under microgravity, demonstrating that thermal capillary flow and added vibration effectively control the direction and stability of fluid interfaces, which offers new control strategies for fluid management under microgravity [5]. Govindan and Dreyer examined the stability of liquid interfaces during filling processes under microgravity through experiments in the Bremen Drop Tower, discovering that liquid interfaces exhibit stable but non-constant characteristics at different volumetric flow rates, with interface stability related to critical flow rates [6]. Chen investigated the phenomenon of liquid rising between plates through capillary action under microgravity, showing that the rising height of the liquid between plates is related to the angle between plates, the dynamic contact angle between the liquid and the plate walls, and the viscous resistance of the fluid [7].

In terms of advancements and innovations in fluid management devices, J. Hartwig provided a detailed historical review of PMDs used for propellant acquisition under low gravity conditions, covering the design concepts, basic fluid physics, and operating principles of PMDs, which offers valuable references for the future PMD designs [8]. Chung studied the effects of using coating and pulse-flow techniques for cryogenic spray quenching on simulated propellant tank walls under microgravity, showing significant improvements in cooling efficiency and reduced propellant consumption [9]. Minai and Kuzmich explored the optimal design of lateral PMDs through Computational Fluid Dynamics (CFD) methods to enhance the energy characteristics of launch vehicles and reduce propellant residuals [10]. Chato conducted ground tests on the reflux capabilities of Low-temperature Propellant Liquid Acquisition Devices (LADs), aiding in understanding the working principles of LADs under microgravity and their sensitivity to the flow rate and the tank internal pressure [11]. Alipour proposed a design method for PMD systems under zero gravity conditions, laying a new theoretical foundation for the design of PMDs, especially in handling liquid and gas separations and ensuring the propellant supply [12]. Motooka explored the application of porous metals in the gas–liquid equilibrium propulsion systems of small spacecraft, demonstrating through experiments that porous metals can effectively manage liquid propellants [13].

Ensuring continuous supply to propulsion systems will enhance the performance and reliability of small spacecraft. Baeten established a coupled membrane–liquid dynamics model, providing a new design and analysis tool for liquid storage and management systems on spacecraft, with a particular focus on simulating liquid behavior under microgravity and laying a theoretical foundation for the future designs of fluid management systems on spacecraft [14]. Shukla investigated the phase separation performance of Screen Channel Liquid Acquisition Devices (LADs) under microgravity, offering significant theoretical and experimental evidence for designing efficient liquid supply systems that are used in the microgravity environment [15].

Recent research has also begun to focus on the issues of interface sloshing under the influence of minor disturbances. Bourdelle proposed a new design of a model-oriented controller for addressing the issues of propellant sloshing on spacecraft and enhancing spacecraft stability and precision during attitude adjustments and other operations, which is crucial for increasing the success rate of space missions [16]. Liu used the Smoothed Particle Hydrodynamics (SPH) method to simulate liquid sloshing in spacecraft tanks and its impact on spacecraft separation, offering a new perspective for understanding and controlling liquid sloshing [17]. Hou studied the impact of longitudinal excitation on the sloshing behavior of liquid hydrogen in spacecraft tanks under microgravity, providing important reference information for the design of spacecraft thermal management systems, especially in handling the storage and transmission of liquid hydrogen under microgravity [18]. Hu et al. analyzed the fluid sloshing motion in blade-type PMDs with and without anti-slosh baffles under microgravity through numerical simulations and microgravity experiments, finding that anti-slosh baffles significantly improve fluid stability [19]. Liu et al. conducted experimental research on the liquid sloshing behavior in blade-type surface tension tanks used in high-orbit satellites, offering important guidance for the designs of spacecraft structures and control systems [20]. Khoshnood et al. proposed a mechanical model to simulate the fuel slosh dynamics and its impact on the stability and control of spacecraft, aiming to simplify the analysis process and reduce computational workload [21]. Fries et al. explored the modeling of fluid motion in spacecraft propellant tanks, particularly focusing on how fluid sloshing affects the rigid body motion of spacecraft and how it can be controlled through reaction control systems [22]. Dumitrache and Deleanu conducted ANSYS CFX simulations of the behavior of two fluids (seawater and air) inside a spherical cabin, performed fluid–structure interaction analysis, and emphasized the importance of evaluating liquid sloshing effects [23]. Yu proposed a numerical method for the linear sloshing problem of inviscid incompressible liquids considering surface tension effects, particularly relevant in low gravity environments [24]. Leiter et al. discussed the challenges of spacecraft dynamics during orbital correction maneuvers under the effects of liquid sloshing and proposed a robust control methodology to ensure the stability and performance of spacecraft under multi-tank sloshing disturbances [25].

Previous research has primarily focused on the behavior of fluids and vibration mechanisms in microgravity environments. However, during spacecraft maneuvers such as orbit changes, braking, and docking, maintaining high levels of microgravity is not feasible. These maneuvers cause significant changes in liquid distribution and center of mass, which can severely affect the success rate of orbit changes and docking control. Past research does not address these conditions, leaving a scientific gap in understanding liquid distribution and center of mass changes under varying gravity environments. Furthermore, the issue of interface reconstruction in tanks with central columns under these varying gravity conditions has not been comprehensively studied. In this study, an Experimental Unit (EU) is designed, and experiments under microgravity conditions are conducted in the Chinese Space Station (CSS), and the formation of a gas–liquid interface during the filling process is observed. Based on these observations, the distribution of the gas–liquid interface around the central column under different gravity conditions is theoretically derived. To accurately and quickly predict the interface, a numerical calculation program based on the shooting method is developed, which allows the calculation of interface contour through the input volume. The results by the Volume of Fluid (VoF) method are compared with the experiment results in space, and the derived bubble distribution through the VoF method is used to verify the accuracy of theoretical derivations. Based on this, a detailed analysis of the changes in liquid properties within the tank under different gravity conditions is conducted.

2. Filling Experiment for the Tank with a Central Column on CSS

Surface tension tanks are the most extensively utilized satellite propellant tanks worldwide, leveraging surface tension to manage liquid transport and gas–liquid separation, thereby supplying non-gaseous propellant to engines or thrusters. Plate-type tanks, as the predominant variant of surface tension tanks globally, cater to a variety of flow requirements and are adaptable to different microgravity conditions, which are particularly favorable for large satellite platforms with relatively low microgravity environments. They represent the current trajectory of surface tension tank development. To investigate the positioning and static-dynamic characteristics of liquids in microgravity environments and the filling performance of plate-type tanks under such conditions, experiments involving the filling of plate-type tanks with a central column were conducted aboard the Chinese Space Station.

The filling of a tank with a central column (TA) is conducted in space. The designed cross-section of the tank in the space station, depicted in Figure 1, features a spherical shell with a central column and four baffles. The tank’s solvent volume is approximately 108 mL. Its configurations are detailed in

Table 1. Description of TA.

TA	Information
Shape	Spherical
Size	Inner radius 30 mm
Components of PMD	4 guide vanes and 4 small liquid-storage vanes
Type of management	Fully managed
Width of gaps between the guide vanes and the model wall	Gradually reduced from 1.4 mm at the gas port to 0.8 mm at the liquid port.

, where the wall material is transparent plastic, and the PMD is made of titanium alloy. The tank is filled with dyed 10 CST silicone oil, whose properties are provided in

Table 2. Fluid properties (25 °C).

Fluid	μ (kg/(m*s))	ρ (kg/m3)	σ (N/m)	ν (10−6 m2/s)
Air	1.789 × 10−5	1.225	/	1.460 × 10−5
SF10	0.00935	935	0.0201	10

.

During the space experiment, the TA is placed within a comprehensive experimental system. It is initially filled with gas, then the filling liquid is slowly pumped into the TA from the bottom inlet, and the liquid will be expelled through the top outlet. When the filling volume of the liquid is small, it accumulates at the inlet. As the volume of liquid gradually increases due to the small contact angle of the silicone oil, it rapidly ascends along the wall baffle to the outlet and eventually closes off to form a spherical shell-like liquid layer. Figure 2 shows the liquid surfaces at 8 s and 40 s during the filling process.

Through programming, the grayscale values of images from the tank without liquid are subtracted from those of the processed images to obtain pixel information of the liquid interface. Subsequently, the liquid interface is marked in red and denoised to a certain extent. The bubble shapes obtained by the program are shown in Figure 2, where the internal red curve represents the liquid–gas interface, and the external red curve represents the solid-liquid interface. Due to the contact angle of silicone oil on PMMA walls being 0–10°, a noticeable high-curvature bending occurs when the bubble contour approaches the central column.

3. Small Bond Number Theory

In the tank experiments on CSS, we have observed the formation of annular bubbles during the filling process, a phenomenon not previously seen in earlier experiments. Once the acceleration due to gravity surpasses a certain level, the annular bubbles rise to the top, forming an interface that completely separates the gas and liquid phases. However, in aerospace engineering, tanks on orbiting spacecraft frequently encounter various levels of microgravity, necessitating an analysis of their morphological interfaces under different gravity levels. To facilitate a better theoretical analysis, the situation inside the spherical tank is simplified by removing the baffles and storage plates, leaving only the cylindrical central column, as shown in Figure 3. This figure represents a simplified cross-sectional view of the annular bubbles inside the tank and the forces that act on them.

Figure 3 below depicts a cross-sectional view of the liquid distribution inside the tank when the Bond number is less than 1. The silver part represents the tank shell, the blue part represents the liquid in the tank, and the white part represents the bubbles in the tank. The tank has a cylindrical central column with a radius that is constant everywhere and equal to r₁. The z-axis is the vertical axis of symmetry of the tank, and the r-axis is the horizontal line at the tank’s lowest point, with the two axes intersecting at the origin.

Since the liquid surface is axially symmetric about the central axis of the tank, the analysis is conducted only on the right half of the cross-sectional diagram. In the diagram, $\theta$ represents the contact angle of the gas–liquid interface in the tank, $\phi$ represents the angle between the tangent of the profile curve of the liquid surface and the horizontal direction, and $\alpha$ represents the angle between the tank wall and the horizontal line. Point A is where the bubble contour intersects with the central column, with coordinates (r₁, z₁), and Point B is where the bubble contour intersects with the tank wall, with coordinates (r₂, z₂).

In space, the equilibrium free surface of the bubble within the tank satisfies the differential Equation (1),

$$\begin{array}{c}\frac{1}{r}\frac{d}{dr}\frac{rdz/dr}{\sqrt{1+{\left(\frac{dz}{dr}\right)}^2}}-\frac{\rho gz}{\sigma }=\frac{\mathsf{\Delta }P}{\sigma }=A\end{array}$$

(1)

In the equation, σ represents the surface tension of the liquid, ΔP is the pressure difference between the liquid and the gas, g is the local gravitational constant, $\rho$ is the density of the liquid, and A is a constant to be determined.

The term with dz/dr in the equation can be simplified through relation (2).

$$\frac{dz/dr}{\sqrt{1+{(dz/dr)}^2}}=\frac{dz}{\sqrt{d{z}^2+d{r}^2}}=\sin \phi$$

(2)

Then, multiplying Equation (1) by r and integrating with respect to r leads to Equation (3) as follows:

$$r\sin \phi -\frac{\rho g}{\sigma }\int zrdr=\frac{A}{2}{r}^2+B$$

(3)

The gradient of the contour of the liquid surface can be written using Equation (4) as follows:

$$\frac{\mathrm{d}z}{\mathrm{d}r}=\tan \phi =\frac{\sin \phi }{{(1-{\sin }^2\phi )}^{0.5}}$$

(4)

By combining Equations (3) and (4), a differential equation for z as a function of r can be obtained. Solving this differential equation yields the equation for the gas–liquid interface. However, due to the complexity of the equation, it is not directly solvable; thus, perturbation methods are used for the solution. Below is the process of solving the gas–liquid interface equation using perturbation methods.

Perturbation Method Solution Process

Before using the perturbation method, it is necessary to first nondimensionalize Equations (3) and (4) as shown in Equation (5):

$$x=\frac{r}{r_2}\hspace{1em}y=\frac{z}{r_2}$$

(5)

where $A^{\prime }$ and $B^{\prime }$ are the nondimensionalized parameters, $\epsilon =\frac{\rho g}{\sigma }{r_2}^2$, is the nondimensional Bond number.

This leads to Equation (6):

$$\left\{\begin{array}{l}\frac{dx}{dy}=\tan \phi =\frac{\sin \phi }{{(1-{\sin }^2\phi )}^{0.5}}\\ x\sin \phi -\epsilon \int yxdx=\frac{A^{\prime }}{2}{x}^2+B^{\prime }\end{array}\right.$$

(6)

Let $\epsilon$ be the perturbation parameter, and expand $\phi$, $y$, $A^{\prime }$, $B^{\prime }$ respectively to obtain the expanded forms.

$$\left.\begin{array}{c}y=y_0+\epsilon y_1\\ \phi ={\phi }_0+\epsilon {\phi }_1\\ A^{\prime }=A_0^{\prime }+\epsilon A_1^{\prime }\\ B^{\prime }=B_0^{\prime }+\epsilon B_1^{\prime }\end{array}\right\}$$

(7)

Introducing the order of perturbation (Equation (7)) into the system of Equation (6), we obtain the zeroth-order equation

$$\left\{\begin{array}{l}x\sin {\phi }_0=\frac{A_0^{\prime }}{2}{x}^2+B_0^{\prime }\\ \frac{\mathrm{d}{y}_0}{\mathrm{d}x}=\tan {\phi }_0=\frac{\sin {\phi }_0}{{(1-{\sin }^2{\phi }_0)}^{0.5}}\end{array}\right.$$

(8)

with the following boundary conditions:

$$\begin{array}{c}x=\frac{r_1}{r_2},{\phi }_0=\theta -\pi /2\\ x=1,{\phi }_0=\alpha +\theta \hspace{1em}{y}_0=f\left(r_2\right)\end{array}$$

(9)

Based on the calculations in reference [26], the solution is known to be as follows:

$$y_0=[c\frac{r_1}{r_2}F(\phi ,k)+E(\phi ,k)]$$

(10)

$$A_0^{\prime }=2r_2\frac{r_1\sin \left(\theta -\frac{\pi }{2}\right)-r_2\sin (\alpha +\theta )}{r_1^2-r_2^2}$$

(11)

$$B_0^{\prime }=\frac{-r_1^2{r}_2^2\sin (\alpha +\theta )-r_1{r}_2^3\sin (\pi /2-\theta )}{r_1^2-r_2^2}$$

(12)

$${\phi }_0=\arcsin (\frac{A_0^{\prime }}{2}x-\frac{B_0^{\prime }}{x})$$

(13)

where $F(\phi ,k)$ and $E(\phi ,k)$ are elliptic integrals of the first and second kinds, respectively, and

$$c=\frac{r_2\sin \left(\theta -\frac{\pi }{2}\right)-r_1\sin (\alpha +\theta )}{r_1\sin \left(\theta -\frac{\pi }{2}\right)-r_2\sin (\alpha +\theta )}$$

(14)

$$k^2=(r_2^2-c^2{r}_1^2)/r_2^2$$

(15)

Using the zeroth-order solution, a first-order equation can be obtained as follows:

$$\left\{\begin{array}{c}{\phi }_{1}x\cos {\phi }_0+{\int }_{\frac{r_1}{r_2}}^x{y}_{0}tdt=\frac{A_1^{\prime }}{2}{x}^2+B_1^{\prime }\\ \frac{d{y}_1}{dx}=\frac{{\phi }_1}{{\cos }^2{\phi }_0}\end{array}\right.$$

(16)

To deal with the more complex integral terms in the expressions, we incorporate the formula of $y_0$ obtained earlier along with Equation (16) into the integration, yielding the following formula

$$\begin{array}{c}V(x)={\int }_{\frac{r_1}{r_2}}^x{y}_{0}t\mathrm{d}t=\frac{2\pi }{3{r_2}^2}[(2c^2{r}_1^2+3c{r}_1{r}_2+2r_2^2)E(\phi ,k)\\ -c^2{r}_1^{2}F(\phi ,k)+{(r_2^2-a^2{x}^2)}^{\frac{1}{2}}{(a^2{x}^2-c^2{r}_1^2)}^{\frac{1}{2}}]-\frac{2\pi x^2}{r_2}[c{r}_{1}F(\phi ,k)+r_{2}F(\phi ,k)]\end{array}$$

(17)

with the following boundary conditions:

$$\left.\begin{array}{c}x=\frac{r_1}{r_2},{\phi }_1=0\\ x=1,{\phi }_1=0\hspace{1em}{y}_1=0\end{array}\right\}$$

(18)

By incorporating the boundary conditions, the values of the two parameters can be determined as follows:

$$A_1^{\prime }=2\frac{V_0}{r_2(r_1^2-r_2^2)}$$

(19)

$$B_1^{\prime }=\frac{-r_2{V}_0}{(r_1^2-r_2^2)}$$

(20)

where

$$\begin{array}{c}{V}_0={\int }_{\frac{r_1}{r_2}}^1{y}_{0}tdt=\frac{2\pi }{3{r_2}^2}[(2c^2{r}_1^2+3c{r}_1{r}_2+2r_2^2)E(\phi ,k)\\ -c^2{r}_1^{2}F(\phi ,k)+{(r_2^2-r_1^2)}^{\frac{1}{2}}{(r_1^2-c^2{r}_1^2)}^{\frac{1}{2}}]-\frac{2\pi r_1^2}{{r_2}^3}[c{r}_{1}F(\phi ,k)+r_{2}F(\phi ,k)]\end{array}$$

(21)

After substituting the parameters, we obtain the following:

$${\phi }_1(x)=\frac{1}{x\cos {\phi }_0}[\frac{r_2^2-r^2}{r_2^2-r_1^2}{V}_0-V(x)]$$

(22)

Then, by substituting ${\phi }_1$ into Equation (16), we obtain the following expression:

$$y_1={\displaystyle \underset{\frac{r_2}{r_1}}{\overset{x}{\int }}\frac{1}{x{\cos }^3{\phi }_0}[\frac{r_2^2-r^2}{r_2^2-r_1^2}{V}_0-V(x)]}dx$$

(23)

Due to the complexity of integrating $y_1$, the numerical solution can only be obtained through numerical integration algorithms.

4. Theory of Big Bond Number

When the Bond number exceeds 1, the liquid surface tends to flatten so that the influence of surface tension is only significant near the central column and the tank walls. In this situation, the boundary layer theory can be used to solve the equations. It is assumed that there exists a core region covering the majority of the liquid surface, within which $\frac{dx}{dy}=\tan \phi$ is very small, and the boundary layer regions near the tank walls and the central column, where $\phi$ rapidly increases to the specified boundary values. The values of parameters are determined by matching the core and boundary layer solutions in their transition area.

In the core region ($\left|\phi \right|\ll 1$), due to the very small value of $\left|\phi \right|$, the Young–Laplace equation can be written in the following form:

$$\frac{1}{x_c}\frac{d}{d{x}_c}\left\{x_c\frac{d{y}_c}{d{x}_c}\left[1+O\left({\phi }^2\right)\right]\right\}-\epsilon y_c-A’=0$$

(24)

where the subscript c represents the core region.

By introducing boundary condition (9), the solution to this equation is the following:

$$y_c=\frac{A_i}{\epsilon }\left[I_0\left({\epsilon }^{\frac{1}{2}}{x}_c\right)-1\right]\left[1+O\left({\phi }^2\right)\right]+\frac{A_k}{\epsilon }{K}_0\left({\epsilon }^{\frac{1}{2}}{x}_c\right)\left[1+O\left({\phi }^2\right)\right]+h$$

(25)

where $I_0$ and $K_0$ are the Bessel functions of the first and second kind, respectively, and h is to be determined by comparison with the real boundary conditions after solving the equations, with $A_i$ and $A_k$ having the following:

$$A_i+A_k=A’$$

(26)

In the boundary layer region, referring to the derivation in reference [27], set

$$\beta ={\epsilon }^{-\frac{1}{2}}$$

(27)

The boundary layer at the tank wall has a thickness of

$$x_w=\beta (1-r)$$

(28)

Now, using this as a variable, we reconstruct the Young–Laplace equation in the vicinity of the tank wall as follows:

$$\begin{array}{c}\frac{1}{x_w}\frac{d}{d{x}_w}\frac{x_{w}d{y}_w/d{x}_w}{\sqrt{1+{\left(\frac{d{y}_w}{d{x}_w}\right)}^2}}-\epsilon y_w=A_w\end{array}$$

(29)

where the subscript W represents the boundary layer region near the tank wall.

By introducing $\phi ={\tan }^{-1}(d{y}_w/d{x}_w)$, Equation (29) can be transformed into the following:

$$\left.\begin{array}{l}\frac{d{y}_w\left(\phi \right)}{d\phi }=\frac{\sin \phi }{\epsilon y_w-\left(\sin \phi \right)/x_w+A_w}\\ \frac{d{x}_w\left(\phi \right)}{d\phi }=\frac{\cos \phi }{\epsilon y_w-\left(\sin \phi \right)/x_w+A_w}\end{array}\right\},(0<\phi <\frac{1}{2}\pi -\theta )$$

(30)

Next, expanding $y_w$ and $x_w$ using perturbation theory results in

$$\left\{\begin{array}{l}{y}_w=h+\beta y_{w1}+{\beta }^2{y}_{w2}\\ x_w=1-\beta x_{w1}-{\beta }^2{x}_{w2}\end{array}\right.$$

(31)

with corresponding boundary conditions as follows:

$$\left.\begin{array}{l}\phi =\alpha -\theta \hspace{1em},\hspace{1em}{x}_{w1}\left(\phi \right)=x_{w2}\left(\phi \right)=0\\ \phi =0\hspace{1em},\hspace{1em}{y}_{w1}\left(\phi \right)=y_{w2}\left(\phi \right)=0\end{array}\right\}$$

(32)

Next, we solve the first-order equation. During this process, since the value of $A_w$ has a higher order of smallness compared to other quantities, the first-order equation can be simplified to the following form:

$$\left\{\begin{array}{l}d{y}_{w1}/d\phi =\sin \phi /y_{w1}\\ d{x}_{w1}/d\phi =-\cos \phi /y_{w1}\end{array}\right.$$

(33)

Under the boundary conditions, based on the calculations in reference [27], the solution is known to be as follows:

$$y_{w1}=2\sin \frac{1}{2}\phi$$

(34)

$$x_{w1}=\log \left[\tan \frac{1}{4}{\phi }_1/\tan \frac{1}{4}\phi \right]+2\cos \frac{1}{2}{\phi }_1-2\cos \frac{1}{2}\phi$$

(35)

where ${\phi }_1=\alpha -\theta$.

Similarly, the second-order equations are as follows:

$$\left\{\begin{array}{l}\frac{d{y}_{w2}}{d\psi }=-\frac{\sin \phi }{{y_{w1}}^2}{y}_{w2}+\frac{{\sin }^2\phi }{{y_{w2}}^2}\\ \frac{d{x}_{w2}}{d\phi }=\frac{\cos \phi }{{y_{w1}}^2}{y}_{w2}-\frac{\sin \phi \cos \phi }{{y_{w1}}^2}\end{array}\right.$$

(36)

After introducing the boundary conditions, based on the calculations in reference [27], the solution is known to be the following:

$$y_{w2}=\frac{2}{3}\frac{1-{\cos }^3\frac{1}{2}\phi }{\sin \frac{1}{2}\phi }$$

(37)

$$\begin{array}{c}{x}_{w2}=-\frac{2}{3}{\sin }^2\frac{1}{2}{\phi }_1-\frac{1}{6}{(1+\cos \frac{1}{2}{\phi }_1)}^{-1}+\frac{1}{2}\log \frac{\tan \frac{1}{4}{\phi }_1}{\tan \frac{1}{4}\phi }\\ +\frac{2}{3}{\sin }^2\frac{1}{2}\phi +\frac{1}{6}{(1+\cos \frac{1}{2}\phi )}^{-1}\end{array}$$

(38)

Next, the boundary layer solution and the core solution are matched to determine the values of parameters. For the boundary layer solution, as $\phi \to 0$, $y_w\to 0$ and $x_w\gg 1$, so after simplification by substituting $x_w$ with $\phi$, we can obtain the following:

$$\begin{array}{c}{y}_w=4\beta \{1+\frac{1}{2}[1-x_w(\phi )]+\beta [\frac{19}{12}-\cos \frac{1}{2}{\phi }_1-\frac{3}{3}{\sin }^2\frac{1}{2}{\phi }_1-\frac{1}{6}{(1+\cos \frac{1}{2}{\phi }_1)}^{-1}]+o(\beta )\}\\ \times \exp \left[-(1-x_w(\phi \left)\right)/\beta +\log \tan \frac{1}{4}{\phi }_1-2(1-\cos \frac{1}{2}{\phi }_1)\right]\end{array}$$

(39)

For the core solution, let $y_c=y_w$ and $x_c=x_w$, and at this point, $x_c$ approaches 1, and the value of $x_c/\beta$ is large enough; therefore, the expansion can be written as follows:

$$y_c={\left(2\pi \right)}^{-\frac{1}{2}}{A}_i{\beta }^{\frac{5}{2}}\exp (x_c/\beta )[1+\frac{1}{2}(1-x_c)+\frac{1}{8}\beta +o(\beta )]$$

(40)

After matching, $A_i$ can be obtained as follows:

$$A_i=e^{-1/\beta }{\beta }^{-\frac{3}{2}}{A}_{i0}(1+\beta A_{i1}+o(\beta ))$$

(41)

where

$$A_{i0}=4{\left(2\pi \right)}^{\frac{1}{2}}\tan \frac{1}{4}{\phi }_1\exp \{-2[1-\cos \frac{1}{2}{\phi }_1]\}$$

(42)

$${\lambda }_1=\frac{3}{24}-\cos \frac{1}{2}{\psi }_1-\frac{2}{3}{\sin }^2\frac{1}{2}{\psi }_1-\frac{1}{6}{[1+\cos (\frac{1}{2}{\psi }_1)]}^{-1}$$

(43)

For the boundary layer at the central column, we assume its thickness is the following:

$$x_p=\beta (r_1+r)$$

(44)

where the subscript p represents the boundary layer near the central column.

Similar to the situation in the boundary layer at the tank wall, the Young–Laplace equation is written in the same form as Equations (33) and (36). By introducing the boundary conditions

$$\left.\begin{array}{l}\phi =\theta -\frac{1}{2}\pi \hspace{1em},\hspace{1em}{x}_{p1}\left(\phi \right)=x_{p2}\left(\phi \right)=0\\ \phi =0\hspace{1em},\hspace{1em}{y}_{p1}\left(\phi \right)=y_{p2}\left(\phi \right)=0\end{array}\right\}$$

(45)

the equations are similarly solved, and we obtain the following:

$$y_{p1}=2\sin \frac{1}{2}\phi$$

(46)

$$x_{p1}=\log \left[\tan \frac{1}{4}{\phi }_2/\tan \frac{1}{4}\phi \right]+2\cos \frac{1}{2}{\phi }_2-2\cos \frac{1}{2}\phi$$

(47)

$$y_{p2}=\frac{2}{3}\frac{1-{\cos }^3\frac{1}{2}\phi }{\sin \frac{1}{2}\phi }$$

(48)

$$\begin{array}{c}{x}_{p2}=-\frac{2}{3}{\sin }^2\frac{1}{2}{\phi }_2-\frac{1}{6}{(1+\cos \frac{1}{2}{\phi }_2)}^{-1}+\frac{1}{2}\log \frac{\tan \frac{1}{4}{\phi }_2}{\tan \frac{1}{4}\phi }\\ +\frac{2}{3}{\sin }^2\frac{1}{2}\phi +\frac{1}{6}{(1+\cos \frac{1}{2}\phi )}^{-1}\end{array}$$

(49)

where ${\phi }_2=\frac{1}{2}\pi -\theta$.

Now, we match the boundary layer solution with the core solution. For the boundary layer solution at this location, we perform the same operation as for the boundary layer solution at the wall; it can be written as follows:

$$\begin{array}{c}{y}_w=4\beta \{1+\beta [\frac{1}{12}-\cos \frac{1}{2}{\phi }_1-\frac{2}{3}{\sin }^2\frac{1}{2}{\phi }_1-\frac{1}{6}{(1+\cos \frac{1}{2}{\phi }_1)}^{-1}]+o(\beta )\}\\ \times \exp \left[\log \tan \frac{1}{4}{\phi }_1-2(1-\cos \frac{1}{2}{\phi }_1)\right]\left[\frac{1}{\beta }(x_c-\frac{r_1}{r_2})\right]\end{array}$$

(50)

For the core solution, let $y_p=y_c$ and $x_p=x_c$. Near the central column, x approaches $\frac{r_1}{r_2}$. Since we are discussing a spherical tank and since the fill level is essentially between 30 and 70%, it is assumed that $\frac{r_1}{r_2}\to 0$. When $\epsilon$ < 100, the core solution can be expanded and simplified as follows:

$$y_c=A_k\beta \ln {\left(\frac{2\beta r_2}{r_1}\right)}^{-1}\left[\frac{1}{\beta }(x_c-\frac{r_1}{r_2})\right]$$

(51)

From this, we obtain the following:

$$A_k=\ln \left(\frac{2\beta r_2}{r_1}\right)A_{k_0}\left[1+\beta A_{k_1}+o\left(\beta \right)\right]$$

(52)

where

$$A_{k_0}=\tan \frac{1}{4}{\phi }_1\exp (2-2\cos \frac{1}{2}{\phi }_1)$$

(53)

$$A_{k_1}=[\frac{1}{12}-\cos \frac{1}{2}{\psi }_1-\frac{2}{3}{\sin }^2\frac{1}{2}{\psi }_1-\frac{1}{6}{(1+\cos \frac{1}{2}{\psi }_1)}^{-1}]$$

(54)

when $\epsilon$ > 100, it can be approximated that $A_k$ = 0.

After introducing the boundary layer solution into the formulas, we can obtain the following:

$$h=\frac{z_3}{r_2}-y_w(\alpha -\theta )$$

(55)

Regardless of the method used, the formula for calculating the volume V of the liquid in the tank can be obtained as follows:

$$V=2\pi r_2^3{\displaystyle \underset{\frac{r_1}{r_2}}{\overset{1}{\int }}x[y-f(x)]dy}$$

(56)

5. Numerical Simulation

5.1. Numerical Method

To further verify the accuracy of the theory, commercial software FLUENT is used for simulation calculations. The simulation is based on a simplified grid model that includes the central column, tank walls, and the space in between. Due to symmetry, a quarter-sphere mesh model is utilized, with a radius of 50 mm, as shown in Figure 4, where the surface generates a hexahedral mesh on this model, with a total of 124,500 elements. The maximum element size is 2.4 × 10⁻⁴ m. The average mass per unit was recorded at 0.73, with the mean skewness and maximum skewness measured at 0.11 and 0.53, respectively. In addition, considering the drastic changes in the liquid surface near the contact surfaces, a boundary layer method is used to refine the mesh near the surface of the central column and the tank walls.

In the simulation, multiple sets of calculations are carried out by varying the values of contact angle and liquid volume to simulate different operating conditions. To reduce computational effort, the liquid surface is pre-set as a flat plane for the calculations. The initialization model is shown in Figure 5, where the red part represents the gas phase and the blue part represents the liquid phase.

The motion of the fluid is described using the discrete form of the Navier–Stokes (N-S) equations, with its momentum equation shown as follows:

$$\frac{\partial }{\partial t}(\rho \overrightarrow{\nu })+\mathsf{\nabla }\cdot (\rho \overrightarrow{\nu }\overrightarrow{\nu })=-\mathsf{\nabla }p+\mathsf{\nabla }\cdot [\mu (\mathsf{\nabla }\overrightarrow{\nu }+\mathsf{\nabla }{\overrightarrow{\nu }}^T)]+\rho \overrightarrow{g}+\overrightarrow{F}$$

(57)

where $\overrightarrow{g}$ is the local gravitational constant, $\overrightarrow{F}$ represents the body force, and $\rho$ and $\overrightarrow{\nu }$ denote the fluid density and velocity vector respectively. The values of $\rho$ and $\mu$ are not constant across different elements but are determined by the phase field, depending on the properties of all phases passing through the control volume, and their formulae are provided further below.

The mass conservation equation employs a phase-field model, described by the Volume of Fluid (VOF) model, which depicts the behavior of gas–liquid two-phase flow. This model tracks the interface between phases by solving a continuity equation for the volume fraction of one (or more) phase(s). The following equation describes the evolution of the gas–liquid interface tracked by the VOF method in this model:

$$\frac{1}{{\rho }_q}[\frac{\partial }{\partial t}({\alpha }_q{\rho }_q)+\mathsf{\nabla }\cdot ({\alpha }_q{\rho }_q{\overrightarrow{\nu }}_q)=S_{{\alpha }_q}+\sum _n^{p=1}(\begin{array}{cc}{\dot{m}}_{pq}& -{\dot{m}}_{qp}\end{array})]$$

(58)

This is a discrete equation where ${\dot{m}}_{qp}$ represents the mass transfer from phase q to another phase p and ${\dot{m}}_{pq}$ denotes the mass transfer from phase p to phase q. The subscripts q and p represent the two material phases in the simulation. $S$ stands for source terms, t for time, and $\alpha$ for the fraction of a material in a local cell, with

$${\alpha }_q+{\alpha }_p=1$$

(59)

To solve the flow field, the model solves a single momentum equation across the entire control volume, sharing the resulting velocity field among the phases.

After the phase fractions have been determined in each cell, equations for $\rho$ and $\mu$ in Equation (58) can be solved, with their calculation formulae given below:

$$\rho =\sum {\alpha }_q{\rho }_q$$

(60)

$$\mu =\sum {\alpha }_q{\mu }_q$$

(61)

Substituting these numerical values back into Formula (57) allows for the calculation of the corresponding velocities.

The description of surface tension employs the continuum surface force (CSF) model, which interprets surface tension as a continuous three-dimensional effect on the interface rather than as a boundary condition on the interface. The surface tension effects are represented as body forces on the surface using the divergence theorem and are reflected in the momentum equation, with its formula being the following:

$$F_{\nu ol}={\sigma }_{pq}\frac{\rho {\kappa }_q\nabla {\alpha }_q}{\frac{1}{2}\left({\rho }_q+{\rho }_p\right)}$$

(62)

where ${\kappa }_q$ is the curvature of the fluid phase boundary, ${\sigma }_{pq}$ is the surface tension between the two fluid phases, and $\nabla$ is the gradient operator. The values of $\rho$ are determined by Formula (60).

After the control equations are defined, the properties of the two fluids listed in Table 2 are substituted for fluid phases p and q, respectively, with the boundary conditions on the left and right cross-sections set to be symmetric and the contact angles on the surface of the central column and the tank walls are set to their respective values for solution. Utilizing the aforementioned model, with the time step set to 10⁻⁴ s, the problem is solved with the SIMPLE algorithm until the liquid surface becomes stable, which yields the condition of the liquid surface.

5.2. Comparison with Theoretical Results

Figure 6 displays 2D cross-sectional views of the gas–liquid interface simulation results under different gravity conditions, with a contact angle of 90 degrees. In the figure, the red part represents the gas phase, and the blue part represents the liquid phase.

When the gravity level is within the range of 0.01–10 (m/s²) (Bond number between 0.1 and 100), it can be observed that there are significant fluctuations at the gas–liquid interface. At this time, the effect of surface tension on the liquid distribution is more pronounced. When gravity exceeds 1 (m/s²), the liquid surface tends to be flat, showing significant changes only near the central column and the tank walls. At this point, the influence of gravity is dominant. This change at the interface also verifies the application of the two different calculation methods discussed earlier.

From the numerical calculation results, a gas–liquid isosurface with a phase ratio of 0.3 to 0.7 is extracted. With the calculation program, the intersections of the liquid surface with the surface of the central column and the tank wall, A and B, and the lowest point C of the liquid surface are identified in the figure. We measure the distances of these three points from the tank’s lowest point in the z-direction and the distances of points C and B from the symmetry axis, which correspond to the coordinates of points A, B, and C in the mathematical model. By substituting the coordinates of point A into the mathematical model and solving the model with a program written based on it, we obtain theoretical results for the coordinates of points B and C. Comparing these with the numerical simulation results, as shown in

Table A1. Comparison between Theoretical and Numerical Results.

g (m/s2)	Bond	Contact Angle	V (mm3)	r2	z2	z1 (mm)		Error Rate	z3 (mm)		Error Rate
						Numerical Results	Theoretical Results		Numerical Results	Theoretical Results
0.05	5	50	2 × 105	49.9	53	46.7	46.3	0.86%	35	35.4	1.14%
1	100	50	2.2 × 105	50	51	50.8	51	0.39%	48.9	48.6	0.61%
10	1000	50	2.5 × 105	50	50.3	50.3	50.3	0.02%	49.2	49.4	0.41%
0.05	5	90	2 × 105	44	59	42.9	43.1	0.47%	37	38	2.70%
1	100	90	2.2 × 105	49.8	52	51.5	52	0.97%	5.1	6.7	2.37%
10	1000	90	2.5 × 105	50	50.7	50.7	50.7	0.02%	48.7	48.6	0.21%
0.05	3.8	50	0.8 × 105	43.8	25.4	27.57	27.01	2.03%	23.25	23.2	0.22%
1	72	50	1 × 105	42.8	24.2	26.33	26.74	1.56%	23.25	23.53	1.20%
10	720	50	1 × 105	43	24.46	25.37	25.23	0.55%	23.97	24.16	0.79%
0.05	4.3	90	0.8 × 105	46.6	32	26.32	26.66	1.29%	19.19	19.68	2.55%
1	81	90	1 × 105	44.8	27.8	28.9	29.5	2.08%	23.5	22.9	2.55%
10	770	90	1 × 105	43.6	25.4	25.9	26	0.39%	23.98	23.9	0.33%

, the accuracy is above 95%. On this basis, by overlaying the theoretically calculated contour with increased transparency on the corresponding cloud map by numerical simulation, as shown in Figure 7, it can be seen that the contours of the liquid surfaces basically coincide. Considering the errors in numerical simulations, we can say that the theoretical results are in good agreement with reality.

6. Further Analysis

6.1. Procedures for Calculating Contours

Two types of programs for predicting liquid surface morphologies have been developed based on the theoretical model. The first program requires inputs of the vertical coordinate of point A and the level of gravity (assuming the inherent properties of the two phases are known). Then, as mentioned earlier, it directly solves the Young–Laplace equation to obtain the corresponding liquid surface contour. Figure 8 illustrates the comparison of surface contours generated by this program under different contact angles and gravity conditions.

The second program requires inputs of the liquid volume and gravity level (with the inherent properties of the two phases known). It then uses the shooting method to determine the coordinates of point A on the liquid surface at that volume of liquid. Similar to the first program, it solves the Young–Laplace equation to obtain the corresponding liquid surface contour. Figure 9 illustrates the comparison of surface contours generated by this program under different contact angles and gravity conditions.

The method presented in this paper is based on the Young–Laplace equation and a simple numerical algorithm, allowing for quick and accurate calculation of liquid surface profiles under any conditions ranging from microgravity to normal gravity. Compared to traditional numerical algorithms, its advantages include significant time savings and adaptability to different gravity conditions. However, a drawback of this method is that it is limited to predicting static and quasi-static situations, and further research is needed to address the prediction of oscillations caused by rapid changes.

6.2. Analyses of Relevant Parameters

Using the developed numerical calculation program, analyses of various surface properties of the liquid inside the tank are conducted. Figure 10, Figure 11, Figure 12, Figure 13 and Figure 14 represent the analyses of various properties of a liquid with a volume of x under different gravity conditions in the tank. In these analyses, H₁ denotes the z-coordinate of the intersection point of the liquid surface with the tank wall, H₂ represents the z-coordinate of the intersection point of the liquid surface with the surface of the central column, Δh is the height difference between the highest and lowest points on the profile of the liquid surface, and S is the area of the gas–liquid interface of the liquid surface, with its calculation formula being the following:

$$S=4\pi {\displaystyle {\int }_{{\Gamma }_1}r}ds=4\pi {{\displaystyle {\int }_{r_2}^{r_1}r\left[1+{\left(\frac{dz}{dr}\right)}^2\right]}}^{0.5}dr$$

(63)

Z_C represents the center of mass of the liquid, with its calculation formula being the following:

$$z_c=\frac{2\pi {r_2}^4{\displaystyle \underset{\frac{r_1}{r_2}}{\overset{1}{\int }}x{[y-f(x)]}^{2}dy}}{V}$$

(64)

when the Bond number is between 1 and 10, significant changes start to occur in various parameters. As the Bond number increases further, the liquid surface tends to become flat, and the interface properties at different contact angles almost converge to the same value. This change marks the transition from a surface tension-dominated interface to a gravity-dominated interface, which aligns with outcomes from numerical simulations. The variations in the levels H₁ and H₂ of liquid surfaces with differing contact angles are notably complex, as Figure 10 and Figure 11 represent. Surfaces with larger contact angles exhibit a descending curve, while those with smaller contact angles show an ascending trend. This occurs because, at low Bond numbers, surface tension is the predominant force. Smaller contact angles facilitate easier ascension along the spherical wall surfaces, which are inherently inclined. Consequently, at low Bond numbers, a liquid surface with a 90-degree contact angle forms a concave shape, which transitions to a flat plane as the Bond number increases, resulting in a rise in H₁. In Figure 12, the height difference of the liquid surface is selected to better illustrate the impact of the contact angle on the liquid surface rise. Under low Bond number conditions, a larger contact angle results in a higher liquid surface rise, with a noticeable difference compared to smaller contact angles. However, this distinction diminishes as the Bond number increases. Figure 13 shows the variation in the area of the gas–liquid interface, a metric that directly determines the influence of surface effects on liquid distribution according to the Young–Laplace equation. This figure visually demonstrates the differences in surface effects of liquids with varying contact angles under low Bond number conditions. Figure 14 calculates the change in the liquid’s center of mass, which is a critical parameter for researchers studying PMD (Propellant Management Device) issues. Liquids with smaller contact angles exhibit better stability in their center of mass.

6.3. Comparison with Other Works

Concus [27] also conducted theoretical calculations of cylindrical liquid surfaces under different Bond numbers, obtaining similar liquid surface profiles and images of characteristic parameters varying with the Bond number. However, his study was limited to a simplistic wall environment, specifically a cylindrical tube with no variations along the z-axis and without a central column. This means that the images produced by Concus are simpler, and the trends in the characteristic parameters are relatively straightforward. Because his research did not account for the influence of different wall shapes on the height of the liquid surface rise, his results cannot be directly applied to most of the current spacecraft tanks.

Chen S.T. [28], also based on the Young–Laplace equation, studied the distribution of liquid in slits within the context of PMD issues, gaining a preliminary understanding of static liquid surfaces in tanks. Chen S.Y. [29] focused on the same engineering background and examined the gas–liquid interface in spherical tanks without a central column under microgravity conditions. They analyzed the contribution of different tank wall shapes and contact angles to the liquid surface rise. However, their studies did not further explore liquid surfaces under varying gravity conditions, nor did they consider the commonly used central column structure in modern plate tanks.

In summary, this study targets the most widely used central column tank structure in current PMD devices, encompassing most gravity conditions that might occur in contemporary aerospace engineering, making it more suitable for practical engineering applications.

7. Conclusions

This paper delves into the gas–liquid interface observed within a tank containing a central column during microgravity experiments conducted in the Chinese Space Station and carries out a detailed study of the state of the gas–liquid interface under different gravity conditions. Based on the Young–Laplace equation and by utilizing perturbation methods, the contour of the gas–liquid interface at Bond numbers between 0.1 and 1 is derived. For the Bond numbers greater than 1, we note that by simplifying the Young–Laplace equation and combining perturbation methods with the boundary layer theory, the contour of the gas–liquid interface at Bond numbers between 1 and 1000 is obtained. The behavior of the liquid under different gravity conditions is simulated using the VOF model. Simulation results are compared with theoretical results from perturbation methods, with an accuracy exceeding 95%.

Two numerical calculation programs based on the perturbation theory are developed. The first program calculates the contour of the gas–liquid interface by inputting the coordinates of the intersection point of the profile of gas–liquid interface with the surface of the central column; the second program determines the contour of the gas–liquid interface by inputting the volume of liquid within the tank. Moreover, based on these programs, an analysis of various properties of the liquid inside the tank under different gravity conditions is conducted, revealing the trends in the behavior of the liquid within the tank as gravity changes.

These findings hold significant importance for the management of tanks in aerospace engineering, such as dealing with the impact of oscillations generated during spacecraft launches and orbital changes on the management of tanks.