## Help: IMEX Scheme for Incompressible Navier-Stokes Fails Without Incremental Form

General
1

Hi all,

I’m trying to solve the incompressible Navier-Stokes equations using an IMEX (implicit-explicit) time-stepping scheme in NGSolve, and I’m encountering issues.

Instead of using the incremental form recommended in the official tutorials, I attempted to directly solve the full system at each time step without explicitly splitting the nonlinear convection term. In other words, I assembled the full system as:

Mun+1+dt⋅(Aun+1+C(un+1))=rhsM u^{n+1} + dt \cdot \left( A u^{n+1} + C(u^{n+1}) \right) = rhsMun+1+dt⋅(Aun+1+C(un+1))=rhs

where C(u)C(u)C(u) denotes the convection term. I also tried moving the convection term to the right-hand side explicitly, but the system either fails to converge or produces unstable, oscillatory results.

In contrast, the incremental approach used in NGSolve examples—e.g., forming the right-hand side as (r.mat - dt * conv.mat) * gfu.vec—seems much more stable.

:question: I’d like to ask:

  • Is it necessary to use the incremental form when using IMEX schemes with Navier-Stokes?
  • Is there a recommended way to linearize or stabilize the full system if solving it directly?
  • Are there any working examples that implement full system assembly (including convection) successfully?

Any suggestions or insights would be greatly appreciated.

Best regards,
777

from ngsolve import *
from ngsolve.webgui import Draw
from netgen.occ import *
from netgen.webgui import Draw as DrawGeo
import matplotlib
import numpy as np
import matplotlib.pyplot as plt

shape = Rectangle(2,0.41).Circle(0.2,0.2,0.05).Reverse().Face()  # 先绘制一个矩形，再在(0.2,0.2)位置绘制圆，并Reverse使其为interface
shape.edges.name = "cyl" # 圆形表面
shape.edges.Max(X).name = "outlet"
shape.edges.Min(X).name = "inlet"
shape.edges.Max(Y).name = "wall"
shape.edges.Min(Y).name = "wall"
mesh = Mesh(OCCGeometry(shape, dim=2).GenerateMesh(maxh = 0.07)).Curve(3) # 三阶多项式拟合曲线
n = specialcf.normal(mesh.dim)

# # 绘制物理域及计算域
# DrawGeo(shape); 
# Draw(mesh); # OCCGeometry会根据绘制图像复杂性生成不同粗细的网格拟合（如圆附近网格会更细）

nu = 0.001
k = 2 # 空间阶数
dt = 0.001 # 定义时间步长
t = 0
i = 0
time_vals, drag_x_vals, drag_y_vals = [],[],[]

V = VectorH1(mesh, order = k, dirichlet = "wall|cyl|inlet")
Q = H1(mesh, order = k-1) # 未指明的边界默认为自然边界条件（在弱形式中会体现为零Neumann边界条件）
X = V * Q # 定义速度场与压力场的耦合空间

gfu = GridFunction(X) # 设置GridFunction作为解域，用于存储解作为后续初值使用
velocity, pressure = gfu.components # X由V，P组成，gfu.components即为[V，P]
uin = CF((1.5*4*y*(0.41-y)/(0.41**2),0)) # CF为系数矩阵
velocity.Set(uin, definedon = mesh.Boundaries("inlet"))

drag_x_test = GridFunction(X)
drag_x_test.components[0].Set(CoefficientFunction((-20.0,0)), definedon=mesh.Boundaries("cyl"))
drag_y_test = GridFunction(X)
drag_y_test.components[0].Set(CoefficientFunction((0,-20.0)), definedon=mesh.Boundaries("cyl"))

# # 分别绘制速度场和压力场处置
# Draw(velocity);
# Draw(pressure);

# 设置Trial和Test函数
(u,p),(v,q) = X.TnT()

a = BilinearForm((nu*InnerProduct(grad(u),grad(v)) - div(v)*p - div(u)*q)*dx).Assemble()
f = LinearForm(X).Assemble() # 默认为0

res = f.vec - a.mat*gfu.vec
inv_stokes = a.mat.Inverse(X.FreeDofs())
gfu.vec.data += inv_stokes * res

# # 绘制stokes的解（NS的初值）
# sceneu = Draw(velocity,mesh,"u")
# scenep = Draw(pressure,mesh,"p")

mstar = BilinearForm(InnerProduct(u,v)*dx + dt*(nu*InnerProduct(grad(u),grad(v)) - div(v)*p - div(u)*q - (Grad(u) * u) * v)*dx).Assemble()
inv = mstar.mat.Inverse(X.FreeDofs())
conv = BilinearForm(X, nonassemble = True) # 对流项并非双线性，此处仅是利用Bilinear的函数性质运算（即作为函数），适用于增量处理情况，也可以通过牛顿迭代等方法处理为线性项并在每次迭代中Assemble处理
conv += (Grad(u) * u) * v * dx
r = BilinearForm(InnerProduct(u,v)*dx).Assemble()
f = LinearForm(X)

t = 0
i = 0
tend = 1
gfut = GridFunction(gfu.space,multidim = 0)
gfut.AddMultiDimComponent(gfu.vec)
while t < tend-0.5*dt:
    print ("\rt=", t, end="")
    conv.Assemble()
    f = (r.mat - dt * conv.mat) * gfu.vec
    # f.vec.data -= Integrate(0.5 * InnerProduct(velocity, velocity) * InnerProduct(velocity, n), mesh, BND)
    gfu.vec.data = inv * f
    
    t = t + dt
    i += 1
    time_vals.append(t)
    drag_x_vals.append(InnerProduct(res, drag_x_test.vec))
    drag_y_vals.append(InnerProduct(res, drag_y_test.vec))
    if i % 50 == 0:
        gfut.AddMultiDimComponent(gfu.vec)

Draw (gfut.components[0], mesh, interpolate_multidim=True, animate=True,
      min=0, max=1.9, autoscale=False, vectors = True);
Draw (gfut.components[1], mesh, interpolate_multidim=True, animate=True,
      min=-0.5, max=1, autoscale=False);

# Plot drag force over time
plt.plot(time_vals, drag_x_vals)
plt.xlabel('time'); plt.ylabel('drag'); plt.title('drag'); plt.grid(True)
plt.show()   
# Plot lift force over time
plt.plot(time_vals, drag_y_vals)
plt.xlabel('time'); plt.ylabel('lift'); plt.title('lift'); plt.grid(True)
plt.show()