The Finite Element Numerical Investigation of Free Surface Newtonian and Non-Newtonian Fluid Flows in the Rectangular Tanks

In this paper, we develop a new computational framework to investigate the sloshing free surface flow of Newtonian and non-Newtonian fluids in the rectangular tanks. We simulate the flow via a two-phase model and employ the fixed unstructured mesh in the computation to avoid the mesh distortion and reconstruction. As for the solution of Navier–Stokes equation, we utilize the SUPG finite element method based on the splitting scheme. The same order interpolation functions are then used for velocity and pressure. Moreover, the moving interface is captured via the concise level set method. We take advantage of the implicit discontinuous Galerkin method to handle the solution of level set and its reinitialization equations. A mass correction technique is also added to ensure the mass conservation property. The dam break-free surface flow is simulated firstly to demonstrate the validity of our mathematical model. In addition, the sloshing Newtonian fluid in the tank with flat and rough bottoms is considered to illustrate the feasibility and robustness of our computational scheme. Finally, the development of free surface for non-Newtonian fluid is also studied in the two tanks, and the influence of power-law index on the sloshing fluid flow is analyzed.


Introduction
e sloshing free surface flow in a rectangular tank occurs in many areas, such as oil exploitation, ocean engineering, and hydraulic engineering, and it often involves Newtonian and non-Newtonian fluids. e free surface will move under the gravity, and the profile of interface front changes with time evolution. In recent years, the numerical investigation of this problem has become a research hotspot and received a great deal of attention [1][2][3][4][5].
Some authors have used a mesh-free method to study the sloshing free surface fluid flow. In 2012, Shao et al. [6] have proposed an improved smoothed particle hydrodynamics (SPH) to simulate the liquid sloshing. In their calculation, the Reynolds averaged turbulence model is coupled with the SPH method. ey also modified the scheme to achieve smoother pressure field and utilized additional algorithm to treat the solid wall boundary condition. In 2014, Cao et al. [7] found a more appropriate kernel function in the SPH simulation. ey also considered dummy particles and a new way to handle the moving boundary. In 2020, Fu et al. [8] presented a semi-Lagrangian meshless framework to study the sloshing phenomenon. e localized radial basis function collection method is employed to get velocity, and then, the free surface elevation is computed via second-order explicit Runge-Kutta method. e semi-Lagrangian algorithm constrains the lateral motion of inner collocation points. Generally speaking, it is difficult to handle the boundary condition and the irregular computational domain for the mesh-free method. In addition, the research about sloshing non-Newtonian fluid flow based on the particle method is very few.
In the past several decades, the standard particle finite element method (PFEM) [9][10][11] is proposed and has been used for the simulation of free surface fluid flow. e general way is to solve the Navier-Stokes equation via FEM on grid. Since the boundary is moving, it is tracked by a purely Lagrangian method according to the velocity. In 2014, Oñate et al. [12] developed a Lagrangian elemental FEM to study the free surface flow. e equations are integrated over the elements, which is the same as in the classical FEM. In each remeshing step, the original elements will be discarded and a new triangulation will be generated. In 2020, Franci et al. [13] have proposed another Lagrangian nodal integration method to simulate the free surface flow. e definitions of variables are all at nodes, and the integrals are performed over nodal patches. In their approach, the results will be less affected by the remeshing operations. In these methods, when the mesh is distorted, the reconstruction operation is necessary to guarantee the mesh quality. However, this process is difficult and time-consuming, especially for the complex computational geometry.
In this paper, we aim to research the sloshing fluid flow via a two-phase flow model with the fixed computational mesh. We employ the splitting scheme to decouple the unified Navier-Stokes equations. After that, the elliptic subequations are solved via standard FEM for the high efficiency.
e streamline upwind/Petrov-Galerkin (SUPG) method [14,15] is utilized to solve the hyperbolic subequation for the stability. e same order interpolation functions are used for the velocity and pressure. As for the capturing of free surface, the simple and efficient level set method [16,17] is utilized. We employ the implicit discontinuous Galerkin method [18,19] to handle the level set and its reinitialization equations for the accuracy and stability. With a view to the main drawback of the level set method [20], a correction technique [21,22] is also added to preserve the mass conservation property. In the numerical modelling, we consider both the sloshing Newtonian and non-Newtonian fluid flows in the rectangular tank with flat and rough bottoms. According to our computational algorithm, there is no need for the mesh movement or reconstruction and it is easy to deal with boundary conditions in the irregular domain. e paper is organized as follows: in the next section, we introduce the two-phase flow model and the power-law model of the non-Newtonian fluid; in Section 3, we describe the splitting scheme and temporal and spatial discretizations of our computational algorithm in details; in Section 4, the sloshing Newtonian and non-Newtonian fluid flows in the tank with flat and rough bottoms are all studied to illustrate the validity, feasibility, and robustness of the computational scheme. We also compare the development of free surface for different cases and analyze the influence of power-law index on it; finally, in Section 5, we give the summary of the concluding remarks.

Two-Phase Flow Mathematical Model
As for the free surface flow of sloshing liquid in a rectangular tank, it mainly consists of liquid and air. erefore, it is natural to regard this problem as a two-phase flow. In our mathematical model, the level set equation is employed to capture the free interface, and the governing equations are written in a unified form. e details are shown as follows.

e Level Set Equation.
e zero contour of level set function ϕ represents the free interface Γ in the level set method. According to the definition, the absolute value of ϕ denotes the distance from the point to the interface |ϕ(x)| � min In the liquid area, the value of level set function ϕ is less than zero and it is bigger than zero in the air region. e level set equation is written as follows [23]: As for the incompressible flow, the above equation is rewritten in the following conservation form: In order to keep the signed distance function property of the level set function, the reinitialization procedure is essential.
e conservative form of level set reinitialization equation with initial condition is given in the following equation [23]: . φ 0 is the initial level set function needed to be reinitialized, and τ represents the artificial time. ε equals to 1.2 h, and h is the maximum element diameter.

A Unified Form of Governing Equation.
As for the incompressible liquid and gas, the governing equations for each fluid are both Navier-Stokes equation. e main distinctions are the density and viscosity of these two fluids. For the convenience of computation, we intend to unify the governing equations as one formulation in the whole computational region via the level set function. e unified momentum equation and the continuity equation are written as where p is the pressure. u � (u, v), and u and v are the components of velocity in the horizontal and vertical directions, respectively. In addition, f � (0, −g), g means the gravitational acceleration, and there is no external force in the horizontal direction. In the whole region, the formulations of density and viscosity based on the level set function are given in the following equations [22,24]: where ρ l and ρ g represent the density of liquid and gas and μ l and μ g denote the viscosity of liquid and gas. Τhe definition of H(ϕ), the smoothed Heaviside function [25], is given in the following equation: where ε is selected as 1.2 h in the calculation.

Power-Law Model.
According to the rheological theory, the viscosity of non-Newtonian fluid will be affected by the velocity gradient. We take the power-law model [26,27] to describe a nonlinear relation between shear stress and the rate of deformation. e strain rate tensor is denoted as d � (∇u + (∇u) T )/2, and its magnitude is c � ������ � 2(d: d). us, in terms of power-law model, the viscosity is where n represents the power-law index. When n � 1, n < 1, and n > 1, the model represents the Newtonian fluid, the pseudoplastic fluid, and the dilatant fluid, respectively.

The Numerical Algorithm
As for the modelling of sloshing free surface flow, it mainly involves the numerical solution of governing equations for the flow field and the capture of moving interface. We utilize the SUPG method based on a splitting scheme to handle the unified Navier-Stokes equations. In addition, the level set and its reinitialization equations are solved via implicit discontinuous Galerkin method.

Temporal Discretization and Splitting Scheme.
As for the temporal discretization of unified Navier-Stokes equations, the simplest Euler algorithm is used to discretize the term with time derivative. e nonlinear convective term is treated explicitly. Moreover, the pressure and viscous terms are both managed implicitly. And then, it gets the discretized momentum equation in time: According to the splitting scheme [28,29], we can obtain three subequations. We first neglect the terms about pressure and viscous in equation (8) and introduce an intermediate velocity u to replace u n+1 . It is then able to get the hyperbolic subequation: After that, we bring in another intermediate velocity and achieve an equation about pressure: We take divergence on both sides of equation (10) and apply the velocity divergence free condition for the intermediate velocity u to get the Poisson equation of pressure: And then, the intermediate velocity u is able to calculate based on the following equation: In the last stage, it is able to acquire the Helmholtz equation about velocity: When we plus the subequations of (9), (10), and (13) together, we are able to recover equation (8).

Preliminary for the Spatial Discretization.
Before the following introduction of spatial discretizations, we need to do some preparatory works. Let us assume that the computational domain Ω is divided into N nonoverlapping small conformal triangles and h represents the maximum element diameter. us, the area Ω is able to be approximated by Ω h � ∪ N i�1 Ω k . e continuous and discontinuous spaces are defined as where P k (Ω i ) represents the polynomials of degree less than or equal to k on Ω i . C 2 k (Ω h ) and D 2 k (Ω h ) are the symbols of vector version.
In the discontinuous space, we take the "+" superscript and "−" superscript to represent the information in the exterior and the interior of an element, respectively. e average operator [30] is defined as It is effective for both scalar and vector. Furthermore, the jump operators [30] for the scalar and vector are also given as According to the average and jump operators, the popular Lax-Friedrichs numerical flux [30] is able to be written as follows: Mathematical Problems in Engineering 3.3. Spatial Discretizations of the Subequations. After the temporal discretization and splitting scheme, we receive several semidiscrete subequations. e first subequation of equation (9) belongs to the hyperbolic equation, and we employ the SUPG method to discretize it in space. Let us find u h in the space of C 2 k (Ω h ), and the test function of w h + τu n h · ∇w h is multiplied on the both sides of equation (9). erefore, we are able to get the following spatial discretization forms of the two component equations of equation (9): e other subequations are the elliptic equation. us, we employ the standard FEM to solve them. e spatial discretization form of equation (11) is shown in the following equation: e FEM discretization of the two component equations of equation (12) is that In the same way, the spatial discretizations of the two component equations of equation (13) are displayed as

Solver for Level Set and Its Reinitialization Equations.
Furthermore, we employ the implicit discontinuous Galerkin method for the solution of level set and its reinitialization equations. We take the solution of level set equation as an example to explain the full process. e implicit temporal discretization of equation (2) is written as follows: As for the discontinuous Galerkin spatial discretization, we multiply the test function l h on the both sides of equation (24) and also conduct the integration by parts twice for the convective term. e spatial discretization form of the convective term is that where ϕ h ′ is in the D k (Ω h ) space.
As for the numerical flux of (u h ′ ϕ h ′ ) * , we choose the Lax-Friedrichs flux in equation (17) to substitute it: us, the final spatial discretization of equation (24) is written as In the same way, the fully discretization scheme of equation (3) is written as e time step is selected according to the following equation [22,23]: where h denotes the minimum grid size and N represents the degree of polynomials. In order to keep the mass conservation property, we add a mass correction step in our computational framework. Let S represent the mass of fluid obtained via a numerical scheme. Se is the exact mass of fluid, and L denotes the length of free surface. en, we can calculate the value of c according to the equality of c � (Se − S)/L. After that, we are able to modify the value of level set function ϕ to ϕ − c. e details of mass correction technique are found in the literature [21,22].

Computational Process.
After the description of the temporal and spatial discretizations, we give the outline of computational strategy in this part. e details are shown as follows: (1) Initialize the level set function, velocity, and pressure (2) For n � 0, ... , N (t � 0 ⟶ T) (a) Solve equations (18) and (19) to receive the intermediate velocity u h and v h (b) Deal with equation (20) to achieve the pressure field (c) Handle equations (21) and (22) to get the intermediate velocity u h and v h (d) Solve the level set and its reinitialization equations to update the moving interface front (e) Take advantage of the technique to make mass correction (f ) If the fluid is non-Newtonian, it needs to update its viscosity according to the power-law model (g) Output the results at this step

Validity of Mathematical Model.
In this part, we consider the classical dam break-free surface flow to validate our mathematical model. A schematic diagram is given in Figure 1. In the beginning, the water is sustained in the region of (0,1) × (0,1) and a � 1 m. e nonslip boundary conditions are used for the velocity on the solid walls. We set pressure as zero on the upper boundary. e computational mesh is shown in Figure 2, and there are 32088 elements. All the simulations are executed on the Dell desktop computer with i7-9700 CPU @ 3.00 GHz, and the Mathematical Problems in Engineering 5 CPU time is 21860 s. We have compared our numerical results with the experimental data in Figures 3 and 4 [31,32]. e profile of free surface, the position of surge front, and the remaining water column height on the left wall are all identical with the results in the literature [31,32], which demonstrate the validity of our mathematical model.

e Sloshing Free Surface Flow.
In this section, we study the sloshing free surface fluid flows in different tanks. e initial geometries of the sloshing problem are depicted in       and the radius is 0.1 B. We will consider the Newtonian (water) and non-Newtonian fluids (pseudoplastic fluid and the dilatant fluid) in the following simulations, and the firstorder polynomials are employed.

e Sloshing Free Surface of Newtonian Fluid.
We first investigate the water sloshing in the tank with smooth bottom. On the solid walls, we take nonslip boundary conditions for the velocity. e pressure is set as zero on the upper boundary. In the computation, we have employed three successively refined Mesh1, Mesh2, and Mesh3. e minimum gird size is 1/120, 1/150, and 1/180, and the total triangular elements are 10002, 13590, and 22778. e unstructured Mesh2 is exhibited in Figure 6. e time step Δt is selected as 0.001.
According to the initial geometry of the problem, the fluid will go down on the left wall due to the gravity. In the meanwhile, the fluid will go up along the right wall because of the extrusion from the left side. When the fluid reaches the highest point on the right wall and the lowest point on the left wall, it will then go down on the right wall and run up on the left wall. When the fluid gets the highest point on the left wall and the lowest point on the right wall, this is one circle. After that, the fluid will continue to flow periodically.
In Figures 7(a) and 7(b), the profile of free surface and the value of pressure along y � 0.1 at t � 0.88 s is depicted. e red, blue, and black lines represent the results from Mesh1, Mesh2, and Mesh3, respectively. From these two figures, our numerical algorithm shows good mesh convergence. e maximum mass deviations on three meshes are also given in Table 1 Figure 9. It is observed that the values of pressure are sensible with the corresponding shape of interface front. We also provide the other author's results [13] about the distribution of pressure at different time instants in Figure 10. e profiles of free surface and distributions of pressure are all identical with the reports in the literature [12,13], which illustrate the validity of our computational strategy. e maximum value of pressure in our simulation is about 2850, and it is 2600 for the report in the literature. e relative error is about 8.9%. e case of initial nonlinear free surface, i.e., sine function, is also simulated. e profiles of free surface at several different time instants (a) t � 0, (b) t � 0.8 s, (c) t � 1.44 s, and (d) t � 1.84 s are displayed in Figure 11. e development of free surface is reasonable, and there is no numerical oscillation, which illustrates the good ability of our computational scheme.
In addition, we also consider the water sloshing in the second tank with rough bottom. We have employed three successively refined Mesh1, Mesh2, and Mesh3 with 9200, 12639, and 20778 elements, respectively. e minimum grid size is 1/120, 1/150, and 1/180. e computational time step Δt is 0.001. In Figure 12   respectively. e red, blue, and black lines represent the results from Mesh1, Mesh2, and Mesh3, which demonstrate the good mesh convergence of the computational method. In Table 2, the maximum mass deviations on three refined meshes are given. As for the simulations in the tank with rough bottom, the CPU time is 5400 s, 7560 s, and 12588 s on Mesh1, Mesh2, and Mesh3, respectively. e profiles of moving interface at different time instants are also shown in Figure 14. e fluid on the left wall at t � 0.88 s is much higher than the result in Figure 8(a) that is mainly because of the convex bottom of this tank. e distributions of pressure at different time instants are also depicted in Figure 15. Although the computational geometry is more complex with some singular points, there is no numerical oscillation for the interface and pressure, which demonstrates the robustness of our numerical scheme.

e Sloshing Free Surface of Pseudoplastic Fluid.
In this part, we consider the sloshing of pseudoplastic fluid in the two tanks. e boundary conditions and the time step are in accordance with those in Section 4.2.1. e computational Mesh2 of two tanks is utilized in the following simulations. As for the power-law index in equation (7), we take it as n � 0.75.
is means the viscosity of this non-Newtonian fluid will decrease with the increase in shear rate. e free surface in the first tank at t � 0.88 s, 2.06 s, 3.2 s, and 5.0 s is shown in Figure 16. We are able to find that the fluid on the left wall at t � 0.88 s is much higher than the same case in Section 4.2.1. e free surface at different times is coarser than those in Section 4.2.1. ese phenomena make sense because of the shear-thinning property of the fluid. e CPU time is 8580 s for this simulation on Mesh2.
We also show the development of free surface for the pseudoplastic fluid in the second tank. e fluid on the left wall at t � 0.88 s is also much higher than the result in Figure 16. ere is no instability for the evolution of pseudoplastic fluid in the second tank with complex geometry. e CPU time is 7990 s for this simulation on Mesh2.

e Sloshing Free Surface of Dilatant Fluid.
In this part, we continue to study the sloshing of another non-Newtonian fluid, i.e., the dilatant fluid, in the two tanks. e    boundary conditions and the time step are in accordance with those in Section 4.2.1. e computational Mesh2 in Section 4.2.1 for two tanks is employed in the following simulations. e power-law index of n in equation (7) is chosen as 1.5. is means the viscosity of this non-Newtonian fluid will increase with the enlargement of shear rate.
Comparing with the results in Figures 8 and 17, it is able to know that, with the increase in n, the sloshing amplitude of free surface becomes smaller and the profile of interface seems more compact and smoother. e vertical velocity along y � 0.1 at t � 0.33 s in the first and second tanks is depicted in Figures 18(a) and 18(b), respectively. e red, blue, and black lines represent the velocity of the pseudoplastic fluid, the Newtonian fluid, and the dilatant fluid, respectively. In the vicinity of the boundary, the velocity of the pseudoplastic fluid is sharpest due to its shear-thinning property. e velocity boundary layer of dilatant fluid is wider than the other two fluids, which is mainly because its viscosity tends to be bigger with the increase in shear rate. ese findings are identical with the reports in the literature [26]. e CPU time is 8580 s for this simulation on Mesh2. e sloshing free surface of dilatant fluid in the second tank at four different time instants is also given in Figure 19. e fluid on the left wall at t � 0.88 s in the second tank is higher than that in the first tank. e free surface looks smoother with the comparison of the results in Figures 14 and 20. e CPU time is 7990 s for this simulation on Mesh2.

Conclusions
In this paper, we have proposed a SUPG finite element coupled with implicit discontinous Galerkin method to simulate the sloshing free surface fluid flow in two rectangular tanks. We have simulated this problem via a twophase model to avoid the mesh movement or reconstruction. e Newtonian and non-Newtonian fluids sloshing in the tank with smooth and convex rough bottom have all been studied. e profiles of free surface and the distributions of pressure for the Newtonian sloshing problem are all identical with the results in the literature, which shows the validity of our coupled method. ere exists no numerical oscillation for all the cases, which illustrate the robustness of the numerical algorithm. In addition, the simulation results show good mesh convergence of our combined approach for all the cases. e good mass conservation property is also demonstrated in the whole flow process.
For the same tank, the shake amplitude of the non-Newtonian fluid with n � 0.75 is biggest and it is smallest for the non-Newtonian fluid with n � 1.5. e interface front is smooth and compact for the non-Newtonian fluid of n � 1.5, and it becomes coarser for the non-Newtonian fluid with n � 0.75. It is also able to find that, in the vicinity of boundary, the velocity of the pseudoplastic fluid (n � 0.75) is sharpest and it is smoothest for the dilatant fluid (n � 1.5).

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e author declares that there are no conflicts of interest.