CFD中文网

@李东岳 在 势流算法问题 中说：

推导太长没细看，如果都是1 -2 1，那你的矩阵不相容无法求解。

老师，我推导出来的系数矩阵是这个，上面是我写错了
\begin{equation}
\begin{bmatrix} 1& -1& 0& 0\-1& 2& -1& 0\ 0& -1& 2&-1 \ 0& 0& -1& 1\ \end{bmatrix}
\end{equation}

@李东岳 在 势流算法问题 中说：

推导太长没细看，如果都是1 -2 1，那你的矩阵不相容无法求解。

我和您推导的结果最大的区别在于系数矩阵的那个3和源项，这个3是咋求出来的呢？

@李东岳 在 势流算法问题 中说：

网格2和3，离散出来的系数为啥不一样

老师，是我笔误，2和3网格系数都是1 -2 1，推导过程没错。

@李东岳 在 势流算法问题 中说：

你得矩阵看起来不相容。你用的都是零梯度边界条件么

第一个矩阵是用进口固定梯度，出口固定值条件求出来的，第二个矩阵是用入口零梯度，出口固定值，加上一个源项求出来的。求出来这俩也不等价啊

@李东岳 在 势流算法问题 中说：

哪个第一个图？

另外，按照 http://dyfluid.com/potentialFoam.html 上面写的，目前我们要证明的是，链接里面的方程10（进口固定梯度，出口固定值），与方程11（进口零法向梯度，出口固定值，附加源项）的等价。考虑4网格系统，方程10与11离散出来的矩阵都是：

\begin{equation}
\begin{bmatrix}
1 & -1 & & \\
-1 & 2 & -1 & \\
& -1 & 2 & -1 \\
& & -1 & 3 \\
\end{bmatrix}
\begin{bmatrix}
\Phi_1 \\
\Phi_2 \\
\Phi_3 \\
\Phi_4 \\
\end{bmatrix} \equiv
\begin{bmatrix}
u_{inlet} \\
0 \\
0 \\
v \\
\end{bmatrix}
\end{equation}

老师，方程(11)和进口固定梯度，出口固定值得到的矩阵我又手算了一遍。对于4网格系统:

1. 进口固定梯度，出口固定值

边界条件:

cell 1:

cell 2:

同理cell 3:

cell 4:

矩阵为:

2. 方程(11)

边界条件:

cell 1:

cell 2:

同理cell 3:

cell 4:

矩阵为:

和您算的系数矩阵第四行不一致，和等号右边的第一行和第四行也不一致，但计算过程我检查了好几遍，也没发现问题，问题出在哪儿了？

谢谢岳子老师。

@Afr1yne
1.我知道是湍流模型，我想在openfoam的代码中找到它
3.我了解过FvVectorMatrix，说它存储了Ax=b的相关信息。我也能猜到不是连等，我只是想看看最底层的定义，但我不知道在哪里找。
在你想标红的字前后加上`

@Afr1yne 在 PISO算法讨论 中说：

@yjc95

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

谢谢您的回复。
1.“pisoFoam用湍流模型把粘性项模化turbulence->divDevSigma(U)”具体是咋模化的呢？
2.“fvModels.source(U)是源项”是哪个源项？是体积力、电场力这种？
3.假设我打开了动量预测，这俩==该咋理解呢？

@李东岳
@李东岳
老师，这是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++语法特点吗？

@yjc95 在 PISO算法讨论 中说：

@李东岳

1. 方程(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 在 PISO算法讨论 中说：

@李东岳

1. 方程(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，因为我把这一行注释之后报错了。

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

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

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

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

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

好的

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

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

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

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

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

我就想问该怎么化简呀0.0

@李东岳

1. 方程(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的解释。

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

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)是啥意思呢？

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

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

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