编写自己的求解器

OpenFOAM

---

创建于:20200621

---

0.0 目的

如果不写自己的求解器，那么学这个软件将毫无意义。而写自己的求解器，常见的，无外乎几种目的：

• 不是通用方程。比如多孔介质流动，里面的速度压力方程其实是个能量方程（darcy可由伯努利方程导出）；
• 系数不是常数，而是根据变量变化的。比如超临界二氧化碳换热，传热系数、比热等跟随温度变化明显；
• 存在源项。比如蒸发沸腾，相函数存在源项，而且源项是温度的函数；
• 要额外计算一些通用求解器不知道的量，比如气固流动的沉积问题，沉积量可能需要单独计算；

这些需求里面，其实除了第一条是需要改动求解器框子的，剩下的无非是在一般计算的中间插入一些计算步骤罢了。所以下面写的时候我先从最简单的求解器（简单方程、没有源项、常系数）开始写起，然后单独处理上述若干种不同情况。

ps 还有一种高级情况需要修改求解器，那就是不同的离散格式。这个我不太懂，感觉比较高端，以后再慢慢考虑吧。。
pps 这个笔记局限在不可压粘性流体里面。可压缩和欧拉方程我基本不懂。。。

0.1 顺序

所谓顺序，就是一个求解器从头到尾要干的事情。初学者对这个往往没概念导致看求解器代码的时候云山雾罩的。这里写清楚一点：

• 第一步是创建一个网格，实际上是根据生成的网格创建一个对象。一般不需要自己搞，只要一句话：

#include "createMesh.H"

• 接下来是往这个网格上放变量和系数，就是指定每个cv里的变量是多少。在求解器里也只是一句话：

#include "createFields.H"

但实际上没那么easy，因为具体要读取什么变量和系数每个求解器不太一样，有的时候所以要自己修改一些内容。

• 第三步就是离散方程组啦。具体怎么离散，取决于方程是怎么求解的。就我目前的了解，常见的方程无非是两种：压力速度耦合方程组，和被动标量方程。
  
  • 前者就是质量守恒+ns方程，一般可以写成：
    $$\nabla\cdot\hat u=0\\ \frac{\partial u}{\partial t} + \hat u\nabla\cdot \hat u=-\frac1\rho\nabla p+ \nabla\cdot(\nu\nabla\hat u))=-\frac1\rho\nabla p+ \nu\Delta\hat u$$
    这个是用来求解p和u的，虽然是非常基础的方程，但是处理起来其实有点难。。那些最难学的部分，什么simple族算法，piso，棋盘分布，rei-chow插值还有各类乱七八糟的对流项离散格式，基本都是用在这个上边了。因为这个比较成熟一般没啥改动，这个笔记就先掠过这部分，单独到simple族算法的时候再说。至于现在，我们只需要知道这个方程组的任务就是为下面的求解提供一个速度场（或者phi场，按照of的说法）就足够了。
  • 剩下的就是我比较喜欢的被动标量方程了。很多方程都是这个类别，比如能量方程（T），vof方法的相方程，level-set方法的f方程，rans族里面的各种湍动能什么的方程，这些不影响流动的，单纯只是被速度输运的量都可以归到被动标量方程里。这类方程的特点是：只依赖于速度，所以求解前需要知道速度；不作用于速度，所以不需要和ns方程联合求解。举个例子，在多孔介质的流体输运中有如下方程：
    $$\frac{\partial \theta}{\partial t}=\nabla\cdot(\kappa\nabla h)\\ \theta = g(h)=h$$
    虽然写起来是两个方程，但是第二个方程带入1式，真正求解的就是一个关于h的方程罢了（z是高度，已知）。这里面左边是瞬态项，右边是扩散项，连个速度都没有真是太好了。我们接下来就从这个方程出发试着写一个自己的求解器。
    简单起见，首先我们假定输运系数（渗透率）是常数。
    对于上述第二类方程，在of中直接调用各类离散函数就能完成。离散的结果是个特定的对象，一般是 fvVectorMatrix或者fvScalarMatrix，调用的函数就是fvm::ddt等。
• 接下来是求解。说到底只要一个solve函数就能搞定，例如：

pEqn.solve();

或者

fvVectorMatrix UEqn(...);
solve(UEqn == -fvc::grad(p));

，但是of毕竟是个通用计算框架，里面有大量的考虑松弛、非正交修正的语句。目前这个笔记只涉及最简单的矩形结构化网格，多余的因素尽量扣除。

• 最后就是输出了，调用语句：

runTime.write();

即可。

1. 创建一个空的求解器

当然可以从一个空文件从第一个字符开始敲，但是of提供了很多生成摸板的工具。19年初去学软件的时候就有所耳闻，但是如何创建一个空的求解器还是今天看了一个教程才知道的，那就是使用命令：

foamNewApp mySolver

这样就可以创建一个路径，里面有一个叫做mySolver.C的文件，文件里有写好的文件头和主函数的框子等等。这个源文件外头还有一个Make路径，Make里面也已经有了现成的files和options 文件了。
值得注意的是，这句命令最好在applications目录下运行，这样新生成的solver会和系统的求解器放在一个路径下，避免新手搞不定编译的问题。

2. 创建网格对象

在新建的空求解器中加入一句话：

#include "createMesh.H"

即可。这句话应该是将求解器路径下的/constant/polyMesh下的网格文件读入程序，编程一个对象。

3. 读取初场和系数：creatFields.H

接下来就是读取初场和系数。在上一句下面加一句：

#include "createFields.H"

当然这个要自己写，原因前面已经写过。所以我们在.C文件的同一个路径下新建一个creatFields.H文件。里面的内容，按照我们的示例方程也就是

$$\frac{\partial \theta}{\partial t}=\nabla\cdot(\kappa\nabla(h+z))\\ \theta = g(h)$$
来写。那么我们要读入的变量就是 h，要读入的系数就是 $\kappa$
变量和系数的读入（和读出，因为有些量，不只是变量是要输出来给人看的）在of里通过对象注册机来完成。新手完全没必要琢磨这是个啥，只要知道变量和系数是这样写入写出的就vans了。

读取系数

先读$\kappa$，仿照scalarTransportFoam的写法：

Info<< "Reading transportProperties\n" << endl;
IOdictionary transportProperties
(
    IOobject
    (
        "transportProperties",
        runTime.constant(),
        mesh,
        IOobject::MUST_READ_IF_MODIFIED,
        IOobject::NO_WRITE
    )
);

这一段原封不动抄进去就行。到这里是创建了一个IOdictionary的对象，叫做TransportProperties。创建的时候就用IoObject来初始化。如何初始化呢？首先从从/constant路径下的叫做"transportProperties"的文件里读入数据，这分别是IoObject里面前两行的意思。剩下的就是把这个东西注册（也就是写到）mesh上，mesh是第二步里创建出来的东西。
ok，现在我们有一个对象叫做transportProperties，这个对象是IOdictionary类的，也就拥有了这个类下的各种函数。那么我们要用的就是查找具体系数的功能。虽然这个例子里只有一个系数，看上去多次一举。但是真正的求解器会有若干个系数，所以要分别处理。写入$\kappa$的语句是：

Info<< "Reading diffusivity DT\n" << endl;
dimensionedScalar kappa
(
    transportProperties.lookup("kappa")
);

可以看到，我们创建了一个对象叫做kappa，是一个dimensionedScalar类的对象。这个对象初始化的时候值从哪来？使用lookup函数从文件里找来。这个函数从哪里来？从我们上一步创建的transportProperties里来。上一步我们说，transportProperties是个IOdictionary类的对象，lookup函数也是从这个类里来的。

读取初场

ok，接下来搞初场。我们的方程只有一个变量，也就只有一个初场。这个初场放在./0路径下面，就叫做h。读取的方式跟上面差不多。上面我们读取系数的时候创建了一个系数属于的对象，那么现在我们也要创建一个变量属于的对象。此处是h，是个标量所以我们要创建一个标量场volScalarField。在of里，写成：

Info<< "Reading field T\n" << endl;
volScalarField h
(
    IOobject
    (
        "h",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);

关于为什么这里有两个mesh，以及对象注册机的写详细说明，参阅资料1，2和3。从我个人的感觉来说，这部分属于软件架构和数据结构的内容，涉及的是如何把软件写的高效，并不涉及太多cfd的内容，如果仅仅是在of的框架下编写自己的求解器或许不需要太过在意。。。。

4.方程离散

接下来的任务是创建一个离散的方程组。方程组在of中一般以fvVectorMatrix或者fvScalarMatrix的类的对象存在。前者是用来储存速度这类矢量，压力、温度、vof相函数等标量储存在后者。此处我们创建一个fvScalarMatrix的对象hEqn：

 fvScalarMatrix hEqn
 (
    fvm::ddt(h)
    == fvm::laplacian(kappa, h)
);

可以看见这里的写法很类似与真实的数学写法。
常见的微分算符，其实就是瞬态微分，散度和laplace，分别写作ddt,div,laplacian。这三者又分别属于fvm和fvc两个命名空间，fvm是隐式，fvc是显式，取决与需要的离散格式。
当然还有一个要注意的地方是时间循环。在of里没有真正的稳态计算（？）所有的问题都是瞬态的，只是在稳态问题中可以放大时间步让方程迅速收敛到稳态而已（此处存疑）。所以在写稳态求解器的时候一样要写上瞬态项。
既然有瞬态项就必须有步进策略，而在of中步进一般由simple算法的框架来完成，其实这种简单求解器完全不需要simple算法，此处只是为了通用性而调用simple而已。
具体的语句是，在此前，就使用mesh生成一个叫做simple的对象，属于simpleControl类，语句是：

simpleControl simple(mesh);

接下来，在上述离散步骤的外面需要套一个循环：

while (simple.loop(runTime))
{
}

5.方程求解

构建好方程组之后，只需要一个简单的solve就能搞定求解：

hEqn.solve();

6. 结果输出

每个时间步都需要把结果输出来，所以加上一句：

runTime.write();

整个求解器就算初具雏形了。看一下整体的样子：

#include "fvCFD.H"// 必备头文件
#include "fvOptions.H"// 与源项有关，目前不涉及可以删掉
#include "simpleControl.H"// 时间推进有关
int main()// 主函数
{
    #include "setRootCaseLists.H"//必备
    #include "createTime.H"//创建时间对象
    #include "createMesh.H"//创建mesh 后续频繁使用
    simpleControl simple(mesh);//创建simple对象，用于控制时间推进
    #include "createFields.H"//创建场，读取初场和系数
    Info<< "\nCalculating h transport\n" << endl;//提示计算开始
    #include "CourantNo.H"// 计算库朗数
    while (simple.loop(runTime))// 时间循环开始
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;//输出当前时间
        while (simple.correctNonOrthogonal())//非正交修正，本求解器不涉及可以删掉
        {
            fvScalarMatrix hEqn//创建线性方程组
            (
                fvm::ddt(h)
              == fvm::laplacian(kappa, h)
            );
            hEqn.relax();//松弛，可以删掉
           // fvOptions.constrain(hEqn);//与源项有关，目前不涉及，可以删掉
            hEqn.solve();//求解
            //fvOptions.correct(h);//同上
        }
        runTime.write();//输出
    }
    Info<< "End\n" << endl;//提示计算结束
    return 0;
}

接下来要做的是把源代码编程一个可执行程序。所有的solver都是可执行程序。所以需要编译。编译通过wmake即可，之前的Make文件夹就是为这个准备的。

6.5 一个插曲

实际上示例的方程还是有点麻烦，因为$\theta$实际上是：

$$\theta = \frac{a}{b-h^c}$$
因为上述被动标量方程只能处理单个变量，所以应该写成：
$$\frac{\partial \theta}{\partial t}=\frac{\partial \theta}{\partial h}\frac{\partial h}{\partial t}=\kappa\Delta h$$
这样子有2种处理方法：

• 一种是直接把${\partial \theta}/{\partial h}$除过去，这样就变成了变系数问题
• 第二种就比较暴躁，直接扯出两个方程，一个是已知h求$\theta$，另一个是用求得的$\theta$反求h。这样写起来比较简单，所以就写了一个试试：

fvScalarMatrix thetaEqn
        (
            fvm::ddt(theta)
        );
        solve(thetaEqn == fvm::laplacian(kappa,h));//求解第一个方程得到theta
        h = -1*pow(-1600000+341532/(theta-0.075),0.252525);//代入求解h

没想到竟然支持这么狂野的表达式，竟然编译通过了。。。
等后面网格画好跑一下试试。

7. 编译

使用wmake，这个of自带的编译系统即可。因为之前我们使用foamNewApp命令生成了Make，内容都不用修改。直接在在当前求解器目录下：

wmake

就vans了。
这部分有个很强的教程可以参考。

8. 运行求解器

求解的前提是准备网格+初场。我不太喜欢blockMesh，又学不太会snappyHexMesh，所以网格都用商软来画，画完转换一下。这部分的内容见另一个笔记：使用ICEM生成的网格进行计算。

9. 变系数问题

10. 源项的处理

fvOptions