Stream function formulation for Navier-Stokes equations in 2D - boundary conditions

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.

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$.

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.