简介这份文献是面向风工程与计算流体力学初学者的OpenFOAM基础培训讲义适合刚接触开源CFD、希望建立软件整体认识的读者。压缩包仅含1个PDF文件约162KB篇幅轻量便于快速浏览与课堂配套学习。已有555人学习。内容从OpenFOAM的开发背景与版本演进讲起梳理其面向对象C类库、有限体积离散、非结构化多面体网格、大规模并行计算等核心能力并延伸到自动网格运动、网格拓扑变化、拉格朗日粒子追踪、有限元与有限面积离散等扩展功能。讲稿还介绍元类、容器类、字典、实时选择表与分层设计理念可帮助读者理解求解器框架、代码重用和用户定制思路为后续阅读源码、搭建风工程算例和排查网格与并行设置问题提供入门线索。1. 风工程为什么用 OpenFOAM一次阻塞比返工接手一个屋面风压复核任务时模型是别人拷过来的网格看着很规整simpleFoam收敛得也漂亮。结果迎风面中部 Cp 算到 1.3屋角吸力峰值只有 -1.1跟风洞报告对不上。查了两天问题出在计算域宽度只留了 3 倍建筑高度阻塞比接近 12%气流被硬挤过去整个压力场都是歪的。这类返工在风工程里太常见了OpenFOAM 本身不难难的是把风洞里的物理约束翻译成命令行里的参数。风工程 OpenFOAM 基础培训要解决的正是这件事把 OpenFOAM 当作一台计算风洞来用而不是当通用 CFD 软件来用。它适合结构风荷载、建筑群风环境、行人高度风舒适度这类以钝体绕流为核心的问题也适合已经会用商业软件、但说不清湍流入口和 y 到底该怎么定的人。2. 求解器与湍流模型怎么挑从 simpleFoam 到 LES 的选型边界2.1 不可压求解器分工simpleFoam、pisoFoam、pimpleFoam风工程里绝大多数算例的风速远低于 0.3 马赫空气按不可压处理密度当成常数。真正要选的是时间项怎么处理。simpleFoam走 SIMPLE 算法把时间项设成steadyState一轮迭代直接求解平均流场成本最低方案比选、多风向角扫掠、参数敏感性分析都靠它。它的局限也很明确稳态假设下没有涡脱落钝体尾流里那些周期性的东西会被抹平成一个看起来还行的定常解。pisoFoam走 PISO 算法真正的时间推进配合 LES 或 DES 用。它要求 Courant 数压在 1 以下时间步长往往在 1e-4 秒量级跑十几秒物理时间就是几十万步成本是稳态的几十到上百倍。值不值得取决于你要的是平均风压还是峰值吸力和脉动风压谱——后者只能靠它。pimpleFoam走 PIMPLEPISO 和 SIMPLE 的混合允许 Courant 数大于 1可以只做一次外迭代校正。非定常 RANS、带热浮力的街区通风、以及稳态不收敛想试试非定常的过渡场景我一般先拿它试水。求解器时间项风工程典型任务相对成本simpleFoamsteadyState平均风压、平均风速场、方案比选1pimpleFoam瞬态可大 Courant非定常 RANS、热浮力、过渡态1030pisoFoam瞬态Courant 1LES/DES涡脱落、脉动风压50200potentialFoam无给稳态求解器一个像样的初场 1较新的发行版把求解器合并成了统一的foamRun -solver incompressibleFluid功能一致只是入口不同。选型逻辑不变先问自己要平均量还是瞬时量。2.2 湍流模型k-ε 稳、k-ω SST 贴壁、LES 抓涡标准 k-ε 是最省事的选择。它对自由剪切层和远场比较宽容收敛性好工程上算整体风压分布、街区平均风速场够用。缺点写在教科书里也写在实践里钝体背风侧的分离区会被过度预测再附点往后拖背压偏低。某跨平台系统做过方柱绕流的对比标准 k-ε 给出的再附长度比实验长约 20%realizable k-ε 能收到 10% 以内。k-ω SST 在近壁区用 k-ω远场切回 k-ε分离预测比 k-ε 准是钝体绕流里性价比最高的一档。代价是对网格更挑想要它发挥优势壁面第一层要压到 y ≈ 1如果加壁面函数则要保持 y 在 30100最怕停在 530 的缓冲层里两头不靠。LES 是另一档。它能直接解析大尺度涡给出峰值吸力和脉动风压谱屋面围护结构设计、行人高度阵风这些都靠它。但网格量按 Re^1.8 左右的趋势涨一栋 60 米高的楼按 1:100 缩尺做 LES网格轻松过千万还得跑够统计时间。我的习惯是先稳态拿平均场定方案再挑一两个关键风向做 LES 精细核算不全程上 LES。提示湍流模型选了 k-ε 就别指望背风面 Cp 精确到小数点后两位用它做相对比较A 方案和 B 方案谁更差比做绝对值标定靠谱得多。2.3 最小算例骨架0/ constant/ system/ 三个目录OpenFOAM 的算例结构是死的风工程算例也一样。摸清这三个目录剩下的都是填内容。# 从自带算例拷一个不可压模板改成自己的 cp -r $FOAM_TUTORIALS/incompressible/simpleFoam/pitzDaily ./caseWind cd caseWind # 目录骨架 # 0/ U p k epsilon nut —— 初始场与边界条件 # constant/ transportProperties —— 运动黏度 nu # turbulenceProperties —— 湍流模型开关 # polyMesh/ —— blockMesh 产出的背景网格 # triSurface/ —— 建筑 STL # system/ blockMeshDict snappyHexMeshDict # controlDict fvSchemes fvSolutiontransportProperties里只关心nu20 摄氏度空气取 1.5e-5 m²/s缩尺模型不改 nu改的是特征尺寸和风速。turbulenceProperties里simulationType RAS配RASProperties选模型LES则换LESProperties。fvSchemes和fvSolution是数值层面的东西稳态钝体绕流有一套相对固定的写法下面会贴。真正决定结果对不对的是0/目录边界条件写错后面调什么都没用。3. 计算域与网格把阻塞比和 y 先在纸上算清楚3.1 域尺寸的硬约束上游 5H、下游 15H 与 3% 阻塞比域尺寸不是凭感觉画的。上游一般留 5HH 为建筑高度让来流有空间发展下游留 15H尾流有足够长度恢复否则出口处还有回流pressureInletOutletVelocity会被逼着往计算域里灌气残差就开始抖侧面单侧 5H、顶部 56H保证顶部不出现明显的加速。阻塞比是更硬的约束。风洞实验的通行门槛是迎风投影面积不超过试验段截面的 3%数值模拟可以放宽到 5%但超过 5% 之后迎风面 Cp 会系统性偏高屋角吸力偏低误差可以到 20% 以上。# 由建筑尺寸反算最小计算域宽度保证阻塞比不超标 H, B, D 60.0, 40.0, 30.0 # 建筑高、迎风宽、进深单位 m blockage 0.03 # 风洞 3%纯数值可放宽到 0.05 A_front B * H # 迎风面投影面积 A_cross A_front / blockage # 所需最小横截面积 H_domain 6.0 * H # 域高取 6H W_domain A_cross / H_domain # 由面积反算域宽 print(f迎风面积 {A_front:.0f} m2截面至少 {A_cross:.0f} m2) print(f域高 {H_domain:.0f} m 时总宽至少 {W_domain:.0f} m单侧 {W_domain/2:.0f} m) print(f建议域范围 x:[-{5*H:.0f}, {15*H:.0f}] y:±{W_domain/2:.0f} z:[0, {H_domain:.0f}])代入 60 米高、40 米宽的建筑3% 要求总宽至少 222 米单侧 111 米。这个数字很多人第一眼觉得夸张但它就是风洞准则换算过来的结果。实际做的时候我会先按 5% 估一版跑完看迎风面 Cp 是否明显高于 1.0驻点理论上限偏高就加宽。注意域加宽不是免费的网格量按面积涨。宽域配粗网格和窄域配细网格前者通常更准。3.2 blockMesh 背景网格关键行与 checkMesh 判据blockMesh只负责铺一个六面体背景网格建筑由snappyHexMesh挖出来。所以背景网格不用太细均匀一点、正交性好一点就行真正的加密交给后面。// system/blockMeshDict 节选H 60 m 的算例 vertices ( (-300 -150 0) // 0 ( 900 -150 0) // 1 ( 900 150 0) // 2 (-300 150 0) // 3 (-300 -150 360) // 4 ( 900 -150 360) // 5 ( 900 150 360) // 6 (-300 150 360) // 7 ); blocks ( hex (0 1 2 3 4 5 6 7) (120 60 60) simpleGrading (1 1 3) // z 向渐变地面附近密、高空疏 );simpleGrading (1 1 3)的意思是 x、y 方向均匀z 方向最后一层与第一层厚度比为 3。大气边界层的速度梯度集中在地面附近这样铺能在不增加网格量的前提下把地面分辨率提上去。三个方向都用单一渐变会限制加密位置真要做建筑附近密、远处疏得多块拼或者干脆靠snappyHexMesh的refinementRegions补。boundary段里定义inlet、outlet、ground、top、side五个 patch类型先都写patch具体边界条件在0/里再定。生成完立刻跑checkMeshblockMesh log.blockMesh 21 checkMesh log.checkMesh 21 grep -E non-orthogonality|skewness|aspect ratio log.checkMesh三个数字要盯最大非正交角最好小于 65 度平均小于 40 度最大偏斜小于 4长宽比在边界层里大是正常的但背景网格层面别超过 1000。这三个指标失控后面snappyHexMesh贴出来的层也会跟着歪。3.3 snappyHexMesh 贴体加密三层加密与首层厚度snappyHexMesh是风工程网格的主力。流程是surfaceFeatureExtract提特征线snappyHexMesh挖洞、加密、贴层最后checkMesh复核。// system/snappyHexMeshDict 关键段 castellatedMesh true; snap true; addLayers true; geometry { building.stl { type triSurfaceMesh; name building; } } castellatedMeshControls { maxLocalCells 4000000; maxGlobalCells 20000000; minRefinementCells 10; nCellsBetweenLevels 3; // 相邻加密级之间至少隔 3 层防跳变 resolveFeatureAngle 30; features ( { file building.eMesh; level 3; } ); refinementSurfaces { building { level (2 3); } // 最小 2 级、最大 3 级 } locationInMesh (0 0 200); // 取一个肯定在流体域内的点 allowFreeStandingZoneFaces true; } addLayersControls { layers { building { nSurfaceLayers 5; } ground { nSurfaceLayers 3; } } expansionRatio 1.2; // 相邻层厚度比别超过 1.3 firstLayerThickness 0.002; // 首层 2 mm由 y 反算得来 minThickness 0.0004; // 最小层厚低于它就塌陷 featureAngle 60; nGrow 0; }level (2 3)表示最小加密 2 级、特征线附近 3 级背景网格 5 米的话最细能到 0.6 米。nCellsBetweenLevels 3是个容易被忽略但很关键的数加密级之间跳得太快非正交角会飙snap阶段就容易出负体积。firstLayerThickness是绝对厚度单位米必须由 y 反算不能拍脑袋。minThickness一般取首层的 0.10.2 倍太小了不起作用太大会直接导致层塌陷。surfaceFeatureExtract log.surfaceFeatureExtract 21 snappyHexMesh -overwrite log.snappyHexMesh 21 checkMesh log.checkMesh 21 grep -i collapsed\|negative log.snappyHexMesh-overwrite省掉一次改名但出错就回不去了第一次跑建议不加。日志里出现collapsed或negative volume先看minThickness和nSurfaceLayers绝大多数是这两个数给大了。3.4 第一层高度反算一张 y 对照表和两个公式y 的定义是 y uτ·y/ν其中 uτ 是壁面摩擦速度y 是首层高度ν 是运动黏度。工程上 uτ 可以用 0.05·Uref 粗估够用了。按 ν 1.5e-5 m²/s 算Uref (m/s)uτ 估计 (m/s)y1 首层y30 首层y100 首层50.250.06 mm1.8 mm6 mm100.500.03 mm0.9 mm3 mm150.750.02 mm0.6 mm2 mm201.000.015 mm0.45 mm1.5 mm表里最扎眼的是一列比一列小。想全解析壁面y ≈ 110 m/s 风速下首层只要 0.03 毫米一栋 60 米的楼配上这个首层边界层网格量能顶到天上去。这就是为什么工程 RANS 普遍走壁面函数路线把首层放在 y 30100取个 13 毫米层数 5 层、膨胀比 1.2。# 由目标 y 反算首层高度反过来也行 nu, Uref, yplus 1.5e-5, 10.0, 50.0 u_tau 0.05 * Uref y_first yplus * nu / u_tau print(fu_tau{u_tau:.3f} m/s, y{yplus} 对应首层 {y_first*1000:.2f} mm) # 已知首层验算实际 y y_given 0.002 print(f首层 {y_given*1000:.1f} mm 实际 y {u_tau*y_given/nu:.0f})算完之后别急着信跑起来用yPlusfunctionObject 输出一遍实际分布。模型表面的 y 是随位置变的迎风面驻点附近低、屋角高能差三倍。看的是绝大部分面积落在哪个区间不是平均值。4. 大气边界层入口怎么写速度剖面、粗糙度与出口配置4.1 幂律剖面的 codedFixedValue 写法风洞里的来流不是均匀风是带速度梯度的大气边界层。直接给fixedValue uniform (10 0 0)是最常见的错误驻点压力勉强对顶部加速区和分离剪切层全错因为整个流场没有剪切的驱动力。幂律是最简单的一种写法用codedFixedValue直接按面心坐标算inlet { type codedFixedValue; value uniform (0 0 0); name powerLawInlet; code #{ const scalar Uref 8.0; // 参考高度风速 m/s const scalar Zref 10.0; // 参考高度 m const scalar alpha 0.16; // 地貌指数 const vector dir(1, 0, 0); // 来流方向 const fvPatch p patch(); const vectorField Cf p.Cf(); // 每个面心的坐标 vectorField U *this; // 当前边界上的 U forAll(Cf, i) { const scalar z max(Cf[i].z(), 0.0); U[i] dir * Uref * pow(z / Zref, alpha); } #}; }三个参数里alpha决定剖面形状开阔水面 0.12开阔草地 0.140.16郊区低层建筑 0.20密集城区 0.250.30 都有人用。Uref和Zref是配套的别把 10 米参考风速配 1 米参考高度那算出来的剖面整体偏大一个量级。缩尺模型要按缩尺比换算几何按米建模、风速也按对应的速度比换算nu保持真实值不变。提示codedFixedValue每次求值都要重跑一遍循环边界面上万面时开销可观。面数特别大的话改用timeVaryingMappedFixedValue读一张预先写好的剖面表更省。4.2 k 与 epsilon 入口湍流强度 10% 是怎么来的只给速度剖面湍流是零边界层是层流化的分离剪切层起不来模型的侧面和背风面会明显失真。入口必须同时给 k 和 epsilon或者 omega。k 由湍流强度反算k 1.5·(U·I)²。风洞大气边界层的湍流强度在参考高度附近大致 8%15%取 10% 是常见起点。epsilon 由混合长度反算混合长度取边界层厚度的 0.070.1 倍。// 0/k inlet { type turbulentIntensityKineticEnergyInlet; intensity 0.10; // 10% 湍流强度 value uniform 1.2; // 8 m/s、10% 时约 1.2 m2/s2 } // 0/epsilon inlet { type turbulentMixingLengthDissipationRateInlet; mixingLength 20.0; // 约 0.1 倍边界层厚度 m value uniform 0.05; }intensity和value的关系是 k 1.5(U·I)²求解器会用value作为初值、按intensity随时间更新。两个数不一致会导致第一步的 k 场有个跳变虽然能收敛但日志难看。用 k-ε 有个绕不开的问题幂律速度剖面配常数 k、epsilon在流场里不是平衡态往下游走剖面会自己衰减模型位置上根本到不了目标风速。较新的发行版内置了atmBoundaryLayer系列边界条件按对数律同时给 U、k、epsilon能保持平衡。字段名各分支略有差异落到自己环境里先用源码或foamSearch确认一遍再写。风洞里测了剖面的老老实实用timeVaryingMappedFixedValue读表那是最保真的做法。4.3 地面粗糙度roughnessLength 与壁面函数的配对地面是 wall但不能当光滑壁面处理。城市下垫面的粗糙度直接影响底部风速和街谷里的流动。做法是给地面的nut用粗糙壁面函数把空气动力学粗糙度长度 z0 换算成等效砂粒粗糙度 Ks// 0/nut ground { type nutkRoughWallFunction; Ks uniform 1.0; // 等效砂粒粗糙度由 z0 换算 Cs uniform 0.5; // 粗糙度常数 value uniform 0; }换算式是 Ks ≈ 9.793·z0/Cs。z0 取 0.05 米郊区/低层城区Cs 取 0.5得 Ks ≈ 1.0 米。这个值看着大但它不是几何尺寸是一个等效参数别拿去跟首层高度比。配套的还有0/k地面用kqRWallFunctionvalue 设 0 即可0/epsilon地面用epsilonWallFunction。k 和 epsilon 的壁面处理没配对最常见的表现是残差降不下去、地面附近速度出现非物理的加速。如果是 LES地面就没这么简单了得用nutUSpaldingWallFunction或者干脆把首层压到 y ≈ 1 全解析。后者成本高通常只在关键区域模型周边这么做。4.4 出口、顶部、侧面的边界类型与回流风险出口的经典组合是速度用pressureInletOutletVelocity压力用fixedValue uniform 0。前者在出流时是零梯度一旦局部出现回流就切回固定值比单纯的zeroGradient稳得多是尾流算例的标配。顶面和侧面用slip速度压力zeroGradient。用slip而不是 wall是为了模拟风洞试验段的自由滑移边界避免侧壁面额外生成边界层把流场往中间挤。有人喜欢用symmetry在纯速度场上和slip等价看习惯。// 0/U 的出口段 outlet { type pressureInletOutletVelocity; value uniform (0 0 0); } // 0/p 的出口段 outlet { type fixedValue; value uniform 0; }出口离模型太近时回流会一直顶到模型背风面表现为背压忽高忽低、残差周期性抖动。这时候要么把出口往下游推要么在出口前加一段缓变的多孔区工程上少用不推荐最简单的还是回到 3.1 节把下游长度补到 15H。5. 避坑排查五个让算例翻车的现场5.1 y 停在 530 缓冲层壁面函数和解析壁面两头不靠现象模型表面 Cp 分布形状对但整体偏移迎风面偏高、屋角吸力偏低yPlus输出显示表面大部分区域落在 530。原因这一层既不能当对数律区用壁面函数也没细到能解析黏性底层两种壁面处理都在这里失效。nutkWallFunction在对数律区外推出来的壁面剪应力会系统性偏大。解决二选一不要骑墙。走壁面函数就把firstLayerThickness加到 y 50 左右层数保持 5 层以上、膨胀比 1.2走全解析就把首层压到 y 1同时确认nut用nutLowReWallFunction或对应的低雷诺数处理。改完一定重新跑一遍yPlus别看设置值。5.2 阻塞比超标迎风面 Cp 超过 1.0现象驻点附近 Cp 算到 1.11.3理论上限是 1.0怎么调湍流模型都下不来。原因计算域截面太小气流绕流时被挤着加速参考静压和来流速度都对不上。这是纯几何问题跟数值格式无关。解决按 3.1 节的脚本重算域宽把阻塞比压到 5% 以内对标风洞就压到 3%。加宽之后网格总数会涨可以同步把远处网格放粗总网格量不一定增加多少。改完先跑一个空域不放建筑验证入口给剖面出口截面的平均速度应该和入口基本一致偏差超过 3% 说明域还不够。5.3 稳态算钝体不收敛残差在 1e-3 附近画正弦现象simpleFoam前 500 步残差降得挺顺之后开始周期性振荡监测量某个测点的 p也跟着抖再跑 2000 步还是这样。原因钝体绕流本身有涡脱落物理上就是非定常的。强行求稳态解等于让迭代在几个近似解之间跳。不是数值没调好是物理假设不成立。解决三个方向。一是接受它用fieldAverage对最后一段做时间平均得到准平均场工程上够用二是换pimpleFoam做非定常 RANS时间步按涡脱落周期的 1/100 取三是降松弛因子、开伪瞬态pseudoTransient或localEuler把振荡压下去。我一般先试第三个压不住再上第二个。5.4 入口只给均匀风整体 Cp 偏差不大分离区全错现象迎风面 Cp 和风洞对得上顶部和侧面分离区的吸力明显偏弱尾流长度偏短。原因均匀入口没有速度梯度靠的是模型自身的钝体分离产生剪切缺少大气边界层自带的那份湍动能。同时入口 k 和 epsilon 如果给了个常数量级也往往不对k 给小了。解决按 4.1 和 4.2 补上剖面和湍流入口。补完之后一定要验证在模型位置放一个空的竖直采样线不放建筑看速度剖面是否保持你给的形状对数律剖面在下游衰减超过 5% 就说明边界条件不守恒得换atmBoundaryLayer系列或者加一段前置发展段。5.5 snappyHexMesh 层塌陷表面出现一片秃斑现象模型表面局部没有边界层网格checkMesh报collapsed或者层数明显少于设定值非正交角在那一带飙到 80 度以上。原因minThickness本文还有配套的精品资源点击获取
简介这份文献是面向风工程与计算流体力学初学者的OpenFOAM基础培训讲义适合刚接触开源CFD、希望建立软件整体认识的读者。压缩包仅含1个PDF文件约162KB篇幅轻量便于快速浏览与课堂配套学习。已有555人学习。内容从OpenFOAM的开发背景与版本演进讲起梳理其面向对象C类库、有限体积离散、非结构化多面体网格、大规模并行计算等核心能力并延伸到自动网格运动、网格拓扑变化、拉格朗日粒子追踪、有限元与有限面积离散等扩展功能。讲稿还介绍元类、容器类、字典、实时选择表与分层设计理念可帮助读者理解求解器框架、代码重用和用户定制思路为后续阅读源码、搭建风工程算例和排查网格与并行设置问题提供入门线索。1. 风工程为什么用 OpenFOAM一次阻塞比返工接手一个屋面风压复核任务时模型是别人拷过来的网格看着很规整simpleFoam收敛得也漂亮。结果迎风面中部 Cp 算到 1.3屋角吸力峰值只有 -1.1跟风洞报告对不上。查了两天问题出在计算域宽度只留了 3 倍建筑高度阻塞比接近 12%气流被硬挤过去整个压力场都是歪的。这类返工在风工程里太常见了OpenFOAM 本身不难难的是把风洞里的物理约束翻译成命令行里的参数。风工程 OpenFOAM 基础培训要解决的正是这件事把 OpenFOAM 当作一台计算风洞来用而不是当通用 CFD 软件来用。它适合结构风荷载、建筑群风环境、行人高度风舒适度这类以钝体绕流为核心的问题也适合已经会用商业软件、但说不清湍流入口和 y 到底该怎么定的人。2. 求解器与湍流模型怎么挑从 simpleFoam 到 LES 的选型边界2.1 不可压求解器分工simpleFoam、pisoFoam、pimpleFoam风工程里绝大多数算例的风速远低于 0.3 马赫空气按不可压处理密度当成常数。真正要选的是时间项怎么处理。simpleFoam走 SIMPLE 算法把时间项设成steadyState一轮迭代直接求解平均流场成本最低方案比选、多风向角扫掠、参数敏感性分析都靠它。它的局限也很明确稳态假设下没有涡脱落钝体尾流里那些周期性的东西会被抹平成一个看起来还行的定常解。pisoFoam走 PISO 算法真正的时间推进配合 LES 或 DES 用。它要求 Courant 数压在 1 以下时间步长往往在 1e-4 秒量级跑十几秒物理时间就是几十万步成本是稳态的几十到上百倍。值不值得取决于你要的是平均风压还是峰值吸力和脉动风压谱——后者只能靠它。pimpleFoam走 PIMPLEPISO 和 SIMPLE 的混合允许 Courant 数大于 1可以只做一次外迭代校正。非定常 RANS、带热浮力的街区通风、以及稳态不收敛想试试非定常的过渡场景我一般先拿它试水。求解器时间项风工程典型任务相对成本simpleFoamsteadyState平均风压、平均风速场、方案比选1pimpleFoam瞬态可大 Courant非定常 RANS、热浮力、过渡态1030pisoFoam瞬态Courant 1LES/DES涡脱落、脉动风压50200potentialFoam无给稳态求解器一个像样的初场 1较新的发行版把求解器合并成了统一的foamRun -solver incompressibleFluid功能一致只是入口不同。选型逻辑不变先问自己要平均量还是瞬时量。2.2 湍流模型k-ε 稳、k-ω SST 贴壁、LES 抓涡标准 k-ε 是最省事的选择。它对自由剪切层和远场比较宽容收敛性好工程上算整体风压分布、街区平均风速场够用。缺点写在教科书里也写在实践里钝体背风侧的分离区会被过度预测再附点往后拖背压偏低。某跨平台系统做过方柱绕流的对比标准 k-ε 给出的再附长度比实验长约 20%realizable k-ε 能收到 10% 以内。k-ω SST 在近壁区用 k-ω远场切回 k-ε分离预测比 k-ε 准是钝体绕流里性价比最高的一档。代价是对网格更挑想要它发挥优势壁面第一层要压到 y ≈ 1如果加壁面函数则要保持 y 在 30100最怕停在 530 的缓冲层里两头不靠。LES 是另一档。它能直接解析大尺度涡给出峰值吸力和脉动风压谱屋面围护结构设计、行人高度阵风这些都靠它。但网格量按 Re^1.8 左右的趋势涨一栋 60 米高的楼按 1:100 缩尺做 LES网格轻松过千万还得跑够统计时间。我的习惯是先稳态拿平均场定方案再挑一两个关键风向做 LES 精细核算不全程上 LES。提示湍流模型选了 k-ε 就别指望背风面 Cp 精确到小数点后两位用它做相对比较A 方案和 B 方案谁更差比做绝对值标定靠谱得多。2.3 最小算例骨架0/ constant/ system/ 三个目录OpenFOAM 的算例结构是死的风工程算例也一样。摸清这三个目录剩下的都是填内容。# 从自带算例拷一个不可压模板改成自己的 cp -r $FOAM_TUTORIALS/incompressible/simpleFoam/pitzDaily ./caseWind cd caseWind # 目录骨架 # 0/ U p k epsilon nut —— 初始场与边界条件 # constant/ transportProperties —— 运动黏度 nu # turbulenceProperties —— 湍流模型开关 # polyMesh/ —— blockMesh 产出的背景网格 # triSurface/ —— 建筑 STL # system/ blockMeshDict snappyHexMeshDict # controlDict fvSchemes fvSolutiontransportProperties里只关心nu20 摄氏度空气取 1.5e-5 m²/s缩尺模型不改 nu改的是特征尺寸和风速。turbulenceProperties里simulationType RAS配RASProperties选模型LES则换LESProperties。fvSchemes和fvSolution是数值层面的东西稳态钝体绕流有一套相对固定的写法下面会贴。真正决定结果对不对的是0/目录边界条件写错后面调什么都没用。3. 计算域与网格把阻塞比和 y 先在纸上算清楚3.1 域尺寸的硬约束上游 5H、下游 15H 与 3% 阻塞比域尺寸不是凭感觉画的。上游一般留 5HH 为建筑高度让来流有空间发展下游留 15H尾流有足够长度恢复否则出口处还有回流pressureInletOutletVelocity会被逼着往计算域里灌气残差就开始抖侧面单侧 5H、顶部 56H保证顶部不出现明显的加速。阻塞比是更硬的约束。风洞实验的通行门槛是迎风投影面积不超过试验段截面的 3%数值模拟可以放宽到 5%但超过 5% 之后迎风面 Cp 会系统性偏高屋角吸力偏低误差可以到 20% 以上。# 由建筑尺寸反算最小计算域宽度保证阻塞比不超标 H, B, D 60.0, 40.0, 30.0 # 建筑高、迎风宽、进深单位 m blockage 0.03 # 风洞 3%纯数值可放宽到 0.05 A_front B * H # 迎风面投影面积 A_cross A_front / blockage # 所需最小横截面积 H_domain 6.0 * H # 域高取 6H W_domain A_cross / H_domain # 由面积反算域宽 print(f迎风面积 {A_front:.0f} m2截面至少 {A_cross:.0f} m2) print(f域高 {H_domain:.0f} m 时总宽至少 {W_domain:.0f} m单侧 {W_domain/2:.0f} m) print(f建议域范围 x:[-{5*H:.0f}, {15*H:.0f}] y:±{W_domain/2:.0f} z:[0, {H_domain:.0f}])代入 60 米高、40 米宽的建筑3% 要求总宽至少 222 米单侧 111 米。这个数字很多人第一眼觉得夸张但它就是风洞准则换算过来的结果。实际做的时候我会先按 5% 估一版跑完看迎风面 Cp 是否明显高于 1.0驻点理论上限偏高就加宽。注意域加宽不是免费的网格量按面积涨。宽域配粗网格和窄域配细网格前者通常更准。3.2 blockMesh 背景网格关键行与 checkMesh 判据blockMesh只负责铺一个六面体背景网格建筑由snappyHexMesh挖出来。所以背景网格不用太细均匀一点、正交性好一点就行真正的加密交给后面。// system/blockMeshDict 节选H 60 m 的算例 vertices ( (-300 -150 0) // 0 ( 900 -150 0) // 1 ( 900 150 0) // 2 (-300 150 0) // 3 (-300 -150 360) // 4 ( 900 -150 360) // 5 ( 900 150 360) // 6 (-300 150 360) // 7 ); blocks ( hex (0 1 2 3 4 5 6 7) (120 60 60) simpleGrading (1 1 3) // z 向渐变地面附近密、高空疏 );simpleGrading (1 1 3)的意思是 x、y 方向均匀z 方向最后一层与第一层厚度比为 3。大气边界层的速度梯度集中在地面附近这样铺能在不增加网格量的前提下把地面分辨率提上去。三个方向都用单一渐变会限制加密位置真要做建筑附近密、远处疏得多块拼或者干脆靠snappyHexMesh的refinementRegions补。boundary段里定义inlet、outlet、ground、top、side五个 patch类型先都写patch具体边界条件在0/里再定。生成完立刻跑checkMeshblockMesh log.blockMesh 21 checkMesh log.checkMesh 21 grep -E non-orthogonality|skewness|aspect ratio log.checkMesh三个数字要盯最大非正交角最好小于 65 度平均小于 40 度最大偏斜小于 4长宽比在边界层里大是正常的但背景网格层面别超过 1000。这三个指标失控后面snappyHexMesh贴出来的层也会跟着歪。3.3 snappyHexMesh 贴体加密三层加密与首层厚度snappyHexMesh是风工程网格的主力。流程是surfaceFeatureExtract提特征线snappyHexMesh挖洞、加密、贴层最后checkMesh复核。// system/snappyHexMeshDict 关键段 castellatedMesh true; snap true; addLayers true; geometry { building.stl { type triSurfaceMesh; name building; } } castellatedMeshControls { maxLocalCells 4000000; maxGlobalCells 20000000; minRefinementCells 10; nCellsBetweenLevels 3; // 相邻加密级之间至少隔 3 层防跳变 resolveFeatureAngle 30; features ( { file building.eMesh; level 3; } ); refinementSurfaces { building { level (2 3); } // 最小 2 级、最大 3 级 } locationInMesh (0 0 200); // 取一个肯定在流体域内的点 allowFreeStandingZoneFaces true; } addLayersControls { layers { building { nSurfaceLayers 5; } ground { nSurfaceLayers 3; } } expansionRatio 1.2; // 相邻层厚度比别超过 1.3 firstLayerThickness 0.002; // 首层 2 mm由 y 反算得来 minThickness 0.0004; // 最小层厚低于它就塌陷 featureAngle 60; nGrow 0; }level (2 3)表示最小加密 2 级、特征线附近 3 级背景网格 5 米的话最细能到 0.6 米。nCellsBetweenLevels 3是个容易被忽略但很关键的数加密级之间跳得太快非正交角会飙snap阶段就容易出负体积。firstLayerThickness是绝对厚度单位米必须由 y 反算不能拍脑袋。minThickness一般取首层的 0.10.2 倍太小了不起作用太大会直接导致层塌陷。surfaceFeatureExtract log.surfaceFeatureExtract 21 snappyHexMesh -overwrite log.snappyHexMesh 21 checkMesh log.checkMesh 21 grep -i collapsed\|negative log.snappyHexMesh-overwrite省掉一次改名但出错就回不去了第一次跑建议不加。日志里出现collapsed或negative volume先看minThickness和nSurfaceLayers绝大多数是这两个数给大了。3.4 第一层高度反算一张 y 对照表和两个公式y 的定义是 y uτ·y/ν其中 uτ 是壁面摩擦速度y 是首层高度ν 是运动黏度。工程上 uτ 可以用 0.05·Uref 粗估够用了。按 ν 1.5e-5 m²/s 算Uref (m/s)uτ 估计 (m/s)y1 首层y30 首层y100 首层50.250.06 mm1.8 mm6 mm100.500.03 mm0.9 mm3 mm150.750.02 mm0.6 mm2 mm201.000.015 mm0.45 mm1.5 mm表里最扎眼的是一列比一列小。想全解析壁面y ≈ 110 m/s 风速下首层只要 0.03 毫米一栋 60 米的楼配上这个首层边界层网格量能顶到天上去。这就是为什么工程 RANS 普遍走壁面函数路线把首层放在 y 30100取个 13 毫米层数 5 层、膨胀比 1.2。# 由目标 y 反算首层高度反过来也行 nu, Uref, yplus 1.5e-5, 10.0, 50.0 u_tau 0.05 * Uref y_first yplus * nu / u_tau print(fu_tau{u_tau:.3f} m/s, y{yplus} 对应首层 {y_first*1000:.2f} mm) # 已知首层验算实际 y y_given 0.002 print(f首层 {y_given*1000:.1f} mm 实际 y {u_tau*y_given/nu:.0f})算完之后别急着信跑起来用yPlusfunctionObject 输出一遍实际分布。模型表面的 y 是随位置变的迎风面驻点附近低、屋角高能差三倍。看的是绝大部分面积落在哪个区间不是平均值。4. 大气边界层入口怎么写速度剖面、粗糙度与出口配置4.1 幂律剖面的 codedFixedValue 写法风洞里的来流不是均匀风是带速度梯度的大气边界层。直接给fixedValue uniform (10 0 0)是最常见的错误驻点压力勉强对顶部加速区和分离剪切层全错因为整个流场没有剪切的驱动力。幂律是最简单的一种写法用codedFixedValue直接按面心坐标算inlet { type codedFixedValue; value uniform (0 0 0); name powerLawInlet; code #{ const scalar Uref 8.0; // 参考高度风速 m/s const scalar Zref 10.0; // 参考高度 m const scalar alpha 0.16; // 地貌指数 const vector dir(1, 0, 0); // 来流方向 const fvPatch p patch(); const vectorField Cf p.Cf(); // 每个面心的坐标 vectorField U *this; // 当前边界上的 U forAll(Cf, i) { const scalar z max(Cf[i].z(), 0.0); U[i] dir * Uref * pow(z / Zref, alpha); } #}; }三个参数里alpha决定剖面形状开阔水面 0.12开阔草地 0.140.16郊区低层建筑 0.20密集城区 0.250.30 都有人用。Uref和Zref是配套的别把 10 米参考风速配 1 米参考高度那算出来的剖面整体偏大一个量级。缩尺模型要按缩尺比换算几何按米建模、风速也按对应的速度比换算nu保持真实值不变。提示codedFixedValue每次求值都要重跑一遍循环边界面上万面时开销可观。面数特别大的话改用timeVaryingMappedFixedValue读一张预先写好的剖面表更省。4.2 k 与 epsilon 入口湍流强度 10% 是怎么来的只给速度剖面湍流是零边界层是层流化的分离剪切层起不来模型的侧面和背风面会明显失真。入口必须同时给 k 和 epsilon或者 omega。k 由湍流强度反算k 1.5·(U·I)²。风洞大气边界层的湍流强度在参考高度附近大致 8%15%取 10% 是常见起点。epsilon 由混合长度反算混合长度取边界层厚度的 0.070.1 倍。// 0/k inlet { type turbulentIntensityKineticEnergyInlet; intensity 0.10; // 10% 湍流强度 value uniform 1.2; // 8 m/s、10% 时约 1.2 m2/s2 } // 0/epsilon inlet { type turbulentMixingLengthDissipationRateInlet; mixingLength 20.0; // 约 0.1 倍边界层厚度 m value uniform 0.05; }intensity和value的关系是 k 1.5(U·I)²求解器会用value作为初值、按intensity随时间更新。两个数不一致会导致第一步的 k 场有个跳变虽然能收敛但日志难看。用 k-ε 有个绕不开的问题幂律速度剖面配常数 k、epsilon在流场里不是平衡态往下游走剖面会自己衰减模型位置上根本到不了目标风速。较新的发行版内置了atmBoundaryLayer系列边界条件按对数律同时给 U、k、epsilon能保持平衡。字段名各分支略有差异落到自己环境里先用源码或foamSearch确认一遍再写。风洞里测了剖面的老老实实用timeVaryingMappedFixedValue读表那是最保真的做法。4.3 地面粗糙度roughnessLength 与壁面函数的配对地面是 wall但不能当光滑壁面处理。城市下垫面的粗糙度直接影响底部风速和街谷里的流动。做法是给地面的nut用粗糙壁面函数把空气动力学粗糙度长度 z0 换算成等效砂粒粗糙度 Ks// 0/nut ground { type nutkRoughWallFunction; Ks uniform 1.0; // 等效砂粒粗糙度由 z0 换算 Cs uniform 0.5; // 粗糙度常数 value uniform 0; }换算式是 Ks ≈ 9.793·z0/Cs。z0 取 0.05 米郊区/低层城区Cs 取 0.5得 Ks ≈ 1.0 米。这个值看着大但它不是几何尺寸是一个等效参数别拿去跟首层高度比。配套的还有0/k地面用kqRWallFunctionvalue 设 0 即可0/epsilon地面用epsilonWallFunction。k 和 epsilon 的壁面处理没配对最常见的表现是残差降不下去、地面附近速度出现非物理的加速。如果是 LES地面就没这么简单了得用nutUSpaldingWallFunction或者干脆把首层压到 y ≈ 1 全解析。后者成本高通常只在关键区域模型周边这么做。4.4 出口、顶部、侧面的边界类型与回流风险出口的经典组合是速度用pressureInletOutletVelocity压力用fixedValue uniform 0。前者在出流时是零梯度一旦局部出现回流就切回固定值比单纯的zeroGradient稳得多是尾流算例的标配。顶面和侧面用slip速度压力zeroGradient。用slip而不是 wall是为了模拟风洞试验段的自由滑移边界避免侧壁面额外生成边界层把流场往中间挤。有人喜欢用symmetry在纯速度场上和slip等价看习惯。// 0/U 的出口段 outlet { type pressureInletOutletVelocity; value uniform (0 0 0); } // 0/p 的出口段 outlet { type fixedValue; value uniform 0; }出口离模型太近时回流会一直顶到模型背风面表现为背压忽高忽低、残差周期性抖动。这时候要么把出口往下游推要么在出口前加一段缓变的多孔区工程上少用不推荐最简单的还是回到 3.1 节把下游长度补到 15H。5. 避坑排查五个让算例翻车的现场5.1 y 停在 530 缓冲层壁面函数和解析壁面两头不靠现象模型表面 Cp 分布形状对但整体偏移迎风面偏高、屋角吸力偏低yPlus输出显示表面大部分区域落在 530。原因这一层既不能当对数律区用壁面函数也没细到能解析黏性底层两种壁面处理都在这里失效。nutkWallFunction在对数律区外推出来的壁面剪应力会系统性偏大。解决二选一不要骑墙。走壁面函数就把firstLayerThickness加到 y 50 左右层数保持 5 层以上、膨胀比 1.2走全解析就把首层压到 y 1同时确认nut用nutLowReWallFunction或对应的低雷诺数处理。改完一定重新跑一遍yPlus别看设置值。5.2 阻塞比超标迎风面 Cp 超过 1.0现象驻点附近 Cp 算到 1.11.3理论上限是 1.0怎么调湍流模型都下不来。原因计算域截面太小气流绕流时被挤着加速参考静压和来流速度都对不上。这是纯几何问题跟数值格式无关。解决按 3.1 节的脚本重算域宽把阻塞比压到 5% 以内对标风洞就压到 3%。加宽之后网格总数会涨可以同步把远处网格放粗总网格量不一定增加多少。改完先跑一个空域不放建筑验证入口给剖面出口截面的平均速度应该和入口基本一致偏差超过 3% 说明域还不够。5.3 稳态算钝体不收敛残差在 1e-3 附近画正弦现象simpleFoam前 500 步残差降得挺顺之后开始周期性振荡监测量某个测点的 p也跟着抖再跑 2000 步还是这样。原因钝体绕流本身有涡脱落物理上就是非定常的。强行求稳态解等于让迭代在几个近似解之间跳。不是数值没调好是物理假设不成立。解决三个方向。一是接受它用fieldAverage对最后一段做时间平均得到准平均场工程上够用二是换pimpleFoam做非定常 RANS时间步按涡脱落周期的 1/100 取三是降松弛因子、开伪瞬态pseudoTransient或localEuler把振荡压下去。我一般先试第三个压不住再上第二个。5.4 入口只给均匀风整体 Cp 偏差不大分离区全错现象迎风面 Cp 和风洞对得上顶部和侧面分离区的吸力明显偏弱尾流长度偏短。原因均匀入口没有速度梯度靠的是模型自身的钝体分离产生剪切缺少大气边界层自带的那份湍动能。同时入口 k 和 epsilon 如果给了个常数量级也往往不对k 给小了。解决按 4.1 和 4.2 补上剖面和湍流入口。补完之后一定要验证在模型位置放一个空的竖直采样线不放建筑看速度剖面是否保持你给的形状对数律剖面在下游衰减超过 5% 就说明边界条件不守恒得换atmBoundaryLayer系列或者加一段前置发展段。5.5 snappyHexMesh 层塌陷表面出现一片秃斑现象模型表面局部没有边界层网格checkMesh报collapsed或者层数明显少于设定值非正交角在那一带飙到 80 度以上。原因minThickness本文还有配套的精品资源点击获取
