岳子写OpenFOAM代码20多年了（算上我导师的10多年）。发现一个问题，就是CFD公式需要经常查，代码也经常需要复制粘贴。
高级的程序员，往往需要最简朴素的复制粘贴。
那么从哪里复制？岳子的CFD代码，都是自己总结的。每次都来这里复制。现在把这些代码都公开给大家。
希望大家复制的高兴，复制的开心，一起冲最高级的CFD程序员！
随时间变化的速度进口
inlet
{
type uniformFixedValue;
uniformValue coded;
name pulse;
codeInclude
#{
#include "mathematicalConstants.H"
#};
code
#{
return vector
(
0,
0.5*(1 - cos(constant::mathematical::twoPi*x)),
0
);
#};
}
边界条件中调用fvc、mesh、字典文件、时间步长等相关信息
INLET
{
type codedFixedValue;
value uniform (10 0 0);
name linearTBC1;
codeInclude
#{
#include "fvCFD.H"
#};
codeOptions
#{
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude
#};
codeLibs
#{
-lmeshTools \
-lfiniteVolume
#};
code
#{
// 调用mesh
const fvMesh& mesh = this->patch().boundaryMesh().mesh();
// 调用字典文件
dictionary C = mesh.lookupObject<dictionary>("physicalProperties");
scalar test(readScalar(C.lookup("test")));
const fvPatch& inletPatch = this->patch();
// 调用时间步长
const scalar deltaT = mesh.time().deltaT().value();
const volVectorField& U = mesh.lookupObject<volVectorField>("U");
// 调用fvc
volTensorField gradU = fvc::grad(U);
// 计算inlet的面积
label patchID = mesh.boundaryMesh().findPatchID("inlet");
const polyPatch& myPatch = mesh.boundaryMesh()[patchID];
scalar patchArea = 0.0;
forAll(myPatch, faceI)
{
patchArea += mesh.magSf().boundaryField()[patchID][faceI];
}
// 访问U的值
volVectorField& U = mesh.lookupObjectRef<volVectorField>("U");
#};
}
OpenFOAM中基本CFD变量的声明
//需要在程序中植入
const scalarField& V = mesh.V(); // 网格体积
const surfaceVectorField& Sf = mesh.Sf(); // 网格面矢量
const surfaceVectorField& Cf = mesh.Cf(); // 网格面心位置
const surfaceScalarField& magSf = mesh.magSf(); // 网格面积的模
const surfaceVectorField& N = Sf/magSf; // 网格面法向量
const label& nCells = mesh.nCells(); // 网格单元数量
const label& nInternalFaces = mesh.nInternalFaces(); // 网格内部面数量
const meshCellZones& cellZones = mesh.cellZones(); // 网格cellZone编号
wordList zoneNames = mesh.cellZones().names(); // 网格cellZone的名字
label zoneI = mesh.cellZones().whichZone(celli); // 判断celli在哪个cellZone里
label i = cellZones[0][i]; // 网格cellZone[0]的真实网格编号
const fvPatchList& patches = mesh.boundary(); // 边界网格，是一系列patch的集合
const polyBoundaryMesh& boundaryMesh = mesh.boundaryMesh(); // 边界网格，是一系列patch的集合
const label patchID = boundaryMesh.findPatchID("inlet"); // inlet边界的ID
volScalarField S(magSqr(symm(tgradU()))); // 形变率的双点积
fvMesh mesh1 // 创建网格1
(
IOobject
(
"mesh1",
runTime.timeName(),
runTime,
IOobject::MUST_READ
)
);
fvMesh mesh2 // 创建网格2
(
IOobject
(
"mesh2",
runTime.timeName(),
runTime,
IOobject::MUST_READ
)
);
IOdictionary physicalProperties // 声明字典文件
(
IOobject
(
"physicalProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
);
bool training =
readBool(physicalProperties.lookup("training")); // 去字典文件中读取关键词
scalar RR(readScalar(generalProperties.lookup("RR"))); // 去字典文件中读取关键词
const volScalarField& T // 调取引用场
(
mesh().lookupObject<volScalarField>("T")
);
volScalarField T // 体心场
(
IOobject
(
"T",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
volScalarField B
(
IOobject
(
"B",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar(dimless, 0.0),
A.boundaryField().types() // 把A的边界条件赋值给B
);
volScalarField B
(
IOobject
(
"B",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar(dimless, 0.0),
zeroGradientFvPatchScalarField::typeName // B的边界条件均为零梯度
);
tmp<volScalarField> deltaLambdaT // tmp体心场
(
volScalarField::New
(
"deltaLambdaT",
mesh,
dimensionedScalar(dimless, 1.0)
)
);
surfaceScalarField phi //面心场
(
IOobject
(
"phi",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar(dimless, 0.0)
);
dimensionedScalar DT //带单位的值
(
"DT",
dimensionSet(0,0,0,0,0),
1.0
);
volScalarField::Internal TInternal("TInternal", T); //不带边界的体心场
PtrList<volScalarField> nutlist(BATCH); //场数组
forAll(nutlist, i)
{
nutlist.set
(
i,
new volScalarField
(
IOobject
(
"nutlist" + name(i),
mesh.time().timeName(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
),
mesh
)
);
}
volScalarField y(wallDist::New(mesh).y()); // 壁面距离，需要#include "wallDist.H"
U.primitiveFieldRef() = ... // 修改U的内部场的值
U.ref() = ... // 修改U的内部场的值
U.boundaryFieldRef() = ... // 修改U的边界场的值
bound(S, dimensionedScalar(S.dimensions(), 0)); // 对S场做有界处理
volScalarField& alphat = // 强制修改const常量
const_cast<volScalarField&>
(
mesh().lookupObject<volScalarField>("alphat")
);
List<label> markedCell; // 声明一个数组并赋值
for (int i = 0; i < 10; i++)
{
markedCell.append(i);
}
IOField<scalar> utau
(
IOobject
(
"utau",
runTime.constant(),
"../postProcessing", // 将场输出到postProcessing文件夹中
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
scalarField(totalFSize,0.0)
);
forAll(mesh.boundary(), patchi) // 判断壁面的边界类型
{ // 需要#include "wallFvPatch.H"
if (isA<wallFvPatch>(mesh.boundary()[patchi]))
{
}
if
(
isA<fixedValueFvPatchScalarField>
(
T.boundaryField()[patchi]
)
)
{
}
if
(
isA<processorPolyPatch>
(
mesh.boundaryMesn()[patchi]
)
)
{
}
}
禁用log输出
SolverPerformance<scalar>::debug = 0;
更改变量近壁网格体心值
label patchID = mesh.boundaryMesh().findPatchID("wall");
const scalarField TfaceCell
= T.boundaryField()[patchID].patchInternalField();
k.boundaryFieldRef()[patchID] == sqrt(TfaceCell);
矩阵操作
scalarSquareMatrix squareMatrix(3, Zero);
squareMatrix(0, 0) = 4;
squareMatrix(0, 1) = 12;
squareMatrix(0, 2) = -16;
squareMatrix(1, 0) = 12;
squareMatrix(1, 1) = 37;
squareMatrix(1, 2) = -43;
squareMatrix(2, 0) = -16;
squareMatrix(2, 1) = -43;
squareMatrix(2, 2) = 98;
const scalarField S(3, 1);
LUscalarMatrix L(squareMatrix);
scalarSquareMatrix inv(3); //矩阵求逆
L.inv(inv);
scalarField x(3, Zero);
scalarField Mx = squareMatrix*x; //矩阵乘以向量
scalarField x(L.solve(S)); //计算L \cdot x = S
监控代码计算时间
#include <chrono>
auto start = std::chrono::steady_clock::now();
// Functions here
auto end = std::chrono::steady_clock::now();
auto diff = end - start;
Info<< "Calculate nodes and weights, using "
<< std::chrono::duration <double, std::milli> (diff).count()
<< " ms" << endl;
显性、隐性离散
fvm::ddt(T) \\隐性时间项离散
fvm::laplacian(T) \\隐性拉普拉斯项离散
fvm::div(phi, T) \\隐性对流项离散
fvm::Sp(coeff, T) \\隐性源项离散
fvm::SuSp(coeff, T) \\依据coeff的符号进行隐性或显性源项离散
fvc::ddt(T) \\显性时间项离散
fvc::laplacian(T) \\显性拉普拉斯项离散
fvc::div(phi, T) \\显性对流项离散
fvc::grad(T) \\显性梯度项离散
在算例层面进行后处理操作：
// 输出多相界面、界面高度要对libwaves.so进行包含
libs
(
"libwaves.so"
);
// functions需要写入到算例的controlDict里
functions
{
// 输出等值面
#includeFunc isoSurface(isoField=alpha.water, isoValue=0.5, fields=())
// 输出波动速度场，注意，每个时间步文件夹下的UMean的数值要正确
uPrime
{
type subtract;
libs ("libfieldFunctionObjects.so");
fields (U UMean);
result uPrime;
executeControl writeTime;
writeControl writeTime;
}
// 输出波动速度场自身的张量积
uPrime2
{
type multiply;
libs ("libfieldFunctionObjects.so");
fields (uPrime uPrime);
result uPrime2;
executeControl writeTime;
writeControl writeTime;
}
// 输出空气龄
age
{
libs ("libfieldFunctionObjects.so");
type age;
diffusion on;
writeControl writeTime;
executeControl writeTime;
}
// 输出舒适度
comfort
{
libs ("libfieldFunctionObjects.so");
type comfort;
clothing 0.5;
metabolicRate 1.2;
extWork 0;
relHumidity 60;
writeControl writeTime;
executeControl writeTime;
}
// 输出热通量
wallHeatFlux1
{
type wallHeatFlux;
libs ("libfieldFunctionObjects.so");
region fluid;
patches (".*Wall");
}
// 输出形变率
dymmy
{
type coded;
libs ("libutilityFunctionObjects.so");
writeControl writeTime;
executeControl writeTime;
code
#{
#};
codeExecute
#{
const volVectorField& U
(
mesh().lookupObject<volVectorField>("U")
);
volScalarField strainRate(sqrt(2.0)*mag(symm(fvc::grad(U))));
strainRate.write();
#};
}
// 更改deltaT
setDeltaT
{
type coded;
libs ("libutilityFunctionObjects.so");
codeExecute
#{
const Time& runTime = mesh().time();
if (runTime.userTimeValue() >= -15.0)
{
const_cast<Time&>(runTime).setDeltaT
(
runTime.userTimeToTime(0.025)
);
}
#};
}
// 输出缓存场
cacheTemporaryObjects(kEpsilon:G);
#includeFunc writeObjects(kEpsilon:G)
// 输出失踪粒子
particles
{
libs ("liblagrangianFunctionObjects.so");
type particles;
}
// 输出场的最大值最小值
minMaxp
{
type fieldMinMax;
functionObjectLibs ("libfieldFunctionObjects.so");
fields
(
U
);
location no;
writeControl timeStep;
writeInterval 1;
}
// 输出力
forces
{
type forceCoeffs;
libs ("libforces.so");
writeControl timeStep;
writeInterval 1;
patches ("motorBike.*");
rho rhoInf; // Indicates incompressible
log true;
rhoInf 1; // Redundant for incompressible
liftDir (0 0 1);
dragDir (1 0 0);
CofR (0.72 0 0); // Axle midpoint on ground
pitchAxis (0 1 0);
magUInf 20;
lRef 1.42; // Wheelbase length
Aref 0.75; // Estimated
}
// 输出流线
#includeFunc streamlinesLine
(
funcName=streamlines, start=(0 0.5 0), end=(9 0.5 0), nPoints=24, U
)
// 输出某个场的体积分
volFieldValue1
{
type volFieldValue;
libs ("libfieldFunctionObjects.so");
writeControl writeTime;
log yes;
writeFields no;
regionType all;
name outlet;
operation volIntegrate;//min, max, etc.
//weightField phi;
fields
(
alpha
);
}
// 输出某个场的体平均
volFieldValue1
{
type volFieldValue;
libs ("libfieldFunctionObjects.so");
writeControl writeTime;
log yes;
writeFields no;
regionType all;
name outlet;
operation volAverage;//min, max, etc.
fields
(
T
);
}
// 将某个场的坐标系进行变更
cartesianToCylindrical
{
type cylindrical;
libs ("libfieldFunctionObjects.so");
origin (0 0 0);
axis (1 0 0);
field U;
writeControl outputTime;
writeInterval 1;
}
#includeFunc patchFlowRate(patch=outlet1) // 计算边界流率
#includeFunc faceZoneFlowRate(name=fz1) // 计算某个faceZone的流率
#includeFunc fieldAverage(U, p, prime2Mean = yes) // 输出场的时间平均
#includeFunc residuals // 输出残差
#includeFunc mag(UPrime2Mean) // 输出场的模
#includeFunc components(U) // 输出速度分量
#includeFunc MachNo // 输出马赫数
#includeFunc yPlus // 输出y+
#includeFunc Q // 输出Q
#includeFunc vorticity // 输出涡量
#includeFunc Lambda2 // 输出lambda
#includeFunc reconstruct(phi) // 输出重组的速度场
#includeFunc scalarTransport(T, alphal=1, alphat=1) // 输出s标量传输
#includeFunc patchAverage(patch=outlet, fields=(p U)) // 输出边界积分值
}
全部代码，请阅读原文
