PISO算法讨论

李东岳

122页画红线那个，我不觉得是他写的那样。OpenFOAM里面的PISO同样略去了临点。

yjc95

@李东岳老师，这就是我想请教您的。我看您在用户手册上写到过，您跟的Issa讨论过PISO算法，那么openfoam中的PISO是Issa提出的那样吗？还是PISO算法发展到现在已经有很多变体了。

李东岳

我之前见过weller说过：至少rhoPimpleFoam这个求解器，是Henry Weller完全自己推出来的。没参考任何资料。

我看过1985年原生的PISO算法，现在有点忘了，我只记得OpenFOAM的pisoFoam的算法。大略看一眼的话，OpenFOAM跟原生的PISO有区别，你看原生的PISO湍流模型就知道。

李东岳

如果这本书是你买的，你可以给他发邮件问一下他书里面说的，OpenFOAM中的pisoFoam，哪里体现了并没有略去临点的思想。

yjc95

@李东岳好的老师。再向您请教个问题，icoFoam里面，fvSolution下定义了pRefCell和pRefValue，在求解方程中用到的p是已经减去这个值的吗？

yjc95

@yjc95 老师，我看源代码解决这个问题了，是先减去参考点的值再求解压力方程。

yjc95

@李东岳老师，还有几个问题向您请教。

1.

$ \triangledown p $离散后为啥要除以体积？

2.

这三段代码对应的应该是（3）、（4）、（5）都除以$ \Delta t$ 吧？

3.

icoFoam.C源代码中，

        if (piso.momentumPredictor())
        {
            solve(UEqn == -fvc::grad(p));
        }

// Non-orthogonal pressure corrector loop
            while (piso.correctNonOrthogonal())

中piso.momentumPredictor()和piso.correctNonOrthogonal()是在哪个文件中定义的？constant和system两个文件夹下没找到。

4.

icoFoam.C源代码中，HbyA的定义为啥是HbyA(constrainHbyA(rAU*UEqn.H(), U, p))，constrainHbyA啥意思？

5.

组建phiHbyA用到的 fvc::flux(HbyA) + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)是啥意思呢？

李东岳

方程就是这样离散的
已更新
fvSolution
参考无痛苦笔记
ddtCorr确保时间步更小的时候结果正确，算法我还没看

yjc95

@李东岳

方程(27)等号右边等于
\begin{equation}
\Delta t \int_{}^{}\frac{1}{A}\bigtriangledown pd\textbf{S}
\end{equation}

\begin{equation}
=\Delta t\sum \frac{1}{A_{f}}(\bigtriangledown p)_{f}\cdot \textbf{S}
\end{equation}
面上压力梯度咋就变成了方程(26)矩形框的那部分？

3.cavity案例的fvSolution里p-preconditioner DIC对应piso.momentumPredictor，因为有DIC这一设置，所以if为真，执行动量预测步。同样在fvSolution里，PISO-nNonOrthogonalCorrectors 0对应的应该是piso.correctNonOrthogonal，但值为0，while语句不执行，就不会求解压力方程了。但运行该案例，是有求解压力的。

4.无痛苦笔记中没有关于constrainHbyA的解释。

李东岳

你得方程2，里面的$\nabla p$还需要继续化简

yjc95

@李东岳在 PISO算法讨论中说：

你得方程2，里面的$\nabla p$还需要继续化简

我就想问该怎么化简呀0.0

李东岳

化简成26行括号里面的形式

yjc95

@李东岳在 PISO算法讨论中说：

化简成26行括号里面的形式

这样化简的原理是啥呀？没想明白

李东岳

$\nabla p=\frac{1}{V}\int \nabla p \rd V$ 是基本的离散法则。我向你推荐算法课 http://dyfluid.com/class.html

yjc95

@李东岳在 PISO算法讨论中说：

$\nabla p=\frac{1}{V}\int \nabla p \rd V$ 是基本的离散法则。我向你推荐算法课 http://dyfluid.com/class.html

好的

yjc95

@李东岳在 PISO算法讨论中说：

$\nabla p=\frac{1}{V}\int \nabla p \rd V$ 是基本的离散法则。我向你推荐算法课 http://dyfluid.com/class.html

老师，我突然想明白了，$\bigtriangledown p$是定义在体心上的，当然就等于控制体内各处压力梯度的值之和除以体积了。脑子一下子短路了(T＿T)

yjc95

@yjc95 在 PISO算法讨论中说：

@李东岳

方程(27)等号右边等于
\begin{equation}
\Delta t \int_{}^{}\frac{1}{A}\bigtriangledown pd\textbf{S}
\end{equation}

\begin{equation}
=\Delta t\sum \frac{1}{A_{f}}(\bigtriangledown p)_{f}\cdot \textbf{S}
\end{equation}
面上压力梯度咋就变成了方程(26)矩形框的那部分？

3.cavity案例的fvSolution里p-preconditioner DIC对应piso.momentumPredictor，因为有DIC这一设置，所以if为真，执行动量预测步。同样在fvSolution里，PISO-nNonOrthogonalCorrectors 0对应的应该是piso.correctNonOrthogonal，但值为0，while语句不执行，就不会求解压力方程了。但运行该案例，是有求解压力的。

4.无痛苦笔记中没有关于constrainHbyA的解释。

老师，问题3，piso.correctNonOrthogonal对应的是fvSolution下的PISO-nCorrectors，因为我把这一行注释之后报错了。

yjc95

@yjc95 在 PISO算法讨论中说：

@李东岳

方程(27)等号右边等于
\begin{equation}
\Delta t \int_{}^{}\frac{1}{A}\bigtriangledown pd\textbf{S}
\end{equation}

\begin{equation}
=\Delta t\sum \frac{1}{A_{f}}(\bigtriangledown p)_{f}\cdot \textbf{S}
\end{equation}
面上压力梯度咋就变成了方程(26)矩形框的那部分？

3.cavity案例的fvSolution里p-preconditioner DIC对应piso.momentumPredictor，因为有DIC这一设置，所以if为真，执行动量预测步。同样在fvSolution里，PISO-nNonOrthogonalCorrectors 0对应的应该是piso.correctNonOrthogonal，但值为0，while语句不执行，就不会求解压力方程了。但运行该案例，是有求解压力的。

4.无痛苦笔记中没有关于constrainHbyA的解释。

老师，问题4，我在您最新的无痛苦笔记里找到了

yjc95

@李东岳
 @李东岳
老师，这是pisoFoam里的速度方程

    fvm::ddt(U) + fvm::div(phi , U)
+  MRF.DDt(U)
+  turbulence->divDevSigma(U)
==
    fvModels.source(U)

solve(UEqn == -fvc::grad(p))

这是icoFoam里的速度方程

    fvm::ddt(U) 
+ fvm::div(phi , U)
-  fvm::laplacian(nu , U)

solve(UEqn == -fvc::grad(p))

比较不难发现，pisoFoam里的速度方程多了MRF.DDt(U)，turbulence->divDevSigma(U)
和== fvModels.source(U)，少了fvm::laplacian(nu , U)。

MRF.DDt(U)是和旋转结构有关，在cavityRAS这个案例中用不到。

1.

turbulence->divDevSigma(U)应该是这个案例中用到的K-epsilon方程，这个方程的具体代码在哪个文件里可以看到呢？fvModels.source(U)就是fvm::laplacian(nu , U)吗？

2.

第一段代码中已经有==了，第二段代码又有一个==，从逻辑上讲，应该将fvModels.source(U)移到等号左边，再等于-fvc::grad(p)，这是C++语法特点吗？

Afr1yne

@yjc95

icoFoam是不能用湍流模型的，fvm::laplacian(nu , U)是动量方程右边的粘性项。pisoFoam用湍流模型把粘性项模化turbulence->divDevSigma(U)。fvModels.source(U)是源项。
2.solve(UEqn == -fvc::grad(p))是momentumPredictor里面的，不开的话是用不上的，我感觉你是不是把它理解成求解了？

CFD中文网

PISO算法讨论

1.

2.

3.

4.

5.

1.

2.