@random_ran 您好，想咨询一下，这么多的网格您是如何进行后处理的呢？paraview 会崩溃的吧？

@yu_tian 很简单，你在decomposPar 之前renumber, 然后把生成的网格复制到constant 中替换，再进行分块和并行计算。

@李东岳 非常抱歉，这个硬盘我还真不知道。我用的是UCL的集群。

@李东岳 我看教程有pvserver 可以并行。但是有时候不知道学校的管理员是不是给封了，以前速度还可以，现在用起来就很揪心。集群上模拟最大的问题就是网格多了以后数据文件都比较大都是100G+，没个4T的硬盘真不太敢做模拟。后处理和下载数据都是心病。

目前使用OpenFoam 在集群上做模拟，网格数量级在千万级别。每个时间步的文件都非常大，对于这样的情况，如何能够有效的进行后处理，希望各位朋友提供帮助和分享一下经验。

@samuel-tu 为什么感觉PIMPLE才是真正的PISO算法，而PISO算法只是进行了多次压力和速度修正，在OpenFOAM中并没有判断一个时间步内收敛与否。从流程图上看也没有区别，https://marinecfd.xyz/post/piso-algorithm/ 这个帖子还对PISO的细节进行了分析。

@李东岳 我下面仔细推导了一下：
首先 rAU(1.0/UEqn.A());，这个公式的计算得到的结果是：

对角系数

这是A()函数：

 template<class Type>
 Foam::tmp<Foam::volScalarField> Foam::fvMatrix<Type>::A() const
 {
     tmp<volScalarField> tAphi
     (
         volScalarField::New
         (
             "A("+psi_.name()+')',
             psi_.mesh(),
             dimensions_/psi_.dimensions()/dimVol,
             extrapolatedCalculatedFvPatchScalarField::typeName
         )
     );
 
     tAphi.ref().primitiveFieldRef() = D()/psi_.mesh().V();
     tAphi.ref().correctBoundaryConditions();
 
     return tAphi;
 }
 

这是D()对角系数的平均化处理

 template<class Type>
 Foam::tmp<Foam::scalarField> Foam::fvMatrix<Type>::D() const
 {
     tmp<scalarField> tdiag(new scalarField(diag()));
     addCmptAvBoundaryDiag(tdiag.ref());
     return tdiag;
 }

边界对对角系数的影响：

 template<class Type>
 void Foam::fvMatrix<Type>::addCmptAvBoundaryDiag(scalarField& diag) const
 {
     forAll(internalCoeffs_, patchi)
     {
         addToInternalField
         (
             lduAddr().patchAddr(patchi),
             cmptAv(internalCoeffs_[patchi]),
             diag
         );
     }
 }

从上面的代码可以得到：
$$
A_p=\frac{\bar{D}}{\Delta V}
$$
其中：
$$
\bar{D}=average(a_p')+diag
$$
式子中average(Ap')表示的是边界对主对角线系数的影响的平均；diag表示的是内部面离散的主对角线系数。

周围系数作为源项：

H()

 template<class Type>
 Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>>
 Foam::fvMatrix<Type>::H() const
 {
     tmp<GeometricField<Type, fvPatchField, volMesh>> tHphi
     (
         GeometricField<Type, fvPatchField, volMesh>::New
         (
             "H("+psi_.name()+')',
             psi_.mesh(),
             dimensions_/dimVol,
             extrapolatedCalculatedFvPatchScalarField::typeName
         )
     );
     GeometricField<Type, fvPatchField, volMesh>& Hphi = tHphi.ref();
 
     // Loop over field components
     for (direction cmpt=0; cmpt<Type::nComponents; cmpt++)
     {
         scalarField psiCmpt(psi_.primitiveField().component(cmpt));
 
         scalarField boundaryDiagCmpt(psi_.size(), 0.0);
         addBoundaryDiag(boundaryDiagCmpt, cmpt);
         boundaryDiagCmpt.negate();
         addCmptAvBoundaryDiag(boundaryDiagCmpt);
 
         Hphi.primitiveFieldRef().replace(cmpt, boundaryDiagCmpt*psiCmpt);
     }
 
     Hphi.primitiveFieldRef() += lduMatrix::H(psi_.primitiveField()) + source_;
     addBoundarySource(Hphi.primitiveFieldRef());
 
     Hphi.primitiveFieldRef() /= psi_.mesh().V();
     Hphi.correctBoundaryConditions();
 
     typename Type::labelType validComponents
     (
         psi_.mesh().template validComponents<Type>()
     );
 
     for (direction cmpt=0; cmpt<Type::nComponents; cmpt++)
     {
         if (validComponents[cmpt] == -1)
         {
             Hphi.replace
             (
                 cmpt,
                 dimensionedScalar(Hphi.dimensions(), 0)
             );
         }
     }
 
     return tHphi;
 }

 template<class Type>
 void Foam::fvMatrix<Type>::addBoundaryDiag
 (
     scalarField& diag,
     const direction solveCmpt
 ) const
 {
     forAll(internalCoeffs_, patchi)
     {
         addToInternalField
         (
             lduAddr().patchAddr(patchi),
             internalCoeffs_[patchi].component(solveCmpt),
             diag
         );
     }
 }

这里面的求解过程包含了：
$$
H=\frac{\left [-\sum a_{\mathbf{N}}\mathbf{U_N} + (\mathbf{b}+\mathbf{b'})+(average(a_p')-a_p')\mathbf{U_C}\right]}{\Delta V}
$$

最后：

$$
HbyA= \frac{H}{A}=\frac{\Delta V}{average(a_p')+diag}\frac{\left [-\sum a_{\mathbf{N}}\mathbf{U_N} + (\mathbf{b}+\mathbf{b'})+(average(a_p')-a_p')\mathbf{U_C}\right]}{\Delta V}
$$
$$
HbyA= \frac{H}{A}=\frac{\left [-\sum a_{\mathbf{N}}\mathbf{U_N} + (\mathbf{b}+\mathbf{b'})+(average(a_p')-a_p')\mathbf{U_C}\right]}{average(a_p')+diag}
$$

@李东岳 如果只是离散内部面是相同的，但是将边界的影响引入到了对角线上的Ap系数后，对角线系数就不一样了。比如Wall边界条件，里面涉及到了应力的分解：

@李东岳 李老师，您好。volScalarField rAU(1.0/UEqn.A());这个重理论上分析应该是矢量才对，但是细节上的操作确实做成了体标量场。

 template<class Type>
 Foam::tmp<Foam::volScalarField> Foam::fvMatrix<Type>::A() const
 {
     tmp<volScalarField> tAphi
     (
         volScalarField::New
         (
             "A("+psi_.name()+')',
             psi_.mesh(),
             dimensions_/psi_.dimensions()/dimVol,
             extrapolatedCalculatedFvPatchScalarField::typeName
         )
     );
 
     tAphi.ref().primitiveFieldRef() = D()/psi_.mesh().V();
     tAphi.ref().correctBoundaryConditions();
 
     return tAphi;
 }

$ A=\frac{A_p}{\Delta V}$
在计算系数Diag的过程中，将边界对主对角线的系数进行了平均更新：

 template<class Type>
 Foam::tmp<Foam::scalarField> Foam::fvMatrix<Type>::D() const
 {
     tmp<scalarField> tdiag(new scalarField(diag()));
     //这里进行了平均处理并修正了主对角线系数。
     addCmptAvBoundaryDiag(tdiag.ref());
     return tdiag;
 }

为什么这样操作呢？这样会使得系数矩阵不一样了!

volScalarField rAU(1.0/UEqn.A());

按照道理说这个应该是矢量场才对，因为三个方向volVectorField。

winsway_zero@foxmail.com

最近在学习李老师刚发表的文章：twoWayGPBEFoam: An open-source Eulerian QBMM solver for monokinetic bubbly flows。文章后面附有FOAM的code,不过在编译的时候出现了一个小问题。

就是fvMesh中没有solver这个成员函数（openFOAM-7.0），希望能够得到李老师帮助，根据code学习一下该算法。谢谢，谢谢~

你可以读取openFOAM处理好的网格文件啊，它的网格文件还是很有规律的，至少我现在已经实现读取和处理openFOAM的网格数据了，自己在编写求解器。

大家好，我现在用FLuent 做一个二维i微通道里液滴运动的模拟，连续相为油，液滴为水。在初始化时液滴通道里patch了一个水的区域，连续相初始有速度。使用VOF 模型计算，但是在最终计算的过程中，发现液滴的大小被压缩或者扩张了。【像下面图像展示的那样，液滴逐渐被拉长了，似乎不再满足质量守恒】实际情况下不应该出现这种现象的。请问大家这是什么原因呢？

@TINGHAIK 完全可以，没有问题的。

gmsh可以，这个是开源的网格软件。

@东岳 谢谢东岳老师哈

@Samuel-Tu 大佬，期待你搞一个附带公式的流程图

@东岳
主要是这个哈

@东岳 岳哥，帮忙看看呗，我写的公式是否正确，特别是梯度的算数平均那里哈~