Investigation of liquid sloshing characteristics in tanks using an OpenFOAM solver with arbitrary excitation

Liquid sloshing in containers occurs in different fields, such as spacecraft, aircraft, ships, and tank cars. The liquid contained in huge tanks usually sloshes and poses great risks to the equipment. For example, liquid sloshing has significant effects on the management and stability of spacecraft because sloshing forces and torques can interact with the control system through a feedback loop.¹ In addition, to meet the requirement of a higher transportation capacity, the tank volume of large civil aircraft is gradually increasing. Large sloshing caused by braking, rolling, or other situations may result in structural failure and safety problems.² Therefore, a number of studies have been conducted on the characteristics and suppression of liquid sloshing under different conditions.^2–4 

Efforts to explore the dynamic behaviors of liquid sloshing can be classified into three aspects: theoretical studies,^5,6 experiments,⁷ and numerical simulations.^8,9 Abramson⁵ published an analysis of the dynamic behavior of liquids in moving containers, with applications to propellants in space vehicle fuel tanks. Hasheminejad and Soleimani¹⁰ developed theoretical formulations based on linearized potential theory and confirmed that the length and fill depth of a container could strongly affect the sloshing frequency. In experimental studies, Akyıldız et al.¹¹ investigated the liquid sloshing characteristics in a cylindrical tank with various fill levels and ring baffles under the excitation of roll motion. Xue et al.¹² experimentally investigated the effectiveness of four types of baffles in suppressing pressure under a wide range of forcing frequencies. Kim et al.¹³ compared the statistical analysis results of the extremal probability and exceedance probability distributions and found the effectiveness of baffles linked to a spring system in reducing the sloshing magnitude. Liu et al.¹⁴ conducted a series of experiments on the sloshing characteristics of three immiscible liquids in a rectangular tank. The results showed some unique characteristics of three-layer liquid sloshing. The requirement of a complex testing system and high-precision sensors in experimental investigations makes this inefficient. Nevertheless, theoretical and experimental studies have provided feasible methods for analyzing the dynamic behavior of liquid sloshing in containers. For severe sloshing behaviors under arbitrary accelerations, it is difficult for theoretical and experimental methods to sort out due to the high nonlinearity.

With the development of Computational Fluid Dynamics (CFD) technology, simulations using CFD methods have provided a strong impetus for research on liquid sloshing. For free surface flow simulation, the Volume of Fluid (VOF),¹⁵ level-set,^16,17 and front-tracking¹⁸ are regular methods, among which the VOF method can guarantee mass conservation. Liu and Lin¹⁹ used the second-order accurate VOF method to track the distorted and broken free surface under arbitrary six DOF (Degree of Freedom) external excitations. Goudarzi and Sabbagh-Yazdi²⁰ investigated nonlinear sloshing under both harmonic and earthquake induced excitations using the VOF method. They discussed the practical limitations of the linear solution for evaluating the response of seismically excited liquids. Cao et al.²¹ utilized the VOF method to demonstrate the sloshing dynamics of a two-layer liquid under horizontal sinusoidal excitation. It was observed that the interface displacement was mainly influenced by the lower liquid depth, whereas the free-surface displacement primarily depends on the total liquid depth. Kang et al.²² evaluated the nonlinear hydrodynamic pressures on the inner wall of rigid cylindrical tanks resulting from liquid sloshing because of the simultaneous action of two horizontal and one vertical ground motion component. The moving mesh technique was used, and a transport equation considering the deformation control volumes was employed. Wang et al.²³ used the open-source code OpenFOAM to study the interaction between the sloshing flow and random porous structure, and a harmonic roll excitation was considered. For cryogenic liquids, Chung et al.²⁴ performed a multiphase sloshing simulation using the VOF method for liquefied hydrogen in a type C cargo tank, focusing on the impact pressure without considering the thermal effect. Li et al.²⁵ simulated the thermal dynamic behavior of liquid hydrogen in a storage tank for trailers using the VOF method incorporated with a phase-change model to describe the evaporation phenomenon of hydrogen. Nevertheless, the simulations were mostly conducted under regular sloshing conditions, and long-term arbitrary sloshing conditions were rarely considered.

In addition to the traditional CFD method, some new simulation approaches have potential application values in the study of tank sloshing. Gándara et al.²⁶ analyzed the two-dimensional sloshing of water in a partially filled stepped-based tank using an Arbitrary Lagrangian–Eulerian (ALE) method. Jiao et al.²⁷ used the Smoothed Particle Hydrodynamics (SPH) method to simulate the motion of two side-by-side LNG (Liquefied Natural Gas) ships coupled with tank sloshing in regular waves because of the advantage of SPH in dealing with free surface and fluid-structure coupling problems. Cai et al.²⁸ compared the traditional CFD, SPH, and ALE methods for the evaluation of rigid body force in liquid sloshing problems and showed that all the numerical approaches failed to accurately predict the maximum absolute rigid body force while the wave was sharply separated from the main domain owing to the limited resolution. Considering the interface breakup and limitations of the interface-tracking method relying on the mesh resolution, a coupled VOF and Lagrangian particle tracking solver was developed by Heinrich and Schwarze²⁹ to capture drops at small scales, which is an appropriate way to reveal more detailed interface structures of liquid sloshing in tanks.

Previous studies focusing on liquid sloshing were usually solved under regular excitations described by functions or simplified excitation changes, which are not applicable for long-term arbitrary accelerations. For arbitrary accelerations induced by waves, a simulation with fluid-structure coupling can be adopted; however, it is not applicable to spacecraft or aircraft. As for the mesh motion approach, the moving velocity should be derived; additional velocity correction is required in the computation, and it is difficult to track an arbitrarily moving object in post-processing. Consequently, the present study aims to propose a general solver that can directly load arbitrary acceleration data for the simulation of severe liquid sloshing under various conditions.

In the present study, a solver in OpenFOAM is implemented for liquid sloshing to remedy the inconvenience of existing solvers. This solver is based on the interFoam solver and uses the VOF method to capture the free surface. Without mesh motion, an interface is developed to read the data of arbitrary accelerations set in the sheet. In this way, any external excitation, including gravitational acceleration and translational and rotational inertia forces, can be considered. Moreover, this solver is also implemented with algorithms for solving the energy equation and phase change for cryogenic liquids, such as liquid oxygen and LNG in tanks.

The present solver directly solves liquid sloshing by using a time-varying gravity $g⃗t$ in 3D in the momentum equation and implements an algorithm for reading the acceleration data input from users. In this solver, the PIMPLE loop is used for pressure-velocity coupling. Before the PIMPLE loop, an algorithm to refresh the gravity field for every time step is incorporated so that an arbitrary external excitation can be included. The interpolation algorithms are listed in Table I.

The new solver is named tankInterFoam, which is specific to arbitrary liquid sloshing in tanks. A flow chart of the solver is shown in Fig. 1, including the three main parts. The extra code in the green and yellow frames is developed for loading the equivalent external excitation by reading the imported acceleration data and for coupling with the heat and mass transfer at each time step. Compared with the mesh motion method, the algorithm should move the mesh before the iteration in every time step and calculate the relative velocity throughout the entire computational domain in each iteration, thereby consuming additional computational resources. Second, the mesh motion method needs to input the velocity of the object at every time step, but the velocity is usually unknown in actual processes and the acceleration can be easily measured. Therefore, the mesh motion method is insufficient for conditions in which the acceleration is arbitrary and the unknown object velocity needs to be integrated. Conversely, the acceleration values can be directly applied to the present solver, and only the recalculation of the body force is required at every time step, making the simulation more effective and suitable for arbitrary sloshing conditions.

The present simulation sets the excitation amplitude and frequency to 0.01 m and 0.94 ω₁, respectively, based on the experimental configuration, where ω₁ is the first-mode natural frequency equal to 6.403 rad/s. Three probes are located at distances of 65, 105, and 185 mm above the bottom of the tank to monitor the time series of the dynamic impact pressure. A schematic of the simulated tank and mesh configuration is shown in Fig. 2.

As shown in Fig. 2, a rectangular tank with dimensions of 570 × 310 × 570 mm³ is constructed. The boundary conditions for all six faces of the tank are set as no-slip walls. The entire computational domain is uniformly divided into 50 × 25 × 50 square cells, and the total number of cells is 62 500. The water is initially set at the bottom of the tank, and the depth is 180 mm. Other numerical configurations and physical attributes are listed in Table II. For the turbulence model, the simulation is conducted using the laminar model, the accuracy of which is validated by Cai et al.²⁸ for similar sloshing cases as it obtains the minimum impulse overprediction among the investigated turbulence models using a relatively coarse mesh and requires considerably less calculation time.²⁸ 

Figure 3 compares the free surface structures captured at the same instants in the experiment and simulation, showing good consistency. The slosh wave breaking process occurs close to the tank wall, and lots of bubbles around the impact zone are trapped when the sloshing wave impacts on the tank wall. The results show that the time-varying gravity field can replace the tank motion and capture the liquid sloshing characteristics well under external excitation.

The time series of the dynamic impact pressures detected by the three probes are shown in Fig. 4. Two peaks can be found in the hump. According to Jiang et al.,³⁶ the first peak is caused by the initial impact of the sloshing liquid on the wall, and the second peak is mainly produced by the following water in the case of gravity after the previous impact. The comparison also shows good consistency, although differences are found in the amplitude.

The evaporation and condensation of cryogenic liquids in the container is important. In this section, the effect of sloshing on the evaporation and condensation of the liquid is investigated. To test the effectiveness of the phase-change component of the tankInterFoam solver, two imaginary simulations with different excitation amplitudes are conducted. The frequency is 0.94ω₁, and the amplitudes are set to 0.005 and 0.001 m for large and small sloshing, respectively. The saturation temperature of LNG is set to 111 K. The initial temperature of LNG is set to 90 K, and the wall temperature is set to 115 K. In this case, the LNG evaporates owing to the heat transfer between the liquid and the wall. The simulation parameters and physical properties for this case are listed in Table III. The computational domain and mesh configuration are the same as those shown in Fig. 2. The walls are also set as no-slip walls, and the laminar model is employed.

The time series of the volume-averaged evaporation rates under different sloshing conditions is shown in Fig. 5. At the beginning, the evaporation rate under low-amplitude sloshing is larger than that under high-amplitude sloshing. However, the evaporation rate under higher-amplitude sloshing increases faster, and at ∼80 s, the evaporation rate under high amplitude sloshing becomes larger than that under low amplitude sloshing. The reason for this phenomenon is discussed in Sec. III C, which shows the evolution of the temperature field.

Figure 6 shows the evolution of the temperature field under different sloshing conditions. We chose two instants at ∼30 and 80 s for comparison. The temperature first increases on the wall and free surface, and then the high temperature enters the liquid. As shown, the heat transfer under low-amplitude sloshing is slower than that under high-amplitude sloshing, so the heat concentrates around the interface, and it is found that the temperature around the center of the interface under low-amplitude sloshing is higher. In this situation, some of the parts may reach the saturated temperature and start evaporating under low-amplitude sloshing, whereas the temperature under high-amplitude sloshing reaches the saturated temperature later. At ∼80 s, the temperature in most regions is higher under larger sloshing and reaches the saturated temperature so that the evaporation rate is higher under higher-amplitude sloshing after 80 s, as seen in the curves in Fig. 5.

In engineering, baffles are typically used to inhibit liquid sloshing in containers. In this section, a regular cylindrical tank with a baffle set in the middle is utilized to investigate the characteristics of liquid sloshing in a tank with a baffle. The width of the baffle W is 50 mm, and the thickness of the baffle h is 10 mm. In this section, the developed solver is used in a mild sloshing situation to investigate the effects of the baffle on liquid sloshing. A cylindrical tank with a diameter D = 916 mm and a height H = 2000 mm is taken into consideration, and the structure of the tank with a baffle is shown in Fig. 7(a). The computational domain is divided into hexahedral cells as shown in Fig. 7(b). The cell size is ∼10 mm in the entire domain, and the total cell number is ∼1 400 000. Because the main focus is on the change in the liquid surface near the baffle during sloshing, the computational mesh is refined above and below the baffle. Other numerical configurations and physical attributes are listed in Table IV.

Five conditions with different liquid heights are taken into account, with different distances from the free surface to the center of the baffle (d_s). Liquid sloshing is conducted with an initial condition in which the free surface is settled under a 0.1g acceleration in the X direction and a −1.0g acceleration in the Z direction. Sloshing is simulated under a −1.0g acceleration in the Z direction. The dimensionless distances from the free surface to the center of the baffle (d_s/R₀) in the present work are 0, 0.05, 0.1, 0.15, and 0.2, where R₀ is the radius of the tank.

From the comparison of the damping ratios among the experiment, simulation, and empirical formula in Fig. 8(b), it can be observed that the simulated results confirm well with the experimental results for different conditions. The empirical results show a different trend when d_s/R₀ ≤ 0.1, as they fail to consider the overall factors. From the experimental results, it can be found that the damping ratio around d_s/R₀ = 0.075 reaches a maximum value. The damping ratios obtained by the present simulation are slightly larger than the experimental results when d_s/R₀ ≤ 0.1, which may be due to the idealization in the simulation and the error in calculating the damping ratio. It is interesting that the damping ratio reaches a maximum value not when the height of liquid is equal to the height of the baffle but when the free surface is slightly higher than the baffle.

To analyze the reason for the above-mentioned deviation, the free surface evolutions simulated for d_s/R₀ = 0.0 and d_s/R₀ = 0.1 in the present work are shown in Fig. 9. In the figure, the motion of the free surface is marked by arrows, where the red arrows indicate that free surface movement is restricted by the baffle. When d_s/R₀ = 0.0, it can be found that the free surface on the right-hand side is moving up freely without the effect of the baffle, while the free surface on the left-hand side is affected by the baffle from 20 to 21 s. When the baffle is immersed in the liquid when d_s/R₀ = 0.1, as shown in Fig. 9(b), the effect of the baffle always acts during the entire sloshing cycle. Therefore, it can be inferred that when the baffle is slightly immersed in the liquid and is within the intense effect region of the baffle, the damping ratio is the maximum. As the liquid height further increases and the effect of the baffle weakens, the damping ratio decreases. Therefore, the maximum damping ratio occurs when the baffle is slightly immersed in the liquid, which rectifies the empirical result that the maximum damping ratio occurs when the liquid height is equal to the baffle height.

The liquid sloshing frequencies under different conditions after FFT spectrum analysis are shown in Fig. 10. It can be observed that when d_s/R₀ > 0.1, the frequencies are ∼0.305 Hz and change slightly. By decreasing the distance from the free surface to the baffle, the frequency increases, and it is ∼0.3389 Hz when d_s/R₀ = 0.0. For d_s/R₀ = 0.05, although the damping ratio is larger, the frequency is smaller compared with that of d_s/R₀ = 0.0. For d_s/R₀ > 0.05, the sloshing is mainly controlled by the resistance of the baffle, and with a decrease in the resistance of the baffle, the frequency decreases with decreasing the damping ratio. For d_s/R₀ = 0.0, the effect of the baffle is slightly different from that of the other conditions, and the free surface is mainly controlled by the inner circle of the baffle. Therefore, the frequency of d_s/R₀ = 0.0 is larger than that of d_s/R₀ = 0.05, although the damping ratio is smaller than that of d_s/R₀ = 0.05.

The impact of baffle thickness h is further investigated in the present work. Five different baffle thicknesses are selected and simulated considering d_s/R₀ = 0. The evolutions of the sloshing forces in the X direction are shown in Fig. 11(a), and the damping ratios are shown in Fig. 11(b). According to the simulation results, as the baffle thickness increases, the baffle damping ratio decreases in an approximately linear manner. The damping ratio for a baffle thickness of 5 mm is ∼0.044, which is 4.7% higher than that for a thickness of 10 mm. The damping ratio for a baffle thickness of 15 mm is about 0.039, which is 7.7% lower than that for a baffle thickness of 10 mm. Therefore, it can be concluded that the baffle thickness is inversely related to the damping ratio, suggesting that the baffle thickness can be designed to be as thin as possible while meeting the strength requirement.

In spacecraft or during the voyage of ships, the fuel in the container may slosh violently, which cannot be described by theoretical analysis. Using the tankInterFoam solver, arbitrary input data, including mutational overloads, can be easily input by users. In this section, a regular cylindrical container containing LNG is used for the simulation. The diameter and height of the container are 0.916 and 2.0 m, and the initial fill rate is 50%. This simulation imagines a scenario that is similar to the situation of rocket launching and separation. The time series of accelerations in three directions are shown in Fig. 12.

From 0 to 20 s, the increase in acceleration in the Z direction corresponds to the process of axial acceleration, such as launching. At 20 s, the acceleration in the Z direction suddenly disappears, indicating the proximity of the main thruster. Subsequently, the acceleration in the Z direction changes direction, and the acceleration in the X direction becomes the main external force owing to the change in attitude. In this process, the liquid in the tank should violently roll and cause highly nonlinear sloshing forces, which should be simulated.

With the tankInterFoam solver, these complicated accelerations can be easily imported without functions. Figure 13 shows the simulated structures of the free surface under drastically changing acceleration. In the simulation, the liquid in the tank is assumed to be LNG with a density of 0.422 kg/m³ and a temperature of 111 K. As the acceleration in the Z direction is released and the acceleration in the X direction increases, the liquid rolls upward along the left-hand side wall. Then, the liquid concentrates in the upper left corner of the tank after several cycles of large-amplitude sloshing. At 23 s, we can find that the free surface is greatly fragmentated owing to the impact on the wall, and at 26 s, the free surface is also fragmented owing to the restoration of liquid after the first excessive sloshing.

Figure 13(b) shows the temperature distribution on the free surface. At the first moment of ∼20 s, the liquid over the walls has a higher temperature because mild sloshing causes the liquid around the wall to have a more adequate heat exchange. When fierce motion occurs, the free surface sinks into the liquid and breaks up. The liquid carrying a high temperature cools quickly, and the free surface shows a relatively consistent heat distribution. At the end of the violent sloshing, the free surface tends to stabilize, and the liquid near the wall rises again.

Figure 14 shows the free-surface evolution at different instants to determine the influence of the baffle on violent sloshing. At t = 21 and 22 s, a significant delay in the free surface change can be observed with the effect of the baffle. At t = 27 s, the free surface is also relatively leveled off because the baffle can also restrain sloshing in the Z direction.

Figure 15 shows the evolution of the forces on the container wall monitored during violent sloshing in containers with and without a baffle. The trend of force evolution conforms to the change in acceleration because the force mainly acts in the same direction as the acceleration. As the acceleration mainly changes in the Z direction, the fluctuation in the Z-direction force is relatively violent. The force values also coincide with the theoretical values according to Newton’s second law. Generally, as seen in the time-domain diagram from t = 0 s to t = 20 s, the relatively steady loads cause liquid sloshing regularly in all directions in the tank without the baffle. Conversely, regular sloshing is suppressed in the tank with the baffle owing to its relatively high damping ratio. After 20 s, violent sloshing in the Z direction is found slightly delayed in the tank with a baffle. Because the Y direction sloshing is suppressed from t = 0 s to t = 20 s, the sloshing in the Y direction after 20 s is also disappeared in the tank with a baffle as it is induced from the regular sloshing from 0 to 20 s.

To figure out the influence of the baffle on liquid sloshing in the container, a comparison of the force evolution in tanks with and without a baffle is shown in Fig. 16. The curves are enlarged to show the detailed force evolution during the regular and violent sloshing stages. From 0 to 20 s, due to the slight change in the overload direction, the liquid in the tank without baffle sloshes regularly and the amplitude gradually increases, which may cause equipment oscillation. When the baffle is set at the position of the free surface, liquid sloshing is effectively suppressed.

As seen from the force evolution curves from 19 to 35 s, the red dashed box indicates that the sloshing in the Y direction, which continues from the regular sloshing, is weakened. For sloshing in the X and Z directions, it can be found that the forces in the X and Z directions are enlarged first as the baffle’s interaction with the liquid and the sloshing is delayed by the baffle, but the inhibition effect is negligible when the free surface is nearly perpendicular to the baffle.

In the present work, a solver for solving liquid sloshing problems in tanks is developed within the open source CFD library OpenFOAM v2206. The developed solver tankInterFoam is modified from the interFoam solver by implementing an interpolating algorithm for arbitrary external excitations and a heat and mass transfer algorithm for cryogenic liquids. The developed solver provides an easy method to simulate liquid sloshing in tanks and can be widely applied to reveal the characteristics of liquid sloshing in tanks under arbitrary excitations. Simulations under different conditions have been conducted based on the solver, and some conclusions can be drawn as follows:

• The characteristics of the liquid in the tank under sinusoidal excitation are used to verify the morphological characteristics of the free surface. The amplitude and frequency characteristics of the forces on the wall are compared with the experimental results, showing remarkable consistency of the solver.

• The energy equation with phase change is incorporated into the solver, and the liquid evaporation and condensation in the LNG storage tank are simulated under two sloshing conditions with different amplitudes, showing that the evaporation rate is lower at the beginning and increases much faster in the later stage under higher amplitude sloshing.

• From the simulations with different distances between the interface and the baffle, the damping ratio is found to reach a maximum value when the baffle is slightly immersed in the liquid at d_s/R₀ = 0.075 in the present work and then decreases with increasing the distance. Through the simulation of violent sloshing under arbitrary acceleration, it is found that the baffle delays the peak force and allows sloshing to subside more quickly.

• From the simulations with different baffle thicknesses from 5 to 15 mm, the damping ratio is found to decrease with increasing the baffle thickness approximately linearly. In this study, the damping ratio for a baffle thickness of 5 mm is 4.7% higher than that for a thickness of 10 mm. The damping ratio for a baffle thickness of 15 mm is ∼7.7% lower than that for a baffle thickness of 10 mm.

The supplementary material encompasses the solver code for liquid sloshing in tanks supporting arbitrary excitations developed using OpenFOAM.

This research was funded by the National Natural Science Foundation of China (Grant Nos. 52372366 and U23A20669) and the Key Research and Development Program of Zhejiang Province (Grant No. 2024C01168).

The authors have no conflicts to disclose.

Linmin Li: Investigation (equal); Methodology (equal); Software (equal); Writing – original draft (equal). Bohan Shen: Formal analysis (equal); Visualization (equal). Zuchao Zhu: Funding acquisition (equal); Project administration (equal). Qile Ren: Funding acquisition (equal); Project administration (equal). Junhao Zhang: Funding acquisition (equal); Project administration (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.