发布一个基于PLIC-VOF
方法的两相流求解器,详见 https://github.com/daidezhi/interPlicFoam
队长别开枪
帖子
-
-
@izumi 关于数值计算和CFD中的松弛技术,写了一点,对你的算例有用没用凑合看看吧
-
@飞火流星jyj 不好意思,这两天太忙没顾得上,刚刚把代码发上来。链接地址:
meshInfoFoam
,它会将体单元的体积、中心和面单元的面积、中心、面积矢量、单位矢量全部写入初始文件夹0
,具体结果可以参考链接meshInfoFoam output
。代码里还测试了一些涉及单个网格元素操作的成员函数。 -
//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;
-
- 对流项
div(rhoPhi, U)
和div(phi, nuTilda)
使用高阶格式; - 网格,做一下网格尺寸无关性实验。
- 对流项
-
@yhdthu
不可压流动中驱动流体流动的是压力梯度,一般对流场中的一个特定位置设定参考位置和参考压力值。如果包含第一类压力边界条件,使用p_rgh进行设置更为方便。而且多相流中,使用p_rgh能使得计算更加稳定和准确,非得使用p的话,就得像fluent使用诸如body-force-weighted等特殊的压力插值格式计算压力梯度了。在不同的流体交界面附近,压力p是连续的,但是压力梯度是非连续的,因为两边的流体密度不一样。
参考高度的不同只会影响p的值,不会影响其梯度的值,因为它们之间相差一个梯度为零的等值压力场。
最后,体积力项包含在p和p_rgh的转换关系里,详细参考东岳流体文档 http://dyfluid.com/interFoam.html ,公式(22)-(25)。 -
@yhdthu
连续方程
动量插值
-
较新版本的Tecplot (例如2015版本)已经支持直接读取OpenFOAM计算结果。
-
if (Pstream::parRun()) { if (Pstream::master()) { Info << "balabala..." << endl; } } else { Info << "balabala..." << endl; }
这样串行并行就都OK了,可能需要
#include "Pstream.H"
。希望能帮到你。 -
@金石为开 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求解器开发 -
@李东岳
我觉得关键点在于进行当前时间步Rhie-chow动量插值时上一时间步界面速度$\mathbf{U}_{f}^{n}$的计算,如果没有ddtcorr()
,$\mathbf{U}_{f}^{n}$只是简单通过相邻单元基于动量方程系数矩阵主对角元素加权平均得到,如果加了修正项,那么$\mathbf{U}_{f}^{n}$使用的就是上一时间步里基于动量插值得到的结果,所以在phiHbyA
的计算里需要修正界面的体积通量,下面的代码是EulerDdtScheme.C
527~545行(OpenFOAM 2.4.0)中的一部分fluxFieldType phiCorr ( phi.oldTime() - (mesh().Sf() & fvc::interpolate(U.oldTime())) );
- (mesh().Sf() & fvc::interpolate(U.oldTime()))
用于消去计算HbyA
时加上的体中心处的旧值。修正的原因是即便体中心处的速度场满足连续方程,基于动量方程系数矩阵主对角元素加权平均往面中心插值,时间步长越小,越倾向于简单算术平均,再加上非结构不规则网格等因素,很难再保证插值后的面速度场仍然保持0散度条件。
个人浅见,欢迎批评指正。 -
@Irisch 你
maxCo
和maxAlphaCo
设置的多少?网格是不是均匀分布的?最小尺寸是多少?代数方程组选择的哪种求解器?收敛残值设置的多少?影响求解时间的因素太多了,还有,按你说的你分了24个区,为啥用96核?建议你把整个算例和log
日志贴出来 -
@史浩
UEqn.H()
是关于U
的一个函数,因此每次迭代更新U
以后也会更新phiHbyA
。 -
@guohuiqun 搞好了,更新到Github页面了。
-
@mohui Splitting 2D and 3D yields two papers.
-
@东岳 不过代数方法最大的问题还是在
sharpness
,界面处的VOF场数值弥散严重,但是优点是够光滑,最近有时在想能不能把两个结合起来。 -
@东岳 可以先用GetData (http://getdata-graph-digitizer.com/download.php) 得到图片里曲线的数据,然后使用excel进行FFT运算 (http://www.stem2.org/je/Excel_FFT_Instructions.pdf)。
-
plus版本里的
ccmToFoam
很好,还支持多域网格的导入,我用的重叠网格都是ccm生成再转OpenFOAM格式的。其实无论hex-core poly mesh
还是general poly mesh
都是使用非结构网格的格式存储的,就数据格式而言没有本质区别。不嫌麻烦的话可以再安装一个plus版本,用它转换网格,再用你自己中意的版本进行计算。其实OpenFOAM各发行版各有优点,有时需要多版本切换,把活干了。
interPlicFoam
OF使用SIMPLE计算10步报错停止,SIMPLEC成功迭代收敛的原因
网格体积
网格体积
openfoam SA模型计算空翼出现Cl和Cd算不准的问题?
关于p_rgh,p,rho*gh的讨论
关于p_rgh,p,rho*gh的讨论
有谁知道怎么在openfoam3.0.1上编译foamToTecplot360?求教。
Info没有输出问题
关于基本求解器中laplacianFoam热传导明明有时间导数为什么还是用稳态的simple算法?
icoFoam的一些细节问题
有关Openfoam在超算并行运算的问题
重新看icoFoam
alphaInitializerFoam
找到vof中interface的位置
请问 vof方法中,surface iso-value为什么取值为0.5
如何从一个图上后处理出周期?鼓泡床直接模拟
讨论下与openfoam匹配的画网格软件