改進(jìn)最小二乘法的非結(jié)構(gòu)網(wǎng)格梯度重構(gòu)算法

肖藝，明平劍

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

非結(jié)構(gòu)網(wǎng)格有限體積法的梯度重構(gòu)是計(jì)算流 體力學(xué)(computational fluid dynamics，CFD)的必要環(huán)節(jié)，其準(zhǔn)確性和健壯性至關(guān)重要，是CFD研究的重點(diǎn)之一。由于非結(jié)構(gòu)化網(wǎng)格高階精度格式還存在諸多問(wèn)題，目前應(yīng)用最為廣泛的方法仍然是具有二階空間精度有限體積法。對(duì)于二階空間精度的有限體積法，需要進(jìn)行分片恒定梯度重構(gòu)。對(duì)于高雷諾數(shù)的粘性流動(dòng)模擬，通常會(huì)遇到大長(zhǎng)寬比的彎曲網(wǎng)格，且邊界層單元類型比較復(fù)雜，復(fù)合層材料導(dǎo)熱中也會(huì)遇到此類網(wǎng)格。因此，非結(jié)構(gòu)網(wǎng)格有限體積法的梯度重構(gòu)仍然是有重要意義的研究課題。Diskin等[1]對(duì)非結(jié)構(gòu)網(wǎng)格有限體積法梯度重構(gòu)算法作了全面系統(tǒng)的歸類，并且分析了各種梯度重構(gòu)算法在二維大長(zhǎng)寬比網(wǎng)格上的理論截?cái)嗾`差和數(shù)值收斂精度。Shim等[2]提出了格林-高斯加權(quán)的最小二乘梯度重構(gòu)方法(weight least squares(green-gauss),WLSQ(G))，該方法的權(quán)函數(shù)的構(gòu)造綜合了距離反比加權(quán)最小二乘法(weight least squares，WLSQ)和格林-高斯法(green-gauss,G-G)，并且考慮了單元中心之間的方位影響。Mishrik等[3]以二維四邊形結(jié)構(gòu)化網(wǎng)格為例，從理論上分析了格心型最小二乘法的截?cái)嗾`差與網(wǎng)格幾何因素的相關(guān)性。張帆等[4]提出了基于格點(diǎn)的距離反比加權(quán)最小二乘梯度重構(gòu)法(vertex-based weight least squares,VWLSQ(1))，即基于節(jié)點(diǎn)加權(quán)最小二乘梯度重構(gòu)法。對(duì)于邊界單元梯度重構(gòu)，也有專門(mén)的研究[5-7]。

本文基于哈爾濱工程大學(xué)自主研發(fā)的大型通用輸運(yùn)方程數(shù)值分析軟件通用數(shù)值計(jì)算平臺(tái)，針對(duì)大長(zhǎng)寬比彎曲網(wǎng)格，改進(jìn) WLSQ(G)得到了 BWLSQ(G)。

1 非結(jié)構(gòu)網(wǎng)格有限體積法梯度重構(gòu)

1.1 格心型有限體積法

格心型有限體積法的變量均儲(chǔ)存在單元中心，考慮穩(wěn)態(tài)熱傳導(dǎo)控制方程：

·(kT)=0

(1)

式中：k為導(dǎo)熱系數(shù)；T為溫度。對(duì)于格心型有限體積法，利用高斯公式在每個(gè)單元內(nèi)積分為：

(2)

對(duì)于非結(jié)構(gòu)化網(wǎng)格：

(3)

1.2 最小二乘法

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

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

(4)

(5)

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

(6)

方程組(5)還可通過(guò)正規(guī)化或格萊姆-施密特正交化[5]求解，文獻(xiàn)[5]中指出，正規(guī)化法重構(gòu)大長(zhǎng)寬比的網(wǎng)格單元的梯度時(shí)會(huì)出現(xiàn)矩陣病態(tài)問(wèn)題，進(jìn)而導(dǎo)致梯度重構(gòu)誤差較大。本文中的最小二乘問(wèn)題均通過(guò)求解式(6)最小值求解。

計(jì)算模板單元的選取是最小二乘法的核心問(wèn)題，以圖1中心單元為例，梯度重構(gòu)的計(jì)算模板單元一般分為共面型模板(圖1中的深色單元)和共點(diǎn)型模板(圖1中的所有單元)。當(dāng)網(wǎng)格質(zhì)量較差時(shí)，引入更多單元的共點(diǎn)型模板能提高梯度重構(gòu)精度，但也會(huì)導(dǎo)致計(jì)算量增加。對(duì)于二維的大長(zhǎng)寬比網(wǎng)格，文獻(xiàn)[1]中分析了WLSQ(n)法在共面單元模板和共點(diǎn)單元模板的理論截?cái)嗾`差，指出當(dāng)計(jì)算模板單元最大距離比L=max{rCj}/min{rCj}與網(wǎng)格長(zhǎng)寬比同量級(jí)時(shí)(見(jiàn)圖2)，WLSQ(n)能保持梯度重構(gòu)精度；而當(dāng)計(jì)算模板單元距離比的最大值遠(yuǎn)小于網(wǎng)格長(zhǎng)單元比時(shí)，WLSQ(n)梯度重構(gòu)精度下降。因此，計(jì)算模板單元的選取直接影響最小二乘法梯度重構(gòu)的精度。

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

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

1.3 格林-高斯法

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

(7)

式中：φf(shuō)K為單元的第k個(gè)面的面心變量值；ΩC為單元的體積。

2 改進(jìn)的加權(quán)最小二乘梯度重構(gòu)

2.1 高斯加權(quán)最小二乘法

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

(8)

式中：αk為線性插值系數(shù)?？紤]加權(quán)最小二乘法，求解式(6)最小值得：

(9)

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

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

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

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

(10)

即高斯加權(quán)最小二乘法的權(quán)函數(shù)為：

(11)

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

(12)

式中g(shù)i只與網(wǎng)格幾何關(guān)系有關(guān)。圖3為格林-高斯權(quán)重示意圖。

圖3 格林-高斯權(quán)重示意Fig.3 Schematic of Green-Gauss weight

2.2 高斯加權(quán)最小二乘法理論截?cái)嗾`差分析

2.2.1 大長(zhǎng)寬比彎曲四邊形網(wǎng)格

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

A=Rhq/hr>>1

(13)

網(wǎng)格的變形程度為：

(14)

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

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

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

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

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

(15)

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

ω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)法對(duì)圖2所示網(wǎng)格內(nèi)部單元的梯度重構(gòu)誤差為Ο(hθ)，為一階收斂精度。

2.2.2 大長(zhǎng)寬比彎曲三角形網(wǎng)格

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

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

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

由式(11)計(jì)算各點(diǎn)權(quán)重為：

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

注意到Rhθ<

(19)

故WLSQ(G)法對(duì)圖5所示網(wǎng)格內(nèi)部單元的梯度重構(gòu)誤差為ε=Ο(1)，為零階收斂精度。

2.2.3 矩形計(jì)算區(qū)域內(nèi)部單元

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

由式(11)計(jì)算各點(diǎn)權(quán)重為：

式(12)得：

(20)

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

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

故WLSQ(G)法對(duì)圖6所示網(wǎng)格內(nèi)部單元的梯度重構(gòu)誤差為Ο(hx)，為一階收斂精度。

2.3 改進(jìn)邊界單元計(jì)算模板

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

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

3 數(shù)值算例與分析

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

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

表4 網(wǎng)格尺寸Table 4 Mesh size

求解單層壁圓筒穩(wěn)態(tài)導(dǎo)熱問(wèn)題，分別在粗、中等、細(xì)3種網(wǎng)格尺寸的三角形/四邊形網(wǎng)格測(cè)試。本章以中等規(guī)格三角形/四邊網(wǎng)格為例，分析不同梯度重構(gòu)算法的計(jì)算結(jié)果?？刂品匠毯瓦吔鐥l件為：

(21)

式(20)的解析解為：

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

(22)

從圖9(a)可看到，對(duì)于圓環(huán)形計(jì)算域四邊形網(wǎng)格算例中的內(nèi)部單元，G-G法和VWLSQ(1)法的計(jì)算結(jié)果相近，且比WLSQ(G)法誤差整體小1個(gè)量級(jí)。G-G和VWLSQ(1)法的誤差沿徑向增大，而WLSQ(G)的誤差沿徑向減小。從圖9(b)可看到，LSQ在1.0≤r≤1.1內(nèi)的誤差最大達(dá)到0.99。從圖10可看到，對(duì)于邊界單元，BWLSQ(G)、G-G、VWLSQ(1)的誤差相近，均為一階收斂精度；WLSQ(G)為二階收斂精度，但其誤差比BWLSQ(G)、G-G、VWLSQ(1)高2個(gè)量級(jí)；LSQ為零階收斂精度。分析計(jì)算模板單元的幾何關(guān)系，基于格心的最小二乘法的計(jì)算模板如圖11所示，其格心梯度由其計(jì)算模板單元內(nèi)插得到，且計(jì)算模板的最大距離比L=Ο(A)，進(jìn)而導(dǎo)致LSQ計(jì)算誤差較大。

圖9 圓環(huán)域四邊形網(wǎng)格梯度重構(gòu)誤差分布Fig.9 Distribution of gradient reconstruction error of quadrilateral mesh in ring domain

從圖12中可以看到，對(duì)于圓環(huán)域三角形網(wǎng)格內(nèi)部單元，G-G法誤差隨徑向增大，在1.0～1.1內(nèi)，其誤差為0.01～0.034，但其誤差在r=1.1附近最大；BWLSQ(G)和WLSQ(G)計(jì)算結(jié)果相同，誤差隨徑向減小；VWLSQ(1)法和LSQ法誤差在r=1附近均大于0.90，且均隨徑向減??；在r=1處，WLSQ(G)和BWLSQ(G)最大誤差為0.2，符合2.2.2節(jié)中理論截?cái)嗾`差分析得到的ε=Ο(1)，主要由徑向梯度誤差導(dǎo)致。從圖13中可知，對(duì)于邊界單元，BWLSQ(G) 為一階精度，其誤差最小，且比WLSQ(G)誤差小2個(gè)量級(jí)；G-G法誤差次之，為一階精度；WLSQ(G)為一階精度，誤差較大；LSQ和VWLSQ(1)均為零階精度，且誤差較大。分析計(jì)算模板單元的幾何關(guān)系，基于格點(diǎn)的距離反比權(quán)重最小二乘法(VWLSQ(1))的計(jì)算模板單元如圖14所示，當(dāng)網(wǎng)格長(zhǎng)寬比較大時(shí)，其格點(diǎn)梯度由外插值得到，且計(jì)算模板的最大距離比L≈2<

圖10 圓環(huán)域四邊形網(wǎng)格邊界單元梯度重構(gòu)誤差Fig.10 Gradient reconstruction error of boundary element of quadrilateral mesh in ring domain

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

圖12 圓環(huán)域三角形網(wǎng)格梯度重構(gòu)誤差分布 Fig.12 Distribution of gradient reconstruction Error of triangle mesh in ring domain

對(duì)圓環(huán)穩(wěn)態(tài)導(dǎo)熱算例，從圖15和圖16可知在四邊形網(wǎng)格上，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個(gè)量級(jí)。從圖18可以看出，LSQ計(jì)算發(fā)散；G-G和VWLSQ(1)殘差衰減最快，BWLSQ(G)次之；當(dāng)殘差收斂時(shí)，BWLSQ(G)的所用迭代步數(shù)為WLSQ(G)的1/3，故針對(duì)三角形網(wǎng)格，BWLSQ(G)能有效的提高WLSQ(G)的計(jì)算效率。分析可知，對(duì)于四邊形網(wǎng)格，BWLSQ(G)和WLSQ(G)的邊界待求單元的計(jì)算模板均滿足L=Ο(A)，因此在四邊形網(wǎng)格上的計(jì)算收斂性效率相同。對(duì)于三角形網(wǎng)格，WLSQ(G)的邊界待求單元的計(jì)算模板為L(zhǎng)=Ο(1)<

圖13 圓環(huán)域三角形網(wǎng)格邊界單元梯度重構(gòu)誤差Fig.13 Gradient reconstruction error of boundary elements of triangle mesh in ring domain

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

圖15 圓環(huán)穩(wěn)態(tài)導(dǎo)熱四邊形網(wǎng)格誤差分布Fig.15 Error distribution of quadrilateral mesh of steady state heat conduction in ring domain

圖16 圓環(huán)穩(wěn)態(tài)導(dǎo)熱誤差四邊形網(wǎng)格殘差收斂Fig.16 Convergence of residual of steady state heat conduction in ring domain

圖17 圓環(huán)穩(wěn)態(tài)導(dǎo)熱三角形網(wǎng)格誤差分Fig.17 Error distribution of triangle mesh of steady state heat conduction inring domain

圖18 圓環(huán)穩(wěn)態(tài)導(dǎo)熱三角形網(wǎng)格殘差收斂Fig.18 Convergence of residual of triangle mesh of steady state heat conduction in ring domain

表5 BWLSQ(G)在圓環(huán)穩(wěn)態(tài)導(dǎo)熱三角形網(wǎng)格的計(jì)算結(jié)果Table 5 Calculation results of BWLSQ(G) of triangular mesh of steady-state heat conduction in the ring domain

表6 BWLSQ(G)在圓環(huán)穩(wěn)態(tài)導(dǎo)熱四邊形網(wǎng)格的計(jì)算結(jié)果Table 6 Calculation results of BWLSQ(G) of quadrilateral mesh of steady-state heat conduction in the ring domain

4 結(jié)論

1)對(duì)于大長(zhǎng)寬比彎曲三角形/四邊形網(wǎng)格，G-G法的計(jì)算結(jié)果要好于加權(quán)重的最小二乘法，且在邊界單元為一階收斂精度。

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

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

下一步工作需深入分析計(jì)算模板單元幾何關(guān)系和權(quán)重對(duì)最小二乘插值系數(shù)矩陣的影響，并且應(yīng)用于更復(fù)雜的實(shí)際問(wèn)題中測(cè)試其適用性。