OpenFOAM风工程实战：算例骨架、网格生成与湍流校核

简介：这份 PDF 为风工程 OpenFOAM 基础培训中的 OpenFOAM 简介讲义，面向刚接触开源 CFD、希望理解软件结构与设计理念的高校师生和工程技术人员。内容从 OpenFOAM 的历史沿革讲起，梳理其与帝国理工学院 Gosman 小组、Charlie Hill 博士论文及 Henry Weller、Hrvoje Jasak 工作的关联，并重点说明有限体积离散、非结构化多面体网格、面向对象类库封装、大规模并行计算、自动网格运动与拉格朗日粒子追踪等能力。讲义还通过 solve 方程表达式、张量代数与场代数分层、容器类与字典机制等示例，呈现代码重用、模块化和用户定制的实现思路，并简列航空航天、汽车、能源环境等应用方向。包内共 1 个 PDF 文件，约 162KB，属轻量文献资料；已有 555 人学习，适合用作 OpenFOAM 入门概念梳理和后续查证线索。

1. 风工程里为什么先啃 OpenFOAM：一份讲义的第一页决定后面三天

很多人拿到一份叫「风工程 OpenFOAM 基础培训」的讲义，翻到「OpenFOAM 简介」那几页，看到终端、字典文件、编译命令，第一反应是先跳过，直接找算例抄。这恰恰是最容易把后续几天工作白费的地方。OpenFOAM 不是打开就能算的图形软件，它对一套建筑风环境算例的要求是：计算域尺寸、地面粗糙度、入口风剖面、湍流模型、壁面函数、离散格式六件事必须自洽。任何一环给错，残差照样往下掉，云图照样花花绿绿，结论却是错的。这几页简介真正要交代的，是这套自洽关系怎么落地。做建筑风压、行人高度风速、小区风环境评估的人，以及从图形界面软件转过来、想自己搭一条可复用流程的人，都绕不开它。

2. OpenFOAM 的算例骨架与求解器选择：0 / constant / system 各管什么

OpenFOAM 的每个算例在磁盘上都是一棵固定形状的目录树，风工程算例也不例外。先把这棵树的职责分清楚，比急着改参数有用得多。根目录下有三个文件夹：0、constant、system，还有一个跑起来之后才会生成的processor*或者时间步目录。新手最常见的错误是分不清某个参数该写进哪个文件，于是反复重启计算。

2.1 从 0 目录开始：U、p、k、epsilon 四个文件的边界条件长什么样

0目录放的是初始场和边界条件，一个物理量一个文件。风工程稳态计算至少需要U（速度）、p（压力）、k（湍动能）、epsilon或omega（耗散率或比耗散率）、nut（湍流运动黏度）。文件结构三段式：dimensions写量纲，internalField写内部初值，boundaryField逐个 patch 写边界类型。

U的典型写法如下，注意入口用fixedValue、出口用inletOutlet、建筑与地面用noSlip、顶部与侧面用slip：

// 0/U dimensions [0 1 -1 0 0 0 0]; internalField uniform (0 0 0); boundaryField { inlet { type fixedValue; value uniform (10 0 0); // 先给常数跑通，再换成风剖面 } outlet { type inletOutlet; // 回流时按入口值处理，防止出口倒灌发散 inletValue uniform (0 0 0); value uniform (0 0 0); } ground { type noSlip; } building { type noSlip; } top { type slip; } // 顶部不能设 wall，否则相当于加了盖 sides { type slip; } }

p这一块有一个风工程里几乎人人踩过一次的点：不可压求解器里的p是运动压力，也就是压力除以密度，量纲是[0 2 -2 0 0 0 0]，单位是 m²/s²。很多人按静压 101325 去设定出口值，结果整个压力场被抬到一个毫无意义的量级，云图看起来还挺正常。正确做法是出口fixedValue给 0，其余zeroGradient，后面算风压系数时直接Cp = p / (0.5 * Uref²)，不用再除密度。

// 0/p dimensions [0 2 -2 0 0 0 0]; internalField uniform 0; boundaryField { inlet { type zeroGradient; } outlet { type fixedValue; value uniform 0; } ground { type zeroGradient; } building { type zeroGradient; } top { type zeroGradient; } sides { type zeroGradient; } }

k和epsilon入口值需要按湍流强度和特征尺度反推，第 4 章展开。nut在壁面上要选壁面函数类型，比如nutkWallFunction或nutUSpaldingWallFunction，前者对 y+ 落在对数律区更稳，后者对粗糙壁面和过渡区更宽容。

2.2 constant 与 system：网格、物性、离散格式、求解控制分开放

constant装的是随时间不变的东西：polyMesh/里的网格（points、faces、owner、neighbour、boundary五个文件），运动黏度nu，以及湍流模型配置。湍流模型字典在不同发行版里名字不统一，常见的是turbulenceProperties、RASProperties或momentumTransport，内容大同小异：

// constant/turbulenceProperties（名字随发行版可能不同） simulationType RAS; RAS { RASModel kEpsilon; turbulence on; printCoeffs on; }

system装的是求解控制：controlDict管时间推进与输出，fvSchemes管离散格式，fvSolution管线性求解器与松弛因子，blockMeshDict、snappyHexMeshDict、decomposeParDict也都放在这里。风工程稳态算例的controlDict关键项就这么几行，其中deltaT在稳态求解器里只是伪时间步，不影响物理结果，但会影响残差下降的平滑度：

application simpleFoam; startFrom startTime; startTime 0; stopAt endTime; endTime 2000; deltaT 1; writeControl timeStep; writeInterval 100; purgeWrite 3; // 只留最近几个时间步，省磁盘

fvSchemes里对风工程影响最大的是divSchemes。一上来全用linear看起来精度高，钝体绕流很容易震荡；全用upwind稳定但数值耗散大，尾流被抹平，风压系数偏低。折中做法是速度用linearUpwind，湍流量用upwind：

divSchemes { default none; div(phi,U) Gauss linearUpwind grad(U); div(phi,k) Gauss upwind; div(phi,epsilon) Gauss upwind; div((nuEff*dev2(T(grad(U))))) Gauss linear; } laplacianSchemes { default Gauss linear corrected; } snGradSchemes { default corrected; }

fvSolution里，SIMPLE段要设nNonOrthogonalCorrectors，网格非正交度大就设 2；relaxationFactors里压力压到 0.3，速度 0.7，湍流量 0.7，这是风工程稳态计算比较通用的起手值。较新的发行版把松弛因子分成fields和equations两组，压力归前者，其余归后者。

2.3 simpleFoam、pimpleFoam、pisoFoam 怎么挑：风工程的三条判断线

求解器选错，后面所有参数都是白调。判断线其实就三条：要不要时间项、要不要浮力、要不要求解能量方程。

绝大多数建筑风环境评估做的是平均量，simpleFoam就够。要算旋涡脱落、峰值负压、脉动特性，得上瞬态；网格尺度到厘米级、时间步到毫秒级，机时成本立刻上涨一个数量级，所以先想清楚报告里到底要出哪个数。定稳态还是瞬态这件事，通常要在算之前定下来，算完再改就是重来。

用potentialFoam给初值是性价比很高的一步：先跑几十步势流，把U和p大致铺开，再切到simpleFoam，能省掉前面几百步的摸索期。做法是potentialFoam读完同一个0目录、把结果写回0或者某个时间步，然后改startFrom latestTime。

3. 把一栋建筑塞进计算域：blockMesh 与 snappyHexMesh 的实操参数

风工程网格分两步走：先用blockMesh造一个规整的长方体背景网格，再用snappyHexMesh把建筑 STL 表面“切”进去并加边界层。这套流程稳定、可脚本化，是 OpenFOAM 做钝体绕流绕不开的一步。

3.1 计算域多大才不算白算：阻塞率与上下游距离的取法

计算域尺寸不是随便定的。以建筑特征高度 H 为基准，常见取法是：上游 3H～5H，下游 10H～15H，侧面各 5H，顶部 5H～6H。下游要留长是因为尾流恢复慢，取短了出口边界会反过来影响建筑背风面的压力。

阻塞率是另一个硬指标：建筑迎风面投影面积除以计算域迎风截面积，一般控制在 3% 以内。阻塞率超了，来流被“挤”着加速，建筑表面风压整体偏高，而且这种偏差不会随网格加密而消失。一个 100 m 高的塔楼，计算域宽度和高度至少要做到 400 m 量级，这不是浪费网格，是保证结果可解释。

3.2 blockMeshDict 写背景网格：vertices、blocks、boundary 三块

blockMeshDict的结构是：convertToMeters定单位，vertices列 8 个角点，blocks声明六面体及其三向网格数，boundary给六个面命名并指定 patch 类型。风工程里通常还会用blocks的分段写法，让建筑附近密、远场疏：

convertToMeters 1; vertices ( (-150 -100 0) (450 -100 0) // 上游 150，下游 450 (450 200 0) (-150 200 0) (-150 -100 400) (450 -100 400) (450 200 400) (-150 200 400) ); blocks ( hex (0 1 2 3 4 5 6 7) (60 40 40) // 三向网格数 simpleGrading (5 5 8) // 三向渐变：近地面、近建筑加密 ); boundary ( inlet { type patch; faces ((0 4 7 3)); } outlet { type patch; faces ((1 2 6 5)); } ground { type wall; faces ((0 3 2 1)); } top { type patch; faces ((4 5 6 7)); } sides { type patch; faces ((0 1 5 4) (3 7 6 2)); } );

simpleGrading的参数是每个方向末段与首段的网格尺寸比，大于 1 表示从起点向终点变疏。上面(5 5 8)的意思是水平两向从上游向下游放大 5 倍，竖向从地面向顶部放大 8 倍。竖向不要放太大，顶部网格太粗会让顶部边界的滑移条件传递出数值误差。写完执行：

blockMesh

终端会打印网格点数、单元数和边界 patch 列表，先确认单元数在预期量级，再往下走。这一层网格不用太细，它只是给snappyHexMesh提供骨架，真正的分辨率由后面加密级别决定。

3.3 snappyHexMesh 三步走：castellated、snap、addLayers 的关键参数

snappyHexMesh一次运行里其实跑了三个阶段，每个阶段一个控制块，全部在system/snappyHexMeshDict里。建筑 STL 放在constant/triSurface/下，命名以.stl结尾，字典里引用时去掉后缀。

第一阶段 castellated，把背景网格按表面附近逐级细化并“挖掉”建筑内部：

castellatedMesh true; snap true; addLayers true; geometry { building.stl { type triSurfaceMesh; name building; } } castellatedMeshControls { maxLocalCells 2000000; maxGlobalCells 8000000; minRefinementCells 10; nCellsBetweenLevels 3; // 相邻加密级之间至少隔 3 层，防止跳级 resolveFeatureAngle 30; // 曲面夹角超过 30 度就加密，屋檐、棱角靠它 locationInMesh (0 0 50); // 这个点必须在流体域内，在建筑内部就全挖没了 features ( { file "building.eMesh"; level 3; } ); refinementSurfaces { building { level (3 4); // 最小 3 级，最大 4 级 patchInfo { type wall; } } } refinementRegions { buildingRegion { mode inside; levels ((1e-6 5)); // 建筑附近再补一层 } } }

locationInMesh是最容易翻车的一项：它决定“哪一侧是流体”。点选在建筑内部，整个域会被判成固体，网格挖空，后面所有求解都跑不起来。

第二阶段 snap，把阶梯状的表面顶点吸附到 STL 上，这决定建筑表面是否光滑：

snapControls { nSmoothPatch 3; tolerance 2.0; // 相对局部网格尺寸的吸附容差 nSolveIter 30; nRelaxIter 5; nFeatureSnapIter 10; implicitFeatureSnap false; }

第三阶段 addLayers，给壁面加边界层网格，直接影响壁面剪切和风压精度：

addLayersControls { relativeSizes true; layers { building { nSurfaceLayers 5; } } expansionRatio 1.2; // 相邻层厚比，超过 1.3 容易畸变 finalLayerThickness 0.3; // 相对第一层之外的末层厚度 minThickness 0.1; // 低于此值就不加了，避免负体积 featureAngle 60; slipFeatureAngle 30; nGrow 0; nRelaxIter 5; }

加层失败是常态，不是例外。终端会打印“X out of Y faces extruded successfully”，挤出成功率低于 90% 就要回头降nSurfaceLayers或者调featureAngle。

3.4 checkMesh 输出怎么读：四条不能越过的红线

网格跑完必须过一遍checkMesh，四项指标要盯：

checkMesh -allGeometry -allTopology | tee log.checkMesh

把输出重定向到文件是血泪经验：出问题时回头翻日志，比凭记忆猜快得多。checkMesh报Failed 1 mesh checks不代表一定要重做，要看具体是哪个，非正交角 66 度通常能靠nNonOrthogonalCorrectors 2扛过去，偏斜度超 6 就老老实实回去改snapControls。

4. 大气边界层入口与湍流模型：风工程最容易糊弄过去的两组参数

同样的网格，入口风剖面给常数还是给指数律，建筑顶部的风速能差一倍；湍流模型选k-epsilon还是k-omega SST，背风面的负压系数能差 30%。这两组参数是风工程 OpenFOAM 计算里最值得花时间的地方。

4.1 指数律风剖面：codedFixedValue 逐面给风速

真实大气边界层里，平均风速随高度按指数律增长：U(z) = Uref * (z / zref)^alpha。alpha是地面粗糙度指数，开阔场地常见 0.11～0.16，郊区 0.16～0.22，城市中心 0.22～0.33。这个值给错，等于换了一个场地在算。

入口用codedFixedValue逐面赋值，代码在字典里直接编译：

// 0/U 的 inlet 段 inlet { type codedFixedValue; value uniform (0 0 0); name ablInlet; code #{ const fvPatch& p = patch(); const vectorField& Cf = p.Cf(); // 每个面心的坐标 vectorField& U = *this; const scalar Uref = 10.0; // 参考高度风速 const scalar zref = 10.0; // 参考高度 const scalar alpha = 0.22; // 郊区地貌 forAll(Cf, faceI) { const scalar z = max(Cf[faceI].z(), 0.5); U[faceI] = vector(Uref*pow(z/zref, alpha), 0, 0); } #}; }

这段代码的逻辑是遍历入口 patch 的每个面心，读它的 z 坐标，按幂律算出该面的速度矢量。max(..., 0.5)是为了防止最底层网格面心 z 接近 0 时幂函数把速度压到 0。参数调整只改三行：Uref、zref、alpha，改完不需要重新编译整个求解器，codedFixedValue会自动重编。如果第一次用不确定接口，先用uniform (10 0 0)跑通，再换代码，能省掉排查编译错误的时间。

4.2 k 与 epsilon 的入口值：湍流强度怎么反推

入口湍动能按定义反推：k = 1.5 * (U * I)²，I是湍流强度。风工程里常按高度给一个衰减型剖面，简化为常数也能接受，但不要给成 0.001 这种量级，太小的入口湍流会让尾流恢复得极慢。

耗散率按epsilon = Cmu^0.75 * k^1.5 / l给，Cmu取 0.09，l是湍流积分尺度，钝体绕流常取 0.07 倍的特征长度。用k-omega SST时对应omega = k^0.5 / (Cmu^0.25 * l)。这两个式子不难，难的是别把单位搞混，k的量纲是[0 2 -2 0 0 0 0]，和p一样，属于运动量。

4.3 k-epsilon 与 k-omega SST 的选型：附着流、分离流、钝体绕流

做建筑表面平均风压，k-omega SST通常是更稳妥的起点；做大范围场地、几十栋楼的流场，realizable k-epsilon收敛快、机时省。两边都要壁面函数配合，见下一节。

4.4 壁面函数与 y+：地面这一刀切在哪

y+ 是无量纲壁面距离，决定壁面第一层网格落在黏性底层、过渡区还是对数律区。用壁面函数时，第一层网格要落在 30 < y+ < 300；想解析到壁面，y+ 要压到 1 以下，代价是第一层网格厚度降到毫米量级，单元数翻几倍。

风工程地面通常是粗糙壁面，用nutkWallFunction加粗糙度参数是常见做法：在nut的壁面段里给Ks（等效砂粒粗糙度）和Cs，Ks与地面粗糙度指数对应，城市地貌的Ks能到 1 m 量级。给了粗糙度，第一层网格就不用做到毫米级，y+ 落在对数律区即可。

// 0/nut 的 ground 段 ground { type nutkWallFunction; value uniform 0; Ks uniform 0.5; // 等效砂粒粗糙度，按地貌取 Cs uniform 0.5; // 粗糙度常数，一般 0.5～1.0 }

建筑表面一般按光滑壁面处理，Ks给 0 或很小值。侧面和顶部必须是slip或symmetry，如果误写成wall，等于给计算域加了个盖子，来流被压着加速，屋顶风速会明显偏高，而且残差曲线看不出异常，只能靠对比发现。

5. 排查清单：跑不出像样结果时先翻这五条

风工程算例出问题，八成不是求解器本身，而是边界、网格或参数之间不自洽。下面五条按出现频率排，每条给出可以直接执行的处置顺序。

5.1 残差很漂亮，建筑表面 Cp 却出现 ±5 的怪值

现象：U、p残差都降到 1e-5 以下，但建筑表面风压系数出现正负 5 以上的尖刺，位置集中在檐口和角部。

原因：snappyHexMesh生成的表面存在台阶状细节，局部偏斜度大，面法向压力梯度计算失真；也可能是 layers 挤出失败，壁面处网格尺寸突变。

解决：先跑checkMesh -allGeometry看最大偏斜度，超 4 就回snappyHexMeshDict调resolveFeatureAngle（降到 20 度让棱角加密）；再看 layers 挤出成功率，低于 90% 时把nSurfaceLayers从 5 降到 3；最后确认fvSchemes里snGradSchemes是corrected而不是uncorrected。

5.2 第一步就发散，U 残差直接飙到 1 以上

现象：第一个时间步残差就跳到 1 量级，几轮之后浮点溢出退出。

原因：入口速度、湍流参数和压力边界互相打架。最常见的是入口p设了fixedValue、出口也设了fixedValue，压力场被两端夹死；或者epsilon入口给了个极大的值，湍流黏性直接爆掉。

解决：入口p用zeroGradient，出口p用fixedValue 0，全流场只有一个压力参考点；k和epsilon先给常数、量级正常（k在 0.1～1 量级，epsilon在 0.001～0.1 量级），跑通再换剖面；松弛因子先压到p 0.2 / U 0.5，稳定几百步后再放开。

5.3 屋顶风速比来流还大，出口像堵了墙

现象：建筑顶部和两侧的风速接近或超过入口参考风速，流线在顶部明显上翘。

原因：顶部或侧面 patch 被误设成wall、noSlip，或者blockMeshDict里boundary的面顶点顺序写反，导致法向朝内。还有一种情况是出口用了zeroGradient而不是inletOutlet，回流时无法正确出流。

解决：先foamToVTK看U的边界矢量方向，确认slip生效；出口U改成inletOutlet；计算域顶部高度从 5H 加到 6H，阻塞率重算一遍，控制在 3% 以内。

5.4 并行结果和串行对不上

现象：4 核并行跑出来的建筑表面平均风压，比单核跑低 10% 以上。

原因：decomposePar之后改了0目录的边界条件却没有重新分区，processor*里还是旧的0；或者decomposeParDict里numberOfSubdomains和mpirun -np的核数不一致，多出来的进程在算空网格。

解决：改完任何0或constant里的内容，decomposePar -force重跑一次；并行前用reconstructPar -latestTime做一次串并行对照，取一个监测点比较，偏差超过 1% 就说明分区或边界没同步。

decomposePar -force mpirun -np 4 simpleFoam -parallel | tee log.simpleFoam reconstructPar -latestTime

5.5 跑了五千步不收敛，监测量一直小幅震荡

现象：残差卡在 1e-3 上下反复，建筑表面的力系数和某点风速在 ±2% 内来回摆，看不出趋势。

原因：SIMPLE的松弛因子偏大，或者网格非正交度较高但nNonOrthogonalCorrectors只给了 0；也可能是尾流本身存在低频非定常，稳态求解器本来就不该追求残差单调下降。

解决：nNonOrthogonalCorrectors加到 2，压力松弛降到 0.2；用functionObject监控力系数和特征点速度，看最后 500 步的波动幅度，稳定在 1% 以内就可以停，不必强求残差到 1e-6。真怀疑非定常，就换pimpleFoam算几百个时间步看尾流是否周期性摆动。

6. 结果可信度怎么验：残差、监测量、网格无关性三层校核

算完不是结束，是开始怀疑的时候。残差只是数值收敛的旁证，不能当作物理正确的证据。我一般的顺序是：先看监测量是否稳定，再做网格无关性，最后和规范或风洞数据对一对。这三层里任何一层不过，前面的云图都不该往报告里放。

6.1 从残差到力系数：用 function object 盯住特征量

在controlDict里挂一个forces函数对象，每若干步输出一次建筑整体力和力矩，比盯残差有用得多：

functions { buildingForces { type forces; libs ("libforces.so"); patches (building); rho rhoInf; rhoInf 1.225; CofR (0 0 0); writeControl timeStep; writeInterval 10; } }

patches必须是建筑壁面 patch 名，CofR是力矩参考点，一般取建筑底面中心。输出在postProcessing/forces/下，取最后几百步做平均，看阻力和升力系数的漂移量，1% 以内算稳。

6.2 网格无关性与 Cp 换算：一张表看清差异

网格无关性至少做三套，单元数按 1.3～1.5 倍递增，比较同一个特征量。记录格式如下，数值是示意，实际要填自己的结果：

偏差降到 5% 以内、再加密变化不明显，那套网格就可以定稿。风压系数按Cp = p / (0.5 * Uref²)算，注意这里的p是运动压力，别再多除一次密度，这是从图形界面软件转过来的人最容易犯的错。

换算完和荷载规范或风洞试验数据对一对，重点看迎风面中心和背风面中心两个点，角部和檐口允许有更大偏差，因为那些位置对网格和湍流模型都极敏感。

我踩过的最大一个坑，是花两天调一张漂亮的流线图，最后发现locationInMesh点在了建筑内部，整个计算域其实算的是一团静止空气。从那以后我的习惯是：snappyHexMesh跑完先看checkMesh的单元数，再截一张 z 中截面的网格图，确认建筑是空的、外面是实的。这个动作十秒钟，能省掉两天。希望帮到你。

本文还有配套的精品资源，点击获取