流場特性驅(qū)動的亞松弛因子自動獲取方法

摘要：為解決傳統(tǒng)亞松弛因子選取過度依賴計算經(jīng)驗并且沒有考慮流場本身特性的問題，提出了一種流場特性驅(qū)動的亞松弛因子自動獲取方法。針對壓力耦合方程組的半隱式方法（SIMPLE）系列算法，基于流場速度梯度給出了用于表征速度分布的表征量，在此基礎(chǔ)上利用非線性函數(shù)對該表征量進(jìn)行轉(zhuǎn)化，導(dǎo)出了速度亞松弛因子的計算表達(dá)式，并且在頂蓋驅(qū)動流、后臺階流和水箱流3個經(jīng)典算例下對比分析了固定亞松弛因子與自動獲取方法在一階迎風(fēng)格式、對流項二階迎風(fēng)插值（QUICK）格式和保證隱定性的二階插分格式下的計算表現(xiàn)。結(jié)果表明：所提自動獲取方法在各算例和離散格式中的計算時間均與固定亞松弛因子下的最短計算時間相當(dāng)，有些工況下甚至小于固定亞松弛因子的最短計算時間；當(dāng)后臺階流算例采用30×180和30×210的網(wǎng)格，QUICK格式下固定亞松弛因子為0.1、0.5、0.9時，計算均不收斂，改用所提自動獲取方法后計算收斂；采用自動獲取方法的計算殘差最終可以達(dá)到10^－14量級，與收斂情況下的固定亞松弛因子量級相同。研究為進(jìn)一步優(yōu)化亞松弛因子自動獲取方法提供了參考。

關(guān)鍵詞：SIMPLE算法；有限體積法；高精度格式；亞松弛因子

中圖分類號：TK124 文獻(xiàn)標(biāo)志碼：A

DOI：10.7652/xjtuxb202504004 文章編號：0253-987X（2025）04-0040-10

Flow Field Characteristics-Driven Automatic Under-Relaxation Factor Acquisition Method

KANG Jie， FAN Hua， WU Dongyin

（School of Energy and Power Engineering， Xi’an Jiaotong University， Xi’an 710049， China）

Abstract：To address the issue of traditional under-relaxation factor selection overly reliant on computational experience and overlooking inherent flow field characteristics， a method for automatically acquiring under-relaxation factors driven by flow field characteristics is proposed. For the semi-implicit method for pressure-linked equations （SIMPLE） series algorithm， a quantity representing the velocity distribution based on the velocity gradient of the flow field is provided. Utilizing a nonlinear function to transform this quantity， an expression for calculating the velocity under-relaxation factor is derived. Comparative analyses of fixed under-relaxation factors and the automatic acquisition method are conducted in three classic cases of lid-driven flow， backward-facing step flow， and tank flow under first-order upwind， quadratic upwind interpolation of convective kinematics （QUICK）， and stability-guaranteed second-order difference schemes. The results show that the proposed automatic acquisition method yields computation times comparable to the shortest computation times under fixed under-relaxation factors in all cases and discretization schemes， and in some scenarios， even shorter than the shortest computation time under fixed under-relaxation factors. In the case of the backward-facing step flow with grid sizes of 30×180 and 30×210 under the QUICK scheme， computations do not converge with fixed under-relaxation factors set at 0.1， 0.5， and 0.9， but converge when using the proposed automatic acquisition method. The residual values obtained using the automatic acquisition method can reach the order of 10^－14， equivalent to the order of fixed under-relaxation factors under convergence conditions. This research provides a reference for further optimizing automatic under-relaxation factor acquisition methods.

Keywords：SIMPLE algorithm; finite volume method; high-order scheme; under-relaxation factors

目前，在流動傳熱計算過程中，離散方程壓力與速度的求解方法可以分為聯(lián)立直接求解和分離式求解兩類。壓力耦合方程組的半隱式方法（SIMPLE）系列算法是分離式求解中的一種經(jīng)典算法，該算法及其衍生算法在工程實踐中得到了廣泛的應(yīng)用^［1-7］。SIMPLE算法在預(yù)估-校正的過程中引入了亞松弛因子，但是目前關(guān)于亞松弛因子的選取并沒有一個最優(yōu)的方法，在流動傳熱問題的計算中大多依賴以往的計算經(jīng)驗。

許多學(xué)者圍繞著如何選取最優(yōu)的亞松弛因子進(jìn)行了研究。在固定亞松弛因子方面，Patankar^［8］提出將速度和壓力亞松弛因子分別設(shè)為0.5和0.8，以確保大多數(shù)情況下的迭代收斂性。Barron等^［9］進(jìn)一步分析了亞松弛因子在湍流計算中的作用，并推薦將速度和壓力亞松弛因子分別設(shè)置在0.8～0.9和0.2的范圍內(nèi)。Marek等^［10］提出了一種基于殘差收斂與發(fā)散趨勢的混合松弛格式。Gopinath等^［11］提出了一種基于正交陣列的亞松弛因子選擇方法。

在自動獲取亞松弛因子的過程中，有許多基于殘差的方法。Chatwani等^［12］從動量方程殘差入手，提出了一種自動獲取壓力亞松弛因子的方法，但是這種方法在提高SIMPLE算法計算效率方面的效果不夠理想。范曉偉等^［13］基于離散方程殘差平方和最小的原理，建立了一種壓力亞松弛因子的自動生成方法，并且通過計算驗證了該方法的可靠性。在亞松弛因子自動獲取的研究中，模糊控制方法也被引入。Liu等^［14］采用模糊控制方法，為SIMPLER算法提出了一種自動獲取最佳速度亞松弛因子的策略，該方法在速度亞松弛因子較小時的加速效果顯著。王艷寧等^［15］將這一方法擴(kuò)展到三維流動問題的求解中，并且通過經(jīng)典算例的計算，證實了模糊控制方法的有效性。除此之外，Morii等^［16］基于速度修正值，通過亞松弛因子補(bǔ)償了在修正速度時所舍棄的相鄰節(jié)點速度修正值的影響，從而提出了SOAR算法。馬淼等^［17］單獨對壓力的亞松弛因子加上了自適應(yīng)調(diào)整，一定程度上節(jié)省了計算時間。

上述自動獲取方法大多著眼于數(shù)學(xué)上的推演，而忽略了流場本身的性質(zhì)。陶文銓團(tuán)隊^［18］提出可以將離散方程中周圍鄰點對應(yīng)的系數(shù)拆分為對流部分和擴(kuò)散部分，通過平衡對流部分和擴(kuò)散部分，從而推導(dǎo)得出亞松弛因子的計算表達(dá)式。

本文從流場本身的流動特性出發(fā)，提出了一種亞松弛因子的自動獲取方法，并且采用SIMPLE算法編寫求解器，通過對頂蓋驅(qū)動流、后臺階流和水箱流這3個經(jīng)典算例進(jìn)行計算，對比研究了固定亞松弛因子和自動獲取亞松弛因子在不同算例中的表現(xiàn)，驗證了該自動獲取方法的實用性。

1 SIMPLE算法與離散格式

一般在處理流動傳熱問題時，連續(xù)性方程、動量守恒方程和能量守恒方程可以統(tǒng)一寫成

（ρφ） t+ div （ρUφ）= div （Γφ grad φ）+Sφ（1）

式中：ρ為密度， kg/m 3；t為時間， s ；φ為廣義變量，在連續(xù)性方程中為常數(shù)1，在動量守恒方程中為速度，在能量守恒方程中為溫度；U為速度矢量；Γφ為廣義擴(kuò)散系數(shù)；Sφ為廣義源項。

有限體積法基于相鄰控制體之間的關(guān)聯(lián)，將控制方程進(jìn)行離散，最終可以得到控制方程的離散方程。如圖1所示，以控制體O為基準(zhǔn)，W和W′分別代表控制體O西面的第1個和第2個控制體體心，E和E′分別代表控制體O東面的第1個和第2個控制體體心，w和e分別代表控制體O的西界面和東界面。后文中，各變量采用大寫角標(biāo)W′、W、O、E、E′代表該變量為相應(yīng)控制體體心處變量，各變量采用小寫角標(biāo)w、e代表該變量為相應(yīng)控制體界面處變量。

在SIMPLE算法中，以二維穩(wěn)態(tài)流動中的x方向動量方程為例，引入亞松弛因子之后，對應(yīng)的離散方程變?yōu)?/p>

aOαuO=∑IaIuI+b+（pO-pE） Δ y+aOα（1－α）u ^old O

（2）

式中：a為離散方程的系數(shù)；α為速度u的亞松弛因子；u為水平方向的速度， m/s ；I為與控制體O相鄰方向的控制體體心，二維網(wǎng)格中I∈{W，E，N，S}；b為源項；p為壓力， Pa ； Δ y為控制體高度， m ；u ^old O為上一輪迭代的速度， m/s 。

離散方程中系數(shù) a 的計算，與所采用的離散格式息息相關(guān)，如果離散格式選擇不當(dāng)，也會影響計算的收斂性^［19］。在交錯網(wǎng)格^［20］的基礎(chǔ)上，可以通過不同的離散格式計算離散方程中的系數(shù)。

1.1 一階迎風(fēng)格式

在一階迎風(fēng)格式中，僅考慮控制體上游方向的物理量φ對控制體的影響。以速度控制體為例，東界面上的速度ue可以表示為

ue=uO，uegt;0uE，uelt;0（3）

1.2 對流項二階迎風(fēng)插值（QUICK）格式

QUICK格式是一種改進(jìn)的迎風(fēng)格式^［21］。在QUICK格式中，不僅考慮了控制體上游方向的物理量φ對控制體的影響，同時也考慮了控制體下游方向的物理量φ對控制體的影響。東界面上的速度ue可以表示為

ue=uO+uE2－18Q（4）

Q=uE－2uO+uW，uegt;0

uO－2uE+uE′，uelt;0

1.3 保證穩(wěn)定性的二階插分（SGSD）格式

在描述界面上的物理量時，由于中心差分格式和二階迎風(fēng)格式存在各自的優(yōu)缺點，所以有學(xué)者提出可以將這兩種格式結(jié)合起來，這就是SCSD格式^［22］。以交錯網(wǎng)格下的速度控制體為例，在SCSD格式中，東界面上的物理量可以表示為

φ ^SCSD e=βφ ^CD e+（1－β）φ ^SUD e（5）

式中：φ ^CD e為中心差分格式計算出的物理量東界面值；φ ^SUD e為二階迎風(fēng)格式計算出的物理量東界面值；β為加權(quán)系數(shù)，0≤β≤1。

在SCSD格式的基礎(chǔ)之上，SGSD格式^［23］提供了一種β的計算方法，可以很好地提高計算的穩(wěn)定性。 β 計算式如下

β=22+P Δ （6）

式中：P Δ 為Pe數(shù)。

2 自動獲取方法的建立

在固定亞松弛因子方法中，全流場采用固定的亞松弛因子進(jìn)行計算，固定值的選取無法通過統(tǒng)一的公式計算，只能通過試算并結(jié)合以往經(jīng)驗來確定。在這種傳統(tǒng)方法下，亞松弛因子沒有考慮流場的速度分布情況，也無法針對流場的不同區(qū)域采取因地制宜的亞松弛策略。

為了提高計算效率，需要建立一種可以因地制宜調(diào)整每個控制體亞松弛因子的方法。首先，需要尋找一個指標(biāo)來表征每個控制體所處位置的速度分布情況。以x方向動量方程為例，在引入亞松弛因子之后，將式（2）中的∑IaIuI展開，得到

∑IaIuI=aEuE+aWuW+aNuN+aSuS

（7）

式中：系數(shù)aE、aW、aN和aS根據(jù)具體的離散格式計算得到。在系數(shù)的計算的過程中，會涉及到Fe、Fw、Fn和Fs

Fe=ρue Δ y

Fw=ρuw Δ y

Fn=ρvn Δ x

Fs=ρvs Δ x（8）

式中：Fe、Fw、Fn和Fs代表了通過控制體界面流入（或流出）控制體的流量。由于ue、uw、vn和vs分別代表控制體界面上的速度，所以流場在x方向的速度梯度表征量δx可以用系數(shù)a的表達(dá)式進(jìn)行表征

δx=aE－aWaO（9）

同理，對于 y 方向，有

δy=aN－aSaO（10）

每個控制體速度梯度的表征量 δ 為

δ=aE－aWaO+aN－aSaO（11）

根據(jù)文獻(xiàn)［21］，離散方程的系數(shù)滿足

aO=aE+aW+aN+aS（12）

根據(jù)絕對值不等式，可以求得 δ 滿足

0≤δ≤aE+aW+aN+aSaO（13）

0≤δ≤1（14）

雖然 δ 表征了流動的速度梯度，但是無法直接當(dāng)作亞松弛因子進(jìn)行使用。這是因為根據(jù)式（11），在對流作用較小的區(qū)域，δ取值可能達(dá)到10^－3量級，如果將其直接作為亞松弛因子使用，反而會嚴(yán)重拖慢計算速度。為了避免此類極端值帶來的負(fù)面影響，可以考慮引入合適的數(shù)學(xué)表達(dá)式對 δ 的表達(dá)式進(jìn)行鈍化處理。本文推薦選取指數(shù)函數(shù)作為映射函數(shù)，建立了δ［0，1］fα auto ［0.1，0.9］的映射關(guān)系，以便控制亞松弛因子處于合適的范圍。指數(shù)函數(shù)的一般表達(dá)式為

f（δ）=A－B e ^Cδ（15）

式中：A、B和C均為常系數(shù)。

將 f（δ） 映射到［0.1，0.9］這一區(qū)間，有

α auto =0.1+0.8f（δ）－f（0）f（1）－f（0）（16）

α auto =0.1+0.81－ e ^Cδ1－ e C（17）

經(jīng)過化簡可以看到，化簡結(jié)果與系數(shù)A和B無關(guān)，僅與C有關(guān)。本文選取C=－1.2，此時式（17）變?yōu)?/p>

α auto =0.1+0.81－ e ^－1.2δ1－ e ^－1.2（18）

引入非線性函數(shù)之后，由于α auto 與δ是正相關(guān)的，所以α auto 同樣可以用來表征控制體所在位置的速度梯度。另一方面，引入非線性函數(shù)后，在δ的變化過程中，即便δ=0，也有α auto =0.1，這樣就避免了極端的δ導(dǎo)致的計算過慢問題。經(jīng)過鈍化之后，α auto 可以根據(jù)流場的速度梯度調(diào)整亞松弛因子的大小，因地制宜地選取與流場流速分布相適應(yīng)的亞松弛因子，以提高計算的速度。綜上所述可知，式（18）可以作為一種較好的速度亞松弛因子自動獲取方法。

3 算例介紹

為了檢驗該自動獲取方法在實際計算中的表現(xiàn)，本文選取了3個經(jīng)典算例，包括頂蓋驅(qū)動流、后臺階流和水箱流，通過與傳統(tǒng)的固定亞松弛因子方法進(jìn)行對比，來驗證自動獲取方法的實用性。

3.1 頂蓋驅(qū)動流

頂蓋驅(qū)動流是驗證流體動力學(xué)問題計算方法的一個基準(zhǔn)案例。如圖2所示，該模型由一個二維方腔組成，方腔的頂蓋具有切向速度，其他3個壁面速度為0。此運動引起了一種流動現(xiàn)象，其特征是在空腔中心有一個大渦，在角落有一些小渦，而雷諾數(shù)的大小影響著流動中渦的大小和數(shù)量。

本文取頂蓋驅(qū)動流的尺寸為1m×1m。頂蓋的水平速度u top =1 m/s ，豎直方向速度v top =0，其余兩條邊界為無滑移壁面，計算域被劃分為40×40和80×80兩套網(wǎng)格。該模型中，雷諾數(shù)定義如下

Re1=u top L top ν（19）

式中：L top 為頂蓋的長度， m ；ν為運動黏性系數(shù)， m2/s 。

3.2 后臺階流

如圖3所示，在后臺階流中，流體先流過一個狹窄的入口流道，接著進(jìn)入橫截面突然變大的流道。流體在流經(jīng)臺階處的時候，會形成漩渦。文獻(xiàn)［24］表明，隨著 Re 數(shù)的增大，在流道的頂部也會出現(xiàn)小渦。

本文取后臺階流的尺寸h=1 m 。入口速度u in =1 m/s ，豎直方向速度v in =0，出口處為充分發(fā)展邊界，其余壁面為無滑移壁面，計算域被劃分為30×120的網(wǎng)格。該模型中，雷諾數(shù)定義如下

Re2=u in Hν（20）

式中：u in 為入口速度， m/s ；H為充分發(fā)展段的管道高度，m。

3.3 水箱流

在水箱流模型中，水流從水箱的右下方流入，經(jīng)過一個未封閉的腔體之后，從水箱的左上方流出。由于水流在流動的過程中，速度的方向發(fā)生了改變，所以在腔體內(nèi)也會形成漩渦，當(dāng) Re 數(shù)較大時，還會在水箱的角落形成小渦。

本文取水箱流的尺寸為1m×1m。入口速度v′ in =1 m/s ，水平方向速度u′ in =0。壁面均為無滑移壁面，計算域被劃分為40×40和80×80兩套網(wǎng)格。該模型中雷諾數(shù)定義如下

Re3=v′ in Lν（21）

式中：L為水箱的長度， m 。

4 計算結(jié)果及分析

本文基于SIMPLE算法編寫求解器，分別使用一階迎風(fēng)格式、QUICK格式和SGSD格式進(jìn)行離散，對頂蓋驅(qū)動流、后臺階流和水箱流算例進(jìn)行計算。在計算過程中對每個離散方程分別計算殘差，并認(rèn)為當(dāng)殘差滿足一定條件時，計算收斂。計算收斂條件為

max {Rm sum ，Ru sum ，Rv sum }lt;10^－6（22）

式中：Rm sum 為連續(xù)性方程的殘差；Ru sum 為X方向動量方程的殘差；Rv sum 為Y方向動量方程的殘差。

為了驗證計算結(jié)果的準(zhǔn)確性，首先將計算結(jié)果與3個經(jīng)典算例的已有研究結(jié)果^［24-26］進(jìn)行比較分析。在此基礎(chǔ)之上，為了全面評估自動獲取亞松弛因子方法相較于傳統(tǒng)固定亞松弛因子方法的優(yōu)勢，本文不僅記錄了自動獲取方法在各算例中的表現(xiàn)，同時也系統(tǒng)地統(tǒng)計了在亞松弛因子分別固定為0.1、0.2、0.3、0.4、0.5、0.6、0.7、0.8、0.9、1.0時的計算結(jié)果。本文的計算時間均為CPU時間（s），數(shù)字精確到小數(shù)點后3位。為了排除壓力松弛因子的干擾，本文所有計算的壓力亞松弛因子統(tǒng)一取0.2^［9］。

4.1 頂蓋驅(qū)動流的計算結(jié)果

當(dāng)Re1=1000、網(wǎng)格數(shù)為40×40時，采用一階迎風(fēng)格式、QUICK格式和SGSD格式分別進(jìn)行計算，計算時間如圖4所示。

圖4中部分亞松弛因子的固定值未呈現(xiàn)出對應(yīng)數(shù)據(jù)點，表示亞松弛因子在該固定值下，計算不收斂。下文相同，不再贅述。從圖4可以看出，在一階迎風(fēng)格式下，亞松弛因子取固定值時均能計算收斂，而在QUICK格式和SGSD格式下，當(dāng)亞松弛因子的取值過大時，則會出現(xiàn)不收斂的現(xiàn)象。

當(dāng) Re 1=1500、網(wǎng)格數(shù)為80×80時，3種離散格式所對應(yīng)的計算時間如圖5所示。

從圖4和圖5可以看出，在頂蓋驅(qū)動流算例中，總體上隨著亞松弛因子的增大，計算時間逐漸減小，但是在高階格式（QUICK和SGSD格式）中，當(dāng)亞松弛因子取值較大時，會出現(xiàn)計算不收斂的現(xiàn)象。當(dāng)采用自動獲取方法之后，各種離散格式下的計算均收斂，并且計算時間基本與采用固定亞松弛因子所需要的最短時間相當(dāng)。

4.2 后臺階流的計算結(jié)果

在后臺階流中，當(dāng) Re 2=500、網(wǎng)格數(shù)為30×120時，采用3種離散格式分別進(jìn)行計算，計算時間見圖6?？梢钥闯觯涸赒UICK格式下，當(dāng)亞松弛因子取值過小時，會出現(xiàn)發(fā)散的現(xiàn)象；在SGSD格式下，當(dāng)亞松弛因子取值過大時，也會出現(xiàn)發(fā)散現(xiàn)象。

在網(wǎng)格數(shù)不變的情況下，將 Re 2變?yōu)?000，3種離散格式所對應(yīng)的計算時間如圖7所示。

從圖6和圖7可以看出，在后臺階流算例中，總體上隨著亞松弛因子的增大，計算時間先減小后增大，在高階格式中，無論是亞松弛因子過小還是過大，都有可能出現(xiàn)計算不收斂的現(xiàn)象。當(dāng)采用自動獲取方法后，各種離散格式下的計算均收斂，并且計算時間基本與固定亞松弛因子所需要的最短時間相當(dāng)。

值得一提的是，在圖7中，QUICK格式下所有采用固定亞松弛因子的計算全部發(fā)散，但是在改用自動獲取方法之后，計算依然以較快的速度收斂，計算速度甚至追上了一階迎風(fēng)格式。

4.3 水箱流的計算結(jié)果

在水箱流中，當(dāng)Re3=500、網(wǎng)格數(shù)為40×40時，3種離散格式對應(yīng)的計算時間如圖8所示。

當(dāng) Re 3=1000、網(wǎng)格數(shù)為80×80時，3種離散格式所對應(yīng)的計算時間如圖9所示。

從圖8和圖9可以看出，在水箱流算例的兩種 Re 3數(shù)下，采用自動獲取方法相比于采用固定亞松弛因子，計算速度均體現(xiàn)出了顯著的優(yōu)勢。

4.4 網(wǎng)格數(shù)對計算的影響

在工程應(yīng)用中，選擇合適的網(wǎng)格數(shù)至關(guān)重要。從計算精度角度來講：網(wǎng)格數(shù)越多，對物理場的刻畫就越精確；網(wǎng)格數(shù)越少，對物理場的刻畫誤差就越大。從計算時間角度來講，在離散方程組的求解過程中，網(wǎng)格數(shù)直接決定了系數(shù)矩陣的大小，所以在其他條件相同的情況下，網(wǎng)格數(shù)越多，計算量越大，計算時間也就越長。為了進(jìn)一步驗證本文提出的自動獲取方法對不同網(wǎng)格數(shù)的適應(yīng)情況，對每種算例均選取了5組網(wǎng)格數(shù)，觀察自動獲取方法的計算表現(xiàn)。各算例對應(yīng)的網(wǎng)格數(shù)如表1所示。

在此基礎(chǔ)之上，本文選取了0.1、0.5、0.9這3個固定亞松弛因子作為對照組，來分析固定亞松弛因子方法的收斂性表現(xiàn)。根據(jù)4.1～4.3小節(jié)的計算結(jié)果可知：在頂蓋驅(qū)動流和水箱流中，0.1代表了收斂性最高的固定值，0.9代表了計算速度最快的固定值；在后臺階流中，0.5代表了計算速度最快的固定值。各算例的計算時間如圖10～12所示。

如圖10所示，在頂蓋驅(qū)動流中，總體上各個離散格式的計算收斂性均較好。在SGSD格式下，當(dāng)亞松弛因子較大時出現(xiàn)了計算不收斂的情況，相同情況下自動獲取方法的收斂情況較好，同時保持著較快的收斂速度。

如圖11所示，在后臺階流中，隨著網(wǎng)格數(shù)的變化，QUICK格式和SGSD格式均出現(xiàn)了計算不收斂的情況：QUICK格式在30×120和30×150網(wǎng)格下，當(dāng)亞松弛因子取0.1時，計算發(fā)散，在30×180和30×210網(wǎng)格下，3個固定亞松弛因子對應(yīng)的計算均不收斂；在SGSD格式下，當(dāng)亞松弛因子取0.9時，在5種網(wǎng)格數(shù)下的計算均發(fā)散。

采用自動獲取方法后，QUICK格式下的計算收斂性得到了顯著的改善，5種網(wǎng)格數(shù)下的計算均能夠收斂。通過與圖11（a）和圖11（c）的對比可知，在QUICK格式下，自動獲取方法依然保持了較快的計算速度。在一階迎風(fēng)格式和SGSD格式中，自動獲取方法的計算時間也均與亞松弛因子取0.5時的計算時間相當(dāng)，體現(xiàn)出了較快的計算速度。

如圖12所示，在水箱流中，隨著網(wǎng)格數(shù)的變化，各個格式的收斂性均較好，僅在計算時間上存在差別，自動獲取方法在各個離散格式下均保持著較快的計算速度。

對3個算例進(jìn)行橫向比較可以發(fā)現(xiàn)，在后臺階流中，固定亞松弛因子的計算收斂性最差。然而，當(dāng)采用自動獲取亞松弛因子的方法時，計算收斂性顯著提高，優(yōu)勢明顯，這一結(jié)果與后臺階流的流場分布有關(guān)。在后臺階流的流場中：入口段存在一個狹長的層流區(qū)域，流體速度的方向保持相對穩(wěn)定，變化不大；在靠近出口的充分發(fā)展段，流場也表現(xiàn)出較為平穩(wěn)的特性，流體速度在此區(qū)域基本保持不變。后臺階流的回流區(qū)主要集中在流道截面積突然擴(kuò)大的位置，這一區(qū)域內(nèi)流體速度的大小和方向都會經(jīng)歷顯著的變化。相較之下，頂蓋驅(qū)動流和水箱流的流場則沒有界限清晰的流速變化平緩區(qū)域和劇烈區(qū)域。這兩種算例的流場均由許多大小不一的漩渦組成，流速在這些區(qū)域內(nèi)持續(xù)變化。

前文已經(jīng)提及，在固定亞松弛因子方法中，亞松弛因子沒有考慮流場的速度分布情況，也無法針對流場的不同區(qū)域采取因地制宜的亞松弛策略。所以，在速度變化平緩區(qū)域和劇烈區(qū)域界限明顯的后臺階流中，固定亞松弛因子的這一缺點表現(xiàn)最為明顯，當(dāng)采用高階格式（QUICK格式和SGSD格式）時，出現(xiàn)了大面積計算不收斂的現(xiàn)象。與之形成鮮明對比的是，自動獲取亞松弛因子的方法可以精準(zhǔn)地捕獲速度變化劇烈的回流區(qū)，因地制宜地采取合適的亞松弛因子進(jìn)行計算，提高了計算速度。

4.5 殘差的變化趨勢

基于4.4小節(jié)對計算收斂性的討論，本文對殘差隨著迭代數(shù)的變化趨勢進(jìn)行研究。本文依然沿用了0.1、0.5、0.9這3個固定亞松弛因子，統(tǒng)計了在各算例下，殘差隨著迭代數(shù)的變化趨勢，并將其與自動獲取方法下的殘差曲線進(jìn)行了對比。為了完整展現(xiàn)殘差的最終狀態(tài)，在進(jìn)行計算時，刪除了程序里關(guān)于計算收斂的判斷語句，人為觀察殘差不再下降時停止計算，結(jié)果如圖13所示。

從圖13可以看出，在固定亞松弛因子的方法下，如果想要實現(xiàn)較好的收斂性，只能犧牲計算效率，減小亞松弛因子。相較之下，采用本文提出的自動獲取方法后，計算的收斂性顯著提高。同時，在3個算例中，自動獲取方法均表現(xiàn)出了較好的計算速度。在計算收斂的情況下，例如圖13（c）中，采用自動獲取方法或者固定亞松弛因子，殘差最終都能達(dá)到10^－14量級。

5 結(jié) 論

本文以流場速度梯度為依據(jù)，提出了一種亞松弛因子的自動獲取方法。通過3個經(jīng)典算例以及3種離散格式進(jìn)行計算，得到的結(jié)論如下。

（1）與固定亞松弛因子相比，本文提出的自動獲取方法在不同算例中均體現(xiàn)出了較快的計算速度，各工況下計算時間均與采用固定亞松弛因子的最短計算時間相當(dāng)，有些工況下甚至小于固定亞松弛因子的最短計算時間。

（2）自動獲取方法適用于多種離散格式。當(dāng)采用固定亞松弛因子時，如果亞松弛因子選擇不合理，會出現(xiàn)計算不收斂的現(xiàn)象，而采用自動獲取方法后，計算在收斂的同時可以得到理想的收斂速度。

（3）隨著網(wǎng)格數(shù)的增加，各個算例的計算時間均有不同程度增加。除此之外，網(wǎng)格數(shù)的改變對QUICK格式和SGSD格式下的后臺階流沖擊較大，造成了大面積計算不收斂的現(xiàn)象，而采用自動獲取方法可以顯著地改善計算不收斂的現(xiàn)象，有效地提高了計算收斂性。

（4）在計算收斂的情況下，本文的自動獲取方法可以達(dá)到與固定亞松弛因子相同的殘差量級。

（5）本文是以SIMPLE算法下的速度亞松弛因子為例提出的自動獲取方法。后續(xù)可以圍繞著適用于多種衍生算法的統(tǒng)一自動獲取方法以及速度壓力耦合的自動獲取方法開展進(jìn)一步研究。

參考文獻(xiàn)：

［1］宣領(lǐng)寬， 龔京風(fēng)， 張文平， 等. 用于非均質(zhì)復(fù)合材料應(yīng)力分析的交錯網(wǎng)格有限體積法 ［J］. 西安交通大學(xué)學(xué)報， 2014， 48（3）： 121-127.

XUAN Lingkuan， GONG Jingfeng， ZHANG Wenping， et al. Staggered cell-vertex finite volume method for analyzing stress in heterogeneous composite materials ［J］. Journal of Xi’an Jiaotong University， 2014， 48（3）： 121-127.

［2］VAN DOORMAAL J P， RAITHBY G D. Enhancements of the simple method for predicting incompressible fluid flows ［J］. Numerical Heat Transfer， 1984， 7（2）： 147-163.

［3］PATANKAR S V， SPALDING D B. A calculation procedure for heat， mass and momentum transfer in three-dimensional parabolic flows ［J］. International Journal of Heat and Mass Transfer， 1972， 15（10）： 1787-1806.

［4］王鋒， 孫中國， 劉啟新， 等. Y型微混合器內(nèi)流流場移動粒子半隱式法數(shù)值分析 ［J］. 西安交通大學(xué)學(xué)報， 2021， 55（5）： 162-170.

WANG Feng， SUN Zhongguo， LIU Qixin， et al. Numerical analysis of the flow field in a Y-type micromixer using moving particle semi-implicit method ［J］. Journal of Xi’an Jiaotong University， 2021， 55（5）： 162-170.

［5］李冬， 柳云釗， 陳鑫， 等. 酚醛樹脂基復(fù)合材料傳熱-熱解模型構(gòu)建與試驗研究 ［J］. 西安交通大學(xué)學(xué)報， 2024， 58（8）： 49-59.

LI Dong， LIU Yunzhao， CHEN Xin， et al. Research on construction and experimental verification of heat transfer-pyrolysis model for phenolic resin-based composite materials ［J］. Journal of Xi’an Jiaotong University， 2024， 58（8）： 49-59.

［6］王長安， 侯育杰， 劉成昌， 等. 燃煤鏈條爐摻燒半焦NO x 生成特性及配風(fēng)優(yōu)化數(shù)值模擬研究 ［J］. 西安交通大學(xué)學(xué)報， 2022， 56（4）： 1-12.

WANG Changan， HOU Yujie， LIU Chengchang， et al. Numerical simulation study on NO x "emission characteristics and the optimization of air induction of coal-fired chain boiler blended with semi-coke ［J］. Journal of Xi’an Jiaotong University， 2022， 56（4）： 1-12.

［7］鄧雅軍， 張琳， 孫東亮， 等. IDEAL算法在OpenFOAM界面兩相流求解器中的實現(xiàn)與性能分析 ［J］. 西安交通大學(xué)學(xué)報， 2020， 54（12）： 147-157.

DENG Yajun， ZHANG Lin， SUN Dongliang， et al. Implementation and performance analysis of IDEAL algorithm for unsteady two-phase flow problems in OpenFOAM ［J］. Journal of Xi’an Jiaotong University， 2020， 54（12）： 147-157.

［8］PATANKAR S V. Numerical heat transfer and fluid flow ［M］. Boca Raton， FL， USA： CRC Press， 1980.

［9］BARRON R M， SALEHI NEYSHABOURI A A. Effects of under-relaxation factors on turbulent flow simulations ［J］. International Journal for Numerical Methods in Fluids， 2003， 42（8）： 923-928.

［10］MAREK R， STRAUB J. Hybrid relaxation： a technique to enhance the rate of convergence of iterative algorithms ［J］. Numerical Heat Transfer： Part BFundamentals， 1993， 23（4）： 483-497.

［11］GOPINATH R， GANESAN V. Orthogonal arrays： An introduction and their application in optimizing underrelaxation factors in ‘simple’-based algorithm ［J］. International Journal for Numerical Methods in Fluids， 1992， 14（6）： 665-680.

［12］CHATWANI AV， TURAN A. Technical note improved pressure-velocity coupling algorithm based on minimization of global residual norm ［J］. Numerical Heat Transfer： Part B Fundamentals， 1991， 20（1）： 115-123.

［13］范曉偉， 陰建民， 陳鐘頎. SIMPLE算法中壓力松馳因子的自動生成 ［J］. 鄭州紡織工學(xué)院學(xué)報， 1997， 8（3）： 7-12.

FAN Xiaowei， YIN Jianmin， CHEN Zhongqi. The auto production of pressure relaxation factor in simple algorithm ［J］. Journal of Zhengzhou Textile Institute， 1997， 8（3）： 7-12.

［14］LIU Xunliang， TAO Wenquan， ZHENG Ping， et al. Control of convergence in a computational fluid dynamic simulation using fuzzy logic ［J］. Science in China： Series E Technological Science， 2002， 45（5）： 495-502.

［15］王艷寧， 孫東亮， 苗政， 等. 模糊控制技術(shù)在SIMPLER算法中的應(yīng)用及求解性能分析 ［J］. 西安交通大學(xué)學(xué)報， 2016， 50（1）： 78-84.

WANG Yanning， SUN Dongliang， MIAO Zheng， et al. Application of a fuzzy control method in SIMPLER algorithm and its solving performance analysis ［J］. Journal of Xi’an Jiaotong University， 2016， 50（1）： 78-84.

［16］MORII T， VIEROW K. The soar method for automatically optimizing simple relaxation factors ［J］. Numerical Heat Transfer： Part B Fundamentals， 2000， 38（3）： 309-332.

［17］馬淼， 李春光， 景何仿. 確定SIMPLE算法中壓力松弛因子的自適應(yīng)方法 ［J］. 甘肅科學(xué)學(xué)報， 2011， 23（4）： 102-105.

MA Miao， LI Chunguang， JING Hefang. The adaptive method for fixing the pressure relaxation factor of SIMPLE algorithm ［J］. Journal of Gansu Sciences， 2011， 23（4）： 102-105.

［18］YOU Wei， LI Zengyao， TAO Wenquan. A general self-adaptive under-relaxation strategy for fast and robust convergence of iterative calculation of incompressible flow ［J］. Numerical Heat Transfer： Part B Fundamentals， 2020， 77（4）： 299-310.

［19］宇波， 焦開拓， 陳宇杰， 等. 有限容積法對流項離散格式的穩(wěn)定性與有界性研究 ［J］. 西安交通大學(xué)學(xué)報， 2023， 57（4）： 107-121.

YU Bo， JIAO Kaituo， CHEN Yujie， et al. Study on the stability and boundedness of discrete schemes for convection term based on finite volume method ［J］. Journal of Xi’an Jiaotong University， 2023， 57（4）： 107-121.

［20］陶文銓. 數(shù)值傳熱學(xué) ［M］. 2版. 西安： 西安交通大學(xué)出版社， 2001.

［21］王福軍. 計算流體動力學(xué)分析 ［M］. 北京： 清華大學(xué)出版社， 2004.

［22］NI Mingjiu， TAO Wenquan， WANG Shangjin. Stability-controllable second-order difference scheme for convection term ［J］. Journal of Thermal Science， 1998， 7（2）： 119-130.

［23］LI Z Y， TAO W Q. A new stability-guaranteed second-order difference scheme ［J］. Numerical Heat Transfer： Part B Fundamentals， 2002， 42（4）： 349-365.

［24］ERTURK E. Numerical solutions of 2-D steady incompressible flow over a backward-facing step： part I high Reynolds number solutions ［J］. Computers amp; Fluids， 2008， 37（6）： 633-655.

［25］GHIA U， GHIA K N， SHIN C T. High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method ［J］. Journal of Computational Physics， 1982， 48（3）： 387-411.

［26］MAZHAR Z. A procedure for the treatment of the velocity-pressure coupling problem in incompressible fluid flow ［J］. Numerical Heat Transfer： Part B Fundamentals， 2001， 39（1）： 91-100.

（編輯 陶晴）