改進最小二乘法的非結構網格梯度重構算法

肖藝，明平劍

(哈爾濱工程大學 動力與能源工程學院，黑龍江 哈爾濱 150001)

非結構網格有限體積法的梯度重構是計算流 體力學(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 格心型有限體積法

格心型有限體積法的變量均儲存在單元中心，考慮穩態熱傳導控制方程：

·(kT)=0

(1)

式中：k為導熱系數；T為溫度。對于格心型有限體積法，利用高斯公式在每個單元內積分為：

(2)

對于非結構化網格：

(3)

1.2 最小二乘法

假定相連單元中心之間的變量線性分布為:

φF=φC+φCr·(rF-rC)

(4)

(5)

方程(5)為超定方程，可通過解函數L的最小值求解：

(6)

方程組(5)還可通過正規化或格萊姆-施密特正交化[5]求解，文獻[5]中指出，正規化法重構大長寬比的網格單元的梯度時會出現矩陣病態問題，進而導致梯度重構誤差較大。本文中的最小二乘問題均通過求解式(6)最小值求解。

計算模板單元的選取是最小二乘法的核心問題，以圖1中心單元為例，梯度重構的計算模板單元一般分為共面型模板(圖1中的深色單元)和共點型模板(圖1中的所有單元)。當網格質量較差時，引入更多單元的共點型模板能提高梯度重構精度，但也會導致計算量增加。對于二維的大長寬比網格，文獻[1]中分析了WLSQ(n)法在共面單元模板和共點單元模板的理論截斷誤差，指出當計算模板單元最大距離比L=max{rCj}/min{rCj}與網格長寬比同量級時(見圖2)，WLSQ(n)能保持梯度重構精度；而當計算模板單元距離比的最大值遠小于網格長單元比時，WLSQ(n)梯度重構精度下降。因此，計算模板單元的選取直接影響最小二乘法梯度重構的精度。

圖1 最小二乘梯度重構單元模板Fig.1 Least squares gradient reconstruction stencil

圖2 插值模板單元最大距離比Fig.2 Maximum distance ratio of interpolation stencil

1.3 格林-高斯法

格林-高斯法是應用最廣泛的空間梯度重構方法，由散度定理，單元中心變量空間梯度重構為：

(7)

式中：φfK為單元的第k個面的面心變量值；ΩC為單元的體積。

2 改進的加權最小二乘梯度重構

2.1 高斯加權最小二乘法

距離反比權重是加權的最小二乘法最常用的一種，但其與單元中心之間的方位無關。為了單元中心方位的影響，Shima[2]提出了高斯權重最小二乘法(WLSQ(G))，考慮格林-高斯法，式(6)可進一步改寫為:

(8)

式中：αk為線性插值系數。考慮加權最小二乘法，求解式(6)最小值得：

(9)

式中：rck=rk-rc；Δrck和lck分別為單元C和相連單元k之間的距離和單位距離矢量。

對比式(8)和(9)，在網格近似正交的情況下，可以認為nfk≈lck，僅考察單元k的權重，因此忽略常數系數1/Vo和M-1，可得:

ωk=αkSfk/Δrck=αWkSfk/Δrck

為在一維情況下達到二階精度，取αWk為：

(10)

即高斯加權最小二乘法的權函數為：

(11)

式中：rc為單元C的格心坐標向量；rfk、nfk分為單元C的第k個面的面心坐標向量、外法線單位面積矢量；rk為單元C的第k個相連單元的格心坐標向量。WLSQ(G)法的最終形式為:

(12)

式中gi只與網格幾何關系有關。圖3為格林-高斯權重示意圖。

圖3 格林-高斯權重示意Fig.3 Schematic of Green-Gauss weight

2.2 高斯加權最小二乘法理論截斷誤差分析

2.2.1 大長寬比彎曲四邊形網格

考慮圖4所示的大寬比彎曲四邊形網格，網格在坐標(r,θ)的空間步長分別為hr和hθ,網格長寬為：

A=Rhq/hr>>1

(13)

網格的變形程度為：

(14)

圖4中WLSQ(G)法重構單元O梯度的計算模板單元包括的N、S、W、E單元，以O為原點建立直角坐標系(e,n)，各單元中心在坐標系(r,θ)和(n,e)的坐標如表1所示。

圖4 大長寬比彎曲四邊形網格示意Fig.4 Schematic of large aspect ratio curved quadrilateral mesh

表1 四邊形網格模板單元中心在不同坐標系中的坐標值Table 1 Coordinate values of cell center of the stencil of quadrilateral mesh in different coordinate systems

在以O為原點的直角坐標系中，一般函數f(n,e)可由O、N、S、W、E共5個點通過線性最小二乘擬合。線性擬合函數fr(e,n)為：

fr(e,n)=fo+eeφ+nnφ

(15)

其中fo=f(0,0)，eφ和nφ通過求解式(12)的得到。由式(11)計算各點的權重為：

ωN=ωS=2Rtan(hθ/2)/hr≈Rhθ/hr=A

(16)

ωW=ωE=hr/[2Rsin(hθ/2)]≈hr/Rhθ=1/A

(17)

解式(12)得：

(18)

故WLSQ(G)法對圖2所示網格內部單元的梯度重構誤差為Ο(hθ)，為一階收斂精度。

2.2.2 大長寬比彎曲三角形網格

將圖4中的大寬比彎曲四邊形網格劃分為三角網格，如圖5所示。圖5中WLSQ(G)法重構單元O梯度的計算模板單元包括的1、2、3單元，以0為原點建立直角坐標系(e,n)，各單元中心在坐標系(r,θ)和(n,e)的坐標如表2所示。

圖5 大長寬比彎曲的三角形網格示意Fig.5 Schematic of large aspect ratio curved triangular mesh diagram

表2 三角形網格模板單元中心在不同坐標系中的坐標值Table 2 Coordinate values of cell center of the stencil of triangular mesh in different coordinate systems

由式(11)計算各點權重為：

ω1=ω2=3，ω3=3hr/(2Rhθ).

注意到Rhθ<

(19)

故WLSQ(G)法對圖5所示網格內部單元的梯度重構誤差為ε=Ο(1)，為零階收斂精度。

2.2.3 矩形計算區域內部單元

考慮圖6所示的大長寬比三角形網格，網格x和y方向的空間步長分別為hx和hy，網格長寬比A=hx/hy>>1。圖4中，WLSQ(G)法重構單元O梯度的計算模板單元為1、2、3單元，以O為原點建立直角坐標系(x,y)，各單元中心在坐標系(n,e)的坐標如表3所示。

由式(11)計算各點權重為：

式(12)得：

(20)

圖6 矩形域大長寬比三角形網格示意Fig.6 Schematic of large aspect ratio triangle mesh of rectangle domain

表3 模板單元中心在不同坐標系中的坐標值Table 3 Coordinate values of cell center of the stencil in different coordinate systems

故WLSQ(G)法對圖6所示網格內部單元的梯度重構誤差為Ο(hx)，為一階收斂精度。

2.3 改進邊界單元計算模板

對圖7所示的大長寬比彎曲三角形網格邊界單元O，WLSQ(G)法的計算模板單元由單元1、單元3、邊界2組成。由2.2.2中對三角形網格分析可知，共面型計算模板L=Ο(1)<

圖7 邊界單元梯度重構模板Fig.7 Gradient reconstruction stencil of boundary element

3 數值算例與分析

考慮網格的彎曲、大長寬比，分別對WLSQ(G)、VWLSQ(1)、WLSQ(G)、G-G，LSQ法在圖8所示的圓形域網格進行測試。采用了四邊形網格和三角形網格，其中三角形網格是在四邊形網格的基礎上有序的按對角線剖分為2個三角形得到。為測試不同梯度重構方法在邊界單元的數值收斂精度，按文獻[2]中對圖8所示的網格逐漸細化(分為別粗、中等、細3種網格尺寸)，網格尺寸見表4。測試函數為f=r2=x2+y2，相對誤差ε=|1-φ/(2r)|。

圖8 圓環域網格示意(中等網格尺寸)Fig.8 Schematic of ring domain mesh(medium mesh size)

表4 網格尺寸Table 4 Mesh size

求解單層壁圓筒穩態導熱問題，分別在粗、中等、細3種網格尺寸的三角形/四邊形網格測試。本章以中等規格三角形/四邊網格為例，分析不同梯度重構算法的計算結果。控制方程和邊界條件為：

(21)

式(20)的解析解為：

T*=T1+(T2-T1)ln(r/r1)/ln(r2/r1)

(22)

從圖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計算誤差較大。

圖9 圓環域四邊形網格梯度重構誤差分布Fig.9 Distribution of gradient reconstruction error of quadrilateral mesh in ring domain

從圖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<

圖10 圓環域四邊形網格邊界單元梯度重構誤差Fig.10 Gradient reconstruction error of boundary element of quadrilateral mesh in ring domain

圖11 大長寬比彎曲四邊形網格基于格心的最小二乘模板Fig.11 Cell-centered least squares stencil of large aspect ratio curved quadrilateral mesh

圖12 圓環域三角形網格梯度重構誤差分布 Fig.12 Distribution of gradient reconstruction Error of triangle mesh in ring domain

對圓環穩態導熱算例，從圖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)<

圖13 圓環域三角形網格邊界單元梯度重構誤差Fig.13 Gradient reconstruction error of boundary elements of triangle mesh in ring domain

圖14 大長寬比彎曲三角形網格基于格點的最小二乘模板Fig.14 Vertex-centered least squares stencil of large aspect ratio curved quadrilateral mesh

圖15 圓環穩態導熱四邊形網格誤差分布Fig.15 Error distribution of quadrilateral mesh of steady state heat conduction in ring domain

圖16 圓環穩態導熱誤差四邊形網格殘差收斂Fig.16 Convergence of residual of steady state heat conduction in ring domain

圖17 圓環穩態導熱三角形網格誤差分Fig.17 Error distribution of triangle mesh of steady state heat conduction inring domain

圖18 圓環穩態導熱三角形網格殘差收斂Fig.18 Convergence of residual of triangle mesh of steady state heat conduction in ring domain

表5 BWLSQ(G)在圓環穩態導熱三角形網格的計算結果Table 5 Calculation results of BWLSQ(G) of triangular mesh of steady-state heat conduction in the ring domain

表6 BWLSQ(G)在圓環穩態導熱四邊形網格的計算結果Table 6 Calculation results of BWLSQ(G) of quadrilateral mesh of steady-state heat conduction in the ring domain

4 結論

1)對于大長寬比彎曲三角形/四邊形網格，G-G法的計算結果要好于加權重的最小二乘法，且在邊界單元為一階收斂精度。

2)針對大長寬比彎曲三角形/四邊形網格的邊界單元，BWLSQ(G)法能有效減小邊界單元梯度重構誤差，達到一階精度；

3)針對四邊形網格，BWLSQ(G)計算效率與WLSQ(G)相同；針對三角形網格，BWLSQ(G)能有效提高WLSQ(G)的計算效率。

下一步工作需深入分析計算模板單元幾何關系和權重對最小二乘插值系數矩陣的影響，并且應用于更復雜的實際問題中測試其適用性。