Simulation of incompressible viscous flow using finite element method

. This study focused on simulating incompressible viscous flow using the finite element method. This study used velocity and pressure as unknowns known as primitive variable formulations. Simulation of incompressible fluid flow poses numerical challenges due to the presence of nonlinear convective terms in Navier-Stokes equations and the incompressible nature of the fluid. If the connection between velocities and pressure is not discretized correctly, the stable and convergent velocities might be gained, but the obtained pressure will be oscillatory. To avoid these difficulties, continuous quadratic and additional cubic bubble functions will be used for the velocity field and linear functions for the pressure field. This kind of discretization satisfies the Ladyzhenskaya-Babuška -Brezzi (LBB) stability condition. Two cases of different Reynolds numbers were used to test the formulation's effectiveness. In the case of Reynolds number 0.12, no vortices were formed, suggesting that the flow is primarily governed by fluid friction, and fluid inertia has minimal effect. In the case of Reynolds number 120, the vortex formation, which is known as Von Kármán vortex street, appeared. These results concluded that the formulation using the finite element method is correct.


Introduction
The study of incompressible fluid flow plays an important role in many applications of civil engineering such as the distribution of pollutants, water quality flow, floodplain inundation, and flowing streams of water in sea waves.Incompressible fluid flows can be classified based on various criteria such as internal or external flows, laminar or turbulent flows, and single-phase or multiphase flows [1].Almost every type of flow is governed by the same equations, namely the Navier-Stokes equations [2].Incompressible fluid flow can be studied in three approaches, namely experimental, analytical, and numerical approaches.This work employed a numerical method to study incompressible flow.Computational simulations are pivotal in understanding complex incompressible fluid flow phenomena since the simulations can mimic real and ideal conditions and give comprehensive information at a relatively low cost in a short period.Among various numerical methods such as finite difference and finite volume methods, the finite element method (FEM) has proven to be a valuable tool, since FEM can easily handle complex regions and complicated boundary conditions in real-world problems.The FEM is a generalization of the classical Ritz and Galerkin weighted-residual methods.However, the FEM method for analyzing incompressible fluid flows has some challenges, which do not occur in analyzing solid mechanics problems.In solid mechanics, the classical Ritz and Galerkin methods are the same as the problem of minimizing an unconstrained convex quadratic functional.Meanwhile, Galerkin weightedresidual formulations of Navier-Stokes equations are not equivalent to the minimization of unconstrained convex quadratic functionals.Other challenges of the simulation of incompressible fluid flows are the presence of nonlinear convective terms in the Navier-Stokes equations and the incompressibility of the fluid.These problems lead to stability issues, i.e. the convective terms will create spurious oscillations (wiggles) solutions.Therefore, the standard Galerkin formulations have to be modified.
This study used the Galerkin mixed method, in that velocities and pressure variables were approximated as fundamental unknowns.The word mixed indicates that velocities variables are used together with the force-like variable, pressure, in a single formulation.Through this formulation, the original minimization problem converts to a saddle-point problem.The sufficient condition for a saddle-point problem to have a unique solution is Ladyzhenskaya-Babuška-Brezzi (LBB) condition.The LBB condition guarantees that the discretization of a saddle point problem is stable.To satisfy the LBB condition, this study used higher-order velocities elements and lower-order pressure elements.
Since the Galerkin mixed method is not robust for numerical modeling of incompressible viscous flow past the rigid surface of an immersed body or along a fixed impermeable wall.To make the Galerkin method robust, bubble functions are introduced.
Bubble functions are functions located on the elements' interiors and vanish on the element boundaries.With the help of the bubble function, the solution to the boundary value problem will be broken down into the sum of the coarse-scale and fine-scale solutions.The classical Galerkin finite element will obtain the coarse-scale solution, meanwhile, the bubble function will take care of the fine-scale solution.
This study will employ the finite element code which was developed by Matuttis and his coworker of the Department of Mechanical and Control Engineering, University of Electro-Communications (UEC), Tokyo [3][4][5][6][7].The code employed a second-order polynomial for velocities and a linear function for pressure.Taking into account the bubble function the second-order element becomes the third-order element [7].This development eliminates the need for a genuine third-order element, streamlining the computational process while ensuring LBB stability conditions.
To demonstrate the formulations' effectiveness, the code will be used to simulate the incompressible viscous flow of different velocities past a sphere.This viscous flow will create a thin boundary layer adjacent to the sphere's surface; then a trailing wake will exist in the flow behind the cylinder.
Specifically, the focus lies on simulating a twodimensional Von Kármán vortex street, a well-known flow phenomenon characterized by the formation of alternating vortices in the wake of a body.This phenomenon has been extensively studied and validated in various flow situations [3].The Von Kármán vortex street will be simulated under different flow conditions, i.e. by varying Reynolds numbers (Re).

Underlying basic equations
Two-dimensional flow can be modeled under the following assumptions: 1.One of the dimensions (z-direction) is very long, and no flow exists along that direction.

The velocities in the other two directions do not vary
with the z-directions.With these assumptions, the governing basic equation for such incompressible flow can be summarized as follows.Comprehensive treatment can be found in standard continuum mechanics books such as [8][9][10].

Mass equation
The conservation of mass can be written as Equation 1.
where  is the fluid density,   and   are the fluid velocities in the -direction and -direction, respectively.In the case of incompressible flow, the conservation of mass Equation 1 reduces to Equation 2.
Equation 2 states that no volume change during the deformation of incompressible fluid.

Momentum equation
Equations 3-4 shows the conservation of linear momentum.
where   ,   and   are the Cartesian components of the total stress tensor , while   and   are the body forces in the -direction and -direction, respectively.

Constitutive equation
Equations 5-6 describe the constitutive equations for viscous fluid.
where  is the viscosity,  is the pressure, and   ,   and   are the Cartesian components of the viscous stress tensor .

Galerkin formulation
The discretization of the Navier-Stokes equations will be done by formulating the weak form of Equations The functions   and   are called the shape functions.
The importance of FEM is that it provides a continuous solution function over the whole computational domain.This is achieved by associating a shape function with each node of each element.In fluid dynamics, selection of the shape function is constrained by a condition, namely the Ladyzhenskaya-Babuška-Brezzi (LBB) stability condition.The LBB states that the order of elements for the velocities must be one order higher than for the pressures, as the Navier-Stokes equations also contains the velocities must be one order higher than for the pressures [5].
Gresho and Sani [11] presented options of different order of polynomials to be used as shape functions, some of which are shown in Fig. 1.To take the LBB condition into account, third order elements are necessary.However, increasing the number of elements between particles would increase the computational effort considerable.The P2+ combines the second order P2 element and the bubble function of the third order P3 element.Thus, using the P2+ element for the velocity and P1 element for the pressure (P2+ P1 element) still guarantees the LBB stability without the need for a full third order P3 element.
The barycentric coordinates of point  can be obtained by drawing a line connecting  and each of the triangle's corner points , , and .Each barycentric coordinate corresponds to the fraction of area of the sub-triangle (Equation 15).Since all sub-triangles combined is the area of the whole triangle, the sum of all barycentric coordinates is   +   +   = 1 (Fig. 2).

Methodology
The simulation is done using the code that was built on MATLAB programming environment.The code used in this research is a modified version of the original FEM code by Matuttis Laboratory, University of Electro-Communications, Tokyo.

Geometry and mesh generation
Two channels with different geometries are configured for the simulation.The first channel is configured as a microfluidic channel with a length of 8 mm, a width of 1 mm.A particle with a diameter of 0.1 mm is placed in the centerline of the channel acting as an obstacle to the flow.The second channel is similarly constructed but with bigger scales, which are the length of 500 mm, width of 200 mm and 100 mm diameter particle.The characteristic length  ℎ of the system is based on the diameter of the particle.Therefore, the first channel has  ℎ = 0.1  while the second channel has  ℎ = 100 .The geometry is discretized using a graded mesh technique (Fig. 3) to create a computational grid for the simulation [8].The mesh is refined in the region near the particle and to the rear of the particle where vortices are expected to appear.The use of graded mesh allows to capture the flow features accurately while giving reasonable computing time.

Fluid properties
The fluid used in the simulation is water.The properties of water, such as density and viscosity, are specified to accurately model the flow behavior.Fluid density is set to  = 1000 / 3 and the dynamic viscosity is set to  = 10 −3  ⋅ .

Numerical solver
The Newton-Raphson method is employed as the numerical solver to solve the discretized Navier-Stokes equations.The Newton-Raphson method is an iterative technique that allows for the efficient solution of nonlinear systems of equations.Since the Navier-Stokes equations are not possible to solve analytically, an iterative approach to calculate the flow field at each timestep is needed.This is possible with the Newton-Raphson method, an iterative method to calculate the roots of equations where an initial guess  0 is refined through iterations  until it approaches (  ) = 0 +  with small enough residual  (Equation 16) [12].

Simulation parameter
The simulation is started with a low Reynolds number  (Equation 17) and is increased over time until it reached maximum   .This is done to ensure that the onset of the vortex is produced properly.The adjustment of  is done by changing the prescribed inflow velocity  0 (Table 1).
Table 1.Variation of inflow velocity  0 over time.
The simulation is performed with specific time step sizes and total simulation time (Equation 18).The time step size is computed via time-adaptive backward differentiation formula of second order (BDF2) predictorcorrector [6,10].To control the timestep size , the deviation between   and   .
The indices  − 1 and  + 1 denote the previous and new timestep size.The maximum allowed timestep in the simulation is set to   = 0.02 seconds.The total simulation time is limited to 120 seconds.

Results
After the simulation, post-processing and analysis of the results are conducted.The flow field variables, such as velocity and pressure, are extracted from the simulation data.Quiver plots are used to visualize the velocity field.While the vortices can be seen from the arrows of quiver plot (Fig. 4), the vorticity can be better visualized by calculating the curl which then is plotted as a colormap (Fig. 5).Trail of vortices is apparent in the larger channel with  = 120.The formation of vortices indicates that the flow is in the transient state where the inertia of fluid has substantial effect to the flow.The obstacle in the channel caused disturbance to the flow which initiated the formation of vortices.Note that while there are vortices, the flow is still in the laminar region which is confirmed by the shapes of vortices that are relatively regular shaped and not turbulent.
On the contrary, there are no vortices formed in small channel with  = 0.12 although the flow is also being obstructed by a particle.The low Reynolds number ( < 1) indicates that the flow is largely govern by the fluid friction and the fluid inertia has minimal to no effect.This means that even when the flow is disturbed by the obstacle, the fluid inertia is not strong enough to start the formation of vortices.

Conclusion
In conclusion, the simulation results have successfully verified the accuracy of the developed Finite Element Method (FEM) code for incompressible fluid flow.The presence of vortices in the channel with a Reynolds number (Re) of 120 indicates the influence of fluid inertia on the transient flow state, while the absence of vortices in the small channel with a Re of 0.12 highlights the dominance of fluid friction.These findings align well with existing observations, validating the FEM code's capability to effectively model complex flow phenomena.

Fig. 2 .
Fig. 2. Velocity field of the fluid flow near the particle at maximum  for each channel.

Fig. 3 .
Fig. 3. Vorticity field of the fluid flow showing trails of vortices.
0 is imposed to represent the incoming flow.At the outlet, a pressure boundary condition is prescribed.The top and bottom walls of the channel are considered as no-slip boundaries, where the velocity is set to zero.
2.5 Closure of the initial boundary value problemAppropriate boundary and initial conditions should be prescribed to close the problem.Here the primitive variables are velocities and pressure and the Dirichlet boundary conditions are given in the means that a prescribed velocity profile 2 and 7-8.To do that, Equation 7 is multiplied by the test function  1 , then is integrated over the element Ω  .The integration by parts gives Equations 9-11.