利用OpenFOAM开发一个求解波动方程的的求解器。

波动方程表达式为：

式中，c为波速，h为波高。

1 创建文件

利用foamNewApp快速创建文件结构。

run
foamNewApp demo13
cd demo13
touch createFields.H

文件结构如下所示。

2 源代码

• 头文件createFields.H中实现物理量的读取

// 需要从字典文件中读取变量c，并且定义待求标量h
 
Info << "读取transportProperties文件" << endl;
 
IOdictionary transportProperties
(
    IOobject
    (
        "transportProperties", //文件名
        runTime.constant(),    //文件位置
        mesh,
        IOobject::MUST_READ_IF_MODIFIED,
        IOobject::NO_WRITE
     )
);
 
// 定义标量C，该变量从字典文件中读取
dimensionedScalar C
(
    "C",         //字典中的关键字
    dimVelocity, // 速度量纲，也可以写成dimensionSet(0,1,-1,0,0)
    transportProperties
);
 
Info << "读取标量场h" << endl;
// 定义标量场h
volScalarField h
(
    IOobject
    (
        "h",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);

• 源文件demo13.C

#include "fvCFD.H"
int main(int argc, char *argv[])
{
    // 给个程序使用说明
    argList::addNote("本程序用于计算给定网格下的波动方程");
 
    #include "setRootCase.H"
    #include "createTime.H"
 
    // 添加fvMesh对象mesh
    #include "createMesh.H"
    // 包括前面定义的createFields.H头文件
    #include "createFields.H"
 
    Info << nl << "开始时间迭代计算" << endl;
    while (runTime.loop())
    {
        Info << nl << "Time = " << runTime.timeName() << endl;
        Info << "求解物理场h" << endl;
        // 控制方程
        fvScalarMatrix hEqn
        (
            // sqr函数为计算C的平方
            fvm::d2dt2(h) == fvm::laplacian(sqr(C), h)
        );
        hEqn.solve(); // 求解方程
 
        runTime.write(); // 写出数据
        Info << "Execution Time = " << runTime.elapsedClockTime() << "s" << endl;
    }
 
    // * * * * * * * * * * * * * * * * * * * * //
    Info << nl;
    runTime.printExecutionTime(Info);
 
    Info << "Endn" << endl;
    return 0;
}

3 测试案例

采用二维模型，几何尺寸为2x2m，网格尺寸0.02 m，计算域中心位置有一个0.2x0.2m的区域，将h的值设置为1，其他区域h=0。

这里挑重要的文件进行说明，细节可直接查看案例文件。

• system/fvSolution文件中需要添加h的求解方式

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSolution;
}
// * * * * * * * * * * * * * * //
solvers
{
 
    "h|hFinal"
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-08;
        relTol          0.01;
    }
 
}

• setFieldsDict字典文件进行初始化

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      setFieldsDict;
}
// * * * * * * * * * * * * //
 
defaultFieldValues
(
    volScalarFieldValue h 0
);
 
regions
(
    boxToCell
    {
        box (-0.1 -0.1 -1) (0.1 0.1 1);
 
        fieldValues
        (
            volScalarFieldValue h 1
        );
    }
 
);

• transportProperties文件中指定波速

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      transportProperties;
}
// * * * * * * * * * * * * * * * * * * //
// 定义波速C
C                   [0 1 -1 0 0 0 0]  0.05;

• 0/h文件指定待求量h的初始值与边界值

FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      h;
}
// * * * * * * * * * * * * //
dimensions      [0 1 0 0 0 0 0];
internalField   uniform 0;
boundaryField
{
    topBot
    {
        type            empty;
    }
 
    north
    {
        type fixedValue;
        value uniform 0;
    }
 
    south
    {
        type fixedValue;
        value uniform 0;
    }
 
    east
    {
        type fixedValue;
        value uniform 0;
    }
 
    west
    {
        type fixedValue;
        value uniform 0;
    }
}

计算结果如下图所示。

（完毕）