# Navier-Stokes Stream function-vorticity formulation

Hey there!
I am trying to implement the 2D lid-driven cavity problem in the stream-function formulation.

Follows a minimal working example of what I coded until now:

from fenics import *
##############################################################################
## PDE Parameters Definition
n = 32
Re= 100
nu = 1.0/Re
mesh = UnitSquareMesh(n,n)
V = FunctionSpace(mesh,'CG',2)
## Boundary Conditions
def boundary(x, on_boundary):
return on_boundary
## Homogeneus Boundary Conditions
bc0 = DirichletBC(V, Constant(0.0), boundary)
## Defining PDE
psi = TrialFunction(V)
phi = TestFunction(V)

'''

Hi,
The articles that you proposed ain’t freely accessible, can write down the formula you’re referencing to?
For the boundary condition you want to impose Neumann bc?

Hi Yoann,
Yes, I want to impose Neumann BC!

Ok thanks,
I have some questions before answering. The equation (5.10) is your PDE? The \Delta is a Laplacian or \nabla ?
It’s important to have a clear formulation of your PDE. For example you can rewrite (5.10) with \nabla\psi instead of \partial_{x_1}\psi and \partial_{x_2}\psi.
Finally are you sure that the Neumann bc and the dirichlet bc are compatible?

Dear Yoann,

if I understand well, \Delta is the Laplacian (in another place of the paper (5.10) is called ‘‘Nonlinear Biharmonic Problem’’).

Yes, problem (5.10) with boundary conditions (6.7) and (6.8) is the problem I want to implement (I realize just now that it is slightly different from the problem in (6.9) ).

I am not an expert in PDE and for this reason, I decided to post this topic.

I believe that conditions (6.7) and (6.8) are not compatible as they are written. Since this should be equivalent to the lid-driven cavity problem I believe that (6.7) should be imposed just on three sides of the square and (6.8) should be imposed on one side of the square. I am not sure of this and I would appreciate knowing your opinion.

Thank you once more!

Stefano

You may also want to check the Biharmonic demo since the weak form 5.16, if you want to implement as it is, requires you to have your approximation in H^2(\Omega), which again would require at least C^1-continuous spaces. These are difficult to construct compared to C^0 functions. Even though FIAT as the celebrated Argyris triangle, as far as I know, it is not supported in ffc.

Dear bhaveshshrimali,
Yes, I think only FIAT supports the Argyris triangle, ffc doesn’t; so you cannot do any meaningful calculations using it. I have tested this in dolfin-2018.1.0; and I don’t think dolfin-2019.1.0 would either. You can check by typing FunctionSpace(mesh, "Argyris", 5).
Firedrake has full support for Argyris as I understand but even though it does, imposing Dirichlet boundary conditions might not be straightforward for C^1 continuous elements. You can check the paper, they have examples and code as supplements. If you are looking to use Argyris, you can also check a lightweight library scikit-fem. It also has an example on using Argyris.
I haven’t used tiGAR, but it does seem to support C^1 splines on quads/hex @kamensky might be the best person to help you with it.