[MAE]
[Project#05:
Lid Driven 
Cavity] [542] 
[Engineering Applications of Computational Fluid Dynamics] 
[ 3rd April
2011]

By Mandeep Singh 

Person # 3721 2672 
  

Contents
1  Introduction ................................................................................................................................... 3 

1.1  Streamfunction ...................................................................................................................... 3 

1.2  Vorticity .................................................................................................................................. 3 

2  Problem definition / Problem statement ...................................................................................... 3 

3  Method of Solution : ...................................................................................................................... 5 

4  Discussion of results ....................................................................................................................... 7 

4.1  Plots for vorticity and convergence of the solution for various mesh sizes .......................... 7 

4.2  Plots for the Iteration for PSOR iteration as a function of time for various relaxation factor
  9 

4.3  Contour plots for the Stream function and Vector plots for the velocities for different 
Reynolds number ............................................................................................................................. 14 

5  Summary and Conclusion ............................................................................................................ 21 

6  Appendix ...................................................................................................................................... 22 

6.1  Matlab Codes written for solving the iterations .................................................................. 22 

6.2  Boundary Condition Calculations for vorticity ..................................................................... 27 

6.2.1  Boundary 1 (Left Hand Side) ........................................................................................ 27 

6.2.2  Boundary 2 (Right Hand Side) ...................................................................................... 28 

6.2.3  Boundary 4 (Top edge) ................................................................................................. 28 

6.2.4  Boundary 3 (Bottom edge) ........................................................................................... 29 

7  References ................................................................................................................................... 30 

2 | P a g e  

 
1 Introduction
 

This project deals with the solution of vorticity, stream function and velocity fields in 
a  laminar  incompressible  flow.  Navier  Stokes  equations  are  used  to  calculate  the 
solution of equations. First we will understand the terms vorticity , stream function. 

1.1 Streamfunction
 

The streamfunction represents the two dimensional position representation flow which can 
be utilized to calculate the stream lines or the trajectories of the steady state flow. The first 
derivative  of  the  stream  function  give  the  fluild  parcels’s  velocity  and  second  derivative 
gives  the  accelerations.    A  continuous  interconnected  stream  function  gives  the 
stereamlines  for  a  given  snapshot  in  time.  Understanding  the  location  of  streamlines  is 
critical  in  studying  the  flow  pattern  for  engineering  application.  The  streamfunction  for  a 
given domain can be solved by the streamfunction equation for one instance in time if the 
initial position and magnitude of vorticity is known. 

1.2 Vorticity
Vorticity is a physical quantity in fluid dynamics in general for a fluid parcel gives a measure 
of  its  localized  rotation.  Numerically  vorticity  is  defined  as  the  divergence  of  the  velocity 
field  ×  . If  the vorticity is known everywhere in  the flow the stream function(and hence 
the velocity components) is determined by solving a Poisson equation. 

2 Problem definition / Problem statement

 

In this project we are given a lid driven square cavity in which the fluid flows over the top 
edge  with  a  velocity  (u  =  1)  in  x‐direction.  On  rest  of  the  boundaries  u  =  0.  The  lid  driven 
cavity has the non‐dimensional length of 1 x 1. Also given that the velocity component v = 0 , 

stream function   0  on all boundaries. The scheme of the project is shown in fig. 1. 

3 | P a g e  

 
 Figure 1 (Lid Driven Cavity) 

We need to solve the given system for the vorticity, stream function and velocity fields. The 
convergence  criteria for  the stream function for all the cases is  given by  ε  = 1  x 10‐10.  The 
following cases needs to be analyzed for this project  

Table 1(Cases to be evaluated) 

S. No  Condition to Explore  Parameter given 

1  Plot ω a  function of time at X = 0.5 and Y = 0 for grid  Relaxation Factor = 1 
size 5 x 5, 9 x 9, 17 x 17  for all cases. Re =10 

2  Plot number of PSOR iteration as  a  function of time for  Varying Relaxation 
grid size 5 x 5, 9 x 9, 17 x 17, Determining Optimum  Factor = 1.0  to 1.5 for 
Relaxation Factor for the grid refinement and estimating  all cases.Re=10 
the relaxation factor for 35 x 35 

3  Contour Plot of stream function and vector plot of  Cases for Re=100 , 
velocities at t = 2 for grid size 17x17  1000 , 10000 , 50000 

4 | P a g e  

 
3 Method of Solution :
 

The  driven  cavity  problem  is  a  classical  problem  that  has  wall  boundaries  surrounding  the 
entire computational region.  

 In this problem we will utilize the Neumann and Dirichlet’s boundary conditions to start our 
iteration.  Whole  system  is  divided  into  grid  of  different  sizes  as  required.  The  problem 
assume incompressible viscous flow in the cavity is driven by the uniform translation of the 
moving  upper  lid.  We  utilized  the  vorticity‐stream  function  method  to  solve  the  driven 
cavity problem. 

First we calculated plugged the boundary conditions to our system. 

For the Left , right , top and bottom boundary v = 0 ,   0 . For top boundary u = 1 and for 

rest of the boundaries u  = 0. The value of     is evaluated  over the stationary walls using 

 
u  and  v    which yields the value of     f (x)    and       f ( y )   which can 
y x
only be satisfied for   const  c  0  

Then  we  calculate  the  vorticity  for  the  boundaries  and  for  the  inner  nodes.  For  the 
boundaries , vorticity  ω  can be defined as  

For Left and right boundaries  

 2
   
x 2

For Top and bottom boundaries  

 2
   
y 2

We utilized the Taylors series formulation ( See Appendix for detail) to approximate the 
vorticity on the boundaries.  

For the inner nodes vorticity is given by relation 

5 | P a g e  

 
  uin1, j  uin1, j uin, j 1  uin, j 1  t  in1, j  2in, j  in1, j in, j 1  2in, j  in, j 1 
in, j 1  ijn  t      
 2x 2y  Re  x 2 y 2 
 

For the Stream function the values are evaluated using Gauss‐ Seidel algorithm – Point 
successive over relaxation (PSOR).  

Gauss Seidel algorithm using PSOR method for the stream function is given as  

 ik, j 1  (1  ) ijk 


2(1   )
2

 ik1, j  ik11,j   2 ( ik, j 1  ik, j 11 )  ik, j x 2 
              (2) 

Here the function (Stream function   ) is calculated for the iteration K+1 and is computed 

will we get the convergence,    is the relaxation parameter which converges the solution at 
a faster rate hence increasing the computation efficiency. We have used Gauss‐Seidel PSOR 
method with five point grid. The grid used is shown in the figure 2 below 

Figure 2(Grid Point arrangement for Gauss‐Seidel algorithm ‐ PSOR) 

Also  note  that  the  grid  arrangement  shown  above  implies  that  the  corner  nodes  are  not 
required in the boundary conditions to calculate the values of the inner nodes. In case we 
want  to  consider  to  calculate  the  corner  node  values  ,  we  need  to  apply  the  ghost  point 
method.  

The overall solution procedure for the given system can be summarized as  

1. Specify the geometry and the properties like Reynolds number, length, width , grid size 
etc 

6 | P a g e  

 
 2. Specify the initial conditions 
3. Determine  t  (non – dimensional)  
Now  t   is calculated using the stability conditions similar to what applies in the FTCS 
method for the courant number c. 

Here  t    = min (convection time t c  , diffusion time t d ) 

 2 
Now for the convection time step t   
 Re(u  v )   
2 2

 
 
Now for the diffusion time step t d   Re 
  1 1 

 2 2  2    
  x y  
 
  

in, j 1
4. Now we solve for the vorticity        

5. Solve for stream function  in, j 1  

6. Solve for    u in, j 1 and  vin, j 1         

7. Repeat step 3 – 6 until desired time or steady state is achieved. 

We can also calculate the pressure in the post processing using the Navier‐Stokes equation 
for momentum . 

4 Discussion of results

4.1 Plots for vorticity and convergence of the solution for various mesh
sizes
 

Here we plotted the vorticity  ω  as a function of time at X = 0.5 and Y = 0  for the maximum 
non‐dimensional time of  2.0. Here we have used the Relaxation factor to be 1.0.  

7 | P a g e  

 
  
Figure 3 (Plot for comparison of the Vorticity as a function of time t for various mesh sizes at location X = 0.5 and Y = 0) 

We  can  observe  form  the  plot  that  the  value  of  vorticity    ω    has  increases  by  refining  the 
mesh size. There is a significant improvement in the vorticity when we change the mesh size 
from 5 x 5 to 9 x 9.  

The  value  of  final  vorticity  after  time  2  is  0.62  compared  to  0.68  when  the  grid  is  refined 
from 5 x 5 to 9 x 9. Hence there is a change of approximately 32 % in the value of vorticity. 
When we refine the mesh from 9 x 9 to 17 x 17, we observe a change in the value of 0.68 to 
0.7 which is a change of 3.5% in the value of vorticity. 

Since  for  these  iterations,  we  used  the  relaxation  factor  to  be  1.0,  the  grid  refinement 
increases the computational time of the solution.  

Table 2(Table Showing the Computational Time for the various mesh sizes) 

S. No  Grid Size  Computational time for Re = 10 , Relaxation Factor = 1.0 

1  5 x 5  0.055439 seconds 

8 | P a g e  

 
 2  9 x 9  1.216796 seconds 

3  17 x 17  66.039699 seconds 

Here by estimating the time for the computation we can see that the convergence time for 
the  grid  size  5  x  5  comes  out  the  least  i.e.,  0.05  s  and  for  17  x  17,  it  comes  out  to  be  the 
highest  i.e.,  1  min  6s  which  is  almost  .  Though  these  values  depend  on  the  processor  and 
ram , but gives us the brief estimation of the performance of the convergence time for the 
various  mesh  grid  sizes.Hence  we  observed  that  the  convergence  occurs  by  refining  the 
mesh size but the time for the convergence increases since by increasing the mesh size , we 
have to calculate the convergence at more nodes. Also the accuracy increases upto certain 
level of mesh size. 

The  computational  time  is  represented  by  the  area  of  the  plot  in  the  following  chart  for 
comparison 

 
Figure 4 (Time estimation for the various grid sizes) 

4.2 Plots for the Iteration for PSOR iteration as a function of time for
various relaxation factor
 

9 | P a g e  

 
 In this case we changed the value of the relaxation factor for the various grid sizes i.e., 5x5, 
9x9,  17x17.  We  are  required  to  find  the  optimum  value  of  the  relaxation  factor  for  the 
various grids. 

For 5x5 grid size , we obtain the following plot for PSOR iterations as a function of time 

 
Figure 5 ( Figure showing the PSOR iterations for Stream function as a function of time for various relaxation factors for 
Grid size 5x5) 

For 5 x 5 grid size , we can see that the value of the relaxation factor comes out to be 1.2 . 
We varied the value of relaxation factor from 1 to 1.5 . 

Now for the grid size 9 x 9, we get the following fig 6.  

10 | P a g e  

 
  

Figure 6(Figure showing the PSOR iterations for Stream function as a function of time for various relaxation factors for 
Grid size 9x9) 

Here we can see that in the beginning the relaxation factor for the grid size comes out to be 
1.5 in the beginning but at the time approaches the value of 2 the optimum relaxation factor 
comes out to be 1.4 becomes optimum in the end which gives us the minimum number of 
iterations (18) compared to 22 iteration. Since if we consider the weighted average of the 
values of iteration , we can see that the optimum value of the relaxation factor comes out to 
be 1.5 

Now for the grid size of 17 x 17 we evaluated the following plot  

11 | P a g e  

 
  

Figure 7(Figure showing the PSOR iterations for Stream function as a function of time for various relaxation factors for 
Grid size 17x17) 

Hence for the grid size of 17 x 17, we can observe the value of relaxation factor comes out 
to be 1.7 

Hence we observed that the relaxation factor for the grid sizes 5x5, 9x9, 17x17 are 1.2 , 1.5, 
1.7 respectively. We utilized extrapolation for the relaxation factor we obtained the value of 
the relaxation factor for the grid size 35 x 35 to be 1.9 

12 | P a g e  

 
  
Figure 8(Relaxation Factor Extrapolation for estimating at grid size 35x35) 

As the grid size is refined, the trend shows that the optimum value of the relaxation factor 
increases. The value of the relaxation factor for the grid size 35 x 35 comes out to be 1.9. 
since  the  grid  is  refined  ,  the  number  of  computations  per  loop  increases  which  results  in 
increase in the relaxation factor.  

Hence since as the grid was refined we observed more time is required  for the computation 
hence by extrapolating we can approximate the relaxation factor the value for 35 x 35 grid, 
which can really be helpful for any study to know the approximate value of the relaxation 
factor before hand and hence reducing the computational time and cost. 

Hence we observed that the optimum value of the relaxation factor increases as we refine 
the grid. As far as the trend for the number of PSOR iterations /time step are concerned , we 
see that for a given relaxation factor , at the time t approaches 2 , the number of iterations 
required for the convergence decreases. Also as we increase the relaxation factor upto the 
optimum  value  the  iteration  per  unit  time  decreases  and  the  iterations  again  increases 
beyond the optimum relaxation factor. This increases the computational time and cost. As 
the grid is refined we require more iterations for getting the convergence in more number 
of nodes. Hence to converge the solution we have to increase the relaxation parameter. 

13 | P a g e  

 
  

4.3 Contour plots for the Stream function and Vector plots for the
velocities for different Reynolds number
For the case 3 , we are required to plot the contour for the stream function and the vector 
plot  of  velocities  at  t  =  2.0  for  different  values  of  the  Reynolds  number  i.e.,  100,  1000, 
10,000 and 50,000. 

For higher Re values the primary vortex shifts more to the centre.  For the grid size of 17 x 
17, Re = 100 , 1000 , 10000, 50000, and Relaxation factor = 1.5, the computation time is = 14 
seconds, 1.5 min , 13 min , 1 hr respectively which depend on the computer configuration 
These  values  gives  us  approximation  and  tells  the  advantage  PSOR  which  provide  the 
relaxation parameter   

Figure 9 and 10 shows the contour and velocity plots for the Reynolds number 100. Here we 
can see the primary vortex shifted towards right.The center of the vortex comes at around 
0.65 in the X direction and 0.8 in the Y‐direction. 

14 | P a g e  

 
  
Figure 9(Figure showing the stream function contour for the Grid size of 17x17 for the Reynolds number 100) 

The velocity vector plots shows the vector component of the instantaneous velocities of the 
resultant  of  u  and  v  velocities.  The  vector  plot  shows  the  velocity  direction  and  we  can 
clearly see the  circulation area near the top of the plot. This area is more clearly visualized 
using the stream function in the contour plot. 

15 | P a g e  

 
  
Figure 10(Figure showing the vector plot for the Grid size of 17x17 for the Reynolds number 100) 

For the Reynolds number of 1000 we can see that the primary vortex shifts to the right. Also the 
center of the rotation has shifted by 0.1 up and is at 0.9 compared to 0.8 as in the previous case.  

 
Figure 11(Figure showing the stream function contour for the Grid size of 17x17 for the Reynolds number 1000) 

16 | P a g e  

 
  
Figure 12(Figure showing the vector plot for the Grid size of 17x17 for the Reynolds number 100) 

Fig.  12  shows  the  velocity  plot  for  the  Reynolds  number  1000  and  we  can  see  that  the 
velocity  components  vortex  have  shifted  towards  the  right.  At  the  top  the  velocity  vector 
component have only the u component since v = 0 at the top hence the velocity vector has 
the  dimensionless  magnitude  of  1.  Also  the  center  of  rotation  has  shifted  to  the  up 
compared to the case where Reynolds number was 100.  

17 | P a g e  

 
  
Figure 13(Figure showing the stream function contour for the Grid size of 17x17 for the Reynolds number 10,000) 

 
Figure 14(Figure showing the vector plot for the Grid size of 17x17 for the Reynolds number 10,000) 

18 | P a g e  

 
 Figure  13  and  14  shows  the  variation  in  the  streamfunction  and  the  velocity  vector  at  the 
top. The center of the primary vortex has shifted to the center again. In this case we have 
changed the Reynolds number to 10000. That means now it is a turbulence model and this is 
what the change in the solution. The time step is reduced drastically for the high values of 
the Reynolds number. May be because of the high speed due to the change in the speed of 
the fluid , our method is unable to map the velocities properly for the given grid size. Hence 
we may require to refine our mesh for the proper convergence of the solution. 

Similarly we can observe that the primary vortex formed comes to the center of the cavity 
and  it  seems  like  our  scheme  does  not  work  for  very  high  Reynolds  number  for  the  given 
size of the grid. Fig 16 shows the velocity vector plot showing the vortex at the center. And 
the values of the velocities are too less compared to 1 at the top hence the velocity vector 
arrows appear very small. 

 
Figure 15(Figure showing the stream function contour for the Grid size of 17x17 for the Reynolds number 50,000) 

19 | P a g e  

 
  
Figure 16(Figure showing the vector plot for the Grid size of 17x17 for the Reynolds number 50,000) 

So we have observed that the stream lines shiftes to the right and the center od the vorticity 
liftes up when we changes the Reynolds number for 100 to 1000.  This shift occurs since the 
Re is proportional to the velocity but our velocity is constant at the top and u becomes more 
prone than v. We have also assumed the flow to be incompressible hence there is no rate of 
change in density. And we know that Re also depends on the viscosity and the dependency 
of time step on Reynolds number is bringing this change in the solution. When we increases 
the Re to 10000 or 50000, we can see the system becomes turbulent. Hence we may require 
to  use  the  turbulence  model  turbulence models are used to solve for the mean flow
behaviour and calculate the statistics of the fluctuations. Hence our Solution seems to be
unrealistic at and beyong 10000 Re 

 
 
 
 
 

20 | P a g e  

 
5 Summary and Conclusion
 

In this project we studied lid driven cavity base model for the varying grid size , relaxation 
factor and Reynolds number. Here we have observed that by varying the size of the grid the 
computational time increases and out convergence time increases since we have to find the 
convergence at more nodes. The results accuracy is enhanced by increasing the size of the 
mesh but beyond cetain mesh size the variation in accuracy is not phenomenal and it just 
adds computational time and cost to the solution convergence. The computational cost can 
be reduced by using the optimal value of the relaxation parameter. Also we have seen that 
the  relaxation  parameter  increases  by  increasing  the  mesh  size.  This  happens  because  we 
have more nodes to solve but due to refining of the mesh size we are able to use increases 
relaxation parameter. We estimated the value of the relaxation parameter for the 35 x 35 
grid using the value of the previous mesh size which can reduce our computational time and 
cost using such analogies. We also observed that by increasing the Reynolds number from 
100  to  1000,  the  center  of  primary  vortex  shifts  towards  the  right  but  at  high  Reynolds 
number like 10000, our model no longer looks physically realistic. This may be because we 
may need turbulent model scheme to solve such high Reynold number.   

21 | P a g e  

 
6 Appendix
6.1 Matlab Codes written for solving the iterations
 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Name : Mandeep Singh %
% Project 5 (Lid Driven Cavity) %
% Assigned: 04/19/11 Due: 05/03/11 %
% MAE 542 Engineering Applications of Computational Fluid Dynamics %
% The standard setup solves a lid driven cavity problem %
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%------------------------------------------------------------------%
%%

clear all

clc

%% Given
% on the right, left and bottom of the lid driven cavity
% u = v = psi(stream function) = 0
% On the top this cavity u = 1 ; v = psi(stream function) = 0
%% -----------------------------------------------------------------------
%% INPUT
tic
Re = 50000; % Reynolds Number
eps=10^-10; % Covergence criteria
om= 1.5; % Omega the Relaxation Factor
tot = 2 ; % Total Time

% dimension
Lx = 1; % Length in X-dir.
Ly = 1; % Length in Y-dir.

% Mesh Size
nx=17; %Space steps number in x-direction
ny=17; %Space steps number in y-direction

dx= Lx/(nx-1);
dy= Ly/(ny-1);
b=dx/dy; % Beta for PSOR

% Velocity over the surface

u=zeros(nx,ny);
v=zeros(nx,ny);
psi=zeros(nx,ny);
o=zeros(nx,ny);

22 | P a g e  

 
% Left Boundary
v(:,1)=0;
u(:,1)=0;
psi(:,1)=0;

% Right Boundary
v(:,ny)=0;
u(:,ny)=0;
psi(:,ny)=0;

% Bottom Boundary
v(nx,:)=0;
u(nx,:)=0;
psi(nx,:)=0;

% Top Boundary
v(1,:)=0;
u(1,:)=1;
psi(1,:)=0;

psiKP1=zeros(nx,ny);

counter = 0;
factorial = 1;
y=0;

while counter <= 2

% Calculating the time period
y=y+1;
for i=1:nx
for j=1:ny
dtc(i,j)=2/(Re*(u(i,j)^2+v(i,j)^2));
end
end
dtd=Re/(2*((1/dx^2)+(1/dy^2)));
dt=min(dtd,(min(dtc)));
dt = min(dt);

counter=counter+dt;
p(y)=counter;

for i=1:nx
for j=1:ny

%For Calculating the omega Vorticity

if j==1
% Left BC
o(i,1)=2*(psi(i,1)-psi(i,2))/(dy^2);
o1=o;
elseif j==ny
% Right BC
o(i,ny)=2*(psi(i,ny)-psi(i,ny-1))/(dy^2);
o1=o;

23 | P a g e  

 
 elseif i==nx
% Bottom BC
o(nx,j)=2*(psi(nx,j)-psi(nx-1,j))/(dx^2);
o1=o ;
elseif i==1
% Top Boundary
o(1,j)=2*((psi(1,j)-psi(2,j))/(dx^2)-((1/dx)));
o1=o;
end
end
end

for i=2:nx-1
for j=2:ny-1
o1(i,j)=o(i,j)+dt*(((u(i,j+1)*o(i,j+1)-u(i,j-1)*o(i,j-
1))/(2*dy))+...
((v(i+1,j)*o(i+1,j)-v(i-1,j)*o(i-1,j))/(2*dx)))+...
((dt/Re)*(((o(i,j+1)-2*o(i,j)+o(i,j-
1))/(dy^2))+...
((o(i+1,j)-2*o(i,j)+o(i-1,j))/(dx^2))));
end
end
o=o1;

for i= 1:y
JP(y)=o1(nx,(ny+1)/2);
end

for kl=1:100000
psiold=psi;

for i=2:nx-1
for j=2:ny-1

psi(i,j)=((1-
om)*psi(i,j))+(om*(1/(2*(1+b^2))))*(psi(i,j+1)+psi(i,j-1)+...
(b^2)*(psi(i+1,j)+psi(i-1,j))+o1(i,j)*dy^2);
% IN THE ABOVE LINE o(i,j) is the source term
end
end

m=0;
for i=2:nx-1
for j=2:ny-1
error(i,j)=(psi(i,j)-psiold(i,j))^2;
m=m+error(i,j);
end
end
jp3=sqrt(m);

jp=jp3;
if(jp3 < 10^-10)
break;

24 | P a g e  

 
 end
end
JP2(y)=kl;

for i=2:nx-1
for j=2:ny-1

u(i,j)=-(psi(i+1,j)-psi(i-1,j))/(2*dy);
v(i,j)=(psi(i,j+1)-psi(i,j-1))/(2*dx);
end
end

end

toc

psifinter = flipud(psi);
psifinal = fliplr(psifinter);

U = fliplr(u);
you = flipud(U);
V = fliplr(v);
vee = flipud(V);
x=linspace(0,1,nx);
y=linspace(0,1,ny);
quiver(x,y,you,vee);figure(gcf)
ylabel({'y-axis'});
xlabel({'x-axis'});
title({'Vector Plot for the Stream Function for Re = 10000, Grid Size 17 x
17, \Omega 1.5'});

figure (2)
contourf(x,y,psifinal,'DisplayName','psi');figure(gcf)
ylabel({'y-axis'});
xlabel({'x-axis'});
title({'Contour Plot for the Stream Function for Re = 10000, Grid Size 17
x 17, \Omega 1.5'});

% figure (1)
% plot(p,JP,'-ro')
% hold on
% xlabel('time','fontWeight','bold','fontSize',10);
% ylabel('\omega - vorticity','fontWeight','bold','fontSize',10);
% title('\omega as a function of time for Re = 10',...
% 'fontWeight',...
% 'bold','fontsize',10);
% grid on

25 | P a g e  

 
% figure (2)
% plot(p,JP2,'-b')
% hold on
% xlabel('time - t','fontWeight','bold','fontSize',10);
% ylabel('Iterations','fontWeight','bold','fontSize',10);
% title('\Omega - Number of iterations as a function of time for Re = 10 ,
Grid Size 17 x 17',...
% 'fontWeight',...
% 'bold','fontsize',10);
% grid on

26 | P a g e  

 
  
 

6.2 Boundary Condition Calculations for vorticity

We utilize the Taylor series expansion for calculating the value of vorticity at the boundaries 
of the cavity 

6.2.1 Boundary 1 (Left Hand Side)

 

For the left hand side wall of the cavity, 

 
Figure 17 

 
For approximating the stream function on the second node in the x‐direction, we used taylor series 
expansion which yield the expression 
 
  2 x 2
 2, j   1, j  x  2   (x 3 )
 
x i , j x i , j 2
 
But from definition and boundary condition we know that, 

v0
For  the left boundary x i, j
,  

 2
 i , j
Also  x 2        
i, j

Hence using these values in the taylor series expansion, we get 
 
 2( 1, j   2, j ) 
i , j    
 x 2 

27 | P a g e  

 
  
 

6.2.2 Boundary 2 (Right Hand Side)

Similar  to  the  left  boundary  condition  ,  we  calculate  the  right  side  baoundary  condition 
using the Taylor series expansion 

 
Figure 18 

 
 
 
 
The stream function for the second last node on the right hand side can be written using the Taylor 
series expansion 
 
  2 x 2
 NX 1, j   NX 1, j  x  2   (x 3 )
x NX , j x NX , j 2
 
 
Now from definition and boundary conditions
 
 
  2
v0  i , j
For  the Right hand side boundary  x i, j
, Also  x 2        
i, j

By using the above relation in the Taylor series expansion oof the stream function , we get the 
following relation for vorticity 
 
 2( NX , j   NX 1, j ) 
i , j   
 x 2   
 

6.2.3 Boundary 4 (Top edge)

 

28 | P a g e  

 
  
Figure 19 

Using the Taylor series expansion as done before, we get 
 
  2 y 2
 i , NY 1   i , NY  y  2   (y 3 )
 
y i , NY y i , NY
2
 
  2
 u 1  i , NY
For  the Top boundary  y , Also  y 2        
i, j i, j

 
Hence we obtain the relation given below for the vorticity to be 
 
 2( i , NY   i , NY 1 ) u 
i , j    
 y 2 y   
 

6.2.4 Boundary 3 (Bottom edge)

 
 

 
Figure 20 

Similarly for the bottom boundary conditions and Taylor series expansion we can obtain the vorticity 
to be 

 2( i ,1   i , 2 ) 
i , j    
 y 2 

29 | P a g e  

 
 

7 References
1.  Lecture notes by Prof. Desjardin 
2. Hoffman & Chian, Computational Fluid Dynamics Vol‐I 

30 | P a g e