Numerical Solution of the Navier-Stokes Equations for Incompressible Fluid Flow by Crank-Nicolson Implicit Scheme

The Navier-Stokes (N-S) equations for incompressible fluid flow comprise of a system of four nonlinear equations with five flow fields such as pressure P, density ρ and three velocity components u, v, and w. The system of equations is generally complex due to the fact that it is nonlinear and a mixture of the three classes of partial differential equations (PDEs) each with distinct solution methods. The N-S equations fully describe the unsteady fluid flow behaviour of laminar and turbulent types. Previous studies have shown existence of general solutions of fluid flow models but little has been done on numerical solution for velocity of flow in N-S equation of incompressible fluid flow by Crank-Nicolson implicit scheme. In practice, real fluid flows are compressible due to the inevitable variations in density caused by temperature changes and other physical factors. Numerical approximations of the general system of Navier-Stokes equations were made to develop numerical solution model for incompressible fluid flow. Adequate solutions of the latter produce numerical solutions applicable in numerical simulation of fluid flows useful in engineering and science. Non-dimensionalization of variables involved was done. Crank-Nicolson (C.N) implicit scheme was implemented to discretize partial derivatives and appropriate approximation made at the boundaries yielded a linear system of N-S equations model. The linear numerical system was then expressed in matrix form for computation of velocity field by Computational fluid dynamics (CFD) approach using MATLAB software. Numerical results for velocity field in two dimensional space, u(x,y,t) and v(x,y,t) generated in uniform 32×32 grids points of the square flow domains, 0≤x≤1.0 and 0≤y≤1.0 were presented in three dimensional figures. Results showed that the velocity in two dimensional space does not change suddenly for any change in spatial levels, x and y. Therefore, C-N implicit Scheme applied to solve the N-S equations for fluid flow is consistent.


Introduction
The system of Navier-Stokes (N-S) equations for incompressible fluid flow fully describes the dynamics of Newtonian fluid. It is basically the statement of the principles of conservation of momentum of fluid flow in differential and integral forms. It is in essence the Newton's Second law of motion of fluids as by a system on nonlinear PDEs. Conventionally, each of the three classes of second order PDE namely, elliptic, parabolic and hyperbolic have distinct numerical solution methods. The general form of N-S equation is a complex system due to the fact that it consists of nonlinear PDEs in 3D and a mixture of the three classes of second order PDEs. The unsteady viscous term renders them parabolic type, the diffusion term make them hyperbolic type and the source or sink of forces render them elliptic (Poisson's equation). Pressure term coupled with velocity field for instance, is a function of space and time and determined as part of general solution of the incompressible momentum equation. The mass conservation equations can easily and readily be solved by any of the methods applicable to differential equations. Therefore, solution of the system of momentum conservation equations was the goal.
Numerical solutions for a desired fluid flow problem require modification of the general system, numerical approximation and application of appropriate initial and boundary condition. This develops numerical solution model of flow fluid for simulation of variety of flow situations of the problem.
No fluid can be considered perfectly incompressible even in carefully controlled flow situations since minute variations in density occur due to inevitable variations in temperature and many other factors [19]. However, the incompressible fluid flow model is accurate and necessary approximation in the sense that it contains certain distinctive mathematical and physical properties (eg. pressure, density, velocity) suitable and applicable for modeling a variety of flow situations that require special computational techniques. Taking heat conduction in the fluid flow to be constant, the following general form N-S equations can fully describe the motion of a viscous incompressible fluid flows in some bounded domain Ωϵℛ , n = 2,3.
Where, the three velocity components in Cartesian coordinates systemV = , , , gradient operator ∇= + + ! and the Laplace operator ∇ = " " + " " . The Pressure gradients in the system provide the degree of freedom needed to ensure that the condition of continuity equation is satisfied by the flow velocity. Therefore, pressure field coupled with the equation of continuity is determined as part of the general solution of the momentum equation of the incompressible fluid flow [2].
Practically, lamina (smooth and orderly) flow if any occur only a small section of the flow domain of many boundary layer flows. Fluid flows past a surface exhibit considerable turbulence (disorderly flow) [16]. Reynolds number, Re being the ratio of inertial force to viscous force determines the nature of flow. Certain range of Re gives the transition from lamina to turbulent flow types.
Generally, the solutions for a specific practical flow problem desired to model a flow phenomenon require rigorous, relevant or reasonable approximations of general system of fluid flow equations to arrive at flow models that satisfy the conditions of the desired flow situation.
The solutions for Burger's equation which is model of the incompressible N-S equation in 1D without equations of continuity and pressure functions using modified C-N scheme, showed consistency of the numerical solution method [7,8,18]. Comparing the numerical results for the of the problem with those of pure C-N implicit scheme applying uniform grid with constant Re at 4000 and very small time stepping (∆t) at 0.001 [8], showed better consistency and rate of convergence in the modified C-N scheme than pure C-N scheme. Though, the model problem did not focus on the principles of conservation.
Linearization done for the Burger's Equation by Newton's Methods and numerical solution by Gauss's elimination with partial pivoting achieved convergent numerical results of velocity field [18]. Similarly, results from the solutions for viscosity and velocity observed values close to those from exact solution were [7].
Numerical solution for time-dependent Navier-Stokes models in curvilinear coordinates with its three velocity and vorticity components discretized by second order central difference in space and third order semi-implicit Runge-Kutta in time produced highly efficient, accurate and consistent results [11].
The C-N implicit scheme combined with time efficient ADI combined using iterative methods produces even more efficient numerical solution model [15]. The central difference approximations of the C-N implicit scheme applied on initial value problems (IVP) of the quasi linear parabolic PDEs and hyperbolic PDEs produce more accurate numerical results than forward difference approximations with decreased steps sizes [4,17].
Initial conditions and Von Neumann boundary conditions must be applied [8] to develop a complete and accurate numerical solution model. Neumann boundary conditions were preferably implemented in. Numerical solution of particularly the flow problem because it produces accurate and reliable numerical results [14].
Numerical solution for N-S equations with the Dirichlet boundary conditions exist in a continuous flow domain (well-posedness) and the solution converge to suitable weak solutions [5]. The Results for velocity-pressure solution solutions %, & for the 3D N-S equations in a bounded domain '∶=Ω × 0, * proved possible extensions to numerical approximation by C-N schemes. But lack of regularity in & made it difficult to prove global existence of solutions.

Incompressible Fluid Flow Fields
Fluids in motion are deformed, translated from one position to another and rotated at different magnitude and direction during the flow. The nonlinear PDEs of fluid flow problem are arrived at by the coupling of the primitive variables involved in the practical flow situation.

Velocity Gradient Tensor
Viscous fluids flow is basically caused by resistance to shear stresses. The translational motion of the fluid elements produces velocity gradient tensor, a second order tensor, given as 3 × 3 obtained as follows. The diagonal elements in the velocity gradient tensor denote the rate at which the fluid flow may undergo viscous strain during flow while the off-diagonal components show changes in the position of fluid elements in the directions normal to the direction of flow velocity components (shear motions). The rate of deformation tensor is equal to the velocity gradients and its transpose [1].
In this case, the off-diagonal elements are symmetrical (equal) and decomposed as follows; The extreme right parenthetical parts denote rotation of the fluid during flow with the angular velocity vector ω. These produce vorticity due to Coriolis forces which in essence cause deflections on the fluid elements in the direction normal to the direction of entire fluid flow. They yield pure shear, rotation or even a combination of both depending on the flow situation.
Therefore, the fluid in motion possesses four physical features of flow velocity namely: pure translation, vorticity, dilation (or contraction) and symmetric rate-of-shear tensor. Translation is specified by velocity components, (u, v, w) in some time t while the rate of circulation which gives vorticity specified by twice angular velocity components ω , ω and ω ! given by; This results in viscous stress which is linearly proportional to (twice) the strain rate and vorticity, ζ is twice the angular velocity ω [3]. The rate of strain of fluid flowing is therefore a second order tensor.
The 3 × 3 matrix can further be expanded as follows.
The Superposition of the four features of fluid flow, one or more of which may be zero at any point of the flow determine the fluid flow behavior.
Equation (6) can be decomposed further to obtain the sum of the rate -of-strain tensor (symmetric) and the rotation tensor (skew-symmetric) as follows.

Inertial Force on the Fluid Flow
For a fixed control volume (constant unit dimensions) of incompressible fluid undergoing translational motion in equation (2), the control mass is constant during flow since there is neither source nor sink of mass during the flow. The respective transition of mass flux components in the x, y and z directions during fluid flow in equation (5) The change in mass of the fluid in the control volume in motion during flow in the direction of x parallel to the dxdy plane, for instance, is given by; Therefore, total change in mass flux of fluid flow through the control volume is zero. There is neither sink nor source of mass during the flow of incompressible fluid. Therefore, the change in mass is given by; The unsteady change in mass flux in each of the corresponding component direction during the fluid flow is due to dilation or contraction of the fluid. This is in essence material derivatives (derivatives of mass flux by chain rule) of mass flux in each component directions. They are given as; The equation for steady incompressible fluid flow, -V 2 = 0 when ρ is constant, is the divergence-free in the velocity vector fields given by: This yields the equation of continuity of flow. Similarly, the change in momentum fluxes in x, y and z directions during translational motion given by; The change in momentum of the fluid in the x direction parallel to the dxdy plane of control volume is given by: Similar expressions for changes can be obtained for pressure gradient and the viscous stress components [10]. The change momenta of fluid flow along the directions of x, y and z past their respective orthogonal planes of the control volume become; Therefore, the respective rate of change in momenta of the fluid flow after a time t gives the respective acceleration field, a of the flow. They comprise of the steady and advection terms denoting the inertial effect in the fluid flow equation. They are obtained by material derivatives of the momentum flux components given by; The system yields momentum equations which basically signify the inertial forces per unit area on control surface. The corresponding components give the inertial forces when multiplied by density (mass per unit volume) given as;

Viscous Forces and Dynamical Balances
The vorticity term described in equation (5) produces considerable viscous effects (forces) due to fluid strain and circulation. The result is a second order tensor (specified by surfaces and forces), viscous stress tensor b given by; For incompressible fluid, the normal stresses, c on control surfaces are isotropic (equal independent of direction), σ = σ = σ ! and each is given by: In particular, viscous stress of fluid flux normal to the x direction through the orthogonal plane parallel to the xz plane σ for instance, is given by; Where, μ is coefficient of dynamic viscosity and η is the coefficient of elastic viscosity of the fluid in motion.
This produces normal force on the control surface when multiplied by the control surface area of the corresponding orthogonal plane. Thus, the corresponding normal surface forces exerted on the control surfaces of the fluid in the x, yand z directions are given by; The average of the equal normal stresses for a continuous flow of incompressible fluid is given by; The equation of continuity for incompressible fluid obtained in (12) due to pressure effects implies that the average normal viscous stress is equal to negative pressure.
It then follows that the Stokes hypothesis condition of Newton's viscosity assumption for Newtonian fluid holds [12]. This is given by; Similarly, the symmetrical shear stresses in (2.8) are given as follows.
They are directly opposite off-diagonal components (shear stress) which act to deform the fluid diagonally (shear strain). They yield corresponding viscous forces in the fluid in motion given by; Viscous fluid flows at high Reynolds numbers Re (the ratio of inertial force to viscous force) depicts quick variation of vertical velocity u(y) (velocity varying with depth y from solid surface). The stress tensor, maps linearly the surface force on the surface of the deformable continuum to the unit normal on the surface. Therefore, the resultant of viscous (surface) forces is the algebraic sum of forces due to the normal and shear stresses given by; Thus, the total surface force due to the normal and shear stresses shown in the system of equations (27) is given by;

Methodology
The procedures incurred solution of the problem were as follows.

Formulation of Conservative Fluid Flow Equations
The conservative system of Navier-Stokes equations for incompressible fluid flow constitutes system of fluid equations that satisfy the conservation principles such of mass and momentum during fluid flow.
The significant viscous forces together with body forces cause imbalance between forces in the fluid resulting fluid motion (inertial effects). Dynamical balance exist when the net force due to algebraic sum of viscous forces is equal to inertial force (acceleration of the fluid). This is the Newton's law of motion of the fluid [6]. This is given by a system of flow equations in 3D as follows; Boussinessq approximation implemented accordingly develops the desired flow model without necessarily changing the conservation laws [3]. The density is treated as constant except in the gravity term of body force. The Reynolds number Re becomes the key parameter that determines the nature of the flow.

Coriolis Effects and Dynamical Balance
The Coriolis force due to rotation of the fully immersed sphere in the direction of flow velocity deflects the fluid in the direction perpendicular to the direction of rotation. This yields angular acceleration.
Necessary approximation is achieved by assuming geostrophic balance; a balance between Coriolis forces and the pressure gradient The momentum conservation equations for the boundary layer of fluid flow past a rotating earth forms the following system of 2D equation of fluid motion.
Where, ν = r V is the kinetic viscosity.
With the addition of the continuity equation, this system of equation comprise of three equations with three dependent variables, u, v and P.

Conservative Fluid Flow Equations
Conservative systems of Navier-Stokes equations of incompressible fluid flow are basically the statements of principles of conservation of mass and momentum of fluids in integrals forms and converting into differential forms.
By Gauss's divergence theorem, the volume integral of the divergence of mass flux, ρu equals the surface integral of the mass flux normal to control surface [6]. This is given by; where, the minus sign indicates that • is an outward normal and where, S u, t is the algebraic sum (sources or sinks) of forces that produce inertial force.
The fluids circulating at a velocity u about the curve C of the rotating sphere is such that the rate of fluid circulation provide vorticity which is twice the angular velocity ƒ= curl "3…. Thus, the momentum conservation equation in x direction in equation (29)  The conservative equation of motion is given by; Where, the u = u x, y, z is the velocity component in x direction.
Pressure field in the entire flow domain adjust instantaneously to the velocity perturbations during circulation. It results from effects of normal viscous stresses in the momentum conservation equations. The divergence of the velocity field deviates from zero when the fluids circulate around a curved surface.
Pressure distribution is altered by the net force due to horizontally rotating sphere. Large pressure difference at this point produces vorticity of the fluid. This phenomenon is increased as Re increases causing significant shear stresses within the boundary layers. The momentum conservation equations for fluid flow were approximated to the following system of 2D equations of fluid flow. ‡ + ‡ | This shows that the momentum fluxes are caused by the pressure gradients and viscous effects.

Numerical Solutions
The following approaches were applied in numerical solution of the flow problem.

Linearization
The dependent and independent variables in the incompressible Navier-Stokes equations governing the flow problem were non-dimensionalized then normalized to order of unity using the corresponding scaling variables. Pressure gradients were scaled with density and velocity. The variables and parameters occurring in the flow equations were scaled using the following scales: Where, ρ and μ are constants. The conservative boundary layer equation for fluid flow past a rotating sphere in 2D was then scaled as follows; The system of nonlinear equation (41) is discretized by Crank-Nicolson finite difference scheme to produce a system of linear algebraic equations [4]. The resulting system algebraic of equations is given as follows. (43) Fluid Flow by Crank-Nicolson Implicit Scheme

Neumann Boundary Conditions
The initial and boundary conditions for u x, y, t and v x, y, t were implemented numerical solutions [14]. The analytical solution the system of Navier-Stokes equations (39) at any point x, y, t is given by the following equations (53)

Numerical Solution Model
Using forward finite difference to approximate equations (44-53) gives the following respectively. Neumann's boundary conditions was used to approximate U W±q,¡ In equations (55)-(57) ω = n or n + (64)

Numerical Computation and Results
The C-N scheme developed for the N-S equations of fluid flow model formed a series of linear algebraic equations that could be expressed in matrix form. The numerical results can therefore be generated using MATLAB software.
The data used to obtain the results were k 0.001, h 0.1, l 0.1 and Re 4000 using bounded (square) computational domainΩ Ñ x, y : 0 Ò x Ò 1, 0 Ò y Ò 1. Numerical computations were performed using uniform grid, with a mesh width∆x ∆y 0.1. The results are presented in 32 ) 32 grid points graphically in three-dimensional figures.

Discussion of Results
The 3-D figures clearly show that the solutions are not changing suddenly for change in x and y hence the results for CN Scheme developed are consistent. Figures in two-dimensional space presented showed the values of the absolute errors at the corresponding points along the x-axis and y-axis are the same.
Graphs of u x, y, t and v x, y, t against Ø is presented and is exactly the same as that of u x, y, t and v x, y, t against Ù. This shows that the C-N scheme applied on the N-S equation of fluid flow problem produce stable solutions

Conclusion
The Crank-Nicholson scheme was developed and the numerical solutions for velocity fields, u x, y, t and v x, y, t of the system of Navier stokes equation for incompressible fluid flow were presented. The graphs for the velocity fields were plotted against Ø and Ù respectively using MATLAB for the C-N scheme developed. The schemes developed formed a system of linear equations that could be expressed in matrix form and therefore MATLAB software was required to generate the numerical results.
The 3-D graphical figures clearly show that the solutions are not changing suddenly for change in x and y hence the results for CN Scheme developed are consistent.
The scheme was also consistent because the results were not changing suddenly for small change in space hence convergent.

Suggestions for Further Research
This study has not explored all aspects of fluid dynamics. The compressibility analysis and Rossby number (ratio of nonlinear acceleration to) regarding rotational flows. Also the study pressure effects on the flows were by narrowing pressure gradients due to the viscous stresses.
Further explorations can be done on effects of other external forces (gravity, Lorentz, centripetal forces), magnetohydrodynamics and geographical features of atmospherics fluid flows.