buoyantBoussinesqSimpleFoam代码解析及使用心得
Boussinesq
在对有浮升力存在的换热问题进行求解时，常常使用Boussinesq假设将问题简化，在保持不可压求解相对便利性的同时更在一定程度上缓解方程组的非线性。随着计算机算力的不断提高，基于可压流对浮升力作用下的换热问题进行求解已不再昂贵。然而，👉🏼一方面，在很多领域，Boussinesq假设的沿用久经考验，成为业界传统或被写入手册（Best Practice Guidelines），因而并没有移花接木的必要；👉🏼另一方面，在开展一些新工作时，需要在等温流动计算结果的基础上增加对换热问题的考虑，而前期研究很有可能是建立在对不可压流动的求解上的，若为了考虑密度差带来的影响而将问题全盘转移则会大费周章。
为了较快得到结果，稳态解一般是首先考虑的，然而浮力流动往往并不呈现显著的稳态特性，在很多问题中很难通过稳态求解达到收敛，虽然Boussinesq假设可以在一定程度上改善问题的收敛性，不过对于该类问题，即使是收敛解也很有可能存在误导（misleading）。在此类问题的稳态求解中，经常有为了使问题收敛而不得不使用并不符合物理规律（non-physical）的边界条件来伺候求解器的无奈境况。很多时候使用自认为最接近物理实际的边界条件却造成发散；而换成其他看似离题的边界条件却能让问题很好地收敛并与实验数据有较好的吻合，这种感受着实一言难尽😂。尽管过程充满 狗血，我们还是得迎难而上，毕竟对该类问题进行稳态求解仍然是常规的研究思路。
本期内容着眼于对有浮升力存在时换热问题的稳态求解在OpenFOAM中的实现，围绕buoyantBoussinesqSimpleFoam求解器，从控制方程、迭代思路、代码实现和求解计算四个方面入手。
👇🏼
1
控制方程
首先看Boussinesq假设，这一假设包含三部分：
1 | 流体粘性耗散忽略不计，因而能量方程不含由此带来的“耗散函数”源项Φ(dissipation function)； |
2 | 除密度外其它物性为常数； |
3 | 密度的影响只在动量方程中与体积力有关的项中出现，其它项中仍为常数。 |
以冷流体温度为参考温度Τref（Reference Temperature [K]），其密度为参考密度ρ0[kg/m3]，同时假设密度随温度线性变化，有
方程1
ρ=ρ0[1-β(Τ-Τref)]
其中，β为体胀系数（Thermal Expansion Coefficient [K-1]）。对于理想气体，β值一般取参考温度的倒数，即β≈1/Τref，所以在室温下，空气的体胀系数可取3.33×10-3，对于15K以内的温差，计算误差小于1%（Ferziger and Peric 2001, p.15)。对于其他流体，更一般的定义为：
基于Boussinesq假设的不可压连续性方程、动量方程、与能量方程如下（爱因斯坦运算符）：
方程2
方程3
方程4
根据需要还可添加湍流模型，此处省略。
2
迭代思路
buoyantBoussinesqSimpleFoam使用大家耳熟能详的“压力耦合方程组的半隐式方法“，即SIMPLE算法（Semi-Implicit Method for Pressure Linked Equations）进行迭代求解。SIMPLE算法具体流程不在此详述，而在随后的代码中注释。提前说明的是：
速度修正值方程由动量离散方程推得，并略去了邻点速度修正值
的影响，这意味着每一个格点的速度修正值
视作与压力修正值
有关，格点间的u’不相制约，因而无需建立代数方程组，这就是所谓的“显式求解”。而对于压力修正值方程，相邻格点的影响不作忽略，格点间p’相制约，需要通过建立代数方程组进行整场求解，所谓“隐式求解”，这对于动量离散方程中u*也是一样。速度修正值方程显式求解，动量离散方程和压力修正值方程隐式求解，这就是所谓的“半隐”，将在代码中体现。该求解器采用了压力变换（Pressure Shift）的思路，可在压力修正值方程中将浮升力的影响同样纳入考虑，介绍如下。
在无需采用Boussinesq假设的可压动量方程中，等号右边为（注意此处P为大写，单位N/m2）：
式5
引入参考压力Prgh，
方程6
其中href,k代表压力参考点的坐标，令位差hk=(x-href)k
将方程6回代到式5，涉及对如下压力梯度的计算：
方程7
整理后，得到，
式8
相应的，对于采用Boussinesq假设的不可压动量方程，方程3变形为，
方程9
其中为ρk为折算密度（effective kinematic density）
方程10
动量方程的半离散化，形式为：
方程11
两边同除，有
方程12
其中的源项为，
方程13
修正后的速度up也满足动量方程，由
相应的速度修正值
，
方程14
对求散度并令其为0（满足连续性方程），并结合
得到
方程15
按SIMPLE的半隐思想，略去上式中相邻节点速度修正值的影响（第二项），得到最终的压力修正值，
方程16
下面将结合OpenFOAM-6的代码进行解析。
求解器源码位于$FOAM_SOLVERS/heatTransfer/buoyantBoussinesqSimpleFoam目录。包含的文件如下：
文件名 | 功能 |
buoyantBoussinesqSimpleFoam.C | 主程序 |
creatFields.H | 初始化各物理量场 |
readTransportProperties.H | 读取物性参数 |
UEqn.H | 解动量离散方程 |
TEqn.H | 解温度离散方程 |
pEqn.H | 解压力修正值方程和速度修正值方程 |
ℹ以下大量代码，可左右滑动查看 | |
buoyantBoussinesqSimpleFoam.C
#include"fvCFD.H"//有限体积法头文件
#include"singlePhaseTransportModel.H"//声明singlePhaseTransportModel类，基于粘性的单相流动输运模型
#include"turbulentTransportModel.H"//湍流模型的命名空间
#include"noRadiation.H"//声明noRadiation类，“不考虑辐射的辐射模型”
#include"fvOptions.H"//声明option类，提供一系列成员函数用于执行对源项的操作，例如T方程中即将用到的constrain()和correct()
#include"simpleControl.H"//声明simpleControl类，提供一系列成员函数用于执行simple循环指令,例如下面的loop()
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
intmain(int argc, char *argv[])
{
#include"postProcess.H"//声明后处理函数，以便实时调用
#include"setRootCaseLists.H"//检查case文件结构
#include"createTime.H"//根据controlDict中的设置，确认当前时间步，并定义runTime对象
#include"createMesh.H"//定义前述fvMesh类的对象mesh
#include"createControl.H"//定义前述simpleControl类的对象simple
#include"createFields.H"//初始化各物理量场，见下一节
#include"initContinuityErrs.H"//定义并初始化连续性误差
turbulence->validate(); //在开始求解前，确认已正确计算湍流粘度，此为部分湍流模型的功能
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (simple.loop(runTime)) //执行simple循环
{
Info<< "Time = " << runTime.timeName() << nl << endl;
{
#include"UEqn.H"//解动量离散方程
#include"TEqn.H"//解温度离散方程
#include"pEqn.H"//解压力修正值方程和速度修正值方程
}
laminarTransport.correct(); //修正运动粘度nu
turbulence->correct(); //解湍流变量离散方程
runTime.write(); //写数据到磁盘
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;
return0;
}
creatFields.H
Info<< "Reading thermophysical properties\n" << endl;
Info<< "Reading field T\n" << endl;
volScalarField T //通过volScalarField类的构造函数，定义一个对象T
( //该构造函数有两个参数
IOobject //一个是对IOobject构造函数的调用
( //该构造函数有五个参数
"T", //变量名称
runTime.timeName(), //由runTime对象的成员函数timeName()获取当前时间步，与时间关联
mesh, //网格对象mesh，与空间关联
IOobject::MUST_READ, //读写属性1：必读，因而需要用户在case中提供T文件
IOobject::AUTO_WRITE //读写属性2：自动写
),
mesh //另一个参数是对象mesh，通过mesh对象间接读取数据进行初始化
);
Info<< "Reading field p_rgh\n" << endl;
volScalarField p_rgh //类似的，定义标量场p_rgh，必读，自动写，需要用户提供p文件
(
IOobject
(
"p_rgh",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading field U\n" << endl;
volVectorField U //类似的，定义向量场U，必读，自动写，需要用户提供U文件
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
#include"createPhi.H"//调用头文件createPhi.H
/* 头文件中的语句为
surfaceScalarField phi //定义面向量场phi, 即面通量，自动写，无需要用户提供phi文件
(
IOobject
(
"phi",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
fvc::flux(U)
); */
#include"readTransportProperties.H"//调用头文件 readTransportProperties.H
//其中使用singlePhaseTransportModel类的构造函数定义了laminarTransport对象，
Info<< "Creating turbulence model\n" << endl;
autoPtr<incompressible::turbulenceModel> turbulence //使用智能指针实例化一个对象turbulence
( //用incompressible::turbulenceModel类的构造函数，进行初始化
incompressible::turbulenceModel::New(U, phi, laminarTransport)
);
// Kinematic density for buoyancy force
volScalarField rhok //定义标量场rhok，默认读写属性，无需提供相应文件
(
IOobject
(
"rhok",
runTime.timeName(),
mesh
),
1.0 - beta*(T - TRef) //rhok值由表达式（10）进行更新
);
// kinematic turbulent thermal thermal conductivity m2/s
Info<< "Reading field alphat\n" << endl;
volScalarField alphat //定义标量场alphat，必读，自动写，需用户提供alphat文件
(
IOobject
(
"alphat",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
#include"readGravitationalAcceleration.H"//从case文件中读取指定的重力加速度
#include"readhRef.H"//调用头文件readhRef.H
/*执行的语句为
Info<< "\nReading hRef" << endl;
uniformDimensionedScalarField hRef //读取一个均一的带单位的标量场hRef，同式(6)
(
IOobject
(
"hRef",
runTime.constant(), //不随时间变化
mesh,
IOobject::READ_IF_PRESENT, //可指定，也可不给
IOobject::NO_WRITE //不会自动生成文件
),
dimensionedScalar(dimLength, 0) //传一个值为0，单位为长度单位的标量作初始化
);
*/
#include"gh.H"//调用头文件gh.H
/*执行的语句为
dimensionedScalar ghRef(- mag(g)*hRef); // 定义标量场 ghRef = -|g|×hRef
volScalarField gh("gh", (g & mesh.C()) - ghRef); //gh 见式(7)，体心值，
surfaceScalarField ghf("ghf", (g & mesh.Cf()) - ghRef); //gh 见式（7），面心值
// & 为重载的运算符，代表点乘，mesh.C()为体心位置矢，mesh.Cf()为面心位置矢
*/
volScalarField p //标量场p,不需要提供p文件
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
p_rgh + rhok*gh //用p_rgh, rhok和gh进行初始化
);
label pRefCell = 0; //默认压力参照点
scalar pRefValue = 0.0;
setRefCell
(
p,
p_rgh,
simple.dict(),
pRefCell,
pRefValue
);
if (p_rgh.needReference())
{
p += dimensionedScalar //整体加减，以使得默认参照点处的压力等于默认参考压力
(
"p",
p.dimensions(),
pRefValue - getRefCellValue(p, pRefCell)
);
}
mesh.setFluxRequired(p_rgh.name());
#include"createIncompressibleRadiationModel.H"//辐射模型
#include"createFvOptions.H"//定义option类的对象fvOptions
readTransportProperties.H
singlePhaseTransportModel laminarTransport(U, phi); //定义了前述singlePhaseTransportModel类的对象laminarTransport
// Thermal expansion coefficient [1/K]
dimensionedScalar beta //定义一个有单位的标量beta
(
"beta",
dimless/dimTemperature, //单位为 无量纲/温度单位 即1/K
laminarTransport //从上述laminarTransport对象中读取，laminarTransport会从transportProperties字典中读取
);
// Reference temperature [K]
dimensionedScalar TRef("TRef", dimTemperature, laminarTransport);
// Laminar Prandtl number
dimensionedScalar Pr("Pr", dimless, laminarTransport);
// Turbulent Prandtl number
dimensionedScalar Prt("Prt", dimless, laminarTransport);
UEqn.H
tmp<fvVectorMatrix> UEqn
(
fvm::div(phi, U) //对速度求散度(隐式加入等号左侧系数矩阵)
+ turbulence->divDevReff(U) //见式（3）第二项
//对偏应力(deviatoric stress)求散度（divergence）
//其中的R代表雷诺时均，eff代表总的扩散系数,
//即运动粘度nu及湍流扩散系数nu_t的和nu_eff，
==
fvOptions(U) //源项
);
UEqn.relax(); //速度亚松弛
fvOptions.constrain(UEqn); //通过源项实现特定区域速度恒定(若指定)
if (simple.momentumPredictor()) //若fvSolution文件中SIMPLE下的momentumPredictor设为YES
{
solve
(
UEqn
==
fvc::reconstruct //由以下的面心值得到体心值，显式加入等号右侧
(
(
- ghf*fvc::snGrad(rhok) //ghf代表gh的面心值，再与 “折算密度” rhok的梯度相乘
- fvc::snGrad(p_rgh) //p_rgh的梯度
)*mesh.magSf() //乘以面积 同式（9）
)
);
fvOptions.correct(U); //通过源项实现速度限值修正(若指定)
}
TEqn.H
{
alphat = turbulence->nut()/Prt; //求湍流热扩散系数 alpha_t = nut / Prt
alphat.correctBoundaryConditions(); //更新边界
volScalarField alphaEff("alphaEff", turbulence->nu()/Pr + alphat); //求总热扩散系数
fvScalarMatrix TEqn //温度方程同(4)
(
fvm::div(phi, T)
- fvm::laplacian(alphaEff, T)
==
radiation->ST(rhoCpRef, T) //实际上没有考虑辐射，因为头文件中默认选取的辐射模型为无
+ fvOptions(T) //源项
);
TEqn.relax(); //温度亚松弛
fvOptions.constrain(TEqn); //通过源项实现指定区域温度恒定(若设置)
TEqn.solve(); //解温度方程
radiation->correct(); //实际上什么都不做
fvOptions.correct(T); //通过源项实现温度限值修正(若指定)
rhok = 1.0 - beta*(T - TRef); //根据新的温度场更新rhok
}
pEqn.H
{
volScalarField rAU("rAU", 1.0/UEqn.A()); // rAU = 1/A
surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU)); //rAUf = (1/A)_f
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh)); // HbyA = H / A 见(12)
UEqn.clear();
surfaceScalarField phig(-rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf());
//phig = b/A = - gh * grad(rhok) / A 见（12）
surfaceScalarField phiHbyA
(
"phiHbyA",
fvc::flux(HbyA)
); //phiHbyA = (HbyA)_f 见（16）
adjustPhi(phiHbyA, U, p_rgh);
phiHbyA += phig; // 把 b/A 并入 H/A
// Update the pressure BCs to ensure flux consistency
constrainPressure(p_rgh, U, phiHbyA, rAUf);
while (simple.correctNonOrthogonal())
{
fvScalarMatrix p_rghEqn
(
fvm::laplacian(rAUf, p_rgh)== fvc::div(phiHbyA) //参考压力修正值方程 同（16）
);
p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
p_rghEqn.solve(); //求解 参考压力修正值方程
if (simple.finalNonOrthogonalIter())
{
// 计算面通量
phi = phiHbyA - p_rghEqn.flux();
p_rgh.relax(); //压力方程亚松弛
// 速度修正值方程，可见下式中速度由参考压力直接得到，没有求解线性方程组，所谓显式修正
U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
U.correctBoundaryConditions();
fvOptions.correct(U);
}
}
#include"continuityErrs.H"
p = p_rgh + rhok*gh; //更新p
if (p_rgh.needReference())
{
p += dimensionedScalar
(
"p",
p.dimensions(),
pRefValue - getRefCellValue(p, pRefCell)
);
p_rgh = p - rhok*gh;
}
}
边界条件
OpenFOAM安装文件中自带了两个算例，位于$FOAM_TUTORIALS/heatTransfer/buoyantBoussinesqSimpleFoam/路径。
遗憾的是这两个算例都为密闭空间自然对流，而对于有开口存在的问题，提供以下参考。
一般对于单纯的速度-压力耦合问题，OpenFOAM对边界条件组合的评价如下：
inlet | outlet | Stability |
flowRateVelocity for U | fixedValue for p | Exellent |
totalPressure for p | totalPressure for p | Very Good |
totalPressure for p | fixedValue for p | Good |
fixedValue for p | fixedValue for p | Poor |
以上都为第一类边界条件。 对于不可压流动，由于速度和压力的强耦合，在同一边界上，速度和压力不能同时都给第一类边界条件。对于出口边界，通常都是给压力设置第一类边界条件，而对速度给第二类边界条件，具体选取则涉及对回流的考虑。 OpenFOAM对出口边界条件的选取给予如下建议： | ||
Type | Condition | Stability |
Outflow | zeroGradient for U | Unstable if flow reveres |
Blocked | inletOutlet for U | Good, but unphysical |
Return flow I | pressureInletOutleVelicity for U fixedValue for p | Good |
Return flow II | pressureInletOutleVelocity for U totalPressure for p | Very Good |
| ||
而在具体应用时，这些边界条件没有绝对的优劣之分，越是复杂的问题越是视阴晴而定。举例来说，在一些回流存在但并不占优的问题中巧妙地使用充分发展假设可以让问题简化而达到收敛，而有时充分照顾到回流求解反而不收敛，所谓给脸不要敛😑。
对于有浮升力存在的换热问题，问题显得稍有不同，因为采用了压力变换的思路，引出了prgh，与速度的耦合关系有相应的变动。研读OpenFOAM自带的算例可以看到，对prgh 采用的是fixedFluxPressure边界，所求得的压力梯度满足动量方程，有助于计算的稳定。很多的反馈表明在使用其他出口边界条件都会发散的情况下，使用fixedFluxPressure边界能使问题更好地收敛。这一边界条件同样适用于壁面。
对于温度场，需要指定T和αt的边界条件，在流域内部，αt由湍流模型算得，在壁面则需要指定壁面函数。
边界类型 | T | alphat |
入口 | fixedValue | calculated |
出口 | totalTemperature | calculated |
壁面 | 按需选择fixedValue/zeroGradient/fixedGradient | alphatJayatillekeWallFunction |
| ||
1 | 如果需要指定局部热源：通过topoSet命令选择相关区域，然后在fvOptions中设置scalarSemiImplicitSource |
2 | 如果需要指定局部恒温：通过topoSet命令选择相关区域，然后在fvOptions中设置fixedTemperatureConstraint |
3 | 如果需要设置温度上下限：limitTemperature |
⚠值得注意的是，从热力学的角度出发，对于一个合理的稳态解，能量进出应是匹配的，这在设置边界条件时应该给予考虑。如果是有内热源的问题，对于闭口系可设置与热源强度相当的热流边界，对于开口系也要考虑热力平衡。
差分格式
在求解过程中，会有各种问题造成求解发散。除了结合物理实际使用合适的边界条件外，还应使用合适的差分格式和松弛因子以保证稳定性。
👉🏼对于差分格式，CFD理论告诉我们高阶格式会带来over-shoot，这在温度场的计算中有直观的体现，温度常常溢出应有的范围；有时如果初始场设置不合理，使用高阶格式也会让求解在一开始就发散。在不妥协阶数的前提下，可以把二阶迎风换为带限的差分格式如limitedLinear和vanLeer看是否有好转，进一步的妥协是采用不会造成over-shoot的一阶迎风(upwind)，这一格式具有假扩散的特性，精度低，但能够帮助排查问题究竟是不是出在差分格式上，如果还是不稳定甚至发散，则检查边界条件、松弛因子以及思考对应的物理问题是否能发展到稳态。
👉🏼对于难伺候的问题，可以先算得等温流动的初始流场，而后设置所需的温度边界条件，用一阶迎风让温度场参与计算，残差稳定或收敛后再对温度使用limitedLinear或vanLeer格式，并使用limitTemperature源项对温度范围进行限制。
松弛因子
对于松弛因子，SIMPLE算法的修正思路明确对速度和压力的亚松弛处理提出要求；同时，由于流场与温度场的耦合，在循环体内，ρk 是先使用速度场的预测值计算温度场得到的，而后再进入压力方程参与计算，所以在温度场求解中也应做亚松弛处理。对于湍流模型也是类似。在实践中，可以先从这样的松弛因子开始：
u | p_rgh | T | turbulence parameters |
0.7 | 0.3 | 0.7 | 0.7 |
如果发现某个量不稳定或发散，则可适当调小它的松弛因子。松弛因子只是保证求解器的稳定性，但并不能让由不合理的边界条件或物理描述带来的发散得到解决，更多时候还是要多从其他设置上找原因。
5
总结
本文介绍了OpenFOAM中用于求解有浮升力存在时不可压流换热问题的求解器buoyantBoussinesqFoam，列出了理论推导的关键步骤，介绍了SIMPLE算法的迭代思路，对程序的大部分语句给出了注释及理论对照，最后分享了一些求解器使用上的经验和技巧，希望对大家有所帮助。由于本人水平和时间有限，难免有错漏荒唐之处，还请不吝指正。
