ofsp2025-08
参考:
- https://github.com/UnnamedMoose/BasicOpenFOAMProgrammingTutorials
- https://www.topcfd.cn/simulation/solve/openfoam/openfoam-program/
- https://www.tfd.chalmers.se/~hani/kurser/OS_CFD/
- https://github.com/ParticulateFlow/OSCCAR-doc/blob/master/openFoamUserManual_PFM.pdf
- https://www.youtube.com/watch?v=KB9HhggUi_E&ab_channel=UCLOpenFOAMWorkshop
- http://dyfluid.com/#
- https://ss1.xrea.com/penguinitis.g1.xrea.com/study/OpenFOAM/index.html
感谢原作者们的无私引路和宝贵工作。
计算流体力学的通用基本方程为
求解的波动方程如下
我们依然用 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
无需修改,我们仍然使用求解控制的设置(步长步数等)decomposedParDict
和PDRblockMeshDict
不用管,本文用不到
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
的时候要注意坐标点 - 设置
boxToCell
中box
的两个角点坐标 - 注意,
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
。
虽然系列讨论追求给初学者带来“连续”“逐深”“有方向”“不间断”的学习体验,但这是避不开的一环。