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

Authors: Puyang Gao
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

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.


---

## Body

## 1. Introduction

The 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. The 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–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. They 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. They 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. The 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. The 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–11] is proposed and has been used for the simulation of free surface fluid flow. The 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. The 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. The 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. The streamline upwind/Petrov–Galerkin (SUPG) method [14, 15] is utilized to solve the hyperbolic subequation for the stability. The 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.The 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 Section3, 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.


## 2. 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. Therefore, 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. The details are shown as follows.


### 2.1. The Level Set Equation

The 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=minxΓ∈Γx−xΓ. In the liquid area, the value of level set function ϕ is less than zero and it is bigger than zero in the air region.The level set equation is written as follows [23]:(1)∂ϕ∂t+u⋅∇ϕ=0.As for the incompressible flow, the above equation is rewritten in the following conservation form:(2)∂ϕ∂t+∇⋅uϕ=0.In order to keep the signed distance function property of the level set function, the reinitialization procedure is essential. The conservative form of level set reinitialization equation with initial condition is given in the following equation [23]:(3)∂φ∂τ+∇⋅wφ=signφ0+∇⋅wφ,φx,0=φ0x,where x  =  (x, y), w=signφ0∇φ/∇φ, and signφ0=φ0/φ02+ε2. φ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.


### 2.2. A Unified Form of Governing Equation

As for the incompressible liquid and gas, the governing equations for each fluid are both Navier–Stokes equation. The 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. The unified momentum equation and the continuity equation are written as(4)ρϕ∂u∂t+ρϕu⋅∇u=−∇p+μϕ∇2un+1+f,∇⋅u=0,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]:(5)ρϕ=ρl+ρg−ρlHϕ,μϕ=μl+μg−μlHϕ,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:(6)Hϕ=0,ϕ<−ε,121+ϕε+1πsinπϕε,−ε≤ϕ≤ε,1,ϕ>ε,where ε is selected as 1.2 h in the calculation.


### 2.3. 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. The strain rate tensor is denoted as d=∇u+∇uT/2, and its magnitude is γ=2d:d. Thus, in terms of power-law model, the viscosity is(7)μl=μlγn−1,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.


## 2.1. The Level Set Equation

The 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=minxΓ∈Γx−xΓ. In the liquid area, the value of level set function ϕ is less than zero and it is bigger than zero in the air region.The level set equation is written as follows [23]:(1)∂ϕ∂t+u⋅∇ϕ=0.As for the incompressible flow, the above equation is rewritten in the following conservation form:(2)∂ϕ∂t+∇⋅uϕ=0.In order to keep the signed distance function property of the level set function, the reinitialization procedure is essential. The conservative form of level set reinitialization equation with initial condition is given in the following equation [23]:(3)∂φ∂τ+∇⋅wφ=signφ0+∇⋅wφ,φx,0=φ0x,where x  =  (x, y), w=signφ0∇φ/∇φ, and signφ0=φ0/φ02+ε2. φ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.


## 2.2. A Unified Form of Governing Equation

As for the incompressible liquid and gas, the governing equations for each fluid are both Navier–Stokes equation. The 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. The unified momentum equation and the continuity equation are written as(4)ρϕ∂u∂t+ρϕu⋅∇u=−∇p+μϕ∇2un+1+f,∇⋅u=0,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]:(5)ρϕ=ρl+ρg−ρlHϕ,μϕ=μl+μg−μlHϕ,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:(6)Hϕ=0,ϕ<−ε,121+ϕε+1πsinπϕε,−ε≤ϕ≤ε,1,ϕ>ε,where ε is selected as 1.2 h in the calculation.


## 2.3. 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. The strain rate tensor is denoted as d=∇u+∇uT/2, and its magnitude is γ=2d:d. Thus, in terms of power-law model, the viscosity is(7)μl=μlγn−1,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.


## 3. 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.


### 3.1. 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. The 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:(8)ρϕun+1−unΔt+ρϕun⋅∇un=−∇pn+1+μϕ∇2un+1+f.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 un+1. It is then able to get the hyperbolic subequation:(9)ρϕu˜−unΔt+ρϕun⋅∇un=0.After that, we bring in another intermediate velocity and achieve an equation about pressure:(10)ρϕu˜˜−u˜Δt=−∇pn+1.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:(11)∇2pn+1=ρϕΔt∇⋅u˜.And then, the intermediate velocityu˜˜ is able to calculate based on the following equation:(12)u˜˜=u˜−Δt∇pn+1ρϕ.In the last stage, it is able to acquire the Helmholtz equation about velocity:(13)−μϕ∇2un+1+ρϕΔtun+1=ρϕΔtu˜˜+f.When we plus the subequations of (9), (10), and (13) together, we are able to recover equation (8).


### 3.2. 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 intoN nonoverlapping small conformal triangles and h represents the maximum element diameter. Thus, the area Ω is able to be approximated by Ωh=∪i=1NΩk. The continuous and discontinuous spaces are defined as(14)CkΩh=v∈C0Ω:v|Ωi∈PkΩi,∀Ωi∈Ωh,DkΩh=l∈L2Ω:l|Ωi∈PkΩi,∀Ωi∈Ωh,where Pk (Ωi) represents the polynomials of degree less than or equal to k on Ωi. Ck2Ωh and Dk2Ω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. The average operator [30] is defined as(15)u=12u−+u+.It is effective for both scalar and vector. Furthermore, the jump operators [30] for the scalar and vector are also given as(16)u=n^−u−+n^+u+,u=n^−⋅u−+n^+⋅u+.According to the average and jump operators, the popular Lax–Friedrichs numerical flux [30] is able to be written as follows:(17)f∗u−,u+=fu+λ2u,λ≥maxs∈u−,u+∂fs∂u.


### 3.3. Spatial Discretizations of the Subequations

After the temporal discretization and splitting scheme, we receive several semidiscrete subequations. The 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 Ck2Ωh, and the test function of wh+τuhn⋅∇wh is multiplied on the both sides of equation (9). Therefore, we are able to get the following spatial discretization forms of the two component equations of equation (9):(18)∫Ωhρϕu˜h−uhnΔtw1h+τuhn⋅∇w1hdΩh−∫Ωhρϕuhn⋅∇uhn:w1h+τuhn⋅∇w1hdΩh=0,∀w1h∈CkΩh,(19)∫Ωhρϕv˜h−vhnΔtw2h+τuhn⋅∇w2hdΩh−∫Ωhρϕuhn⋅∇vhn:w2h+τuhn⋅∇w2hdΩh=0,∀w2h∈CkΩh,where uh=uh,vh, wh=w1h,w2h, and τ=h/2uh.The other subequations are the elliptic equation. Thus, we employ the standard FEM to solve them. The spatial discretization form of equation (11) is shown in the following equation:(20)−∫Ωh∇phn+1⋅∇w1hdΩh+∫∂Ωh∂phn+1∂nw1hds=ρϕΔt∫Ωh∇⋅u˜w1hdΩh,∀w1h∈CkΩh.The FEM discretization of the two component equations of equation (12) is that(21)∫Ωhu˜˜hw1hdΩh=∫Ωhu˜h−pxn+1Δtρϕw1hdΩh,(22)∫Ωhv˜˜hw1hdΩh=∫Ωhv˜h−pyn+1Δtρϕw1hdΩh.In the same way, the spatial discretizations of the two component equations of equation (13) are displayed as(23)∫Ωhμϕ∇uhn+1⋅∇w1hdΩh+∫ΩhρϕΔtuhn+1w1hdΩh=∫ΩhρϕΔtu˜˜hw1hdΩh,∀w1h∈CkΩh,.∫Ωhμϕ∇vhn+1⋅∇w2hdΩh+∫ΩhρϕΔtvhn+1w2hdΩh=∫ΩhρϕΔtv˜˜h−gw2hdΩh,∀w2h∈CkΩh.


### 3.4. 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. The implicit temporal discretization of equation (2) is written as follows:(24)ϕn+1−ϕndt+∇⋅uϕn+1=0.As for the discontinuous Galerkin spatial discretization, we multiply the test functionlh on the both sides of equation (24) and also conduct the integration by parts twice for the convective term. The spatial discretization form of the convective term is that(25)∫Ωi∇⋅uh′ϕhlhdΩi=∫Ωi∇⋅uh′ϕh′lhdΩi−∫∂Ωin⋅uh′ϕh′lhds+∫∂Ωin⋅uh′ϕh′∗lhds,lh∈DkΩi,where ϕh′ is in the DkΩh space.As for the numerical flux ofuh′ϕh′∗, we choose the Lax–Friedrichs flux in equation (17) to substitute it:(26)∇⋅uh′ϕh′,lhΩi=∇⋅uh′ϕh′,lhΩi−12uh′ϕh′,lh∂Ωi+λ2n⋅ϕh′,lh∂Ωi,lh∈DkΩi.Thus, the final spatial discretization of equation (24) is written as(27)ϕh′n+1−ϕh′ndt,lhΩi+∇⋅uh′ϕh′n+1,lhΩi=0,lh∈DkΩi.In the same way, the fully discretization scheme of equation (3) is written as(28)ϕh′n+1−ϕh′ndt,lhΩi+∇⋅wh′ϕh′n+1,lhΩi−12wh′ϕh′n+1,lh∂Ωi+λ2n⋅ϕh′n+1,lh∂Ωi=signφ0+∇⋅wh′n+1ϕh′n+1,lhΩi,lh∈DkΩi.The time step is selected according to the following equation [22, 23]:(29)Δt≤h2N+11u,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. LetS 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. Then, 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. The details of mass correction technique are found in the literature [21, 22].


### 3.5. Computational Process

After the description of the temporal and spatial discretizations, we give the outline of computational strategy in this part. The details are shown as follows:(1)
Initialize the level set function, velocity, and pressure(2)
Forn=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(3)
Output the final results and terminate the calculation


## 3.1. 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. The 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:(8)ρϕun+1−unΔt+ρϕun⋅∇un=−∇pn+1+μϕ∇2un+1+f.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 un+1. It is then able to get the hyperbolic subequation:(9)ρϕu˜−unΔt+ρϕun⋅∇un=0.After that, we bring in another intermediate velocity and achieve an equation about pressure:(10)ρϕu˜˜−u˜Δt=−∇pn+1.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:(11)∇2pn+1=ρϕΔt∇⋅u˜.And then, the intermediate velocityu˜˜ is able to calculate based on the following equation:(12)u˜˜=u˜−Δt∇pn+1ρϕ.In the last stage, it is able to acquire the Helmholtz equation about velocity:(13)−μϕ∇2un+1+ρϕΔtun+1=ρϕΔtu˜˜+f.When we plus the subequations of (9), (10), and (13) together, we are able to recover equation (8).


## 3.2. 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 intoN nonoverlapping small conformal triangles and h represents the maximum element diameter. Thus, the area Ω is able to be approximated by Ωh=∪i=1NΩk. The continuous and discontinuous spaces are defined as(14)CkΩh=v∈C0Ω:v|Ωi∈PkΩi,∀Ωi∈Ωh,DkΩh=l∈L2Ω:l|Ωi∈PkΩi,∀Ωi∈Ωh,where Pk (Ωi) represents the polynomials of degree less than or equal to k on Ωi. Ck2Ωh and Dk2Ω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. The average operator [30] is defined as(15)u=12u−+u+.It is effective for both scalar and vector. Furthermore, the jump operators [30] for the scalar and vector are also given as(16)u=n^−u−+n^+u+,u=n^−⋅u−+n^+⋅u+.According to the average and jump operators, the popular Lax–Friedrichs numerical flux [30] is able to be written as follows:(17)f∗u−,u+=fu+λ2u,λ≥maxs∈u−,u+∂fs∂u.


## 3.3. Spatial Discretizations of the Subequations

After the temporal discretization and splitting scheme, we receive several semidiscrete subequations. The 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 Ck2Ωh, and the test function of wh+τuhn⋅∇wh is multiplied on the both sides of equation (9). Therefore, we are able to get the following spatial discretization forms of the two component equations of equation (9):(18)∫Ωhρϕu˜h−uhnΔtw1h+τuhn⋅∇w1hdΩh−∫Ωhρϕuhn⋅∇uhn:w1h+τuhn⋅∇w1hdΩh=0,∀w1h∈CkΩh,(19)∫Ωhρϕv˜h−vhnΔtw2h+τuhn⋅∇w2hdΩh−∫Ωhρϕuhn⋅∇vhn:w2h+τuhn⋅∇w2hdΩh=0,∀w2h∈CkΩh,where uh=uh,vh, wh=w1h,w2h, and τ=h/2uh.The other subequations are the elliptic equation. Thus, we employ the standard FEM to solve them. The spatial discretization form of equation (11) is shown in the following equation:(20)−∫Ωh∇phn+1⋅∇w1hdΩh+∫∂Ωh∂phn+1∂nw1hds=ρϕΔt∫Ωh∇⋅u˜w1hdΩh,∀w1h∈CkΩh.The FEM discretization of the two component equations of equation (12) is that(21)∫Ωhu˜˜hw1hdΩh=∫Ωhu˜h−pxn+1Δtρϕw1hdΩh,(22)∫Ωhv˜˜hw1hdΩh=∫Ωhv˜h−pyn+1Δtρϕw1hdΩh.In the same way, the spatial discretizations of the two component equations of equation (13) are displayed as(23)∫Ωhμϕ∇uhn+1⋅∇w1hdΩh+∫ΩhρϕΔtuhn+1w1hdΩh=∫ΩhρϕΔtu˜˜hw1hdΩh,∀w1h∈CkΩh,.∫Ωhμϕ∇vhn+1⋅∇w2hdΩh+∫ΩhρϕΔtvhn+1w2hdΩh=∫ΩhρϕΔtv˜˜h−gw2hdΩh,∀w2h∈CkΩh.


## 3.4. 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. The implicit temporal discretization of equation (2) is written as follows:(24)ϕn+1−ϕndt+∇⋅uϕn+1=0.As for the discontinuous Galerkin spatial discretization, we multiply the test functionlh on the both sides of equation (24) and also conduct the integration by parts twice for the convective term. The spatial discretization form of the convective term is that(25)∫Ωi∇⋅uh′ϕhlhdΩi=∫Ωi∇⋅uh′ϕh′lhdΩi−∫∂Ωin⋅uh′ϕh′lhds+∫∂Ωin⋅uh′ϕh′∗lhds,lh∈DkΩi,where ϕh′ is in the DkΩh space.As for the numerical flux ofuh′ϕh′∗, we choose the Lax–Friedrichs flux in equation (17) to substitute it:(26)∇⋅uh′ϕh′,lhΩi=∇⋅uh′ϕh′,lhΩi−12uh′ϕh′,lh∂Ωi+λ2n⋅ϕh′,lh∂Ωi,lh∈DkΩi.Thus, the final spatial discretization of equation (24) is written as(27)ϕh′n+1−ϕh′ndt,lhΩi+∇⋅uh′ϕh′n+1,lhΩi=0,lh∈DkΩi.In the same way, the fully discretization scheme of equation (3) is written as(28)ϕh′n+1−ϕh′ndt,lhΩi+∇⋅wh′ϕh′n+1,lhΩi−12wh′ϕh′n+1,lh∂Ωi+λ2n⋅ϕh′n+1,lh∂Ωi=signφ0+∇⋅wh′n+1ϕh′n+1,lhΩi,lh∈DkΩi.The time step is selected according to the following equation [22, 23]:(29)Δt≤h2N+11u,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. LetS 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. Then, 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. The details of mass correction technique are found in the literature [21, 22].


## 3.5. Computational Process

After the description of the temporal and spatial discretizations, we give the outline of computational strategy in this part. The details are shown as follows:(1)
Initialize the level set function, velocity, and pressure(2)
Forn=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(3)
Output the final results and terminate the calculation


## 4. Numerical Results

### 4.1. 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 Figure1. In the beginning, the water is sustained in the region of (0,1) × (0,1) and a = 1 m. The nonslip boundary conditions are used for the velocity on the solid walls. We set pressure as zero on the upper boundary. The 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 CPU time is 21860 s. We have compared our numerical results with the experimental data in Figures 3 and 4 [31, 32]. The 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.Figure 1
The schematic diagram of dam break-free surface flow.Figure 2
The computational mesh.Figure 3
The free surface at (a)t = 0.1 s, (b) t = 0.2 s, (c) t = 0.4 s, and (d) t = 0.6 s (experimental data in the left column and numerical simulations in the right column).
(a)(b)(c)(d)Figure 4
Comparison between the numerical results on Mesh2 and the experimental data: (a) surge front position; (b) remaining water column height.
(a)(b)


### 4.2. The Sloshing Free Surface Flow

In this section, we study the sloshing free surface fluid flows in different tanks. The initial geometries of the sloshing problem are depicted in Figure5. The height of fluid on the left and right walls is H1 = 0.35 m and H2 = 0.15 m, respectively. The width of the tank is B = 0.5 m, and the height is H = 0.6 m. In the second tank, the rough bottom consists of five convex semicircles 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 first-order polynomials are employed.Figure 5
The initial geometry of the sloshing problem in the (a) first tank and (b) second tank.
(a)(b)


#### 4.2.1. The 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. The pressure is set as zero on the upper boundary. In the computation, we have employed three successively refined Mesh1, Mesh2, and Mesh3. The minimum gird size is 1/120, 1/150, and 1/180, and the total triangular elements are 10002, 13590, and 22778. The unstructured Mesh2 is exhibited in Figure6. The time step Δt is selected as 0.001.Figure 6
The computational Mesh2 for the first tank.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 Figures7(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. The 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. The maximum mass deviations on three meshes are also given in Table 1. The maximum mass error decreases with the increase in mesh elements. The results illustrate the good mass conservation property. As for the simulations in the tank with flat bottom, the CPU time is 6180 s, 8040 s, and 13507 s on Mesh1, Mesh2, and Mesh3, respectively.Figure 7
(a) The profile of free surface; (b) the pressure alongy = 0.1 at t = 0.88 s on three meshes (red, blue, and black lines denote results from Mesh1, Mesh2, and Mesh3).
(a)(b)Table 1
The maximum mass deviation on different meshes.
Mesh typeMesh1Mesh2Mesh3Maximum mass deviation2.2%1.3%0.8%The development of free surface is displayed in Figure8 at four different times: 0.88 s, 2.06 s, 3.2 s, and 5.0 s. In the entire tank, the contours of pressure at 0.88 s, 2.06 s, 3.2 s, and 5.0 s are exhibited in 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. The 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. The maximum value of pressure in our simulation is about 2850, and it is 2600 for the report in the literature. The relative error is about 8.9%.Figure 8
Development of free surface at several different time instants: (a) 0.88 s; (b) 2.06 s; (c) 3.2 s; (d) 5.0 s.
(a)(b)(c)(d)Figure 9
The pressure contours at four-time instants: (a) 0.88 s; (b) 2.06 s; (c) 3.2 s; (d) 5.0 s.
(a)(b)(c)(d)Figure 10
The distribution of pressure at different times [13]: (a) 0.88 s; (b) 2.06 s; (c) 5.0 s.
(a)(b)(c)The case of initial nonlinear free surface, i.e., sine function, is also simulated. The 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. The development of free surface is reasonable, and there is no numerical oscillation, which illustrates the good ability of our computational scheme.Figure 11
The profile of interface front at (a)t = 0, (b) t = 0.8 s, (c) t = 1.44 s, and (d) t = 1.84 s.
(a)(b)(c)(d)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. The minimum grid size is 1/120, 1/150, and 1/180. The computational time stepΔt is 0.001. In Figure 12, the unstructured Mesh2 is shown.Figure 12
The computational Mesh2 for the second tank.The interface front and the value of pressure alongy = 0.1 at t = 0.88 s are depicted in Figures 13(a) and 13(b), respectively. The 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.Figure 13
(a) The profile of free surface; (b) the pressure alongy = 0.1 at t = 0.88 s on three meshes (red, blue, and black lines denote results from Mesh1, Mesh2, and Mesh3).
(a)(b)Table 2
The maximum mass deviation on different meshes.
Mesh typeMesh1Mesh2Mesh3Maximum mass deviation2.4%1.5%0.9%The profiles of moving interface at different time instants are also shown in Figure14. The 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. The 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.Figure 14
The development of free surface at several different times: (a) 0.88 s; (b) 2.06 s; (c) 3.2 s; (d) 5.0 s.
(a)(b)(c)(d)Figure 15
The pressure contours at four time instants: (a) 0.88 s; (b) 2.06 s; (c) 3.2 s; (d) 5.0 s.
(a)(b)(c)(d)


#### 4.2.2. The Sloshing Free Surface of Pseudoplastic Fluid

In this part, we consider the sloshing of pseudoplastic fluid in the two tanks. The boundary conditions and the time step are in accordance with those inSection 4.2.1. The 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. This means the viscosity of this non-Newtonian fluid will decrease with the increase in shear rate.The free surface in the first tank att = 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. The free surface at different times is coarser than those in Section 4.2.1. These phenomena make sense because of the shear-thinning property of the fluid. The CPU time is 8580 s for this simulation on Mesh2.Figure 16
Free surface of pseudoplastic fluid at different times: (a)t = 0.88 s, (b) t = 2.06 s, (c) t = 3.2 s, and (d) t = 5.0 s in the first tank.
(a)(b)(c)(d)We also show the development of free surface for the pseudoplastic fluid in the second tank. The fluid on the left wall att = 0.88 s is also much higher than the result in Figure 16. There is no instability for the evolution of pseudoplastic fluid in the second tank with complex geometry. The CPU time is 7990 s for this simulation on Mesh2.


#### 4.2.3. The 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. The boundary conditions and the time step are in accordance with those inSection 4.2.1. The computational Mesh2 in Section 4.2.1 for two tanks is employed in the following simulations. The power-law index of n in equation (7) is chosen as 1.5. This 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.Figure 17
Free surface of dilatant fluid at different times: (a)t = 0.88 s, (b) t = 2.06 s, (c) t = 3.2 s, and (d) t = 5.0 s in the first tank.
(a)(b)(c)(d)The vertical velocity alongy = 0.1 at t = 0.33 s in the first and second tanks is depicted in Figures 18(a) and 18(b), respectively. The 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. The 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. These findings are identical with the reports in the literature [26]. The CPU time is 8580 s for this simulation on Mesh2.Figure 18
The vertical velocity of the pseudoplastic fluid (red line), Newtonian fluid (blue line), and the dilatant fluid (black line) in the (a) first tank and (b) second tank.
(a)(b)The sloshing free surface of dilatant fluid in the second tank at four different time instants is also given in Figure19. The fluid on the left wall at t = 0.88 s in the second tank is higher than that in the first tank. The free surface looks smoother with the comparison of the results in Figures 14 and 20. The CPU time is 7990 s for this simulation on Mesh2.Figure 19
Free surface of dilatant fluid at different times: (a)t = 0.88 s, (b) t = 2.06 s, (c) t = 3.2 s, and (d) t = 5.0 s in the second tank.
(a)(b)(c)(d)Figure 20
Free surface of pseudoplastic fluid at different times: (a)t = 0.88 s, (b) t = 2.06 s, (c) t = 3.2 s, and (d) t = 5.0 s in the second tank.
(a)(b)(c)(d)


## 4.1. 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 Figure1. In the beginning, the water is sustained in the region of (0,1) × (0,1) and a = 1 m. The nonslip boundary conditions are used for the velocity on the solid walls. We set pressure as zero on the upper boundary. The 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 CPU time is 21860 s. We have compared our numerical results with the experimental data in Figures 3 and 4 [31, 32]. The 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.Figure 1
The schematic diagram of dam break-free surface flow.Figure 2
The computational mesh.Figure 3
The free surface at (a)t = 0.1 s, (b) t = 0.2 s, (c) t = 0.4 s, and (d) t = 0.6 s (experimental data in the left column and numerical simulations in the right column).
(a)(b)(c)(d)Figure 4
Comparison between the numerical results on Mesh2 and the experimental data: (a) surge front position; (b) remaining water column height.
(a)(b)


## 4.2. The Sloshing Free Surface Flow

In this section, we study the sloshing free surface fluid flows in different tanks. The initial geometries of the sloshing problem are depicted in Figure5. The height of fluid on the left and right walls is H1 = 0.35 m and H2 = 0.15 m, respectively. The width of the tank is B = 0.5 m, and the height is H = 0.6 m. In the second tank, the rough bottom consists of five convex semicircles 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 first-order polynomials are employed.Figure 5
The initial geometry of the sloshing problem in the (a) first tank and (b) second tank.
(a)(b)


### 4.2.1. The 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. The pressure is set as zero on the upper boundary. In the computation, we have employed three successively refined Mesh1, Mesh2, and Mesh3. The minimum gird size is 1/120, 1/150, and 1/180, and the total triangular elements are 10002, 13590, and 22778. The unstructured Mesh2 is exhibited in Figure6. The time step Δt is selected as 0.001.Figure 6
The computational Mesh2 for the first tank.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 Figures7(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. The 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. The maximum mass deviations on three meshes are also given in Table 1. The maximum mass error decreases with the increase in mesh elements. The results illustrate the good mass conservation property. As for the simulations in the tank with flat bottom, the CPU time is 6180 s, 8040 s, and 13507 s on Mesh1, Mesh2, and Mesh3, respectively.Figure 7
(a) The profile of free surface; (b) the pressure alongy = 0.1 at t = 0.88 s on three meshes (red, blue, and black lines denote results from Mesh1, Mesh2, and Mesh3).
(a)(b)Table 1
The maximum mass deviation on different meshes.
Mesh typeMesh1Mesh2Mesh3Maximum mass deviation2.2%1.3%0.8%The development of free surface is displayed in Figure8 at four different times: 0.88 s, 2.06 s, 3.2 s, and 5.0 s. In the entire tank, the contours of pressure at 0.88 s, 2.06 s, 3.2 s, and 5.0 s are exhibited in 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. The 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. The maximum value of pressure in our simulation is about 2850, and it is 2600 for the report in the literature. The relative error is about 8.9%.Figure 8
Development of free surface at several different time instants: (a) 0.88 s; (b) 2.06 s; (c) 3.2 s; (d) 5.0 s.
(a)(b)(c)(d)Figure 9
The pressure contours at four-time instants: (a) 0.88 s; (b) 2.06 s; (c) 3.2 s; (d) 5.0 s.
(a)(b)(c)(d)Figure 10
The distribution of pressure at different times [13]: (a) 0.88 s; (b) 2.06 s; (c) 5.0 s.
(a)(b)(c)The case of initial nonlinear free surface, i.e., sine function, is also simulated. The 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. The development of free surface is reasonable, and there is no numerical oscillation, which illustrates the good ability of our computational scheme.Figure 11
The profile of interface front at (a)t = 0, (b) t = 0.8 s, (c) t = 1.44 s, and (d) t = 1.84 s.
(a)(b)(c)(d)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. The minimum grid size is 1/120, 1/150, and 1/180. The computational time stepΔt is 0.001. In Figure 12, the unstructured Mesh2 is shown.Figure 12
The computational Mesh2 for the second tank.The interface front and the value of pressure alongy = 0.1 at t = 0.88 s are depicted in Figures 13(a) and 13(b), respectively. The 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.Figure 13
(a) The profile of free surface; (b) the pressure alongy = 0.1 at t = 0.88 s on three meshes (red, blue, and black lines denote results from Mesh1, Mesh2, and Mesh3).
(a)(b)Table 2
The maximum mass deviation on different meshes.
Mesh typeMesh1Mesh2Mesh3Maximum mass deviation2.4%1.5%0.9%The profiles of moving interface at different time instants are also shown in Figure14. The 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. The 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.Figure 14
The development of free surface at several different times: (a) 0.88 s; (b) 2.06 s; (c) 3.2 s; (d) 5.0 s.
(a)(b)(c)(d)Figure 15
The pressure contours at four time instants: (a) 0.88 s; (b) 2.06 s; (c) 3.2 s; (d) 5.0 s.
(a)(b)(c)(d)


### 4.2.2. The Sloshing Free Surface of Pseudoplastic Fluid

In this part, we consider the sloshing of pseudoplastic fluid in the two tanks. The boundary conditions and the time step are in accordance with those inSection 4.2.1. The 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. This means the viscosity of this non-Newtonian fluid will decrease with the increase in shear rate.The free surface in the first tank att = 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. The free surface at different times is coarser than those in Section 4.2.1. These phenomena make sense because of the shear-thinning property of the fluid. The CPU time is 8580 s for this simulation on Mesh2.Figure 16
Free surface of pseudoplastic fluid at different times: (a)t = 0.88 s, (b) t = 2.06 s, (c) t = 3.2 s, and (d) t = 5.0 s in the first tank.
(a)(b)(c)(d)We also show the development of free surface for the pseudoplastic fluid in the second tank. The fluid on the left wall att = 0.88 s is also much higher than the result in Figure 16. There is no instability for the evolution of pseudoplastic fluid in the second tank with complex geometry. The CPU time is 7990 s for this simulation on Mesh2.


### 4.2.3. The 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. The boundary conditions and the time step are in accordance with those inSection 4.2.1. The computational Mesh2 in Section 4.2.1 for two tanks is employed in the following simulations. The power-law index of n in equation (7) is chosen as 1.5. This 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.Figure 17
Free surface of dilatant fluid at different times: (a)t = 0.88 s, (b) t = 2.06 s, (c) t = 3.2 s, and (d) t = 5.0 s in the first tank.
(a)(b)(c)(d)The vertical velocity alongy = 0.1 at t = 0.33 s in the first and second tanks is depicted in Figures 18(a) and 18(b), respectively. The 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. The 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. These findings are identical with the reports in the literature [26]. The CPU time is 8580 s for this simulation on Mesh2.Figure 18
The vertical velocity of the pseudoplastic fluid (red line), Newtonian fluid (blue line), and the dilatant fluid (black line) in the (a) first tank and (b) second tank.
(a)(b)The sloshing free surface of dilatant fluid in the second tank at four different time instants is also given in Figure19. The fluid on the left wall at t = 0.88 s in the second tank is higher than that in the first tank. The free surface looks smoother with the comparison of the results in Figures 14 and 20. The CPU time is 7990 s for this simulation on Mesh2.Figure 19
Free surface of dilatant fluid at different times: (a)t = 0.88 s, (b) t = 2.06 s, (c) t = 3.2 s, and (d) t = 5.0 s in the second tank.
(a)(b)(c)(d)Figure 20
Free surface of pseudoplastic fluid at different times: (a)t = 0.88 s, (b) t = 2.06 s, (c) t = 3.2 s, and (d) t = 5.0 s in the second tank.
(a)(b)(c)(d)


## 4.2.1. The 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. The pressure is set as zero on the upper boundary. In the computation, we have employed three successively refined Mesh1, Mesh2, and Mesh3. The minimum gird size is 1/120, 1/150, and 1/180, and the total triangular elements are 10002, 13590, and 22778. The unstructured Mesh2 is exhibited in Figure6. The time step Δt is selected as 0.001.Figure 6
The computational Mesh2 for the first tank.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 Figures7(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. The 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. The maximum mass deviations on three meshes are also given in Table 1. The maximum mass error decreases with the increase in mesh elements. The results illustrate the good mass conservation property. As for the simulations in the tank with flat bottom, the CPU time is 6180 s, 8040 s, and 13507 s on Mesh1, Mesh2, and Mesh3, respectively.Figure 7
(a) The profile of free surface; (b) the pressure alongy = 0.1 at t = 0.88 s on three meshes (red, blue, and black lines denote results from Mesh1, Mesh2, and Mesh3).
(a)(b)Table 1
The maximum mass deviation on different meshes.
Mesh typeMesh1Mesh2Mesh3Maximum mass deviation2.2%1.3%0.8%The development of free surface is displayed in Figure8 at four different times: 0.88 s, 2.06 s, 3.2 s, and 5.0 s. In the entire tank, the contours of pressure at 0.88 s, 2.06 s, 3.2 s, and 5.0 s are exhibited in 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. The 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. The maximum value of pressure in our simulation is about 2850, and it is 2600 for the report in the literature. The relative error is about 8.9%.Figure 8
Development of free surface at several different time instants: (a) 0.88 s; (b) 2.06 s; (c) 3.2 s; (d) 5.0 s.
(a)(b)(c)(d)Figure 9
The pressure contours at four-time instants: (a) 0.88 s; (b) 2.06 s; (c) 3.2 s; (d) 5.0 s.
(a)(b)(c)(d)Figure 10
The distribution of pressure at different times [13]: (a) 0.88 s; (b) 2.06 s; (c) 5.0 s.
(a)(b)(c)The case of initial nonlinear free surface, i.e., sine function, is also simulated. The 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. The development of free surface is reasonable, and there is no numerical oscillation, which illustrates the good ability of our computational scheme.Figure 11
The profile of interface front at (a)t = 0, (b) t = 0.8 s, (c) t = 1.44 s, and (d) t = 1.84 s.
(a)(b)(c)(d)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. The minimum grid size is 1/120, 1/150, and 1/180. The computational time stepΔt is 0.001. In Figure 12, the unstructured Mesh2 is shown.Figure 12
The computational Mesh2 for the second tank.The interface front and the value of pressure alongy = 0.1 at t = 0.88 s are depicted in Figures 13(a) and 13(b), respectively. The 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.Figure 13
(a) The profile of free surface; (b) the pressure alongy = 0.1 at t = 0.88 s on three meshes (red, blue, and black lines denote results from Mesh1, Mesh2, and Mesh3).
(a)(b)Table 2
The maximum mass deviation on different meshes.
Mesh typeMesh1Mesh2Mesh3Maximum mass deviation2.4%1.5%0.9%The profiles of moving interface at different time instants are also shown in Figure14. The 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. The 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.Figure 14
The development of free surface at several different times: (a) 0.88 s; (b) 2.06 s; (c) 3.2 s; (d) 5.0 s.
(a)(b)(c)(d)Figure 15
The pressure contours at four time instants: (a) 0.88 s; (b) 2.06 s; (c) 3.2 s; (d) 5.0 s.
(a)(b)(c)(d)


## 4.2.2. The Sloshing Free Surface of Pseudoplastic Fluid

In this part, we consider the sloshing of pseudoplastic fluid in the two tanks. The boundary conditions and the time step are in accordance with those inSection 4.2.1. The 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. This means the viscosity of this non-Newtonian fluid will decrease with the increase in shear rate.The free surface in the first tank att = 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. The free surface at different times is coarser than those in Section 4.2.1. These phenomena make sense because of the shear-thinning property of the fluid. The CPU time is 8580 s for this simulation on Mesh2.Figure 16
Free surface of pseudoplastic fluid at different times: (a)t = 0.88 s, (b) t = 2.06 s, (c) t = 3.2 s, and (d) t = 5.0 s in the first tank.
(a)(b)(c)(d)We also show the development of free surface for the pseudoplastic fluid in the second tank. The fluid on the left wall att = 0.88 s is also much higher than the result in Figure 16. There is no instability for the evolution of pseudoplastic fluid in the second tank with complex geometry. The CPU time is 7990 s for this simulation on Mesh2.


## 4.2.3. The 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. The boundary conditions and the time step are in accordance with those inSection 4.2.1. The computational Mesh2 in Section 4.2.1 for two tanks is employed in the following simulations. The power-law index of n in equation (7) is chosen as 1.5. This 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.Figure 17
Free surface of dilatant fluid at different times: (a)t = 0.88 s, (b) t = 2.06 s, (c) t = 3.2 s, and (d) t = 5.0 s in the first tank.
(a)(b)(c)(d)The vertical velocity alongy = 0.1 at t = 0.33 s in the first and second tanks is depicted in Figures 18(a) and 18(b), respectively. The 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. The 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. These findings are identical with the reports in the literature [26]. The CPU time is 8580 s for this simulation on Mesh2.Figure 18
The vertical velocity of the pseudoplastic fluid (red line), Newtonian fluid (blue line), and the dilatant fluid (black line) in the (a) first tank and (b) second tank.
(a)(b)The sloshing free surface of dilatant fluid in the second tank at four different time instants is also given in Figure19. The fluid on the left wall at t = 0.88 s in the second tank is higher than that in the first tank. The free surface looks smoother with the comparison of the results in Figures 14 and 20. The CPU time is 7990 s for this simulation on Mesh2.Figure 19
Free surface of dilatant fluid at different times: (a)t = 0.88 s, (b) t = 2.06 s, (c) t = 3.2 s, and (d) t = 5.0 s in the second tank.
(a)(b)(c)(d)Figure 20
Free surface of pseudoplastic fluid at different times: (a)t = 0.88 s, (b) t = 2.06 s, (c) t = 3.2 s, and (d) t = 5.0 s in the second tank.
(a)(b)(c)(d)


## 5. 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 two-phase model to avoid the mesh movement or reconstruction. The Newtonian and non-Newtonian fluids sloshing in the tank with smooth and convex rough bottom have all been studied.The 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. There 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. The good mass conservation property is also demonstrated in the whole flow process.For the same tank, the shake amplitude of the non-Newtonian fluid withn = 0.75 is biggest and it is smallest for the non-Newtonian fluid with n = 1.5. The 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).



---

*Source: 7548192-2020-11-02.xml*