Tutorial 2 – Numerics SIMPLE Scheme

Posted on July 13, 2017 - Tutorials

This tutorial is intended for the full version of the toolbox.

Problem definition

In this tutorial we show how to implement the SIMPLE algorithm for pressure-velocity coupling problem, which occurs in the solution of incompressible flows. The original SIMPLE algorithm (Semi Implicit Method for Pressure Linked Equations) proposed in [1] has been mainly intended for finite volume discretization technique. However, we will show that the SIMPLE algorithm is basically an algorithm to solve a system of linear equations written in matrix form and resulting from the discretization of a Stokes (or Navier-Stokes) problem. The sole discretization method (finite volume or finite element) is not meaningful whatsoever for the application of the original SIMPLE procedure.

The algorithm assumes that the discretized equations can be written in matrix form as follows:

    $$A\vec{u}+Bp=\vec{f}$$

    $$C\vec{u}=0$$

In the above u stands for all velocity unknowns, $p$ for pressure, the first equation discretizes momentum equations and the second one the divergence constraint. C is typically equal to $B^{T}$ . This is the usual matrix form of equations resulting from the either finite volume or finite element discretization of the Stokes or Navier-Stokes problem. The resulting matrix

    $$\left(\begin{array}{cc} A & B\\ C & 0 \end{array}\right)$$

is regarded in literature as the Uzawa problem matrix or saddle point problem matrix

SIMPLE algorithm steps

The SIMPLE algorithm works as follows:

Set an initial pressure guess $p^{*}$ and velocity field guess $v^{*}$ and solve the momentum equations:

    $$A\vec{u^{*}}=\vec{f}-Bp^{*}$$

Let us now assume that the correct velocity field will be given by

    $$v=v^{*}+v'$$

and the correct pressure field by

    $$p=p^{*}+p'$$

Thus, equations for the correct velocity and pressure fields read:

    $$A(\vec{u^{*}}+\vec{u'})+B(p^{*}+p')=\vec{f}$$

    $$C(\vec{u^{*}}+\vec{u'})=0$$

Substratction of equations for the guessed fields from the equations for the correct fields, yields equations for pressure and velocity corrections:

    $$A\vec{u'}+Bp'=\vec{f}$$

    $$C\vec{u'}=0$$

Let us now make the main approximation of the SIMPLE algorithm. At this stage it is assumed that the equation for the velocity correction can be approximating by neglecting off-diagonal terms in the momentum equations and assume that the velocity correction equation can be written as

    $$D\vec{u'}+Bp'=\vec{f}$$

where

    $$D=diag(A)$$

Hence, velocity correction can be casted as

    $$\vec{u'}=D^{-1}(\vec{f}-Bp')$$

Let us now request that the corrected velocity field satisfies divergence constraint. To do so, we substitute velocity correction $\vec{u'}$ into following equation

    $$C(\vec{u^{*}}+\vec{u'})=0$$

which yields

    $$CD^{-1}Bp'=CD^{-1}\vec{f}+C\vec{u^{*}}$$

The above equation is the final equation for the pressure correction. Matrix $D$ is diagonal, so easily invertible, matrix multiplications are also easy to perform and the resulting matrix $CD^{-1}B$ is symmetric and positive definite, thus algorithms as conjugate gradient method can be successfully applied to solve this problem. We solve now the eqaution for pressure correction and in the next iteration we solve momentum equation for $\vec{u^{*}}$ with the new pressure field $p^{*}$ . The following procedure converges to the correct solution provided that we apply underrelaxation of pressure and velocity changes. For a Stokes problem which we consider here, it is enough to introduce underrelaxation for pressure, i.e.

    $$p=p^{*}+\alpha_{p}p'$$

For a nonlinear problem like full Navier-Stokes problem, also momentum equations must be underrelaxed.

Note here that it is wrong to underrelax the velocity correction as $\vec{v}=\vec{v^{*}}+\alpha_{v}\vec{v'}$ , but it is necessary to underrelax momentum equation. For details of implementing the SIMPLE procedure for the whole Navier-Stokes problem, refer to [2], [3] or [4].

Code

Let us now write the whole function which implements all the above steps. Below, $u$ will stand for the whole solution vector containing both velocity and pressure unknowns, $SM$ for the whole system matrix, $x0$ will contain initial velocity and pressure guess and $nVnodes$ denotes total number of velocity nodes in the compuational grid.

function u = simpleAlgorithm(SM, F, x0, nVnodes)

Let us first split the whole SM matrix into A, B and C matrices

A = SM(1:2*nVnodes,1:2*nVnodes);
B = SM(1:2*nVnodes,2*nVnodes+1:end);
C = SM(2*nVnodes+1:end,1:2*nVnodes);

Now, retrieve RHS vector of the momentum equations and initial guess of the pressure field

Fu = F(1:2*nVnodes);
p_star = x0(2*nVnodes+1:end);

Prepare also other matrices that will be needed in the algorithm and do not change from one iteration to another in case of a Stokes problem

D = diag(A);
invD = diag(1./D);
CinvDB = C*invD*B;

Define underrelaxation parameter

alpha_p = 0.02;

Start iteration loop

for iter = 1:100
    u_star = A\(Fu-B*p_star);
    p_prime = CinvDB\(C*u_star);
    p = p_star + alpha_p*p_prime;
    p_star = p;
end

Output final solution

u = [u_star; p];

Close function

end

Solving channel flow

In this section we test our algorithm by solving a simple channel flow with parabolic inlet velocity profile and two no-slip walls. We do not explain the steps to be made to solve a fluid flow problem in QuickerSim CFD Toolbox for MATLAB®. The first tutorial included in the documentation teaches all these basics.

clc;
clear;

% Fluid properties
nu = 0.1;

% Load mesh
load('channelMesh.mat');

% Convert mesh to second order
[p,e,t] = convertMeshToSecondOrder(p,e,t);

% Define inlet velocity profile as a vector (not a function handle)
velProf = [(6*p(2,:).*(1-p(2,:)))', zeros(size(p,2),1)];

% Generate indices and init solution
indices = generateIndices2D(p,t);
u = initSolution(p,t,[0 0],0);

% Assemble Stokes matrix for 2D
[SM, F] = assembleStokesMatrix2D(p, e, t, nu);%, [1 0]);

% Impose boundary conditions
[SM, F] = imposeCfdBoundaryCondition2D(p,e,t,SM,F,1,'wall');
[SM, F] = imposeCfdBoundaryCondition2D(p,e,t,SM,F,3,'wall');
[SM, F] = imposeCfdBoundaryCondition2D(p,e,t,SM,F,4,'inlet',velProf);

% Solve flow field by using SIMPLE algorithm (size(p,2) returns number of nodes)
u = simpleAlgorithm(SM, F, zeros(size(SM,1),1), size(p,2));

% Refine pressure data to the second order grid
pressure = generatePressureData(u,p,t);

% Display horizontal velocity field
figure(1)
displaySolution2D(p,t,u(indices.indu),'x-velocity')


% Display pressure field
figure(2)
displaySolution2D(p,t,pressure,'pressure');

References

[1] Patankar, S. V. and Spalding, D.B. (1972), “A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows”, Int. J. of Heat and Mass Transfer, Volume 15, Issue 10, October 1972, Pages 1787-1806.

[2] The SIMPLE Algorithm, Lecture 32 by Prof. Jayathi Y. Murthy from Purdue University, http://engineering.purdue.edu/ME608/webpage/course-organization.html (access on 28.11.2015)

[3] The SIMPLE Algorithm (Cont’d), Lecture 33 by Prof. Jayathi Y. Murthy from Purdue University, http://engineering.purdue.edu/ME608/webpage/l33.pdf (access on 28.11.2015)

[4] H. Versteeg, W. Malalasekera, “An Introduction to Computational Fluid Dynamics: The Finite Volume Method (2nd Edition).

 

Tags