1. Introduction
Nowadays the sloshing of fuel inside the tank of vehicles is a very important research issue. Many studies have been done in order to find a methodology able to predict the oscillation of the liquid inside the tank. As is well known, the sloshing phenomenon can be defined as the motion of the free surface of a liquid in a partially filled tank. The fluid motion is due to external forces and to the interaction with the gas phase. Consequently, the effects of sloshing on the containers are noise, inertial loads and structural vibrations and potentially is can also have effects on the stability of the vehicle [
1]. The relevance of the implications on the vehicle depends on several aspects. First of all, it depends on the type of motion, the amplitude and the frequency of excitation, but also on the geometrical aspects like the tank geometry, the fill levels and the liquid properties. The motion of the vehicle during the travelling is another important aspect that has to be taken into account in the sloshing analysis because with high acceleration, cornering and braking fluid inside the tank is forced to move. However, the prediction models need to include the interaction between the liquid and the gas by integrating multiphase applications.
In this paper the sloshing is analyzed looking in particular at fuel tanks, but this phenomenon is observed and studied for the all vehicle containers. The proposed numerical methodology can be applied to study sloshing in all the tanks inside a car. It has been observed that, in automotive applications, the sloshing does not affect the vehicle stability during travel; however, it may cause, as said, unwanted acoustic noise issues in the passenger compartment, structural vibrations, problems with low fluid (fuel, oil, water, etc.) level management and vapor emissions. Therefore, the individuation of a robust numerical model able to predict the free surface motion of a liquid is of great interest to the scientific community. For this reason, many researchers have investigated sloshing phenomena with modeling and experimental approaches.
Experimental analysis has been carried out by Pal et al. [
2] to study the sloshing in a cylindrical container with and without baffles. Aus der Wiesche [
3] also studied the sloshing in automotive fuel tanks experimentally, analyzing in particular the noise deriving from the phenomenon. He established a correlation between the simulated pressure fluctuations and the recorded slosh noise within the sloshing liquid. Abramson [
4] and Ibrahim [
5] have also analyzed the sloshing using both experimental and modeling techniques for different containers shapes.
Numerical models have been investigated by Liu and Lin [
6] and Popov et al. [
7]. In particular, Liu and Lin [
6] analyzed a nonlinear liquid sloshing with broken free surfaces while Popov et al. [
7] investigated the effect on the fluid motion of curvature and acceleration in rectangular containers. They observed that the most intense sloshing occurs when the fill level is approximately 30–60%. Ubbink [
8] developed an accurate Computational Fluid Dynamics (CFD) methodology able to predict the flow behavior of immiscible fluids. Fluids (water and air) are separated by a well-defined interface and are subject to an external solicitation to reply to a vehicle’s movement on a bend. The choice of the water instead of gasoline/diesel is strategic because of the safety, the simplicity of disposal and the unlimited availability at no cost. However, in the future the same Experimental Fluid Dynamics (EFD) and Computational Fluid Dynamics (CFD) analyses will be performed with fuels.
Even if the study concerns the analysis of the fluid dynamics of two immiscible fluids, it considers the strong cohesion forces between their molecules. In fact, when fluids are miscible, the study can be approached considering the surface tension (resistance to mixing). A negative value of the surface tension indicates no resistance to mixing, as underlined by Batchelor [
9]. As well known, the flow behavior of immiscible fluids depends on chemical reactions or phase changes due to high temperatures.
The flow of immiscible fluids can be classified into three groups based on the interfacial structures and topographical distributions of the phases, namely segregated flows, transitional or mixed flows and dispersed flows [
10]. The classification is explained by considering a closed tank-filled with a liquid and a gas. The container is forced to oscillate with low amplitudes and frequencies. The sloshing in this case falls in the first class where two phases remain separated with a single well-defined interface. If the frequency and amplitude are increased the phenomenon is classified as the second class. In this case there is a mixing or translating of flow because the waves become unstable and break [
8]. The mixing and the translating motion effect on the fluid and, as consequence, typically the interface between liquid and gas breaks up. Therefore, in this regime two small bubbles of gas are trapped in the liquid.
The last category is the third class. It defines the flow behavior when the fluid is shaken violently with a dispersion of the flows. The gas is suspended as small bubbles within the liquid. Rebouillat and Liksonov [
11] presented a review analyzing issues related to the solid–fluid interaction in tanks partially filled by liquid. Also, in this review the sloshing has been studied by means numerical approaches able to predict the phenomenon.
Another important aspect is the study of the fluid–structure interaction. This aspect has been investigated more recently using exclusively numerical methods. These numerical solutions concern the primarily finite-elements and finite-differences methods applied to Euler or Lagrangian fluid and solid domains. Li et al. [
12] solve the dynamic behavior of sloshing liquids in a moving tank by adopting the material point method (MPM). The MPM is a numerical scheme able to evaluate the impact pressure based on a contact algorithm over background grids. Authors validated the model by comparing simulation results with experimental data. Oxtobyet al. [
13] developed a semi-implicit volume of fluid free-surface modelling approach. The model demonstrated to be a good tool to study flow problems with a violent free-surface motion. A hybrid-unstructured edge-based vertex-center finite volume discretization with an entirely matrix free solution methodology was adopted. The high resolution artificial compressive (HiRAC) volume-of-fluid method captures the free surface in violent flow regimes. In this research numerical and experimental data were compared in a regime of violent sloshing.
The proposed research fits well in this scenario already described. In fact, it has the objective to introduce an approach to fully describe the sloshing phenomenon inside a fuel tank that combines experimental and numerical approaches. A preliminary model tank has been designed and then prototyped. This container has been realized with two faces of Plexiglas in order to visually capture the fluid oscillations during the tests. A test bench has been realized by the company Fiat Chrysler Automobiles (F.C.A.) [
14]. The “8-DOF bench” integrates a high-resolution camera with a resolution of 1280 × 720 pixels and a shot of 25 frames per second. An accurate three-dimensional CFD model has been built up starting from the designed tank geometry [
15]. The model has demonstrated to predict correctly the interfacial structure that is exposed to large deformations, merging or breakup. A mesh sensitivity analysis has been done allowing for an excellent accuracy (in comparison with experimental tests) and for computational efficiency. By the end, 3D pictures of both tests and numerical model have been compared for each water level and instant in the range [7 ÷ 12] s with a time step of 0.25 s. Then, pictures have been post-processed using Matlab
® (R2017 a of the company MathWorks Inc., Natick, MA, USA) in order to compare the center of gravity and the free surfaces. For this reason, the images have been reduced in two dimensions. The University of Naples “Federico II” and FCA Automobiles are still working in the numerical model to predict not only the sloshing phenomena but also the tank filling and the evaporation of the fuel into the tank itself. A similar study has also been carried out with the Flow-3D code [
16].
2. Experimental Set-Up
Before describing the developed numerical model to predict the sloshing phenomenon inside a fuel tank, the tank has been tested in collaboration with F.C.A.|Fiat Chrysler Automobiles (Pomigliano d’Arco, Naples, Italy). The bench is shown in
Figure 1. It has been designed by Moog Inc. (Elma, NY, USA) flowing the specifications provided by F.C.A. The bench is a hexapod with six independent actuators (with a total stroke of 600 mm) arranged in three triangles. Triangles are connected on the base and the platforms, thus allowing all six DOFs. The top platform is linked to a tilt table with two additional actuators (with a total stroke of 400 mm). This solution allows one to extend the pitch and roll envelope, thus the name of “8-DOF bench” [
15].
The hexapod design has been conceived based a small payload flight or driving simulators. There are ball screw types of actuators for both hexapod and tilt table. This solution allows a fast response, smooth run, low noise emission and low maintenance. The bench is equipped with an absolute encoder, to continuously monitor the position of the tank, and with a mechanical brake, integrated in the DC servo motor. The maximum excursion of the fuel tank allowed by the system is more than 400 mm in every direction, reaching peak accelerations of about 1 g (~10 m/s2). Pitch and roll rotation excursions, including the tilt table, reach over 50 deg. With the described test bench, it is possible to reproduce combined sway and surge accelerations of 1 g peak amplitude down to a frequency of about 0.7 Hz. This cuts out much of the low frequency dynamic that is very relevant for automotive applications. The tilt table allows one to simulate lower frequency accelerations using pitch and roll inclinations of the tank assembly. In this way, it is possible to reproduce experimentally the overall effect on fluid motion of very long accelerations, and changes like curves, speeding up and braking.
As clearly shown in
Figure 1b, the tank (highlighted in the black rectangle) has been located inside a rigid steel frame, representing installation conditions in the vehicle. The frame is then fixed on the top platform of the tilt table.
A Moog Replication software package (The integrated Test Suite, Moog Inc., Elma, NY, USA) manages the test bench. The software allows to define, view and edit the command signals, to operate and monitor the test, to operate signal identification and matching. The Moog Replication software has been specifically configured for the “8-DOF bench”. The bench, of course, includes many sensors. In fact, there are:
- -
XsensMTi G-700 ™ (Xsens, Enschede, Netherlands), measuring all six components of acceleration, plus magnetometer and GPS.
- -
Gyroscopes and tridimensional magnetometers make up the employed inertial platform. Gyroscopes that are also able to record rolling, pitching and yaw.
Even if the bench is able to obtain the rolling, pitching and yaw, in this paper the comparison between experimental and numerical data has been done only for the translation. Equations for the translation follow a harmonic law reported below, where the frequency is of 0.7 Hz and 0.5 Hz:
The research is still in progress therefore the effects on the sloshing due to rotation are currently under investigation.
Figure 1 shows the F.C.A. company test bench. In particular, in
Figure 1a a 3D drawing of the bench is shown, where the automotive tank under investigation has been located during tests. The tank is clearly shown in the black rectangle in
Figure 1b where, with the ellipse at the point (I the high-resolution camera located in front of the tank) has been highlighted. The camera has been rigidly fixed on the tilt table over the platform with an adjustable rod. Of course, it is necessary to fix the digital camera in the way to film any test performed at the bench and at the same time to maintain the position of the tank. The camera, as shown in
Figure 1c, is interposed between the tank and an opaque black panel in order to reduce reflection of the Plexiglas panel of the tank and, consequently, improve the video quality. The camera used has a resolution of 1280 × 720 pixels and shoots 25 frames per second. In
Figure 1c a zoomed view of the tank is shown. Point II in figure is the fill sensor, III is the air vent while IV is the exhaust valve of the tank.
Dynamic tests have been performed on a tank properly designed for the study made by high-density polyethylene and Plexiglas (see
Figure 1c) and it was fixed to the Moog test bench by means of two L-shaped steel profiles to ensure a high rigidity of the system during any test condition. Starting from the central section of a multilayer plastic tank of a F.C.A. commercial vehicle the tank polyethylene shape has been drawn. The designed and tested tank has exactly the same section of a production automotive tank thus two dimensions, width and height, are the same. The third dimension (depth) is lower (of 100 mm, about ¼ or ⅕ lower than the original tank). The depth has been modified introducing two faces in Plexiglas to capture the fluid oscillations during the tests. Therefore, the prototyped tank is 870 mm in length, 265 mm in height and 100 mm in depth and front and back transparent plates have 5 mm of thickness. It has been designed not only to observe the sloshing phenomenon but also to guarantee resistance, rigidity and waterproofness. The tank, during the test, is completely closed; therefore, no gas or liquid can enter or exit from the system.
As already mentioned, the working fluid is simply water colored by a dark blue food colorant. The choice of the water as a working fluid has been done to easily monitor the sloshing phenomenon. The choice is also strategic because of the safety, the simplicity of disposal and the unlimited availability at no cost. Therefore, as a multiphase application, the properties of both phases are listed in
Table 1.
Tests have been done for three different free surface levels H. The tank has been filled at:
H = 165 mm;
H = 101 mm;
H = 65 mm.
In
Figure 2 the three levels are shown. Numerical simulations have been done for each fluid level in the tank and results have been compared with tests with a high agreement.
In order to compare tests with the numerical results images post-process is necessary. All the pictures captured by the camera highlighted in
Figure 1b have been analyzed using Matlab
®. The post-processing phase has been done in order to obtain clear images usable to obtain important information like the center of gravity and the free surface.
4. Numerical Model Description
From the numerical point of view, the adopted methodology consists of a three-dimensional CDF approach developed using the commercial code PumpLinx
® (V4.6.3, Simerics Inc., Bellevue, WA, USA) [
18]. The numerical model has been built up including new capabilities to correctly predict the flow behavior of two immiscible fluids, separated by a well-defined interface. Thus, the model consists of two immiscible fluids in movement separated by an interface. In particular, in this section, all equations governing the phenomena are presented separately.
The model mathematically describes the fluid flow solving the equation of conservation of mass and momentum. Furthermore, to describe correctly the fluid sloshing inside the tank, additional properties, such as viscosity, surface tension and compressibility, have been implemented as input for the simulations.
Considering a flow quantity
and a control volume reported in
Figure 5 [
15,
17,
18,
19], the Gauss’s theorem can be applied as reported into the equation:
where
indicates the time,
is the convective flux over the boundary due to the motion of the fluid,
is the fluid velocity,
is the flux over the boundary conditions due to diffusion,
is the volume of the control volume,
is the internal source,
the source at the boundary. In the analyzed case the
and the
(the internal source and the source at the boundary) are the gravity and the displacement law variation (function of the frequency) acting on the fluid volume of the tank.
If the flow quantity is the mass, the transport equation can be defined as in (3) assuming no sources such as chemical reactions or phase changes in the system. Equation (3) is also known as the continuity condition. It states that, for the incompressible case, the mass of a fluid in a closed domain can only be changed by flow across the boundaries:
Before to describe the built model, it is important to deeply describe the methodology adopted presenting the equations solved during simulation. Particular attention must be reserved to the definition of the model for the turbulence. In this application, good results have been achieved using the Re-Normalization Group (RNG)
k −
ε model for the turbulence. In fact, for the specific problem it has proven to be really accurate. The standard k − ε model, used for the simulations presented in this paper, follows the two equations reported below [
20,
21,
22,
23,
24,
25,
26,
27]:
with
= 1.42,
= 1.68,
= 1,
= 1.3,
is the fluid velocity and
is surface motion velocity while
and
are the turbulent kinetic energy and the turbulent kinetic energy dissipation rate Prandtl numbers. As said in
Section 2 of the paper, the designed tank has been tested for three different liquid levels. The container is filled by water and air therefore simulations need to solve a multiphase flow. This type of flows is better solved with explicit numerical schemes, especially if an accurate prediction of the free surface is required. To achieve a stable solution, these schemes require the Courant-Friedrichs-Lewy condition to be satisfied. This condition is defined by Equation (6):
The magnitude is the Courant Number and it represents the distance that the interface between two components would travel, across a mesh cell, during one time-step. The Courant Number controls the sub time-steps size in the calculation of the transport of the components for the multiphase module. In the code PumpLinx®, the Courant Number is set to a default value of 0.34, with a lower number providing more accuracy but, on the other hand, longer simulation time.
The calculation of a sub time-step, Δtsub, is performed only for the cells in the domain that have a mixture of components. Then the smallest of those calculations is used as the sub time-step, Δtsub, for updating the transport of the components.
As mentioned above, an explicit method preserves the shape of the interface between components and it is more accurate than the implicit one. A Blending Factor (
BF) is used in conjunction with the higher order interpolation schemes (e.g., Central and 2nd Order Upwind) to help stabilizing the convergence by including the 1st order Upwind scheme using the formula reported below [
17,
25]:
Higher values of the Blending Factor make the solution more stable, but potentially at the expense of accuracy. Therefore, a typical value falls in the range from [0 ÷ 0.5]. In this study, it is zero.
First of all, before running the simulation, the gravity force (in this case, the direction is along the
z-axis) has been included in the model. In thermodynamics, surface tension is dealt with on a macroscopic level based on the statistical behavior of interfaces [
3]. Sabersky and Acosta [
20] noticed that it is better to consider the microscopic level to understand the mechanism of this phenomenon. In fact, many other molecules surround a liquid molecule and away from the interface.
The surface tension coefficient σ is defined by Yuan [
21] and by Massey [
22]. It is the amount of work necessary to create a unit area of free surface. It always exists for any pair of fluids, and its magnitude is determined by the nature of the fluids. For immiscible fluids, the value is always positive and for miscible fluids such as water and alcohol, it is negative [
9].
An important aspect of surface tension is that it creates a pressure jump Δp across a curved surface. The following equation correlates the pressure jump across a curvature surface with the surface tensor:
where
and
are the principal radii of curvature of the surface,
is the mean curvature,
is the surface tension coefficient and
is the pressure (the larger pressure
is on the concave side of the curved surface).
In order to achieve the final goal of the study, some assumptions have been fixed. First of all, the flow inside the tank is incompressible, Newtonian and homogenous with a density and viscosity μ constant. At last, the temperature dependence has been not considered. Furthermore, the is considered equal to the tank velocity on the tank walls; therefore, there are non-slip conditions on the walls. As said, the analysis has been done considering two immiscible fluids subject to external forces. The flow of immiscible fluids can be classified into three classes; these classes have been already described in the introduction section. The first class concerns the sloshing in this case where the two phases remain separated with a single well-defined interface with a motion with low frequency and amplitude. Thus, this study falls in the first one because oscillations are gently with a low amplitude and frequency. The same modeling approach can be used also in the other two classes of sloshing already defined in the introduction (classes two and three).
The experimental setup has been replied using a three-dimensional CFD numerical approach [
15,
21,
22,
23,
24,
25,
26]. Therefore, starting from the real CAD geometry of the prototype tank, the fluid volume has been extracted and then meshed. The fluid volume of the tank is shown in
Figure 6b. This volume has been extracted from the real fuel tank in
Figure 6a.
As a study of two immiscible fluids, before importing the geometry in the simulation code for the grid generation, it is important to define in advance the initial (static) free surface for each simulation. The liquid level has been already indicated with the letter H. Experimental tests have been done for three H levels:
- -
H = 65 mm;
- -
H = 101 mm;
- -
H = 165 mm.
In
Figure 7 the free surface in the tank drawing at H = 101 mm is presented. Each dimension has been taken with respect to the coordinate (0, 0, 0), center of the drawing. The bold red line in
Figure 7 represents the considered free surface of one of the three analyzed cases.
The extracted fluid volume in
Figure 6b has been meshed using a body-fitted binary tree approach. This grid type is accurate and efficient because the parent-child tree architecture allows for an expandable data structure with reduced memory storage. In this architecture, the binary refinement is optimal for transitioning between different length scales and resolutions within the model. The majority of cells are cubes, which is the optimum cell type in terms of orthogonallity, aspect ratio, and skewness thereby reducing the influence of numerical errors and improving speed and accuracy. It is important to underline that since the grid is created from a volume, it can tolerate inaccurate CAD surfaces with small gaps and overlaps.
The entire grid of the tank fluid volume is shown in the
Figure 8 with also a zoom view of the grid in a particular area around edges. With the described grid generation approach, it has been demonstrated [
15] that it is possible to reach good accuracy and low computational time also by increasing the grid density in the boundary layer and on the surface. However, in case of multiphase simulation, it is better to have a mesh as uniform as possible. The zoomed view in
Figure 8 shows that in this particular application the generated mesh almost is uniform, especially around edges. To better simulate the sloshing phenomenon, the grid has been subdivided and cut in the regions of high curvature and small details.
Due to the complex analysis it has been selected the automatic refinement algorithm of PumpLinx® for the mesh generation. With this approach, it is possible to assign a maximum and a minimum cells size (and, of course, the cell size at the walls) and the mesh generator automatically adapts the refinement strategy and the levels of refinement as a function of the prescribed values and the geometry provided. This is a proprietary methodology developed by Simerics Inc. (Simerics Inc., Bellevue, WA, USA), the code developer.
A maximum cell size parameter of 3 mm has been chosen, whereas the minimum cell size parameter has been fixed at 0.5 mm. The maximum cell size is a parameter used to build the grid which controls the size of the cells while the minimum cell size is a parameter used to limit minimum dimension of cells. This means that, no cell in the volume can have a cell side smaller than the minimum cell size. At last, the cell size on surfaces also needs to be fixed. This parameter is used to control the size of the cells for all surfaces of a mesh volume.
As a dynamic phenomenon, simulations are transient and have been run on a multi-core Windows® 64-bit PC (CPU with two Intel(R) Xeon(R) E5-2640 v2 processors operating at @2.00 GHz, equipped with 192 Gb RAM). Before showing the numerical results, it is important to underline that, as normally done, a mesh sensitivity analysis has been performed.
Figure 9 shows three of the generated meshes. The first grid called “Mesh 1” consists of 2 million cells. The maximum cell size and the cell size on surfaces have been set to 2 mm while the minimum cell size is of 0.5 mm. The grid called “Mesh 2” consists of 1 million cells with a maximum cell size and the cell size on surfaces set to 3 mm and with the same value of the “Mesh 1” for the minimum cell size. The grid called “Mesh 3” has only 50.000 cells. Also in this case the minimum cell size is of 0.5 mm while the other two parameters are higher than in other cases.
Simulations have been run for the three meshes with the same boundary conditions and results with 1 million and 2 million cells have demonstrated the same level of agreement with the experimental data. However, the model with 2 M cells has higher computational time than the one with 1 M cells. At the same way, even if the computational time of the mesh with only 50.000 cells is lower than for others, results are too far from the experimental data. For these reasons, in this study the mesh size of one million cells has been employed; it shows excellent accuracy and low computational time.
Having built the mesh of the model, it is important, before running simulations, to set all the input parameters. In this case, the fluid level is one of them. It can be set writing the Equations (11) and (12), reported in the following paragraphs, in the expression editor panel of the code. It automatically splits the fluid volume of the tank, in
Figure 7, into two phase
ao and
wo where
ao is the air and
wo is the water. Equations (11) and (12), therefore, define the percentage of filling of the tank. In this research, three levels of filling the tank have been studies (called Test 1, 2 and 3) better described in the next section (
Section 5). The free surface of the water, for each filling level, evolves during simulation as function of the time, as done during the tests. Both, tests and simulations, last 12 s. In particular, the air level is defined as in the Equation (9) instead the water level in the Equation (10):
Equations (9) and (10) are referred to the free surface level of 101 mm already shown in
Figure 7. There is no special mesh treatment at the interface between the air and the liquid. The mesh is uniform across the flow field because the liquid interfaces in this simulation can potentially reach everywhere in the domain and change very quickly through time. A surface tension of 0.0073 N/m has been applied on the water. As known, the surface tension depends on the combined properties of the fluids on either side of the interface, but for the current release, the surface tension is limited to associating it to only one component, even for multiple fluids.
As done during tests, an assigned displacement along the direction y for the entire tank has been assigned as input for the model. The displacement law variation has been already shown in the paper
Section 2 (Equation (1)). The displacement law variation is, therefore, function of the frequency (0.5 Hz and 0.7 Hz) and time. In order to replay the experimental data, simulations have been run for 17 s with 340 steps all saved for comparing numerical results with tests in the validation phase.