在本节中，通过在 interFoam 求解器中添加一个新的 PDE 来实现一个新的求解器，以模拟两个不可压缩的不混溶流体相。 PDE 及其支持代码的目的是展示如何执行实现，而不是模拟真实的物理传输过程。流体界面上和跨流体界面的传输现象建模需要仔细和严格的推导。在大多数情况下，该等式源于 compressibleInterFoam 求解器中可用的实现。在那里，由于使用了耦合压力、温度和密度的状态方程，因此考虑了跨界面的热传递。为简单起见，热力学状态方程被省略，假设为层流。

9.3.1 附加模型方程

作为层流两相流中非定常被动标量输运模型的示例，使用以下等式 $\frac{\partial ρT}{\partial t} + \nabla \cdot (ρ U T) - \nabla \cdot (D_{e ff} \nabla T) = 0 (9.1)$ 其中 $T$ 表示温度， $ρ$ 是密度， $U$ 是速度， $D_{e ff}$ 是相混合物的导热系数。

interFoam 求解器在单场公式中实现两相 Navier-Stokes 方程，以模拟两个不可混溶的不可压缩流体相的两相流。该模型通过引入额外的标量（体积分数场）来区分流体相，将两种不混溶流体相的流动描述为单个连续体的流动。关于两相流建模的更多信息可以在 [8] 的书中找到。

任何具有介于两者之间的 α 值的单元都被视为界面单元，其材料特性加权为两个主要相的混合物，因为 α 定义为体积分数 $α = \frac{V _{1}}{V} (9.2)$ 其中 V1 是总体积为 V 的单元中phase 1 所占的体积。然后使用体积分数场对单个连续体的特性进行建模，形成混合物量，例如混合物粘度 $ν = ν_{1} α_{1} + (1 - α_{1}) ν_{2} (9.3)$ 或混合密度: $ρ = ρ_{1} α_{1} + (1 - α_{1}) ρ_{2} (9.4)$ 其中 ν1, ν2 和 ρ1, ρ2 分别是两个流体相的运动粘度和密度。方程 9.1中Deff 中的有效热传导系数在本例中以类似方式建模，使用 α1 场 $D_{e ff} = \frac{α k _{1}}{C _{v 1}} + \frac{( 1 - α ) k _{2}}{C _{v 2}} (9.5)$ 这里， $k_{1} 、 k_{2}$ 和 $C_{v 1} 、 C_{v 2}$ 代表各个流体相的传导系数和比热容。

由方程 9.5 建模的热传导系数基于体积分数场，其方式与其他相属性类似。以这种方式使用体积分数物理场值，为此示例将物理属性的恒定值分配给两个流体相的主体区域。流体界面的区域将具有区间 [0, 1] 内的 α1 值，因此将定义相应相属性的两个恒定值之间的过渡区域。过渡区轮廓的形状将由描述混合特性的模型的性质决定，如传导系数（见方程 9.5）。这种方法如何准确地应用于跨移动界面的热传递的实际物理问题对于描述如何在 OpenFOAM 中执行求解器修改并不重要，因此不在范围之内。

9.3.2 准备求解器进行修改

在更改现有求解器的源代码之前，必须创建求解器应用程序的新副本。将 interFoam 求解器的源代码目录复制到个人应用程序目录并重命名：

cp -r $FOAM_SOLVERS/multiphase/interFoam/ $FOAM_RUN/../applications/
cd $FOAM_RUN/../applications/
mv interFoam heatTransferTwoPhaseSolver
cd heatTransferTwoPhaseSolver
mv interFoam.C heatTransferTwoPhaseSolver.C

提示：每个 OpenFOAM 版本都为用户工作目录导出 $FOAM_RUN 变量，但在 OpenFOAM 安装期间不会生成它。如果此目录不可用，可以使用 mkdir -p $FOAM_RUN 创建它。

在接下来的步骤中，必须删除所有不需要修改的文件。这样做是为了防止代码重复，应始终避免。如果代码重复，那么维护求解器就会出现问题，因为随着 OpenFOAM 的发展，求解器文件需要更新。这就是为什么在实施求解器时对现有文件的更改应保持在最低限度。

提示：在开发新的求解器时，将共享文件中引入的更改保持在最低限度。在重用现有文件的同时引入其他文件并将它们包含到新的求解器应用程序中，使新求解器更易于维护，因为重用现有文件将在新的 OpenFOAM 版本中获取这些文件中的更改。

应该从目录中清除所有未更改的文件：

rm -rf interMixingFoam/ overInterDyMFoam/
rm alphaSuSp.H createFields.H correctPhi.H initCorrectPhi.H rhofs.H pEqn.H UEqn.H

提示：如果从一开始就知道必须修改哪些求解器文件，则只需复制这些文件和 Make 目录。

由于求解器源代码文件已重命名，因此必须修改 Make/files 以反映更改。 Make/files 配置文件应仅包含此处显示的行：

heatTransferTwoPhaseSolver.C
EXE = $(FOAM_USER_APPBIN)/heatTransferTwoPhaseSolver

让 wmake 构建系统知道要从重命名的求解器源代码文件中构建并安装在用户应用程序二进制目录中的可执行应用程序。因此，这将指示 wmake 编译文件 heatTransferTwoPhaseSolver.C 及其依赖项，并创建一个名为 heatTransferTwoPhaseSolver 的可执行求解器。此时调用 wmake 会失败：我们已经删除了方程文件、场初始化等，最终得到

ls
heatTransferTwoPhaseSolver.C Make

为了通知新求解器重新使用来自 interFoam 的文件和来自 $FO A M_{S} O L V ERS / m u lt i p ha se / V o F 的共享文件， M ak e / o pt i o n s 必须修改为 ‘‘‘ ba s h I NC = - I$ (FOAM_SOLVERS)/multiphase/interFoam
-I $(FO A M_{S} O L V ERS) / m u lt i p ha se / V o F ...‘‘‘ 这样， w mak e 被指示首先查看 in t er F o am 求解器的文件夹以查找与 in t er F o am 共享的文件，然后查看 ‘$ FOAM_SOLVERS/multiphase/VoF `以查找与 VoF 求解器系列共享的文件。

此时，可以使用 wmake 构建新的 heatTrasnferTwoPhaseSolver，它的行为与 interFoam 完全相同。

9.3.3 添加温度场和热传导系数

必须从transportProperties字典中为每个相查找两个新的所需材料相属性：热传导系数$k$和比热容 $C_v$。此外，需要初始化新的温度 $T$ 和有效传热系数 $D_{eff} $。初始化被放置在一个新的初始化文件 createHeatTransferFields.H 中，连同温度和传热系数的初始化，如下面的代码段所示。代码段所示的新 createHeatTransferFields.H 应插入到新求解器中 createFields 之后.H 文件，在 immiscibleIncompressibleTwoPhaseMixture 被初始化之后。类型 dimensionedScalar 的热容量将从 transportProperties 字典文件中每个阶段的子字典中加载。完成所有字典查找和声明后，新模型方程 9.1 将被添加到求解算法中。

const dictionary& phase1dict = mixture.subDict ("phase1");
const dictionary& phase2dict = mixture.subDict ("phase2");
auto k1 (phase1dict.get<dimensionedScalar> ("k"));
auto Cv1 (phase1dict.get<dimensionedScalar> ("Cv"));
auto k2 (phase2dict.get<dimensionedScalar> ("k"));
auto Cv2 (phase2dict.get<dimensionedScalar> ("Cv"));
volScalarField T
(
    IOobject
    (
        "T",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);
volScalarField Deff
(
    "Deff",
    (alpha1*k1/Cv1 + (1.0 - alpha1)*k2/Cv2)
);

9.3.4 编写温度方程

温度方程放置在单独的文件 TEqn.H 中，然后该文件包含在求解器应用程序中。 TEqn.H 文件中定义的代码实现了两个操作。它计算系数 Deff 并求解方程 9.1 给出的被动温度传递 PDE。

Deff 系数计算为基于单元中心 alpha1 值（公式 9.5）的线性加权平均值。对于这个例子，计算 Deff 系数的方法与 compressibleInterFoam 求解器中使用的方法相同。由于 Deff 系数会因单元而异，因此它被实现为 volScalarField，而不是单个标量。在 createHeatTransferFields.H 中初始化扩散系数字段后，可以如下所示实现传输方程。

Deff == alpha1*k1/Cv1 + (1.0 - alpha1)*k2/Cv2);
solve
(
    fvm::ddt(rho, T)
    + fvm::div(rhoPhi, T)
    - fvm::laplacian(Deff, T)
);

现在 TEqn.H 已完成，需要将其添加到在主求解器应用程序文件 heatTransferTwoPhaseSolver.C 中实现的求解算法中. TEqn.H 应该插入到内部 PIMPLE 循环的上方，与下面代码段所示的动量方程相邻。

while (pimple.loop())
{
    ...
    mixture.correct();
    if (pimple.frozenFlow())
    {
        continue;
    }
    #include "TEqn.H"
    #include "UEqn.H"
    ...
}

这样就完成了将热传递方程添加到 interFoam 求解器的求解算法所需的源代码修改。应该从构建过程生成的旧文件中清理求解器目录，然后编译求解器应用程序：

wclean
wmake

完成新求解器后，需要设置与新求解器兼容的模拟案例。

9.3.5 设置案例

提示：最终案例设置可以在案例存储库中找到，位于 chapter09/2DheatXferTest 文件夹中。

创建新求解器后，现在需要新的初始条件、边界条件、材料属性和求解器控制输入参数来处理传热计算。在0文件夹下，温度场以T文件的形式添加，其必须具有合适的量纲：

FoamFile
{
    version        2.0;
    format         ascii;
    class         volScalarField;
    location    "0";
    object         T;
}
dimensions [0 0 0 1 0 0 0];

应该设置边界条件，在这个例子中我们简单地使用零梯度边界条件。要为新的温度场设置初始条件，需要编辑系统目录下的 setFieldsDict 文件。在本例中，我们将上升气泡的温度设置为高于其环境的值：

defaultFieldValues
(
    volScalarFieldValue alpha1 0
    volVectorFieldValue U (0 0 0)
    volScalarFieldValue T 300
);
regions
(
    sphereToCell
    {
        centre (1 1 0);
        radius 0.2;
        fieldValues
        (
            volScalarFieldValue alpha1 1
            volScalarFieldValue T 500
        );
    }
);

还应修改constant目录中的 transportProperties 文件，为每个相添加 k 和 Cv 系数值。我们任意选择了这些值：

e1
transportModel Newtonian;
nu        nu         [ 0 2 -1 0 0 0 0 ] 1e-02;
rho        rho         [ 1 -3 0 0 0 0 0 ] 100;
k        k         [ 1 1 -3 -1 0 0 0] 100;
...

e2
transportModel Newtonian;
nu        nu          [ 0 2 -1 0 0 0 0 ] 1e-02;
rho        rho         [ 1 -3 0 0 0 0 0 ] 100;
k         k         [ 1 1 -3 -1 0 0 0] 100;
Cv        cv         [ 0 2 -2 -1 0 0 0] 100;
...

由于向求解器添加了新的 PDE，因此应在 system/fvSchemes 中为其微分算子选择离散化方案，以防未提供或不适用默认选项。在此示例中，默认 CDS（高斯线性）方案将替换为迎风方案：

divSchemes
{
    default                Gauss linear;
    div(rho*phi,U)         Gauss limitedLinearV 1;
    div(rho*phi,T)         Gauss upwind;

最后，应设置涉及温度方程离散化得到的线性系统解的参数。将TFinal 添加到 fvSolution 文件中，例如

TFinal
{
    solver                PBiCG;
    preconditioner          DILU;
    tolerance            1e-06;
    relTol                0;
}

这些是设置传热方程求解器及其参数的选项。这应该是运行带有新模型方程的新求解器所需的最后一个案例配置更改。求解器可以在案例目录中运行，并且可以分析模拟的结果。

9.3.6 执行求解器

为了运行案例，请确保示例代码库的库和应用程序代码通过调用自动编译:

./Allwmake

在示例代码存储库的主目录中。要启动求解器，只需在仿真案例目录中执行以下命令：

heatTrasferTwoPhaseSolver

图 9.1 显示了用于模拟二维上升气泡的体积分数场 α 和温度场 T 的分布。气泡温度设置为高于周围流体的温度，这将导致气泡在上升时向其周围释放热通量。

不同几何形状的热传导存在精确解，例如实心圆柱体或板。为一个被热静止空气包围的圆柱体或两个接触的板创建一个模拟案例，并验证新的求解器。解是否收敛到精确解？

图 9.1 2D 上升气泡模拟情况下两个选定模拟时间的温度分布。