Full Dynamic Model for Liquid Sloshing Simulation in Cylindrical Tank Shape

Omar Noui¹, Mohamed Bouazara^1*, Marc J. Richard²

¹Department of Applied Science, University of Quebec at Chicoutimi, Chicoutimi, Canada.
²Department of Mechanical Engineering, Laval University, Quebec, Canada.

DOI: 10.4236/wjm.2024.144004                

Abstract

This study presents a comprehensive full dynamic model designed for simulating liquid sloshing behavior within cylindrical tank structures. The model employs a discretization approach, representing the liquid as a network of interconnected spring-damper-mass systems. Key aspects include the adaptation of liquid discretization techniques to cylindrical lateral cross-sections and the calculation of stiffness and damping coefficients. External forces, simulating various vehicle maneuvers, are also integrated into the model. The resulting system of equations is solved using Maple Software with the Runge-Kutta-Fehlberg method. This model enables accurate prediction of liquid displacement and pressure forces, offering valuable insights for tank design and fluid dynamics applications. Ongoing refinement aims to broaden its applicability across different liquid types and tank geometries.

Keywords

Fluid-Structure Interaction, Equivalent Mechanical Model, Liquid Discretization, Spring-Mass Model, Spring-Mass Network, Liquid Simulation

1. Introduction

Over the past few years, there has been a considerable focus on studying the stability of tank-trucks, with particular attention given to the impact of liquid motion within the tank on its dynamic behavior. Numerous analytical and numerical models have been developed to assess both static and dynamic forces exerted by the liquid and their consequent effects on the vehicle’s behavior. However, these sophisticated models necessitate complex mathematical and computational tools, as well as expensive simulation software. Conversely, experimental studies on this subject have proven to be costly and delicate. Consequently, some researchers have been motivated to devise new and simpler methods for studying liquid motion in partially filled tanks. The primary approach involves substituting the movement of the liquid with the oscillations of mechanical models, such as spring-mass systems or pendulums. Several 2D mechanical models have been created to simulate liquid sloshing in tanks [1] [2] [3] [4]. Nevertheless, these models overlook various crucial parameters, particularly in the case of portable tanks. Forces acting on the tank can originate from any direction, leading to an underestimation of liquid motion and pressure forces applied to the tank walls, potentially resulting in design flaws. This inherent limitation renders plan mechanical models ineffective for accurately simulating liquid motion in practical scenarios.

Mechanical models find applications across various fields for simulating complex phenomena. Specifically, in graphic computer animation, spring-mass systems are commonly employed due to their simplicity and ease of implementation [5] [6] [7]. Furthermore, utilizing these models in deformable body simulations offers a combination of a geometric representation that aligns with the object’s topology and a physical representation that models internal and external interactions.

This study focuses on developing a new full dynamical model to simulate the overall motion of liquid in a moving horizontal cylindrical tank, irrespective of the filling rate and liquid type. The model can simulate lateral, longitudinal, and vertical displacements, as well as evaluate pressure forces applied to the walls. The central concept of the 3D model involves representing the liquid as a mesh of spring-mass-damper systems. The process begins by dividing the liquid in the tank into multiple masses along each axis. The movement of each mass is simulated by displacing its center of mass, which forms the nodes of the mesh. Flexible edges connecting adjacent nodes incorporate parallel springs and dampers.

This study outlines the discretization method for the liquid based on the cylindrical tank shape. It details the computation of masses and initial coordinates for each node, followed by obtaining the global equations of motion for the developed model. Subsequently, steps are taken to determine other parameters of the model, such as the stiffness of the springs and damping coefficients of the dampers. Finally, the equations are solved for various conditions, and the results are compared with existing literature. Whether for the liquid displacement or for the pressure forces applied on the tank walls, the results are well comparable to the results obtained by more complex numerical models. Thus, this model allows to save calculations time while ensuring precision in the computation.

2. Liquid Discretization

There are several conventional tank shapes for road vehicles, with circular, elliptical, and modified oval lateral cross-sections being the most commonly utilized. The selection of these shapes is not arbitrary; rather, elliptical and modified oval sections are chosen due to their lower center of mass compared to cylindrical tanks. This design choice contributes to enhanced vehicle stability and a reduced rollover threshold.

However, it’s important to note that these shapes come with trade-offs. While elliptical and modified oval sections lower the center of mass, providing stability benefits, they also result in a wider free surface of the liquid, leaving more space for sloshing. This may potentially increase the rollover threshold, working against the desired behavior of the vehicle. Despite these considerations, cylindrical tanks remain the most prevalent in the industry.

This paper focuses on studying the cylindrical tank, where the lateral walls are assumed to be straight without any camber. The cylindrical shape continues to be widely used, and in this study, the specific characteristics and dynamics of such tanks are explored. Figure 1 below illustrates the circular lateral cross-section and the chosen reference frame for this investigation:

Figure 1. Circular lateral cross-section of the tank.

Each section is defined by its center O' and its radius R. The length of the tank is noted by L. Equation of the section is given by the formula:

$y^2+{\left(z-R\right)}^2=R^2$ (1)

The first step in liquid discretization is to compute the height $Z_h$ of the free surface, according to the filling rate $\epsilon$ ($0<\epsilon <1$ ). In this case, its value was obtained by solving the following equation:

${\displaystyle {\int }_0^{Z_h}f\left(z\right)\text{d}z}=\frac{\epsilon }{1-\epsilon }\left[{\displaystyle {\int }_{Z_h}^{2R}f\left(z\right)\text{d}z}\right]$ (2)

where $f\left(z\right)=\sqrt{R^2-{\left(z-R\right)}^2}$ .

In this study, we adapted the liquid discretization model outlined in a previous study [8] to accommodate a circular lateral cross-section shape. The objective of this model is to partition the liquid shape along each axis, resulting in numerous liquid portions referred to as particles. The centroids of these particles serve as the nodes of the primary model. Notably, the liquid’s length (X-axis) is divided by M, its width (Y-axis) by 2N, and its height (Z-axis) by P.

The subsequent step involves determining the mass and initial coordinates of the center of mass for each particle. To accomplish this, various scenarios were investigated corresponding to the selected discretization method. The calculation of masses and coordinates for each particle depends on its precise location within the tank. Detailed procedures and formulas can be found in our previous study [8].

Following this, we presented two examples of the liquid discretization model. The first example (Figure 2 and Figure 3) illustrates a cylindrical tank filled to

Figure 2. Example 1 for liquid discretization (2D).

Figure 3. Example 1 for liquid discretization (3D).

70% with M = 15, N = 4, and P = 7. For these parameters, we obtained 840 nodes, as depicted in the subsequent figures.

The second example (Figure 4 and Figure 5) shows a 50% filled cylindrical tank with M = 30, N = 6 and P = 10. We obtain 7560 nodes for these parameters as follows:

Figure 4. Example 2 for liquid discretization (2D).

Figure 5. Example 2 for liquid discretization (3D)

3. Stiffness and Damping Coefficients

To simulate liquid sloshing, the flexible edges connecting the nodes are composed of a parallel spring and damper. Indeed, each node relates to six edges, two edges along each axis except those on the free surface; they are connected to 5 edges. Centers of mass of the particles contacting the tank walls are connected to the structure using the same type of edges. To compute the spring stiffnesses along the X-axis and Y-axis, we utilized and adjusted the formula proposed by Dodge [9], which calculates the natural frequencies of each level of the liquid in a rectangular tank

${\omega }_n^2=\pi \left(2n-1\right)\frac{g}{a}\tanh \left[\pi \left(2n-1\right)\frac{h}{a}\right]$ (3)

In adapting the formula to a cylindrical shape, where “a” represents the tank width, “h” signifies the liquid height, and “g” denotes the acceleration due to gravity, we derive the natural frequency for each level of the liquid, allowing us to calculate the global stiffness associated with each level. We assume that the springs at each level are in series and that the stiffness increases as the spring approaches the tank walls.

Subsequently, to compute the spring stiffnesses along the Z-axis, we establish that the maximum vertical displacement remains fixed for each filling rate, with the values obtained from literature [10]. This assumption facilitates the calculation of the global vertical spring stiffness. Moreover, we presume that the springs at each level of the liquid are in parallel, with diminishing significance of movement as the liquid approaches the bottom of the tank. Consequently, the global vertical stiffness is distributed considering these assumptions.

For detailed procedures and formulas to compute all stiffnesses, please refer to our previous study [11].

Furthermore, the damping coefficients are assumed to be equal to one hundred times the liquid’s dynamic viscosity coefficient, denoted by η. Given that the new full dynamical model represents the liquid using a set of solid bodies, higher damping coefficients are required to restrain the movement of these bodies effectively. Through experimentation with various values, it was observed that results were more accurate when η was multiplied by 100. Consequently, the damping coefficient of each node, denoted as G(i,j,k), is computed using the following formula:

$c_{1,i,j,k}=100\eta$ (4)

Alternatively, in light of all the assumptions made for computing stiffnesses, we employed a calibration factor to multiply longitudinal and lateral stiffnesses. The magnitude of this factor is contingent upon the number of nodes utilized in the discretization model. Specifically, the greater the number of nodes in the model, the lesser the necessity for a calibration factor.

4. Equations of Motion

To formulate the equations of motion for the new model, the analysis began with an examination of the forces acting on each node. The motion of an individual node is governed by three equations. Each node experiences three types of forces: the stresses from the springs (five or six forces), the damping forces (five or six forces), and external forces. It’s important to note that the initial length of each spring corresponds to the initial distance between its adjacent nodes. External forces can simulate various vehicle maneuvers, such as lane changes, turning, or braking.

The global system of equations for the model is derived by assembling the equations for each node. This global system of equations can be represented as follows:

$\left[ℳ\right]\left(\ddot{r}-a_{ext}\right)+A-T=0$ (5)

where $\left[ℳ\right]$ represents the mass matrix, $\ddot{r}$ is the vector of nodes accelerations, $T$ is the vector of all spring stresses, $A$ is the vector of damping forces and $a_{ext}$ is the vector of external acceleration.

Subsequently, Maple Software [12] is utilized to solve the global system of equations employing the Runge-Kutta-Fehlberg method, RKF45 [14]. The coordinates of the liquid’s center of mass and the forces generated by its movement are computed based on the displacements and velocities of each node obtained by solving the equations.

5. Results

5.1. Lateral Sloshing

In this section, some results in terms of lateral liquid displacement are given for a 50% and 70% filled cylindrical tank with a length X = 7.5 m, a radius r = 1.02 m. The liquid used in this simulation is the domestic oil with a density $\rho =960\text{kg}\cdot {\text{m}}^{-3}$ and a dynamic viscosity $\eta =0.048\text{kg}\cdot {\text{m}}^{-1}\cdot {\text{s}}^{-1}$ . The selected values for the discretizing model are M = 12, N = 2 and P = 5. Therefore, the dynamical model consists of 240 nodes and 390 links. The simulation consists of displaying the model response to lateral forces applied on each node caused by movement in a curve. This acceleration is defined as follows:

$\{\begin{array}{ll}{a}_y=0\hfill & t\le t_1\hfill \\ a_y=-\frac{A}{2}\cos \left(2\left(t-t_1\right)+\pi \right)+t_2\hfill & t_1<t\le t_2\hfill \\ a_y=A\hfill & t<t_2\hfill \\ A=3\text{m}/{\text{s}}^{\text{2}},t_1=0.1\text{s},t_2=0.6\text{s}\hfill & \hfill \end{array}$ (6)

Firstly, the results of liquid simulation are shown for a 50% filled tank. Figure 6 shows the lateral movement of the nodes in the middle of the tank according to time t where $1\text{s}\le t\le 5.5\text{s}$ with a step of 0.5s. The blue line represents the liquid free surface. This curve is a natural cubic spline obtained by a kriging method [13]. In this case, the interpolation points are the upper nodes of the model. The purpose is to visualise the deformation of the liquid free surface and the waves generated by the liquid sloshing.

The following figures (Figure 7 and Figure 8) show lateral and vertical displacement of the liquid center of mass for $0\text{s}\le t\le 8\text{s}$ due to the lateral acceleration, As shown, the maximum lateral displacement of liquid center of mass is 11.4 cm. Then, the lateral displacement seems to be stabilizing around 9 cm. The maximum vertical displacement of the liquid center of mass from its initial position is 4.7 cm then the vertical displacement seems to be stabilizing around 2 cm.

The following graphics show lateral and vertical pressure forces generated by liquid sloshing using the liquid center of mass acceleration. Figure 9 shows that

Figure 6. Lateral liquid sloshing of 50% filled tank for $1\text{s}\le t\le 5.5\text{s}$ .

Figure 7. Lateral displacement of the liquid in a 50% filling rate.

Figure 8. Vertical displacement of the liquid in a 50% filling rate.

Figure 9. Lateral pressure forces generated by liquid sloshing for 50% filling rate.

lateral force generated by liquid sloshing can reach a maximum of 37 KN and it decreases successively by the effect of the dynamic viscosity. Also, Figure 10 shows that vertical force generated by liquid sloshing can reach a maximum of

Figure 10. Vertical pressure forces generated by liquid sloshing for 50% filling rate.

42 KN. Vertical pressure forces are not negligible, and it’s important to take them into account.

Secondly, the same results of liquid simulation are shown for a 70% filled tank. Figure 11 shows the lateral movement of the nodes in the middle of the tank according to time t where $1\text{s}\le t\le 5.5\text{s}$ with a step of 0.5 s. Comparing results of lateral liquid sloshing between 50% and 70% filling rates, it is found that it is less deformation of the free surface in the case of 70% filling rate. This is related to the smaller area of the free surface and a heavier liquid charge. It is possible to make the same observation by comparing lateral and vertical displacement of the global center of mass of the liquid.

The following figures show lateral and vertical displacements of the liquid center of mass for $0\text{s}\le t\le 8\text{s}$ due to the lateral acceleration in the case of 70% filling rate. As shown in Figure 12, the maximum lateral displacement of the liquid center of mass is 8.8 cm. Then, the lateral shift seems to be stabilizing around 6.4 cm, about 2.6 cm less than the 50% filling rate case. For, the vertical shift (Figure 13), the maximum displacement of the liquid center of mass from its initial position is about 4 cm then the vertical displacement seems to be stabilizing around 2.5 cm, about 0.5 cm higher than the 50% filling rate case, which is negligible. It is concluded that the lateral shift of the liquid charge is higher when the tank is 50% filled.

The following graphics show lateral and vertical pressure forces generated by liquid sloshing using the liquid center of mass acceleration where the tank is 70% filled. Figure 14 shows that lateral force generated by liquid sloshing can reach a maximum of 47 KN. We note that lateral pressure forces increase by 21% compared to the 50% filling rate case. Also, Figure 15 shows that vertical pressure forces can reach a maximum of about 52 KN which is 19% higher than vertical pressure forces generated with a 50% filling rate. This results from the fact that the quantity of transported load is higher in the case of 70% filling rate. Despite that the center of mass displacement is smaller when the tank is 70% filled, pressure forces are greater because of the mass of the liquid in motion.

Figure 11. Lateral liquid sloshing of 70% filled tank for $1\text{s}\le t\le 5.5\text{s}$ .

5.2. Longitudinal Sloshing

In this simulation, a longitudinal force is applied to the liquid in a 50% filled cylindrical tank. Dimensions of the tank are X = 4.5 m and r = 1.05 m. The liquid used in this simulation is the domestic oil with a density $\rho =960\text{kg}\cdot {\text{m}}^{-3}$ and a dynamic viscosity $\eta =0.048\text{kg}\cdot {\text{m}}^{-1}\cdot {\text{s}}^{-1}$ . The selected values for the discretizing model are M = 16, N = 2 and P = 3. Therefore, the dynamical model consists of

Figure 12. Lateral displacement of the liquid in a 70% filling rate.

Figure 13. Vertical displacement of the liquid in a 70% filling rate.

Figure 14. Lateral pressure forces generated by liquid sloshing for 70% filling rate

Figure 15. Vertical pressure forces generated by liquid sloshing for 70% filling rate.

192 nodes and 255 links. Longitudinal deceleration simulating the vehicle braking is given by the following formula:

$\{\begin{array}{ll}{a}_x=0\hfill & t\le t_1\hfill \\ a_x=-A\cos \left(\pi \left(t+t_2-2t_1\right)\right)\hfill & t_1<t\le t_2\hfill \\ a_x=A\hfill & t>t_2\hfill \\ A=2.5\text{m}/{\text{s}}^{\text{2}}\text{or}4.5\text{m}/{\text{s}}^{\text{2}},t_1=1\text{s},t_2=1.5\text{s}\hfill & \hfill \end{array}$ (7)

To study the effect of braking forces on the liquid movement, two acceleration amplitudes are compared in this section, 2.5 m/s² and 4.5 m/s². The following figures (Figure 16 and Figure 17) show X and Z coordinates of the liquid center of mass. These coordinates are obtained by using the instantaneous coordinates of the model nodes. We notice that a higher acceleration amplitude makes higher longitudinal and vertical displacements, as for the amplitude of oscillations.

Figure 16. Longitudinal displacement of the liquid center of mass. (a) 2.5 m/s²; (b) 4.5 m/s².

Figure 17. Vertical displacement of the liquid center of mass. (a) 2.5 m/s²; (b) 4.5 m/s².

Figure 18. Shape of the liquid free surface (A = 2.5 m/s²).

Figure 18 and Figure 19 provide a 3D overview of the liquid free surface for the two acceleration amplitudes. This surface is obtained by creating 3D natural cubic splines using a multidimensional kriging method [14]. Interpolation

Figure 19. Shape of the liquid free surface (A = 4.5 m/s²).

points are the nodes of the upper nodes of the new model. The images show the shape of the free surface between t = 1 s and t = 8 s at each second. It is clearly seen that the free surface is deformed in the front of the tank because of longitudinal forces generated by the movement of the other nodes. Moreover, a higher amplitude of acceleration generates a greater deformation of the free surface.

Figure 20 and Figure 21 depict the longitudinal and vertical pressure forces, respectively, resulting from the movement of the liquid center of mass. Notably, it can be observed that pressure forces escalate with the initial deceleration. In both directions, the amplitude of the force more than doubled when the deceleration increased from 2.5 m/s² to 4.5 m/s².

In conclusion, the simulations presented in this section represent preliminary results that successfully validate the effectiveness of the new 3D dynamical model. The obtained results exhibit a close resemblance to the literature concerning liquid displacements and pressure forces generated by liquid motion, particularly when comparing filling rates of 50% and 70% [15]. Moreover, the new model stands out by offering a visual representation of liquid motion through the modeling of the free surface using a cubic spline, providing an advantage over other

Figure 20. Longitudinal force of the liquid center of mass. (a) 2.5 m/s²; (b) 4.5 m/s².

Figure 21. Vertical displacement of the liquid center of mass. (a) 2.5 m/s²; (b) 4.5 m/s²

equivalent mechanical models. It is noteworthy that the computational efficiency of this new model is commendable, as the presented results were achieved in less than ten minutes using only a laptop for computation.

6. Conclusions

In this study, we present a novel full dynamical model tailored for simulating liquid behavior within a cylindrical portable tank. The model entails discretizing the liquid and portraying it as a network of spring-damper-mass systems. Initially, we formulate a specific liquid discretization model to determine the nodes of the dynamical model. Each neighboring pair of nodes is linked by a parallel spring and damper. Subsequently, we incorporate certain assumptions to compute the stiffness and damping coefficients of the model.

Following this, we develop the global system of equations of motion. We then present preliminary results to validate the model, comparing it to existing literature. The results of the liquid simulation closely align with the literature, demonstrating high accuracy concerning liquid displacement and pressure forces generated by load shifts. Moreover, the new full model provides a visually realistic representation of the liquid motion.

In future work, we will present and discuss additional results, incorporating modifications to different parameters such as the type of liquid, vehicle maneuvers, and the number of nodes. This ongoing exploration aims to further refine and expand the applicability of the proposed full dynamical model to develop new software in this domain.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Abramson, H.N. (1966) The Dynamic Behavior of Liquids in Moving Containers. NASA SP-106.
[2]	Ranganathan, R., Ying, Y. and Miles, J.B. (1993) Analysis of Fluid Slosh in Partially Filled Tanks and Their Impact on the Directional Response of Tank Vehicles. SAE Technical Paper 932942, 39-45. https://doi.org/10.4271/932942
[3]	El Damatty, A.A. and Sweedan, A.M.I. (2006) Equivalent Mechanical Analog for Dynamic Analysis of Pure Conical Tanks. Thin-Walled Structures, 44, 429-440. https://doi.org/10.1016/j.tws.2006.03.016
[4]	Noui, O. (2012). Étude et modélisation des citernes en aluminium pour les véhicules routiers. Mémoire de maîtrise, Université du Québec à Chicoutimi. https://doi.org/10.1522/030341785
[5]	Provot, X. (1995) Deformation Constraints in a Mass-Spring Model to Describe Rigid Cloth Behavior. National Institute for Research in Computer Science and Control (France).
[6]	Balci, M. and Foroosh, H. (2006) Real-Time 3D Fire Simulation Using a Spring-Mass Model. 2006 12th International Multi-Media Modelling Conference, Beijing, 4-6 January 2006, 8. https://doi.org/10.1109/mmmc.2006.1651309
[7]	Mesit, J., Guha, R. and Chaudhry, S. (2007) 3D Soft Body Simulation Using Mass-Spring System with Internal Pressure Force and Simplified Implicit Integration. Journal of Computers, 2, 34-43.
[8]	Bouazara, M., Noui, O. and Richard, M.J. (2014). 3D Dynamical Model for Liquid Motion Simulation in a Partially Filled Tank. Volume 4A: Dynamics, Vibration, and Control, Montreal, 14-20 November 2014. https://doi.org/10.1115/imece2014-36211
[9]	Dodge, F.T. (2008) The New Dynamic Behaviour of Liquids in Moving Containers. Southwest Research Institute (SRI), San Antonio, Texas.
[10]	Kang, X. (2001) Optimal Tank Design and Directional Dynamic Analysis of Liquid Cargo Vehicles under Steering and Braking. Ph.D. Thesis, Concordia University Canada.
[11]	Noui, O., Bouazara, M. and Richard, M.J. (2016) Study of Liquid Motion and Pressure Forces Applied on the Walls of Partially Filled Moving Tank. Volume 4A: Dynamics, Vibration, and Control, Phoenix, 11-17 November 2016. https://doi.org/10.1115/imece2016-65170
[12]	Maple (2016) Maplesoft, a Division of Waterloo Maple Inc., Waterloo, Ontario.
[13]	Fehlberg, E. (1969) Low-Order Classical Runge-Kutta Formulas with Step Size Control and Their Application to Some Heat Transfer Problems. NASA Technical Report 315.
[14]	Fortin, A. (2015) Analyse numérique pour ingénieurs. 5th Edition, Presses Internationales Polytechnique.
[15]	Toumi, M., Bouazara, M. and Richard, M.J. (2009) Impact of Liquid Sloshing on the Behaviour of Vehicles Carrying Liquid Cargo. European Journal of Mechanics-A/Solids, 28, 1026-1034. https://doi.org/10.1016/j.euromechsol.2009.04.004