OpenFOAM二次开发教程07第一个二次开发实战——自定义标量输运求解器版本与事实声明标量输运方程的 OpenFOAM 写法fvm::ddt/fvm::div/fvm::laplacian/fvm::Sp/fvm::Su以本机$FOAM_SRC/finiteVolume/finiteVolume/fvm/下的算子实现与官方 Doxygen 为准。fvScalarMatrix是fvMatrixscalar的常用别名形式求解接口见官方 DoxygenfvMatrix.HSolverPerformanceType solve(...)。源项的物理含义与系数均为教学示例值不代表任何标准规定也不对应任何真实装置。验证方法残差轨迹、守恒性检查为通用工程做法具体判据阈值需按问题尺度自定。一句话结论一个实用的标量输运求解器只需写出fvm::ddt(T) fvm::div(phi, T) - fvm::laplacian(DT, T) fvm::Sp(S, T) Su这一条方程并solve()把耗散型源项用fvm::Sp隐式化是让求解器在强源项下仍稳定的关键一步。〇、本篇要解决的认知问题Q1标量输运方程在 OpenFOAM 里怎么写每一项用哪个算子、隐式还是显式Q2源项怎么处理才既有物理意义又稳定fvm::Sp与Su如何配合Q3加了新场T之后除了代码还必须在算例里补哪些配置漏了会报什么错Q4怎么判断我算的标量场是对的守恒性怎么检查Q5这个求解器可以直接改成哪些实用模型温度场、组分浓度、被动标量、湍流标量一、机制解析1.1 标量输运方程为什么它是二次开发的Hello World标量输运方程是 CFD 里最通用的形式之一∂T/∂t ∇·(U T) - ∇·(D ∇T) S │ │ │ │ 时间项 对流项 扩散项 源项它之所以是二次开发的第一个真正有用的例子是因为温度场、组分浓度、盐度、被动标量、湍动能等大量物理量都写成这个形式。学会它你就掌握了 80% 的加一个新方程的场景。为什么这对你重要很多工程需求本质就是给现有流动加一个标量场——比如在已有流场上算温度分布“算示踪剂浓度”“算某种添加剂的扩散”。这些都不需要重写流场求解器只需要写这个标量方程并耦合进现有求解流程。1.2 每一项的算子选择与隐式/显式决策方程项OpenFOAM 写法隐式/显式决策理由时间项fvm::ddt(T)隐式时间推进必须隐式否则时间步长受限于 CFL对流项fvm::div(phi, T)隐式对流是输运主角隐式处理保证稳定性phi视为已知面通量扩散项fvm::laplacian(DT, T)隐式扩散系数DT视为已知若DT依赖T则需外迭代源项fvm::Sp(S, T) Su半隐式负系数部分隐式增强对角占优常数部分显式反直觉的默认值陷阱phi面通量在不可压缩求解器里不是从T推出来的而是在流场求解阶段就已经算好、且必须满足连续性∇·phi 0。如果phi不守恒不满足连续性标量场必然不守恒——这是温度场总量漂移最常见也最隐蔽的根因排查时要先查流场而不是查标量方程。1.3 源项的半隐式处理从Sp/Su到物理直觉把源项写成关于 T 的线性形式是 CFD 的标准做法S(T) ≈ Sp · T Su Sp 与 Su 本身与 T 无关可以是场其中Sp通常取负耗散Su是常数部分生成。放回方程fvScalarMatrixTEqn(fvm::ddt(T)fvm::div(phi,T)-fvm::laplacian(DT,T)fvm::Sp(Sp,T)Su// Sp·T 隐式进矩阵Su 显式进右端);为什么Sp取负就稳定一个负的Sp相当于增加矩阵对角而矩阵求解的迭代收敛条件正是对角占优。所以任何随 T 增大而增大消耗的物理都应该写成负Sp。经验法则源项拆分的判据是能否写出线性形式。如果源项强烈非线性如 Arrhenius 反应速率仍可用负系数部分隐式化 其余显式配合外迭代逐步更新——这正是半隐式能稳定的原因。1.4 加一个新场必须同步改三处配置这是新手最常漏的一环。新增T场后必须在算例里补三处位置补什么漏了的后果0/T场文件含dimensions、internalField、每个 patch 的typeMUST_READ报 “cannot find file”system/fvSolution的solversT的线性求解器条目运行期报找不到 T 的 solver 设置system/fvSchemes若T用到新项如div(phi,T)、laplacian(DT,T)的格式运行期报 “keyword … undefined”若用了default noneconstant/transportProperties视实现扩散系数等物性运行期报字典缺项记住口诀加场 加文件 加 solver 加可能格式。三处齐了才能跑。1.5 怎么判断算对了三个层次的验证数值层残差是否单调下降到容差以下SolverPerformance判据。守恒层无源、封闭边界时标量的总量应守恒体积分近似不变。这是最有力的物理判据。极限层把扩散系数设为极大 → 场应趋于均匀把速度设为零 → 只剩扩散可对照解析或已知行为。铁律 7 的落地任何自定义求解器都必须通过这三个层次中的至少两个才能用于工程。跑通了和算对了是两回事。二、完整代码与逐行剖析代码 2-1myScalarTransportFoam.C/*---------------------------------------------------------------------------*\ myScalarTransportFoam.C —— 带对流、扩散与半隐式源项的标量输运求解器 方程∂T/∂t div(phi,T) - laplacian(DT,T) Sp*T Su 依赖需要一个已知的不可压缩流场U, phi作为输入示例假设其已存在/已解出 \*---------------------------------------------------------------------------*/#includefvCFD.Hintmain(intargc,char*argv[]){#includesetRootCase.H#includecreateTime.H#includecreateMesh.H#includecreateFields.H// 本求解器的场U、phi、T、DT、Sp、Su#includeinitContinuityErrs.HInfo\nStarting scalar transport time loop\nendl;while(runTime.loop()){InfoTime runTime.timeName()nlendl;// ---- 若流场需同步求解应在此处先解 U本篇假设 phi 已知且满足连续性----// ---- 装配标量输运方程 ----fvScalarMatrixTEqn(fvm::ddt(T)fvm::div(phi,T)// 对流phi 为面通量面心场-fvm::laplacian(DT,T)// 扩散DT 为扩散系数此处设为常量场fvm::Sp(Sp,T)// 半隐式源项耗散部分进矩阵对角Su// 显式源项常数/生成部分);// ---- 求解并打印收敛统计 ----SolverPerformancescalarperfTEqn.solve();InfoT: nIter perf.nIterations() initial perf.initialResidual() final perf.finalResidual() converged perf.converged()nlendl;// ---- 边界同步内部场已更新边界必须重新求值 ----T.correctBoundaryConditions();// ---- 守恒性监控统计 T 的体积分用于守恒性判据----// T*mesh.V() 给出每个单元的“T x 体积”sum 得到总量constscalar totalTgSum(T.primitiveField()*mesh.V());Infointegral(T) totalTnlendl;runTime.write();InfoExecutionTime runTime.elapsedCpuTime() s\nendl;}InfoEnd\nendl;return0;}逐行剖析方程的四项一一对应 §一.2 的表格fvm::Sp(Sp, T) Su是半隐式源项的标准写法。注意Sp与Su都可以是场volScalarField这让空间上非均匀的源变得容易表达。SolverPerformancescalar perf TEqn.solve();与第 06 篇一致打印残差是纪律不是可选项。T.correctBoundaryConditions()再次强调解完方程后边界必须同步。gSum(T.primitiveField() * mesh.V())守恒性监控的核武器。gSum是 OpenFOAM 的并行归约求和跨所有 processor 求和mesh.V()是单元体积。把每一时刻的integral(T)打印出来就能看出总量是否漂移——若漂移明显且边界无源说明phi不守恒或边界设置不当。perf.converged()不要忽略这个布尔量。当它长期为假说明容差没达到结果不可信。代码 2-2createFields.H新增场与系数/*---------------------------------------------------------------------------*\ createFields.H —— myScalarTransportFoam 的场集合 \*---------------------------------------------------------------------------*/// ---- 流场本求解器假设 phi 已知由前置流场求解或初值提供----InfoReading field U\nendl;volVectorFieldU(IOobject(U,runTime.timeName(),mesh,IOobject::MUST_READ,IOobject::AUTO_WRITE),mesh);#includecreatePhi.H// 由 U 与网格构造面通量 phi面心场// ---- 待求标量场 T ----InfoReading field T\nendl;volScalarFieldT(IOobject(T,runTime.timeName(),mesh,IOobject::MUST_READ,IOobject::AUTO_WRITE),mesh);// ---- 扩散系数 DT演示用常量场真实问题可从 transportProperties 读取 ----InfoReading diffusivity DT\nendl;volScalarFieldDT(IOobject(DT,runTime.timeName(),mesh,IOobject::NO_READ,IOobject::NO_WRITE),mesh,dimensionedScalar(DT,dimViscosity,1e-5)// 示例值请按真实物性与单位设定);// ---- 源项系数 Sp隐式通常取负与 Su显式----// 这里用“常量场”演示真实问题可由 cellZone、温度依赖关系或外部数据驱动。volScalarFieldSp(IOobject(Sp,runTime.timeName(),mesh,IOobject::NO_READ,IOobject::NO_WRITE),mesh,dimensionedScalar(Sp,dimless/dimTime,-0.1)// 示例值负系数保证对角占优);volScalarFieldSu(IOobject(Su,runTime.timeName(),mesh,IOobject::NO_READ,IOobject::NO_WRITE),mesh,dimensionedScalar(Su,dimless/dimTime,0.0)// 示例值默认无生成项);逐行剖析U用MUST_READ本求解器把流场当输入所以必须有0/U。若流场也要在本求解器里解需要在时间循环里加上速度-压力耦合第 20 篇的耦合插件会涉及。DT/Sp/Su用NO_READ 构造函数给初值这是程序内部生成的辅助场的标准做法——它们不来自磁盘而是由代码定义。NO_WRITE表示不必写回磁盘除非你要后处理它们。维度用dimensionedScalar指定dimViscosity运动黏度量纲用于DTSp/Su用无量纲/时间即 1/s是因为方程里它们与∂T/∂t同级。这些维度必须与方程各项自洽否则装配时报维度不匹配第 06 篇报错 3-3。反直觉点DT设为1e-5只是示例。真实扩散系数必须从constant/transportProperties或热物性字典读取第 08 篇写死在代码里是工程禁忌——它让物理参数无法通过配置调整。代码 2-3算例配置补丁system/fvSchemes与system/fvSolution需新增的部分// ---------- system/fvSchemes 需确保包含 ---------- ddtSchemes { default Euler; // 时间项一阶隐式欧拉可换 backward 提高精度 } divSchemes { default none; div(phi,T) Gauss limitedLinear 1; // 标量对流限幅线性有界性较好 } laplacianSchemes { default Gauss linear corrected;// 扩散含非正交修正 } // ---------- system/fvSolution 需新增 T 的求解器 ---------- solvers { T { solver PBiCGStab; // 非对称对流项使矩阵一般非对称 preconditioner DILU; tolerance 1e-8; relTol 0.01; } (U|k|epsilon) // 若同时解流场/湍流可用正则一次配置多个 { solver smoothSolver; smoother symGaussSeidel; tolerance 1e-6; relTol 0.1; } }逐行剖析div(phi,T)用限幅格式如limitedLinear标量对流最容易出现非物理过冲/欠冲限幅格式能保证有界性是工程默认选择。经验法则标量输运优先用有界格式除非你确知需要高精度且网格足够好。PBiCGStabDILU非对称矩阵的常规组合。为什么标量方程通常非对称因为对流项破坏了对称性。选错成PCG对称求解器不是不能用而是收敛慢、可能不收敛。(U|k|epsilon)这种正则式子字典OpenFOAM 支持用正则表达式一次为多个场配置求解器这是大型算例里省事的标准写法正则匹配的场名必须真的存在。注意这里没有给DT/Sp/Su配 solver它们是NO_READ的辅助场不作为方程求解对象。只有被solve()的场才需要 solver 条目——这是判断该不该加 solver的判据。代码 2-4守恒性与残差验证脚本#!/bin/sh# verify_scalar.sh —— 运行并抽取“残差 T 总量”轨迹做守恒性检查# 用法CASE_SRC含 0/U 与 0/T 的算例 sh verify_scalar.shset-euCASE_SRC${CASE_SRC:?请用 CASE_SRC... 指定一个已建网格、含 0/U 与 0/T 的算例}WORK$PWD/_scalar_caserm-rf$WORK;cp-r$CASE_SRC$WORKcd$WORKmyScalarTransportFoam-case.log.myScalarTransportFoam21||{echo[FAIL] 求解器退出非零见 log.myScalarTransportFoam;exit1;}echo 残差与总量轨迹 grep-Einitial |integral\(T\) log.myScalarTransportFoamecho 判据 1残差是否收敛 ifgrep-qconverged 1log.myScalarTransportFoam;thenecho[OK] 至少一步报告 converged 1elseecho[WARN] 未见 converged 1请检查 fvSolution 容差与 fvSchemes 格式fiecho 判据 2总量漂移相对变化# 取第一个与最后一个 integral(T) 值算相对变化first$(grepintegral(T) log.myScalarTransportFoam|head-1|awk{print $NF})last$(grepintegral(T) log.myScalarTransportFoam|tail-1|awk{print $NF})awk-va$first-vb$lastBEGIN{ if (a 0) { print [INFO] 初始总量为 0无法计算相对漂移; exit } d (b - a) / a; if (d 0) d -d; printf 初始%.6g 末态%.6g 相对漂移%.4f%%\n, a, b, d*100; if (d 0.01) print [OK] 相对漂移 1%无源且封闭边界下的期望行为; else print [WARN] 相对漂移 1%优先排查 phi 是否满足连续性div(phi)≈0; }逐行剖析grep -q converged 1把是否收敛过作为硬判据之一直接来自SolverPerformance。用awk算相对漂移而非拍脑袋守恒性判据必须量化第 04 篇练习规范。漂移大时的第一排查方向phi不守恒div(phi)不为零。这是最容易被忽略、但最致命的根因——标量方程的守恒性继承自流场的连续性。依然在官方算例副本上跑铁律 7。三、常见报错与排查报错 3-1-- FOAM FATAL IO ERROR: cannot find file .../0/T。现象启动即报找不到 T。根因新增了场但没在0/下加文件§一.4 的三处配置只做了代码。解法从同类官方算例复制一个T文件修改dimensions、internalField与各 patch 的type特别检查dimensions是否与代码里T的预期维度一致。报错 3-2-- FOAM FATAL IO ERROR: keyword solver ... T is undefined找不到 T 的求解器。现象进入求解时报缺 solver 配置。根因fvSolution的solvers块里没有T条目。解法加入T子字典配solver/preconditioner/tolerance/relTol。判据凡是被solve()的场都必须有 solver 条目NO_READ的辅助场不需要。报错 3-3-- FOAM FATAL IO ERROR: keyword div(phi,T) is undefined in fvSchemes。现象装配对流项时找不到格式。根因divSchemes里用了default none且没有为div(phi,T)指定格式注意括号内表达式必须与代码一致。解法在divSchemes中加div(phi,T) 格式;。注意场名大小写与括号内写法要和代码完全一致。报错 3-4T 场出现负值或过冲温度跑到物理上不可能的区间。现象结果里T出现非物理极值。根因对流格式无界如用了Gauss linear纯中心差分或源项显式处理导致过冲。解法对流改用有界格式limitedLinear、vanLeer、limitedLinearV等具体名称以官方文档为准把耗散型源项改为fvm::Sp隐式必要时减小时间步deltaT。报错 3-5-- FOAM FATAL ERROR: dimension mismatch ...。现象装配方程时维度冲突。根因DT、Sp、Su的维度与方程其他项不匹配例如把Sp配成了无量纲而非 1/时间。解法逐一打印各量维度对账让Sp/Su的维度与∂T/∂t同级即量纲/时间DT的量纲使laplacian(DT,T)与∂T/∂t同级。这是纯粹的维度代数画在纸上两分钟就能推清楚比试错快得多。四、动手练习练习 1编译与跑通把代码 2-1、2-2 与Make/{files,options}配好编译并在一个含0/U、0/T的算例副本上运行。判定wmake无 error日志出现T: nIter ...与integral(T) ...以End结束。练习 2守恒性验证用代码 2-4 计算相对漂移。判定在无源Su0、Sp不改变总量且边界封闭的设置下相对漂移 1%若超阈值能指出优先排查 phi 的连续性。练习 3源项稳定性实验把Sp从-0.1改为0.1正系数 源项随 T 增大而增大观察残差与数值稳定性。判定能观察到稳定性变差迭代增多、可能发散或出现异常值并解释正 Sp 破坏对角占优的机制。练习 4格式影响把div(phi,T)的格式从Gauss limitedLinear 1换成Gauss linear中心差分在强对流设置下对比结果。判定能观察到非物理过冲/欠冲增加例如 T 超出边界值的范围并能用有界性 vs 精度解释取舍。练习 5思考题无标准答案把这个求解器改造成温度场求解器即T就是温度。列出需要修改的三处。验证要点(a)T的dimensions是否为温度量纲(b)DT是否应改为热扩散系数并考虑是否需除以密度与比热© 是否考虑源项量纲变为温度/时间第 08 篇将给出物性字典的规范读法。五、小结与下一篇预告本篇是你从会写求解器到会写有用的求解器的转折点。四条要点方程四项对应四个算子ddt/div/laplacian/源项源项半隐式化负系数进fvm::Sp、常数进Su是稳定性的关键加场必改三处配置场文件、solver 条目、格式验证要分三层残差、守恒性、极限行为并用gSum(T*mesh.V())把守恒性变成可量化的数字。第 08 篇《物性与热物性模型扩展》将解决本篇留下的一个工程隐患DT不应该写死在代码里。我们要解剖thermophysicalProperties字典与hConst/hPolynomial/janaf三类热力学模型搞清楚物性字典如何变成运行期类实例并给出自定义物性的规范做法——那是把求解器从玩具变成工程工具的必经一步。本篇认知问题回显FAQQ1标量输运方程在 OpenFOAM 里怎么写A写成fvm::ddt(T) fvm::div(phi, T) - fvm::laplacian(DT, T) fvm::Sp(Sp, T) Su其中 T 是待求标量场phi 是面通量面心场 surfaceScalarFieldDT 是扩散系数Sp/Su 是源项的线性化系数。ddt 为时间项、div 为对流项、laplacian 为扩散项前三项都用 fvm:: 隐式处理源项用 fvm::Sp隐式负系数部分加 Su显式常数部分。装配为 fvScalarMatrix 后调用 solve()。Q2源项怎么处理才既物理又稳定A把源项线性化为 S(T) ≈ Sp·T Su其中 Sp 与 T 无关可为场。耗散或消耗型物理使 Sp 为负写进 fvm::Sp(Sp, T) 隐式进矩阵对角增强对角占优、提升稳定性常数或生成部分放 Su 显式进右端。强烈非线性源项如 Arrhenius 速率也可用同样方式线性化配合外迭代逐步更新。判据是源项能否写成关于 T 的线性形式负系数部分隐式化是关键。Q3加了新场 T 之后必须在算例里补哪些配置A三处。第一是 0/T 场文件含 dimensions、internalField 与各 patch 的 type否则 MUST_READ 报 cannot find file。第二是 system/fvSolution 的 solvers 块中为 T 增加线性求解器条目solver/preconditioner/tolerance/relTol否则运行期报缺 solver 设置被 solve() 的场必须有 solver 条目而 NO_READ 的辅助场如 DT/Sp/Su不需要。第三是 system/fvSchemes 中补齐 T 用到的项格式如 div(phi,T)、laplacian(DT,T)若用了 default none 则遗漏会立即报 keyword undefined。Q4怎么判断标量场算得对不对A分三个层次验证。数值层看残差是否单调下降到 fvSolution 设定的容差以下由 SolverPerformance 的 nIterations/initialResidual/finalResidual/converged 判定。守恒层在无源且边界封闭时检查标量总量是否守恒用 gSum(T.primitiveField()*mesh.V()) 计算体积分并比较首末相对漂移工程判据常取 1%若漂移过大首要排查 phi 是否满足连续性div(phi)≈0因为标量守恒性继承自流场连续性。极限层把扩散系数设极大看场是否趋于均匀或把速度设为零对比纯扩散行为。Q5这个求解器可以改成哪些实用模型A标量输运方程形式通用可直接改造为温度场T 为温度DT 改为热扩散并考虑密度与比热、组分浓度T 为质量或摩尔分数DT 为组分扩散系数、被动标量与示踪剂、盐度或其他输运量也可作为湍流标量如 k、omega 类方程的基础结构。改造时需同步调整 T 的量纲、扩散系数定义与源项量纲并让辅助场从配置字典读取而非写死在代码里。