OpenFOAM二次开发教程（06）：有限体积离散与方程装配——fvMatrix 与 fvm/fvc 算子

OpenFOAM二次开发教程（06）：有限体积离散与方程装配——fvMatrix 与 fvm/fvc 算子

版本与事实声明

• fvMatrix的求解接口见官方 DoxygenfvMatrix.H：提供SolverPerformance<Type> solve(fvMatrix<Type>&, const word&)等重载，“Solve returning the solution statistics given convergence tolerance”。
• fvm（隐式有限体积算子）与fvc（显式有限体积算子）两族函数的可用清单以本机$FOAM_SRC/finiteVolume/finiteVolume/fvm|fvc/与官方 Doxygen 为准。
• SolverPerformance是求解器返回值类型（含初始残差、最终残差、迭代次数等统计），字段名以官方源码为准。
• 本文示例为教学用 Poisson 方程求解器，物理设置与边界值均为示例，不代表任何标准规定。

一句话结论：在 OpenFOAM 里"求解一个方程"分两步——用fvm::（隐式）与fvc::（显式）算子把偏微分方程装配成fvMatrix，再调用solve()得到SolverPerformance统计；判据不是"有没有解出来"，而是"残差是否降到fvSolution设定的容差以下"。

〇、本篇要解决的认知问题

• Q1：为什么 OpenFOAM 不像传统 CFD 代码那样"写个线性系统再交给求解器"，而是搞出fvMatrix这一层？
• Q2：fvm::与fvc::到底差在哪？什么时候必须用隐式、什么时候只能用显式？
• Q3：solve()返回的SolverPerformance里有什么？怎么用它判断"这次迭代算不算成功"？
• Q4：fvm::Sp与fvm::Su有什么区别？为什么隐式处理源项能显著改善稳定性？
• Q5：为什么fvSchemes里的格式选择会直接影响fvMatrix的系数？两者是怎么联动的？

一、机制解析

1.1 从偏微分方程到 fvMatrix：OpenFOAM 的"装配"哲学

传统 CFD 代码的写法是：先离散、拿到系数矩阵的稀疏结构、再调用线性求解器（如 PETSc、hypre）。OpenFOAM 的写法不同——它把"离散"和"矩阵装配"合并成一个可由方程表达式驱动的过程：

偏微分方程（人类写法） ∂T/∂t + ∇·(U T) - ∇·(ν ∇T) = S │ │ 用 fvm:: / fvc:: 写成 C++ 表达式 ▼ fvMatrix<T> 方程 ( fvm::ddt(T) + fvm::div(phi, T) - fvm::laplacian(nu, T) == fvm::Sp(a, T) + Su ) │ │ solve() 触发：按 fvSchemes 的格式生成系数 → 交给 fvSolution 指定的求解器 ▼ SolverPerformance （残差、迭代次数等统计）

为什么这样设计：它让"方程"成为代码里的一等公民。改一项离散格式，只需改fvSchemes（配置），不用改代码；换一个线性求解器，只需改fvSolution。这就是"配置驱动"在算法层的体现，也是本系列第 01 篇决策表里"能靠字典解决的不写 C++"的深层原因。

1.2 fvm 与 fvc：隐式与显式的分工

两族算子的区别只有一个词：这个项进不进矩阵。

核心规则：方程左边（要解的量）用fvm::；当作"已知量"用fvc::。

常见组合示例：

// 稳态标量输运方程：隐式处理对流与扩散，显式给出源项（Su）fvScalarMatrixTEqn(fvm::div(phi,T)// 对流：隐式（通量 phi 视为已知）-fvm::laplacian(DT,T)// 扩散：隐式（扩散系数 DT 视为已知）==Su// 源项：显式（直接作为右端项）);TEqn.solve();// 触发装配系数 + 调用 fvSolution 中的线性求解器

反直觉点：fvm::div(phi, T)里的phi本身通常是用fvc::interpolate(U) & mesh.Sf()或createPhi.H得到的显式量。也就是说，同一个方程里可以混用两族算子——"哪些项进矩阵"是你的物理/数值决策，不是语法限制。

1.3 fvm::Sp 与 fvm::Su：源项的隐式化

源项处理是二次开发中最常被低估的技术点。

• fvm::Su(Su_value, T)：把源项显式加入右端项（Su_value按给定值计算）。实现简单，但对大源项不稳定。
• fvm::Sp(Sp_value, T)：把源项隐式处理为"系数 × 待求量"，即把它写进矩阵对角。当Sp_value的符号使其增强对角占优时，稳定性显著提升。

为什么这对你重要：任何"耗散型"物理（如化学反应消耗、辐射吸收、多孔介质阻力）都应该优先隐式化。经验法则：如果源项可以写成"负系数 × 待求量 + 常数"，那么把负系数部分放进fvm::Sp、常数部分放Su，通常能大幅改善收敛——这条经验在第 13 篇的fvModels（如semiImplicitSource，即"半隐式源项"）里被官方直接采用。

1.4 solve() 与 SolverPerformance：怎么判定"算成功"

fvMatrix.H提供的求解接口签名为：

template<classType>SolverPerformance<Type>solve(fvMatrix<Type>&,constword&);

官方文档对它的一句说明是"Solve returning the solution statistics given convergence tolerance"——返回的是"求解统计"，判据是收敛容差。SolverPerformance内含（字段名以官方源码为准）：

• 初始残差 / 最终残差（initial、final residual）；
• 迭代次数（nIterations）；
• 是否收敛（converged 标志）；
• 求解器名与场名。

工程含义：solve()不会因为"没收敛"而抛异常，它只是把结果告诉你。你必须自己判断。这就是为什么大量求解器代码里会写：

SolverPerformance<scalar>perf=TEqn.solve();// 或者更常见的写法：solve() 后由外层检查残差与最大迭代数

最佳实践：在自建求解器里，把每步的SolverPerformance打印/记录到日志（Info或写入 CSV），让"收敛轨迹"成为可审计数据。第 15 篇性能优化正是靠这份轨迹来判断"是网格问题还是求解器设置问题"。

1.5 fvSchemes 与 fvMatrix 的联动

这是最容易让初学者困惑的一环：fvMatrix的系数并不是写死的，而是由fvSchemes在装配时决定的批量选择。

具体地：

• fvm::div(phi, T)装配时，会去fvSchemes的divSchemes里找div(phi,T)对应的格式（如Gauss linearUpwind grad(T)）；
• fvm::laplacian(DT, T)会去laplacianSchemes找对应格式（如Gauss linear corrected）；
• fvm::ddt(T)会去ddtSchemes（如Euler、backward）。

推论：同一份 C++ 求解器代码，配上不同的fvSchemes，装配出的矩阵完全不同。所以"结果不对"时，先别怀疑代码——先核对fvSchemes。这是第 04 篇"分层排查"纪律在第 06 篇的具体落地。

二、完整代码与逐行剖析

代码 2-1：手写 Poisson 方程求解器myPoissonFoam.C

/*---------------------------------------------------------------------------*\ myPoissonFoam.C —— 教学用 Poisson 方程求解器 求解：laplacian(T) = f （稳态扩散/势场问题的最小形式） 用途：演示 fvMatrix 装配 + solve() + SolverPerformance 的完整闭环。 \*---------------------------------------------------------------------------*/#include"fvCFD.H"// fvMesh、volField、fvm/fvc、fvMatrix、solve 等// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //intmain(intargc,char*argv[]){#include"setRootCase.H"#include"createTime.H"#include"createMesh.H"// ---- 1. 读入待求场 T（Poisson 方程的解）----Info<<"Reading field T\n"<<endl;volScalarFieldT(IOobject("T",runTime.timeName(),mesh,IOobject::MUST_READ,IOobject::AUTO_WRITE),mesh);// ---- 2. 定义一个扩散系数 DT（无量纲或与方程一致，示例值）----// dimensionedScalar 用 T.dimensions() 派生维度，避免手写指数出错constdimensionedScalarDT("DT",dimless,1.0);// ---- 3. 跑一个“伪时间循环”，让控制权交给 runTime 与写出机制 ----// 对稳态问题，可在循环内反复解方程直至收敛；这里演示最小闭环。while(runTime.loop()){Info<<"Time = "<<runTime.timeName()<<nl<<endl;// ---- 4. 装配 Poisson 方程：laplacian(T) == 0（示例：无源 Poisson/Laplace）----// 注意：fvm::laplacian 是隐式项，会进入矩阵对角与邻点系数；// 右端项若为场（如源项），可用 fvc:: 计算后放在 == 右侧。fvScalarMatrixTEqn(-fvm::laplacian(DT,T)// 取负号使矩阵对角为正定（符号约定，见剖析）);// ---- 5. 求解并获取统计信息 ----SolverPerformance<scalar>perf=TEqn.solve();// ---- 6. 显式检查收敛性：solve() 不会替代你做判断 ----// perf 中提供初始/最终残差与迭代次数（字段名以官方源码为准）Info<<"T solve: "<<"nIter = "<<perf.nIterations()<<" initial = "<<perf.initialResidual()<<" final = "<<perf.finalResidual()<<" converged = "<<perf.converged()<<nl<<endl;// ---- 7. 边界重新求值：内部场变了，边界必须同步（第 05 篇纪律）----T.correctBoundaryConditions();// ---- 8. 写出 ----runTime.write();Info<<"ExecutionTime = "<<runTime.elapsedCpuTime()<<" s\n"<<endl;}Info<<"End\n"<<endl;return0;}

逐行剖析：

• DT用dimless（无量纲）是示例选择：真实问题里 DT 的维度必须与方程一致。用T.dimensions()派生或明确写出dimensionSet都行，但要保证一致性（第 05 篇的纪律）。
• fvm::laplacian(DT, T)是隐式项：它把"每个单元的 T 与其邻点 T 的线性组合"变成矩阵系数。这一行就是"装配"的核心动作。
• 取负号- fvm::laplacian(DT, T)：这是矩阵符号约定的结果。OpenFOAM 的矩阵装配遵循"对角占优、系数符号统一"的内部约定，不同方程形式下需要调整符号使矩阵良态。实践建议：直接对照官方同类求解器的符号写法，不要凭直觉试——这是最容易"看起来对、其实矩阵已病态"的地方。
• SolverPerformance<scalar> perf = TEqn.solve();：官方签名即SolverPerformance<Type> solve(fvMatrix<Type>&, const word&)（Impl版本还有带word的重载）。返回值含收敛统计。
• perf.nIterations()/initialResidual()/finalResidual()/converged()：把"是否收敛"变成可打印、可审计的数值。这是二次开发最该养成的习惯：不要相信"跑了没报错"，要相信残差数字。
• T.correctBoundaryConditions()：这行不能省。解完方程后内部场更新了，边界必须按各自类型重新计算，否则内部/边界不自洽（第 05 篇"两套数据"的直接后果）。
• runTime.write()：按controlDict的写出设置落盘。稳态问题可以只用少量迭代步，靠executeControl之类的函数对象做收敛监控。

代码 2-2：配套的system/fvSchemes与system/fvSolution（关键块）

// ---------- system/fvSchemes（节选）---------- ddtSchemes { default Euler; } gradSchemes { default Gauss linear; } divSchemes { default none; } // 本方程无对流项，保持 none 即可 laplacianSchemes{ default Gauss linear corrected; } // Poisson 主角：拉普拉斯格式 interpolationSchemes { default linear; } snGradSchemes { default corrected; } // ---------- system/fvSolution（节选）---------- solvers { T { // 拉普拉斯离散得到的矩阵通常对称，用 PCG + DIC 是常规选择 solver PCG; preconditioner DIC; tolerance 1e-8; // 绝对容差 relTol 0.01; // 相对容差：早期迭代放宽以省时间 } }

逐行剖析：

• laplacianSchemes的corrected是"含非正交修正"：网格非正交时必须用（否则结果有误），代价是多一次非正交循环。判据：checkMesh报告的非正交角越大，越需要修正。
• divSchemes { default none; }：没有对流项就显式声明 none，这样如果误加了fvm::div会立刻报"缺格式"，而不是静默用错格式。这是第 04 篇"显式声明"纪律的延续。
• T的求解器选PCG+DIC：对称矩阵配对称求解器。如果错配成非对称求解器（如PBiCGStab），不是不能算，而是收敛慢、内存高——选错不报错但代价很大，这正是第 15 篇要优化的对象。
• tolerance与relTol的分工：relTol是"相对初始残差的下降比例"，用于前几步迭代；tolerance是绝对下限。经验法则：relTol用于"省早期迭代时间"，tolerance用于"保证最终精度"。

代码 2-3：收敛性验证脚本（POSIX Shell）

#!/bin/sh# verify_poisson.sh —— 在官方算例副本上验证 myPoissonFoam 的残差下降行为# 用法：CASE_SRC=<含 0/T 的算例路径> sh verify_poisson.shset-euCASE_SRC="${CASE_SRC:?请用 CASE_SRC=... 指定一个含 0/T 的算例（可用官方传热/势场类算例副本）}"WORK="$PWD/_poisson_case"rm-rf"$WORK";cp-r"$CASE_SRC""$WORK"cd"$WORK"blockMesh>log.blockMesh2>&1||{echo"[FAIL] blockMesh 失败";exit1;}myPoissonFoam-case.>log.myPoissonFoam2>&1||{echo"[FAIL] 求解器失败";exit1;}echo"== 残差轨迹（每步一行）=="grep-E"T solve:"log.myPoissonFoam||trueecho"== 判定 =="# 判据 1：出现求解统计行grep-q"nIter"log.myPoissonFoam&&echo"[OK] 已输出求解统计（含迭代次数与残差）"# 判据 2：日志以 End 正常结束grep-q"^End"log.myPoissonFoam&&echo"[OK] 求解器正常结束"# 判据 3：至少产生一个时间目录（说明 runTime.write 生效）n=$(ls-d[0-9]*2>/dev/null|grep-v'^0$'|wc-l)["$n"-ge1]&&echo"[OK] 产生了$n个非初始时间目录"||echo"[WARN] 未产生新时间目录，请检查写设置"

逐行剖析：

• 用grep -E "T solve:"抽取每步残差行：这就是第 15 篇性能优化所需的"残差轨迹"原始数据，本篇先把它落下来。
• 两个硬判据（有统计行、以End结束）+ 一个软判据（有时间目录）：区分"必须通过"与"视配置而定"，避免把合理情况误判为失败。
• 依然在官方算例副本上验证（铁律 7），不改动原算例。

三、常见报错与排查

报错 3-1：--> FOAM FATAL ERROR: ... keyword laplacian(DT,T) is undefined in dictionary .../system/fvSchemes。
现象：装配时找不到拉普拉斯项的格式。根因：fvSchemes的laplacianSchemes用了default none又没有为该具体项指定格式，或格式键拼写与代码中的表达式不完全一致（laplacian(DT,T)的括号内容必须匹配）。解法：核对报错给出的键名，在laplacianSchemes中补齐对应条目；注意键名里的变量名要与代码一致。

报错 3-2：求解"收敛"了但残差没有下降（initial与final几乎相等）。
现象：SolverPerformance显示迭代次数很低、残差几乎没降。根因：容差设置过宽（tolerance/relTol太大），或方程装配有问题（例如把所有项都放到右端显式处理，导致矩阵几乎只有对角）。解法：先调紧tolerance（如1e-8）验证；若残差仍不降，检查方程装配——确保待求项用了fvm::而非fvc::。

报错 3-3：--> FOAM FATAL ERROR: ... dimension mismatch ...。
现象：装配方程时报维度冲突。根因：方程各项维度不一致（例如fvm::laplacian(DT, T)中DT的维度与T、网格尺度的组合无法得到与其它项相同的维度）。解法：打印各项维度对账；优先让系数从被作用场派生维度，不要手写dimensionSet指数。

报错 3-4：迭代次数爆炸（几百上千次）甚至发散。
现象：nIter极大，或残差上升到NaN。根因：矩阵病态——常见于显式处理了本应隐式的项（大源项用了Su而非Sp）、泊松方程缺少参照（纯 Neumann 问题解不唯一）、或fvSchemes的对流格式在强对流下不稳定。解法：把大源项改隐式（fvm::Sp）；纯 Neumann 问题在fvSolution里设pRefCell/pRefValue固定参照；强对流改用更稳的格式或降低时间步。

报错 3-5：编译报no matching function for call to ...（算子调用不匹配）。
现象：fvm::xxx(...)参数类型不对。根因：算子的参数类型有严格要求（例如fvm::div期望面通量surfaceScalarField，你传了体心速度volVectorField）。解法：打开官方 Doxygen 找到该算子的签名，或看官方求解器源码里的实际调用（铁律 1）。最稳的做法：直接抄官方求解器的同类表达式，再改名字。

四、动手练习

• 练习 1（最小闭环）：用代码 2-1 与 2-2 编译并运行myPoissonFoam。判定：wmake无 error；日志出现T solve:行且含nIter；日志以End结束。
• 练习 2（残差实验）：把fvSolution中T的tolerance分别设为1e-3、1e-6、1e-10，各运行一次。判定：能观察到nIterations随tolerance收紧而增大（或finalResidual随之变小）；能用自己的话解释"为什么容差与迭代次数是权衡关系"。
• 练习 3（隐式源项对比）：用fvm::Sp与显式Su两种方式实现同一个"耗散型"源项，分别运行。判定：能观察到隐式版本迭代次数更少或更稳定；能说明"为什么负系数源项隐式化能增强对角占优"。
• 练习 4（格式联动）：在fvSchemes中把laplacianSchemes从Gauss linear corrected改为Gauss linear uncorrected（非正交网格下）。判定：能观察到结果或迭代行为发生变化（对非正交较强的网格更明显），并能说明"为什么 uncorrected 在非正交网格上会引入误差"（此结论需结合checkMesh输出的非正交角度量）。
• 练习 5（思考题，无标准答案）：给定方程∂T/∂t + ∇·(U T) = ∇·(ν ∇T)，写出每一项在 OpenFOAM 中的算子写法并说明隐式/显式选择。验证要点：(a) 时间导数与对流、扩散是否用fvm::；(b) 通量phi是否用面心场并（通常）显式计算；© 是否在fvSchemes中补齐ddt/div/laplacian三类格式（第 07 篇给出完整实现）。

五、小结与下一篇预告

本篇揭示了 OpenFOAM 最核心的机制：方程即代码，装配即离散。四条要点请记住——fvm::进矩阵、fvc::不算系数；fvm::Sp隐式化源项能显著改善稳定性；solve()返回SolverPerformance，收敛要自己判断、要留残差轨迹；fvMatrix的系数由fvSchemes决定，“结果不对先查格式”。

第 07 篇《第一个二次开发实战：自定义标量输运求解器》将把前六篇的知识合成一个真正有用的东西：一个带对流、扩散与源项的标量输运求解器myScalarTransportFoam，并在官方算例上验证守恒性与残差行为——那是你第一次写出"物理正确"的自定义求解器。

本篇认知问题回显（FAQ）

Q1：为什么 OpenFOAM 要有 fvMatrix 这一层，而不是直接给线性系统？
A：fvMatrix 把"离散"与"矩阵装配"合并为一个由方程表达式驱动的过程。你用人话写出方程（如 - fvm::laplacian(DT, T)），solve() 触发时按 fvSchemes 选择离散格式生成系数，再交给 fvSolution 指定的线性求解器。这样方程成为代码里的一等公民：换离散格式只改 fvSchemes 配置，换线性求解器只改 fvSolution，不必改代码，也避免了手写稀疏矩阵结构的繁琐与出错。

Q2：fvm:: 与 fvc:: 的区别是什么？
A：fvm:: 是隐式有限体积算子，把该项离散进 fvMatrix（成为矩阵系数），返回值是 fvMatrix，稳定性高，典型项有 fvm::ddt、fvm::div、fvm::laplacian、fvm::Sp。fvc:: 是显式有限体积算子，把该项直接计算成场（不成为系数），返回 volField 或 surfaceField，稳定性较低，典型项有 fvc::grad、fvc::div、fvc::interpolate。核心规则是：方程左边要解的量用 fvm::，当作已知量用的用 fvc::；同一方程中可混用两族算子。

Q3：solve() 返回的 SolverPerformance 里有什么，怎么用？
A：官方 fvMatrix.H 提供 SolverPerformance solve(fvMatrix&, const word&)，说明为"Solve returning the solution statistics given convergence tolerance"，即返回求解统计。统计包含初始残差、最终残差、迭代次数、是否收敛标志、求解器名与场名等（字段名以官方源码为准）。关键在于 solve() 不会因未收敛而抛异常，必须由你判断：把 nIterations/initialResidual/finalResidual/converged 打印或记录成残差轨迹，才能判定这次迭代是否成功。

Q4：fvm::Sp 与 fvm::Su 的区别？为什么隐式源项更稳定？
A：fvm::Su 把源项显式加到方程右端项（按给定值计算），实现简单但对大源项不稳定。fvm::Sp 把源项隐式处理为"系数乘待求量"，写进矩阵对角。当该系数为负（耗散/消耗型物理，如化学反应消耗、辐射吸收、多孔介质阻力）时，隐式化会增强矩阵对角占优，从而显著改善稳定性与收敛速度。经验做法是把源项写成"负系数×待求量 + 常数"，负系数部分进 fvm::Sp、常数部分进 fvm::Su，官方的 semiImplicitSource 即半隐式源项思路。

Q5：fvSchemes 的格式选择如何影响 fvMatrix？
A：fvMatrix 的系数不是写死的，而是 solve() 时按 fvSchemes 中的格式动态生成的。fvm::ddt(T) 查 ddtSchemes，fvm::div(phi,T) 查 divSchemes，fvm::laplacian(DT,T) 查 laplacianSchemes，梯度与插值查 gradSchemes/interpolationSchemes/snGradSchemes。因此同一份 C++ 求解器代码配上不同 fvSchemes 会装配出完全不同的矩阵，精度与稳定性随之变化。实践含义是"结果不对先核对 fvSchemes 再怀疑代码"；若某类项没有对流项，可用 default none 使遗漏立即报错而非静默用错格式。