非结构网格有限体积法的梯度重构是计算流体力学(computational fluid dynamics，CFD)的必要环节，其准确性和健壮性至关重要，是CFD研究的重点之一。由于非结构化网格高阶精度格式还存在诸多问题，目前应用最为广泛的方法仍然是具有二阶空间精度有限体积法。对于二阶空间精度的有限体积法，需要进行分片恒定梯度重构。对于高雷诺数的粘性流动模拟，通常会遇到大长宽比的弯曲网格，且边界层单元类型比较复杂，复合层材料导热中也会遇到此类网格。因此，非结构网格有限体积法的梯度重构仍然是有重要意义的研究课题。Diskin等^[1]对非结构网格有限体积法梯度重构算法作了全面系统的归类，并且分析了各种梯度重构算法在二维大长宽比网格上的理论截断误差和数值收敛精度。Shim等^[2]提出了格林-高斯加权的最小二乘梯度重构方法(weight least squares(green-gauss), WLSQ(G))，该方法的权函数的构造综合了距离反比加权最小二乘法(weight least squares，WLSQ)和格林-高斯法(green-gauss, G-G)，并且考虑了单元中心之间的方位影响。Mishrik等^[3]以二维四边形结构化网格为例，从理论上分析了格心型最小二乘法的截断误差与网格几何因素的相关性。张帆等^[4]提出了基于格点的距离反比加权最小二乘梯度重构法(vertex-based weight least squares, VWLSQ(1))，即基于节点加权最小二乘梯度重构法。对于边界单元梯度重构，也有专门的研究^[5-7]。

本文基于哈尔滨工程大学自主研发的大型通用输运方程数值分析软件通用数值计算平台，针对大长宽比弯曲网格，改进WLSQ(G)得到了BWLSQ(G)。

1 非结构网格有限体积法梯度重构 1.1 格心型有限体积法

格心型有限体积法的变量均储存在单元中心，考虑稳态热传导控制方程：

$ \nabla {\rm{ }}\cdot{\rm{ }}(k\nabla T) = 0 $

(1)

式中：k为导热系数；T为温度。对于格心型有限体积法，利用高斯公式在每个单元内积分为：

$ \begin{array}{l} \int_\varOmega {\nabla {\rm{ }}\cdot{\rm{ }}(k\nabla T){\rm{d}}\mathit{\boldsymbol{V}}} = \oint\limits_\varOmega {k(\nabla T)\cdot\mathit{\boldsymbol{n}}{\rm{d}}S} {\rm{ }} = \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \oint_s {k\frac{{\partial T}}{{\partial n}}{\rm{d}}S} {\rm{ }} = \sum\limits_{f = 1}^{nf} {{k_f}{{\left( {\frac{{\partial T}}{{\partial n}}} \right)}_f}{S_f}} \end{array} $

(2)

式中：S_f为面积矢量；n为单位面积矢量；V为体积；$ \frac{{\partial T}}{{\partial n}}$为T沿n方向温度梯度。

对于非结构化网格：

$ {\left( {\frac{{\partial T}}{{\partial n}}} \right)_f} = {\left( {\nabla T} \right)_f}\cdot{\kern 1pt} {\kern 1pt} e + {\left( {\nabla T} \right)_f}\cdot{\kern 1pt} {\kern 1pt} \tau $

(3)

式中：$ {\left( {\nabla T} \right)_f} \cdot e$离散为两侧单元变量的线性关系；$ {\left( {\nabla T} \right)_f}\cdot\tau $则显式处理；$ {\left( {\nabla T} \right)_f}$由界面两侧单元梯度线性插值得到。但由于没有统一的笛卡尔坐标，差分方程不再适用于求解离散后的空间梯度项。为了克服这一缺点，最小二乘插值(least squares, LSQ)^[8-9]和格林-高斯法^[10-11]被用来重构有限体积法非结构化网格的空间梯度。Wang等^[12]研究了非结构网格最小二乘法和格林-高斯法的稳定性和数值收敛精度。

1.2 最小二乘法

假定相连单元中心之间的变量线性分布为:

$ {\phi _F} = {\phi _C} + {\rm{ }}\nabla {\phi _C}^r{\kern 1pt} {\kern 1pt} \cdot{\kern 1pt} {\kern 1pt} {\kern 1pt} ({\mathit{\boldsymbol{r}}_F} - {\mathit{\boldsymbol{r}}_C}) $

(4)

式中：ϕ_C和ϕ_F分别为带求梯度单元和与其相连单元的变量值；r_F和r_C分别为单元中心向量，其中$ \nabla {\phi _C}^r = \left( {{{\left( {\frac{{\partial \phi }}{{\partial x}}} \right)}_{{_C}}}\;{{\left( {\frac{{\partial \phi }}{{\partial y}}} \right)}_{{_C}}}} \right)$为拟合梯度。对所有相连单元列出式(3)，可得到方程组：

$ {\left[ {\begin{array}{*{20}{c}} {\Delta {x_{C{F_1}}}}&{\Delta {y_{C{F_1}}}}\\ \vdots & \vdots \\ {\Delta {x_{C{F_N}}}}&{\Delta {y_{C{F_N}}}} \end{array}} \right]_{N \times 2}}{\left[ {\begin{array}{*{20}{c}} {{{\left( {\frac{{\partial \phi }}{{\partial x}}} \right)}_C}}\\ {{{\left( {\frac{{\partial \phi }}{{\partial y}}} \right)}_C}} \end{array}} \right]_{2 \times 1}} = {\left[ {\begin{array}{*{20}{c}} \begin{array}{c} {\phi _{{F_1}}} - {\phi _C}\\ \vdots \end{array}\\ {{\phi _{{F_N}}} - {\phi _C}} \end{array}} \right]_{N \times 1}} $

(5)

方程(5)为超定方程，可通过解函数L的最小值求解：

$ L{\rm{ }} = \sum\limits_{k = 1}^{N(C)} {{{[{\omega _k}({\phi _{{F_k}}} - {\rm{ }}{\phi _C} - \nabla {\phi _C}\cdot{\rm{ }}\Delta {r_{C{F_k}}})]}^2}} $

(6)

式中ω_k为权重，ω_k=1即为LSQ法；ω_k=1/Δr_CFkⁿn=1, 2, 3, 即为加权最小二乘法WLSQ(n)。

方程组(5)还可通过正规化或格莱姆-施密特正交化^[5]求解，文献[5]中指出，正规化法重构大长宽比的网格单元的梯度时会出现矩阵病态问题，进而导致梯度重构误差较大。本文中的最小二乘问题均通过求解式(6)最小值求解。

计算模板单元的选取是最小二乘法的核心问题，以图 1中心单元为例，梯度重构的计算模板单元一般分为共面型模板(图 1中的深色单元)和共点型模板(图 1中的所有单元)。当网格质量较差时，引入更多单元的共点型模板能提高梯度重构精度，但也会导致计算量增加。对于二维的大长宽比网格，文献[1]中分析了WLSQ(n)法在共面单元模板和共点单元模板的理论截断误差，指出当计算模板单元最大距离比L=max{r_Cj}/min{r_Cj}与网格长宽比同量级时(见图 2)，WLSQ(n)能保持梯度重构精度；而当计算模板单元距离比的最大值远小于网格长单元比时，WLSQ(n)梯度重构精度下降。因此，计算模板单元的选取直接影响最小二乘法梯度重构的精度。

1.3 格林-高斯法

格林-高斯法是应用最广泛的空间梯度重构方法，由散度定理，单元中心变量空间梯度重构为：

$ \nabla {\phi _C} = \frac{1}{{{\varOmega _C}}}\sum\limits_{k = 1}^{{N_f}} {{\phi _{{f_k}}}{S_{{f_k}}}} $

(7)

式中：ϕ_fK为单元的第k个面的面心变量值；Ω_C为单元的体积。

2 改进的加权最小二乘梯度重构 2.1 高斯加权最小二乘法

距离反比权重是加权的最小二乘法最常用的一种，但其与单元中心之间的方位无关。为了单元中心方位的影响，Shima^[2]提出了高斯权重最小二乘法(WLSQ(G))，考虑格林-高斯法，式(6)可进一步改写为:

$ \nabla {\phi _C} = \frac{1}{{{\varOmega _C}}}\sum\limits_{k = 1}^{{N_f}} {{\phi _{{f_K}}}{\mathit{\boldsymbol{n}}_{{f_K}}}{S_{{f_K}}}} = \frac{1}{{{\varOmega _C}}}\sum\limits_{k = 1}^{{N_f}} {{\alpha _k}} {\mathit{\boldsymbol{n}}_{{f_k}}}{S_{{f_k}}}\Delta {\phi _k} $

(8)

式中：α_k为线性插值系数。考虑加权最小二乘法，求解式(6)最小值得：

$ \begin{array}{l} \nabla {\phi _C} = {\mathit{\boldsymbol{M}}^{ - 1}}\sum\limits_{k = 1}^{N(C)} {{\omega _k}{r_{ck}}\Delta {\phi _k}} = {\mathit{\boldsymbol{M}}^{ - 1}}\sum\limits_{k = 1}^{N(C)} {{\omega _k}\Delta {r_{ck}}{\mathit{\boldsymbol{l}}_{ck}}\Delta {\phi _k}} \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \mathit{\boldsymbol{M}} = \left[ {\begin{array}{*{20}{l}} {{I_{xx}}}&{{I_{xy}}}\\ {{I_{yx}}}&{{I_{yy}}} \end{array}} \right];{I_{AB}} = \sum\limits_{k = 1}^{N(C)} {{\omega _k}} \Delta A\Delta B \end{array} $

(9)

式中：$ \mathit{\boldsymbol{r}}{_{ck}} = {\rm{ }}\mathit{\boldsymbol{r}}{_k} - {\rm{ }}\mathit{\boldsymbol{r}}{_c} ; \Delta {r_{ck}}$和l_ck分別为单元C和相连单元k之间的距离和单位距离矢量。

对比式(8)和(9)，在网格近似正交的情况下，可以认为n_fk≈ l_ck，仅考察单元k的权重，因此忽略常数系数1/V_o和M^－1，可得:

$ {\omega _k} = {\alpha _k}{S_{{f_k}}}/\Delta {r_{ck}} = {\alpha _{{W_k}}}{S_{{f_k}}}/\Delta {r_{ck}} $

为在一维情况下达到二阶精度，取α_Wk为：

$ {\alpha _{{W_k}}} = {\left[ {\frac{{2({\mathit{\boldsymbol{r}}_{{f_k}}} - {\mathit{\boldsymbol{r}}_c}) \cdot {\mathit{\boldsymbol{n}}_{{f_k}}}}}{{({\mathit{\boldsymbol{r}}_k} - {\mathit{\boldsymbol{r}}_c}) \cdot {\mathit{\boldsymbol{n}}_{{f_k}}}}}} \right]^2} $

(10)

即高斯加权最小二乘法的权函数为：

$ {\omega _k} = {\left[ {\frac{{2({\mathit{\boldsymbol{r}}_{{f_k}}} - {\mathit{\boldsymbol{r}}_c}) \cdot {\mathit{\boldsymbol{n}}_{{f_k}}}}}{{({\mathit{\boldsymbol{r}}_k} - {\mathit{\boldsymbol{r}}_c}) \cdot {\mathit{\boldsymbol{n}}_{{f_k}}}}}} \right]^2}\frac{{{S_{{f_k}}}}}{{\left| {{\mathit{\boldsymbol{r}}_k} - {\mathit{\boldsymbol{r}}_c}} \right|}} $

(11)

式中：r_c为单元C的格心坐标向量；r_fk、n_fk分为单元C的第k个面的面心坐标向量、外法线单位面积矢量；r_k为单元C的第k个相连单元的格心坐标向量。WLSQ(G)法的最终形式为:

$ \nabla \phi = {\mathit{\boldsymbol{M}}^{ - 1}}\sum\limits_{i = 1}^N {{\omega _i}} \Delta {\phi _i}\Delta {\mathit{\boldsymbol{r}}_i} = \sum\limits_{i = 1}^N \Delta {\phi _i}{\mathit{\boldsymbol{g}}_i} $

(12)

式中g_i只与网格几何关系有关。图 3为格林-高斯权重示意图。

2.2 高斯加权最小二乘法理论截断误差分析 2.2.1 大长宽比弯曲四边形网格

考虑图 4所示的大宽比弯曲四边形网格，网格在坐标(r, θ)的空间步长分别为h_r和h_θ, 网格长宽为：

$ A{\rm{ }} = {\rm{ }}R{h_q}/{\rm{ }}{h_r} \gg 1 $

(13)

网格的变形程度为：

$ \Gamma {\rm{ }} = {\rm{ }}R(1 - {\rm{cos}}\theta )/{\rm{ }}{h_r} \approx Rh_\theta ^2/(2{h_r}) = {\rm{ }}A{h_\theta }/2 $

(14)

图 4中WLSQ(G)法重构单元O梯度的计算模板单元包括的N、S、W、E单元，以O为原点建立直角坐标系(e, n)，各单元中心在坐标系(r, θ)和(n, e)的坐标如表 1所示。

在以O为原点的直角坐标系中，一般函数f(n, e)可由O、N、S、W、E共5个点通过线性最小二乘拟合。线性拟合函数f^r(e, n)为：

$ {f^r}(e,n){\rm{ }} = {f_o} + {\rm{ }}e{\nabla _e}\phi {\rm{ }} + {\rm{ }}n{\nabla _n}\phi $

(15)

其中f_o=f(0, 0)，$ {\nabla _e}\phi $和${\nabla _n}\phi $通过求解式(12)的得到。由式(11)计算各点的权重为：

$ {\omega _N} = {\rm{ }}{\omega _S} = {\rm{ }}2R{\rm{tan}}({h_\theta }/2)/{\rm{ }}{h_r} \approx R{h_\theta }/{\rm{ }}{h_r} = A $

(16)

$ {\omega _W} = {\rm{ }}{\omega _E} = {\rm{ }}{h_r}/[2R{\rm{sin}}({h_\theta }/2)]{\rm{ }} \approx {h_r}/{\rm{ }}R{h_\theta } = 1/{\rm{ }}A $

(17)

解式(12)得：

$ \left\{ \begin{array}{l} {\nabla _e}\phi = {f_e} + {\rm{ }}O({h_\theta })\\ {\nabla _n}\phi = {f_n} + {\rm{ }}O({h_\theta }) \end{array} \right. $

(18)

故WLSQ(G)法对图 2所示网格内部单元的梯度重构误差为O(h_θ)，为一阶收敛精度。

2.2.2 大长宽比弯曲三角形网格

将图 4中的大宽比弯曲四边形网格划分为三角网格，如图 5所示。图 5中WLSQ(G)法重构单元O梯度的计算模板单元包括的1、2、3单元，以0为原点建立直角坐标系(e, n)，各单元中心在坐标系(r, θ)和(n, e)的坐标如表 2所示。

由式(11)计算各点权重为：

$ {\omega _1} = {\omega _2} = 3,{\omega _3} = 3{h_r}/(2R{h_\theta }). $

注意到Rh_θ < < h_r, 且R=1，式(12)得：

$ \left\{ \begin{array}{l} {\nabla _e}\phi = {f_e} + {\rm{ }}O({h_\theta })\\ {\nabla _n}\phi = {f_n} + {\rm{ }}O(1) \end{array} \right. $

(19)

故WLSQ(G)法对图 5所示网格内部单元的梯度重构误差为ε=Ο(1)，为零阶收敛精度。

2.2.3 矩形计算区域内部单元

考虑图 6所示的大长宽比三角形网格，网格x和y方向的空间步长分别为h_x和h_y，网格长宽比A=h_x/h_y> >1。图 4中，WLSQ(G)法重构单元O梯度的计算模板单元为1、2、3单元，以O为原点建立直角坐标系(x, y)，各单元中心在坐标系(n, e)的坐标如表 3所示。

由式(11)计算各点权重为：

$ {\omega _1} = 3,{\omega _2} = 3{h_x}/\sqrt {{h_x}^2 + 4{h_{{y^2}}}} , $

$ {\omega _3} = 3{h_x}/\sqrt {4{h_x}^2 + {h_{{y^2}}}} $

式(12)得：

$ \left\{ \begin{array}{l} {\nabla _x}\phi = {f_x} + {\rm{ }}O({h_x})\\ {\nabla _y}\phi = {f_y} + {\rm{ }}O({h_x}) \end{array} \right. $

(20)

故WLSQ(G)法对图 6所示网格内部单元的梯度重构误差为Ο(h_x)，为一阶收敛精度。

2.3 改进边界单元计算模板

对图 7所示的大长宽比弯曲三角形网格边界单元O，WLSQ(G)法的计算模板单元由单元1、单元3、边界2组成。由2.2.2中对三角形网格分析可知，共面型计算模板L=Ο(1) < < A，若变量值只沿径向变化，且沿切向几乎无变化，则由此计算模板重构梯度在径向存在较大误差。为了扩大计算模板单元，本文将边界单元梯度重构的权重取为ω=1/r_iⁿ。对于格心型WLSQ(N)法，只有当计算模板的各单元到中心单元距离的比值的最大值Li=Ο(A)时，WLSQ(N)才能有效减小梯度重构误差^[1]。因此，当引入图 5中的4号单元时，可满足计算模板Li=Ο(A)。为减小WLSQ(G)边界单元梯度重构误差，本文中将边界单元的梯度重构计算单元模板扩大为与边界单元的所有共节点单元(图 7中“+”号单元)，得到BWLSQ(G)法(Bound-WLSQ(G))。

3 数值算例与分析

考虑网格的弯曲、大长宽比，分别对WLSQ(G)、VWLSQ(1)、WLSQ(G)、G-G，LSQ法在图 8所示的圆形域网格进行测试。采用了四边形网格和三角形网格，其中三角形网格是在四边形网格的基础上有序的按对角线剖分为2个三角形得到。为测试不同梯度重构方法在边界单元的数值收敛精度，按文献[2]中对图 8所示的网格逐渐细化(分为别粗、中等、细3种网格尺寸)，网格尺寸见表 4。测试函数为f=r²=x²+y²，相对误差$ \varepsilon = \left| {1 - \nabla \phi /\left( {2r} \right)} \right|$。

求解单层壁圆筒稳态导热问题，分别在粗、中等、细3种网格尺寸的三角形/四边形网格测试。本章以中等规格三角形/四边网格为例，分析不同梯度重构算法的计算结果。控制方程和边界条件为：

$ \left\{ {\begin{array}{*{20}{l}} {\frac{{\rm{d}}}{{{\rm{d}}r}}\left( {r\frac{{{\rm{d}}T}}{{{\rm{d}}r}}} \right) = 0}\\ {r = {r_1},t = {T_1}}\\ {r = {r_2},t = {T_2}} \end{array}} \right. $

(21)

式(20)的解析解为：

$ {T^*} = {T_1} + ({T_2} - {T_1}){\rm{ln}}(r/{r_1})/{\rm{ln}}({r_2}/{r_1}) $

(22)

式中：r₁=1.0, r₂=5.598, T₁=50.0, T₂=10.0, 相对误差为$\varepsilon = \left| {1 - \widetilde T{\rm{ }}/{T^*}} \right|, \widetilde T{\rm{ }} $为数值解。

从图 9(a)可看到，对于圆环形计算域四边形网格算例中的内部单元，G-G法和VWLSQ(1)法的计算结果相近，且比WLSQ(G)法误差整体小1个量级。G-G和VWLSQ(1)法的误差沿径向增大，而WLSQ(G)的误差沿径向减小。从图 9(b)可看到，LSQ在1.0≤r≤1.1内的误差最大达到0.99。从图 10可看到，对于边界单元，BWLSQ(G)、G-G、VWLSQ(1)的误差相近，均为一阶收敛精度；WLSQ(G)为二阶收敛精度，但其误差比BWLSQ(G)、G-G、VWLSQ(1)高2个量级；LSQ为零阶收敛精度。分析计算模板单元的几何关系，基于格心的最小二乘法的计算模板如图 11所示，其格心梯度由其计算模板单元内插得到，且计算模板的最大距离比L=Ο(A)，进而导致LSQ计算误差较大。

从图 12中可以看到，对于圆环域三角形网格内部单元，G-G法误差随径向增大，在1.0~1.1内，其误差为0.01~0.034，但其误差在r=1.1附近最大；BWLSQ(G)和WLSQ(G)计算结果相同，误差随径向减小；VWLSQ(1)法和LSQ法误差在r=1附近均大于0.90，且均随径向减小；在r=1处，WLSQ(G)和BWLSQ(G)最大误差为0.2，符合2.2.2节中理论截断误差分析得到的ε=Ο(1)，主要由径向梯度误差导致。从图 13中可知，对于边界单元，BWLSQ(G)为一阶精度，其误差最小，且比WLSQ(G)误差小2个量级；G-G法误差次之，为一阶精度；WLSQ(G)为一阶精度，误差较大；LSQ和VWLSQ(1)均为零阶精度，且误差较大。分析计算模板单元的几何关系，基于格点的距离反比权重最小二乘法(VWLSQ(1))的计算模板单元如图 14所示，当网格长宽比较大时，其格点梯度由外插值得到，且计算模板的最大距离比L≈2 < < A进而导致基于格点的距离反比权重最小二乘法误差较大；当长宽比随径向减小时，其格点梯度为内插值得到，故其误差随径向减小。

对圆环稳态导热算例，从图 15和图 16可知在四边形网格上，VWLSQ、LSQ、GG、BWLSQ(G)、WLSQ(G)的误差和残差收敛速度相同。从图 17可以看出，VWLSQ(1)整体误差最大，BWLSQ(G)和WLSQ(G)次之，G-G误差最小；在r=1处，VWLSQ(1)的误差要比G-G、BWLSQ(G)和WLSQ(G)高1个量级。从图 18可以看出，LSQ计算发散；G-G和VWLSQ(1)残差衰减最快，BWLSQ(G)次之；当残差收敛时，BWLSQ(G)的所用迭代步数为WLSQ(G)的1/3，故针对三角形网格，BWLSQ(G)能有效的提高WLSQ(G)的计算效率。分析可知，对于四边形网格，BWLSQ(G)和WLSQ(G)的边界待求单元的计算模板均满足L=Ο(A)，因此在四边形网格上的计算收敛性效率相同。对于三角形网格，WLSQ(G)的边界待求单元的计算模板为L=Ο(1) < < A；而BWLSQ(G)的计算模板满足L=Ο(A)，因此，BWLSQ(G)能减小边界单元梯度重构的误差，进而计算收敛效率。如表 5所示，对圆环稳态导热，当BWLSQ(G)边界单元权重为ω=1/r_iⁿ=1/r_i⁴时，对于不同尺寸的三角形网格均能得到收敛解；而当n=1, 2, 3时，计算结果不稳定。如表 6所示，BWLSQ(G)在边界单元权重取不同值时均能收敛。增大距离倒数次数可以减小误差，但次数不能无限取大，其可能会造成最小二乘病态问题，进而导致误差。综合三角形网格和四边形网格测试结果，当边界单元权重取距离倒数四次方时稳定性较好。综合圆环导热测试结果表明，G-G、WLSQ(G)、BWLSQ(G)、VWLSQ(1)计算稳定性较好；未加权的LSQ法计算稳定性较差，并不适用于大长宽比的三角形网格。

4 结论

1) 对于大长宽比弯曲三角形/四边形网格，G-G法的计算结果要好于加权重的最小二乘法，且在边界单元为一阶收敛精度。

2) 针对大长宽比弯曲三角形/四边形网格的边界单元，BWLSQ(G)法能有效减小边界单元梯度重构误差，达到一阶精度；

3) 针对四边形网格，BWLSQ(G)计算效率与WLSQ(G)相同；针对三角形网格，BWLSQ(G)能有效提高WLSQ(G)的计算效率。

下一步工作需深入分析计算模板单元几何关系和权重对最小二乘插值系数矩阵的影响，并且应用于更复杂的实际问题中测试其适用性。