Understanding The Navier-Stokes Solvers

The Navier-Stokes equation explains the motion of viscous fluid using partial differential equations. There is no exact solution for Navier-Stokes equation due to the non-linear set of partial differential equations. The Navier-Stokes equations use the momentum and continuity equations expressed in equation 1.

The first equation of equation 1 states that the continuity does not change and that mass is conserved. The second equation of equation 1 is the momentum equation. The momentum equation is Newton’s second law where mass times acceleration equals forces. The momentum equation has been adapted to flow by including viscosity, pressure and velocity force related to time. Combining the momentum and continuity equations forms non-linear Navier-Stokes equation 2 to solve flow in a viscous, incompressible fluid.

There are many methods to solve Navier-Stokes equations such as the Newton method and Chorin method. The two methods, Newton and Chorin, are used to solve different problems such as steady solver, solutions that are not time dependent and unsteady solver, solutions that change over time. Note that the methods mentioned can be adapted to be a steady or an unsteady solver with the additional time dependent terms, as explored later, however may perform worse to an unsteady solver counterpart. Also, this guide is a simplistic review of different solvers for Navier-Stokes equations and the associated solvers where the purpose is to highlight the solver choice for solving CFD microfluidic problems. For in-depth detail regarding the Navier-Stokes equations used, consult Wikipedia, YouTube (personally I like Fluid Mechanics 101) and all good resources.

Often the decision to use either solver is based on computational resources as well as whether the problem has turbulence. Turbulence can be determined by calculating the Reynolds number, equation 3.

Steady Solver

How Navier-Stokes Newton Method Works

The Newton method aims to solve the Navier-Stokes equation in the non-linear form by investigating iterations of the equation equalling 0, equation 4. The iteration Newton method attempts to solve the velocity by adjusting the velocity. Through the iterations, the solutions should begin to converge upon the accurate velocity value. The Newton method for solving Navier-Stokes is often referred to as a steady solver due to the iterations towards a final solution and the lack of time-dependent variables compared to unsteady solvers. Steady solvers’ lack of time-dependent variables means they are poor at simulating turbulence.

The implementation of the Navier-Stokes Newton method was adapted from the FEniCS demo – Stokes equations with an iterative sovler and Oasis default NSCoupled solver.

Unsteady Solver

How Navier-Stokes Chorin Method Works

Navier-Stokes equation can be adjusted to include a time-dependent variable to describe the fluid, equation 2. The splitting solver works by separating the flow solution using pressure, p, and velocity, u. One of the oldest methods was proposed by Chorin (1968) and Temam (1969), often referred to as Chorin’s method. Sometimes this method is called the “projection method”, equation 5.

The Chorin method works by decoupling the velocity and pressure fields. This is achieved by calculating a tentative velocity, correct for pressure and then correct for velocity to solve the Navier-Stokes equation. The Chorin method solves turbulent flow over time by iterating from a starting assumption based on inlet and outlet domains.

The implementation of the Navier-Stokes Chorin method was adapted from the FEniCS demo – Incompressible Navier-Stokes equations.

How Navier-Stokes Newton Method + Time Dependent Variables Works

As the name suggests, the Newton method, equation 4, has an additional implementation of time-dependent variables, equation 6, based upon equation 2 describing the fluid force at a given time. This method produces an unsteady solver and can speed up computation when compared to the Chorin “splitting method” by attempting to solve both the velocity and pressure at the same time. However, this method can produce errors where a solution cannot be converged upon and an initial solution never arrived upon. This is due to the linearised inertial term which can fail to make the subproblem linear.

Problem Specific Variables

Inflow velocity, boundary domains and other flow properties set up the problem variables.

Parameters	Value
Water Dynamic Viscosity  (kg / m.s)	1.003e-3
Water Density  (kg / m3)	1
Water Velocity (m / s)	1.5
Time Step (s)	1e-4

One difference between using the steady and unsteady was enabling “iterative_solver”. Steady and unsteady solvers for the Navier-Stokes equations in FEniCS may require solver methods optimisation. Solver methods differ between linear and non-linear solvers where different methods, pre-conditioners and tolerances (relative and absolute) need to be tailored to the specific problem. It is recommended to run a direct solver, such as MUMPS, before investigating iterative solvers to test steady solver solutions. Unsteady solvers with iterative solvers are more forgiving in solving the Navier-Stokes equations and generally can be solved with iterative solvers default settings.

Software

BERNAISE, Binary ElectRohydrodyNAmIc SolvEr, is an open-source, free CFD solver created in FEniCS – a solver for partial differential equations (PDEs). BERNAISE is a flexible high-level solver of electrohydrodynamic flows in complex geometries. It is written in Python and built on the FEniCS project, which in turn effectively interfaces to optimized linear algebra backends such as PETSc. The solver is described in Asger Bolet’s, Gaute Linga’s and Joachim Mathiesen’s paper.

This particular version of BERNAISE has been modified by myself to use .XDMF computer aided design (CAD) meshes as well as has been focused towards microfluidic examples. Due to the large size of the simulation, additional computational power may be required.

Step By Step

1. A CAD file can be used as well as a mesh generation by FEniCS / Dolfin. In this example, mesh is generated with FEniCS Dolfin. Documentation can be found here for FEniCS / Dolfin mesh generation. A full description of 2D and 3D CAD conversion is available on the BERNAISE wiki.
2. Assign parameters in the problem python file. In this case, the problem has been called flow_around_cylinder and the problem specific parameters are defined in the problem() function. More example problems can be found in the problems folder. In particular, the boundary labels mentioned in step 1 for assigning boundary properties as well as the Problem Specific Variables mentioned earlier in this blog. It is essential the walls, inlet and outlet are set correctly to solve the problem. Also, depending on which solver was tested, the unsteady solver "basic" and / or "basic_IPCS_Adj" was switched to a steady solver "basicnewton" and vice-versa.
3. Run the problem from the command prompt python sauce.py problem=flow_around_cylinder
4. View the exported .XDMF results in ParaView.

Reviewing The Results

BERNAISE has presented a platform that has implemented steady and unsteady solvers. The difference in the two solutions can be seen when simulating turbulent flow. The solvers were tested with the flow around the cylinder. As expected, turbulence was not shown in the steady solver due to the solution convergence. In comparison, the unsteady solver presented turbulence evolving over time.

Turbulence is not expected due to the Reynolds number below 1000. The Reynolds number is expected to be around 613 ((1 kg / m³ x 1.5 m/s x 0.41 m) / (0.001003 kg / m.s), equation 3.

Steady Solver

The steady solver, Newton method, converged on a solution within 5 seconds and the solution remained the same over time due to the lack of time-dependent variables.

Streamlines and velocity flow presents sharp drop offs from the tail-end of the flow around the cylinder. This is expected due to the flow converging equally around the cylinder but is unrealistic due to the viscosity and expected turbulence created by the cylinder inflow. This can be seen by looking at the differences in the velocity within the tail-end which increases inflow from 0 m/s to approximately 1 m/s before returning to 0 m/s. The changes in flow velocity are unrealistic and suggest turbulent flow. Reviewing the Reynolds equation, equation 3, a laminar flow result would be expected if the channel and cylinder were smaller, flow velocity decreased, fluid density decreased as well as fluid viscosity increased.

Unsteady Solver

The unsteady solver, the Chorin method, calculated the evolution of turbulence over time. Similar to the steady solver, the solution developed a tail after the cylinder initially. However, the solution became turbulent where flow changed preference from above and below the cylinder creating ripples. The colour scale, blue to red indicates the velocity of flow from 0 m/s to 2.4 m/s where the flow was shown to oscillate from higher and lower velocities due to the pressure change over time.

The changes in turbulence are useful for mixing applications or understanding the displacement caused inflow. Reviewing the Reynolds equation, equation 3, the turbulence presented by the unsteady solver can be adjusted by changing the channel and cylinder were smaller, flow velocity decreased, fluid density decreased as well as fluid viscosity increased.

Affecting the Reynolds number

By adjusting the factors affecting the Reynolds number, equation 3, the flow dynamics change dramatically. As the focus of the work going forward is microfluidic related, the model has been shrunk by 1000 times to the micro-scale.

As shown below, no turbulence is measured over time using the unsteady solver with the micro-scaled flow around the cylinder model. This is due to the micrometre linear length, from L = 0.41 m to L = 0.00041 m, leading to little to no turbulence. This is indicated by the new Reynolds number of 0.613 compared to 613 prior to the reduced dimension, equation 3. This suggests that most microfluidic CFD solutions can be solved using a steady solver without time-dependent variables due to the lack of turbulence.

Conclusion

Flow dynamics such as turbulence are presented differently using the unsteady and steady solvers as demonstrated with the flow around the cylinder model. It is fundamental the correct solver is used for analysing the result produced. Often calculating the Reynolds number is an important determination whether the solution requires an unsteady solver. Analysis of a system should take into account these assumptions when using either solver.

Generally, the unsteady solver, Chorin method, is best to visualise laminar and turbulent dynamics within an unknown system accounting for time-dependent and steady states. However, if the dynamics are not time-dependent, a steady state system is preferred due to converging on a solution and, in theory, are less computationally demanding. A steady solver may be optimal when investigating microfluidic models in the future where flow dynamic viscosity dominates over flow inertial forces.

References

1. Ladyzhenskaya, O. A. (1975). Mathematical analysis of navier-stokes equations for incompressible liquids. Annual Review Fluid Mechanics, v, 249–272. https://doi.org/10.1146/annurev.fl.07.010175.001341
2. Fluigent (2021) Basic Properties of Microfluidic Flows. [Accessed on 25th August 2021] Physics of Microfluidics – Basic Properties of Microfluidic Flows | Fluigent
3. C.T. Shin U. Ghia, K.N. Ghia. (1982) High-Resolution for incompressible flow using the Navier-Stokes equations and the multigrid method. J. Comput. Phys., 48:387–411.
4. FreeFEM (2021) Newton Method for the Steady Navier-Stokes Equations. [Accessed on 25th August 2021] Newton Method for the Steady Navier-Stokes equations (freefem.org)
5. A. J. Chorin (1967) A numerical method for solving incompressible viscous flow problems. J. Computational Physics, v. 2, p. 12.
6. Linga, G., Bolet, A., & Mathiesen, J. (2019). Bernaise: A flexible framework for simulating two-phase electrohydrodynamic flows in complex domains. Frontiers in Physics, 7(MAR). https://doi.org/10.3389/fphy.2019.00021

Related Posts