FiPy 2D Navier Stokes 实现：单向导数问题

Question

我正在尝试实现 https://nbviewer.jupyter.org/github/barbagroup/CFDPython/blob/master/lessons/14_Step_11.ipynb 中的示例，但我遇到了一些问题。我认为主要问题是我在边界条件以及方程中定义项方面遇到了一些问题。

偏微分方程是：

初始边界条件为：

我知道我的变量是速度矢量 (u, v) 和压力 (p) 的两个分量。按照示例并使用 FiPy，我将 PDE 编码如下：

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

import fipy
from fipy import *

# Geometry
Lx = 2   # meters
Ly = 2    # meters
nx = 41  # nodes
ny = 41

# Build the mesh:
mesh = Grid2D(Lx=Lx, Ly = Ly, nx=nx, ny=ny)

# Main variable and initial conditions
vx = CellVariable(name="x-velocity", 
                 mesh=mesh,
                 value=0.)
vy = CellVariable(name="y-velocity",
                 mesh=mesh,
                 value=-0.)
v = CellVariable(name='velocity',
                 mesh=mesh, rank = 1.,
                 value=[vx, vy])
p = CellVariable(name = 'pressure',
                 mesh=mesh,
                 value=0.0)

# Boundary conditions
vx.constrain(0, where=mesh.facesBottom)
vy.constrain(0, where=mesh.facesBottom)
vx.constrain(1, where=mesh.facesTop)
vy.constrain(0, where=mesh.facesTop)
p.grad.dot([0, 1]).constrain(0, where=mesh.facesBottom)  # <---- Can this be implemented like this?
p.constrain(0, where=mesh.facesTop)
p.grad.dot([1, 0]).constrain(0, where=mesh.facesLeft)
p.grad.dot([1, 0]).constrain(0, where=mesh.facesRight)

#Equations
nu = 0.1 #
rho = 1.
F = 0.
# Equation definition
eqvx = (TransientTerm(var = vx) == DiffusionTerm(coeff=nu, var = vx) - ConvectionTerm(coeff=v, var=vx) - ConvectionTerm(coeff= [[(1/rho)], [0]], var=p) + F)
eqvy = (TransientTerm(var = vy) == DiffusionTerm(coeff=nu, var = vy) - ConvectionTerm(coeff=v, var=vy) - ConvectionTerm(coeff= [[0], [(1/rho)]], var=p))
eqp = (DiffusionTerm(coeff=1., var = p) == -1 * rho * (vx.grad.dot([1, 0]) ** 2 + (2 * vx.grad.dot([0, 1]) * vy.grad.dot([1, 0])) + vy.grad.dot([0, 1]) ** 2))
eqv = eqvx & eqvy & eqp

# PDESolver hyperparameters

dt = 1e-3   # (s) It should be lower than 0.9 * dx ** 2 / (2 * D)
steps = 100  #
print('Total time: {} seconds'.format(dt*steps))

# Plotting initial conditions

# Solve
vxevol = [np.array(vx.value.reshape((ny, nx)))]
vyevol = [np.array(vy.value.reshape((ny, nx)))]

for step in range(steps):
  eqv.solve(dt=dt)
  v = CellVariable(name='velocity', mesh=mesh, value=[vx, vy])

  sol1 = np.array(vx.value.reshape((ny, nx)))
  sol2 = np.array(vy.value.reshape((ny, nx)))
  vxevol.append(sol1)
  vyevol.append(sol2)

我在第 100 步的结果是这个（与示例中给出的解决方案不一致）：

我认为主要问题在于定义一个特定维度的边界条件（例如dp/dx = 1，dp/dy = 0），以及变量在一维中的导数方程式（在代码中，'eqp'）。

谁能赐教一下？提前致谢！

Answer 1

这是我认为的改进版本。至少结果看起来比较合理。主要变化如下。

使用 linked CFD notebook 中概述的技巧以及压力方程中时间步长的发散。
将矢量速度v更改为面变量，以便我们可以直接使用.divergence。当然可以清理，但是是不同的离散化。我不知道哪个更有效。
修复边界条件。我不确定 p.grad.dot[].constrain 是否在做任何明智的事情。无论如何，零梯度不需要它们，因为这是默认值。
没有求解一个矩阵中的所有方程。一旦您有信心正确地单独解决问题并且您有一个基准可以检查，那么最好这样做。
速度矢量变量在每一步都被重新创建，这意味着这对方程没有影响。 v 现在在循环中明确更新。
在将压力梯度添加到动量方程时不使用 ConvectionTerm。 ConvectionTerm 进行了奇怪的加权，并不是一个简单的区别。从长远来看运行它可能很好用，但在调试时不适合。

这是代码。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from fipy import (
    CellVariable,
    ConvectionTerm,
    Grid2D,
    FaceVariable,
    TransientTerm,
    DiffusionTerm,
    CentralDifferenceConvectionTerm,
    Viewer
)

from fipy.viewers.matplotlibViewer.matplotlibStreamViewer import MatplotlibStreamViewer

# Geometry
Lx = 2   # meters
Ly = 2    # meters
nx = 41  # nodes
ny = 41

# Build the mesh:
mesh = Grid2D(Lx=Lx, Ly = Ly, nx=nx, ny=ny)

# Main variable and initial conditions
vx = CellVariable(name="x-velocity",
                  mesh=mesh,
                  value=0.,
                  hasOld=True)
vy = CellVariable(name="y-velocity",
                  mesh=mesh,
                  value=-0.,
                  hasOld=True)

v = FaceVariable(name='velocity',
                 mesh=mesh, rank = 1)

p = CellVariable(name = 'pressure',
                 mesh=mesh,
                 value=0.0,
                 hasOld=True)

# Boundary conditions

# top

vx.constrain(1, where=mesh.facesTop)
vy.constrain(0, where=mesh.facesTop)
p.constrain(0, where=mesh.facesTop)

# left

vx.constrain(0, where=mesh.facesLeft)
vy.constrain(0, where=mesh.facesLeft)

# right

vx.constrain(0, where=mesh.facesRight)
vy.constrain(0, where=mesh.facesRight)

# bottom

vx.constrain(0, where=mesh.facesBottom)
vy.constrain(0, where=mesh.facesBottom)

#Equations
nu = 0.1
rho = 1.

dt = 0.01

# Equation definition
eqvx = (TransientTerm() == DiffusionTerm(nu) - ConvectionTerm(v) - (1 / rho) * p.grad[0])
eqvy = (TransientTerm() == DiffusionTerm(nu) - ConvectionTerm(v) - (1 / rho) * p.grad[1])
eqp = (DiffusionTerm(1.) == -1 * rho * (v.divergence**2 - v.divergence / dt))


steps = 200
sweeps = 2
print('Total time: {} seconds'.format(dt*steps))
total_time = 0.0

viewer = MatplotlibStreamViewer(v)
pviewer = Viewer(p)
vxviewer = Viewer(vx)
vyviewer = Viewer(vy)

for step in range(steps):
    vx.updateOld()
    vy.updateOld()
    p.updateOld()

    for sweep in range(sweeps):
        res_p = eqp.sweep(var=p, dt=dt)
        res0 = eqvx.sweep(var=vx, dt=dt)
        res1 = eqvy.sweep(var=vy, dt=dt)

        print(f"step: {step}, sweep: {sweep}, res_p: {res_p}, res0: {res0}, res1: {res1}")

        v[0, :] = vx.faceValue
        v[1, :] = vy.faceValue

    viewer.plot()
    pviewer.plot()
    vxviewer.plot()
    vyviewer.plot()

input('end')

FiPy 2D Navier Stokes 实现：单向导数问题

FiPy 2D Navier Stokes implementation: problem with derivatives in one direction

fipy