I'm particularly interested in 2D incompressible version of Navier-Stokes equations that describe flow in some domain of interest:
$$\begin{aligned}\frac{\partial v_x}{\partial t} + \vec{v} \cdot \vec{\nabla} \, v_x &= – \frac{1}{\rho} \frac{\partial p}{\partial x} + \frac{\nu}{\rho} \Delta v_x \\
\frac{\partial v_y}{\partial t} + \vec{v} \cdot \vec{\nabla} \, v_y &= – \frac{1}{\rho} \frac{\partial p}{\partial y} + \frac{\mu}{\rho} \Delta v_y \\
\vec{\nabla} \cdot \vec{v} &= 0 \end{aligned}$$
If we take $\partial/\partial y$ of the first equation and subtract $\partial/\partial x$ of the second, the pressure is eliminated:
$$\frac{\partial \phi}{\partial t} + \frac{\partial \psi}{\partial y} \frac{\partial \phi}{\partial x} – \frac{\partial \psi}{\partial x} \frac{\partial \phi}{\partial y} = \frac{\mu}{\rho} \Delta \phi$$
where $\psi$ is the stream function, $\phi = \Delta \psi$ and velocities follows:
$$v_x = \frac{\partial \psi}{\partial y} \hspace{50pt} v_y = -\frac{\partial \psi}{\partial x}$$
The $\vec{\nabla} \cdot \vec{v} = 0$ equation is satisfied provided the stream function is smooth.
I have an idea how to solve this numerically: in each step we have $\phi$, we can solve the Poisson equation $\Delta \psi = \phi$. Then we make the time step $t \to t + \mathrm{d} t$:
$$\phi_{new} = \phi_{old} + \mathrm{d} t \left( whatever \right)$$
provided we can calculate whatever is in the brackets. Here comes my question: the no-slip boundary condition
$$\left. \vec{v} \right|_{boundary} = 0$$
should be somehow reflected to the function $\psi$, but how? I can force the Dirichlet boundary condition along some normal $\vec{n}$:
$$\left. \vec{n} \cdot \vec{\nabla} \psi \right|_{boundary} = 0$$
but I don't think I can force both $\psi^\prime_x$ and $\psi^\prime_y$ to be zero (or some constant, or whatever I want it to be).
So my question is: how do people usually solve this in this formulation? If I only have the no-penetration condition, it's easy, because that's exactly the Dirichlet condition (where $\vec{n}$ is the vector parallel to the boundary), but no-slip requires both normal and parallel components of gradient to be zero. Is there some subtlety I'm missing?
Moreover, let's imagine the following boundary: square $\left\langle 0, 1 \right\rangle^2$ with a circular hole of some radius less than $1/2$. I understand that I can somehow force the condition that the square boundary is not penetrable, that is, the stream function is constant along the square boundary, let's say $C_1$. Moreover, the circular hole is not penetrable as well, so the stream function is constant along the circle too, let's say $C_2$. Is it true, that $C_1 = C_2$? If not, how can I determine the difference $C_1 – C_2$? (I suspect that the stream function can be adjusted by any global constant, as we always deal with it's derivatives, so what matters is the difference $C_1 – C_2$, not those values themselves)
Is it worth to solve Navier-Stokes equations in 2D like this? It involves a lot of derivatives (a possible source of numerical errors), but I have always had a problem with the original formulation: I couldn't work out the numerical scheme to get the pressure in the next time step (as the equations don't contain any $\partial p/\partial t$ term).
Thank you.
Best Answer
Let's suppose that the boundary is the x-axis. So along the boundary, the stream function is constant. So, $$\frac{\partial^2\psi}{\partial x^2}=0$$ And from the no-slip boundary condition, $$\frac{\partial \psi}{\partial y}=0$$ This can be used to establish the 2nd order finite difference approximation to the value of the vorticity at the boundary:$$\omega=\frac{\partial^2 \psi}{\partial x^2}+\frac{\partial^2 \psi}{\partial y^2}=\frac{2(\psi(I,1)-\psi(I,0))}{(\Delta y)^2}=\omega(I,0)$$ This is used for the boundary condition on $\omega$.
ADDENDUM
First of all, the stream function is known (and constant) at the solid boundaries (because the solid boundaries are stream lines). So it doesn't have to be solved for. It is determined up to an arbitrary constant, and can thus be taken to be zero at one of the boundaries. At the other boundary, the stream function is equal to the volumetric throughput rate per unit width of channel (which is typically known). So you don't need to solve for the stream function at the solid boundaries.
If the tangential derivative $\partial \psi/\partial x=0$ at all locations along the boundary, it's second partial with respect to x must also be equal to zero. This, of course, all follows from the fact that $\psi$ is constant at the boundary.
To integrate the vorticity equation, you need a boundary condition on the vorticity (or at least a 2nd order finite difference approximation to a boundary condition). Just because $\partial \psi/\partial y=0$ does not mean the the second partial of $\psi$ with respect to y is equal to zero at the boundary; this would imply that the vorticity at the boundary is equal to zero, which we know is not correct.
The variable I in the relationships refers to the I'th x grid point. So, back to the boundary condition on vorticity: We have shown so far that, at the boundary, $$\omega=\frac{\partial^2 \psi}{\partial y^2}$$ subject to the constraint that $\partial \psi/\partial y=0$. If we represent these two conditions in 2nd order finite difference form, we obtain: $$\omega(I,0)=\frac{\psi(I,1)-2\psi(I,0)+\psi(I,-1)}{(\Delta y)^2}$$and$$\frac{\psi(I,1)-\psi(I,-1)}{2\Delta y}=0$$If we combine these two finite difference equations, we obtain a 2nd order finite difference approximation to the value of the vorticity at the boundary: $$\omega(I,0)=\frac{2(\psi(I,1)-\psi(I,0))}{(\Delta y)^2}$$ I've successfully used this approach to solving these equations many times.