NS 方程

4. NS 方程#

4.1. 方腔流#

Navier-Stocks 方程:

(4.1)#\[\begin{equation} \left\{ \begin{aligned} &\partial_t u - \mu\Delta u + (u\cdot\nabla)u + \nabla p = f, && {\rm in} \quad \Omega\times(0, T]\\ &\nabla\cdot u = 0, && {\rm in} \quad \Omega\times(0, T] \end{aligned} \right. \end{equation}\]

边界条件除上边界外为滑移边界条件, 上边界速度为 \(u=(1, 0)\). 初边值条件 0.

from firedrake import *

mu = 1
T = 0.25

N_S = 16
N_T = 128

tau = T/N_T
h = 1/N_S

mesh = RectangleMesh(N_S, N_S, 1, 1)

u_0 = as_vector((0, 0))
u_D = as_vector((1, 0))
f = as_vector([0, 0])

firedrake:WARNING OMP_NUM_THREADS is not set or is set to a value greater than 1, we suggest setting OMP_NUM_THREADS=1 to improve performance

4.2. 函数空间#

采用 MINI 元, 即 P1 \(\times\) P1b.

P1b 由 P1 加上 Bubble 组成.

NodalEnrichedElement, EnrichedElement

VectorFunctionSpace 构造向量空间

cell = mesh.ufl_cell()
tdim = cell.topological_dimension()

# Mini element: P1 X P1b
P1 = FiniteElement("CG", cell, 1)
B = FiniteElement("B", cell, tdim+1)
P1b = P1 + B # or P1b = NodalEnrichedElement(P1, B)

V_u = VectorFunctionSpace(mesh, P1b)
V_p = FunctionSpace(mesh, "CG", 1)
V = MixedFunctionSpace([V_u, V_p])

4.3. 弱形式#

线性格式:

(4.2)#\[\begin{equation} \left\{ \begin{aligned} &\frac{1}{\tau}(u^n - u^{n-1}, v) + \mu(\nabla u^n, \nabla v) + \frac{1}{2}(((u^{n-1}\cdot\nabla)u^n, v) - ((u^{n-1}\cdot\nabla)v, u^n)) - (p^n, \nabla\cdot v) = (f^n, v)\\ &(q, \nabla\cdot u^n) = 0 \end{aligned} \right. \end{equation}\]

非线性格式

(4.3)#\[\begin{equation} \left\{ \begin{aligned} &\frac{1}{\tau}(u^n - u^{n-1}, v) + \mu(\nabla u^n, \nabla v) + \frac{1}{2}(((u^n\cdot\nabla)u^n, v) - ((u^n\cdot\nabla)v, u^n)) - (p^n, \nabla\cdot v) = (f^n, v)\\ &(q, \nabla\cdot u^n) = 0 \end{aligned} \right. \end{equation}\]

  • TrialFunctions, TestFunctions:

    以 tuple 返回函数空间中的试验/测试函数,

    主要用于 MixedFunctionSpace.

  • split

    • split: 以索引的方式获取 MixedFunctionSpace 中函数的分量 (保留 UFL 关联信息, 用于定义变分形式)

由于该问题是非线性问题, 我们打算用 NonlinearVariationalSolver 进行求解, 所以下面定义 w 使用了 Function 而不是 TrialFunction/TrialFunctions.

w = Function(V) # u and p
u, p = split(w)

v, q = TestFunctions(V)

w_nm1 = Function(V)
u_nm1, p_nm1 = w_nm1.subfunctions
u_nm1.rename('u_h') # for visualization in paraview
p_nm1.rename('p_h')

F = (
      Constant(1/tau)*inner(u - u_nm1, v)*dx
    + mu*inner(grad(u), grad(v))*dx
    + 1/2*inner(dot(grad(u), u/2), v)*dx
    - 1/2*inner(dot(grad(v), u/2), u)*dx
    - p*div(v)*dx
    + div(u)*q*dx
    - inner(f, v)*dx(domain=mesh)
)

4.4. 定义 Solver#

类似于纯 Neumann 问题, 我们将使用 nullspace 参数.

注意下面混合空间中, 边界条件和 nullspace 的定义.

bc1 = DirichletBC(V.sub(0).sub(0), 0, (1, 2))
bc2 = DirichletBC(V.sub(0).sub(1), 0, (3, 4))
bc3 = DirichletBC(V.sub(0).sub(0), 1, 4)  # upper boundary

nullspace = MixedVectorSpaceBasis(V, [V.sub(0), VectorSpaceBasis(constant=True, comm=mesh.comm)])

problem = NonlinearVariationalProblem(F, w, bcs=[bc1, bc2, bc3])  # F = 0
solver = NonlinearVariationalSolver(problem,
                                    options_prefix='ns',
                                    solver_parameters=None, # {'snes_converged_reason': None, 'snes_max_it': 100},
                                    nullspace=nullspace
                                   )

4.5. 时间循环和保存结果到 pvd 文件#

u_, p_ = w.subfunctions

output = VTKFile('pvd/ns-equation.pvd')

u_nm1.assign(0)
output.write(u_nm1, p_nm1, time=0)

for i in range(N_T):
    t = tau*(i+1)
    
    solver.solve()
    
    u_nm1.assign(u_)
    p_nm1.assign(p_)

    if (i+1)%32 == 0:
        output.write(u_nm1, p_nm1, time=t)

4.6. TODO: 随时间变系数问题#

from firedrake import *
c = Constant(0)

for i in range(5):
    c.assign(i*0.1)
    print(f"i = {i}, c = {c}")

4.7. ParaView 可视化计算结果#

Pipeline 和 Filter

4.7.1. 二维结果 (surf 图)#

Filter: Wrap by scalar

4.7.2. 选择部分区域显示#

View -> Find Data

4.7.3. 并行数据显示各进程区域#

Filter -> Connectivity