基于FVM有限体积算法的非结构化网格稳态通量算法原理及代码实现

FVM算法原理

本文不涉及simple算法，主要是基于有限体积方法对非结构化网格通量的计算，主要参考书籍为The Finite Volume Method in Computational Fluid Dynamics。
在前文中我们介绍了CV-FEM算法实现，这里针对FVM算法介绍，希望能够抛砖引玉，与大家共同学习。
本文涉及的内容为非结构化网格的稳态计算。
FVM基本思想：通过将连续的结构离散为一个个小的单元，并将每一个单元视为一个控制体积，保证在控制体积内满足动量、能量、质量的严格守恒，并将通量、梯度等物理信息储存在单元质心处，通过插值计算内部面的物理信息，通过合理的设置边界条件实现模拟各种流动现象。

扩散项处理

NS方程如式所示，首先我们介绍关于扩散项diffusion term的处理，在非结构化网格结构中，由于相邻网格质心的连线与共享面的法向量具有一定的夹角误差，因此不能达到结构化网格插值计算的二阶精度，需要考虑非正交扩散的影响。

其中，(▽φ)f 表示在单元共享面的梯度值，Ef 是由面矢量Sf分解而来，即Ef + Tf = Sf；这里的Sf大小为面的面积。

上式将扩散项离散为正交贡献、非正交贡献，也称交叉扩散项，有如下性质：

1. 正交贡献项与结构化网格类似，可以通过中心差分计算共享面的梯度▽φ；
2. 交叉扩散项不能直接通过中心差分进行离散，该项需要进行延迟处理，即是利用前一次求出的全场通量φ数值求解梯度，然后作为方程源项迭代求解。

正交贡献项的计算方法

我们可以看到，除Ef以外的量都可以根据单元的几何信息获得，因此该项的计算等价于Ef的确定。

1. 最小校正法
   我们令Ef与Tf正交，使非正交项尽可能的小，则可以得到
   
2. 正交校正法
   这种方法使得无论网格的非正交程度如何，涉及到项的贡献保持与正交网格一致
   
3. 过松弛方法
   这种方法被证明是最稳定的方法，令Tf与Sf正交，当非正交性增大时，涉及到的φC和φF的重要性增加
   

交叉扩散项的计算

对于Tf，我们可以利用Sf - Ef计算得到，关键在于▽φ的计算方法，梯度计算的数值方法有多种，我们这里介绍Green-Gauss(格林-高斯)方法。
本文略去推导过程，数值计算公式为

该式表示，控制体积质心处的梯度等价于单元各面的通量与面矢量的点积/单元体积，我们首先计算共享面的φ值，然后可以通过该式计算得到全场单元质心处的梯度值，最后再去求单元间共享面的梯度值：

其中gc和gf是单元间的插值因子，利用线性插值，有：

其中，f1表示单元表面。
对于梯度计算，同样需要考虑非正交的影响，这里我们利用式所示方法消除非正交影响：

上式是对第一次计算出的结果进行了迭代计算，通过迭代计算消除非正交影响，通常迭代不超过两次。
最后，扩散项可以离散为

需要注意的是：上述项的计算是建立在得到初解的基础上，通过已有的全场φ值迭代计算全场的交叉扩散项，然后放入方程源项，重新求解方程组。

对流项

略去离散过程，对流项可以离散为

其中，mf是质量流量，此处采取了一阶迎风格式(Upwind)。这里 mf 可由式计算
$$m_f=pv{S}_f$$
p表示密度，v表示速度，Sf是面矢量。

边界条件

1. 固定值边界，我们需要在离散方程中调整ac项与bc项，引入

这里，FluxCb、FluxVb是由于边界产生的项。并且FluxVb中的第二项与交叉扩散项类似，也需要按照前文方法进行计算，并延迟处理。

2. 通量边界
   对于通量边界，只需要在源项加入通量影响
   $$b_c=b_c-q{S}_b$$
   其中，q是设置的通量，Sb是边界面的法向量。

代码实现

这里我们的FVM算法继承了前文设置的CVFEM算法的一些函数和信息，可以参阅前文。
下面对代码的数据结构做一些解释：

	class Fvm :public Mesh_cvfem {
	public:
		void findsharedface();
		void init(string cwd);
		void cal(string cwd);
		void cal_Diff(string cwd);
		void cal_Convertion(string cwd);
		void cal_gradient1(vector<double>&, vector<double>&);
		void cal_gradient(vector<double>&, vector<double>&);
		void insert_fai(double, uint64_t);
		//压力速度耦合算法
		void simple();
		void read(string cwd);
		void marge_msh(vector<string>filepath);
		double error_Fvm(vector<double>&, vector<double>&);
		void Write(vector<double>& x, string cwd);
		inline uint64_t faceid(uint64_t, uint64_t, uint64_t);
		//virtual void write() override;
		//
	private:
		vector<double> solver_equtionGaussSeidel(vector<vector<double>>& A, vector<double>& b) override;//高斯赛德尔迭代
		void pgrad(Point&, int, int, Point&, Point&, double);
		//vector<Element_cvfem>mesh_eles;//单元
		//vector<Point>mesh_nodes;//节点
		//vector<vector<int>>mesh_nei;//邻居
		//vector<vector<int>>mesh_neihaxi;//快速检索邻居节点哈希表
		//vector<vector<int>>mesh_bj;//边界三角形网格
		//vector<Point>mesh_gradient;//单元的梯度
		map<uint64_t, vector<int>>sharedface;
		map<uint64_t, double>sharedfai;
		vector<Point>Grad;
		vector<double>diff_k;

		vector<vector<gradcon>>Gradelements;
		vector<gradvalue>grad;

		map<uint64_t, int>fvm_bj;
	};
}

这里我们新增了几个变量用来对内部面和外部面进行一些判断，其中内部面即是遍历所有单元面出现两次的面，边界面即是只出现一次。
为了唯一的给定每个面一个ID，考虑三角面由三个点构成A,B,C，其中A,B,C点的序号为id_A,id_B,id_c, 有：
$$1<i{d}_A<numnodes$$
$$1<i{d}_B<numnodes$$
$$1<i{d}_C<numnodes$$
对于每个面的id，通过将面的三个点的节点序号从小到大顺序排列，并将其映射到立方体结构中，有唯一的id
$$(i{d}_A,i{d}_B,i{d}_C)=sort(i{d}_A,i{d}_B,i{d}_C);$$
$$id=i{d}_A+i{d}_B\ast numnodes+i{d}_C\ast numnodes\ast numnodes$$
我们使用shareface存储所有面的信息，通过面的id检索，对应位置shareface[id]存储了该面周围的两个或一个单元的编号。
同样的，sharefai存储了前一次获得的所有面的φ值，以便求全场梯度迭代使用。
题主开发的FVM代码
FVM-code

理论验证

与CVFEM类似，题主对同样的理论案例进行了验证，理论解析解可以在上篇中找到，这里只放一个验证的对比结果：
横坐标是所有控制体积质心到原点距离，纵坐标是φ值；

考虑对流-扩散无源变扩散系数的情况下，进口φ = 1.0，出口φ = 0.0；

困惑

题主在计算过程中，首先不考虑非正交影响拼凑出方程组后计算得到初解，然后根据初解计算交叉扩散项，并考虑了梯度的非正交影响，进行两次迭代计算全局梯度，将交叉扩散项放入源项拼装出新的方程组重新计算解，此时题主迭代一次得到的结果吻合度很好(上图对比)，多次迭代反而出现结果震荡越来越大，无法收敛，有没有知道原因的大佬解释下。