发布一个基于PLIC-VOF
方法的两相流求解器,详见 https://github.com/daidezhi/interPlicFoam

Best posts made by 队长别开枪
-
interPlicFoam
-
RE: OF使用SIMPLE计算10步报错停止,SIMPLEC成功迭代收敛的原因
@izumi 关于数值计算和CFD中的松弛技术,写了一点,对你的算例有用没用凑合看看吧
-
RE: 网格体积
@飞火流星jyj 不好意思,这两天太忙没顾得上,刚刚把代码发上来。链接地址:
meshInfoFoam
,它会将体单元的体积、中心和面单元的面积、中心、面积矢量、单位矢量全部写入初始文件夹0
,具体结果可以参考链接meshInfoFoam output
。代码里还测试了一些涉及单个网格元素操作的成员函数。 -
RE: 网格体积
//class: fvMesh Info << "\n-Class: fvMesh--------" << endl; //- Return the object registry - resolve conflict polyMesh/lduMesh. // Type: virtual const objectRegistry & Info << mesh.thisDb() << endl; //- Return reference to name. // Type: const word & Info << mesh.name() << endl; //- Return reference to boundary mesh. // Type: const fvBoundaryMesh & mesh.boundary(); //- Internal face owner. // Type: const labelUList & labelList owners(mesh.owner()); //- Internal face neighbour. // Type: const labelUList & labelList neighbours(mesh.neighbour()); //- Return cell volumes. // Type: const DimensionedField< scalar, volMesh > & Info << mesh.V() << endl; //- Return old-time cell volumes. // Type: const DimensionedField< scalar, volMesh > & Info << mesh.V0() << endl; //- Return old-old-time cell volumes. // Type: const DimensionedField< scalar, volMesh > & Info << mesh.V00() << endl; //- Return sub-cycle cell volumes. // Type: tmp< DimensionedField< scalar, volMesh > > Info << mesh.Vsc() << endl; //- Return sub-cycle old-time cell volumes. // Type: tmp< DimensionedField< scalar, volMesh > > Info << mesh.Vsc0() << endl; //- Return cell face area vectors. // Type: const surfaceVectorField & Info << mesh.Sf() << endl; //- Return cell face area magnitudes. // Type: const surfaceScalarField & Info << mesh.magSf() << endl; //- Return cell face motion fluxes. // Type: const surfaceScalarField & Info << mesh.phi() << endl; //- Return cell centres as volVectorField. // Type: const volVectorField & Info << mesh.C() << endl; //- Return face centres as surfaceVectorField. // Type: const surfaceVectorField & Info << mesh.Cf() << endl; //- Return face deltas as surfaceVectorField. // Type: tmp< surfaceVectorField > Info << mesh.delta() << endl; Info << "----------------------\n" << endl; // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // //class: fvMesh pointField points(mesh.points()); faceList faces(mesh.faces()); cellList cells(mesh.cells()); // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // //Class: point = vector Info << "\n-Class: point---------" << endl; point &pt(points[3]); Info << "pt = " << pt << endl; //- Return x component Info << "pt.x() = " << pt.x() << endl; //- Return y component Info << "pt.y() = " << pt.y() << endl; //- Return z component Info << "pt.z() = " << pt.z() << endl; Info << "----------------------\n" << endl; // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // //class: edge Info << "\n-Class: edges----------" << endl; edge eg(faces[1].faceEdge(1)), eg_(faces[1].faceEdge(3)); label pt_1(eg.start()), &pt_2(eg.end()); Info << "eg = " << eg << endl; //- Return start vertex label Info << "eg.start() = " << pt_1 << endl; //- Return end vertex label Info << "eg.end() = " << pt_2 << endl; //- Given one vertex, return the other Info << "eg.otherVertex(eg.end()) = " << eg.otherVertex(pt_2) << endl; //- Return common vertex // - -1: no common vertex Info << "eg.commonVertex(eg_) = " << eg.commonVertex(eg_) << endl; //- Return reverse edge Info << "eg.reverseEdge() = " << eg.reverseEdge() << endl; //- Return centre (centroid) Info << "eg.centre(points) = " << eg.centre(points) << endl; //- Return the vector (end - start) Info << "eg.vec(points) = " << eg.vec(points) << endl; //- Return scalar magnitude Info << "eg.mag(points) = " << eg.mag(points) << endl; //- Return edge line Info << "eg.line(points) = " << eg.line(points) << endl; //- compare edges // Returns: // - 0: different // - +1: identical // - -1: same edge, but different orientation Info << "Foam::edge::compare(eg, eg.reverseEdge()) = " << Foam::edge::compare(eg, eg.reverseEdge()) << endl; Info << "----------------------\n" << endl; // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // //class: face Info << "\n-Class: face----------" << endl; face &fe(faces[4]); Info << "fe = " << fe << endl; //- Return true if the face is empty Info << "fe.empty() = " << fe.empty() << endl; //- Return No. of points corresponding to this face Info << "fe.size() = " << fe.size() << endl; //- Return first point Info << "fe.first() = " << fe.first() << endl; //- Return last point Info << "fe.last() = " << fe.last() << endl; //- Return n-th point Info << "fe.operator[](0) = " << fe.operator[](0) << endl; //- Return the points corresponding to this face Info << "fe.points(points) = " << fe.points(points) << endl; //- Centre point of face Info << "fe.centre(points) = " << fe.centre(points) << endl; //- Calculate average value at centroid of face Info << "fe.average(points, points) = " << fe.average(points, points) << endl; //- Magnitude of face area Info << "fe.mag(points) = " << fe.mag(points) << endl; //- Vector normal; magnitude is equal to area of face Info << "fe.normal(points) = " << fe.normal(points) << endl; //- Return face with reverse direction // The starting points of the original and reverse face are identical. Info << "fe.reverseFace() = " << fe.reverseFace() << endl; //- Which vertex on face (face index given a global index) // returns -1 if not found Info << "fe.which(1966) = " << fe.which(1966) << endl; //- Next vertex on face Info << "fe.nextLabel(1) = " << fe.nextLabel(1) << endl; //- Previous vertex on face Info << "fe.prevLabel(1) = " << fe.prevLabel(1) << endl; //- Return number of edges Info << "fe.nEdges() = " << fe.nEdges() << endl; //- Return edges in face point ordering, // i.e. edges()[0] is edge between [0] and [1] Info << "fe.edges() = " << fe.edges() << endl; //- Return n-th face edge Info << "fe.faceEdge(1) = " << fe.faceEdge(1) << endl; //- Return the edge direction on the face // Returns: // - 0: edge not found on the face // - +1: forward (counter-clockwise) on the face // - -1: reverse (clockwise) on the face Info << "fe.edgeDirection(fe.faceEdge(1)) = " << fe.edgeDirection(fe.faceEdge(1)) << endl; //- compare faces // 0: different // +1: identical // -1: same face, but different orientation Info << "Foam::face::compare(fe, fe) = " << Foam::face::compare(fe, fe) << endl; Info << "----------------------\n" << endl; // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // //class: cell Info << "\n-Class: cell----------" << endl; cell &cl(cells[100]); Info << "cl = " << cl << endl; //- Return true if the cell is empty Info << "cl.empty() = " << cl.empty() << endl; //- Return No. of faces corresponding to this cell Info << "cl.size() = " << cl.size() << endl; //- Return first face Info << "cl.first() = " << cl.first() << endl; //- Return last face Info << "cl.last() = " << cl.last() << endl; //- Return n-th face Info << "cl.operator[](0) = " << cl.operator[](0) << endl; //- Return number of faces Info << "cl.nFaces() = " << cl.nFaces() << endl; //- Return labels of cell vertices Info << "cl.labels(faces) = " << cl.labels(faces) << endl; //- Return the cell vertices Info << "cl.points(faces, points) = " << cl.points(faces, points) << endl; //- Return cell edges Info << "cl.edges(faces) = " << cl.edges(faces) << endl; //- Returns cell centre Info << "cl.centre(points, faces) = " << cl.centre(points, faces) << endl; //- Returns cell volume Info << "cl.mag(points, faces) = " << cl.mag(points, faces) << endl; Info << "----------------------\n" << endl;
-
RE: openfoam SA模型计算空翼出现Cl和Cd算不准的问题?
- 对流项
div(rhoPhi, U)
和div(phi, nuTilda)
使用高阶格式; - 网格,做一下网格尺寸无关性实验。
- 对流项
-
RE: 关于p_rgh,p,rho*gh的讨论
@yhdthu
不可压流动中驱动流体流动的是压力梯度,一般对流场中的一个特定位置设定参考位置和参考压力值。如果包含第一类压力边界条件,使用p_rgh进行设置更为方便。而且多相流中,使用p_rgh能使得计算更加稳定和准确,非得使用p的话,就得像fluent使用诸如body-force-weighted等特殊的压力插值格式计算压力梯度了。在不同的流体交界面附近,压力p是连续的,但是压力梯度是非连续的,因为两边的流体密度不一样。
参考高度的不同只会影响p的值,不会影响其梯度的值,因为它们之间相差一个梯度为零的等值压力场。
最后,体积力项包含在p和p_rgh的转换关系里,详细参考东岳流体文档 http://dyfluid.com/interFoam.html ,公式(22)-(25)。 -
RE: 有谁知道怎么在openfoam3.0.1上编译foamToTecplot360?求教。
较新版本的Tecplot (例如2015版本)已经支持直接读取OpenFOAM计算结果。
-
RE: Info没有输出问题
if (Pstream::parRun()) { if (Pstream::master()) { Info << "balabala..." << endl; } } else { Info << "balabala..." << endl; }
这样串行并行就都OK了,可能需要
#include "Pstream.H"
。希望能帮到你。 -
RE: 关于基本求解器中laplacianFoam热传导明明有时间导数为什么还是用稳态的simple算法?
@金石为开 SIMPLE算法本身是为不可压流动中速度-压力的耦合问题设计出来的,这里温度场的求解没有这类耦合问题,用simple只是为了达到按时间推进和设置非正交修正的目的,在源代码里你把
while (simple.loop())
改成while (runTime.loop())
或者while (runTime.run())
效果都是一样的。如果不需要非正交修正,则只需把while (simple.correctNonOrthogonal()) { solve ( fvm::ddt(T) - fvm::laplacian(DT, T) ); }
改成
solve ( fvm::ddt(T) - fvm::laplacian(DT, T) );
此时我们就没有引入simple的必要了,所以这里使用simple只是为了方便设置非正交修正次数罢了。
这里有一个算例,温度场开始随时间剧烈变化,最后趋于稳定
OpenFOAM求解器开发
Latest posts made by 队长别开枪
-
RE: 界面相变
@hongjiewang 在 界面相变 中说:
通过 https://github.com/daidezhi/alphaInitializerFoam 精确初始化相场后,使用
$\alpha ( 1- \alpha )$
的值进行判断,在$[0,1]$区间其值位于$[0,0.25]$之间,为了剔除相分数十分接近0或者1(体单元接近empty或者full,有些单元由于数值计算的误差累计会出现该为empty的时候其相分数是个十分小的数值,full也有这种可能)的情况,加入阈值判断。代码可以参考// Tolerance const scalar tol(1e-7); const volScalarField test(alpha1*(scalar(1)-alpha1)); forAll(test, cellI) { if (mag(test[cellI]) > tol) { // do something ... } else { // do something else ... } }
-
RE: 基于OpenFOAM做软件开发
@谷柏辰 可以,参考 https://engys.com/products/helyx-os ,其他求解器类似,GUI界面设置网格和参数,后台调用求解器并显示在内嵌终端上。