Least Square Finite Element Method for Viscous Splitting of Unsteady Incompressible Navier–Stokes Equations

In order to solve unsteady incompressible Navier–Stokes (N–S) equations, a new stabilized finite element method, called the viscous-splitting least square FEM, is proposed. In the model, the N–S equations are split into diffusive and convective parts in each time step. The diffusive part is discretized by the backward difference method in time and discretized by the standard Galerkin method in space. The convective part is a first-order nonlinear equation. After the linearization of the nonlinear part by Newton’s method, the convective part is also discretized by the backward difference method in time and discretized by least square scheme in space. C0-type element can be used for interpolation of the velocity and pressure in the present model. Driven cavity flow and flow past a circular cylinder are conducted to validate the present model. Numerical results agree with previous numerical results, and the model has high accuracy and can be used to simulate problems with complex geometry.


Introduction
A stabilized method for the solution of incompressible Navier-Stokes (N-S) equations is used in an important research area in finite element. The stabilization of finite element on the basis of solution techniques consists of two components, one is the stabilization of spatial oscillations caused by the discretization of the convection terms at high Reynolds numbers and the other is the pressure stabilization. To circumvent the spatial oscillations caused by the standard central difference type discretization of the convection term, many efforts have been devoted to develop a nonstandard Galerkin finite element method (FEM). Some of the popular and general stabilization methods include Petrov-Galerkin (Xiao et al., 2013), least-squares finite element method (LSFEM) (Bochev, 1997;Jiang, 1992Jiang, , 1998Jiang and Chang, 1990;Pontaza and Reddy, 2006), and characteristic based split (CBS) (Nithiarasu and Liu, 2006;Zienkiewicz and Codina, 1995).
The second source is associated with the spatial discretization of the incompressible N-S equations by the mixed formulations, which restricts the choice of interpolation spaces for the velocity and pressure fields. The incompatible interpolations that violate the LBB condition may induce spurious spatial oscillations in the resulting pressure field. To avoid the restriction in the interpolation approximation spaces, the stabilization strategies for numerical solutions of the problems have also been extensively investigated in the last decades. The interest in using projection methods to overcome the difficulty in time-dependent viscous incompressible flows started in the late 1960s with the ground breaking work of Chorin (1968). The most attractive feature of projection methods is that only one variable needs to solve a sequence of decoupled equations for the velocity and the pressure at each time step, which makes it very efficient for large scale numerical simulations. In the 1980s, second-order accurate versions of projection method have been proposed, such as Goda (1979) and van Kan (1986).
The idea of the LSFEM is to minimize unconstrained convex least squares functional defined as the sum of the governing equations residuals measured in some norms (mostly ). LSFEM has two major advantages over classic- al Galerkin FEM in the numerical solution of the N-S equations. Firstly, LBB (Brezzi and Fortin, 1991) conditions are satisfied naturally, hence the equal order interpolation for the velocity and pressure can be used. Furthermore, the LS-FEM can generate a symmetric, positive-definite algebraic system of equations that can be robustly solved by the conjugate gradient method (Carey and Jiang, 1986). The straightforward application of the LSFEM to the second-order N-S equations requires -type element (Gunzburger and Bochev, 2009). To avoid the practical difficulties in the implementation of such FEMs and -type element can be used, the second-order N-S equations transform to first-order equations by introducing auxiliary variables. The most popular choice is to introduce the vorticity as a new variable (Ozcelikkale and Sert, 2012). The reduced order system of the second-order N-S equations is the first-order vorticity-velocity-pressure formulation (Nickaeen et al., 2014). There is one additional variable in two-dimensional (2D) case and three additional variables are necessary in three-dimensional (3D) case. The second choice is that the secondorder N-S equations are recast to velocity-stress-pressure formulation by introducing the stress as an independent variable (Pontaza and Reddy, 2006). There are three additional variables in 2D case and six additional variables for 3D case. The third choice is that the second-order N-S equations are recast to the velocity-gradient-velocity-pressure formulation by introducing the velocity gradient (Bochev et al., 1999). There are four additional variables in 2D case and nine additional variables for 3D case.
Splitting algorithm, with a flexible format and good stability, has been developed by Yanenko (1971) in the 1970s. Timmermans et al. (1994) solved the convection-diffusion equation by the T-G-based operator-splitting spectral element method, and the convective equation was solved by the T-G method. Thomas and Majda (1981) developed the error estimates for splitting algorithm, which are uniform in the viscosity as it becomes small for either 2D or 3D fluid flows in all the space. Wang et al. (2011) solved the unsteady incompressible N-S equations with the characteristic-based operator-splitting (CBOS) FEM, which combines the viscous-splitting and CBS algorithms. In the model, the simple explicit characteristic temporal discretization, which involves a local Taylor expansion, is referenced from the CBS algorithm and applied to the discretization of the convective part.
C 0 Simple -type element can be adopted to the first-order vorticity-velocity-pressure formulation of N-S equations. However, new auxiliary variables increase the computational cost, especially in 3D case. In order to improve the computational efficiency, the viscous splitting least square (VSLS) FEM, which is based on some references (Bochev, 1997;Wang et al., 2011), is developed to solve the N-S equations. Driven cavity flow and flow past a circular cylinder are conducted to validate the present model. In Section 2, 2D unsteady incompressible N-S equations are introduced. VSLS FEM is explained in detail in Section 3. In Section 4, the solution process of the present model is described. The two numerical examples are studied using the present model in Section 5. Finally, Section 6 draws some conclusions from the present analysis.

Governing equations
The unsteady viscous incompressible N-S equations in dimensionless forms can be written as follows: ∂u (1) ; is the Cartesian coordinates; is the velocity components in the direction; is the pressure; is the time; is the Reynolds number with as the reference velocity, as the reference velocity, and as the kinematic viscosity; is the external forces.

Viscous-splitting technique
Based on the viscous-splitting technique (Thomas and Majda, 1981), Eqs. (1) and (2) can be split into the diffusive part and the convective part where is the solution of Eq.
(3) at the -th step and the initial value of Eq. (4) at the -th step, is the solution of Eq. (4) at the -th time step and the solution of Eqs. (1) and (2) at the -th step.

Finite element solutions for the diffusive part
The temporal discretization of the diffusive part (3) is performed by the backward difference method, and then we have The weak form of Eq. (5) is established by the Galerkin method Γ Ω δu i where is the boundary of , is the virtual displacements. The viscous term of Eq. (6) is integrated by the parts, and then we have According to the finite element method, the computational domain is divided into finite elements and a basic element is denoted by , thus we can obtain Each element is quadrilateral with four nodes. Let , is the number of the velocity nodes of the basic element and . Substituting Eq. (9) into Eq. (8), and taking giving the FEM equation for the basic element (11) is the permutation operator, . Assembling Eq. (10) over the computational domain gives the FEM equation for the diffusive part (3) Then, Eq. (4) is discretized by the backward difference method in time, thus we have Simplifying Eq. (17), we have L is a first-order partial differential operator, and

R R
where is the residual function. LSFEM is based on minimizing the residual function in a least-squares sense. We construct the least-squares functional (Jiang and Chang, 1990) I u δI = 0 Taking the variation of with respect to and setting leads to the least-squares weak statement The FEM equation for the arbitrary finite element unit is where in which T denotes the transpose. Assembling Eq. (25) over the computational domain gives the FEM equation for Eq. (4) K where is the symmetric and positive definite. Therefore, an iterative technique, such as the conjugate gradient method, can be employed to solve Eq. (29).

Solution process
(2) is taken as the initial value of Eq. (4), which is solved to obtain and ; (3) Processing is performed at a time step, and then return to Step (1).
In the solution procedure, the measure, which is used to determine whether the flow has attained its steady state, can be given by the following expression (Burggraf, 1966)

Num n
where is the number of nodes of the entire computational domain and the superscript denotes the time level. at the coordinate origin. The whole computational domain is divided into elements with four nodes and the grids approaching boundaries are gradually refined.

Computational results at different Reynolds numbers
Re= 100 ×129 when Re=100 and Re=1000, and reaches 257×257 when Re=5000. The number of grids in the present model is only 60×60, which is far smaller than that of Ghia et al. (1982).
The streamlines, vortex and pressure contours for the cavity configurations with the increasing Reynolds numbers from 100 to 5000 are shown in Fig. 3. From the computed results it can be seen that the present method can simulate small structure in flow more effectively and the pressure isobars have no oscillation. Fig. 4 shows the distributions of u-velocity along the vertical middle line and v-velocity along the horizontal middle line under different time steps at Re=1000. Fig. 5 shows the convergence of the iterative procedure under different time steps at Re=1000. Fig. 6 shows the distributions of u-velocity along the vertical middle line and v-velocity along the horizontal middle line under different grid numbers at Re=1000. Fig. 7 shows the convergence of the iterative procedure under different grid numbers at Re=1000.

Test of the mass conservation of the present model
If a model exists to solve the 2D incompressible equations, then the results of the model should satisfy the continuity of the fluid. The continuity will provide a very good mathematical benchmark for checking the results, as suggested by Aydin and Fenner (2001) and Erturk (2009). We use the u-velocity and v-velocity profiles in Fig. 2 to test the accuracy of the results. The integration of these velocity profiles will give plus and minus areas as shown in Fig. 2. The degree to which the plus and minus areas cancel each other, thus the integration gives a value close to zero, will be indicative of the mathematical accuracy of the results. The velocity profiles are integrated using the Newton-Cotes rule to obtain the net volumetric flow rate passing through these sections. The obtained volumetric flow rates are then normalized by a characteristic flow rate 0.5 SHUI Qing-xiang et al. China Ocean Eng., 2018, Vol. 32, No. 4, P. 490-500 Aydin and Fenner, 2001), which is the horizontal rate in the absence of the side-walls (Plane Couette flow) to help quantify the errors. The obtained volumetric flow rate values ( and ) are tabulated in Table 1. The volumetric flow rates in Table 1 are so small that they can be considered as . Therefore, the present model obeys the mass conservation.

Laminar uniform flow past a circular cylinder
Uniform flow past a circular cylinder has been widely investigated both experimentally and numerically. The studies have shown that the flow reaches a steady state with the time evolution at and the vortex shedding phenomenon occurs at . To validate the present model, the laminar flow around a circular cylinder is simulated at different Reynolds numbers.

10D 30D × 18D
The domain consists of a cylinder placed at a distance of from the inlet, where D is the diameter of the cylinder. The distance from the center of the cylinder to the top and bottom sides is equal to 9D. The exit of the domain is placed at a distance of 20D from the center of the cylinder. A rectangular flow field of is divided into 14100 four-node finite elements. The total number of nodal points is 14400 (As shown in Fig. 8). At the inflow boundary, the Dirichlet boundary conditions, and , are enforced. The non-slip boundary condition is used on the surface of the circular cylinder and the two side-walls are free-slip. The zero gradient boundary condition is specified at the outlet boundary. As for the pressure, is adopted at the outlet, while is assumed on the other boundaries, where denotes the outward unit normal vector.
The coefficients of the drag and lift are evaluated as: θ where is the azimuth angle measured from the rear point on the horizontal axis of the cylinder and in the clockwise Table 1 Volumetric flow rates Q 2 through a vertical line and Q 2 through a horizontal line passing through the geometric centre of the cavity f where is the dimensional vortex shedding frequency and determined by a spectral analysis (Fast Fourier Transformation Algorithm, FFT) of the time series of the lift coeffi- For the Reynolds number , the flow reaches a steady state with the time evolution. Fig. 6 shows the distribution of the surface pressure at Re=40. In the figure, solid lines denote the results of the present model, dot markers denote the results of Braza et al. (1986), is the azimuth angle measured from the rear point on the horizontal axis of the cylinder and in the clockwise direction, and is the pressure coefficient and defined as:  SHUI Qing-xiang et al. China Ocean Eng., 2018, Vol. 32, No. 4, P. 490-500 495 ρ p 0 where is the fluid density and is the pressure at the front stagnation point. Fig. 9 shows that the present results agree well with the results of Braza et al. (1986). Meanwhile, the pressure on the front of the cylinder is larger than that on the rear of the cylinder. Re = 40 Fig. 10 shows the velocity vector field in the near region of the cylinder at a steady state for . Fig. 11 shows the characteristic dimensions of flow past a circular cylinder at Re=40. Table 2 lists the geometrical and dynamical parameters of the wake of flow past a circular cylinder

496
SHUI Qing-xiang et al. China Ocean Eng., 2018, Vol. 32, No. 4, P. 490-500 is the length of the recirculation region; (a, b) is the location of the recirculation center. We can see that the present results agree well with the reference values.  Table 3 lists the average drag coefficients , and of the present results at Re=100 and 200. Generally, the present results are well within the range of the results reported by other researchers.

Unsteady flow at Re=100 and 200
Figs. 14 shows the development of the pressure and streamline during a cycle at Re=200. The subplots of Figs. 15a-15e correspond to the time instants of the maximum negative, zero, maximum positive, zero, and maximum negative lift coefficients, respectively. As shown in Figs. 14a and 14e, the upper side of the cylinder experiences the highest pressure in the vortex shedding cycle, whereas the lower side of the cylinder is subject to the lowest pressure and has a larger well-developed attached vortex. Hence, the maximum negative lift force on the circular cylinder is observed. Analogously, the maximum positive lift force can be found in Fig. 14c

498
SHUI Qing-xiang et al. China Ocean Eng., 2018, Vol. 32, No. 4, P. 490-500 distribution of the pressure on the upwind surface of the cylinder has little change and is symmetric in one cycle of vortex shedding. Therefore, the pressure on the upwind surface has a limited effect on the formation of the lift force. Table 4 lists the comparison of the computational efficiency by the two models in 2D case. Numerical tests show that the computing speed of the present model has lower computational efficiency than the VVP-LSFEM. The reason is that the present model requires solving equations with 3*Num and another equation with 2*Num in 2D case and the VVP-LSFEM requires solving only one equation with 4*Num in the same case. However, the VVP-LSFEM requires solving only one equation with 7*Num in 3D case and the present model requires solving equation with 4*Num and another equation with 3*Num in the same case. Therefore, the present model would have higher computational efficiency than the VVP-LSFEM in 3D case.

Conclusions
A new stabilized FEM, called the viscous-splitting least square (VSLS) FEM, has been developed to solve the N-S equations. In each time step, the N-S equations are split into the diffusive and convective parts, which respectively consider the convection dominated characteristics and the diffusive feature by a viscous-splitting algorithm. The standard Galerkin method is adopted to solve the diffusive part, and its results are used as the initial values for the convective part. The convective part is a first-order nonlinear system, so it can be solved directly by the LSFEM without introducing auxiliary variables. The convective part is discretized by the backward difference method in time and discretized by the least square scheme in space after the linearization of the nonlinear part by Newton's method. The resulting matrix is symmetric, positive definite, and solved by the conjugate gradient method. No upwinding or adjustable C 0 parameters are contained in the present model. Compared with the LSFEM by introducing auxiliary variables, given boundary condition would be easier. And moreover,type element can be used for interpolation of the velocity and pressure in the present model.

Re ⩽ 5000
The lid-driven cavity flow problem and the flow past a circular cylinder are adopted to validate the present model. The conclusions are as follows: (1) Solutions of the lid-driven cavity flow are obtained at different Reynolds numbers. At , the results agree well with the available benchmark solutions in the related literature. (2) Through the numerical simulation of the flow past a circular cylinder, it can be seen that the present model can be used to simulate problems with complex geometry.