ui-1/2 = 0.5*(u(i,j)+u(i-1,j))
if ui-1/2 >0 then
(uw)i-1/2 = ui-1/2 wi-1 ; b'=b + $\delta$ Re ui-1/2else
(uw)i-1/2 = ui-1/2 wi ; e'=e + $\delta$Re ui-1/2endif
ui+1/2 = 0.5*(u(i,j)+u(i+1,j))
if uI+1/2 >0 then
(uw)i+1/2 = ui+1/2 wi ; e'=e - $\delta$ Re ui+1/2else
(uw)i+1/2 = ui+1/2 wi+1 ; a'=a - $\delta$ Re ui+1/2endif
The y-direction will be similar except the aspect ratio, , must be included. These modifications to the coefficients must be made before the coefficients are updated for the boundary conditions.
Calculation of pressure
The vorticity-stream function method does not require calculation of pressure to determine the flow field. However, if the force or drag on a body or a conduit is of interest, the pressure must be computed to determine the stress field. Here we will derive the Poisson equation for pressure and determine the boundary conditions using the equations of motion. If the pressure is desired only at the boundary, it may be possible to integrate the pressure gradient determined from the equations of motion. The dimensionless equations of motion for incompressible flow of a Newtonian fluid are as follows.
$$\begin{array}{ccc}\frac{d\mathbf{v}*}{dt*}\hfill & =& \frac{\partial \mathbf{v}*}{\partial t*}+\mathbf{v}*\u2022\nabla *\mathbf{v}*=\frac{\partial \mathbf{v}*}{\partial t*}+\nabla *\u2022(\mathbf{v}*\mathbf{v}*)=\hfill \\ & & -\nabla *P*+\frac{1}{Re}\nabla {*}^{2}\mathbf{v}*\hfill \\ & & \mathrm{where}\hfill \\ \hfill \mathbf{v}*& =& \frac{\mathbf{v}}{U},\phantom{\rule{1.em}{0ex}}\mathbf{x}*=\frac{\mathbf{x}}{L},\phantom{\rule{1.em}{0ex}}t*=\frac{U}{L}t,\phantom{\rule{1.em}{0ex}}P*=\frac{P}{\rho \phantom{\rule{0.166667em}{0ex}}{U}^{2}},\phantom{\rule{1.em}{0ex}}\nabla *=L\phantom{\rule{0.166667em}{0ex}}\nabla \hfill \\ \hfill P& =& p-\rho \phantom{\rule{0.166667em}{0ex}}g\phantom{\rule{0.166667em}{0ex}}z,\phantom{\rule{1.em}{0ex}}Re=\frac{\rho \phantom{\rule{0.166667em}{0ex}}U\phantom{\rule{0.166667em}{0ex}}L}{\mu}\hfill \end{array}$$
Here all coordinates were scaled with respect to the same length scale. The aspect ratio must be included in the final equation if the coordinates are scaled with respect to different length scales. We will drop the * and use
$x$ and
$y$ for the coordinates and
$u$ and
$v$ as the components of velocity.
The following derivation follows that of Hoffmann and Chiang (1993). The equations of motion in 2-D in conservative form is as follows.
$$\begin{array}{c}\frac{\partial u}{\partial t}+\frac{\partial \left({u}^{2}\right)}{\partial x}+\frac{\partial \left(u\phantom{\rule{0.166667em}{0ex}}v\right)}{\partial y}=-\frac{\partial P}{\partial x}+\frac{1}{Re}\left(\frac{{\partial}^{2}u}{\partial {x}^{2}}+\frac{{\partial}^{2}u}{\partial {y}^{2}}\right)\hfill \\ \frac{\partial v}{\partial t}+\frac{\partial \left(u\phantom{\rule{0.166667em}{0ex}}v\right)}{\partial x}+\frac{\partial \left({v}^{2}\right)}{\partial y}=-\frac{\partial P}{\partial y}+\frac{1}{Re}\left(\frac{{\partial}^{2}v}{\partial {x}^{2}}+\frac{{\partial}^{2}v}{\partial {y}^{2}}\right)\hfill \end{array}$$
The Laplacian of pressure is determined by taking the divergence of the equations of motion. We will carry out the derivation step wise by first taking the x-derivative of the x-component of the equations of motion.
$$\frac{\partial}{\partial t}\left(\frac{\partial u}{\partial x}\right)+2\frac{\partial u}{\partial x}\frac{\partial u}{\partial x}+2u\frac{{\partial}^{2}u}{\partial {x}^{2}}+\frac{\partial u}{\partial x}\frac{\partial v}{\partial y}+u\frac{{\partial}^{2}v}{\partial x\partial y}+\frac{\partial v}{\partial x}\frac{\partial u}{\partial y}+v\frac{{\partial}^{2}u}{\partial x\partial y}=-\frac{{\partial}^{2}P}{\partial {x}^{2}}+\frac{1}{Re}\frac{\partial}{\partial x}\left({\nabla}^{2}u\right)$$
Two pair of terms cancels because of the equation of continuity for incompressible flow.
$$\begin{array}{c}\frac{\partial u}{\partial x}\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)=0\hfill \\ u\frac{{\partial}^{2}v}{\partial x\partial y}+u\frac{{\partial}^{2}u}{\partial {x}^{2}}=u\frac{\partial}{\partial x}\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)=0\hfill \\ \mathrm{Thus}\hfill \\ \frac{\partial}{\partial t}\left(\frac{\partial u}{\partial x}\right)+{\left(\frac{\partial u}{\partial x}\right)}^{2}+u\frac{{\partial}^{2}u}{\partial {x}^{2}}+\frac{\partial v}{\partial x}\frac{\partial u}{\partial y}+v\frac{{\partial}^{2}u}{\partial x\partial y}=-\frac{{\partial}^{2}P}{\partial {x}^{2}}+\frac{1}{\mathrm{Re}}\frac{\partial}{\partial x}\left({\nabla}^{2}u\right)\hfill \end{array}$$
Similarly, the y-component of the equations of motion become
$$\frac{\partial}{\partial t}\left(\frac{\partial v}{\partial y}\right)+{\left(\frac{\partial v}{\partial y}\right)}^{2}+v\frac{{\partial}^{2}v}{\partial {y}^{2}}+\frac{\partial u}{\partial y}\frac{\partial v}{\partial x}+u\frac{{\partial}^{2}v}{\partial x\partial y}=-\frac{{\partial}^{2}P}{\partial {y}^{2}}+\frac{1}{Re}\frac{\partial}{\partial y}\left({\nabla}^{2}v\right)$$
The x and y-components of the above equations are now added together and several pairs of terms cancels pair-wise.
$$\begin{array}{c}\frac{\partial}{\partial t}\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)=0\hfill \\ \frac{{\partial}^{2}u}{\partial {x}^{2}}+\frac{{\partial}^{2}v}{\partial x\partial y}=\frac{\partial}{\partial x}\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)=0\hfill \\ \frac{{\partial}^{2}u}{\partial x\partial y}+\frac{{\partial}^{2}v}{\partial {y}^{2}}=\frac{\partial}{\partial y}\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)=0\hfill \\ \frac{\partial}{\partial x}\left(\frac{{\partial}^{2}u}{\partial {x}^{2}}+\frac{{\partial}^{2}u}{\partial {y}^{2}}\right)+\frac{\partial}{\partial y}\left(\frac{{\partial}^{2}v}{\partial {x}^{2}}+\frac{{\partial}^{2}v}{\partial {y}^{2}}\right)=\frac{{\partial}^{2}}{\partial {x}^{2}}\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)+\frac{{\partial}^{2}}{\partial {y}^{2}}\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)=0\hfill \end{array}$$
Therefore the equations reduce to
$${\left(\frac{\partial u}{\partial x}\right)}^{2}+{\left(\frac{\partial v}{\partial y}\right)}^{2}+2\frac{\partial u}{\partial y}\frac{\partial v}{\partial x}=-\left(\frac{{\partial}^{2}P}{\partial {x}^{2}}+\frac{{\partial}^{2}P}{\partial {y}^{2}}\right)$$
The left-hand side can be further reduced by consideration of the continuity equation as follows.
$$\begin{array}{c}{\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)}^{2}={\left(\frac{\partial u}{\partial x}\right)}^{2}+{\left(\frac{\partial v}{\partial y}\right)}^{2}+2\frac{\partial u}{\partial x}\frac{\partial v}{\partial y}=0\hfill \\ \mathrm{from}\phantom{\rule{0.277778em}{0ex}}\mathrm{which}\hfill \\ {\left(\frac{\partial u}{\partial x}\right)}^{2}+{\left(\frac{\partial v}{\partial y}\right)}^{2}=-2\frac{\partial u}{\partial x}\frac{\partial v}{\partial y}\hfill \end{array}$$
Upon substitution, we have
$$-\left(\frac{{\partial}^{2}P}{\partial {x}^{2}}+\frac{{\partial}^{2}P}{\partial {y}^{2}}\right)=2\left(\frac{\partial u}{\partial y}\frac{\partial v}{\partial x}-\frac{\partial u}{\partial x}\frac{\partial v}{\partial y}\right)$$
This equation can be written in terms of the stream function.
$$\frac{{\partial}^{2}P}{\partial {x}^{2}}+\frac{{\partial}^{2}P}{\partial {y}^{2}}=2\left[\frac{{\partial}^{2}\psi}{\partial {x}^{2}}\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}\frac{{\partial}^{2}\psi}{\partial {y}^{2}}-{\left(\frac{{\partial}^{2}\psi}{\partial x\partial y}\right)}^{2}\right]$$