Vorticity Stream Function
Incompressible N-S Equation in 2D
\[\begin{split}\begin{array}{l}
\cfrac{\partial u}{\partial t} +u\cfrac{\partial u}{\partial x}+v\cfrac{\partial u}{\partial y}=
-\cfrac{1}{\rho}\cfrac{\partial p}{\partial x}+\nu\bigg(\cfrac{\partial ^{2}u}{\partial x^{2}}+\cfrac{\partial ^{2}u}{\partial y^{2}}\bigg)\\
\cfrac{\partial v}{\partial t} +u\cfrac{\partial v}{\partial x}+v\cfrac{\partial v}{\partial y}=
-\cfrac{1}{\rho}\cfrac{\partial p}{\partial y}+\nu\bigg(\cfrac{\partial ^{2}v}{\partial x^{2}}+\cfrac{\partial ^{2}v}{\partial y^{2}}\bigg)\\
\end{array}\end{split}\]
Nondimensionalization
\[\begin{split}\begin{array}{l}
\widetilde{u}=\cfrac{u}{U},\widetilde{x}=\cfrac{x}{L},t=t\cfrac{U}{L},
\widetilde{p}=\cfrac{p}{\rho U^{2}}\\
\end{array}\end{split}\]
\[\begin{split}\begin{array}{l}
\cfrac{\partial \widetilde{u}}{\partial \widetilde{t}} +\widetilde{u}\cfrac{\partial \widetilde{u}}{\partial \widetilde{x}}+\widetilde{v}\cfrac{\partial \widetilde{u}}{\partial \widetilde{y}}=
-\cfrac{1}{\rho}\cfrac{\partial p}{\partial \widetilde{x}}+\cfrac{1}{Re}\bigg(\cfrac{\partial ^{2}\widetilde{u}}{\partial \widetilde{x}^{2}}+\cfrac{\partial ^{2}\widetilde{u}}{\partial \widetilde{y}^{2}}\bigg)\\
\cfrac{\partial \widetilde{v}}{\partial \widetilde{t}} +\widetilde{u}\cfrac{\partial \widetilde{v}}{\partial \widetilde{x}}+\widetilde{v}\cfrac{\partial \widetilde{v}}{\partial \widetilde{y}}=
-\cfrac{1}{\rho}\cfrac{\partial p}{\partial \widetilde{y}}+\cfrac{1}{Re}\bigg(\cfrac{\partial ^{2}\widetilde{v}}{\partial \widetilde{x}^{2}}+\cfrac{\partial ^{2}\widetilde{v}}{\partial \widetilde{y}^{2}}\bigg)\\
\end{array}\end{split}\]
The Vorticity Equation
\[\begin{split}\begin{align}
-\cfrac{\partial}{\partial y} &\bigg[\cfrac{\partial {u}}{\partial {t}} +{u}\cfrac{\partial {u}}{\partial {x}}+{v}\cfrac{\partial {u}}{\partial {y}} = -\cfrac{1}{\rho}\cfrac{\partial p}{\partial {x}}+\cfrac{1}{Re}\bigg(\cfrac{\partial ^{2}{u}}{\partial {x}^{2}}+\cfrac{\partial ^{2}{u}}{\partial {y}^{2}}\bigg)\bigg]\\
\cfrac{\partial}{\partial x} &\bigg[\cfrac{\partial {v}}{\partial {t}} +{u}\cfrac{\partial {v}}{\partial {x}}+{v}\cfrac{\partial {v}}{\partial {y}} = -\cfrac{1}{\rho}\cfrac{\partial p}{\partial {y}}+\cfrac{1}{Re}\bigg(\cfrac{\partial ^{2}{v}}{\partial {x}^{2}}+\cfrac{\partial ^{2}{v}}{\partial {y}^{2}}\bigg)\bigg]\\
\end{align}\end{split}\]
\[\begin{split}\begin{array}{l}
\omega=\cfrac{\partial v}{\partial x} -\cfrac{\partial u}{\partial y}\\
-\cfrac{\partial}{\partial y} \bigg(\cfrac{\partial ^{2}{u}}{\partial {x}^{2}}+\cfrac{\partial ^{2}{u}}{\partial {y}^{2}}\bigg)
+\cfrac{\partial}{\partial x}\bigg(\cfrac{\partial ^{2}{v}}{\partial {x}^{2}}+\cfrac{\partial ^{2}{v}}{\partial {y}^{2}}\bigg)
=\bigg(\cfrac{\partial ^{2}{\omega}}{\partial {x}^{2}}+\cfrac{\partial ^{2}{\omega}}{\partial {y}^{2}}\bigg)
\end{array}\end{split}\]
\[\cfrac{\partial {\omega}}{\partial {t}} +{u}\cfrac{\partial {\omega}}{\partial {x}}+{v}\cfrac{\partial {\omega}}{\partial {y}} = \cfrac{1}{Re}\bigg(\cfrac{\partial ^{2}{\omega}}{\partial {x}^{2}}+\cfrac{\partial ^{2}{\omega}}{\partial {y}^{2}}\bigg)\]
The Stream Function Equation
Define the stream function
\[\begin{array}{l}
u=\cfrac{\partial \psi}{\partial y},\quad v=-\cfrac{\partial \psi}{\partial x}
\end{array}\]
which automatically satisfies the incompressibility conditions
\[\cfrac{\partial u}{\partial x} +\cfrac{\partial v}{\partial y}=0\]
by substituting
\[\begin{split}\cfrac{\partial}{\partial x}\cfrac{\partial \psi}{\partial y} -\cfrac{\partial}{\partial y}\cfrac{\partial \psi}{\partial x}=0\\\end{split}\]
Substituting
\[u=\cfrac{\partial \psi}{\partial y},\quad v=-\cfrac{\partial \psi}{\partial x}\]
into the definition of the vorticity
\[\begin{split}\omega=\cfrac{\partial v}{\partial x} -\cfrac{\partial u}{\partial y}\\\end{split}\]
yields
\[\cfrac{\partial ^{2}{\psi}}{\partial {x}^{2}}+\cfrac{\partial ^{2}{\psi}}{\partial {y}^{2}}=-\omega\]
The Navier-Stokes equations in vorticity-stream function form are:
Advection/diffusion equation
\[\begin{split}\cfrac{\partial {\omega}}{\partial {t}} +\cfrac{\partial \psi}{\partial y}\cfrac{\partial {\omega}}{\partial {x}}-\cfrac{\partial \psi}{\partial x}\cfrac{\partial {\omega}}{\partial {y}} = \cfrac{1}{Re}\bigg(\cfrac{\partial ^{2}{\omega}}{\partial {x}^{2}}+\cfrac{\partial ^{2}{\omega}}{\partial {y}^{2}}\bigg)\\\end{split}\]
Elliptic equation
\[\cfrac{\partial ^{2}{\psi}}{\partial {x}^{2}}+\cfrac{\partial ^{2}{\psi}}{\partial {y}^{2}}=-\omega\]
Boundary Conditions for the Stream Function
At the right and the left boundary:
\[\begin{split}\begin{array}{l}
u=0\Rightarrow \cfrac{\partial \psi}{\partial y} =0\\
\Rightarrow \psi=Constant
\end{array}\end{split}\]
At the top and the bottom boundary:
\[\begin{split}\begin{array}{l}
v=0\Rightarrow \cfrac{\partial \psi}{\partial x} =0\\
\Rightarrow \psi=Constant
\end{array}\end{split}\]
Since the boundaries meet, the constant must be the same on all boundaries
\[\psi=Constant\]
The normal velocity is zero since the stream function is a constant on the wall, but the zero tangential velocity must be enforced
At the right and left boundary:
\[\begin{split}v=0\Rightarrow \cfrac{\partial \psi}{\partial x} =0\\\end{split}\]
At the bottom boundary:
\[\begin{split}u=0\Rightarrow \cfrac{\partial \psi}{\partial y} =0\\\end{split}\]
At the top boundary:
\[\begin{split}u=U_{wall}\Rightarrow \cfrac{\partial \psi}{\partial y} =U_{wall}\\\end{split}\]
The wall vorticity must be found from the streamfunction. The stream function is constant on the walls.
At the right and the left boundary:
\[\cfrac{\partial ^{2}{\psi}}{\partial {x}^{2}}+\cfrac{\partial ^{2}{\psi}}{\partial {y}^{2}}=-\omega
\Rightarrow \omega_{wall}=-\cfrac{\partial ^{2}{\psi}}{\partial {x}^{2}}\]
Similarly, at the top and the bottom boundary:
\[\cfrac{\partial ^{2}{\psi}}{\partial {x}^{2}}+\cfrac{\partial ^{2}{\psi}}{\partial {y}^{2}}=-\omega
\Rightarrow \omega_{wall}=-\cfrac{\partial ^{2}{\psi}}{\partial {y}^{2}}\]
boundary
Discretizing the Domain
Uniform mesh (h=constant)
mesh
When using FINITE DIFFERENCE approximations, the values of \(f\) are stored at discrete points and the derivatives of the function are approximated using a Taylor series:
Start by expressing the value of \(f(x+h)\) and \(f(x-h)\) in terms of \(f(x)\)
Finite Difference Approximations
\[\begin{split}\begin{array}{l}
\cfrac{\partial f(x)}{\partial x}=\cfrac{f(x+h)-f(x-h)}{2h}+\cfrac{\partial^{3} f(x)}{\partial x^{3}}\cfrac{h^{2}}{6}+\cdots \\
\cfrac{\partial ^{2}f(x)}{\partial x^{2}}=\cfrac{f(x+h)-2f(x)+f(x-h)}{h^{2}}+\cfrac{\partial^{4} f(x)}{\partial x^{4}}\cfrac{h^{2}}{12}+\cdots \\
\cfrac{\partial f(t)}{\partial t}=\cfrac{f(t+\Delta t)-f(\Delta t)}{\Delta t}+\cfrac{\partial^{2} f(t)}{\partial t^{2}}\cfrac{\Delta t}{2}+\cdots \\
\end{array}\end{split}\]
For a two-dimensional flow discretize the variables on a two-dimensional grid
Laplacian
\[\begin{split}\begin{align}
\cfrac{\partial ^{2}{f}}{\partial {x}^{2}}+\cfrac{\partial ^{2}{f}}{\partial {y}^{2}} & = \cfrac{f_{i+1,j}^{n}-2f_{i,j}^{n}+f_{i-1,j}^{n}}{h^{2}}
+\cfrac{f_{i,j+1}^{n}-2f_{i,j}^{n}+f_{i,j-1}^{n}}{h^{2}} \\ & = \cfrac{f_{i+1,j}^{n}+f_{i-1,j}^{n}+f_{i,j+1}^{n}+f_{i,j-1}^{n}-4f_{i,j}^{n}}{h^{2}}
\end{align}\end{split}\]
\[\begin{split}\cfrac{\partial {\omega}}{\partial {t}} +\cfrac{\partial \psi}{\partial y}\cfrac{\partial {\omega}}{\partial {x}}-\cfrac{\partial \psi}{\partial x}\cfrac{\partial {\omega}}{\partial {y}} = \cfrac{1}{Re}\bigg(\cfrac{\partial ^{2}{\omega}}{\partial {x}^{2}}+\cfrac{\partial ^{2}{\omega}}{\partial {y}^{2}}\bigg)\\\end{split}\]
\[\begin{split}\cfrac{\partial {\omega}}{\partial {t}} =-\cfrac{\partial \psi}{\partial y}\cfrac{\partial {\omega}}{\partial {x}}+\cfrac{\partial \psi}{\partial x}\cfrac{\partial {\omega}}{\partial {y}} + \cfrac{1}{Re}\bigg(\cfrac{\partial ^{2}{\omega}}{\partial {x}^{2}}+\cfrac{\partial ^{2}{\omega}}{\partial {y}^{2}}\bigg)\\\end{split}\]
Using these approximations, the vorticity equation becomes:
\[\begin{split}\begin{align}
\cfrac{{\omega}_{i,j}^{n+1}-{\omega}_{i,j}^{n}}{\Delta{t}} & =
-\bigg(\cfrac{\psi_{i,j+1}^{n}-\psi_{i,j-1}^{n}}{2h}\bigg)\bigg(\cfrac{\omega_{i+1,j}^{n}-\omega_{i-1,j}^{n}}{2h}\bigg)
+\bigg(\cfrac{\psi_{i+1,j}^{n}-\psi_{i-1,j}^{n}}{2h}\bigg)\bigg(\cfrac{\omega_{i,j+1}^{n}-\omega_{i,j-1}^{n}}{2h}\bigg)\\
&+\cfrac{1}{Re}\bigg(\cfrac{\omega_{i+1,j}^{n}+\omega_{i-1,j}^{n}+\omega_{i,j+1}^{n}+\omega_{i,j-1}^{n}-4\omega_{i,j}^{n}}{h^{2}} \bigg)
\end{align}\end{split}\]
The vorticity at the new time is given by:
\[\begin{split}\begin{array}{l}
{\omega}_{i,j}^{n+1}={\omega}_{i,j}^{n}+ {\Delta{t}} \bigg[
-\bigg(\cfrac{\psi_{i,j+1}^{n}-\psi_{i,j-1}^{n}}{2h}\bigg)\bigg(\cfrac{\omega_{i+1,j}^{n}-\omega_{i-1,j}^{n}}{2h}\bigg)\\
+\bigg(\cfrac{\psi_{i+1,j}^{n}-\psi_{i-1,j}^{n}}{2h}\bigg)\bigg(\cfrac{\omega_{i,j+1}^{n}-\omega_{i,j-1}^{n}}{2h}\bigg)\\
+\cfrac{1}{Re}\bigg(\cfrac{\omega_{i+1,j}^{n}+\omega_{i-1,j}^{n}+\omega_{i,j+1}^{n}+\omega_{i,j-1}^{n}-4\omega_{i,j}^{n}}{h^{2}} \bigg)\bigg]
\end{array}\end{split}\]
The stream function equation is:
\[\cfrac{\partial ^{2}\psi}{\partial x^{2}} +
\cfrac{\partial ^{2}\psi}{\partial y^{2}} =-\omega\]
\[\cfrac{\psi_{i+1,j}^{n}+\psi_{i-1,j}^{n}+\psi_{i,j+1}^{n}+\psi_{i,j-1}^{n}-4\psi_{i,j}^{n} }{h^{2}}=-\omega_{i,j}^{n}\]
Discretized Domain
Discretized Domain
Discrete Boundary Condition
Discrete Boundary Condition
\[\psi_{i,j=2}=\psi_{i,j=1}+\cfrac{\partial\psi_{i,j=1}}{\partial y} h+\cfrac{\partial ^{2}\psi_{i,j=1}}{\partial y^{2}}\cfrac{h^{2}}{2} +O(h^{3})\]
Using:
\[\omega_{wall}=-\cfrac{\partial ^{2}\psi_{i,j=1}}{\partial y^{2}} ;\quad U_{wall}=\cfrac{\partial\psi_{i,j=1}}{\partial y}\]
This becomes:
\[\psi_{i,j=2}=\psi_{i,j=1}+U_{wall}h-\omega_{wall}\cfrac{h^{2}}{2} +O(h^{3})\]
Solving for the wall vorticity:
\[\omega_{wall}=(\psi_{i,j=1}-\psi_{i,j=2})\cfrac{2}{h^{2}}+U_{wall}\cfrac{2}{h} +O(h)\]
At the bottom wall (j=1)
\[\omega_{wall}=(\psi_{i,j=1}-\psi_{i,j=2})\cfrac{2}{h^{2}}+U_{wall}\cfrac{2}{h} +O(h)\]
Similarly, at the bottom wall (j=ny):
\[\omega_{wall}=(\psi_{i,j=ny}-\psi_{i,j=ny-1})\cfrac{2}{h^{2}}-U_{wall}\cfrac{2}{h} +O(h)\]
\[\psi_{i,j=ny-1}=\psi_{i,j=ny}+\cfrac{\partial\psi_{i,j=ny}}{\partial y} (-h)+\cfrac{\partial ^{2}\psi_{i,j=ny}}{\partial y^{2}}\cfrac{h^{2}}{2} +O(h^{3})\]
Using:
\[\omega_{wall}=-\cfrac{\partial ^{2}\psi_{i,j=ny}}{\partial y^{2}} ;\quad U_{wall}=\cfrac{\partial\psi_{i,j=ny}}{\partial y}\]
This becomes:
\[\psi_{i,j=ny-1}=\psi_{i,j=ny}+U_{wall}(-h)-\omega_{wall}\cfrac{h^{2}}{2} +O(h^{3})\]
Solving for the wall vorticity:
\[\omega_{wall}=(\psi_{i,j=ny}-\psi_{i,j=ny-1})\cfrac{2}{h^{2}} -U_{wall}\cfrac{2}{h} +O(h)\]
At the left wall (i=1):
\[\psi_{i=2,j}=\psi_{i=1,j}+\cfrac{\partial\psi_{i=1,j}}{\partial x} h+\cfrac{\partial ^{2}\psi_{i=1,j}}{\partial x^{2}}\cfrac{h^{2}}{2} +O(h^{3})\]
Using:
\[\omega_{wall}=-\cfrac{\partial ^{2}\psi_{i=1,j}}{\partial x^{2}} ;\quad 0=\cfrac{\partial\psi_{i=1,j}}{\partial x}\]
This becomes:
\[\psi_{i=2,j}=\psi_{i=1,j}+0h-\omega_{wall}\cfrac{h^{2}}{2} +O(h^{3})\]
Solving for the wall vorticity:
\[\omega_{wall}=(\psi_{i=1,j}-\psi_{i=2,j})\cfrac{2}{h^{2}} +0 \cfrac{2}{h} +O(h)\]
At the right wall (i=nx):
\[\psi_{i=nx-1,j}=\psi_{i=nx,j}+\cfrac{\partial\psi_{i=nx,j}}{\partial x} (-h)+\cfrac{\partial ^{2}\psi_{i=nx,j}}{\partial x^{2}}\cfrac{h^{2}}{2} +O(h^{3})\]
Using:
\[\omega_{wall}=-\cfrac{\partial ^{2}\psi_{i=nx,j}}{\partial x^{2}} ;\quad 0=\cfrac{\partial\psi_{i=nx,j}}{\partial x}\]
Solving for the wall vorticity:
\[\omega_{wall}=(\psi_{i=nx,j}-\psi_{i=nx-1,j})\cfrac{2}{h^{2}} -0\cfrac{2}{h} +O(h)\]
Solving the elliptic equation:
\[\cfrac{\psi_{i+1,j}^{n}+\psi_{i-1,j}^{n}+\psi_{i,j+1}^{n}+\psi_{i,j-1}^{n}-4\psi_{i,j}^{n} }{h^{2}}=-\omega_{i,j}^{n}\]
Rewrite as
\[\psi_{i,j}^{n}=\cfrac{1}{4} \bigg(\psi_{i+1,j}^{n}+\psi_{i-1,j}^{n}+\psi_{i,j+1}^{n}+\psi_{i,j-1}^{n}+\omega_{i,j}^{n}{h^{2}} \bigg)\]
If the grid points are done in order, half of the points have already been updated
Successive Over Relaxation (SOR)
\[\psi_{i,j}^{n+1}=\beta\cfrac{1}{4} \bigg(\psi_{i+1,j}^{n}+\psi_{i-1,j}^{n}+\psi_{i,j+1}^{n}+\psi_{i,j-1}^{n}+\omega_{i,j}^{n}{h^{2}} \bigg)+(1-\beta)\psi_{i,j}^{n}\]
Limitations on the time step
\[\cfrac{v\Delta t}{h^{2}}\le\cfrac{1}{4} \quad\quad \cfrac{(|u|+|v|)\Delta t}{v} \le2\]