OpenFOAM小代码-常用函数
如得到owner和neighbour网格中心的距离？
mesh.surfaceInterpolation::deltaCoeffs() 返回的体心距离的倒数
某个patch的第一个face在boundary中的所有face中的索引
label bFacei = bm[patchi].patch().start() - mesh.nInternalFaces();
//-------------------来自CFD中文网,持续更新中..........---------------------------------
by @李东岳
计算Patch的面积
// 计算某个patch的面积 http://www.cfd-china.com/topic/3498
{
label patchID = mesh.boundaryMesh().findPatchID("inlet");
// Create a polyPatch for looping
const polyPatch& myPatch = mesh.boundaryMesh()[patchID];
// Initialize patchArea
scalar patchArea = 0.0;
// Loop trhough all faces on the polyPatch, adding their magnitude surface
// area vectors
forAll(myPatch, faceI)
{
patchArea += mesh.magSf().boundaryField()[patchID][faceI];
}
}
更改壁面一层网格值
label patchID = mesh.boundaryMesh().findPatchID("wall");
const scalarField& TfaceCell =
T.boundaryField()[patchID].patchInternalField();
k.boundaryFieldRef()[patchID] = sqrt(TfaceCell);
给定patch调用网格
如果p是一个fvPatch
const fvMesh& mesh = p.boundaryMesh().mesh();
类内调用时间步长
假定mesh_是类内可以调用的成员变量。
const scalar deltaT = mesh_.time().deltaT().value();
如果mesh_不可获取，U_可获取：
const scalar deltaT = U_.mesh().time().deltaT().value();
声明一个tmp变量
tmp<volScalarField> deltaLambdaT
(
volScalarField::New
(
"deltaLambdaT",
mesh_,
dimensionedScalar(dimless, 1.0)
)
);
简写边界场
volScalarField::Boundary& meanUb = meanU_.boundaryFieldRef();
forAll(meanUb, patchi)
{
fvPatchScalarField test = meanUb[patchi];
}
判断是否是固定值边界、processor边界
forAll(mesh_.boundary(), P)
{
if
(
isA<fixedValueFvPatchScalarField>
(
T.boundaryField()[P]
)
||
isA<processorPolyPatch>(mesh_.boundaryMesh()[P])
)
{
}
}
codedFixedValue
wall
{
type codedFixedValue;
value uniform 300; //default value
redirectType linearTBC; //name of new BC type
code
#{
const vectorField& Cf = patch().Cf(); // get face center coordinate;
//assume the wall is a vertical wall
scalar ymax = max(Cf&vector(0,1,0)); // `&` is dot product
scalar ymin = min(Cf&vector(0,1,0));
// Info<<"ymax="<<ymax<<",ymin="<<ymin<<nl;
vectorField& vf = *this; // get the boundary field (List of vectors defined at face centers) to fill.
//temporal coordinate, type: scalar
scalar t =this->db().time().value(); // get time
// temporal term
// scalar tt = 1+0.1*sin(5*t);
scalar tt = 1.0;
forAll(Cf,faceI)
{
//spatial coordinates, type: scalar
//scalar x = Cf[faceI].x();
scalar y = Cf[faceI].y();
//scalar z = Cf[faceI].z();
vf[faceI] = ((y-ymin)/(ymax-ymin)*500+300)*tt;
}
#};
//I do not know why I need to add those things
//codeInclude
//#{
// #include "fvCFD.H"
//#};
//codeOptions
//#{
// -I$(LIB_SRC)/finiteVolume/lnInclude
//#};
}
努塞尔数
下列代码可以添加到controlDict中执行
functions
{
coded
{
functionObjectLibs ( "libutilityFunctionObjects.so" );
enabled true;
type coded;
redirectType printMinU;
writeControl timeStep;
writeInterval 1;
codeOptions
#{
-I$(LIB_SRC)/meshTools/lnInclude
#};
codeExecute
#{
const volScalarField& T
(
mesh().lookupObject<volScalarField>("T")
);
const fvPatchList& patches = mesh().boundary();
forAll(patches, patchi)
{
const fvPatch& currPatch = patches[patchi];
if (mesh().time().value() == 30001)
{
if (currPatch.name() == "downwall")
{
fvPatchScalarField nus = T.boundaryField()[patchi];
scalar L = 0.014231;
scalar Tref = 309.1;
scalar Timp = 347.6;
nus = T.boundaryField()[patchi].snGrad()*L/(Timp - Tref);
forAll(T.boundaryField()[patchi], facei)
{
Pout << mesh().C().boundaryField()[patchi][facei].x()/0.014231
<< " " << nus[facei] << nl;
}
}
}
}
#};
}
}
获取某个函数的计算时间
#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;
在某个文件夹下输出scalarField
by @Samuel-Tu
IOField<scalar> utau
(
IOobject
(
"utau",
runTime.constant(),
"../postProcessing",
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
scalarField(totalFSize,0.0)
);
BlockMesh-simpleGrading
处理非均一网格
blocks
(
hex (0 1 2 3 4 5 6 7) (100 300 100)
simpleGrading
(
1 // x-direction expansion ratio
(
(0.2 0.3 4) // 20% y-dir, 30% cells, expansion = 4
(0.6 0.4 1) // 60% y-dir, 40% cells, expansion = 1
(0.2 0.3 0.25) // 20% y-dir, 30% cells, expansion = 0.25 (1/4)
)
3 // z-direction expansion ratio
)
);
输出yplus
functions
{
yplus
{
type yPlus;
functionObjectLibs ("libutilityFunctionObjects.so")；
outputControl outputTime；
enabled true;
}
}
by @Sloan
在流场中添加颗粒失踪
controlDict添加下述内容，同时需要在constant文件夹添加kinematicProperties字典文件
functions
{
tracks
{
libs ("liblagrangianFunctionObjects.so");
type icoUncoupledKinematicCloud;
U U.water;//速度场名字
}
}
在流场中添加欧拉标量传输
controlDict添加下述内容，同时需要在0文件夹下添加s标量场，以及fvScheme、fvSolution都要做适配
functions
{
s
{
type scalarTransport;
libs ("libsolverFunctionObjects.so");
resetOnStartUp no;
field s;
phi phi.water;
}
}
输出场的最大值最小值
minMax
{
type fieldMinMax;
functionObjectLibs ("libfieldFunctionObjects.so");
enabled true;
log true;
write true;
location false;
fields
(
m00
U.water
);
}
对网格进行赋值
forAll(T, C)//遍历网格，T为场
{
T[C] = min(nu[C], 100);//对T的第C个网格，赋值第C个网格的nu与100的最小值
}
声明多个场PtrList
PtrList<surfaceScalarField> abPosCoeff(4);
forAll(abPosCoeff, i)
{
abPosCoeff.set
(
i,
new surfaceScalarField
(
IOobject
(
"abPosCoeff" + name(i),
mesh.time().timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionedScalar(inv(dimLength), 0.0)
)
);
}
补充一个CodeStream实现充分发展的湍流入口，请东岳老师不要建议
dw为到壁面的最短距离
OpenFoam 7：
gasinlet
{
type fixedValue;
value #codeStream
{
codeInclude
#{
#include "fvCFD.H"
#include "DynamicList.H"
#};
codeOptions
#{
-I(LIB_SRC)/meshTools/lnInclude
#};
//libs needed to visualize BC in paraview
codeLibs
#{
-lmeshTools \
-lfiniteVolume
#};
code
#{
const IOdictionary& d = static_cast<const IOdictionary&>
(
dict.parent().parent()
);
const fvMesh& mesh = refCast<const fvMesh>(d.db());
const label id = mesh.boundary().findPatchID("gasinlet");
const fvPatch& patch = mesh.boundary()[id];
vectorField U(patch.size(), vector(0, 0, 0));
const scalar U_inf=93.5;
const scalar LT=0.2;
const scalar visco=2.3183558929634845e-06;
const scalar deltaT=0.37*pow(visco/U_inf,0.2)*pow(LT,0.8);
scalar dw=0.0;
forAll(U, i)
{
const scalar x = patch.Cf()[i][0];
const scalar y = patch.Cf()[i][1];
const scalar z = patch.Cf()[i][2];
DynamicList<scalar> distance(4);
distance.append(mag(y));
distance.append(mag(y-30e-6));
distance.append(mag(z-23e-6));
distance.append(mag(z+23e-6));
dw=min(distance);
distance.clearStorage();
scalar Ux=0.0;
if(deltaT>dw)
{
Ux=U_inf*pow(dw/deltaT,1.0/7.0);
}
else
{
Ux=U_inf;
}
U[i] = vector(Ux, 0, 0);
}
writeEntry(os, "", U);
#};
};
}
@马乔
wordHashSet unMatchedkeys(dict.toc());
forAll(names, nameI)
{
const word& name = names[nameI];
const entry* ePtr = dict.findEntry(name, keyType::REGEX);
if(ePtr)
unMatchedKeys.erase(ePtr->keyWord());
}
class derived
:public base
{
public:
derived() {}
void printInfo()
{
Info << "derived class" << endl;
}
};
int main(int argc, char *argv[])
{
...
derived d;
base::store(&d);
...
}
fvOption小工具
根据GitHub上一个python脚本改写了一个生成fvOption的小工具 https://github.com/Veenxz/fvSchemes
用Foam没有多久，各位大佬可以多提提意见。
# version 2.0
# Based on https://github.com/Carlopasquinucci/fvSchemes
# Generation by Veenxz
# relaesed under license GPL GNU 3.0
steady = True
pseudo_transient = False
precision = 2 # First order 1 or Second order 2
unbounded = False
LUST = False
secondorder = True
maxOrtho = 80
maxSkew = 20
# Header and Footer
h = [
"/*--------------------------------*- C++ -*----------------------------------*|",
"| ========= | |",
"| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |",
"| \\ / O peration | Website: https://openfoam.org |",
"| \\ / A nd | Version: 7 |",
"| \\/ M anipulation | |",
"\*---------------------------------------------------------------------------*/",
"FoamFile", "{", " version 2.0;", " format ascii;",
" class dictionary;", ' location "system";', " object fvSchemes;",
"}",
"// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //",
""
]
footer = '\n// ************************************************************************* //'
nonOrthogonalCorrectors = 1
if (maxOrtho) > 80:
print(
'Warning: mesh is not so nice. Use 3 nonOrthogonalCorrectors in fvSolution file'
)
nonOrthogonalCorrectors = 3
if (maxSkew) > 8:
print('Warning: mesh is not so nice')
if (maxOrtho) > 80:
gradSchemes = (
'{\n default cellLimited Gauss linear 0.5;\n'
' grad(U) faceLimited Gauss linear 1.0;\n}'
)
divSchemes = (
'{\n div(phi,U) Gauss linearUpwind grad(U);\n'
' div(phi,omega) Gauss upwind;\n'
' div(phi,k) Gauss upwind;\n'
' div(phi,e) Gauss upwind;\n'
' div((nuEff*dev(T(grad(U))))) Gauss linear;\n}'
)
laplacianSchemes = (
'{\n default Gauss linear limited 0.333;\n}'
)
snGradSchemes = (
'{\n default Gauss linear limited 0.333;\n}'
)
blending = 0.2
nonOrthogonalCorrectors = 3
if (maxOrtho) > 70:
gradSchemes = (
'{\n default cellLimited Gauss linear 0.5;\n'
' grad(U) cellLimited Gauss linear 1.0;\n}'
)
divSchemes = (
'{\n div(phi,U) Gauss linearUpwind grad(U);\n'
' div(phi,omega) Gauss linearUpwind grad(omega);\n'
' div(phi,k) Gauss linearUpwind grad(k);\n'
' div(phi,e) Gauss linearUpwind grad(e);\n'
' div((nuEff*dev(T(grad(U))))) Gauss linear;\n}'
)
laplacianSchemes = (
'{\n default Gauss linear limited 0.5;\n}'
)
snGradSchemes = (
'{\n default Gauss linear limited 0.5;\n}'
)
blending = 0.5
nonOrthogonalCorrectors = 3
if (maxOrtho) > 60:
gradSchemes = (
'{\n default cellMDLimited Gauss linear 0.5;\n'
' grad(U) cellMDLimited Gauss linear 0.5;\n}'
)
divSchemes = (
'{\n div(phi,U) Gauss linearUpwind grad(U);\n'
' div(phi,omega) Gauss linearUpwind grad(omega);\n'
' div(phi,k) Gauss linearUpwind grad(k);\n'
' div(phi,e) Gauss linearUpwind grad(e);\n'
' div((nuEff*dev(T(grad(U))))) Gauss linear;\n}'
)
laplacianSchemes = (
'{\n default Gauss linear limited 0.777;\n} '
)
snGradSchemes = (
'{\n default Gauss linear limited 0.777;\n} '
)
blending = 0.7
nonOrthogonalCorrectors = 2
if (maxOrtho) > 0:
gradSchemes = (
'{\n default cellMDLimited Gauss linear 0;\n'
' grad(U) cellMDLimited Gauss linear 0.333;\n}'
)
divSchemes = (
'{\n div(phi,U) Gauss linearUpwind grad(U);\n'
' div(phi,omega) Gauss linearUpwind grad(omega);\n'
' div(phi,k) Gauss linearUpwind grad(k);\n'
' div(phi,e) Gauss linearUpwind grad(e);\n'
' div((nuEff*dev(T(grad(U))))) Gauss linear;\n}'
)
laplacianSchemes = (
'{\n default Gauss linear limited 0.95;\n}'
)
snGradSchemes = (
'{\n default Gauss linear limited 0.95;\n}'
)
blending = 0.8
nonOrthogonalCorrectors = 1
if (LUST):
gradSchemes = (
'{\n default Gauss linear;\n'
' grad(U) cellMDLimited leastSquares 1;\n}'
)
divSchemes = (
'{\n div(phi,U) Gauss LUST grad(U);\n'
' div(phi,omega) Gauss LUST grad(omega);\n'
' div(phi,k) Gauss LUST grad(k);\n'
' div(phi,e) Gauss LUST grad(e);\n'
' div((nuEff*dev(T(grad(U))))) Gauss linear;\n}'
)
laplacianSchemes = (
'{\n default Gauss linear corrected;\n}'
)
snGradSchemes = (
'{\n default Gauss linear corrected;\n}'
)
blending = 0.9
nonOrthogonalCorrectors = 1
if (steady):
ddtSchemes = (
'{\n default steadyState;\n}'
)
if (pseudo_transient):
ddtSchemes = (
'{\n default localEuler;\n}'
)
else:
ddtSchemes = (
'{\n default CrankNicolson ' + str(blending) + ' ;\n}'
)
if precision == 1:
ddtSchemes = (
'{\n default Euler;\n}'
)
if (unbounded):
ddtSchemes = (
'{\n default backward;\n}'
)
wallDist = (
"{\n method meshWave;\n}"
)
#open fvSchemes and write inside
f = open("fvSchemes", "w")
for i in h:
f.write(i + "\n")
f.write("ddtSchemes" + "\n")
f.write(ddtSchemes + "\n")
f.write("\n")
f.write("gradSchemes" + "\n")
f.write(gradSchemes + "\n")
f.write("\n")
f.write("divSchemes" + "\n")
f.write(divSchemes + "\n")
f.write("\n")
f.write("laplacianSchemes" + "\n")
f.write(laplacianSchemes + "\n")
f.write("\n")
f.write("snGradSchemes" + "\n")
f.write(snGradSchemes + "\n")
f.write("\n")
f.write("wallDist" + "\n")
f.write(wallDist + "\n")
f.write(footer)
f.close()
print('File fvSchemes created')
浙公网安备 33010602011771号