第 19 章：OpenFOAM® 湍流应用算例（An OpenFOAM® Turbulent Flow Application）

19.1 引言（Introduction）

汽车车身的设计是一件要求很高的工作，需要在品牌风格与消费者审美之间找到妥协，同时还要满足燃油经济性所要求的高效气动特性。在这一情境下使用 CFD 工具是至关重要的：它能帮助研究气动力、粘性效应与湍流边界层的相互作用、流场对车身形状变化的敏感性，以及在各种工况下车辆的阻力系数。本章使用 simpleFoamTurbulent 求解器以及前几章开发并实现的边界条件，分析一个被广泛研究的测试算例——Ahmed 钝体（Ahmed bluff body）。

19.2 Ahmed 钝体（The Ahmed Bluff Body）

在汽车应用中，对阻力贡献最大的因素是车体后方形成的尾流，而对这一尾流的预测是 CFD 中一项困难任务。其原因在于：湍流流动中的流动分离在数值上仍是一个挑战，但分离区的大小极大地影响作用在车体上的预测阻力。因此，对诱导涡流和分离过程的精确模拟是正确预测气动效率的关键。

当代汽车设计包含许多复杂的几何特征，使建模和实验研究都具有挑战性。因此，本章选择 Ahmed 钝体（参考文献 [1, 2]）作为研究对象，它在 Fig. 19.1 中示意性地给出。Fig. 19.1a 给出侧视、正视与俯视图，可从中推断出尺寸；Fig. 19.1b 给出该钝体的三维可视化。尽管几何形状简单，Ahmed 钝体能够发展出三维分离流区域，可研究各种流动现象并与实验数据对比。

Ahmed 钝体是一个广为人知的构型，被广泛用作基准。其后端的斜面几何在侧缘产生一对反向旋转的涡，其强度主要由基底斜面角度决定。该钝体有两种构型，区别在于斜面角度值（25° 与 35°）。本章的模拟考虑 25° 斜面角的构型。

19.3 域与网格（Domain and Mesh）

计算域如 Fig. 19.2 所示，沿 Ahmed 钝体形状的中间截面施加对称条件。利用对称性来缩减计算域，并减轻由于钝体后端涡脱落而预期出现的瞬态行为；这同样有助于增强数值稳定性。

入流与出流边界条件被设置在远离钝体的位置，以最小化与主流区域之间不必要的相互作用，特别是出口与钝体后端流动之间的相互作用。

网格由 snappyHexMesh 生成，它是 OpenFOAM® [3] 软件包自带的一个工具。snappyHexMesh 生成三维网格，包含六面体（hex）与分裂六面体（split-hex）单元，这些单元由立体光刻（STL）格式的三角化表面几何自动生成。网格通过对初始网格逐步细化并将所得 split-hex 网格变形到表面上，从而逐渐贴合表面。可选的阶段会将所得网格回退并插入 cell 层。网格细化级别的设定非常灵活，表面处理对预设的最终网格质量也很稳健。它在并行模式下运行，并在每次迭代中加入负载均衡步骤 [4]。

本章的目的不在于对 snappyHexMesh 应用给出详细介绍，但案例目录 "Ahmed_body/system/snappyHexMeshDict" 中提供了一份带注释的设置文件。Listing 19.1 给出了该文件的一段摘录。摘录中 castellatedMesh true; snap true; addLayers true; 表明三个阶段（castellatedMesh 阶段、snap 阶段、addLayers 阶段）都被启用；几何 geometry 块中定义了一个 triSurfaceMesh 类型的 ahmed_body.stl，命名 ahmed_body；castellatedMeshControls 中 maxLocalCells 300000; 表示在细化前先按权重进行负载均衡的阈值。

网格通过在案例目录下按所示顺序执行以下命令生成：blockMesh—由 system/blockMeshDict 定义 snappyHexMesh 工作所在的总体计算域，该域应在其内部包含 STL 几何；snappyHexMesh—snappyHexMesh 通过三个阶段生成网格，每个阶段的结果分别写入文件夹（1/、2/ 和 3/）；mv 3/polyMesh constant/—将 snappyHexMesh 生成的最终网格移动到 constant 目录中供求解使用；createPatch—overwrite—删除所有空的 patch。

贴体网格由大约 800 000 个多面体单元组成，主要分布在钝体周围两个区域。这使得在保持计算域单元数较少的同时能够很好地分辨尾流区域。钝体附近生成的网格细节如 Fig. 19.3 所示。

在感兴趣区域进行网格单元聚类是常见的做法，它能在恰当分辨流动主要特征的同时降低计算成本。生成的网格完全非结构化，单元为多面体，最多有 15 个面。网格质量最好用 checkMesh 来评估，其输出如 Listing 19.2 所示。Listing 19.2 给出统计信息：points: 806136、faces: 2331063、internal faces: 2292117、cells: 762865、faces per cell: 6.060285896、boundary patches: 10。单元类型按面数分类：hexahedra 744610，prisms 1057，wedges 0，pyramids 0，tet wedges 1，tetrahedra 0，polyhedra 17197；多面体中按面数再分类（4 面 82、5 面 59、6 面 4523、7 面 1563、8 面 410、9 面 6961、10 面 5、11 面 3、12 面 2782、14 面 2、15 面 807）。质量检查输出给出：bounding box (-5 -0.00189057368 0) (10 5 5)，boundary openness (-2.59e-17 -1.40e-15 -4.46e-14) OK，max cell openness 3.96e-16 OK，max aspect ratio 13.98 OK，face area magnitudes OK，cell volumes OK（min 8.48e-09，max 0.0587，total 374.94），non-orthogonality max 62.13、average 5.02 OK，face pyramids OK，max skewness 2.52 OK，coupled point location match (average 0) OK，最终输出 Mesh OK. End。

物理边界条件在 "Ahmed_body/constant/polyMesh/" 目录下的 boundary 文件中设定，如 Listing 19.3 所示。Listing 19.3 给出 10 个 patch 列表：inlet（patch, nFaces 298, startFace 2292117）、outlet（patch, nFaces 298, startFace 2292415）、ffmaxy（patch, nFaces 640, startFace 2292713）、floor（wall, nFaces 13779, startFace 2293353）、top（patch, nFaces 520, startFace 2307132）、ahmed_body_body（wall, nFaces 5750, startFace 2307652）、ahmed_body_body_front_h（wall, nFaces 864, startFace 2313402）、ahmed_body_body_front_v（wall, nFaces 1270, startFace 2314266）、ahmed_body_stilts（wall, nFaces 256, startFace 2315536）、cp（wall, nFaces 15271, startFace 2315792）。

patch inlet 与 outlet 定义为 patch 类型，而 cp 设置为 wall 而非 symmetry plane。为降低流动求解的复杂度，cp patch 将采用 slip 边界条件。代表车体的 "ahmed_body*" patch 定义为 "wall"，其余 patch 定义外部计算域。

19.3.1 初始与边界条件（Initial and Boundary Conditions）

0 目录必须包含以下四个基本文件以进行湍流不可压缩流动模拟：速度 U、压力 p，以及对于 kEpsilon 湍流模型，分别为湍流动能 k 与耗散率 ε。文件 U 定义了速度的边界条件（Listing 19.4）。

U 文件的关键内容如下：dimensions [0 1 -1 0 0 0 0];，internalField uniform (40 0 0);；boundaryField 中 inlet 类型为 surfaceNormalFixedValue，refValue uniform -40（注意 refValue 取负以对应 -40 m/s 法向入流）；outlet 类型为 inletOutlet，inletValue uniform (0 0 0)，value uniform (0 0 0)；ffmaxy 与 top 为 slip；ahmed_body_body、ahmed_body_body_front_h、ahmed_body_body_front_v、ahmed_body_stilts、floor 全部为 noSlipWall（即第 18 章中开发的新边界条件），value uniform (0 0 0)；cp 为 slip。

新的 noSlipWall 边界条件被用来在壁面（包括车身与地面）施加无滑移条件。入口处流体速度设为 40 m/s。沿 cp 平面，使用 slip 壁面边界条件——将粘性力设为零，以模拟一个对称面（patch cp 上由于网格划分造成的轻微几何误差使得无法直接使用 symmetryPlane 边界条件）。zeroGradient 压力边界条件使用零阶外推来计算边界处压力。在出口处，应用 Dirichlet 边界条件以设定参考压力（Listing 19.5）。Listing 19.5 给出压力 boundaryField：inlet 为 zeroGradient；outlet 为 fixedValue，value uniform 0；ffmaxy 为 zeroGradient。

对于 kEpsilon 模型，在入流边界上设定湍流强度与积分长度尺度，而在所有壁面上自动采用标准壁面函数（第 17 章）。k 文件如 Listing 19.6 所示，ε 文件如 Listing 19.7 所示。k 文件中 inlet 类型为 turbulentIntensityKineticEnergyInlet，intensity 0.01，value uniform 0.002；outlet 为 zeroGradient。epsilon 文件中 inlet 类型为 turbulentMixingLengthFrequencyInlet，mixingLength 0.01，phi phi，k k，value uniform 0.002；outlet 为 zeroGradient。

19.3.2 系统文件（Systems Files）

controlDict 文件（Listing 19.8）配置为执行 500 次 SIMPLE 迭代。此外，为了使用针对速度的新 noSlipWall 边界条件，相应的库必须与入流湍流边界条件所需的库一起被包含。这些通过 libs 声明设定。还包含一个特定的运行时后处理配置，以迭代地监测升力与阻力系数的变化，这通过 functions 声明设定。对于外部气动应用而言，比起残差水平，监测钝体上的载荷或力更便于检查收敛性。

controlDict 文件的关键内容：startFrom startTime;，startTime 0;，stopAt endTime;，endTime 500;；libs ("libNoSlip.so" "libincompressibleRASModels.so" "libincompressibleRASModels.so");；functions 中定义 forceCoeffs1，其类型为 forceCoeffs，functionObjectLibs ("libforces.so")，outputControl timeStep，outputInterval 1，log yes，patches ("ahmed_body*")，pName p，UName U，rhoName rhoInf（表示不可压缩），rhoInf 1（对不可压缩冗余），liftDir (0 0 1)，dragDir (1 0 0)，CofR (0.72 0 0)（地面上车轴中点），pitchAxis (0 1 0)，magUInf 40，lRef 1.45（轴距长度），Aref 2.618（估算面积）。

线性求解器规格与松弛因子在 fvSolution 文件中配置，如 Listing 19.9 所示。这里为所有方程选择带不完全分解预条件子的预处理共轭梯度法。对于压力修正方程（pp），由于该方程的拉普拉斯性质，指定对称矩阵求解器。值得指出的是，使用压力修正格式时不需要非正交修正迭代。

fvSolution 文件中 solvers 块的关键内容：pp 用 PCG，preconditioner DIC，tolerance 1e-06，relTol 0.01；"(U|k|epsilon)" 用 PBiCG，preconditioner DILU，tolerance 1e-15，relTol 0.1。SIMPLE { nNonOrthogonalCorrectors 0; }。relaxationFactors 中 pp 0.3，U 0.7，k 0.7，epsilon 0.7。

空间离散设置在 fvSchemes 文件中配置，如 Listing 19.10 所示。对于本模拟，对流项使用一阶迎风格离散，梯度重构使用高斯方法，变量面插值使用线性剖面。

fvSchemes 文件的关键内容：ddtSchemes { default steadyState; }；gradSchemes { default Gauss linear; grad(p) Gauss linear; grad(U) Gauss linear; }；divSchemes 中 default Gauss upwind，div(phi,U) Gauss upwind，div(phi,k) Gauss upwind，div(phi,epsilon) Gauss upwind，div(R) Gauss linear，div((nuEff*dev(T(grad(U))))) Gauss linear；laplacianSchemes 中 default Gauss linear corrected，对 laplacian(nuEff,U)、laplacian((1|A(U)),p)、laplacian(DkEff,k)、laplacian(DepsilonEff,epsilon)、laplacian(DREff,R)、laplacian(DnuTildaEff,nuTilda) 一律 Gauss linear corrected；interpolationSchemes { default linear; interpolate(U) linear; }；snGradSchemes { default corrected; }；fluxRequired { default no; pp ; }。

19.3.3 运行求解器（Running the Solver）

在运行开始时，求解器提供方程残差与力系数信息（Listing 19.11）。基于 controlDict 的设置，力系数也会被写入文件 Ahmed_body/postProcessing/forceCoeffs1/0/forceCoeffs.dat，可用于检查收敛历史。

Listing 19.11 显示第一次迭代（Time = 1）的求解器详细输出：求解 Ux（DILUPBiCG: Initial residual = 1, Final residual = 0.02308834211, No Iterations 2）、求解 Uy（Final residual 0.07259752325, No Iterations 1）、求解 Uz（Final residual 0.05706016776, No Iterations 1）、求解 pp（DICPCG: Final residual 0.009688743115, No Iterations 140）、求解 epsilon（Final residual 2.962105429e-05, No Iterations 1）、求解 k（Final residual 0.000587089458, No Iterations 1）；ExecutionTime = 9.73 s ClockTime = 9 s；forceCoeffs 输出 Cm = 5.971060664, Cd = 18.59033988, Cl = 3.888152128, Cl(f) = 7.915136728, Cl(r) = -4.0269846。

值得指出的是，对于压力修正方程，任何迭代上的残差始终为 1，这是由于 OpenFOAM® 归一化残差的方式（第 10 章中的计算指引）。Fig. 19.4a 给出动量方程残差的收敛历史曲线，Fig. 19.4b 给出升力与阻力系数随模拟进行的变化。Result 明确显示，在 300 次迭代时解已经完全收敛，升力与阻力系数值的变化已变得可忽略。Fig. 19.5 给出钝体周围的静压等值线。由于动压分量的恢复，车身前部压力值显著。等值线还显示，由于流体在曲面上的加速，车身前方位置下游存在低压区。

Fig. 19.6 中钝体后侧出现回流区域（recirculation region），并伴随低速区。Fig. 19.7 中钝体后部形成的尾流在矢量图中清晰可见。该涡流区域的形成是压力损失的原因，是对钝体阻力贡献的主要来源。

19.4 小结（Conclusion）

simpleFoamTurbulent 求解器与前几章开发的边界条件被用于求解 Ahmed 钝体周围的湍流。尽管采用了低阶对流格式与相对粗糙的网格，流动的主要特征仍得到了展示。

本章个人批注

本章是 Moukalled 全书的最后一章，承担"应用算例"的收尾角色。前 18 章建立了从一般输运方程的有限体离散、压力基算法（SIMPLE/SIMPLEC）、湍流模型（标准 k-ε、Realizable k-ε、k-ω SST 等）到 OpenFOAM®/uFVM 边界条件实现的完整知识链条；本章则把这一切打包成"一个真实算例"演示出来。值得注意的是作者选择了 Ahmed 钝体作为基准——这是一个被工业界和学术界广泛接受的几何构型，其气动特性有大量公开实验数据（参考文献 [1] Ahmed 1984、[2] Lienhart 2000）。这意味着本章不仅是一个教程，更是一个"对照实验"：读者可以用自己的结果与实验数据对比来验证其理解。

让我印象深刻的是本章对对称性 + slip 边界条件的组合处理。Ahmed 钝体在中间截面上几何对称，但作者没有直接使用 symmetryPlane 边界条件，而是选择了 cp 壁面配合 slip 边界条件，理由是 "the slight geometric error due to meshing on patch cp does not allow the direct use of a symmetryPlane boundary condition"。这是一个工程上很现实的考虑：snappyHexMesh 生成的网格在 cp 平面附近不可能严格几何对称，强制使用 symmetryPlane 会引入非物理约束；而 slip（零剪切应力）则"放过了"法向速度的不匹配，代价是无法强制两侧速度完全镜像——这个权衡是工程数值方法的一个典型微调。另一个细节是 surfaceNormalFixedValue 在入口处的设置：作者用 refValue uniform -40 而 internalField 是 (40 0 0)，负号源于 surfaceNormalFixedValue 中"参考方向是 patch 的外法向"——这是 OpenFOAM® 边界条件语义的微妙之处。

我还注意到作者在 Listing 19.8 中写"0 0 1"作为 liftDir，而车体的实际 lift 方向应该是 z 方向（垂直于地面）——这里 liftDir (0 0 1) 与 dragDir (1 0 0)、pitchAxis (0 1 0) 共同构成右手系。在 magUInf 40、lRef 1.45、Aref 2.618 这一组参数下，归一化后得到的 Cd = 18.59 在第一次迭代时显然还未收敛到合理量级——这是初始残差为 1 的瞬态显示值，需要在 ~300 步后才能稳定到工业级 Cd ~0.25 范围。书中没有给出最终收敛的 Cd 数值，这多少让我有些意犹未尽；不过从收敛曲线 Fig. 19.4b 来看，作者已经暗示了"300 步后解完全收敛"这一事实——这与书中"500 步"的总配置形成了一致的安全余量。

从全书视角看，ch19 的意义在于它把前 18 章抽象的 FVM 理论与 OpenFOAM® 的具体实现对接到一个"用户视角的真实算例"上。读者只要按书中 Listing 19.1–19.10 抄写 case 文件、运行 blockMesh → snappyHexMesh → simpleFoamTurbulent 就能复现整个工作流。这种"理论 → 工具链 → 算例"的三段式收尾，是教科书把读者从"理解"带到"动手"的常用手法。

与上下章的衔接（一段话）

第 19 章是 Moukalled 全书（Ch. 1–19）的最后一章，承接第 18 章末尾"the next chapter is devoted to detailing the steps needed to solve a turbulent flow problem in OpenFOAM®"的明确预告。第 18 章开发了 noSlipWall 边界条件并展示了 uFVM 的对应实现，本章则在"非正交网格 + 标准 k-ε + SIMPLE 算法 + 多种边界条件类型"的最复杂组合下，把这套工具链落到一个工业级 CFD 基准算例（Ahmed 25° 钝体）上。从全书曲线看：第 1–4 章铺垫 FVM 基础与一般输运方程的数学形式，第 5–12 章建立离散与求解框架，第 13–14 章处理源项与松弛，第 15 章给出不可压缩流的压力基算法（SIMPLE/SIMPLEC），第 16 章拓展到可压缩流，第 17 章引入湍流模型，第 18 章完成边界条件的代码实现，第 19 章则把这一切"打包展示"——构成了 Moukalled 团队"数学 → 数值 → 物理建模 → 软件实现 → 工程算例"的完整闭环。本章引用的 4 篇参考文献（Ahmed 1984、Lienhart 2000、OpenFOAM® 2.3.x、OpenFOAM® User Guide）也都是 OpenFOAM® 社区最常引用的资源——这一选择本身就暗示了本章在"教学算例"与"工程工具"之间的桥梁作用。