基于FiPy实现二维Burgers方程：如何完成方程代码编写？

有资料指出FiPy在处理双曲型PDE方面表现不佳，但相关内容发布于2016年，情况或许已有所改变。我希望实现二维Burgers方程：
$$
\begin{cases}
u_t + (u \partial_x + v \partial_y) u = -\nu \nabla^2 u \
v_t + (u \partial_x + v \partial_y) v = -\nu \nabla^2 v
\end{cases}
$$
其中$u$和$v$是流体流动$\boldsymbol{U}=(u,v)$的x和y分量。为简化操作，我利用链式法则将方程改写为：
$$
\begin{cases}
u_t + \nabla \cdot (\boldsymbol{U} u) = (\nabla \cdot \boldsymbol{U}) u - \nu \nabla^2 u \
v_t + \nabla \cdot (\boldsymbol{U} v) = (\nabla \cdot \boldsymbol{U}) v - \nu \nabla^2 v
\end{cases}
$$
对应到FiPy中的代码框架如下：

u = CellVariable(name='u', mesh=m)
v = CellVariable(name='v', mesh=m)
U = CellVariable(name='U', mesh=m, value=(u, v))
eqn1 = TransientTerm(var=u) + ConvectionTerm(coeff=U, var=u) == DiffusionTerm(coeff=-nu, var=u) + # Not Sure what
eqn2 = TransientTerm(var=v) + ConvectionTerm(coeff=U, var=v) == DiffusionTerm(coeff=-nu, var=v) + # Not Sure what
burgersEq2D = eqn1 & eqn2

现咨询：如何完成上述代码中eqn1和eqn2的编写？

你需要在等式右侧补充$(\nabla \cdot \boldsymbol{U}) u$和$(\nabla \cdot \boldsymbol{U}) v$对应的FiPy表达式，完整代码如下：

import fipy as fp

# 假设已定义网格m和粘度参数nu
u = fp.CellVariable(name='u', mesh=m)
v = fp.CellVariable(name='v', mesh=m)
U = fp.CellVariable(name='U', mesh=m, value=(u, v))
# 计算速度场的散度 ∇·U
div_U = U.divergence

# 完成二维Burgers方程的定义
eqn1 = fp.TransientTerm(var=u) + fp.ConvectionTerm(coeff=U, var=u) == fp.DiffusionTerm(coeff=-nu, var=u) + div_U * u
eqn2 = fp.TransientTerm(var=v) + fp.ConvectionTerm(coeff=U, var=v) == fp.DiffusionTerm(coeff=-nu, var=v) + div_U * v

burgersEq2D = eqn1 & eqn2

核心说明 #

• U.divergence直接计算速度场$\boldsymbol{U}$的散度$\nabla \cdot \boldsymbol{U}$，返回的CellVariable可直接与u、v做乘法，对应改写后方程右侧的$(\nabla \cdot \boldsymbol{U}) u$和$(\nabla \cdot \boldsymbol{U}) v$项。
• ConvectionTerm(coeff=U, var=u)正好对应$\nabla \cdot (\boldsymbol{U} u)$，和改写后的方程左侧完全匹配。
• 针对双曲型PDE的稳定性优化：FiPy支持通过scheme参数为对流项指定迎风格式（如upwind），能有效抑制数值振荡，适合处理含间断解的问题，示例：

fp.ConvectionTerm(coeff=U, var=u, scheme="upwind")

内容的提问来源于stack exchange，提问作者user3166083