跳转至

第 18 章:OpenFOAM® 与 uFVM 中的边界条件(Boundary Conditions in OpenFOAM® and uFVM)

18.1 OpenFOAM® 中的边界条件(Boundary Conditions in OpenFOAM®)

本章开篇即强调,每种边界条件都对应一个物理意义,由数学方程刻画(Dirichlet、Neumann、Robin 等),而在数值方法中必须转化为代数关系。例如入口边界条件同时对速度与压力施加由物理条件给定的数学方程(Dirichlet 与 Neumann),其实现方式取决于所作用的算子(divergence、laplacian、gradient 等)。换言之,Dirichlet/Neumann 等数学分类是边界条件"做什么",而 divergence/laplacian 算子分类是边界条件"作用于哪个方程项"——一个具体边界条件类必须同时回答这两个问题,"Dirichlet 之于 divergence"与"Dirichlet 之于 laplacian"在代数实现上是不对称的,因此一个 BC 类通常要为多个算子提供各自的线性化形式。在 OpenFOAM® 中,几乎所有边界条件的定义都存放在目录 src/finiteVolume/fields/fvPatchFields(Listing 18.1)下,其中主要的边界条件类型分别放在四个子目录(Listing 18.2)中。

第一个子目录 fvPatchField 包含边界条件的一般类定义(即基类)。这个类定义了将被真正边界类继承与使用的主要函数与数据结构,因此所有具体边界条件类都从这里派生——这也是 OpenFOAM® 实现模板化边界条件(同一套接口可用于 scalar / vector / tensor 任意类型)的关键:基类用 template<class Type> 参数化所有 field 类型,因此下面所有具体 BC 都可以基于同一套接口针对标量、矢量、张量分别实例化。第二子目录 basic 包含基本数学定义的边界条件:固定值型(fixedValue,对应 Dirichlet)、零梯度与固定梯度型(zeroGradient / fixedGradient,对应 Neumann)、混合型(mixed,对应 Robin)。该目录下还有一个额外的条目——耦合边界条件(coupled),实现 patch-to-patch 类型条件,把两个边界 patch 耦合在一起;这是与"几何/数学边界条件"不同的另一类——它不指定具体的 Dirichlet/Neumann 值,而是把两个 patch 通过值映射连接,让一边 patch 上的值由另一边 patch 上的值代数决定。第三个子目录 constraint 包含由 coupled 边界类派生的几何类型边界条件,例如周期性边界条件(Fig. 18.1):在这种情况下每个网格 cell 与对应 patch 上的 cell 相关联,使边界 cell 可被当作内部 cell 处理;这是用"几何上的对应关系"把边界条件完全消解到内部方程中的典型做法——边界 cell 拥有完整的内部面通量处理(同样的邻 cell 索引、相邻面的几何信息),只不过其中一个"邻 cell"实际上位于另一侧 patch 上。最后一个子目录 derived 包含从基本 Dirichlet、Neumann、Robin 边界条件派生的所有边界条件;这些派生类本质上只是基本类型的特化——比如把 fixedValue 模板特化成只取 0 的 zeroFixedValue、把 zeroGradient 特化成对称面条件、把 mixed 特化成 inflow/outflow 等工业常用模板。这种 basic/constraint/derived 三段式目录组织,是 OpenFOAM® 把"基础数学类型 / 几何类型 / 工程特化类型"三者解耦的设计体现——新增一类 BC 时,作者只要决定它是"基础/几何/特化"中哪一类,再放到对应子目录即可,而无需触碰核心库。

18.2 边界条件定制(Boundary Condition Customization)

要编写一个新的边界条件,必须理解五个主要函数的作用:updateCoeffs、valueInternalCoeffs、valueBoundaryCoeffs、gradientInternalCoeffs、gradientBoundaryCoeffs。其中 updateCoeffs 负责显式更新边界面形心处的值,每当 patch 字段值需要迭代更新时调用。它本身只负责"显式计算 patch 面形心值"——不涉及矩阵组装。例如在 totalTemperatureFvPatchScalarField 类中(Listing 18.3),updateCoeffs 依据关系 \(T = T_0 - 0.5(c-1)U^2/cR\) 由给定的总温 \(T_0\) 计算静温 \(T\):先用 patch().lookupPatchField<volVectorField, vector>(UName_) 取 patch 上 U 场、用 lookupPatchField<volScalarField, scalar>(psiName_)\(\psi\) 场(用于 \(\gamma = 1/(R\psi)\)),计算 gM1ByG = (\gamma - 1)/\gamma,再用 operator==(T0_/(1.0 + 0.5*psip*gM1ByG*magSqr(Up))) 把结果赋给静温场;函数入口的 if (updated()) return; 保证一次迭代内不重复计算;最后调用 fixedValueFvPatchScalarField::updateCoeffs() 完成基类的状态标记。另一个例子是对整个 patch 设置平均值(Listing 18.4):其思路是把指定均值施加到 patch 所有面上,先取 patch 内场值 Field<Type> newValues(this->patchInternalField()),按面积加权计算当前均值 gSum(magSf * newValues)/gSum(magSf),再用 newValues += (meanValue_ - meanValuePsi) 把每个面加上一个均匀修正量,最后 this->operator==(newValues) 把修正后的值赋给 patch 场;该函数为模板形式,因此可与 scalar、vector、tensor 等多种类型联用。这里的关键观察是:updateCoeffs 只"显式"地写入 patch 值,矩阵贡献完全交给后面的 4 个 *Coeffs 函数——这种"显式 patch 值 + 线性化矩阵系数"的分工模式贯穿 OpenFOAM® 整个边界条件体系。

updateCoeffs() 把面形心处的值显式设置为 Dirichlet 条件的形式,而另外四个函数(valueInternalCoeffs、valueBoundaryCoeffs、gradientInternalCoeffs、gradientBoundaryCoeffs)则用于对边界条件进行线性化以补足 updateCoeffs() 的功能。其中 valueInternalCoeffs 与 valueBoundaryCoeffs 一般用于对 divergence 算子的边界条件进行线性化(因为该算子需要 patch 面上的值);gradientInternalCoeffs 与 gradientBoundaryCoeffs 则用于对 laplace 型算子的边界条件进行线性化(因为该算子需要 patch 面上的梯度)。Table 18.1 给出这一对应关系的小结:divergence 算子用 valueInternalCoeffs / valueBoundaryCoeffs 作为对角 / 源项系数,laplacian 算子用 gradientInternalCoeffs / gradientBoundaryCoeffs。这一划分反映了"散度项关心面处场值、扩散项关心面处梯度"这一物理事实——把系数按算子语义分组,避免边界条件类同时承担两套接口。

为澄清上述原则,本节考察两种边界条件的实现:Neumann(zeroFlux)与 Dirichlet(specifiedValue)。对于正交网格上的 Neumann(zeroFlux)边界条件,由于法向梯度为零,patch 上面的值等于边界 cell 的值(即 \(\phi_b = \phi_C\))。对于 divergence 项(需要 patch 面值),边界 patch 的给定值用 valueInternalCoeffs 与 valueBoundaryCoeffs 写出(分别代表边界 cell 值的线性化部分与其非线性化部分)。例如对零梯度条件(方程 18.1):

\[ \phi_b = \mathrm{FluxC_b}\,\phi_C + \mathrm{FluxV_b} = \mathrm{valueInternalCoeffs}\,\phi_C + \mathrm{valueBoundaryCoeffs} = 1\cdot\phi_C + 0 \]

Listing 18.5 给出了 zeroFlux 条件对 divergence 算子的语法实现:valueInternalCoeffs 返回 tmp<Field<Type>>(new Field<Type>(this->size(), pTraits<Type>::one))(即尺寸等于面数的全 1 向量),valueBoundaryCoeffs 返回 tmp<Field<Type>>(new Field<Type>(this->size(), pTraits<Type>::zero))(即全 0 向量)。注意这里的输入参数是 const tmp<scalarField>&,在 zeroFlux 实现中并未使用——它是 framework 设计的统一接口预留。对于 laplacian 算子,所需的是边界梯度,此时通过梯度线性化(方程 18.2)

\[ \nabla\phi_b = \mathrm{gradientInternalCoeffs}\,\phi_C + \mathrm{gradientBoundaryCoeffs} = 0\cdot\phi_C + 0 \]

设为零(Listing 18.6:gradientInternalCoeffs 与 gradientBoundaryCoeffs 均返回 pTraits::zero)。Table 18.2 汇总了 zeroFlux 条件下这两个系数在 divergence 与 laplacian 下的取值——divergence 时 value 对角系数为 1、源项为 0;laplacian 时 gradient 对角系数与源项均为 0。其物理意义是:因为指定了零通量,边界处的值被设为等于边界 cell 的值,梯度也强制为零。

对于 Dirichlet 边界条件,对系数矩阵的贡献只是方程右端的源项,边界条件不改变对角线。在 OpenFOAM® 中 valueInternalCoeffs 与 valueBoundaryCoeffs 按方程 18.3 定义为:

\[ \phi_b = \mathrm{FluxC_b}\,\phi_C + \mathrm{FluxV_b} = \mathrm{valueInternalCoeffs}\,\phi_C + \mathrm{valueBoundaryCoeffs} = 0\cdot\phi_C + \phi_\mathrm{specified} \]

其物理含义是:边界 patch 上的值是已知的(\(\phi_\mathrm{specified}\)),与内部 cell 值无代数关系,因此对角系数取 0、源项取边界指定值——这正是 Dirichlet 条件"指定边界值而非边界法向梯度"在代数层面的体现。对于 laplacian 算子,边界处的梯度基于 Dirichlet 值,再一次通过梯度线性化写出:

\[ \nabla\phi_b = \mathrm{gradientInternalCoeffs}\,\phi_C + \mathrm{gradientBoundaryCoeffs} = \frac{\phi_C - \phi_b}{\delta} = \left(-\frac{1}{\delta}\right)\phi_C + \frac{\phi_b}{\delta} \]

其中 \(\delta\)patch().deltaCoeffs() 给出的 cell-to-face 距离倒数(\(\delta = 1/|\mathbf{d}|\))。Table 18.3 汇总了 Dirichlet 条件下各系数取值(divergence 时对角系数为 0、源项为边界值;laplacian 时对角系数为 Delta、源项含边界值与 delta)。对应的模板代码见 Listing 18.7:valueInternalCoeffs 返回 pTraits<Type>::zero(对角归零),valueBoundaryCoeffs 返回 *this(即 boundary patch 自身的值——这是 OpenFOAM® 中让"边界值"直接作为源项的常见 trick);gradientInternalCoeffs 返回 -pTraits<Type>::one * this->patch().deltaCoeffs()gradientBoundaryCoeffs 返回 this->patch().deltaCoeffs() * (*this)。对照 Neumann 实现可见:对 divergence 算子,Neumann 把值传给对角系数(1),Dirichlet 把值传给源项(边界值);对 laplacian 算子,Neumann 全归零,Dirichlet 通过 \(\delta\) 把梯度拆为对角 \(-\delta\) + 源项 \(\phi_b \delta\)——这种"算子 × 边界类型"的二维分类决定了每个 BC 类的具体 *Coeffs 形式。

18.3 新边界条件的开发:无滑移壁面条件(Development of a New BC: No Slip Wall Condition)

本节的任务是在 OpenFOAM® 有限体离散框架内为无滑移条件(no slip)定义一种新的边界条件类型(Fig. 18.2)。无滑移条件是流动问题求解中的一个基本边界条件。根据 Newton 粘性定律,粘性流体流过壁面时所受的剪切应力与平行于壁面的速度法向梯度成正比,写为方程 18.4:

\[ \tau_\mathrm{wall} = -\mu\,\frac{\partial \mathbf{v}_\parallel}{\partial (d_\perp)_C} = -\frac{\mu}{(d_\perp)_C}\begin{bmatrix}(1-n_x^2) & -n_y n_x \\ -n_y n_x & (1-n_y^2)\end{bmatrix}\begin{bmatrix}u_C \\ v_C\end{bmatrix} \]

其中下标 \(C\) 指边界 cell 形心处的值,\(\mathbf{v}_\parallel\) 是切向速度分量,\(n = (n_x, n_y)\) 是壁面外法向单位向量。该方程表明无滑移条件是各向异性的——它只依赖于平行于壁面的速度分量与到壁面的法向距离。当壁面与 x 轴平行时(\(n = (0, 1)\)),只有速度的 x 分量影响剪切应力值;事实上由方程 18.4 可得方程 18.5:

\[ \tau_\mathrm{wall} = -\mu\,\frac{\partial \mathbf{v}_\parallel}{\partial (d_\perp)_C} = -\frac{\mu}{(d_\perp)_C}\begin{bmatrix}1 & 0 \\ 0 & 0\end{bmatrix}\begin{bmatrix}u_C \\ v_C\end{bmatrix} \]

显然只有速度 x 分量的对角系数需要注入矩阵——这个"非对角项消失"的退化情况是无滑移条件各向异性矩阵在 x 平行壁面下的特例。OpenFOAM® 中对此条件的一个常见简化是把它当成 Dirichlet 边界条件处理,即 fixedValue (0 0 0)

然而把无滑移条件用 Dirichlet 形式实现会引入一个重要误差:这种做法把剪切应力分布当作各向同性、并使其值受(内部 cell 中)法向速度分量影响。事实上若壁面与 x 方向速度分量平行而采用 Dirichlet 边界条件,则剪切应力被估计为方程 18.6:

\[ \tau_\mathrm{wall} = -\mu\,\frac{\partial \mathbf{v}_\parallel}{\partial (d_\perp)_C} \approx -\mu\,\frac{\partial \mathbf{v}}{\partial (d_\perp)_C} = -\frac{\mu}{(d_\perp)_C}\begin{bmatrix}1 & 0 \\ 0 & 1\end{bmatrix}\begin{bmatrix}u_C \\ v_C\end{bmatrix} \]

与方程 18.4/18.5 比较可知,这种近似会使壁面剪切应力不当地依赖于速度的 y 分量——而对平行 x 轴的壁面,y 分量是法向分量,物理上不应进入切向剪切应力的计算。这里 \(\partial \mathbf{v}/\partial(d_\perp)\)

\[ \partial \mathbf{v}_\parallel/\partial(d_\perp) \]

之差是"切向速度梯度"与"全速度梯度"的差——把全速度梯度代入,意味着法向速度分量 \(v_C\) 也被错误地当成切向贡献处理。以下给出 OpenFOAM® 中无滑移壁面条件的正确实现方法,作者建议在壁面动量离散中始终用这种实现代替 Dirichlet 简化假设。

如前所述,定义新边界条件需要重新定义相应的函数。一个好的起点是拷贝现有边界类,然后修改其算法实现新逻辑。为此可使用 fixedValue 类并针对 vector 类型进行改造。定义新虚类 noSlipWall 的 .H 文件概要见 Listing 18.8:class noSlipWallFvPatchVectorField : public fixedValueFvPatchField<vector>,并以 TypeName("noSlipWall") 注册到运行时选择机制——这样在 case 文件中写 type noSlipWall; 时 OpenFOAM® 就能在工厂中查到该类;类体内只声明两个新函数 gradientInternalCoeffs()gradientBoundaryCoeffs()。值得注意的只是重新定义了 gradientCoeffs 函数,因为 fixedValue 类的 valueCoeffs 与 updateCoeffs 函数已经正确实现——也就是说,updateCoeffs() 仍然只把边界值设为 0,但这次我们要让 gradient 系数体现各向异性。函数声明就绪后,.C 文件的修改从方程 18.1 出发,但写成三维形式(方程 18.7):

\[ \tau_\mathrm{wall} = -\mu\,\frac{\partial \mathbf{v}_\parallel}{\partial (d_\perp)_C} = -\frac{\mu}{(d_\perp)_C}\begin{bmatrix}(1-n_x^2) & -n_y n_x & -n_z n_x \\ -n_y n_x & (1-n_y^2) & -n_z n_y \\ -n_z n_x & -n_z n_y & (1-n_z^2)\end{bmatrix}\begin{bmatrix}u_C \\ v_C \\ w_C\end{bmatrix} \]

方程 18.7 表明矩阵中所有元素都有隐式贡献(它们依赖于内部 cell 值),因此必须用 gradientInternalCoeffs 实现;另一项是壁面速度的切向分量,只依赖于面值,需要用 gradientBoundaryCoeffs 描述。更仔细审视可知,由于动量方程采用分离求解,隐式贡献中的混合项(即 \(-n_y n_x\)\(-n_z n_x\) 等非对角元素)宜存到右端作显式处理。据此方程 18.7 重写为方程 18.8:第一项(对角块 \(\mathrm{diag}(1-n_i^2)\))作为 gradientInternalCoeffs 给出,第二项(含 \(-n_i n_j\) 混合项)作为 gradientBoundaryCoeffs 给出,第三项( \(-\mu\,\mathbf{v}_{\parallel,\mathrm{wall}}/(d_\perp)_C\) )是只依赖面值的额外 gradientBoundaryCoeffs。这种"对角隐式 + 非对角显式"的拆分是 SIMPLE 类分离求解器下对角占优矩阵的标准做法——避免一个 cell 的动量方程依赖多个方向的 cell 值(那样在分离求解中会产生延迟修正问题)。

基于以上公式,gradientInternalCoeffs 函数如 Listing 18.9 所列:先取面法向 vectorField normal = this->patch().nf(),再创建尺寸等于面数的零向量 vectorField impCoeff(this->size(), pTraits<vector>::zero),对每个面 i 由 impCoeff[i][k] = 1 - normal[i][k]^2(k=0,1,2)构造对角矩阵,最后返回 -impCoeff * this->patch().deltaCoeffs()。注意这里返回的是 vectorField 而非标量 field——因为 noSlipWall 是 vector 类型 patch 字段,gradient 系数本身是逐分量的对角张量。gradientBoundaryCoeffs 的计算见 Listing 18.10:先取面法向 normal,计算切向面场 boundTanField = (*this) - ((*this) & normal) * normal(从面场中减去其法向分量),再取 patchInternalField() 存为 intField;最后对每个面 i 按混合项 \(-n_i n_j\) 写出三个分量(其中显式项含 boundTanField 与 intField 的乘积),返回 this->patch().deltaCoeffs() * expCoeff。代码中三个 forAll 循环逐分量写出 expCoeff[i][k] 的形式与方程 18.8 严格对应:第一项是壁面值切向分量 boundTanField[i][k],第二、三项是 intField 与 normal 分量的混合。新边界条件可在任意壁面上使用,只需把 patch 类型定义为 noSlipWall(Listing 18.11:type noSlipWall, value uniform (0 0 0))。

18.4 uFVM 中的无滑移边界条件(The No-Slip Boundary Condition in uFVM)

在 uFVM 中新增边界条件需要为它可能作用的每个 term 都实现一份,不如 OpenFOAM® 模块化,但实现方式相当直接(以无滑移为例即可说明)。需要强调的是,无滑移边界条件是一种"混合"条件:必须同时计算一个 flux(剪切应力)并保证边界速度被设置为指定值(这里为零),因此与 Dirichlet 和 Neumann 都有相似之处。

在 uFVM 中,边界条件信息存储在一个基于 patch 的结构里,由四个数组组成:

  1. thePatchFlux.FLUXC1f:owner cell 的线性化系数
  2. thePatchFlux.FLUXC2f:neighbor cell 的线性化系数
  3. thePatchFlux.FLUXVf:非线性化部分
  4. thePatchFlux.FLUXTf:面处的总 flux

按此分解,边界 flux 写为方程 18.9:

\[ \mathrm{thePatchFlux.FLUXTf} = \mathrm{thePatchFlux.FLUXC1f}\cdot\phi_\mathrm{owner} + \mathrm{thePatchFlux.FLUXC2f}\cdot\phi_\mathrm{boundary} + \mathrm{thePatchFlux.FLUXVf} \]

方程 18.9 展示任意内部面上总 flux 如何按共享该面的 owner 与 neighbor 元素线性化,但对面是 patch 的情况 neighbor 节点实质上是边界节点。对于没有 neighbor 元素定义的边界 face,thePatchFlux.FLUXC2f 始终设为零——这是"边界 face 没有对侧 cell"的直接代数结果;对于内部 face,FLUXC1f 与 FLUXC2f 都非零,从而完整描述了"面通量依赖两侧 cell"的代数关系。值得注意的是,方程 18.9 是 uFVM 中所有 term(扩散项、对流项、源项等)通用的 flux 线性化模式——这意味着 uFVM 用统一的四元组而非 OpenFOAM® 那种"按算子分类的 4 个 *Coeffs 函数"来描述边界条件。对于无滑移条件——其总 flux(即剪切应力)依赖于平行于壁面的速度分量变化——方程 18.10 给出各贡献的具体表达式:

\[ \begin{aligned}\mathrm{theFluxes.FLUXC1f}(\mathrm{iBFaces}) &= \mathrm{area}\cdot\mathrm{TM}\cdot\left(1 - \mathrm{dot}(n_{c0}, n_{c0})\right) \\ \mathrm{theFluxes.FLUXC2f}(\mathrm{iBFaces}) &= 0 \\ \mathrm{theFluxes.FLUXVf}(\mathrm{iBFaces}) &= F_c - \mathrm{area}\cdot\mathrm{TM}\cdot\left(1 - \mathrm{dot}(n_{c0}, n_{c0})\right)\cdot v_c(\mathrm{iOwners}) \\ \mathrm{theFluxes.FLUXTf}(\mathrm{iBFaces}) &= \mathrm{theFluxes.FLUXC1f}(\mathrm{iBFaces})\cdot v_c(\mathrm{iOwners}) \\ &\quad + \mathrm{theFluxes.FLUXC2f}(\mathrm{iBFaces})\cdot v_c(\mathrm{iBElements}) \\ &\quad + \mathrm{theFluxes.FLUXVf}(\mathrm{iBFaces})\end{aligned} \]

其中 \(F_c\) 是按方程 18.11 计算的实际壁面剪切 flux:

\[ F_c = -\mu\,\frac{\mathbf{v}_\parallel}{d_\perp} \]

注意 \(F_c\) 本身正是 OpenFOAM® 方程 18.4 / 18.7 描述的 \(\tau_\mathrm{wall}\),只不过 uFVM 按 term-level 分解后把它单独作为显式修正项而非并入 gradientInternalCoeffs——对照 18.3 节方程 18.8 的三项分解可知,\(F_c\) 对应 OpenFOAM® 那里的第三项 \(-\mu\,\mathbf{v}_{\parallel,\mathrm{wall}}/(d_\perp)_C\) ,而 \(F_c\) 减号之后的部分

\[ -\mathrm{area}\cdot\mathrm{TM}\cdot(1-\mathrm{dot}(n_{c0}, n_{c0}))\cdot v_c(\mathrm{iOwners}) \]

对应 OpenFOAM® 那里 gradientBoundaryCoeffs 中的混合项 \(-n_i n_j v_C\)。两者数学等价,差别只在 C++ 与 MATLAB 的代码组织方式。

完整代码见 Listing 18.12。函数签名为 theFluxes = cfdAssembleStressTermWallNoSlipBC(iPatch, theFluxes, theEquationName, theTerm, iComponent),先取网格、流体 tag 与 boundary 信息(iFaceStart = theBoundary.startFaceiFaceEnd = iFaceStart + numberOfBFaces - 1iBFaces = iFaceStart:iFaceEnd),再由 cfdGetMeshField 取出 term 字段,提取第 iComponent 分量作为 \(v_c\);iComponent 由 1/2/3 决定 e = [1;0;0]/[0;1;0]/[0;0;1](即考察 x / y / z 方向的动量方程);随后取出粘度字段 visc;从网格几何中取 geodiffTfarea、面单位法向 n = [Sf]./[area area area]iOwners;用 cfdGetUCoef 取 TM 系数(即非正交修正等几何因子);随后按 vel_n = (vel·n) nvel_t = vel - vel_n 切/法分解 owner 与 boundary 速度(mag = dot(vel', n')' 取法向分量),并令 vel_wall = velb_t(即把壁面切向速度设为边界 cell 切向速度——在无滑移条件下壁面值 \(\mathbf{v}_{\parallel,\mathrm{wall}} = 0\),但代码仍先按非零切向速度提取,再在 \(F_c\) 项中减去)、nc = n e(外法向与所考察方向单位向量的乘积),

\[ F_c = -\mathrm{area}\cdot\mathrm{TM}\cdot((\mathbf{vel}_t - \mathbf{vel}_\mathrm{wall})\cdot e) \]

;最后按方程 18.10 装填 FLUXC1f、FLUXC2f、FLUXVf、FLUXTf 四个数组。

代码中 iBFaces 与 iBElements 的索引约定值得注意:iElementStart = numberOfElements + iFaceStart - numberOfInteriorFacesiElementEnd = iElementStart + numberOfBFaces - 1——uFVM 把"虚拟的边界元素"放在总 cell 数组(长度 numberOfElements)之后的扩展区域,即把内部 cell 与 boundary element 用同一索引空间编号。这样在通用 flux 装配循环中,owner 与 neighbor 都可以用相同的索引语法 velc(iOwners)velc(iBElements) 访问,无需为边界 face 单独写一版装配代码——这是 uFVM 通过"扩展虚拟元素"实现边界 face 与内部 face 统一处理的小技巧。

另一个几何上的细节是 $1 - \mathrm{dot}(n_{c0}, n_{c0})$ 这个因子:在 2D 下 \(n_{c0} = n_x\)(当 \(e = [1;0]\) 时),于是 \(1 - n_x^2\) 正好是方程 18.4 中的对角元素 \((1-n_x^2)\);在 3D 下 \(\mathrm{dot}(n_{c0}, n_{c0}) = n_i^2\),于是 \(1 - n_i^2\) 对应方程 18.7 中的对角元素 \((1-n_x^2)\)\((1-n_y^2)\)\((1-n_z^2)\)(按 iComponent 决定具体哪个分量)。这种"1 - $n^2$ 因子即对角元素"的设计把方程 18.7 的对角矩阵元素用一行点积即表达,是 uFVM 把 OpenFOAM® 18.3 节"逐分量 for 循环写出"(Listing 18.9/18.10)压缩成"矩阵化点积"的关键简化。从这个对比看 uFVM 在边界条件实现上比 OpenFOAM® 更紧凑(用点积取代循环),但代价是代码"非模板/非类继承"——新增一类 BC 时必须复制一整份函数模板。

18.5 小结(Closure)

源文 18.5 节是一段总结性陈述:明确本章讨论了 OpenFOAM® 中边界条件的实现,并详细列出了添加无滑移边界条件所需的各阶段来说明如何向 OpenFOAM® 引入新边界条件;同时还简要讨论了 uFVM 中的边界条件。下一章专门给出在 OpenFOAM® 中求解湍流问题所需的步骤。

本章个人批注

本章是 Moukalled 团队对前 17 章所有边界条件讨论的一次"集中编译",对应 17 章结尾"the next chapter will discuss the implementation of boundary conditions in OpenFOAM® and uFVM"的承诺。从写作结构看,本章由抽象到具体分为三层:(1) 18.1 节建立目录结构与类继承体系(fvPatchField → basic / constraint / derived),相当于把边界条件按"数学定义 vs 几何约束 vs 派生特化"三种来源分类;(2) 18.2 节给出 5 个核心函数(updateCoeffs / 4 个 *Coeffs)的语义,区分 divergence 与 laplacian 算子下系数取值的差异,并以 Neumann/zeroFlux 与 Dirichlet/specifiedValue 两个最小例子演示;(3) 18.3–18.4 节以"无滑移壁面"为载体,分别给出 OpenFOAM® 与 uFVM 中的完整实现。

让我印象最深的是 18.3 节对 Dirichlet 简化的纠偏。OpenFOAM® 默认把无滑移壁面当 fixedValue(0 0 0) 处理,这意味着方程 18.6 中剪切应力矩阵变成各向同性的 \([1, 0; 0, 1]\),从而错误地把壁面切应力依赖于法向(垂直壁面)速度分量。正确公式(方程 18.4)给出各向异性矩阵 \([1-n_x^2, -n_y n_x; -n_y n_x, 1-n_y^2]\),在壁面平行 x 轴时退化为 \([1, 0; 0, 0]\)——只有切向速度分量真正贡献到剪切应力。这一观察反映出 FVM 框架下"边界条件不只是边界值,还决定矩阵对角与源项"的设计哲学:边界条件直接进入代数方程的系数矩阵,而非仅仅在右端赋一个固定值。值得对应到 ch17 末尾的"壁面函数"讨论:那里通过修改 \(\nu_t\) 的壁面值实现壁面剪切应力(方程 17.77),而 ch18 这里通过自定义 gradientInternalCoeffs 实现同一个壁面剪切应力——两条技术路径殊途同归,但 ch17 用的是"边界值的物理间接修正"、ch18 用的是"边界梯度矩阵的直接定制"。

另一个值得记的细节是 18.4 节 uFVM 与 OpenFOAM® 在设计哲学上的差异。OpenFOAM® 把边界条件实现为可被多个算子共用的"通用 patch 类"(只需实现 5 个 *Coeffs 函数即同时影响 divergence 与 laplacian),而 uFVM 把边界条件按 term 实现——每个 term(散度项、扩散项等)单独一份 BC。这反映了 OpenFOAM® 把"边界 patch 视为一类对象 + 通过模板/继承共享数据结构"的设计取向,与 uFVM"按 term 分别组装系数矩阵"的过程式取向之间的根本差别。本质上 OpenFOAM® 走的是 OOP 路线,uFVM 走的是 MATLAB 风格的过程化路线。方程 18.7 → 18.8 的重写也值得注意:作者把"全矩阵隐式"拆成"对角隐式 + 混合显式右端",原因是动量方程分离求解时把交叉项放右端更便于"延迟修正"或松弛处理——这是 FVM 代码层一个常见但隐蔽的工程权衡。值得一提的是 Listing 18.12 中特意用 vel_wall = velb_t 而不是直接设零,这是为了保留公式在"壁面有切向滑移速度"(如旋转壁面、移动壁面)下的通用性;无滑移只是其特例。

与上下章的衔接(一段话)

第 18 章承接第 17 章末尾"the next chapter will discuss the implementation of boundary conditions in OpenFOAM® and uFVM"的明确预告,是 Moukalled 团队对前 17 章所有边界条件讨论(Dirichlet/Neumann/Robin 数学定义、入口/出口/对称/无滑移壁面/周期性几何条件、湍流壁面函数、SIMPLE/SIMPLEC 算法下压力与速度的边界处理)的"代码层集中编译"。从全书的曲线看:第 5 章给出一般输运方程的边界条件数学形式,第 7–12 章建立离散与求解框架,第 13–14 章处理源项与松弛,第 15–17 章完成压力基算法 + 湍流模型。本章则把所有边界条件从数学到 C++/MATLAB 实现路径串起来,因此是"边界条件"主线的最后一站,也是全书最后一章——下一章将进入湍流求解实例(OpenFOAM® 中的湍流算例),属于"应用/算例"主线而非"基础理论"主线。从更广视角看,ch17 用整章篇幅处理湍流这一物理复杂性,ch18 用整章篇幅处理边界条件这一数值接口复杂性,两者各自代表 FVM 落地过程中两类不可回避的工程问题:物理建模与软件实现。本章的最后一节明示"the next chapter is devoted to detailing the steps needed to solve a turbulent flow problem in OpenFOAM®",构成 Moukalled 团队"理论 → 接口 → 应用"三层结构的最后一跳。