FiPy 2D Navier Stokes 实现：一个方向的导数问题

问题描述

我正在尝试实现 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",value=-0.)
v = CellVariable(name='veLocity',rank = 1.,value=[vx,vy])
p = CellVariable(name = 'pressure',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,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',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').

有人可以启发我吗？提前致谢！

解决方法

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

• 使用 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",value=-0.,hasOld=True)

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

p = CellVariable(name = 'pressure',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')