The Navier-Stokes equation of an incompressible fluid is
$$\begin{align} \nabla \cdot \mathbf{u} & = 0 \tag{1}\label{1}\\ \dfrac{\partial \mathbf{u}}{\partial t} + \left(\mathbf{u} \cdot \nabla\right) \mathbf{u} & = -\nabla p + \mu \cdot \nabla^2 \mathbf{u}\label{2}\tag{2}\end{align}$$
I want to solve the Lid-driven cavity problem, which is a $2$D problem, in a $(1 \times 1)$ square domain.
I read the problems of a numeric fluid's algorithm come mainly from:
- There's no equation for the pressure. One way to bypass it is by using the Projection method.
- When computing the next step (using finite differences for example), $\mathbf{u}$ must satisfy \eqref{1}.
To ignore these problems, I thought about
- Using the variable $s$ (stream line), which satisfies automatically \eqref{1}
- Taking the rotational of \eqref{2}, to get rid of the pressure.
$$u = \dfrac{\partial s}{\partial y} \ \ \ \ \ \ \ \ \ \ \ \ v = -\dfrac{ \partial s}{ \partial x}$$
Then, I get the formulation \eqref{5}! Which is only one equation, only one variable! Should be simpler to solve!
$$\boxed{\nabla^2 \dfrac{\partial s}{\partial t} = \underbrace{\mu \nabla^4 s}_{\text{diffusion}} + \underbrace{\left(\dfrac{\partial s}{\partial x} \dfrac{\partial^3 s}{\partial y^3}-\dfrac{\partial s}{\partial y} \dfrac{\partial^3 s}{\partial x\partial y^2}+\dfrac{\partial s}{\partial x} \dfrac{\partial^3 s}{\partial x^2\partial y}+\dfrac{\partial s}{\partial y} \dfrac{\partial^3 s}{\partial x^3}\right)}_{\text{advection}}}\label{5}\tag{5}$$
The fourth derivative is ugly, but with appropriate elements, it should work.
I tried to implement a numeric solution using FEM on the equation \eqref{5}, but I got terrible results. The oscillates too much in the interior compared with FDM (finite differences method) using the projection method.
Question: Does anyone know why this approach doesn't work? Is it a formulation problem or a numeric problem?
More details:
My domain is a square $ \Omega = \left[0, \ 1\right] \times \left[0, \ 1 \right]$.
The boundary conditions (like shown here) are
$$s(x, y) = 0 \ \ \text{on} \ \partial \Omega$$ $$\dfrac{\partial s}{\partial x} (0, y) = \dfrac{\partial s}{\partial x} (1, y) = \dfrac{\partial s}{\partial y} (x, 0) = 0$$ $$\dfrac{\partial s}{\partial y} (x, 1) = U_{top}(x)$$
I used FEM only on the space, so I transform \eqref{5} into \eqref{6}, to use a time-integrator after that.
$$s(x, \ y, \ t) = \sum_{i} \varphi_i(x, \ y) \cdot S_{i}(t)$$ $$\mathbf{C} \cdot \dot{\mathbf{S}} = \underbrace{\mathbf{D} \cdot \mathbf{S}}_{\text{diffusion}} + \underbrace{\mathbf{S} \cdot \mathbf{A} \cdot \mathbf{S}}_{\text{advection}}\label{6}\tag{6}$$ $$\sum_{j} \mathbf{C}_{ij} \cdot \dot{\mathbf{S}}_{j} = \sum_{j} \mathbf{D}_{ij} \cdot \mathbf{S}_{j} + \sum_{k, \ j} \mathbf{S}_{k} \cdot \mathbf{A}_{kij} \cdot \mathbf{S}_j \ \ \ \ \ \ \forall i\label{7}\tag{7}$$
- $\mathbf{C}$ is a $2$D matrix computed by $$\mathbf{C}_{ij} \equiv \int_{\Omega} \psi_{i} \nabla^2 \varphi_{j} \ d\Omega$$
- $\mathbf{D}$ is a $2$D matrix computed by $$\mathbf{D}_{ij} \equiv \int_{\Omega} \psi_{i} \nabla^4 \varphi_{j} \ d\Omega$$
- $\mathbf{A}$ is a $3$D matrix, which expression is similar, but bigger. Let me know if you need.
I used $2$-order and $3$-order lagrange square elements, so the diffusion term wouldn't be zero. But the results were similar.
The maximum eigenvalue of $\mathbf{C}$ I found was near $1$, then I took the timestep $\Delta t$ of order $0.01$.
My initial state $s_0$ (time $0$) is the stationary state when $\mu = \dfrac{1}{\text{Re}}$ is very big:
$$\nabla^4 s_0 = 0$$
Which solution is here