I am investigating lid-driven cavity flow, demonstrated in the below diagram:
We have a square (two dimensional) domain, with fully Dirichlet conditions for the velocity and fully Neumann conditions for the pressure, with the extra (necessary to to get a unique solution for $p$) constraint that $p\big(t,(0,0)\big)=0$. The governing equations are the incompressible Navier-Stokes equations for isotropic Newtonian fluids:
$$\nabla\cdot \boldsymbol u=0 \tag{1A}$$
$$\partial_t\boldsymbol u+\nabla\cdot(\boldsymbol u\otimes\boldsymbol u)=-\nabla p+\nu\Delta\boldsymbol u\tag{1B}$$
The Reynolds number is defined as $\operatorname{Re}=\frac{UL}{\nu}$. I'll skip the details of the derivation as it is not relevant to the question but, if we recast the equation in non-dimensional form, and let $\operatorname{Re}\to 0$, we get the Stokes equation
$$\Delta\boldsymbol u=\nabla p$$
If we now introduce a streamfunction $\psi$ satisfying $\boldsymbol u=\nabla\times (\psi\boldsymbol e_3)$ and take the curl of the Stokes equation, we get the biharmonic equation,
$$\Delta^2\psi=0\tag{2}$$
Since
$$\boldsymbol u=\begin{bmatrix}\partial_y \psi \\ -\partial_x\psi\end{bmatrix}$$
And recalling our BCs on $\boldsymbol u$, we get the following BCs for $\psi$:
$$(\partial_y\psi)(x,0)=(\partial_y \psi)(0,y)=(\partial_y\psi)(1,y)=0~~,~~(\partial_y\psi)(x,1)=1 \\
(\partial_x \psi)(x,0)=(\partial_x \psi)(x,1)=(\partial_x \psi)(0,y)=(\partial_x \psi)(1,y)=0$$
Because we have a fourth order equation in two variables, these eight boundary conditions should be enough for a solution, I think? I am not sure.
My question:
Is there a good way to analytically or numerically solve $(2)$ with the prescribed BCs? I initially tried separation of variables but quickly saw that this was a dead end because of the mixed term $\frac{\partial^4\psi}{\partial x^2\partial y^2}$. I also tried shoving the equation into Mathematica's NDSolve, but it threw an error, saying "the order of the spatial derivatives may not exceed two", which puzzled me, because I have numerically solved the third-order equation $\partial_t\phi=\partial_x^3\phi$ in Mathematica before. I could possibly code up a purpose-built finite difference solver in Python or FORTRAN but not only would this be very time consuming, I would also be afraid of getting incorrect results. Can I recondition the above equation somehow so that Mathematica will accept it?
Thanks for any advice.
Best Answer
Yes, this exact problem has been studied using the vorticity-streamfunction approach, and the numerical solution has been implemented in Mathematica.