1.OpenFOAM中的源码数据结构-icoFoam为例【未完成】

OpenFOAM中的数据结构-icoFoam为例【未完成】

       撰写博客并整理思路确实能提高工作效率。我本想利用暑假进行一项关于 OpenFOAM 数据结构的源码深入研究，然而日程安排总是源码让我疲于应对各种事务，仿佛永远在追赶时间的源码脚步。近期，源码因大雪封校和周五晚上的源码tcp ip monitor 源码闲暇时光，我终于有了时间来解决这个问题。源码

       我撰写这篇博客的源码目的，是源码希望从数据结构的角度，详细解析 OpenFOAM 如何对一个案例进行运算。源码我将从一个最简单的源码例子出发（1*2*3的网格），通过数字来演示 OpenFOAM 的源码运算过程。本文将分为两部分，源码Cleverqq源码首先分析不带湍流模型的源码 `icoFoam` 例子，之后有空时再探讨带湍流模型的源码 `pisoFoam`。

       在深入代码之前，了解其数学表达式至关重要。关于 `icoFoam` 的数学模型，我主要参考了李东岳博士的论文。`icoFoam` 是一个基于 NS 方程（无湍流项）的简单例子，适用于分析流体动力学基本原理。

       NS 方程描述了动量方程，我们假设忽略压力梯度项，以便进行动量预测。OpenFOAM 使用有限体积法（FVM）对每个项进行体积积分离散化，BW源码最终形成如下方程：

       Vp * dU/dt = Ff - SS * ν * ∇U

       其中，Vp 表示网格单元体积，Ff 为通量，SS 是网格单元面矢量，ν 是动力粘度，U 是速度向量，下标 n 表示当前时间步（已知），r 表示预测时间步（待求）。N 和 P 分别代表相邻网格单元和当前网格单元。

       接着，我们将忽略了的压力梯度项加入方程中，得到：

       Vp * dU/dt = Ff - SS * (p/rho) - SS * ν * ∇U

       这里的libblkid源码 p/rho 表示单位压力，是 OpenFOAM 中定义的压力。

       这实际上是一个线性方程，方程在某一时间 n 上，除了 U^r 未知，其他变量都是已知的，因此可以看作是求解线性方程组（ax=b）的问题。

       接下来是压力泊松方程，我们需要根据 U^r 预测出下一时刻的压力 p^r，以完成循环并满足物理约束。

       循环内，我们通过不断修正 U^r 和 p^r，加入动量方程和连续性方程的libnuma源码物理约束，最终得到满足这些约束的 U^r 和 p^r，再带入下一个时间步进行计算。

       在 `pisoFoam` 中，循环的过程是如何更新速度和压力的呢？

       将动量方程代入连续性方程，可以得到浓缩的 2 合 1 方程式，即压力泊松方程。其中，HbyA 是基于省略了压力项的动量方程计算出的速度预测值，用于更新变量。

       通过 piso 循环，最终目的是在每个循环中得到满足 2 合 1 方程式的预测结果，这样得到的预测结果在物理上是合理的。

       整个流程包括：

       1. 获取第一个压力预测值，通过将动量方程带入压力泊松方程而求得。

       2. 根据第一个压力预测值修正速度。

       3. 循环重复直到收敛。

       在理解了数学模型之后，接下来是观察代码，看看 OpenFOAM 是如何实现这些概念的。在图中，我总结了对 `icoFoam` 代码的理解，提供了一个清晰的视图。跑一个 1*2*3 网格的 case，如自带的 `cavity` case。

       首先调整 `blockMeshDict` 为 1*2*3 的网格。

       运行代码，观察网格形状。

       从 0 到第一个时间步，再到从第一个到第二个时间步。

       1. 从 0 到第一个时间步。

       这段代码的主要目的是将动量方程放入求解器中，通过 `solve` 命令计算出预测速度 U^r。

       疑问在于，`solve(UEqn == -fvc::grad(p));` 之后速度变量 U 直接更新，这在 C++ 中是合法的吗？U 不是 UEqn 的子变量，这行代码没有赋值给 U，怎么就能更新 U？

       经过验证，在 `solve` 语句前后输出 U 的结果确实不同，说明 `solve` 确实修改了外部变量 U。

       这种操作可能涉及指针（pointer）或引用（reference）的使用，以提高代码效率。在 C 语言中，这样的操作是常见的，但在不了解内存管理的情况下可能会引发问题。

       通过复制一个引用实例，观察了 `icoFoam.C` 中的传递方式。从 `UEqn` 到 `solve` 函数的参数传递，可能是通过某种间接引用实现的。

       在 `fvMatrix` 中，包含 A 和部分 b（如 `fvm::ddt(U)`）的结构。`fvMatrix` 分为两部分，一部分用于 A，另一部分用于 b。

       观察到 `fvMatrix` 中确实包含引用变量 psi_，这是 U 的引用位置。`solve` 函数通过创建引用 psi 来更新 U 的值。

       总结了 `fvMatrix` 的结构，确认了 `solve` 函数在更新 U 值时使用了引用或指针。

       在了解了如何在矩阵 A 和 b 中进行操作后，我们接下来需要验证矩阵是如何变化的。

       修改源代码，输出 `solve` 公式中的变量，进行比较。

       对于输出的矩阵和源代码进行了详细的对比分析，确认了矩阵 A 和 b 的结构与预期一致。

       最后，分析了整个 piso 循环过程，从预测速度到修正压力，再到最终得到满足物理约束的 U^r 和 p^r，整个流程在代码中得到了清晰的体现。