🎉 Welcome to Aerosand!.
ofsp2025-08

ofsp2025-08

参考:

感谢原作者们的无私引路和宝贵工作。

前置: OpenFOAM开发编程基础07 第一个求解器 | 𝓐𝓮𝓻𝓸𝓼𝓪𝓷𝓭 (aerosand.cn)

计算流体力学的通用基本方程为

t(ρϕ)+(ρUϕ)=(Γϕ)+Sϕ\frac{\partial}{\partial t}(\rho \phi) + \nabla \cdot (\rho U\phi) = \nabla\cdot(\Gamma\nabla\phi) + S_{\phi}

求解的波动方程如下

2At2=(c2A)\frac{\partial^2 A}{\partial t^2} = \nabla \cdot (c^2 \nabla A)

我们依然用 A 表示待求的物理场。

新建本文的项目文件夹

// terminal
ofsp
mkdir 08_waveSol
cd 08_waveSol

应用准备

// terminal
foamNewApp 08_01_waveSol
cd 08_01_waveSol
cp /home/aerosand{caserun,caseclean} .
cp -r $FOAM_TUTORIALS/incompressible/icoFoam/cavity/cavity .
cp -r debug_case/0 debug_case/0_orig
code .

脚本和说明

略。

求解器

createFields.H

场文件 /userApp/createFields.H 接入 方程所需要的场 A 和参数 c

以后在开发大型求解器的时候,createFields.H 文件会有相当多行代码,创建的场的功能各不相同,读取的字典参数的功能也是各不相同。建议可以在文件开头写一个简单的“目录”。这样子可以通过搜索直达内容分区。

/*
* Contents
*   # Basic fields
*   # Transport properties
*/


/*
* # Basic fields
*/

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


/*
* # Transport properties
*/

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

Info<< "Reading c from transportProperties" << endl;
dimensionedScalar c("c",dimViscosity,transportProperties);

主源码

文件 /userApp/08_01_waveSol.C 其实很简单,就是构建一个时间推进,迭代内,每个时间步求解方程即可。

#include "fvCFD.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
    argList::addNote(
        "This application solves the wave equation over given mesh."
    );

    #include "setRootCase.H"
    #include "createTime.H"

    #include "createMesh.H"

    #include "createFields.H"

    Info<< nl << "Starting time loop..." << endl;

    while (runTime.loop())
    {
        Info<< nl << "Time = " << runTime.timeName() << endl;

        Info<< "Solving h field" << endl;

        fvScalarMatrix hEqn
        (
            fvm::d2dt2(h)
            == fvm::laplacian(sqr(C), h)
        );

        hEqn.solve();

        runTime.write();

        Info<< nl;
        runTime.printExecutionTime(Info);

    }

    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //


    Info<< "End\n" << endl;

    return 0;
}

编译

// termianl
wclean
wmake

调试算例

  • system/blockMeshDict 无需修改,我们仍然使用原几何模型和几何网格
  • system/controlDict 无需修改,我们仍然使用求解控制的设置(步长步数等)
  • decomposedParDictPDRblockMeshDict 不用管,本文用不到

transportProperties

文件 debug_case/constant/transportProperties 内容修改如下

C              0.05;

初始条件

基于初始条件 p 修改得到 A ,内容如下

FoamFile
{
    version     2.0;
    format      ascii;
    arch        "LSB;label=32;scalar=64";
    class       volScalarField;
    location    "0";
    object      A;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 1 0 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    movingWall
    {
        type            fixedValue;
        value           uniform 0;
    }
    fixedWalls
    {
        type            fixedValue;
        value           uniform 0;
    }
    frontAndBack
    {
        type            empty;
    }
}


// ************************************************************************* //

setFieldsDict

使用命令创建一个模板到调试算例下

// terminal
foamGetDict -case debug_case setFieldsDict

修改文件 /debug_case/system/setFieldsDict ,内容如下

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      setFieldsDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

defaultFieldValues
(
    volScalarFieldValue A 0
);

regions
(
    // Set cell values
    // (does zerogradient on boundaries)
    boxToCell
    {
        box (0.04 0.04 0) (0.06 0.06 0.01);

        fieldValues
        (
            volScalarFieldValue A 1
        );
    }
);

// ************************************************************************* //
  • 注意调试算例的 blockMeshDict 中设置了 scalar 0.1; ,几何尺寸有了缩放,所以在设置 box 的时候要注意坐标点
  • 设置 boxToCellbox 的两个角点坐标
  • 注意,setFields 只重设了初始场,区别于场里的

fvSchemes

增加时间项二阶偏导的数值计算离散格式,其他不变

...
d2dt2Schemes
{
    default         Euler;
}

fvSolution

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{
    A
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-06;
        relTol          0.05;
    }

    AFinal
    {
        $p;
        relTol          0;
    }
}

// ************************************************************************* //

也可以写成 A|AFinal ,然后一起设置

controlDict

随意设置计算时间长短等。

脚本

因为额外使用了 setFields 命令,且该命令会更改初始条件,所以需要修改脚本。

caserun 内容如下

#!/bin/sh
cd "${0%/*}" || exit 1                              # Run from this directory
#------------------------------------------------------------------------------

appPath=$(cd `dirname $0`; pwd)                     # Get application path
appName="${appPath##*/}"                            # Get application name

# cd debug_case

# echo $appName

blockMesh -case debug_case | tee debug_case/log.mesh
echo "Meshing done."

cp -r debug_case/0_orig/ debug_case/0/
setFields -case debug_case
echo "Setting fileds done."

$appName -case debug_case | tee debug_case/log.run

caseclean 内容如下

#!/bin/sh
cd "${0%/*}" || exit 1                              # Run from this directory
#------------------------------------------------------------------------------

appPath=$(cd `dirname $0`; pwd)
appName="${casePath##*/}"

cd debug_case

rm -rf log.*
foamCleanTutorials
rm -rf 0/
echo "Cleaning done."

计算

// terminal
./caseclean
./caserun

paraFoam -case debug_case

使用 paraview 可以看到中心点场随着时间的波动变化。

小结

本文讨论求解了一个非定常波动问题。我们通过注释和代码拆出,尝试把求解器结构写的更加清晰易读。首次使用了 OpenFOAM 自带的工具 setFields (本质上也是一种 应用)。

需要注意的是,因为待求的波动方程简单,所以本文求解器对偏微分方程的求解并没有使用任何的算法,只是简单的在每个时间步求解。实际开发中,更多的问题是多个物理场耦合,所以我们需要一些数值算法来求解偏微分方程组,而不是单个偏微分方程。

虽然本系列是开发编程系列,但是在正式进入 OpenFOAM 标准求解器开发前,我们在下一次讨论中不得不讲一下 OpenFOAM 的基础求解算法,主要是 SIMPLE, PISO, PIMPLE

虽然系列讨论追求给初学者带来“连续”“逐深”“有方向”“不间断”的学习体验,但这是避不开的一环。

Last updated on