大家在使用BYRFoam的过程中，有时会有些个性化的需求，比如输出一些自定义的物理量。本篇帖子就举两个例子，教大家如何修改源代码，自定义求解器。

关于自定义求解器

首先说明一点，不建议大家直接修改rdeFoam/source/solver/BYRFoam/中的代码。一个好的习惯是把BYRFoam的代码拷贝一份，新建一个求解器，比如叫cBYRFoam，或者你爱叫啥都可以，基于它进行修改。

  1. 首先在solver/目录下运行:

    1
    2
    3
    cp BYRFoam/ -r cBYRFoam/
    cd cBYRFoam 
    mv BYRFoam.C cBYRFoam.C

  2. 然后修改cBYRFoam/Make/目录里的files文件，内容如下：

    1
    2
    3
    cBYRFoam.C

    EXE = ${FOAM_USER_APPBIN}/cBYRFoam

  3. 之后再对主程序进行正式的修改。

计算并输出 Pmax

前言

如果你想得到爆轰波的数值胞格结构，可以记录一段时间 $[t_0, t_1]$ 内，流场中每个网格点 $(x,y,z)$ 的压强最大值（记为 $P_{max}$ ）:

$$
P_{max} = \max \limits_{t_0 \leq t \leq t_1} P(x,y,z,t)
\tag{1.1}
$$

需要注意它与某一时刻全流场中的压强最大值（记为 $maxP$ ）的区别:

$$
maxP = \max \limits_{(x,y,z) \in \Omega} P(x,y,z,t)
\tag{1.2}
$$

从两者的表达式可以明显看出区别，这里不再赘述。

需要做的修改

  1. 打开createFields.H文件，对照原程序，在下面位置修改代码:

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    volScalarField& he = thermo.he();
    volScalarField& p = thermo.p();

    // new add to create Pmax
    volScalarField Pmax
    (
        IOobject
        (
            "Pmax",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        thermo.p()
    );

  2. 打开cBYRFoam.C文件，对照原程序，在下面位置修改代码:

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    gamma = thermo.gamma();
    c = sqrt(gamma/psi);

    // update Pmax for every mesh point to obtain cellular stucture
    Pmax = max(Pmax,p); 

    runTime.write();

    // Monitor max pressure in flow field at this time
    // Info<< "maxP = " << max(p).value() << endl;

  3. 运行wmake编译cBYRFoam即可。

运行求解器

在算例目录运行cBYRFoam，计算完成后，每个时间步的文件夹里就会有Pmax这个文件。

计算并输出 Qdot

前言

在带多组分化学反应的守恒型 $Euler$ 方程组中，能量方程和组分方程分别为：

$$
\frac{\partial (\rho E)}{\partial t} + \nabla \cdot (\rho \vec u H) = 0
\tag{2.1}
\label{Energy Eq.1}
$$

$$
\frac{\partial (\rho Y_i)}{\partial t} + \nabla \cdot (\rho \vec u Y_i) = \dot \omega _i
\tag{2.2}
$$

其中，$E$ 为单位质量的总内能，$H$ 为单位质量的总焓：

$$
H = E + p/\rho = \sum_{i=1}^N Y_i h_i + u^2/2
\tag{2.3}
$$

我们可以看到，当以总内能 $E$ 为能量变量时，能量方程 $\eqref{Energy Eq.1}$ 是没有源项的。
对于第 $i$ 个组分的比焓 $h_i$，其是温度 $T$ 的函数，定义为：

$$ h_i(T) = h_{0,i} + h_{s,i}(T) = h_{0,i} + \int_{T_{ref}}^{T} C_{p,i} \,\mathrm{d}T \tag{2.4} $$

其中，$h_{0,i}$ 为标准生成焓，定义为：

标准状态下，由自然状态存在的单质通过定压定温反应生成某物质的过程中焓的增量。

对于给定的组分，其标准生成焓 $h_{0,i} = h_{i}(T_{ref})$ 是一个常数，而显焓 $h_{s,i}$ 与温度有关。
我们定义总显焓 $H_s$ 和总显内能 $E_s$ 如下:

$$
\sum_{i=1}^{N} Y_i h_{0,i} = H - H_s = E - E_s
\tag{2.5}
$$

那么能量方程 $\eqref{Energy Eq.1}$ 可变为：

$$
\frac{\partial (\rho E_s)}{\partial t} + \nabla \cdot (\rho \vec u H_s) =
-\sum_{i=1}^N h_{0,i} \left[\frac{\partial (\rho Y_i)}{\partial t} + \nabla \cdot (\rho \vec u Y_i) \right] = -\sum_{i=1}^N h_{0,i} \dot \omega _i
\tag{2.6}
$$

这就是BYRFoam中无粘的能量方程的形式，右边源项代表化学反应放热率（记为 Qdot ），单位 $J/(m^3 \cdot s)$ 。

需要做的修改

  1. 打开createFields.H文件，对照原程序，在下面位置修改代码:

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    Info << "Yi size: " << Y.size() << endl;

    // new add to create Qdot
    volScalarField Qdot
    (
        IOobject
        (
            "Qdot",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        mesh,
        dimensionedScalar("Qdot", dimEnergy/dimTime/dimVolume, 0.0)
    );

  2. 打开cBYRFoam.C文件，对照原程序，在下面位置修改代码:

    1
    2
    3
    4
    5
    6
    7
    8
    // --- Solve energy
    Qdot = reaction->Qdot(); // new add to update Qdot
    solve
    (
        fvm::ddt(rhoE) + fvc::div(EPhi) 
      - fvc::div(sigmaDotU) // this term equals zero in inviscid simulation
     == Qdot
    );

  3. 运行wmake编译cBYRFoam即可。

运行求解器

在算例目录运行cBYRFoam，计算完成后，每个时间步的文件夹里就会有Qdot这个文件。