Created time
Jul 10, 2025 11:10 AM
type
status
date
slug
summary
tags
category
icon
password
Last edited time
Jul 10, 2025 11:14 AM
前言：
Preface:
参考：感谢原作者们的无私引路和宝贵工作。
前置： OpenFOAM开发编程基础06 开发库 | 𝓐𝓮𝓻𝓸𝓼𝓪𝓷𝓭 (aerosand.cn)
通过前面总共七篇讨论，我们已经大概熟悉了 OpenFOAM 中数值应用的架构，熟悉了数值计算的关键元素（时间、网格、场），也熟悉了一些编程的语法、一些工作流程。下面我们使用所掌握的工具，尝试写一个最简单的 OpenFOAM 求解器（不涉及数值算法）。
计算流体力学的通用基本方程为
$$\frac{\partial}{\partial t}(\rho \phi) + \nabla \cdot (\rho U\phi) = \nabla\cdot(\Gamma\nabla\phi) + S_{\phi}$$
为了避免未知量和 OpenFOAM 中的质量通量
phi 误会，我们使用 A 表示待求物理场。下面，我们通过 OpenFOAM 编程求解一个不可压定常无源输运问题，控制方程如下
$$\nabla\cdot(UA) = \nabla\cdot(\Gamma\nabla A)$$
建立本文项目文件夹
应用准备
脚本和说明
略。
求解器
createFields.H
文件
/userApp/createFields.H主源码
方程的构建
上面的主源码咋一看好像也不难理解，但仔细想一想，读者可能会忍不住好奇，方程的构建是怎么样的？
需要强调的是，想要完全理解这些问题，毫无疑问还是要从有限体积法基础开始，深入到矩阵计算中去。强烈建议没有阅读过有限体积法基础的同学还是快速找补一下，避免基础不稳困惑加倍，参考 CFD理论基础02 有限体积法初步 | 𝓐𝓮𝓻𝓸𝓼𝓪𝓷𝓭 (aerosand.cn)。
我们看一下代码中的方程构建，可以按下面方式构建这个方程（结果和上述的主源码没有什么区别）
对流项
对流项的构造是
fvm::div(phi,A) ，这是一个函数，我们查找这个函数为了帮助理解，从代码声明
fvmDiv.H 中摘取几处代码简单讨论如下（不建议继续深究代码，代码挖的太深无益于主要学习的推进）我们了解了此函数的形参类型，明白返回的是
fvMatrix 类型的变量，数学上对应的是和待求量相关的矩阵。扩散项
扩散项的构造是
fvm::laplacian(gamma,A) ，这也是一个函数，我们查找这个函数同样的，摘部分代码如下
我们知道该函数同样返回
fvMatrix 类型，对应数学上和待求量相关的矩阵。构建
离散方程最后形式为
$$Ax = b$$
对流项和扩散项的矩阵相加，构成了待求左侧项，因为没有源项，所以右侧项置空为零，对应的代码就是
求解
构建完成后，我们使用
solve 函数求解线性代数系统暂时不用翻它的源代码。
其他疑惑
为什么散度计算中，使用的是
phi 而不是公式中的 U 呢？此外， fvm 和常见到的 fvc 有什么区别呢？简单理解的话， 求解器中的
fvm 离散所有对系数矩阵有贡献的项（隐式离散），相对比的，fvc 离散和系数矩阵无关的项（显式离散）。但是这么轻飘飘的一句话解释，其实很难让初学者理解。为了讨论上面两个问题，我们需要回到控制方程的偏微分方程数值求解上，从有限体积法开始，一步一步推导，直到矩阵求解，才可以真正明白为什么要使用通量进行计算，为什么会有隐式离散和显式离散的区分。这里先不深究，参见 CFD 基础系列讨论和 OpenFOAM 求解器系列讨论（
ofss）。好在求解语句的代码语法已经十分贴近数学公式，即使带着这些困惑也不要紧，我们仍然可以编写简单的求解器。
编译
调试算例
需要对调试算例做一些修改，以便计算进行。
system/blockMeshDict无需修改，我们仍然使用原几何模型和几何网格
system/controlDict无需修改，我们仍然使用求解控制的设置（步长步数等）
decomposedParDict和PDRblockMeshDict不用管，初学阶段暂时用不到
fvSchemes
- 因为求解涉及到散度和拉普拉斯，所以需要指明两个计算项的离散格式
至于这些格式是在计算过程的哪里起作用的，怎么在起作用的，很难简单解释清楚，还是要参见 CFD 基础系列讨论和 ofss 系列讨论
fvSolution
只计算了 场 A
- 因为求解器只计算了场 A，所以只设置 A 的代数求解器
初始条件
在调试算例的
p 文件上修改文件名称得到 A 文件，内容修改如下U 文件保持不变。运行
终端显式如下
调试算例下多了新的一个时间步文件夹，其中有场 A 的计算结果。
通过 paraview 打开显式结果（以后不再赘述命令）
可以看到一个时间步的计算结果。
没有时间步动画显示，这是当然的，因为这是定常问题，所以我们只计算一个时间步（仅仅是为了不覆盖初始条件）。
小结
我们终于写了第一个求解器，但是还有一些疑惑涉及到有限体积法的数值细节，其他部分的代码和以前讨论里的应用其实没有多少区别。
我们记录一下可能有的疑惑：
fvSchemes文件中对离散格式的指定是怎么反映到求解计算中的？
fvSolution文件中对求解方法的指定是怎么反映到求解计算中的？
这两个问题可以当成黑箱去设置调用，具体细节暂时放一放，需要参考 CFD 基础系列和
ofss 系列，不要担心，到时候会一个一个解决。- 代码实现中为什么使用质量通量而不是直接使用速度？
- 所谓最后
solve的线性代数系统在数学上是怎么组建的？fvm和fvc在其中有什么区别？
这两个可能的疑惑请记下来并在后续学习中逐渐解决。如前所述，好在 OpenFOAM 的语法简明直观，即使带着这些困惑也可以实现简单的求解器。
我们下一篇讨论再求解一个简单的非定常问题。
欢迎留言讨论，反馈建议和意见，赞助打赏。
Feel free to leave comments, feedback, suggestions, opinions, and donations..
Loading...
