CFD中文网

粒子数量n较少，网格数量m较多时，可以对网格单元建立ADT树，然后确定每个粒子位于哪一个网格，复杂度nlog(m)

@wwzhao

是的,OpenFOAM中面的方向定义永远是从小编号指向大编号，Gauss定理中面的方向指向外。也许把这些面分为大面和小面区分一下更容易理解。

对于编号大于$P$的单元$U$,他们之间面上的量为:

$$\phi_ f = \varpi \phi_{P} + (1 - \varpi) \phi_{U}$$

$$\dot { m } _ { f } = ( \rho \mathbf { u } ) _ { f } \cdot \mathbf { S } _ { f }$$

对于编号小于$P$的单元$L$,他们之间面上的量为:

$$\phi_ f = (1 - \varpi) \phi_{P} + \varpi\phi_{L}$$

$$\dot { m } _ { f } =- ( \rho \mathbf { u } ) _ { f } \cdot \mathbf { S } _ { f }$$

$$
\sum _ { f \sim n b ( P ) } ( \rho \mathbf { u } \phi ) _ { f } \cdot \mathbf { S } _ { f }
= \sum _ { N \in L ( P ) } \left( - \dot { m } _ { f } \phi _ { f } \right) + \sum _ { N \in U ( P ) } \dot { m } _ { f } \phi _ { f }
$$
$$
= \sum _ { N \in L ( P ) } - \dot { m } _ { f } \left[ ( 1 - \varpi ) \phi _ { P } + \varpi \phi _ { N } \right] + \sum _ { N \in U ( P ) } \dot { m } _ { f } \left[ \varpi \phi _ { P } + ( 1 - \varpi ) \phi _ { N } \right]
$$

$$
= \left( \sum _ { N \in L ( P) } - \dot { m } _ { f } ( 1 - \varpi ) + \sum _ { N \in U ( P ) } \dot { m } _ { f } \varpi \right) \phi _ { P } + \sum _ { N \in L ( P ) } - \dot { m } _ { f } \varpi \phi _ { N } + \sum _ { N \in U ( P ) } \dot { m } _ { f } ( 1 - \varpi ) \phi _ { N }
$$

这样可以看出lower存的是 $- \varpi _ { f } \dot { m } _ { f } $, upper存的是$\dot { m } _ { f } ( 1 - \varpi ) _ { f } $.并且diag存的是
$$a _ { P } = -[\sum _ { N \in L ( P ) } \dot { m } _ { f } ( 1 - \varpi ) + \sum _ { N\in U ( P ) } (-\dot { m } _ { f } \varpi )]$$

@dzw05 十分感谢！
这个似乎可以完美应用于空气动力学亚跨超和高超问题的计算。自带AUSM+up全速域通量格式和HLLC格式，有LUSGS预处理的GMRES推进方法。貌似不是Jacobian Free的。还有作者的论文。非常好的资源。
我会慢慢测试几个定常的高超算例，确定比较合适的梯度计算方法、限制器等参数。大佬多多指教~

时隔四年了，还是没有有基于Godunov的近似通量和隐式推进方法的开源OpenFOAM application。
已有的替代解决方案是：
1、使用rhoCentralFoam的KT和KNP格式，本质上还是压力基的求解器。可能存在的问题是耗散大以及收敛性差。
2、使用foam-extend的dbnsFoam，里边有气动人常用的HLLC格式和Roe格式，但推进方法只有龙格库塔显式推进。并且Roe格式熵修正处理不好的话会有Cabuncle现象。
3、使用论文Implementation of density-based solver for all speeds in the framework of OpenFOAM作者公开的allSpeedUnsteadFoam,空间格式使用的是带低速预处理的AUSM系列格式。高超声速热化学非平衡计算中很多人喜欢使用AUSM系列的格式。但allSpeedUnsteadFoam依然只有Runge-Kutta推进。

可以按照论文Implementation of density-based implicit LU-SGS solver in the framework of OpenFOAM 的方法实现LUSGS推进方法，工作量其实很小。如果原作者能开源其lusgsFoam最好了。

不过如果想进一步实现其他隐式推进方法，还是需要搞明白blockLduMatrix，并且组装Jacobian矩阵。我想对于foam-extend的大佬来讲，他们已经实现了blockLduMatrix类，以及后续的代数方程预处理及求解器，如果实现密度基的隐式推进方法应该是易如反掌的。不过等了四年依然没有看到有进展。

空闲时间我会自己着手做这件事情。有人同样在做的话可以一起交流合作~

@wwzhao

扩散方程比较容易理解。这里仅以对流方程为例，不考虑边界条件，不考虑网格扭曲的影响。下面是我的理解。

仅考虑对流项
$$\nabla \cdot ( \rho \mathbf { u } \phi ) = 0$$
将对控制单元$O$的体积分转化为面积分，$ f \sim n b ( O )$为控制单元$O$的控制面，$N \sim NB ( O ) $为控制单元$O$的的邻居单元
$$\sum _ { f \sim n b ( O ) } ( \rho \mathbf { u } \phi ) _ { f } \cdot \mathbf { S } _ { f } = 0$$

定义
$$\dot { m } _ { f } = ( \rho \mathbf { u } ) _ { f } \cdot \mathbf { S } _ { f }$$

面心的物理量通过两侧体心物理量加权平均得到：
$$\phi _ { f } = \varpi \phi _ {O} + ( 1 - \varpi ) \phi _ { N }$$
上面公式整理后为：
$$\sum _ { f \sim n b ( O ) } ( \rho \mathbf { u } \phi ) _ { f } \cdot \mathbf { S } _ { f } = \sum _ { f - n b ( O ) } \dot { m } _ { f } \phi _ { f } = a _ { O } \phi _ { O } + \sum _ { N \sim NB ( O ) } a _ { N } \phi _ { N } = 0$$

其中
$$a _ { N } = \dot { m } _ { f } ( 1 - \varpi ) _ { f }$$
$$a _ { O } = \sum _ { f \sim n b ( O ) } \dot { m } _ { f } \varpi = - \sum _ { N \sim NB ( O ) } a _ { N } + \sum _ { f \sim n b ( O ) } \dot { m } _ { f }$$

在gaussConvectionScheme<Type>::fvmDiv中

    fvm.lower() = -weights.primitiveField()*faceFlux.primitiveField();
    fvm.upper() = fvm.lower() + faceFlux.primitiveField();
    fvm.negSumDiag();

可以看出，lower存的是 $- \varpi _ { f } \dot { m } _ { f } $，代表了以Neighbour为单元的代数方程中，Owner的物理量$\phi _ { O }$的贡献 ,upper 存的是$\dot { m } _ { f } ( 1 - \varpi ) _ { f } $，代表了以Owner为单元的代数方程中，Neighbour的物理量$\phi _ { N }$的贡献 。这里所说的Neighbour和Owner，其实是立足于面的，Owner是该面左侧的单元，Neighbour是该面右侧的单元，Owner编号小于Neighbour编号。

至于weights的具体计算，则根据具体的插值格式进行计算，对于Upwind、FROMM、QUICK，TVD等不同格式，分别有不同的计算方法，这里暂不讨论。

这里有个trick，$a _ { N } = \dot { m } _ { f } ( 1 - \varpi ) _ { f }$，为什么对于编号大于Owner的邻居单元，upper 存的是$\dot { m } _ { f } ( 1 - \varpi ) _ { f } $，而对有编号小于Owner的邻居单元，lower存的是 $- \varpi _ { f } \dot { m } _ { f } $呢?这是因为计算面心的系数$\varpi$时在总是认为面左侧是Owner，右侧是Neighbour.

最后是 fvm.negSumDiag()求和计算对角项的系数$a _ { O } $。这里我们并非按照公式$a _ { O } = \sum _ { f \sim n b ( O ) } \dot { m } _ { f } \varpi$，先对所有的cell循环，然后对每个cell的所有面循环，这样计算效率底下，非结构网格程序一般都是对所有的面进行循环，将其贡献分别计入两侧的控制体中。

    for (label face=0; face<l.size(); face++)
    {
        Diag[l[face]] -= Lower[face];
        Diag[u[face]] -= Upper[face];
    }

这里是对列求和，而非对行求和。如果是对行求和的话，得到的是$- \sum _ { N \sim NB ( O ) } a _ { N }$。而对角项的系数$a _ { O } = - \sum _ { N \sim NB ( O ) } a _ { N } + \sum _ { f \sim n b ( O ) } \dot { m } _ { f }$.

那么最后的问题就是为什么对列求和得到的是正确的$a _ { O } = - \sum _ { N \sim NB ( O ) } a _ { N } + \sum _ { f \sim n b ( O ) } \dot { m } _ { f }$，这是巧合么？并不是巧合，而是必然。仔细想想这些系数所代表的含义，

1. lower存的系数代表了控制面右侧的单元(Neighbour)的代数方程中，控制面左侧单元(Owner)的物理量的贡献
2. upper 存的系数代表了控制面左侧单元(Owner)的代数方程中，控制面右侧的单元(Neighbour)的物理量的贡献

另一方面在$\mathbf A\phi=\mathbf b$中

1. 矩阵$A$中的第$i$行的非对角项系数，其实代表了与编号为$i$的单元相邻的单元物理量对$i$的单元控制方程的贡献。
2. 矩阵$A$中的第$j$列的非对角项系数，其实代表了编号为$j$的单元对其相邻的单元控制方程的贡献。

对于对流项，单元$i$对单元$j$的影响不等于单元$j$对单元$i$的影响，但是，编号为$j$的单元对所有与其相邻的单元的贡献的总和，正应该是$j$的单元自身系数的相反数。也正是矩阵$A$中的第$j$列的非对角项系数之和的相反数。

以上理解也只是我个人根据guassConvectionScheme.C文件的自我解读。如果哪里有理解错的地方，还请指正。

如果通过设计方腔形状参数，能够让压力波动（噪声）幅值变小，是否可以发一篇文章？

这个问题困扰我了一天，但我解决了这个问题。
首先，我们应该明白OpenFOAM中并没有存A[][]这样的系数矩阵，存的是diag，lower和upper这些非零的系数。
其次，我们在看一些教程的时候，为了解释lduAddressing中的lower(), upper(), lowerAddr(), and upperAddr(),可能把lower解释成了
the coefficient multiplying $\phi_{owner} $ in the algebraic equation for element neighbour,把upper解释成了
the coefficient multiplying $\phi_{neighbour} $ in the algebraic equation for element owner.
甚至有些教程给出了下面的解释

A[upperAddr[k]][lowerAddr[k]]=lower[k]
A[lowerAddr[k]][upperAddr[k]]=upper[k]

如果这样理解的话，当然negSumDiag是对列求和的，这和我们认知里矩阵的一行代表一个方程相违背。有人强行解释对于扩散方程因为矩阵对称，虽有对行求和对列求和结果一致。

但事实上出现这样认知的原因在于我们对lower和upper的理解出现的错误。Hrvoje Jasak在帖子中也说明了这点
https://www.cfd-online.com/Forums/openfoam-programming-development/72456-matrix-storage-diagonal-coeffients-calculation.html

经过对guassConvectionScheme的探究及对对流方程离散化的推导（原谅我尚未掌握这个论坛上使用公式排版的方法），正确的理解方式应该是：

1. lower存的系数代表了以Neighbour为单元的代数方程中，Owner的物理量的贡献
2. upper 存的系数代表了以Owner为单元的代数方程中，Neighbour的物理量的贡献

并且根据对流方程离散化的推导，我们其实可以知道diag的系数，并非为该行off-diagonal系数的和。而应该为
$$
a_C=-\sum_{f:nb(C)}{a_N}+\sum_{f:nb(C)}{\dot m_f}
$$
基于规则正交网格的推导可以参考
《The Finite Volume Method in Computational Fluid Dynamics An Advanced Introduction with OpenFOAM and Matlab》

如果有时间，我会重新排版这条回复。有任何问题可以和我联系讨论。