INTRODUCTION TO FINITE VOLUME METHOD FOR COMPUTATIONAL FLUID DYNAMICS- II - SOLVING  ONE DIMENSIONAL FLOW FIELD

This notebook has been written in Mathematica by 
Alexandre M. S. Costa
Paulo H. Mareze
Aluizio B. Drugovich
Maringa State University
Maringa PR
Brazil

desouzadacosta@yahoo.com
http://www.geocities.com/desouzadacosta


Version:  23/8/05
More recent versions of this notebook should be available at the web site:
http://www.geocities.com/desouzadacosta/wphd/stabilityfluidiz.zip

Making use of symbolic and numeric capabilities of Mathematica, in this notebook we explore the fundamentals of the finite volume method (FVM). This method is largely employed for solution of computational fluid dynamics (CFD) problems in engineering. Nowadays, There are many commercial CFD packages available. However, for efficient use of these packages, the knowledge of the physics and mathematics behind their solutions are important In this notebook, we discuss the calculation of the flow field procedure and the following concepts related to the FVM: the pressure-velocity coupling, the staggered grid, the SIMPLE method. Finally, we apply all concepts for solving one dimensional flow field.

The application of mathematical modeling and computational simulation to processes and equipment has been almost fundamental in the development of these. Among the computational simulation tools, we detach the CFD (Computational Fluid Dynamics). This tool allows the quantitative prediction of the behavior of equipment, structures and natural phenomena. This prediction are obtained by the following steps:
·	Expressing the relevant physical and chemical laws, besides the required complementary models, in the form of equations which relates pressure, velocity, temperature and concentration values, among others, to points distributed in space and time.
·	Placing the points (i.e., the computational mesh) very close to each other in order to satisfy the continuum hypothesis.
·	Solving the equations systematically and iteratively using error reducing methods
·	Finishing the calculations when convergence is obtained, i.e., when the errors are sufficiently small.

Finally, CFD programs can manipulate a large variety of physical processes in circumstances such as:
·	Fluid flow (e.g. steam, combustion gases, oil),
·	Heat transfer,
·	Phase change,
·	Other simultaneous effects (e.g. thermal stresses, chemical reactions)

The finite volume method is a CFD method for solution of the set of differential equations that describes mass, momentum and energy conservation. The  FVM is related to other approximated numerical methods such as finite differences and finite elements (see Chung).
The finite volume method is characterized by the partition of the spatial domain in a finite number of elementary volumes for which are applied the balances of mass, momentum and energy. The approximated equations for the method can be obtained by two approaches. The first consists in applying balances for the elementary volumes (finite volumes), and the second consists in the integration spatial-temporal of the conservation equations. In this work, we are going to follow the second approach.
The general form of the differential conservation equation written in a cartesian reference frame to be treated by the FVM is:

(∂/(∂t))(ρϕ)+(∂/(∂x))(ρuϕ)+(∂/(∂y))(ρvϕ)+(∂/(∂z))(ρwϕ)=(∂/(∂x))((Γ^ϕ)((∂ϕ)/(∂x)))+(∂/(∂y))((Γ^ϕ)((∂ϕ)/(∂y)))+(∂/(∂z))((Γ^ϕ)((∂ϕ)/(∂z)))+(S^ϕ)

For two dimensional flow in the xy plane , incompressible fluid (density ρ constant), newtonian and constant properties, the variables for Eq. (1) are defined in the following table:

Values for ϕ, (Γ^ϕ), (S^ϕ)
                                             ({{     Balance
      Equation  , ϕ, (Γ^ϕ), (S^ϕ)}, {Mass, 1, 0, 0}, {x Momentum, u, μ, -((∂p)/(∂x))}, {y Momentum, v, μ, -((∂p)/(∂y))}, {Thermal Energy, T, (k/(Subscript[c, p])), (μ/(Subscript[c, p]))Φ}})

u, v, p and T are velocity components for x and y directions, pressure, and temperature respectively. μ , k and cp stands for viscosity, thermal conductivity and specific heat at constant pressure, respectively. We neglect the gravity force and thermal energy generation by viscosity Φ. Integrating the Eq. (1) over a cubic volume of dimensions  Δx, Δy, Δz (see Fig. (1)) and over time (interval between t and t + Δt):

(∫_(V,t))(∂/(∂t))(ρϕ)dVdt+(∫_(V,t))(∂/(∂x))(ρuϕ)dVdt(∫_(V,t))(∂/(∂y))(ρvϕ)dVdt=(∫_(V,t))(∂/(∂x))((Γ^ϕ)((∂ϕ)/(∂x)))dVdt(∫_(V,t))(∂/(∂y))((Γ^ϕ)((∂ϕ)/(∂y)))dVdt(∫_(V,t))(S^ϕ)dVdt

Results

(((Subscript[M, P])(Subscript[ϕ, p])-(Subscript[M, p]^o)(Subscript[ϕ, p]^o))/(Δt))+(Subscript[(OverDot[M]), e])(Subscript[ϕ, e]^θ)-(Subscript[(OverDot[M]), w])(Subscript[ϕ, w]^θ)+(Subscript[(OverDot[M]), n])(Subscript[ϕ, n]^θ)-(Subscript[(OverDot[M]), s])(Subscript[ϕ, s]^θ)=
(Subscript[(Subscript[D, 1])[((∂ϕ)/(∂x))], e]^θ)-(Subscript[(Subscript[D, 1])[((∂ϕ)/(∂x))], w]^θ)+(Subscript[(Subscript[D, 2])[((∂ϕ)/(∂y))], n]^θ)-(Subscript[(Subscript[D, 2])[((∂ϕ)/(∂y))], s]^θ)+((L[(S^ϕ)])^θ)

The process of integration described by Eqs. (2) and (3) is named discretization. The equation  (3) and all similar equations obtained by this type of process are known as discretized equations.
From Eq. (3)

(Subscript[M, p])=ρΔxΔyΔz       ;     (Subscript[M, p])=(ρ^o)ΔxΔyΔz

(Subscript[(OverDot[M]), e])=(Subscript[[ρuΔyΔz], e])    ;    (Subscript[(OverDot[M]), w])=(Subscript[[ρuΔyΔz], w])

(Subscript[(OverDot[M]), n])=(Subscript[[ρvΔxΔz], n])    ;    (Subscript[(OverDot[M]), s])=(Subscript[[ρvΔxΔz], s])

(Subscript[D, 1e])=(Subscript[[(Γ^ϕ)ΔyΔz], e])   ;  (Subscript[D, 1w])=(Subscript[[(Γ^ϕ)ΔyΔz], w])

(Subscript[D, 2n])=(Subscript[[(Γ^ϕ)ΔxΔz], n])   ;  (Subscript[D, 2s])=(Subscript[[(Γ^ϕ)ΔxΔz], s])

From the previous equations the “o” superscript stands for the time t and no superscript stands for the time t + Δt.The superscript “θ” stands for a time between t and t+Δt.In this time the ϕ variable is weighted by θ value (between 0 and 1):

(ϕ^θ) = θ ϕ + (1 - θ)(ϕ^o)

Yet, referring to Eq. (3), the volume face values are weighted by the adjacent volume values:

(Subscript[ϕ, e])=((1/2)+(Subscript[α, e]))(Subscript[ϕ, p])+((1/2)-(Subscript[α, e]))(Subscript[ϕ, E])

(Subscript[ϕ, w])=((1/2)+(Subscript[α, w]))(Subscript[ϕ, W])+((1/2)-(Subscript[α, w]))(Subscript[ϕ, p])

(Subscript[ϕ, n])=((1/2)+(Subscript[α, n]))(Subscript[ϕ, p])+((1/2)-(Subscript[α, n]))(Subscript[ϕ, N])

(Subscript[ϕ, s])=((1/2)+(Subscript[α, s]))(Subscript[ϕ, S])+((1/2)-(Subscript[α, s]))(Subscript[ϕ, p])

Similarly, the face gradient values for ϕ are weighted by the β variable:

(Subscript[[((∂ϕ)/(∂x))], e])=(Subscript[β, e])((((Subscript[ϕ, E])-(Subscript[ϕ, p]))/(Subscript[Δx, e])))

(Subscript[[((∂ϕ)/(∂x))], w])=(Subscript[β, w])((((Subscript[ϕ, p])-(Subscript[ϕ, W]))/(Subscript[Δx, w])))

(Subscript[[((∂ϕ)/(∂x))], n])=(Subscript[β, n])((((Subscript[ϕ, N])-(Subscript[ϕ, p]))/(Subscript[Δy, n])))

(Subscript[[((∂ϕ)/(∂x))], s])=(Subscript[β, s])((((Subscript[ϕ, p])-(Subscript[ϕ, S]))/(Subscript[Δy, s])))

The case for which α = 0 and β = 1 is named central difference scheme.This scheme is commonly employed for pure conduction problems.The case for which β = 0 and α = 0,5 or α = -0,5 is named upwind scheme for positive and negative velocities, respectively.
The last term from equation  corresponds to the last integral on equation . We are going to assume that:

L[(S^ϕ)]=[(Subscript[S, P])((Subscript[ϕ, P])^o) + (Subscript[S, c])] ΔV

Where 
(Subscript[S, P]) and  (Subscript[S, c]) are coefficients defined according to ϕ
ΔV is the volume Δx Δy Δz

Using eqs (4) to (17) and the discretized mass conservation equation (underlined below), equation (3) becomes:

(((Subscript[M, P])(Subscript[ϕ, p]))/(Δt))+(Subscript[ϕ, p]^θ)[(Subscript[(OverDot[M]), e])((1/2)+(Subscript[α, e]))-(Subscript[(OverDot[M]), w])((1/2)-(Subscript[α, w]))+(Subscript[(OverDot[M]), n])((1/2)+(Subscript[α, n]))-(Subscript[(OverDot[M]), s])((1/2)-(Subscript[α, s]))-(Subscript[S, p])ΔV+(Subscript[[(((Subscript[D, 1])β)/(Δx))], e])+(Subscript[[(((Subscript[D, 1])β)/(Δx))], w])+(Subscript[[(((Subscript[D, 2])β)/(Δy))], n])+(Subscript[[(((Subscript[D, 2])β)/(Δy))], s])-(UnderBar[(((Subscript[M, P])-(Subscript[M, p]^o))/(Δt))-(Subscript[(OverDot[M]), e])+(Subscript[(OverDot[M]), w])-(Subscript[(OverDot[M]), n])+(Subscript[(OverDot[M]), s])])]=(Subscript[ϕ, E]^θ)[-((1/2)-(Subscript[α, e]))(Subscript[(OverDot[M]), e])+(Subscript[[(((Subscript[D, 1])β)/(Δx))], e])]+(Subscript[ϕ, W]^θ)[((1/2)+(Subscript[α, w]))(Subscript[(OverDot[M]), w])+(Subscript[[(((Subscript[D, 1])β)/(Δx))], w])]+(Subscript[ϕ, N]^θ)[-((1/2)-(Subscript[α, n]))(Subscript[(OverDot[M]), n])+(Subscript[[(((Subscript[D, 2])β)/(Δy))], n])]+(Subscript[ϕ, S]^θ)[((1/2)+(Subscript[α, s]))(Subscript[(OverDot[M]), s])+(Subscript[[(((Subscript[D, 2])β)/(Δy))], s])]+(((Subscript[M, p]^o)(Subscript[ϕ, p]^o))/(Δt))+(Subscript[S, c])ΔV

Collecting the ϕ coefficients:

(((Subscript[M, P])(Subscript[ϕ, p]))/(Δt))+(Subscript[A, p]^( * ))(Subscript[ϕ, p]^θ)=(Subscript[A, e])(Subscript[ϕ, E]^θ)+(Subscript[A, w])(Subscript[ϕ, W]^θ)+(Subscript[A, n])(Subscript[ϕ, N]^θ)+(Subscript[A, s])(Subscript[ϕ, S]^θ)+(((Subscript[M, p]^o)(Subscript[ϕ, p]^o))/(Δt))+(Subscript[S, c])ΔV

Where

(Subscript[A, p]^( * ))=(Subscript[A, e])+(Subscript[A, w])+(Subscript[A, n])+(Subscript[A, s])-(Subscript[S, p])ΔV-((Subscript[M, P])/(Δt))+((Subscript[M, p]^o)/(Δt))

(Subscript[A, e])=-((1/2)-(Subscript[α, e]))(Subscript[(OverDot[M]), e])+(Subscript[[(((Subscript[D, 1])β)/(Δx))], e])

(Subscript[A, w])=((1/2)+(Subscript[α, w]))(Subscript[(OverDot[M]), w])+(Subscript[[(((Subscript[D, 1])β)/(Δx))], w])

(Subscript[A, n])=-((1/2)-(Subscript[α, n]))(Subscript[(OverDot[M]), n])+(Subscript[[(((Subscript[D, 2])β)/(Δy))], n])

(Subscript[A, s])=((1/2)+(Subscript[α, s]))(Subscript[(OverDot[M]), s])+(Subscript[[(((Subscript[D, 2])β)/(Δy))], s])

The equation (20) for θ = 0 (explicit formulation) becomes:

(((Subscript[M, P])(Subscript[ϕ, p]))/(Δt))=(Subscript[ϕ, p]^o)(((Subscript[M, p]^o)/(Δt))-(Subscript[A, p]^( * )))+∑(Subscript[A, nb])(Subscript[ϕ, NB]^o)+(Subscript[S, c])ΔV

And for  θ = 1 (fully implicit):

(Subscript[A, p])(Subscript[ϕ, p])=∑(Subscript[A, nb])(Subscript[ϕ, NB])  + (((Subscript[M, p]^o)(Subscript[ϕ, p]^o))/(Δt))+(Subscript[S, c])ΔV

Where

(Subscript[A, p])=∑(Subscript[A, nb])-(Subscript[S, p])ΔV+((Subscript[M, p]^o)/(Δt))

The discretized equations for a two dimensional flow corresponding to ϕ variables defined in Table 1,fully implicit formulation Eq (27) and using the constitutive Eq. (33) becomes:

((Subscript[A, p])^u) (Subscript[u, p]) = ((Subscript[A, e])^u) (Subscript[u, E]) + ((Subscript[A, w])^u) (Subscript[u, W])  + ((Subscript[A, n])^u) (Subscript[u, N]) + ((Subscript[A, s])^u) (Subscript[u, S])  - L[(P^u)] ΔV  + ((Subscript[B, p])^(u ))

((Subscript[A, p])^v) (Subscript[v, p]) = ((Subscript[A, e])^v) (Subscript[v, E]) + ((Subscript[A, w])^v) (Subscript[v, W])  + ((Subscript[A, n])^v) (Subscript[v, N]) + ((Subscript[A, s])^v) (Subscript[v, S])  - L[(P^v)] ΔV  + ((Subscript[B, p])^(v ))

((Subscript[A, p])^T) (Subscript[T, p]) = ((Subscript[A, e])^T) (Subscript[T, E]) + ((Subscript[A, w])^T) (Subscript[T, W])  + ((Subscript[A, n])^T) (Subscript[T, N]) + ((Subscript[A, s])^T) (Subscript[T, S])   + ((Subscript[B, p])^(T ))

(((Subscript[M, P])- ((Subscript[M, P])^o))/(Δt))+ (Subscript[(OverDot[M]), e]) - (Subscript[(OverDot[M]), w]) + (Subscript[(OverDot[M]), n]) - (Subscript[(OverDot[M]), s]) = 0

P = P(ρ,T)

The equations (29) to (32) represents, respectively, balances for momentum in the x and y direction, thermal energy and mass.The supercripts u,v,T stands for the ϕ variable considered.Their coefficients are given by eqs (22) to (25) and (28).
It should be noticed that:

L[(P^u)] ΔV = (Subscript[S, c]) ΔV = ∫((∂p)/(∂x))dx dy dz

L[(P^v)] ΔV = (Subscript[S, c]) ΔV = ∫((∂p)/(∂y))dy dx dz

In the terms ((Subscript[B, p])^(u )) and ((Subscript[B, p])^(v ))above are included all the other terms not multiplied by ϕ values.
It should be pointed that (Subscript[S, p]) equal to zero for equations (29) and (30).
This work focus on the solution of incompressible flow that is equivalent for the non-existence of a constitutive equation for the pressure, such as given by Eq (33).This type of flow occurs in a great number of engineering situations.Generally, the numerical solution for the incompressible flow consists on:
Calculating a pressure field that, when inserted in the momentum equations, results in a velocity field that satisfy the mass conservation equation.
One of the largely used methods that fulfill the previous statement is the SIMPLE (Semi-Implicit Method for Pressure-Linked Equations).

Before we discuss the SIMPLE steps, we present the type of variable storage, named staggered, commonly used for the method. For the staggered storage the mass balance equation is discretized using the volume centered on the pressure point (i.e., point P in the Fig 2.(a) below). On the other hand, the discretized momentum equations for x and y directions are formulated using the volumes illustrated in  Figs 2.(b) and 2.(c) below.

(a) Non staggered grid for pressure variable;(b) Staggered grid for u component;(c) Staggered grid for v component.

Steps for the SIMPLE method
1. Guess the initial velocity and pressure fields (u*,v* and P*)
2. Calculate the coefficients for the momentum equation for u and v. Here we are using small caps for referring to staggered volumes such as those given by Figure 2.(b) and 2.(c).
These staggered coefficients are obtained in a similar way to those given by equations (22) to (25) and (28).
3. Solve the momentum equations, using P*, getting u*,v*:

((Subscript[a, e])^u)(SuperStar[(Subscript[u, e])]) =  ∑ ((Subscript[a, nb])^u) (SuperStar[(Subscript[u, nb])])  + (b^u)  + ((SuperStar[(Subscript[p, P])])- (SuperStar[(Subscript[p, E])]))((Subscript[A, e])^a)

((Subscript[a, n])^v)(SuperStar[(Subscript[v, n])]) =  ∑ ((Subscript[a, nb])^v) (SuperStar[(Subscript[v, nb])])  + (b^v)  + ((SuperStar[(Subscript[p, P])])-(SuperStar[(Subscript[p, N])]))((Subscript[A, n])^a)

It should be noticed that the previous equations are similar to Eqs (29) and (30). ((Subscript[A, e])^a)  and  ((Subscript[A, n])^a) are the average east and north volume cross sectional area.(see Fig 1)
4. Calculate the pressure correction P,. We present next how to do this.
The momentum equations, dropping the superscript for convenience, using the correct pressure field are given by:

(Subscript[a, e])(Subscript[u, e]) =  ∑ (Subscript[a, nb]) (Subscript[u, nb])  +  b + ((Subscript[p, P])- (Subscript[p, E]))((Subscript[A, e])^a)

(Subscript[a, n])(Subscript[v, n]) =  ∑ (Subscript[a, nb]) (Subscript[v, nb])  +  b  +((Subscript[p, P])- (Subscript[p, N]))((Subscript[A, n])^a)

Subtracting from the previous equations the equations for the approximated field (Eqs (36) and (37)) :

(Subscript[a, e])((Subscript[u, e])^(,)) =  (UnderBar[∑ (Subscript[a, nb]) ((Subscript[u, nb])^(,))])  + (((Subscript[p, P])^(,))- ((Subscript[p, E])^(,)))((Subscript[A, e])^a)

(Subscript[a, n])((Subscript[v, n])^(,)) =  (UnderBar[∑ (Subscript[a, nb]) ((Subscript[v, nb])^(,))])  + (((Subscript[p, P])^(,))- ((Subscript[p, N])^(,)))((Subscript[A, n])^a)

Neglecting the above underlined terms and solving for the velocity fluctuations:

((Subscript[u, e])^(,)) = (((Subscript[A, e])^a)/(Subscript[a, e]))(((Subscript[p, P])^(,))- ((Subscript[p, E])^(,)))= (Subscript[d, e]) (((Subscript[p, P])^(,))- ((Subscript[p, E])^(,)))

((Subscript[v, n])^(,)) = (((Subscript[A, n])^a)/(Subscript[a, n]))(((Subscript[p, P])^(,))- ((Subscript[p, N])^(,)))= (Subscript[d, n]) (((Subscript[p, P])^(,))- ((Subscript[p, N])^(,)))

Inserting the corrected velocities in the discretized mass balance equation (similar to Eq (32)) for a volume centered on point P results in the pressure correction equation:

(Subscript[a, P]) ((Subscript[p, p])^(,))  = (Subscript[a, E]) ((Subscript[p, E])^(,)) + (Subscript[a, W]) ((Subscript[p, W])^(,)) + (Subscript[a, N]) ((Subscript[p, N])^(,)) + (Subscript[a, S]) ((Subscript[p, S])^(,)) + b

Where

(Subscript[a, E]) = (Subscript[ρ, e]) (Subscript[d, e]) Δy Δz ; (Subscript[a, W]) = (Subscript[ρ, w]) (Subscript[d, w]) Δy Δz ;   (Subscript[a, N]) = (Subscript[ρ, n]) (Subscript[d, n]) Δx Δz ;  (Subscript[a, S]) = (Subscript[ρ, s]) (Subscript[d, s]) Δx Δz ;
(Subscript[a, P]) = (Subscript[a, E]) + (Subscript[a, W]) + (Subscript[a, N]) + (Subscript[a, S ])
b = - ((((Subscript[ρ, P])- ((Subscript[ρ, P])^o)) Δx Δy Δz )/(Δt))-((Subscript[ρ, e])(SuperStar[(Subscript[u, e])]) - (Subscript[ρ, w])(SuperStar[(Subscript[u, w])]))Δy Δz -((Subscript[ρ, n])(SuperStar[(Subscript[v, n])]) - (Subscript[ρ, s])(SuperStar[(Subscript[v, s])]))Δx Δz

5. Correct u* and  v*, obtaining the velocity field that satisfies mass conservation.

(Subscript[u, e])= (SuperStar[(Subscript[u, e])]) + (Subscript[d, e])(((Subscript[P, P])^(,))-((Subscript[P, E])^(,)))

(Subscript[v, n])= (SuperStar[(Subscript[v, n])]) + (Subscript[d, n])(((Subscript[P, P])^(,))-((Subscript[P, N])^(,)))

6. Calculate P using Eq. (48) or Eq. (49)

P =  (SuperStar[P]) + (P^(,))

P =  (SuperStar[P]) + α (P^(,))

7.  Solve conservation equations for other variables, such as temperature, species concentration, etc.

(Subscript[A, P]) (Subscript[T, P]) = (Subscript[A, e]) (Subscript[T, E])  + (Subscript[A, w]) (Subscript[T, W]) +  (Subscript[A, n]) (Subscript[T, N])  + (Subscript[A, s]) (Subscript[T, S])  + ((Subscript[B, p])^(T ))

8.  Make P* = P e returns to step (2) till convergence is acquired.

In the next section, we illustrate the application of the method for a one-dimensional flow.

The example we are going to solve corresponds to proposed problem 6.7 on the book by Patankar.

Considering the one-dimensional incompressible flow through a nozzle shown in the picture

This type of flow can be described by:

(d/dx)(ρuA)=0

(d/dx)(ρuA)u = -A(dp/dx)

where A is the cross sectional area. Equations (51) and  (52) corresponds to differential balances for mass and momentum in the x direction, respectively.
We are going to set the following values for some variables:

ρ = 1; (Subscript[A, A]) = 3; (Subscript[A, B]) = 1; (Subscript[p, 1]) = 28;  (Subscript[p, 3]) = 0;

Futhermore, we are going to guess the initial values for the following variables:

ρuA = 5; (Subscript[u, A]) = (5/3); (Subscript[u, B]) = 5; (Subscript[p, 2]) = 25;

It should be noticed that the discretized momentum equation in the x direction (vide Eq. (52) ) becomes for this simple case:

(ρuA)u(Subscript[|, e]) - (ρuA)u(Subscript[|, w]) = -A((Subscript[P, e]) - (Subscript[P, w]))

Before proceeding, some notes are useful:
The Equation  (53) corresponds to equation (29) with simplifications.
We are going to use the staggered volume for the u component so we are going to rearrange equation (53) in order for being similar to equation (36) . It should be noticed that this rearrangement corresponds to steps 2 and 3 of the SIMPLE method.
The face velocities will be obtained using the UPWIND scheme for u positive, i.e.:

(Subscript[u, e]) = ((1/2)+(1/2))(Subscript[u, p]) + ((1/2)+(1/2))(Subscript[u, E]) =(Subscript[u, p])

(Subscript[u, w]) = ((1/2)+(1/2))(Subscript[u, W]) + ((1/2)-(1/2))(Subscript[u, p]) =(Subscript[u, W])

Using the above equations, the Equation (53) rewritten for the volumes centered on A and on B (see Figure (3)) becomes:

(OverDot[m])(Subscript[U, A]) = (Subscript[A, A])((Subscript[P, 1]) - (Subscript[P, 2]))

(OverDot[m])(Subscript[U, B]) = (OverDot[m])(Subscript[U, A]) + (Subscript[A, B])((Subscript[P, 2]) - (Subscript[P, 3]))

The pressure correction equation for the volume centered on point 2 (corresponding to step 4 of SIMPLE method) is:

(Subscript[a, 2])(Subscript[P, 2]^(´)) = (Subscript[a, 3])(Subscript[P, 3]^(´)) + (Subscript[a, 1])(Subscript[P, 1]^(´)) + b

The values for (Subscript[P, 3]) and (Subscript[P, 1]) values are given, and therefore their corrections  (Subscript[P, 3]^(´)) and (Subscript[P, 1]^(´)) are zero.
For equation (58)  the coefficients (Subscript[a, 2]), (Subscript[a, 3]) and (Subscript[a, 1]) corresponds to (Subscript[a, P]), (Subscript[a, E]) and  (Subscript[a, W]) (given by equations (42) ). Therefore the final correction equation for (Subscript[P, 2]) is given by:

[(Subscript[ρ, e])((Subscript[A, B])/( (OverDot[m])))(Subscript[A, B]) + (Subscript[ρ, w])((Subscript[A, A])/( (OverDot[m])))(Subscript[A, A])](Subscript[P, 2]^(´)) =-[(Subscript[ρ, e])(Subscript[U, B]^( * ))(Subscript[A, B]) + (Subscript[ρ, w])(Subscript[U, A]^( * ))(Subscript[A, A])]

The corrected velocities, corresponding to step 5 of SIMPLE method (Eq.  (46)), are given by:

(Subscript[U, B])= (Subscript[U, B]^( * )) + (((Subscript[A, B])^a)/(Subscript[a, B]))((Subscript[P, 2]^(´))-((Subscript[P, 3]^(´))=0))

(Subscript[U, A])= (Subscript[U, A]^( * )) + (((Subscript[A, A])^a)/(Subscript[a, A]))(((Subscript[P, 1]^(´))=0)- (Subscript[P, 2]^(´)))

Next, we present the Mathematica solution for the problem.

(* known variables *)
(Subscript[A, b])=1 ;
(Subscript[A, a])=3 ;
(Subscript[p, 1])=28 ;
(Subscript[p, 3])=0 ;
ρ = 1 ;
(* guesses *)
(Subscript[p, 2])=25 ;
 (Subscript[U, a ])=5/3 ;
(Subscript[U, b]) = 5 ;  
listP2nounder = {(Subscript[p, 2])} ;
listUanounder= {(Subscript[U, a ])};
listUbnounder = {(Subscript[U, b])} ;
(* loop beginning *)
Do[
ρUA = ρ  (Subscript[U, a ])  (Subscript[A, a])  ; (* step 2 - passo 2 do método simple *)
ClearAll[Uao,Ubo ];
step3  = Solve [{ρUA  Uao == (Subscript[A, a]) ((Subscript[p, 1]) - (Subscript[p, 2])),ρUA Ubo  == ρUA  Uao +(Subscript[A, b])( (Subscript[p, 2])-(Subscript[p, 3]))},{Uao,Ubo}];
(* updating *)
 Uao = Uao  /. step3[[1,1]]; 
Ubo = Ubo  /. step3[[1,2]];
(Subscript[p, 2o]) = (Subscript[p, 2]) ; 
ClearAll[p2l];
step4 = Solve[((ρ  ((Subscript[A, b])/(ρUA))(Subscript[A, b]))+(ρ ((Subscript[A, a])/(ρUA))(Subscript[A, a])))p2l == ρ (Subscript[A, a]) Uao -ρ (Subscript[A, b]) Ubo,p2l] ;
p2l = p2l /. step4[[1]] ;
(Subscript[U, a]) = Uao - ((Subscript[A, a])/(ρUA))p2l  ;
(Subscript[U, b]) = Ubo + ((Subscript[A, b])/(ρUA))p2l  ;
 (Subscript[p, 2]) = (Subscript[p, 2o]) + p2l ;
Print["(Subscript[U, a]) =  ",N[(Subscript[U, a])]," (Subscript[U, b]) =  ",N[(Subscript[U, b])]," (Subscript[p, 2]) =  ",N[(Subscript[p, 2])]];
listP2nounder = Append[listP2nounder,(Subscript[p, 2])] ;
listUanounder= Append[listUanounder,(Subscript[U, a ])];
listUbnounder = Append[listUbnounder ,(Subscript[U, b])];
,{i,1,10}]

(* It should be noticed that the above results shows (Subscript[U, a ])and (Subscript[U, b]) oscilating between two values. This is common characteristic of the SIMPLE method. One alternative for avoiding this non expected behaviour is employing the underrelaxation. The underelaxation procedure is shown in the yellow lines in the cell below *)

(* known *)
(Subscript[A, b])=1 ;
(Subscript[A, a])=3 ;
(Subscript[p, 1])=28 ;
(Subscript[p, 3])=0 ;
ρ = 1 ;
(* guesses *)
(Subscript[p, 2])=25 ;
 (Subscript[U, a ])=5/3 ;
(Subscript[U, b]) = 5 ; 
listP2under = {(Subscript[p, 2])} ;
listUaunder = {(Subscript[U, a ])};
listUbunder = {(Subscript[U, b])} ;
α = 0.5 ; 
(* loop beginning *)
Do[
ρUA = ρ  (Subscript[U, a ])  (Subscript[A, a])  ; (* step 2 *)
ClearAll[Uao,Ubo ];
step3  = Solve [{((ρUA)/α)  Uao == (Subscript[A, a]) ((Subscript[p, 1]) - (Subscript[p, 2])) + ((1-α)/α)ρUA (Subscript[U, a ]),((ρUA)/α) Ubo  == ρUA  Uao +(Subscript[A, b])( (Subscript[p, 2])-(Subscript[p, 3]))+ ((1-α)/α)ρUA (Subscript[U, b])},{Uao,Ubo}];
(* updating *)
 Uao = Uao  /. step3[[1,1]]; 
Ubo = Ubo  /. step3[[1,2]];
(Subscript[p, 2o]) = (Subscript[p, 2]) ; 
ClearAll[p2l];
step4 = Solve[((ρ  ((Subscript[A, b])/(ρUA))(Subscript[A, b]))+(ρ ((Subscript[A, a])/(ρUA))(Subscript[A, a])))p2l == ρ (Subscript[A, a]) Uao -ρ (Subscript[A, b]) Ubo,p2l] ;
p2l = p2l /. step4[[1]] ;
(Subscript[U, a]) = Uao - ((Subscript[A, a])/(ρUA))p2l  ;
(Subscript[U, b]) = Ubo + ((Subscript[A, b])/(ρUA))p2l  ;
 (Subscript[p, 2]) = (Subscript[p, 2o]) + 0.8  p2l ;
Print["(Subscript[U, a]) =  ",N[(Subscript[U, a])]," (Subscript[U, b]) =  ",N[(Subscript[U, b])]," (Subscript[p, 2]) =  ",N[(Subscript[p, 2])]];
listP2under = Append[listP2under,(Subscript[p, 2])] ;
listUaunder = Append[listUaunder,(Subscript[U, a ])];
listUbunder = Append[listUbunder ,(Subscript[U, b])];
,{i,1,10}]

(* Evaluating the next cell, it is presented a comparison between the two cases *)

<<Graphics`MultipleListPlot`
MultipleListPlot[Table[{i,listUanounder[[i]]},{i,1,Length[listUanounder]}],
Table[{i,listUaunder[[i]]},{i,1,Length[listUaunder]}],PlotLegend->{"No Underelaxation","Underelaxation (α=0.5)"},
   PlotLabel -> FontForm["(Subscript[U, a])",
        {"Helvetica-Bold",12}],DisplayFunction->$DisplayFunction,
       Ticks -> Automatic,Axes->None,Frame->True,           
            FrameLabel->{FontForm["iteration number",{"Arial",8}],
                        FontForm["Value",{"Arial",8}] },     
             GridLines->None,
             Background->GrayLevel[0.9],
                 PlotJoined ->True,
                   PlotRange -> All,PlotStyle->{PointSize[0.03], 
                           RGBColor[1,0,.5]} ];
MultipleListPlot[Table[{i,listUbnounder[[i]]},{i,1,Length[listUbnounder]}],
Table[{i,listUbunder[[i]]},{i,1,Length[listUbunder]}],PlotLegend->{"No Underelaxation","Underelaxation (α=0.5)"},
   PlotLabel -> FontForm["(Subscript[U, b])",
        {"Helvetica-Bold",12}],DisplayFunction->$DisplayFunction,
       Ticks -> Automatic,Axes->None,Frame->True,           
            FrameLabel->{FontForm["iteration number",{"Arial",8}],
                        FontForm["Value",{"Arial",8}] },     
             GridLines->None,
             Background->GrayLevel[0.9],
                 PlotJoined ->True,
                   PlotRange -> All,PlotStyle->{PointSize[0.03], 
                           RGBColor[1,0,.5]} ];

ABELL, M. L., BRASELTON, J. P., Mathematica by example, 2nd edition, Academic Press, 1997
CHUNG, T. J., Computational Fluid Dynamics, Cambridge University Press, 2002
INCROPERA, F. P., DeWITT, D. P., Fundamentals of Heat and Mass Transfer, 5th edition, Wiley, 2002
MALISKA, C. R., Heat transfer and Computational Fluid Mechanics (in portuguese), LTC, Brazil, 2004.
PATANKAR, S. V., Numerical heat transfer and Fluid Flow, Hemisphere, 1980.