# Solving the cavity flow problem using the streamfunction-vorticity formulation

In this post I am going to write a (hopefully) simple code in matlab to solve the cavity flow problem using the vorticity stream function formulation. I will the compare the result to the result calculated by the OpenFOAM solver, icoFoam.  Writing your own solvers is fun, rewarding, and is a practice that really cements some of the fundamental knowledge of CFD.

## Contents

1. Introduction to the Cavity Flow Problem
2. Governing Equations
3. Finite Differences and Solver Algorithm
• Understanding and setting the boundary conditions
• Semi-discretization of the governing equations
• Solving the vorticity equation in time
• Solving for the vorticity at the boundaries
• Solving for the stream function
• Putting it all together
4. Results/Comparison to icoFoam
• Re=10
• Re=100
5. Conclusion and some References

## Introduction to the Cavity Flow Problem

The cavity flow problem is described in the following figure. Basically, there is a constant velocity across the top of the cavity which creates a circulating flow inside. To simulate this there is a constant velocity boundary condition applied to the lid, while the other three walls obey the no slip condition. Different Reynolds numbers give different results, so in this post I’ll simulate Re=10, and Re=800. At high Reynolds numbers we expect to see a more interesting result with secondary circulation zones forming in the corners of the cavity.

## Governing Equations

To start let’s remember that in 2D incompressible flow the Navier-Stokes equations can be rewritten to form the following two equations:

$\frac{\partial^2 \psi}{\partial x^2} +\frac{\partial^2 \psi}{\partial y^2} = -\omega$  (1)

$\frac{\partial \omega}{\partial t} + \frac{\partial \psi}{\partial y}\frac{\partial \omega}{\partial x} -\frac{\partial \psi}{\partial x} \frac{\partial \omega}{\partial y} = \nu\left(\frac{\partial^2 \omega}{\partial x^2} +\frac{\partial^2 \omega}{\partial y^2}\right)$    (2)

The two equations above are easily derived and I will probably make a short blog post covering it. For now you can take my word for it… these solve a major complication of the Navier-Stokes equations and that is the pressure-velocity coupling! By making a change of variables from u, v, and p, to $\psi$ and $\omega$ this is now in a simpler form to solve numerically.

## Finite differences and solver algorithm

### Understanding the boundary conditions

Quickly I need to explain the boundary conditions for this problem. They seem simple enough, velocity top and no velocity on the left, bottom and right sides. But remember there is no velocity in our equations! We have stream-function and vorticity. So how do we set our boundary conditions?! Well, the trick is to recognize that the outside can be considered a closed streamline. Meaning that all around the edge it will have a constant value of stream function. What should the value be? It actually doesn’t matter. But most commonly you’ll see this be set at $\psi=0$. This helps simplify some parts of the problem as you’ll see later. And what about vorticity? There is actually no boundary condition for vorticity in fact. It is an output of the simulation. At each iteration of the solver we must calculate the value of the vorticity at the wall either from the stream function or by using the same method for calculating the rest of the vorticity points. I’ll explain how we do this once we get a little further on.

### Semi-discretization of the governing equations

Now, we are going to use finite differences to discretize the equations and create a system of equations that can be solved.  Finite differences are used to approximate the derivatives at each point using the surrounding points that are a finite distance away. This is shown in the next figure. By doing this we create one equation for each node inside the cavity. When all of these equations are put together they form a system of equations that can be solved.

In the stream-function vorticity formulation there are only first and second derivatives in x and y, and a first derivative in time. For these we are going to use the finite difference equations summarized below (Eqs 3-6) which are all 2nd order.

First partial derivatives: (3)  $\frac{\partial Z}{\partial x} = \frac{Z_{i+1,j}-Z_{i-1,j}}{2\Delta x}$  (4)  $\frac{\partial Z}{\partial y} = \frac{Z_{i,j+1}-Z_{i,j-1}}{2\Delta y}$

Second partial derivatives (5) $\frac{\partial^2 Z}{\partial x^2} = \frac{Z_{i+1,j}-2 Z_{i,j}+ Z_{i-1,j}}{\Delta x^2}$ (6) $\frac{\partial^2 Z}{\partial y^2} = \frac{Z_{i+1,j}-2 Z_{i,j}+ Z_{i-1,j}}{\Delta x^2}$

The first step is to insert equations 3-6 into equations 1 and 2. We now have a set of vorticity stream function equations based on finite differences. For simplicity I will solve this problem on a perfectly square grid where $\Delta x = \Delta y = h$. This allows the equations and the coding to be much simpler. We end up with:

(7) $\psi_{i,j}=1/4\left[\psi_{i+1,j}+\psi_{i-1,j}+\psi_{i,j+1}+\psi_{i,j-1}+h^2\omega_{i,j}\right]$ (rearranged for $\psi_{i,j}$ and simplified)

(8) $\frac{\partial \omega}{\partial t} +\frac{\left(\psi_{i,j+1}-\psi_{i,j-1}\right)\left(\omega_{i+1,j}-\omega_{i-1,j}\right)}{4h^2} - \frac{\left(\psi_{i+1,j}-\psi_{i-1,j}\right)\left(\omega_{i,j+1}-\omega_{i,j-1}\right)}{4h^2}= \nu\left(\frac{\omega_{i+1,j}+\omega_{i,j+1}+\omega_{i-1,j}+\omega_{i,j-1}-4\omega_{i,j}}{h^2}\right)$

Next the process can go several ways. We can solve Equation 1 and 2 simultaneously or sequentially, and we could solve Equation 1 first or Equation 2 first. As well we could modify the equations based on some information we already know to be true. For instance, we know ahead of that this problem is a “steady state problem” which means that it doesn’t vary in time. Yet, equation 2 has  a time derivative in it. In one approach, we could get rid of this term from the equation which we know to be zero and write a steady state solver. Alternatively (and the way I chose to do it here) is we use that time derivative to solve the equation in time and as the time gets large (as $t\rightarrow\infty$) the simulation should stop changing and the time derivative would go to zero. Here we are going to solve the equations sequentially starting with the vorticity equation using the time derivative (Equation 2) followed by Equation 1.

### Solving the vorticity equation in time

So to use this  approach we are going to use the most simple time discretization – the explicit Euler scheme (1st order). The easiest way to think about this is we can rearrange Equation 2 to give an equation of the form:

(9) $\frac{\partial \omega}{\partial t}=f(x)$

where f(x) is the function formed by our finite difference equations. We then know from the first order Taylor series that we can say:

(10) $\omega^{k+1}=\omega^{k}+\Delta t\frac{\partial \omega}{\partial t}^{k}$

where: $\frac{\partial \omega}{\partial t} ^k= f_{i,j}^k=\left( -\frac{\left(\psi_{i,j+1}-\psi_{i,j-1}\right)\left(\omega_{i+1,j}-\omega_{i-1,j}\right)}{4h^2} + \frac{\left(\psi_{i+1,j}-\psi_{i-1,j}\right)\left(\omega_{i,j+1}-\omega_{i,j-1}\right)}{4h^2}\right)+ \nu\left(\frac{\omega_{i+1,j}+\omega_{i,j+1}+\omega_{i-1,j}+\omega_{i,j-1}-4\omega_{i,j}}{h^2}\right)$

Here means the time step index. Meaning k is the current time step and k+1 is the next time step. $\Delta t$ is the chosen time step size. This must be chosen based on stability, or a desired time accuracy. Using the equation above we can predict the vorticity throughout the field at the next timestep. This is really the key to the solver we are writing.

### Solving for the vorticity at the boundaries

In the last step I showed how the vorticity can be predicted in time using the explicit Euler formulation. But this is not all that needs to be done to define the vorticity. In the middle points all of the central differences are valid because all of the points have neighbors in all directions. However, at the boundaries things get a little interesting. They don’t have a full set of neighbors! For example on the left there are no nodes for which to use in the central difference equations. There are a couple different ways to do this. The most obvious way is to use one-sided finite differences to formulate the equations at these nodes. However there is an easier way that comes from a simple observation of the Taylor series.

Using the top wall as an example, from Equation 1 we know that:

(11) $\left(\frac{\partial^2 \psi}{\partial y^2}\right)_{i,N} = \omega_{i,1}$

This is because the stream function is constant along the top and therefore $\left(\frac{\partial^2 \psi}{\partial x^2}\right)_{i,1} = 0$ . From the Taylor series we know that

(12) $\psi_{i,N-1}=\psi_{i,N}-\left(\frac{\partial \psi}{\partial y} \right)_{i,N}h +\left(\frac{\partial^2 \psi}{\partial y^2} \frac{h^2}{2}\right)_{i,N} + \textnormal{...}$

There are several properties in the above expansion that we know. For instance, $\left(\frac{\partial \psi }{\partial y}\right)_{i,N}=U_{top}$, and also that $\left(\frac{\partial^2 \psi}{\partial y^2}\right)_{i,N} = \omega_{i,1}$. By subbing these in and rearranging for vorticity we get:

(13) $\omega_{i,N}=2\frac{\psi_{i,N}-\psi_{i,N-1}}{h^2}-\frac{2U_{top}}{h}$

We also know that we set $\psi_{i,N}=0$ so the actual equation for the top boundary becomes:

(14a) $\omega_{i,N}=2\frac{-\psi_{i,N-1}}{h^2}-\frac{2U_{top}}{h}$  (top)

Similarly equation for the other boundaries are:

(14b) $\omega_{i,1}=2\frac{-\psi_{i,2}}{h^2}$  (bottom)

(14c) $\omega_{1,j}=2\frac{-\psi_{2,j}}{h^2}$ (left)

(14d) $\omega_{N,j}=2\frac{-\psi_{N-1,j}}{h^2}$ (right)

### Solving for the streamfunction

After predicting the vorticity at the next time step, we can solve Equation 1 for the corresponding stream function. To solve Equation 1, the way that I have chosen to do it is to rearrange the finite difference equation for the stream-function psi at indices i,j:

(15) $\psi_{i,j}=1/4\left[\psi_{i+1,j}+\psi_{i-1,j}+\psi_{i,j+1}+\psi_{i,j-1}+\omega\right]$

Then it is a matter of using a simple loop structure (probably a for loop) to loop through all the nodes and calculate the predicted psi. Since we are running until we reach steady state, after each time step we should check the value of the time derivative $\frac{\partial \omega}{\partial x}$ and when it reaches an acceptably small value we will call the simulation converged.

### Putting it all together

So to summarize the steps for the solver are:

1. Choose an arbitrary starting point (ie assuming a streamfunction and vorticity at all points)
2. a) Loop through all interior nodes and calculate the time derivative $\partial \omega / \partial t$, and use explicit Euler to guess the vorticity at each point at the next time-step (k+1)
b) At boundary nodes, calculate the vorticity at the boundaries using the simple Taylor series approximation
3. Loop through all interior nodes and solve Equation 1 to predict the stream-function. Keep the stream-function at all boundaries the same ($\psi=0$)
4. Check for convergence (based on some user defined criteria)
5. If converged the simulation is done, you can now post-process and calculate things such as velocity and pressure. If not converged go back to step 2 and continue forward solving in time.

Below is a schematic summary of the solver algorithm I just described:

## Results

To examine the results of my solver I am going to simulate two Reynolds numbers 10, and 800 using both my own solver, and the OpenFoam solver on a comparable grid size. I am assuming the OpenFoam solver to be the superior being since it is written by folks much smarter than I. The OpenFOAM solver that I am using is the solver icoFoam. And coincidentally the OpenFOAM download comes with the cavity problem already set up and all I had to do was hit run! … (After changing the Reynolds number and grid size of course)

### Re = 10

Because of the low Reynolds number I expected both solvers to perform similarly. I compared the output velocity, streamfunction, and vorticity fields. Both of the simulations were done on a 60×60 grid. Take a look:

At this Reynolds number it seems like my solver did a pretty darn good job! The stream function, vorticity, and velocity plots match quite closely. And they look kinda pretty too! I plotted the OpenFOAM output in Matlab so that the images would be plotted the same. Comparing matlab to Paraview just didnt do it for me. Maybe I’ll do a blog about but putting openfoam into Matlab… at a later date.

### Re = 800

Re 800 is supposed to be a much more challenging case. Why? Well for one… the Reynolds number is higher. This means that shear layers are thinner compared to the problem length scale and therefore require more resolution. Additionally, we (or I anyway) know that at these Reynolds numbers reversed flow and separation bubbles are produced in the corners.

For this one each simulation (both icoFoam and the matlab solver I wrote) are solved on a 100×100 grid. The higher resolution is necessary because of the more complex flow. So… lets see how it went!

On thing we notice right away is how much prettier Re 800 is! But again it looks like the solvers match up pretty closely! The vorticity plots are on a log scale. This is a great example of the flow described by Batchelor about high Reynolds number flow with closed streamlines. We have vorticity being convected around the outside of a mostly constant vorticity core.

## Conclusion and some References

If you read this far… thanks for paying attention! This was my first post and first time writing a solver with (maybe) some people watching/reading. If anything is unclear or incorrect (I do make lots of mistakes) comment or something. Here is a link to the matlab code I wrote! Feel free to use it… but if you publish anywhere please reference it: streamvortRe10

Some useful references and alternative sources if you are very interested in this method are:

##### Website pages:

[1] https://www.iist.ac.in/sites/default/files/people/psi-omega.pdf

[2] http://nptel.ac.in/courses/112104030/pdf/lecture24.pdf

[3] http://www.cfm.brown.edu/people/gk/chap11/node1.html

[4] http://www.hindawi.com/journals/isrn/2012/871538/

##### Books

[5] J. C. Tannehill, D. A. Anderson, R. H. Pletcher”Incompressible Navier-Stokes Equations”, Computational Fluid Mechanics and Heat Transfer, 2nd Edition, Taylor and Francis, 1997, pp. 649-659

[6] K. H. Huebner, E. A. Thornton, T. G. Byrom ,”Fluid Mechanics Problems” The Finite Element Method for Engineers, 2nd Edition, John Wiley and Sons, 1995

##### OpenFOAM

[7] OpenCFD, “OpenFOAM- The Open Source Computational Fluid Dynamics (CFD) Toolbox”, http://www.openfoam.com

## 22 thoughts on “Solving the cavity flow problem using the streamfunction-vorticity formulation” Leave a comment ›

1. Ahmed says:

Hi – I am trying to compare the mid-longitudinal u-velocity profile and lateral v-velocity profile to well known benchmarks. But the result is off. Even qualitatively. And I am trying for modest Re = 50 or 100. Has anyone tried this? I will grateful for any insights

2. Pramod says:

Many many thanks for the lucid explanation. Keep writing.!! It’s wonderful

3. mathgra says:

hello

Im master student, as a final project my teacher give a steady flow in a lid-driven cavity by using Naiver-Stokes equations in stream function vorticity form with deltax=deltay=0.1 Re=50,100 stopping criteria for both streamfunctions and vorticity is 0.0001.
initial vorticity is 0 everywhere , stream function zero at boundaries and partial derivative of streamfunction with respect to y is 1. So the only top lid is moving.
We use finite difference for the scheme.
I have to drive streamlines and vorticity equations. I try to write matlab code with help of your code and program runs.First problam it works so slow and second it gives zero as answer.
But I think I missing some point and I am not sure how to put stopping criteria in loop.

Matlab code that I write
clear all
close all
n=11;
x=linspace(0,1,n);
y=linspace(0,1,n);
h=x(2)-x(1);
criteria=0.0001;
Re=50;
[xx,yy]=meshgrid(x,y);
vortnew=zeros(n,n);
psi=zeros(n,n);
U=1;
run=1;
iter=0;

for i=1:n
for j=1:n
psi(i,j,1)=0;
vortnew(i,j,1)=0;
if j==n
psi(i,j,1)=0;
vortnew(i,j,1)=-2/h;
end
end
end

% vorticity equation
if run==1;

while ( 1>criteria )

vortold=vortnew;
psinew=psi;
iter=iter+1;
for i=1:n
for j=1:n
if j~=1 && j~=n && i~=1 && i~=n % Eger boundaryde degilse

f(i,j)=-((psi(i,j+1)-psi(i,j-1))*(vortnew(i+1,j)-vortnew(i-1,j)))/4*(h^2)
+((psi(i+1,j)-psi(i-1,j))*(vortnew(i,j+1)-vortnew(i,j-1)))/4*(h^2)
+(vortnew(i+1,j)-4*vortnew(i,j)+vortnew(i-1,j)+vortnew(i,j+1))/Re*(h^2);

elseif j==1 && i~=1 && i~=n
vortnew(i,j)=2*(psi(i,j)-psi(i,j+1))/h^2; %Bottom
elseif j==n && i~=1 && i~=n
vortnew(i,j)=2*(psi(i,j)-psi(i,j-1))/h^2-2*U/h; % Top
elseif i==1 && j~=1 && j~=n
vortnew(i,j)=2*(psi(i,j)-psi(i+1,j))/h^2; % Left
elseif i==n && j~=1 && j~=n
vortnew(i,j)=2*(psi(i,j)-psi(i-1,j))/h^2; % Right
elseif (i==1 && j==1) %corners equal to neighbors or zero
vortnew(i,j)=0;
elseif (i==n && j==1)
vortnew(i,j)=0;
elseif (i==1 && j==n)
vortnew(i,j)=(vortnew(i+1,j));
elseif (i==n && j==n)
vortnew(i,j)=(vortnew(i-1,j));
end
end
end
end
for i=1:n
for j=1:m
vortold(i,j)=vortnew(i,j)+(1-0.05)*f(i,j);
% To guarantee convergence she give above statement
end
end
%streamfunction equation non-boundary points
for i=1:n
for j=1:n
if j~=1 && j~=n && i~=1 && i~=n
psi(i,j)=0.25*(h^2*vortnew(i,j)+psi(i+1,j)+psi(i-1,j)+psi(i,j+1)+psi(i,j-1));
else
psi(i,j)=0;
end
end
end
end

4. xcvb2 says:

Hi I am trying to increase the number of data point and change the lid speed however I keep getting the error message vortout{s}(:,:) or uout{s}(:,:) has an assignment dimension mismatch I was wondering if you could explain why this is and how I could change the code

5. I see you don’t monetize your website, don’t waste your
traffic, you can earn additional cash every month because
you’ve got hi quality content. If you want to know how to
make extra money, search for: Mrdalekjd methods for \$

6. Mike Strombach says:

Using this code seems to require an existing file ‘lidDrivenCavity.mat’. Can you provide this file?

• Try setting run=1. Then when you run it for the first time it saves your results to lidDrivenCavity.mat. Or you can just comment out that if statement. The idea when I wrote it was that if you set run=0 it meant you just wanted to plots the results you already had.

Let me know if it works. Haven’t looked at this file for a while.

7. AAKASH R says:

Hey,I tried your code , but I have a problem in plotting the secondary vortices. Your code also unlike the picture depicted doesn’t show the secondary vortice.

• The code posted here is for Re 10. You will need to change the code for the new Reynolds number. As well, since the value of the stream function is so low in the secondary vortices you need to change the plotting options in order for them to show.

8. ravi says:

how can i modify the boundary condition and code for back-step flow situation

• Hi Ravi, This would mean modifying the mesh (it’s not square), and the boundary conditions (in-flow, out-flow, wall etc). I think too much to explain in a comment. Perhaps I will do this as a blog post in the future. Thanks for the idea!

9. Kohos says:

Should the second derivative of z wrt y be different? (equation 6)

• Yes! Good catch. The indices for the +/-1 should be the j index. Since that is the index in the y direction.

10. Miquel says:

I’m trying to simulate the cavity-lid problem, but I believe that your equation (7) is not right, I believe “w” should be “w*(Ah^2)”. Please correct me if I’m wrong. Many thanks in advance.

• You are right! I didn’t write out that equation properly. It should be w*h^2 I believe. I have it correctly in the solver process diagram.

• Miquel says:

Ohh!! yes I just see it now, many thanks for your quick reply!

11. Miquel says:

Hi, I’m trying to do a similar study that you have done, but I have a question in your equation (7), I believe the vorticity should be -w*(Ah^2). Could you confirm that this is correct? Many thanks in advance, is a great job.

12. I am not able to capture vortices near the lower wall of cavity. Why is this so?

• It could be because there is an error in your code. Or that you don’t have enough points in the x and y direction. Also it could just be your visualization: the values for the stream function magnitudes are very small in the corner vortices. You could try playing with how you are plotting.

13. Can you simulate for Re=100 and plot streamlines in OpenFoam. My linux is not running and I don’t know how to change Re in Openfoam available for Windows. I want to validate my numerical results. Thanks for quick reply. I appreciate your work which I discovered at the right time.

14. What is the value of psi on moving wall?

• It’s 0! The outside of the cavity is considered a single streamline (psi=0). The fact that the lid is moving is accounted for in the vorticity calculation.

Thanks for reading! And Thanks for the comment! I think I’ll add a diagram to better illustrate the boundary conditions and maybe reorganize the post a bit.

curiosityFluids

This site uses Akismet to reduce spam. Learn how your comment data is processed.