第 9 章 有限体积离散框架：fvMatrix、fvSchemes 与 fvSolution
前面八章铺垫了所有基础设施：字典、注册表、RTS、基元、场、网格、边界条件。本章进入 OpenFOAM 的"心脏"——把连续偏微分方程离散成线性方程组。核心是三个类：
fvMatrix（离散后的稀疏矩阵）、fvSchemes（离散格式配置）、fvSolution（线性求解配置）。理解了本章，你就能看懂任意一个求解器里的fvVectorMatrix UEqn(...)到底在做什么。
9.1 宏观图景：从 PDE 到线性方程
OpenFOAM 求解一个输运方程（以动量方程为例）时，走的是这样一条路径：
连续 PDE
│ 用有限体积法在网格上离散（空间: 面通量; 时间: ddt）
▼
fvMatrix<Type> ← fvm:: 隐式算子组装（对角 + 上/下三角 + 边界系数 + 源项）
│ 映射为通用稀疏格式
▼
lduMatrix ← 与网格无关的"下-对角-上"稀疏矩阵
│ RTS 选择求解器
▼
PCG / PBiCGStab / GAMG ← 迭代求解（第 11 章）
│
▼
新的场值
关键设计是：fvMatrix 直接继承自 lduMatrix（见 9.4），因此从有限体积离散到线性代数求解之间没有额外的转换层——离散器直接产出求解器能吃的矩阵。这是 OpenFOAM "一体化"设计的一个缩影。
9.2 fvSchemes / schemesLookup：离散格式的配置中心
9.2.1 目录表
schemesLookup.H（src/OpenFOAM/matrices/schemes/，第 30–45 行）用一张表列出了 fvSchemes 字典的所有条目类别：
Description
Selector class for finite area/finite volume differencing schemes.
The file will contain these types of entries:
ddtSchemes | ddt | dict | no | none
d2dt2Schemes | d2dt2 (currently only finite-area) | dict | no | none
interpolationSchemes | interpolation | dict | no | linear
divSchemes | div | dict | yes | -
gradSchemes | grad | dict | yes | -
lnGradSchemes | line-normal grad (finite-area) | dict | no | corrected
snGradSchemes | surface-normal grad | dict | no | corrected
laplacianSchemes| laplacian | dict | yes | -
fluxRequired | | dict | no | -
这张表就是 fvSchemes 字典的"语法规范"。要点：
| 条目 | 含义 | 是否有默认 |
|---|---|---|
ddtSchemes | 时间导数格式（Euler/backward/CrankNicolson…） | 无默认 |
interpolationSchemes | 面插值格式 | 默认 linear |
divSchemes | 对流通量/散度格式 | 必需 |
gradSchemes | 梯度格式 | 必需 |
snGradSchemes | 面法向梯度格式 | 默认 corrected |
laplacianSchemes | 拉普拉斯格式 | 必需 |
fluxRequired | 需要通量的场（如 p） | 无 |
fvSchemes 只是 schemesLookup 的一个薄壳（fvSchemes.H 第 58–61 行 class fvSchemes : public schemesLookup），后者才是真正持有 IOdictionary 并做"按名查格式、否则取默认"的地方。内部把每个类别封装成 lookupDetail（schemesLookup.H 第 87–128 行），持有 name_、dict_、default_（一个 ITstream 默认 token 流），并提供 lookup(name)"取具名格式或默认"。
9.2.2 运行期可改
schemesLookup.H 第 52–55 行有个重要说明：
Note
The schemes data are treated as \c READ_MODIFIED even if the
requested \p readOption is nominally MUST_READ or READ_IF_PRESENT.
This allows run-time modifications to behave as expected.
离散格式字典总是按 READ_MODIFIED 处理——即使你以 MUST_READ 打开。这样运行期改 fvSchemes 就能生效。fvSolution 同理（solution.H 第 66–68 行）。这解释了"为什么改松弛因子/换差分格式不用重启程序"。
9.2.3 为什么 fvMesh 要继承 fvSchemes
fvSchemes.H 第 32–33 行：
fvMesh is derived from fvSchemes so that all fields have access to the
fvSchemes from the mesh reference they hold.
因为每个 GeometricField 都持有 mesh 引用，而 mesh 继承自 fvSchemes——于是任何场都能直接 mesh.schemes().ddt(...) 拿到自己的时间格式，无需全局查找。fvSolution 同理。这就是第 8 章 fvMesh 多重继承的深层动机。
9.3 fvSolution / solution：求解器与松弛配置
solution.H（src/OpenFOAM/matrices/solution/，第 30–64 行）说明了 fvSolution 字典的结构：
\par solvers { }
A sub-dictionary listing of linear solver settings. For example,
solvers
{
p
{
solver PCG;
preconditioner DIC;
tolerance 1e-06;
relTol 0.05;
}
}
...
PIMPLE
{
nOuterCorrectors 2;
nCorrectors 1;
nNonOrthogonalCorrectors 0;
}
solution 类的私有数据（solution.H 第 93–130 行）精确对应字典内容：
class solution : public IOdictionary
{
dictionary cache_; // 待缓存临时字段
bool caching_; // 缓存开关
dictionary fieldRelaxDict_; // 场的松弛因子
mutable HashPtrTable<Function1<scalar>> fieldRelaxCache_;
dictionary eqnRelaxDict_; // 方程的松弛因子
mutable HashPtrTable<Function1<scalar>> eqnRelaxCache_;
autoPtr<Function1<scalar>> fieldRelaxDefault_;
autoPtr<Function1<scalar>> eqnRelaxDefault_;
dictionary solvers_; // 各场的线性求解器设置
};
三类配置：
solvers_：每个场（p、U、k…）一套线性求解器设置（solver/preconditioner/tolerance/relTol/maxIter）。第 11 章会看到它如何驱动 RTS 选择PCG/PBiCGStab/GAMG。- 松弛因子：分"场的松弛"（
fieldRelaxDict_）与"方程（矩阵）的松弛"（eqnRelaxDict_，对应fvMatrix::relax()），且支持Function1（时间相关的松弛）。 cache_：cacheTemporaryObjects，即第 3 章提到的临时对象缓存——把某些中间量算一次缓存起来复用。
注意它还持有 HashPtrTable<Function1<scalar>> 缓存：松弛因子可以是时间的函数（Function1），解析一次后缓存复用。
9.4 fvMatrix<Type>：离散后的方程矩阵
9.4.1 类结构与数据
fvMatrix.H（src/finiteVolume/fvMatrices/fvMatrix/，第 118–177 行）：
template<class Type>
class fvMatrix
:
public refCount,
public lduMatrix
{
public:
typedef GeometricField<Type, fvPatchField, volMesh> psiFieldType;
typedef GeometricField<Type, fvsPatchField, surfaceMesh> faceFluxFieldType;
private:
const psiFieldType& psi_; // 被求解的场（引用）
PtrList<fvMatrix<Type>> subMatrices_; // 组装用的子矩阵
bool useImplicit_; // 是否用隐式组装
word lduAssemblyName_; // ldu 组装名
label nMatrix_; // 叠加的矩阵数量
dimensionSet dimensions_; // 维度
Field<Type> source_; // 源项（右端）
FieldField<Field, Type> internalCoeffs_; // 边界的"内部"系数
FieldField<Field, Type> boundaryCoeffs_; // 边界的"源"系数
mutable std::unique_ptr<faceFluxFieldType> faceFluxCorrectionPtr_; // 非正交修正通量
};
逐项解读（这是理解 OpenFOAM 离散的关键）：
psi_：被求解的场，以 const 引用持有（注释说明"求解时才转为非 const"）。这避免了拷贝，也解释了为什么UEqn.solve()必须用const_cast去改psi_。- 继承
lduMatrix：fvMatrix就是一个 LDU 稀疏矩阵。它继承的三组系数（diag/upper/lower，见 9.4.2）来自lduMatrix。 source_：离散方程的右端项（显式部分）。internalCoeffs_/boundaryCoeffs_：边界条件的贡献，类型是FieldField<Field, Type>（第 5 章的"场的场"——每个 patch 一个Field）。这正是第 8 章fvPatchField::valueInternalCoeffs()等函数的返回值去处。faceFluxCorrectionPtr_：非正交修正通量。非正交网格上，拉普拉斯/梯度需要一次"修正"迭代，这个指针缓存修正量（见第 10 章）。subMatrices_/useImplicit_/nMatrix_/lduAssemblyName_：支撑"隐式组装"与"矩阵叠加"，是新版为性能和灵活性引入的机制。
9.4.2 存储格式：为什么叫 "LDU"
lduMatrix 用三个数组存储整个稀疏矩阵（第 11 章会深入，此处先建立直觉）：
| 数组 | 长度 | 含义 |
|---|---|---|
diag_ | nCells | 对角元素（每个单元自己的系数） |
upper_ | nInternalFaces | 上三角元素（面 owner→neighbour 的系数） |
lower_ | nInternalFaces | 下三角元素（面 neighbour→owner 的系数） |
为什么有限体积矩阵天然适配 LDU？ 因为离散方程的耦合只发生在"共享一个面的相邻单元"之间：
A_P φ_P + Σ_{内部面 f} (a_f φ_N) = b_P
每个内部面正好对应矩阵里的一对 off-diagonal 元素（一上一下）。面的 owner 单元对应"行"，neighbour 对应"列"。因此矩阵的非零结构完全由网格拓扑决定，upper_/lower_ 的长度恰好等于内部面数。这就是 OpenFOAM 能用如此紧凑的存储（3 个数组）表达任意多面体网格上离散矩阵的原因。
9.4.3 fvm:: 与 fvc::：隐式 vs 显式
离散算子分为两个命名空间（src/finiteVolume/finiteVolume/fvm/ 与 fvc/）：
fvm::（finite volume method，隐式）：返回tmp<fvMatrix<Type>>。算子把未知量 φ 留在矩阵里（对角/上下三角），组装成一个待解方程。 例如fvm::ddt(rho, U)、fvm::div(phi, U)、fvm::laplacian(mu, U)。fvc::（finite volume calculus，显式）：返回tmp<GeometricField<...>>。算子对已知场求值，得到一个算好的场。 例如fvc::grad(p)、fvc::div(phi)、fvc::laplacian(nu, U)、fvc::interpolate(rho)。
fvm:: 与 fvc:: 的 fvmSup.H 还提供 Su（显式源）与 Sp（隐式源），用于给方程添加源项（第 54–119 行：Su 返回 tmp<fvMatrix<Type>>，Sp 同理）。
典型用法（动量方程，来自求解器）：
fvVectorMatrix UEqn
(
fvm::ddt(rho, U) + fvm::div(phi, U)
+ turbulence->divDevReff(U)
);
UEqn.relax();
solve(UEqn == -fvc::grad(p));
读法：fvm:: 的部分构成矩阵（未知量 U），fvc:: 的部分构成已知的源（-∇p）。 solve(UEqn == ...) 里 == 的作用是把显式项挪到右端并求解。
9.4.4 边界如何"加进"矩阵
fvMatrix 的受保护方法（fvMatrix.H 第 186–242 行）里有一组关键函数：
//- Add patch contribution to internal field
template<class Type2> void addToInternalField(...) const;
template<class Type2> void subtractFromInternalField(...) const;
// Matrix completion functionality
void addBoundaryDiag(scalarField& diag, const direction cmpt) const;
void addCmptAvBoundaryDiag(scalarField& diag) const;
void addBoundarySource(...);
流程是：
- 离散内部面时，写入
diag_/upper_/lower_（来自lduMatrix）； - 遍历边界面时，调用每个 patch 的
valueInternalCoeffs/valueBoundaryCoeffs/gradientInternalCoeffs/gradientBoundaryCoeffs（第 8 章）； - 把这些返回的
Field累加到internalCoeffs_/boundaryCoeffs_； addBoundaryDiag()把"内部系数"加到diag_；addBoundarySource()把"边界系数（源）"加到source_。
这就是"边界条件进入线性方程组"的完整机制。 也正因如此，fvMatrix 必须先知道求解的场（psi_）——因为边界条件需要 psi_ 的 patch 信息。
9.5 求解流程：solve()
9.5.1 分流：segregated vs coupled
fvMatrixSolve.C 第 60–107 行的 solveSegregatedOrCoupled()：
template<class Type>
Foam::SolverPerformance<Type> Foam::fvMatrix<Type>::solveSegregatedOrCoupled
(
const dictionary& solverControls
)
{
...
// Do not solve if maxIter == 0
if (solverControls.getOrDefault<label>("maxIter", -1) == 0)
{
return SolverPerformance<Type>();
}
word type(solverControls.getOrDefault<word>("type", "segregated"));
if (type == "segregated")
{
return solveSegregated(solverControls);
}
else if (type == "coupled")
{
return solveCoupled(solverControls);
}
...
}
- segregated（分离式）：每个场（p、U、k…）各自求解，是传统 SIMPLE/PISO 方式。
- coupled（耦合式）：多个场联立求解（用于强耦合场景）。
求解器参数从 solverControls（即 fvSolution 里该场的子字典）读取：maxIter、solver、preconditioner、tol、relTol 等。
9.5.2 从 fvMatrix 到 lduMatrix::solver
solveSegregated() 内部会把 fvMatrix（作为 lduMatrix）交给一个由 RTS 选择的 lduMatrix::solver（第 4 章的 symMatrix/asymMatrix 表）。选择依据是字典里的 solver 名字（PCG/PBiCGStab/GAMG/smoothSolver）。
因为 fvMatrix 就是 lduMatrix，这个转换是"零成本"的——没有把有限体积矩阵翻译成通用稀疏矩阵的开销。这正是 9.1 图中"没有额外转换层"的含义。
9.5.3 松弛与压力参考
- 松弛：
UEqn.relax()依据fvSolution的relaxationFactors（场松弛与方程松弛）缩放对角与源项，抑制迭代震荡。 - 压力参考：
setComponentReference()（fvMatrixSolve.C第 36–57 行）在封闭域求解压力时，把某个单元/分量的系数与源项做个偏置，消除"压力只定到常数"的奇异性：
internalCoeffs_[patchi][facei].component(cmpt) +=
diag()[psi_.mesh().boundary()[patchi].faceCells()[facei]];
boundaryCoeffs_[patchi][facei].component(cmpt) +=
diag()[psi_.mesh().boundary()[patchi].faceCells()[facei]] * value;
9.5.4 隐式组装 lduPrimitiveMeshAssembly
fvMatrix 里有 useImplicit_ 与 lduAssemblyName_，指向 src/finiteVolume/lduPrimitiveMeshAssembly/。这是把"边界的伪矩阵系数"以更快的方式装配进 LDU 矩阵的机制，用于减少边界处理的重复开销。它仍遵循"网格拓扑 → LDU 结构"的基本框架，只是装配路径更优化。
9.6 小结：三层配置 + 一个矩阵
本章可浓缩为一句口诀：
fvSchemes决定"怎么离散"，fvSolution决定"怎么求解"，fvMatrix装下"离散后的方程"，它本身就是lduMatrix，直接喂给求解器。
在脑子里建立这条链，再看任意求解器代码，就会非常顺：
solve(fvm::ddt(rho, U) + fvm::div(phi, U) + turbulence->divDevReff(U)
==
-fvc::grad(p));
fvm::ddt/div/...按fvSchemes选格式、组装fvMatrix；==把fvc::grad(p)挪成源项；solve()按fvSolution选求解器（RTS）迭代求解；- 边界项经由 patch 的四个 coeffs 函数进入矩阵。
下一章深入这两端的细节：离散算子（ddt/grad/div/laplacian）与插值格式。
本章要点回顾
- 路径：PDE →（有限体积离散）
fvMatrix→（继承）lduMatrix→（RTS）Krylov/GAMG 求解器。 fvSchemes（=schemesLookup）管理ddtSchemes/divSchemes/gradSchemes/laplacianSchemes/interpolationSchemes/snGradSchemes/fluxRequired；格式字典按READ_MODIFIED处理，支持运行期改。fvSolution（=solution）管理solvers{}、场/方程的松弛因子（支持Function1）与cacheTemporaryObjects。fvMatrix<Type>继承refCount + lduMatrix，数据含psi_（引用）、source_、internalCoeffs_/boundaryCoeffs_（边界贡献）、faceFluxCorrectionPtr_（非正交修正）。- LDU 存储（
diag/upper/lower）天然匹配有限体积的"面耦合"，矩阵非零结构由网格拓扑唯一决定。 fvm::返回tmp<fvMatrix<Type>>（隐式、进矩阵），fvc::返回tmp<GeometricField<...>>（显式、算好值）；Su/Sp加源项。- 边界经 patch 的 coeffs 函数进入
diag/source；求解走solveSegregatedOrCoupled→lduMatrix::solver（RTS），并支持松弛与压力参考。
