1 Introduction

Sloshing^[1] is a common phenomenon in nature and industry that usually occurs in containers partially filled with liquid. It involves two or more mutually miscible fluids within a limited space of movement that exhibit highly random and nonlinear characteristics. Sloshing is not only a common physical phenomenon but also one of the key issues in various engineering research fields.

For instance, in shipbuilding, sloshing^[2] is the key phenomenon in ship design. The interior of a ship is equipped with a fuel oil tank, an oil tank, a fresh water tank, and a slop tank liquid cabin. If these hanging liquids in the tanks shake too strongly, the ship may overturn^[3]. By contrast, the problem of liquid sloshing in a storage tank in liquefied natural gas (LNG) and liquefied petroleum gas (LPG) carriers comprises a larger proportion in transportation safety. In the field of aviation, sloshing mainly occurs in aircraft fuel tanks. With aircrafts taking off, landing, or changing states, fuel will always shake. The force and moment produced by sloshing fuel will affect a plane's flight posture^[4].

In a word, sloshing is not only a very common problem in the field of engineering but also the core issue of many key technologies. Solving the problem of sloshing in the process of artificial control will be an important development in these areas.

Although sloshing is very common, its high degree of randomness and nonlinearity yield a theory of sloshing that is far from perfect. In recent years, researchers have focused on studies of fluid-structure interaction, the response of structure, and the elastic tank^[5-6], among other aspects, all of which are fundamental to theoretical research surrounding sloshing. In the study of sloshing theory, the difficulty related to the numerical method is a barrier to most studies. Since the 1870s, various numerical methods have been proposed. With the rapid development of computer technology, the numerical method has gradually become an important means of studying sloshing. The primary methods of solving the governing equation of discrete states are the finite element method, finite difference method, boundary element method, finite volume method, and spectral method^[7]. The key issue in solving the sloshing problem is related to the description of the free liquid surface. Common description methods are the marker and cell (MAC) method, the volume of fluid (VOF) method, and the Level-set method^[8].

This paper aims to accurately simulate and analyze the sloshing problem, to solve some specific sloshing processes with the application of various numerical methods, and to carry out a parametric analysis of the sloshing process. On the basis of the simulation, we use a container plate shape design to reduce the force of side wall hanging based on the optimal control theory and topological optimization.

2 Governing equations and mathematical model

The governing equations are three-dimensional (3D) incompressible Navier-Stokes equations,

(1)

(2)

where v, F_b, p, and μ are the fluid velocity vector, the force vector of unit fluid mass, the pressure, and the coefficient of viscosity, respectively.

The basic idea of using the VOF method to solve fluid motion is to solve the Navier-Stokes equation through application of the finite difference method. Interface tracking is a key step in the VOF method. Among the various ways to carry out interface tracking, the PLIC (interface refactoring method), minimum variance, and ELVIRA^[9] methods can be mentioned. These methods can reflect a relatively realistic interface motion of the fluid. This paper chooses the PLIC method because of its relative technological maturity and because it is the preferred method at present.

The VOF method was proposed by Hirt and Nichols^[10-11] in 1975. Its basic idea is to represent the fluid type by a function, called "volume function" F. F and (1-F), respectively, represent the volume fraction of gas and liquid in the computing area. For a single computing element, the volume fraction takes one of the following three values:

(ⅰ) F=0: on behalf of all for gas in the unit;

(ⅱ) F=1: on behalf of all for liquid in the unit;

(ⅲ) F between 0 and 1: on behalf of the gas-liquid interface present in the unit.

The free surface-sloshing research belongs to the third category mentioned above. The normal direction of the free surface can be determined by the gradient of F after computing the location of the free surface, which will enable us to determine its approximate shape.

The moving boundary method is applied in the process of sloshing calculation. The incentive of sloshing is produced by the boundary, and the calculated result is observed on boundaries with a moving coordinate system.

2.1 Interface tracking

In this paper, the PLIC method is applied to track the interface with a simplified three-step process: the estimation of interface's normal vector, the construction of a two-dimensional surface, and the spread of the interface within the flow. The normal direction of the free surface can be determined by the gradient of F. The free surface can be expressed in the form of a section with the conditions of fluid volume fraction. This section changes over time. To update the new section, one must compute the mass, volume, and momentum flux of each fluid flowing into the adjacent surface. The numerical simulation of flow can be found using a series of calculations. The process of flow with interface can be described either with the Lagrangian or Eulerian methods. The Lagrangian method is used in this paper, and the location of a point in future time can be found from the location of the point in present time.

2.2 The momentum equation and the discrete surface tension

The VOF/PLIC method is used to calculate the volume function. The velocity field is calculated based on the Navier-Stokes equation. The steps followed in the calculation of the velocity field are as follows: first, calculate the temporary speed in two steps; second, solve the pressure by Poisson's equation; and finally, update the velocity field by the pressure value and obtain the velocity field of the next step.

The central difference method is used to discretize the spatial derivatives. In this paper, the grids in the calculation are averaged as physical quantity, such as density and pressure, and are transferred to the interface in the grid. The method is explained as follows:

(3)

However, as some variables such as pressure may change rapidly, depending on specific circumstances, some special methods are needed.

2.3 The moving boundary method

In this paper, the computational area is a swaying rectangle tank in which the motion of fluid is driven by the swaying wall of the tank. In other words, introducing the speed of the same direction and value to the six walls of the tank (i=i_min, i=i_max, j=j_min, j=j_max, k=k_min, and k=k_max) causes the fluid in the tank to slosh.

The research program has gone through the test of the laminar flow and the single disturbance simulation, which proved to be effective. Therefore, the method adopted in this article to simulate the sloshing process is reasonable.

3 Effects on the sloshing process after placing a board in the field

The simplest and most common method to reduce sloshing in the field of engineering^[12] is to insert a board in the fluid. To study this problem, an impermeable board is placed in the liquid storage box, and the ensuing changes in the flow field are analyzed.

In this paper, the physical quantities are dimensionless, and their unit is 1. Sloshing occurs in a relatively flat cuboid liquid storage tank, which oscillates in a periodic motion with a speed of u(t)=0.7cos (50t), driving the fluid in the tank to slosh. An internal boundary condition is added to the sloshing area, namely, a partition set into the center of the flow field with the movement of the outer boundary. The partition is as wide as the computational domain, and higher than the fluid interfaces. The role of the partition is to eliminate the relative speed of fluid flow passing through the wall. The concrete numerical values of the flow parameter are shown in Table 1.

3.1 The kinematic characteristics of flow under a baffle block

When blocked on either side, the fluid sloshes as if exiting in two computational domains at the same time. In Fig. 1, a downward wave is formed first on the right side of the computational domain. On the left side of the computational domain, there is a phenomenon of water driving up, as sloshing occurs on the right side with the subsequent formation of a downward wave on the left side. An upward velocity component appears on the right side of the computational domain.

The simulation result and reality are completely different. In theory, when it happens, sloshing should occur in synchronization on both sides and become identical. We now ask a question: Why does it show an up-down movement in the opposite direction in the left-right domain? It is discovered that the internal boundary conditions of the partition problem outlined in this chapter are different from those of the outer boundary, where the fluid speed through the partition becomes zero in the y- and z-directions and equal to the border in the x-direction. Therefore, we implement artificial changes only to the speed factor in lieu of entirely separating the left-right fluid domains. To verify the accuracy of our calculations, the direction of the sloshing movement is reversed while still appearing out of sync, and it creates a downward wave on the left side of the domain.

Figure 2 depicts the velocity vector change in the process of sloshing for the first cycle results. When sloshing begins, the general trend of the fluid is to hang over the direction of the drive and move to the right. However, on the boundaries of the tank and the baffle plate, a large velocity change is observed. On the left side of the border, most fluid moves downward, and on the right side, most fluid moves upward. That is why the two domains have different surface levels. As the sloshing value on the right side reaches its highest point, the fluid speed decreases (see Fig. 2(b)), and the wave amplitude produced by sloshing is the highest at this moment. Subsequently, the drive sloshing moves to the left, and from then on, so does the fluid. The phenomenon repeats via the above process in reverse.

3.2 The dynamic characteristics of flow under the baffle block

Although a partition board is added in the middle of the fluid domain, fluid movement still has a certain impact on the side walls, as seen in the blue color in Fig. 1. The fluid pressure impact on the side walls is shown in Fig. 3.

To show the pressure changes in flow fields more clearly, the pressure contour is plotted along the flow direction of the vertical section (y=0.5y_max). As seen on the plot, the change in the pressure area is directly proportional to the size of the wave, as shown in Figs. 4(a) and 4(b). In addition, the figure shows cyclic stress variation in the process. In the first half cycle, the low- and high pressure areas increase on the right and left sides, respectively. In the second half of the cycle, both low and high pressure areas on the left and right sides increase, respectively.

During sloshing, an inserted partition board in the middle of the flow field divides the sloshing tank into two flow chambers. At the same time, the partition board can change the sloshing's kinematics and dynamics characteristics. This is clearly observed from the free surface figures, showing that the waves are formed in the side wall and baffle of the internal boundary, which are used to be formed only on the side wall. In addition, compared with the absence of a baffle, the impact pressures on side wall decrease with time.

4 Topology optimization in suppressing sloshing 4.1 The unsteady intelligent optimization and topology optimization methods

In this paper, the solutions of an ordinary differential equation^[13-14] are used for optimization. The optimization goal is to select the target parameters x₁, x₂, x₃, ..., x_n, which make the function value of ω as small as possible. The specific qualifications are as follows:

(4)

In this formula, m₁, m₂, and m₃ are integers that satisfy the relationship 0≤m₁≤ m₂≤ m₃, f_i is a real number, and g_i is a continuous equation.

The method can be used to solve more than one type of problem. For example, when m₁=m₂=m₃, it can solve unconstrained numerical fitting problems; when m₁=0, m₂=1, and m₃>1, it can solve unconstrained minimization problems; when m₁=m₂=0, it can be used to solve the feasibility problem. This kind of optimization method has been applied to many studies, such as those on the topology optimization of biomimetic fishtail shape and optimization of airfoil, and satisfactory results are obtained.

The optimal control problem is like a dynamic optimization problem over time^[15]. The state equations and the initial state of the controlled system are

(5)

(6)

The target is defined as

(7)

For an admissible control u(t)∈ U, t∈ [t₀, t_f], make the system begin from the initial state, and make the following minimum performance indicators:

(8)

The above equation is the optimal control problem, whose result u^* (t) is known as the optimal control or extreme value control. The performance indicators from calculation J^* are known as the optimal performance index.

4.2 The partition topology optimization in the reducing sloshing

In this paper, the optimization method is used in the sloshing analysis of fluid mechanics, combined with the optimal control theory, to optimize the permeable function with a partition board. The permeable function parameter is indicated as Ω. The impact on the side wall of the fluid is regarded as the partition optimization of the objective function. The purpose of our optimization is to find a way to minimize the objective function via adjustment of the grid permeable function parameter around the partition. Thus, the objective function is calculated by the following method:

(9)

In this formula, p is the pressure of the flow field, which is obtained using sloshing flow field calculation.

The partition is divided into 6× 6 grids, and at each grid, a permeable rate value is given. The permeable rate is optimized in a limited amount of feasible region.

The optimal control process of sloshing is as follows: the initial value is Ω = 0. 5, and the fluid is in a semi-permeable state. First, the fluid is left to freely slosh for a period of time (t_start); after the initial motivation effect disappears (t>t_start), every step of the time interval is optimized. To simplify the optimization process, the values of the permeable rate Ω are limited in the range from 0.3 to 0.7. Hence, for the current optimization process, the permeability rate Ω = 0.3 is considered a complete waterproof partition, and as the permeability rate assumes the value of Ω = 0.7, the partition is considered to be completely permeable. Before the next calculation step, the optimization permeability rate is substituted for the original one, and step by step optimization continues until the optimization goal is achieved.

The partition shapes shown in Fig. 5 are representative geometries in the whole optimization control. The third sloshing image in the cycle is chosen as an optimum geometry, as it can avoid unstable effects when sloshing occurs. When the fluid through the baffle rushes to the side wall, the topology optimization of the baffle always finds the shape in Fig. 5(a), namely, fluid is allowed to travel freely in most areas to reduce the impact force as much as possible. In the middle and bottom of the baffle, there are two horizontal bar block areas. Figure 5(b) shows the optimal baffle shapes at late impact on the side walls. At this moment, the shape of the baffle is very chaotic, but there exist some similarities among the shapes, as seen in Fig. 5(c). Most of the area is water-proof, and there will be some semi-permeable areas, as seen in Fig. 5(c). Figure 5(b) can be considered as the transition between two main baffle types. Figure 5(c) shows the optimization shape of the baffle when the fluid leaves the side wall. At this moment, fluid is expelled from the four oblique areas by the baffle. Underneath, the vertical distribution is hollow. Figure 5(d) shows the optimal partition shapes at the transition stage that simply leave from the baffle and are on the way to impact the baffle again. The image is very similar to the results in Fig. 5(a). What makes it different from Fig. 5(a), however, is the presence of a nearly horizontal block of water area, as seen in Fig. 5(a).

Although we only choose the special shape of the partition obtained in a cycle, in fact, the partition shapes change cyclically. Sloshing is a cyclical process, and the partition topology optimization results also change in a cycle.

Next, let us examine what will happen to the value of the resultant force of the impact stress on the side wall under this optimization condition. As seen from Fig. 6, after optimizing the partition, the value of the objective function (the resultant impact force) has obviously been reduced most of the time. Figure 6 proves that partition topology optimization does indeed work.

Once the law of change in the optimal control of the partition shapes is understood, the results obtained are expected to be applied in practice. Since the partition shapes in the actual application could not have been changed, we want to find a way to apply the results to the sloshing reduction of the fixed partition.

In this paper, the value Ω of each point on the partition grid, which occurs in the second sloshing cycle, is averaged as the fixed rate of permeability to the grid. The results are shown in Table 2.

These data are used to obtain the permeability rate chart, as shown in Fig. 7.

According to the area shown in the chart, the partition shape can be drawn roughly, as shown in Fig. 8. The black and white parts of the figure represent the impermeable solid wall and the permeable parts of the partition, respectively.

To test whether one can reduce the impact force on the side wall, sloshing experiments were conducted again by adopting a new permeability rate array. Comparison between the resultant impact force values from experiments and prior optimization values is depicted in Fig. 9. As the diagram shows, the optimized resultant impact force values are significantly reduced. This confirms feasibility of the proposed optimization method.

5 Conclusions

To study the sloshing problem within a cuboid container, the VOF method is used to simulate the sloshing process. The sloshing process is simulated when there is a board settled in the middle of the flow field. Finally, in order to reduce the pressure produced by the sloshing fluid on one side wall, topology optimization and optimal control are used in the design of the board shape.

The results are as follows: the waves are regular during the sloshing, while the surface exhibits chaotic behavior over time. However, a state of dynamic equilibrium can be achieved after adjustment. The conclusion about the reflected fluid by vibration is as follows: the wave size is directly proportional to the size of the wave amplitude, and inversely proportional to the frequency. If we place a board to divide the water in the middle of the flow field, the impact on one side of the wall may increase, while also preserving the sloshing form. The impact on the side wall is taken as the sloshing measure, and Ω denotes the permeability. As for topology optimization, Ω is defined as the independent variable, and the resultant impact on one side wall is defined as the target. After optimization, the results are found to have periodicity and can be separated into several types. Finally, we make an average shape of the board through a general consideration and verify its effect on suppressing sloshing.

In suppressing sloshing, topology optimization is used for the first time and provides a new approach to research in this field of study.