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

波动方程表达式为：

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

1 创建文件

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

runfoamNewApp demo13cd demo13touch 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;// 定义标量场hvolScalarField 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 << "End\n" << 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;}// * * * * * * * * * * * * * * * * * * //// 定义波速CC                   [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;    }}

计算结果如下图所示。

图片

（完毕）

本站仅提供存储服务，所有内容均由用户发布，如发现有害或侵权内容，请点击举报。