发布一个基于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: interPlicFoam
@anubis 试试
gradSchemes
里default Gauss linear;
改成default Gauss pointLinear;
,看看有没有提升。 -
RE: interPlicFoam
@anubis 我们是使用的纯平面,isoAdvector因为使用了iso-surface的概念,使用isoValue去切体单元的每条边后很难保证所有的相交点共面的,所以不得不按照非平面处理。plus版本从2006版本开始也有PLIC算法了,我们自己的PLIC算法植入和他们的区别在于重构算法,最新的版本里我们也使用了自己的时间积分计算方法(正在准备投稿,接收后会开源)。按照目前的测试对比,在溃坝问题里interPlicFOAM和interIsoFoam计算结果非常接近。
-
RE: CCM+到OpenFOAM的网格转换
@王金成
ccmToFoam
本身就在plus版本里,需要你安装libccmio-2.6.1
,编译的时候它会自动检测这个包,检测不到就不会编译了。控制台运行cd $WM_THIRD_PARTY_DIR [ -f libccmio-2.6.1.tar.gz ] || wget ftp://www.daba.lv/pub/TIS/bibliotekas/dazadas/libccmio-2.6.1.tar.gz [ -d libccmio-2.6.1 ] || tar -xzf libccmio-2.6.1.tar.gz ./makeCCMIO lib
source ~/OpenFOAM/OpenFOAM-v2012/etc/bashrc cd $WM_PROJECT_DIR ./Allwmake
ccmToFoam
应该就会安装好了。如果你使用的别的版本就把source ~/OpenFOAM/OpenFOAM-v2012/etc/bashrc
中的版本号替换掉。 -
RE: setFields是否可以在给定的几何域的边界上的网格赋值的时候进行插值
@anubis 把体单元转换为封闭的面网格,使用第三方库求交集,链接中的使用的cork库,效率高但是精度低,不过应该能满足你的要求了。我们有最新的使用CGAL库的版本,而且使用了openmp并行,文章接受后就会开源。