🎉 Welcome to Aerosand!.
Docs
past
ofsp2025
ofsp2025-05

ofsp2025-05

参考：

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

前置： OpenFOAM开发编程基础04 网格类初步 | 𝓐𝓮𝓻𝓸𝓼𝓪𝓷𝓭 (aerosand.cn)

OpenFOAM 的计算绝大部分都是基于物理场的。基于前几次讨论，本文开始讨论 OpenFOAM 中场的基本操作。

为了操作方便，我们可以在 bashrc 中定义快捷命令（参考 OpenFOAM 环境准备 | 𝓐𝓮𝓻𝓸𝓼𝓪𝓷𝓭 (aerosand.cn)）

// bashrc
alias ofsp='cd /home/aerosand/aerosand/ofsp'

建立本文的项目文件夹并进入

// terminal
ofsp
mkdir 05_field
cd 05_field

场的创建

有单位量

dimensionSet

前文用到了 OpenFOAM 的单位预设，其定义在 $FOAM_SRC/OpenFOAM/dimensionSet/dimensionSets.C ，文件中还定义了很多预设单位

回忆 OpenFOAM 的单位设置为 [质量, 长度, 时间, 温度, 摩尔, 电流, 光强] 。

// $FOAM_SRC/OpenFOAM/dimensionSet/dimensionSets.C
...
const dimensionSet dimless;

const dimensionSet dimMass(1, 0, 0, 0, 0, 0, 0);
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0);
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0);
const dimensionSet dimTemperature(0, 0, 0, 1, 0, 0, 0);
const dimensionSet dimMoles(0, 0, 0, 0, 1, 0, 0);
const dimensionSet dimCurrent(0, 0, 0, 0, 0, 1, 0);
const dimensionSet dimLuminousIntensity(0, 0, 0, 0, 0, 0, 1);

const dimensionSet dimArea(sqr(dimLength));
const dimensionSet dimVolume(pow3(dimLength));
const dimensionSet dimVol(dimVolume);

const dimensionSet dimVelocity(dimLength/dimTime);
const dimensionSet dimAcceleration(dimVelocity/dimTime);

const dimensionSet dimDensity(dimMass/dimVolume);
const dimensionSet dimForce(dimMass*dimAcceleration);
const dimensionSet dimEnergy(dimForce*dimLength);
const dimensionSet dimPower(dimEnergy/dimTime);

const dimensionSet dimPressure(dimForce/dimArea);
const dimensionSet dimCompressibility(dimDensity/dimPressure);
const dimensionSet dimGasConstant(dimEnergy/dimMass/dimTemperature);
const dimensionSet dimSpecificHeatCapacity(dimGasConstant);
const dimensionSet dimViscosity(dimArea/dimTime);
const dimensionSet dimDynamicViscosity(dimDensity*dimViscosity);
...

我们可以设置自定义变量的物理单位，如

dimensionSet dimAerosand(1,2,3,4,5,6,7); // just for example

dimensionedType

在 OpenFOAM 代码中可以看到如下类似语句

dimensionedScalar airPressure
(
	"airPressure",
	dimensionSet(1, -1, -2, 0, 0, 0, 0),
	1.0
);

查找源代码如下

// terminal
find $FOAM_SRC -iname dimensionedType.H

找到并摘取其中的构造函数如下

...
	//- Construct from components (name, dimensions, value).
	dimensioned
	(
		const word& name,
		const dimensionSet& dims,
		const Type& val
	);
...

OpenFOAM 中还提供了其他的构造函数，我们可以在不同情况下构造自己的有单位物理量。目前阶段来说，我们不需要继续深挖源代码。点到为止，明白确实可以这样用，并且以后按照这种方式去用就可以了。

当然，dimensionedType.H 也定义了配套的成员函数，比如

	//- Return const reference to value.
	const Type& value() const noexcept { return value_; }

	//- Return non-const reference to value.
	Type& value() noexcept { return value_; }

由此，我们可以便捷的使用 value() 成员函数，例如对上面的 airPressure 取值，再计算最大值

maxP = max(airPressure.value());

注意，OpenFOAM 会检查数学计算式等号两遍的最终单位是否相同，如果不同则会强制报错并中断程序，后续代码会面临并解决这个问题。

volFields

我们会注意到 OpenFOAM 中有典型的语句如下

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

那么 volScalarField 是什么呢？又是如何构造场的呢？

我们讨论一下场是怎么被创建的。

// terminal 
find $FOAM_SRC -iname volFieldsFwd.H

也许有同学疑惑“我怎么知道要查看这个文件”。参考 OpenFOAM开发编程基础03 时间类初步 | 𝓐𝓮𝓻𝓸𝓼𝓪𝓷𝓭 (aerosand.cn) 第一节， 安装 OFextension 插件后，直接在 volScalarField 关键字上右键跳转即可。也可以打开 OpenFOAM Doxygen （OpenFOAM: User Guide: OpenFOAM®: Open source CFD : Documentation），查询相关关键字。除特殊情况，以后不再赘述。

volFieldFwd.H 中有类型定义如下

...
typedef GeometricField<scalar, fvPatchField, volMesh> volScalarField;
typedef GeometricField<vector, fvPatchField, volMesh> volVectorField;
...

可以看到，volFields 包括常见的 volScalarField 和 volVectorField，是体积场量。无论是 volScalarField 还是 volVectorField 其实都是 GeometricField 的模板的不同情况。所以， volFields 只是一种分类的称呼，而不是真正的类。

surfaceFields

同样的，面场量 surfaceFields 包括常见的 surfaceScalarField 和 surfaceVectorField 。

// terminal
find $FOAM_SRC -iname surfaceFieldsFwd.H

surfaceFieldFwd.H 中有类型定义如下

...
typedef GeometricField<scalar, fvsPatchField, surfaceMesh> surfaceScalarField;
typedef GeometricField<vector, fvsPatchField, surfaceMesh> surfaceVectorField;
...

可见，surfaceFields 也是 GeometricField 的模板的不同情况。

surfaceFields 的 surface 都可以写完整，为什么 volFields 的 volume 要简写成 vol 呢？也许是历史代码的原因吧。。。

GeometricField

进一步的，我们查看 GeometricField.H 文件

// terminal
find $FOAM_SRC -iname GeometricField.H

摘取类的主要结构和几个常用构造函数如下

...
template<class Type, template<class> class PatchField, class GeoMesh>
class GeometricField
:
    public DimensionedField<Type, GeoMesh>
{
public:
	...
private:
	...
public:
	...
	// Constructors

	//- Construct given IOobject, mesh, dimensioned<Type> and patch type.
	//  This assigns both dimensions and values.
	//  The name of the dimensioned\<Type\> has no influence.
	GeometricField
	(
		const IOobject& io,
		const Mesh& mesh,
		const dimensioned<Type>& dt,
		const word& patchFieldType = PatchField<Type>::calculatedType()
	);
	...
	//- Construct and read given IOobject
	GeometricField
	(
		const IOobject& io,
		const Mesh& mesh,
		const bool readOldTime = true
	);
	...
	//- Construct as copy resetting IO parameters
	GeometricField
	(
		const IOobject& io,
		const GeometricField<Type, PatchField, GeoMesh>& gf
	);
	...
...

上述构造函数，对应实际的场的构造，分别举例如下

...
	Info<< "Creating field light pressure\n" << endl;
	volScalarField lightP
	(
		IOobject
		(
			"lightP",
			runTime.timeName(),
			mesh,
			IOobject::NO_READ,
			IOobaject::AUTO_WRITE
		),
		mesh,
		dimensionedScalar("lightP",dimPressure,0.0)
		// 以构造的 dimensioned<Type>& 类型的变量初始化
		// 基于名称，单位，初始值
		// find $FOAM_SRC -iname dimensionedType.H
	);
...
	Info<< "Reading field p\n" << endl;
    volScalarField p
    (
        IOobject
        (
            "p",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );
	Info<< "Calculating field mixture density\n" << endl;
	volScalarField rho
	(
		IOobject
		(
			"rho",
			runTime.timeName(),
			mesh,
			IOobject::NO_READ,
			IOobject::AUTO_WRITE
		),
		alpha1*rho1 + alpha2*rho2
		// 传入计算场
	);

这三个例子对应了最常见的场的构造函数，暂时明白这三个即可。其他构造函数以及更深入的代码层面的解析在以后会陆续介绍，无需担心。读者也可以尝试多阅读几个。

场的时间推进

我们以压力场为例，结合场的创建，讨论一下场的时间推进计算。

应用准备 

// terminal 
foamNewApp 05_01_field
cd 05_01_field
cp -r $FOAM_TUTORIALS/incompressible/icoFoam/cavity/cavity debug_case
code .

脚本和说明

略。

单时间步计算

主源码 

#include "fvCFD.H"

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

int main(int argc, char *argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"

    #include "createMesh.H"
    // 常见三个头文件，以后不再赘述

    volScalarField p // 基于 IOobject 和 mesh 创建压力场p
    (
        IOobject // 基于以下参数创建 IOobject
        (
            "p",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE // 根据源代码中的设置自动写入
        ),
        mesh
    );
    // 读取可选：
	//     MUST_READ，必须读取，如果文件不存在则报错
	//     READ_IF_PRESENT，如果存在才读取，如果文件格式错误则报错
	//     NO_READ，不读取
	// 写入可选：
	//     AUTO_WRITE，根据要求自动写入
	//     NO_WRITE，不写入，但是可以通过代码显式写入

    runTime++; // 推进一个时间步
    // 如果省略此行，则不推进时间步，计算结果会写入到初始时间步，覆盖掉初始场

    forAll(p, i) // 遍历压力场
    {
        if (mesh.C()[i][1] < 0.5)
        // 还记得 mesh.C() 返回什么吗？
        // 返回网格中心点坐标的y坐标值（C++从0开始编号）
        {
            p[i] = i*runTime.value(); // 随便给一个可能的计算
        }
    }
    // 完成了整个压力场的计算

    Info<< "Max p = " << max(p).value() << endl; // 直接返回压力场最大值
    Info<< "Min p = " << min(p).value() << endl; // 返回压力场最小值
    Info<< "Average p = " << average(p).value() << endl; // 返回压力场平均值
    Info<< "Sum p = " << sum(p).value() << endl; // 返回压力场总和值

    p.write(); // 显式的代码，要求只写入该时间步的压力场值

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

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

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

    return 0;
}


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

编译运行

注意到主源码中的场的覆写是不区分初始时间步和后续时间步的，而且 foamCleanTutorials 命令不能恢复初始时间步文件夹内的场文件的内容。所以，为了避免初始条件被覆盖，强烈推荐备份初始条件。一旦初始场被覆盖，可以随时复原。操作如下

// terminal
cp -r debug_case/0 debug_case/0_orig

编译运行，结果如下

Create time

Create mesh for time = 0

Max p = 1.995
Min p = 0
Average p = 0.9975
Sum p = 399

ExecutionTime = 0 s  ClockTime = 0 s

End

查看 debug_case/ 文件夹下，多了 0.005/ 文件夹，其中的 0.005/p 文件中写入了计算的场。最后的输出就是最新时间步的场的值。

多时间步计算

如何计算多个时间步呢？

主源码

主源码修改如下

#include "fvCFD.H"

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

int main(int argc, char *argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"

    #include "createMesh.H"

    volScalarField p
    (
        IOobject
        (
            "p",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

    // runTime++;

    for (int i=0; i < 3; ++i) // 计算 3 个时间步
    {
        ++runTime; // 时间推进

        forAll(p, i)
        {
            if (mesh.C()[i][1] < 0.5)
            {
                p[i] = i*runTime.value();
            }
        }
        p.write(); // 写入该时间步
    }

    Info<< "Max p = " << max(p).value() << endl;
    Info<< "Min p = " << min(p).value() << endl;
    Info<< "Average p = " << average(p).value() << endl;
    Info<< "Sum p = " << sum(p).value() << endl;

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

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

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

    return 0;
}


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

关于写成 ++runTime 而不是 runTime++ 的讨论

编译运行

注意要先清理算例，避免错误，以后不再赘述。

// terminal
./caseclean
wclean

编译运行，结果如下

Create time

Create mesh for time = 0

Max p = 5.985
Min p = 0
Average p = 2.9925
Sum p = 1197

ExecutionTime = 0 s  ClockTime = 0 s

End

观察 debug_case/ 文件夹下多了三个时间步的文件夹，如下

|- debug_case/
	|- 0_orig/
		|- p
		|- U
	|- 0/
		|- p
		|- U
	|- 0.005/
		|- p
	|- 0.01/
		|- p
	|- 0.015/
		|- p
	|- constant/
	|- system/

查看各个时间步的压力场值，可以知道终端最后输出的值为最新时间步的计算值。

字典控制计算

即使是多时间步计算，我们也是直接修改主源码的计算时间步个数，每次修改都要重新编译整个应用。如果是大型应用，重新编译会相当麻烦。

在实际工作中，我们希望在应用编译完成后，也可以通过外部文件去控制计算，也就是使用 OpenFOAM 的字典文件 controlDict 控制计算。

主源码

主源码修改如下

#include "fvCFD.H"

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

int main(int argc, char *argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"

    #include "createMesh.H"

    volScalarField p
    (
        IOobject
        (
            "p",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

    while(runTime.loop()) // 时间推进
    {
        forAll(p, i)
        {
            if (mesh.C()[i][1] < 0.5)
            {
                p[i] = i*runTime.value();
            }
        }
		runTime.write(); // 按字典要求写入
    }

    Info<< "Max p = " << max(p).value() << endl;
    Info<< "Min p = " << min(p).value() << endl;
    Info<< "Average p = " << average(p).value() << endl;
    Info<< "Sum p = " << sum(p).value() << endl;

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

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

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

    return 0;
}


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

编译运行

清理调试算例，重新编译运行，结果如下

Create time

Create mesh for time = 0

Max p = 199.5
Min p = 0
Average p = 99.75
Sum p = 39900

ExecutionTime = 0 s  ClockTime = 0 s

End

调试算例 debug_case/ 文件多了几个时间步文件夹，分别是 0.1, 0.2, 0.3, 0.4, 0.5，写入的时间步正好由 system/controlDict 控制。

controlDict 字典的部分参数如下（每 0.005*20 = 0.1 写入一次）

...
deltaT          0.005;

writeControl    timeStep;

writeInterval   20;
...

这些参数都可以随时修改调试，无需重新编译应用。

场的基本操作

下面我们伪计算压力场、温度场、速度场。

  1. 求出网格单元和参考点的最大距离
  2. 根据时间、网格单元向量、参考向量、最大距离这 4 个参数计算压力场
  3. 根据压力梯度通过量纲处理，计算速度场

应用准备 

// terminal
foamNewApp 05_02_field
cd 05_02_field
cp -r $FOAM_TUTORIALS/incompressible/icoFoam/cavity/cavity debug_case

脚本和说明

略。

增加物理场

拷贝过来的调试算例并没有温度场，我们给算例添加初始温度场。

温度和压力一样是标量，所以先拷贝压力场文件再修改。

// terminal
cp -r debug_case/0/p debug_case/0/T
code debug_case/0/T

温度场初始条件 debug_case/0/T，内容如下

FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      T;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

// dimensions      [0 2 -2 0 0 0 0];
dimensions              [0 0 0 1 0 0 0];

// internalField   uniform 0;
internalField           uniform 2023;

boundaryField
{
    movingWall
    {
        // type            zeroGradient;
        type            fixedValue;
        value           uniform 2049;
    }

    fixedWalls
    {
        type            zeroGradient;
    }

    frontAndBack
    {
        type            empty;
    }
}


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

主源码 

#include "fvCFD.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
scalar calculatePressure(scalar t, vector x, vector x0, scalar scale);
// 自定义函数原型声明


int main(int argc, char *argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"

    #include "createMesh.H"

    Info<< "Reading transportProperties\n" << endl;
    IOdictionary transportProperties // 构造字典对象，从字典文件中读入参数
    (
        IOobject
        (
            "transportProperties", // 字典文件名称
            runTime.constant(), // 字典文件路径
            mesh, // 和 mesh 建立关系
            IOobject::MUST_READ_IF_MODIFIED,
            IOobject::NO_WRITE
        )
    );
    dimensionedScalar nu("nu",dimViscosity,transportProperties); // 读取字典文件中的参数

    Info<< "Reading field p\n" << endl;
    volScalarField p // 构造 p 场
    (
        IOobject // 基于 IOobject
        (
            "p", // 场文件的名称
            runTime.timeName(), // 场文件的路径
            mesh, // IOobject 和 mesh 建立关系
            IOobject::MUST_READ, // 必须要给初始条件
            IOobject::AUTO_WRITE
        ),
        mesh // 和 mesh 建立关系
    );

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

    Info<< "Reading field U\n" << endl;
    volVectorField U // 构造 U 场
    (
        IOobject
        (
            "U",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

	// 举例创建空值场
	// 参考上文讨论的 GeometircField 构造函数
    volScalarField zeroScalarField
    (
        IOobject
        (
            "zeroScalarField",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ, // 无需给初始条件的文件
            IOobject::AUTO_WRITE
        ),
        mesh,
        dimensionedScalar("zeroScalarField",dimless,Zero) // 初始化为零
    );

    volVectorField zeroVectorField
    (
        IOobject
        (
            "zeroVectorField",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        mesh,
        dimensionedVector("zeroVectorField",dimless,Zero) // 初始化为零
    );

    volVectorField gradT // 构造 gradT 场
    (
        IOobject
        (
            "gradT",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        fvc::grad(T) // 因为是从 T 场计算而来，所以不需要再次和 mesh 建立关系
    );

    Info<< "Reading/calculating face flux field phi\n" << endl;
    surfaceScalarField phi // 构造 phi 场
    (
        IOobject
        (
            "phi",
            runTime.timeName(),
            mesh,
            IOobject::READ_IF_PRESENT,
            IOobject::AUTO_WRITE
        ),
        fvc::flux(U)
        // fvc::interpolate(U)*mesh.Sf(); // 等同于此行计算
    );
    

    const vector originVector(0.05, 0.05, 0.005); // 构造带初值vector类的对象originVector

    const scalar rFarCell  = // 构造scalar类的对象rFarCell
    max // 求最大值
    (
        mag // 求数值大小
        (
            dimensionedVector("x0", dimLength, originVector) 
            // 以originVector创建dimensionedVector
            - mesh.C() // 减去网格单元中心点坐标，实际上也是向量
        )
    ).value(); // 返回值的大小

    Info<< "\nStarting time loop\n" << endl;

    while (runTime.loop()) // 时间推进
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

        for (label cellI=0; cellI<mesh.C().size(); ++cellI) // 遍历所有网格单元
        {
            p[cellI] = calculatePressure // 通过自定义函数求 p 场
            (
                runTime.time().value(),
                mesh.C()[cellI],
                originVector,
                rFarCell
            );
        }

        U = fvc::grad(p) * dimensionedScalar("tmp",dimTime,1.0); // 求 U 场
        // 额外乘以一个含单位临时量 tmp ，保证两边结果单位相同

        runTime.write(); // 写出计算的场
    }

    Info<< "Calculation done." << endl;

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

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

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

    return 0;
}


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

// 自定义函数的实现
scalar calculatePressure(scalar t, vector x, vector x0, scalar scale)
{
    scalar r(mag(x - x0) / scale); // 构造scalar类的对象r，并初始化

    scalar rR(1.0 / (r + 1e-12)); // 同上

    scalar f(1.0); // 同上

    return Foam::sin(2.0 * Foam::constant::mathematical::pi * f * t) * rR;
    // 返回一个计算值，并没有什么物理意义
}

编译运行，结果如下

Create time

Create mesh for time = 0

Reading transportProperties

Reading field p

Reading field T

Reading field U

Reading/calculating face flux field phi


Starting time loop

Time = 0.005

Time = 0.01

Time = 0.015

Time = 0.02

...


ExecutionTime = 0.01 s  ClockTime = 0 s

End

使用 paraview 将计算结果可视化。注意打开需要路径

// terminal
paraFoam -case debug_case

注意到 debug_case/ 的时间步文件夹下新建了很多场文件，以 debug_case/0.3/ 为例，其文件结构如下

| - 0.3/
	|- uniform/
	|- gradT
	|- p
	|- phi
	|- T
	|- U
	|- zeroScalarField
	|- zeroVectorField

可以查看各个场文件中的数值。

注意：

打开 0.3/phi 可以看到并没有计算更新值，其他非 p,U 文件也没有计算更新值。这是因为场在创建时候给定的计算公式只是初始化，并不会随着时间推进而计算更新。所以，需要在主源码的时间循环中去计算需要的场。

代码整理

在 OpenFOAM 实际应用中，一般的

  • “场的接入”放入单独的文件 createFields.H
  • OpenFOAM 本身提供 createPhi.H
  • 自定义方法也放入单独文件，如 calculatePressure.H

所以该应用的文件结构调整后为

|- userApp/
	|- debug_case/
	|- Make/
		|- files
		|- options
	|- createFields.H
	|- calculatePressure.H
	|- 05_02_field.C

场接入文件 createFields.H 为

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

    dimensionedScalar nu("nu",dimViscosity,transportProperties);

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

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

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

    volScalarField zeroScalarField
    (
        IOobject
        (
            "zeroScalarField",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        mesh,
        dimensionedScalar("zeroScalarField",dimless,Zero)
    );

    volVectorField zeroVectorField
    (
        IOobject
        (
            "zeroVectorField",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        mesh,
        dimensionedVector("zeroVectorField",dimless,Zero)
    );

    volVectorField gradT
    (
        IOobject
        (
            "gradT",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        fvc::grad(T)
    );

    Info<< "Reading/calculating face flux field phi\n" << endl;
    surfaceScalarField phi
    (
        IOobject
        (
            "phi",
            runTime.timeName(),
            mesh,
            IOobject::READ_IF_PRESENT,
            IOobject::AUTO_WRITE
        ),
        fvc::flux(U)
        // fvc::interpolate(U)*mesh.Sf();
    );

自定义方法 calculatePressure.H 为

scalar calculatePressure(scalar t, vector x, vector x0, scalar scale)
{
    scalar r(mag(x - x0) / scale);

    scalar rR(1.0 / (r + 1e-12));

    scalar f(1.0);

    return Foam::sin(2.0 * Foam::constant::mathematical::pi * f * t) * rR;
}

最终主源码整理为

#include "fvCFD.H"

#include  "calculatePressure.H"

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

int main(int argc, char *argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"

    #include "createMesh.H"

    #include "createFields.H"
    
    const vector originVector(0.05, 0.05, 0.005);

    const scalar rFarCell  = 
    max
    (
        mag
        (
            dimensionedVector("x0", dimLength, originVector)
            - mesh.C()
        )
    ).value();

    Info<< "\nStarting time loop\n" << endl;

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

        for (label cellI=0; cellI<mesh.C().size(); ++cellI)
        {
            p[cellI] = calculatePressure
            (
                runTime.time().value(),
                mesh.C()[cellI],
                originVector,
                rFarCell
            );
        }

        U = fvc::grad(p) * dimensionedScalar("tmp",dimTime,1.0);

		phi = fvc::flux(U); // 增加 phi 的计算

        runTime.write();
    }

    Info<< "Calculation done." << endl;

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

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

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

    return 0;
}


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

整理后的这种形式更接近 OpenFOAM 求解器等应用的代码组织形式。

  • 主源码功能划分非常清晰
  • 每个功能都更加易于维护

编译运行结果当然还是一样的。

小结

回顾这 5 篇讨论，我们大概可以感受到数值计算的核心要素所在——时间、网格、物理场。有了这些对象之后，我们便可以组建线性代数方程，进行数值求解。不要着急，我们距离求解器编程越来越近了。

Last updated on
ofsp2025-04ofsp2025-06