OpenFOAM小代码
-
这些是我感觉对初学者有用的一些代码,收藏在此,持续更新。有问题回帖提问,我会对我贡献的代码负责回答。
计算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]; } }
获得Patch的名字
Info << "patch name is " << U.boundaryField()[patch].patch().name() << nl;
网格相关量
const scalarField& V = mesh.V(); // 网格体积 const surfaceVectorField& Sf = mesh.Sf(); // 网格面矢量 const surfaceScalarField& magSf = mesh.magSf(); // 网格面积的模 const surfaceVectorField& N = Sf/magSf; // 网格面法向
更改壁面一层网格值
label patchID = mesh.boundaryMesh().findPatchID("wall"); const scalarField TfaceCell = T.boundaryField()[patchID].patchInternalField(); k.boundaryFieldRef()[patchID] == sqrt(TfaceCell);
#include "wallFvPatch.H" forAll(nu.boundaryField(),patchi) { if (isA<wallFvPatch>(mesh.boundary()[patchi])) { nu.boundaryFieldRef()[patchi] == 0.5*nu.boundaryField()[patchi].patchInternalField(); } }
给定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
inlet { type codedFixedValue; value uniform (0 0 0); //default value name linearTBC; //name of new BC type code #{ //const fvPatch& boundaryPatch = patch(); //const fvBoundaryMesh& boundaryMesh = boundaryPatch.boundaryMesh(); //const fvMesh& mesh = boundaryMesh.mesh(); //const scalar &t = this->db().time().deltaTValue(); scalar inletS = 0.0; const fvPatch& inletPatch = this->patch(); forAll(inletPatch, faceI) { inletS += inletPatch.magSf()[faceI]; } scalar outletS = 0.0019304; scalar alphain = 0.1; vectorField& vf = *this; forAll(vf, i) { vf[i].y() = 0.01*outletS/inletS/alphain; vf[i].x() = 0.0; vf[i].z() = 0.0; } #}; }
努塞尔数
下列代码可以添加到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字典文件
particles { libs ("liblagrangianFunctionObjects.so"); type particles; }
在流场中添加欧拉标量传输
controlDict添加下述内容,同时需要在0文件夹下添加
s
标量场,以及fvScheme、fvSolution都要做适配functions { scalarTransport1 { type scalarTransport; libs ("libsolverFunctionObjects.so"); field s; D 1e-4; fvOptions { s { type semiImplicitSource; active true; selectionMode points; points ( //(0.57 0.001 0.325)//w/h=5 //(0.62 0.001 0.325)//w/h=5 (0.67 0.001 0.325)//w/h=5 //(0.72 0.001 0.325)//w/h=5 //(0.77 0.001 0.325)//w/h=5 ); volumeMode absolute; sources { s { explicit 0.0001; implicit 0; } } } } } }
对网格进行赋值
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) ) ); }
OpenFOAM的二维数组
PtrList<PtrList<scalar>> te(4); forAll(te, i) { te.set(i, new PtrList<scalar>(4)); forAll(te[i], j) { te[i].set(j, new scalar (0.0)); } }
输出最大值最小值
functions { minMaxp { type fieldMinMax; functionObjectLibs ("libfieldFunctionObjects.so"); fields ( U ); location no; writeControl timeStep; writeInterval 1; } }
dimensionedScalar rho(dimDensity, 1.0); dimensionedScalar U(dimVelocity, 1.0); dimensionedScalar Ain(dimArea, 0.000481); dimensionedScalar sin(dimless, 1); dimensionedScalar min = rho*U*Ain*sin; tmp<volScalarField> tDEF ( volScalarField::New ( type(), mesh_, dimensionedScalar(dimless, 0) ) ); tmp<volScalarField> mwj ( volScalarField::New ( type(), mesh_, dimensionedScalar(dimMass/dimTime, 0) ) ); const volScalarField& s = lookupObject<volScalarField>("s"); const volScalarField::Boundary& sBf = s.boundaryField(); volScalarField::Boundary& mwjBf = mwj.ref().boundaryFieldRef(); scalar mw = 0.0; forAll(mesh_.boundary(), patchi) { if (mesh_.boundary()[patchi].name() == "wall") { scalarField Aj = mesh_.boundary()[patchi].magSf(); mwjBf[patchi] = -rho.value()*Aj*sBf[patchi].snGrad()*1e-5; forAll(mesh_.boundaryMesh()[patchi], facei) { mw += mwjBf[patchi][facei]; } } } volScalarField::Boundary& DEFBf = tDEF.ref().boundaryFieldRef(); forAll(mesh_.boundary(), patchi) { if (mesh_.boundary()[patchi].name() == "wall") { scalar area = gSum(mesh_.magSf().boundaryField()[patchi]); Info<< "Wall area is " << area << nl; scalarField Aj = mesh_.boundary()[patchi].magSf(); //DEFBf[patchi] = mwjBf[patchi]/Aj/(mw/A); DEFBf[patchi] = mwjBf[patchi]/Aj; } } return tDEF;
单独一个文件输出直径信息
meanDiameter { type coded; libs ("libutilityFunctionObjects.so"); name error; codeExecute #{ const volScalarField& d = mesh().lookupObject<volScalarField>("d.alpha.oil"); scalar d32 = d.weightedAverage(mesh().V()).value(); if (Pstream::master()) { std::ofstream file; file.open ("d32", std::ofstream::out | std::ofstream::app); file << mesh().time().timeName() << " " << d32 << "\n"; file.close(); } #}; }
拉格朗日输出alpha
cloudFunctions { f1 { type voidFraction; } }
-
补充一个CodeStream实现充分发展的湍流入口,请东岳老师不要建议
$\begin{array}{l}
\delta_{T}>\delta_{w} \rightarrow u=u_{\infty}\left(\frac{d_{w}}{\delta_{T}}\right)^{1 / 7} \
\delta_{T}<\delta_{w} \rightarrow u=u_{\infty}
\end{array}
\delta_{T}=0.37\left(v / u_{\infty}\right)^{0.2} L_{T}^{0.8}
$
dw为到壁面的最短距离OpenFoam 7:
gasinlet { type fixedValue; value #codeStream { codeInclude #{ #include "fvCFD.H" #include "DynamicList.H" #}; codeOptions #{ -I$(LIB_SRC)/finiteVolume/lnInclude \ -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); #}; }; }
-
class base { private: bool ownedByRegistry_; public: template<class Type> static void store(Type* p) { //p->base::ownedByRegistry_ =true; p->ownedByRegistry_ = true; p->printInfo(); } }; class derived :public base { public: derived() {} void printInfo() { Info << "derived class" << endl; } }; int main(int argc, char *argv[]) { ... derived d; base::store(&d); ... }
-
@东岳 在 OpenFOAM小代码 中说:
东岳老师,请教几个小问题:
- 在编译过程中报错:
'mesh_' was not declared
,所以mesh_
是不可以直接使用的,但是我不太理解什么是类内可以调用的成员变量,感觉速度场U,温度场T等都是类内可以调用的成员变量? - 如果
mesh_
不可获取,U_
可获取。上面说mesh_
是类内可以调用的成员变量,然而U_
又是什么? - 声明一个tmp变量,
mesh_
表示什么意思呢?在声明U、T、p等场的时候,mesh_
的位置常常写的是mesh
期待回复,祝好!
- 在编译过程中报错:
-
@Zhong-combustion 你好,我目前也在用LES做加热管流的湍流研究,关于湍流入口的设置看了一些资料,如有兴趣,可以加个qq:1191364394,共同学习探讨。
-
@东岳 东岳老师,没能找出报错原因 ,只能再来麻烦一下您:
INLET1 { type codedFixedValue; value $internalField; name stepInjection; redictType variedValue; code #{ scalar alphapar = this->patch().fluid().phase(); scalar t = this->db().time().value(); if(t >= 0 && t<=5) { alphapar = 0; } else { alphapar = 0.15; } operator==(alphapar) #}; }
-
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')
-
@李东岳 在 OpenFOAM小代码 中说:
IOField<scalar> utau
(
IOobject
(
"utau",
runTime.constant(),
"../postProcessing",
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
scalarField(totalFSize,0.0)
);请问老师,如果我想要在postProcessing中输出一个标量,他只是一个数,并不是场量,我应该怎么定义他的类型。上面这种方法是不是只适应于场量的输出。我想输出的是下面c的数值。
const volScalarField& b = mesh().lookupObject<volScalarField>("alpha.liquid"); scalar c= b.weightedAverage(mesh().V()).value();
-
IOList<scalar> utau ( IOobject ( "utau", runTime.constant(), "../postProcessing", mesh, IOobject::NO_READ, IOobject::AUTO_WRITE ), 1 ); utau[0] = b.weightedAverage(mesh().V()).value();
-
totalLiquid { libs (utilityFunctionObjects); type coded; name totalLiquid; enabled true; writeControl timeStep; writeInterval 1; codeOptions #{ -I$(LIB_SRC)/meshTools/lnInclude #}; codeExecute #{ const volScalarField& b = mesh().lookupObject<volScalarField>("alpha.liquid"); IOList<scalar> liquidFraction ( IOobject ( "liquidFraction", mesh().time().constant(), "../postProcessing", mesh(), IOobject::NO_READ, IOobject::AUTO_WRITE ), 1 ); liquidFraction[0] = b.weightedAverage(mesh().V()).value(); #}; }
东岳老师,我这样添加到controlDict里,运行后postProcessing里并没有出现liquidFraction的值。
-
@hongjiewang 在 OpenFOAM小代码 中说:
totalLiquid { libs (utilityFunctionObjects); type coded; name totalLiquid; enabled true; writeControl timeStep; writeInterval 1; codeOptions #{ -I$(LIB_SRC)/meshTools/lnInclude #}; codeExecute #{ const volScalarField& b = mesh().lookupObject<volScalarField>("alpha.liquid"); IOList<scalar> liquidFraction ( IOobject ( "liquidFraction", mesh().time().constant(), "../postProcessing", mesh(), IOobject::NO_READ, IOobject::AUTO_WRITE ), 1 ); liquidFraction[0] = b.weightedAverage(mesh().V()).value(); #}; }
东岳老师,我这样添加到controlDict里,运行后postProcessing里并没有出现liquidFraction的值。
需要修改两个部分,
1.在下方添加 .write()
2.将codeExecute改为codeWrite
之后就可以得到想要的相含量结果。