Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) – solution in MATLAB®

This work presents a method for the solution of fundamental governing equations of computational fluid dynamics (CFD) using the Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) in MATLAB®. The fundamental governing equations of fluid mechanics are based on three laws of conservation, referred to as the law of conservation of mass, the law of conservation of momentum and the law of conservation of energy. The continuity equation represents the law of conservation of mass, the NavierStokes equations represent the law of conservation of momentum, and the energy equation represents the law of conservation of energy. In SIMPLE, the continuity and Navier-Stokes equations are required to be discretized and solved in a semi-implicit way. This article presents the discretization and method of solution applied to the flow around a 2-D square body. Code is written in MATLAB®. The results show the pressure and velocity fields of the converged solution.

where  represents pressure,  � �⃗ the velocity vector ( +  + ),  ⃗ =    +    +    the body force vector and  is the shear stress tensor as shown in Equation ( 5), The effect of volumetric viscosity is ignored.If the fluid is considered to be a Newtonian fluid, then shear stress is as shown in Equation (6), where  is the dynamic viscosity,   is the kronecker delta, and subscripts ,  and  refer to linear dimensions x, y and z.
The equation based on the conservation of total energy is shown in Equation ( 7 where  =    is internal energy,  2 =  2 +  2 +  2 is the magnitude of velocity vector, ̇ is the energy source,  is the thermal diffusion coefficient and  is temperature.By subtracting the product of momentum equations with their corresponding velocity from the total energy equation, an equation for the conservation of internal energy can be written as shown in Equation (8), where  = − 2 3  and  is the coefficient of dynamic viscosity.
These governing equations can be applied to a physical system in a closed volume.Equations ( 1) to ( 8) are non-linear partial differential equations and difficult to solve analytically.The numerical method of solution involves the discretization in time and space (closed volume) domains.The solution is achieved using iterative numerical techniques such as the Semi-Implicit Method for Pressure-Linked Equations (SIMPLE).
A SIMPLE algorithm is a widely used numerical procedure in computational fluid dynamics (CFD) to solve the fundamental governing equations of fluid mechanics.Prof. Brian Spalding and Suhas Patankar developed the SIMPLE algorithm at Imperial College, London, in the early 1970s [1].Since then, it has been extensively used by many researchers to solve different kinds of fluid flow and heat transfer problems [2][3][4].The literature on computational fluid dynamics discusses the SIMPLE algorithm in detail [5].A modified variant is the SIMPLER algorithm (SIMPLE Revised) introduced by Patankar in 1979 [6].
In order to be solved using numerical techniques, CFD problems need to be discretized in space dimensions [2,[7][8][9][10][11][12][13][14][15].For such numerical problems to be solved, they need to be discretized in nodes and elements [16][17][18].The nodes represent the position in space where parameters are being calculated (for example, pressure, velocities), and the elements define the equations relating to the parameters (for example, continuity and Navier-Stokes) [19].The geometry of the element defines the extent of calculations regarding space dimensions.These space dimensions categorize one-dimensional  This physical behavior of various real-time structures and physical phenomena allows engineers and mathematicians to assume a realistic 3-D situation to be a 2-D or even a 1-D problem.The advantage of reducing the dimensions is that it allows complex equations to be solved much more easily and relatively faster [7].
This paper provides a step by step methodology for SIMPLE, defining the problem, giving the assumptions, degrees of freedom, discretization, and method of solution.The methodology also discusses the flow chart of the method written in MATLAB®.The results section of the paper discusses the pressure and velocity fields obtained from the solution.The paper concludes by discussing the effectiveness and limitations of the given method.

METHODOLOGY
The methodology is divided into two sections: setting up the SIMPLE and the method of solution in MATLAB®.Setting up the SIMPLE includes: assumptions, degrees of freedom, discretization, and pressure and velocity corrections.The method of solution is given in the form of a flowchart as shown in Fig. 4.
This work assumes the CFD problem to be incompressible.This means that the density does not change in time and space and, henceforth, can be considered as a constant value.This assumption is very valid in a range of problems when flow velocities are not too high.In addition, the given CFD problem does not constitute any thermal effects.This assumption allows us to simplify the problem to the extent that there is no need to solve the energy equation.
The domain is defined as 2-D space, as shown in Fig. 2, presenting a 2-D problem of flow around a square.The variables constituting continuity and Navier-Stokes are pressure, xvelocity and y-velocity.With the above-given assumptions, continuity and Navier-Stokes equations can be simplified as shown in Equations ( 9) to (11), where  is the fluid velocity in the x-direction,  is the fluid velocity in the y-direction,  is the fluid pressure,  is the fluid density,  is the coefficient of dynamic viscosity of the fluid, and  is time.

constant velocity at inlet no-slip boundary condition
In order to apply a numerical scheme such as SIMPLE for the solution of Equations ( 9) to (11), each term needs to be discretized.Discretization is the mathematical process of transforming continuous equations into their discrete counterparts.This process is a first step toward making equations suitable for a numerical solution and for implementation on digital machines.In this work, the 2-D domain has been discretized in a staggered grid, as shown in Fig. 3.  10) and ( 11)) can be discretized using the forward difference method as shown in Equations ( 12) and ( 13 where ∆ is the value of the timestep.
Different methods can be used for the discretization of the convective term.The decision as to which method to select is based on the Reynolds number.This non-dimensional number is the ratio of convection and diffusion.If this number is small enough, we can use the central difference method; however, if the Reynolds number is higher, it is better to use the upwind scheme.Convective terms discretized using the central difference method are shown in Equations ( 14) to ( 17 where ∆ and ∆ are the values representing the size of space in the discretized domain. The pressure terms are discretized using the forward difference method, as shown in Equations ( 22) and ( 23 The diffusion terms are discretized using the central difference method, as shown in Equations ( 24) to (27), The discretized terms can be substituted in Navier-Stokes Equations ( 10) and (11), resulting in the following forms, as shown in Equations ( 28) and ( 29 The continuity equation (Equation ( 9)) is discretized using the forward difference method, as shown in Equation (30), where  indicates the extent of convergence.Although its ideal value is zero, a small number is acceptable for a numerical solution.
It is to be noted that Navier-Stokes equations (Equations ( 28) and (29)) are discretized in time domain t.However, the continuity equation (Equation (30)) has been discretized in time domain t+1.This is because the velocities in the new domain are calculated using Navier-Stokes discretized equations (Equations ( 28) and ( 29)) and checked in the discretized continuity equation (Equation (30)).The solution is iterated until the residual reduces to a small number.
Substituting the corrected values in the discretized Navier-Stokes equations (Equations ( 28) and ( 29)) results in Equations ( 34) and ( 35 where  and  can be rearranged in the form of a sparse matrix. The above equation allows solving for residuals in an implicit manner.The obtained pressure corrections can be substituted for velocity corrections (Equations ( 41) and ( 43)).To stabilize the solution, the _ is introduced with a value between 0 to 1.

RESULTS AND DISCUSSION
The pressure and velocity fields are obtained in a 2-D domain, as shown in Fig. 5.The velocity vectors show that the flow was diverted at the corners of the square.Pressure fields demonstrate that the high pressures were developed in front of the square and low pressure behind the square.In the given case, the maximum value of the residual value dropped below 10 -4 .The maximum values of residuals from the continuity equation are plotted against the iteration number, as shown in Fig. 6.

CONCLUSION
The presented study demonstrates that flow around a square in a 2-D domain can be successfully solved using the SIMPLE algorithm.The paper offers an explanation of the SIMPLE algorithm, detailing assumptions, degrees of freedom, discretization, as well as pressure and velocity corrections.The method of solution is given in the form of a flowchart and the code is implemented in MATLAB® and can be accessed from [23].
Computational fluid dynamics (CFD) is based on three basic physical principles: conservation of mass, of momentum and of energy.The governing equations in CFD are based on these conservation principles.The continuity equation is based on conservation of mass, as shown in Equation (1),  is fluid density,  is fluid velocity in the x-direction,  is fluid velocity in the ydirection,  is fluid velocity in the z-direction and  is time.The Navier-Stokes equations based on conservation of momentum are shown in Equation (2) for the x-direction, for the y-direction, ()  + ∇ • � � �⃗ � = − ___________________________________ *Corresponding Author: hassan.a.khawaja@uit.noand in Equation (4) for the z-direction, ()  + ∇ • � � �⃗ � = − (1-D), two-dimensional (2-D) or threedimensional (3-D) elements, as shown in Fig. 1 [20].

Figure 2 .
Figure 2. Domain of a 2-D problem of flow around a square.

Figure 4 .
Figure 4. Flowchart of the code.

Figure 5 .
Figure 5. Pressure and velocity fields around a square in the 2-D domain.

Figure 6 .
Figure 6.Convergence plot: the maximum value of the residual is plotted against the iteration number.