Navier-Stokes treating convection with IMEX scheme

I’m looking to solve navier-Stokes equation with an IMEX method which gives a convective term on the form

c(u_{n+1}, u_n, v) = \int_{\Omega} \nabla u_{n+1} u_n \cdot v d\Omega

as this gives unconditional stability. I have tried to do this with a bilinear form which includes a gridfunction un. Below is is the setup of the spaces and convective part of the bilinear form:

V = VectorH1(mesh,order=k, dirichlet=“wall|cyl|inlet”)
Q = H1(mesh,order=k-1)
X = VQ
(u,p), (v,q) = X.TnT()
gfu = GridFunction(X)
velocity, pressure = gfu.components
uin = CoefficientFunction((1.5
velocity.Set(uin, definedon=mesh.Boundaries(“inlet”))
un = gfu.components[0]
conv = BilinearForm(X)
c = InnerProduct(Grad(u) * un, v)
conv += c

Assembly of conv and updating of un is done while time stepping. This method seems to work for quite a few time steps, until much of the data suddenly become nan values. Any help with figuring out how to do implement this scheme properly would be much appreciated!


if you want to trat the convection semi-implicitly, you need to at the c form to the BilinearForm with all the other implicit terms. As a consequence, you will need to reassemble the matrix and inverse in every time-step.

Best, Henry