旋轉(zhuǎn)坐標系N-S方程無網(wǎng)格/網(wǎng)格混合算法

曹 騁 ，陳紅全 ，張加樂 ，高煜堃

（1.常熟理工學院機械工程學院，江蘇 常熟 215500；2.南京航空航天大學航空學院，南京 210016；3.安徽工業(yè)大學機械工程學院，安徽 馬鞍山 243002）

0 引言

槳葉是渦槳發(fā)動機最為重要的氣動部件之一，其旋轉(zhuǎn)流場也是計算流體力學（Computational fluid dy?namics，CFD）數(shù)值模擬的難點。槳葉氣動特性的優(yōu)劣直接影響發(fā)動機性能。因此，發(fā)展能準確模擬旋轉(zhuǎn)槳葉氣動特性的CFD 算法，對渦槳發(fā)動機設(shè)計與研究具有重要的工程意義。

類似發(fā)動機槳葉這樣的旋轉(zhuǎn)體繞流在工程中廣泛存在，如旋翼飛行器等。傳統(tǒng)的CFD 技術(shù)在模擬這類流場時會遇到2 個難點。一個是目前的CFD 軟件大多基于“網(wǎng)格算法”，涉及計算域網(wǎng)格生成。但對于存在旋轉(zhuǎn)槳葉的渦槳發(fā)動機等多體復(fù)雜氣動外形，要生成合適的整體網(wǎng)格，刻畫出局部間隙等幾何細節(jié)，在網(wǎng)格拓撲約束下并非易事，常需要用到塊結(jié)構(gòu)等特殊處理技術(shù)；另一個是因旋轉(zhuǎn)產(chǎn)生的動邊界問題的處理。目前的求解方法通?；趧泳W(wǎng)格技術(shù)開展，在處理動邊界大位移問題時，常需要結(jié)合網(wǎng)格重構(gòu)、網(wǎng)格重疊等特殊處理技術(shù)，涉及到與角速度相關(guān)的旋轉(zhuǎn)特征物理量的插值運算，這會影響到算法的計算效率與精度。為了突破網(wǎng)格拓撲的制約，Batina等、蒲賽虎、郭曼麗等提出1 類無網(wǎng)格算法，計算域用無網(wǎng)格點離散，既可直接布點，也可用現(xiàn)有的結(jié)構(gòu)或非結(jié)構(gòu)網(wǎng)格點，涉及的空間導(dǎo)數(shù)等物理量僅依據(jù)當?shù)攸c云結(jié)構(gòu)逼近確定，其靈活性更適合處理旋轉(zhuǎn)體等復(fù)雜的幾何外形；為了準確反映旋轉(zhuǎn)角速度等特征物理量，Chen等發(fā)展了一種直接求解旋轉(zhuǎn)坐標系下控制方程的方法。旋轉(zhuǎn)運動網(wǎng)格在旋轉(zhuǎn)坐標系可與物體保持相對靜止，無須反復(fù)重構(gòu)，旋轉(zhuǎn)特征物理量等也無須插值運算；像渦槳發(fā)動機的旋轉(zhuǎn)槳葉、懸停旋翼等非定常繞流，在旋轉(zhuǎn)坐標系下還可被看作準定常，減少了求解的復(fù)雜度?；谶@些特點，該方法已在渦輪機械和旋翼繞流的模擬中獲得廣泛應(yīng)用。但將無網(wǎng)格算法和旋轉(zhuǎn)坐標系下控制方程結(jié)合起來應(yīng)用于求解槳葉這樣的旋轉(zhuǎn)體繞流問題，還鮮有文獻報道。

本文發(fā)展一種基于旋轉(zhuǎn)坐標系的無網(wǎng)格/網(wǎng)格混合算法，用于求解槳葉等旋轉(zhuǎn)物體的黏性繞流問題。基于無網(wǎng)格點云空間導(dǎo)數(shù)逼近等方法離散控制方程，并結(jié)合隱式LU-SGS 算法，發(fā)展形成求解非定常問題的雙時間步隱式格式推進求解的混合算法，選取2 維振蕩翼型非定常繞流、旋轉(zhuǎn)圓柱黏性繞流和模擬發(fā)動機槳葉旋轉(zhuǎn)運動的3 維懸停旋翼繞流進行考核運算，并與文獻計算或試驗結(jié)果進行比較分析。

1 控制方程

2 無網(wǎng)格/網(wǎng)格混合算法

2.1 無網(wǎng)格算法

2.1.1 空間導(dǎo)數(shù)擬合

計算域無網(wǎng)格布點離散后，形成點云結(jié)構(gòu)，點云()由中心點及其周邊衛(wèi)星點組成。空間導(dǎo)數(shù)基于當?shù)攸c云逼近確定。沿用文獻[12]的做法，選用較為典型的加權(quán)最小二乘法，任一函數(shù)在點云中心點的1階擬合導(dǎo)數(shù)為

式中：α、β和γ為只與點云幾何位置信息相關(guān)的無網(wǎng)格系數(shù)；f為中心點與第個衛(wèi)星點連線中點處的函數(shù)近似值。

2.1.2 空間離散

易于發(fā)現(xiàn)，式（1）與慣性坐標系下控制方程的主要區(qū)別在于包含矢通量、和的項及源項。對于前者，用式（3）可得

其中

分裂得到的L與在慣性坐標系下無網(wǎng)格算法對應(yīng)的離散形式一致，可按慣性坐標系已有方法處理，最大程度地減少對已有慣性坐標系下無網(wǎng)格求解程序的改動。W在虛擬交界面上用中心格式確定， |中的速度分量則取為中心點和衛(wèi)星點連線中點處的牽連速度分量。

源項中涉及的密度和速度等物理量直接取中心點處的值即可。式（1）中黏性矢通量部分的求解與慣性坐標系下的相同，其中的湍流黏性系數(shù)采用Spalart-Allmaras湍流模型計算，不再贅述。

2.1.3 時間離散

對式（1）進行空間離散后得到的半離散方程為

式中：R()為殘值，包含了對流矢通量、黏性矢通量和源項經(jīng)空間離散后的部分。

對于非定常流動，式（8）可通過引入虛擬時間和非定常殘值，利用顯式雙時間步在虛擬時間層上推進，以求得下一物理時間步的解W。本文結(jié)合LU-SGS 格式形成隱式雙時間步推進算法。沿用Yoon等的思路，最終得到的2步推進格式為

式中：為對下一物理時間步守恒變量值W的逼近，即當虛擬時間層→時，有→W；()和()為點云中心點的衛(wèi)星點集分割出來的2個子集。

需要注意的是式（12）中包含了源項產(chǎn)生的部分

上述算法將與網(wǎng)格算法集成，結(jié)合第2.3節(jié)中的跨區(qū)交界面處理技術(shù)，發(fā)展形成基于LU-SGS 隱式推進的無網(wǎng)格/網(wǎng)格混合算法。

2.2 邊界條件

以上的速度矢量均為絕對速度矢量。鏡像點處的壓強和密度取為對應(yīng)原像點處的值。

遠場則基于帶源項的特征分析無反射邊界條件確定。

2.3 無網(wǎng)格/網(wǎng)格混合算法實施

計算域采用整體網(wǎng)格和物面附近局部無網(wǎng)格離散形成混合算法，實施的關(guān)鍵是無網(wǎng)格和網(wǎng)格分區(qū)交界面的處理。本文通過引入無網(wǎng)格和網(wǎng)格對偶點，實現(xiàn)了混合算法繞流信息的跨區(qū)交換。為了圖示清晰，以2維為例加以說明。對偶點一般取為交界面處2層網(wǎng)格的格點，如圖1所示。

第1 層稱為無網(wǎng)格對偶點（圖1 空心圓點），第2層稱為網(wǎng)格對偶點（圖1方點）。這些對偶點既被當做點云中的無網(wǎng)格點，也被當作網(wǎng)格單元的格點。作為前者時，純無網(wǎng)格點（圖1 實心圓點）和無網(wǎng)格對偶點周圍點云模板完整（圖1 虛線連接的點云結(jié)構(gòu)），可用無網(wǎng)格算法計算；作為后者時，純網(wǎng)格節(jié)點和網(wǎng)格對偶點周圍網(wǎng)格模板完整（格點格式下），可用格點格式有限體積法求解。2 種算法均將物理量存儲于對偶點處，相互間可直接取用以實現(xiàn)跨區(qū)信息交換，無須插值。

由上述可知，采用原有方法挑選的對偶點，網(wǎng)格部分需要采用格點格式求解。而在求解非定常流動時需要額外的計算修正，涉及到所謂的質(zhì)量矩陣運算。為了避免該問題，本文利用無網(wǎng)格取點的靈活性調(diào)整了對偶點的選取方法。將對偶點選取在網(wǎng)格區(qū)域靠近無網(wǎng)格區(qū)的若干層單元的中心，如圖2 所示（圖中各點含義同圖1）。

圖1 調(diào)整前選取的對偶點

圖2 調(diào)整后選取的對偶點

調(diào)整位置后，存儲于對偶點處的物理量亦可認為存儲于相應(yīng)單元的中心，在計算時可被當作相應(yīng)單元中心處的值直接取用，因此仍然無須插值便可實現(xiàn)跨區(qū)信息交換。但同時，除無網(wǎng)格對偶點所對應(yīng)的單元外（網(wǎng)格模板不完整，但點云模板完整，使用無網(wǎng)格算法計算），其余網(wǎng)格均可統(tǒng)一采用格心格式的有限體積法統(tǒng)一求解，進而避免質(zhì)量矩陣的計算。數(shù)值計算表明，上述處理方法能滿足跨區(qū)交界面數(shù)值通量傳遞的要求，表現(xiàn)為跨區(qū)等值線的連續(xù)（見第3.1節(jié)圖7）。

3 結(jié)果分析

本文對上述發(fā)展的混合算法編制了對應(yīng)的計算程序，并進行了多個典型算例的考核運算。這里先給出NACA 0012翼型的俯仰振蕩非定常動邊界繞流和旋轉(zhuǎn)圓柱黏性繞流算例，再給出旋翼懸停飛行繞流算例。

3.1 NACA 0012翼型俯仰振蕩繞流

NACA 0012翼型俯仰振蕩繞流算例源自AGARD報告，其中翼型的俯仰可視為旋轉(zhuǎn)運動，已被許多文獻用于考核求解旋轉(zhuǎn)體繞流問題。初始時翼型位于平均攻角0.016°處，而后繞25%弦長點作俯仰振蕩，最大振幅對應(yīng)攻角= 2.51°，來流馬赫數(shù)= 0.755，雷諾數(shù)= 5.5×10，試驗振蕩頻率= 62.5 Hz，攻角()隨時間變化為

計算采用的局部無網(wǎng)格區(qū)和網(wǎng)格區(qū)如圖3 所示。遠場邊界落在離物體10倍弦長處，涉及16320個無網(wǎng)格點和4331 個網(wǎng)格單元。為準確捕捉黏性流動的特征，物面附近的無網(wǎng)格區(qū)采用了沿法向的布點方式，第1層無網(wǎng)格點距離物面高度為1×10。

圖3 NACA 0012翼型周圍的無網(wǎng)格區(qū)和網(wǎng)格區(qū)

升力系數(shù)和力矩系數(shù)隨瞬時攻角的變化如圖4、5所示。從圖中可見遲滯環(huán)現(xiàn)象。圖4中的升力系數(shù)相對試驗值整體偏低，但與文獻中的計算值吻合得較好；對于力矩系數(shù)，所給結(jié)果只在個別試驗值處有較大偏差，但整體與試驗值和計算值都符合較好。4個典型時刻翼型表面的壓強系數(shù)分布如圖6所示，相應(yīng)的試驗值也在圖中一并給出，供比較。不同時刻流場的壓強系數(shù)等值線如圖7 所示。展示出翼型俯仰振蕩時，激波等流場特征在翼型上下表面交替演化的情況。

圖4 升力系數(shù)隨攻角變化

圖5 力矩系數(shù)隨攻角變化

圖6 不同振蕩位置翼型表面對應(yīng)壓強系數(shù)

圖7 不同時刻流場的壓強系數(shù)等值線

3.2 旋轉(zhuǎn)圓柱黏性繞流

旋轉(zhuǎn)圓柱黏性繞流在增升減阻、改變旋渦釋放的規(guī)律、降低流動引起的振動等方面具有重要理論價值，歷來為流體力學學者所重視。Coutancea 等也對該流動現(xiàn)象進行了大量觀察和測量，為數(shù)值模擬提供了豐富的參考依據(jù)。這里給出本算例發(fā)展的算法獲得的運算結(jié)果。

設(shè)無窮遠處來流速度為，方向沿軸正向，圓柱直徑為，從靜止開始突然以角速度按逆時針方向旋轉(zhuǎn)，圓柱旋轉(zhuǎn)的角速度= 2。其中為轉(zhuǎn)速比，本文沿用文獻[21]取=0.5。來流條件為=0.1，=200。設(shè)無量綱時間=(，本文計算了∈[0,4]范圍內(nèi)的流場情況。計算采用的無網(wǎng)格區(qū)和網(wǎng)格區(qū)如圖8所示。

圖8 圓柱周圍的無網(wǎng)格區(qū)和網(wǎng)格區(qū)

由于外形規(guī)則，網(wǎng)格區(qū)使用了4500 個結(jié)構(gòu)網(wǎng)格單元，無網(wǎng)格區(qū)則填充了11476 個無網(wǎng)格點，遠場邊界落在離物體10 倍圓柱直徑處。物理時間步長取為總物理時長的1/800，即Δ=/800=/200。

不同無量綱時刻的計算流線和試驗流場照片如圖9 所示。從圖中可見，本文發(fā)展算法所得結(jié)果與試驗結(jié)果基本一致，表現(xiàn)在捕捉的旋渦在強度、位置等特征上與對應(yīng)時刻試驗中觀察到的現(xiàn)象基本吻合。3 個無量綱時刻背風區(qū)上的方向速度分布如圖10 所示。從圖中可見，各時刻的速度分布曲線與試驗值均較為接近。

圖9 各時刻流線與試驗流場比較

圖10 背風區(qū)y = 0 上x 方向速度

3.3 懸停旋翼黏性繞流

懸停旋翼因繞固定軸旋轉(zhuǎn)，相當于發(fā)動機旋轉(zhuǎn)槳葉，因此，給出更具挑戰(zhàn)的3 維Caradonna 旋翼模型懸停飛行算例來模擬發(fā)動機旋轉(zhuǎn)槳葉繞流。該模型由2 片槳葉構(gòu)成，槳葉平面呈長方形，展弦比為6。槳葉剖面為NA?CA0012翼型，沿展向無扭轉(zhuǎn)。計算狀態(tài)取為：槳尖馬赫數(shù)= 0.800（約2268 r/min），槳距角θ= 5°，雷諾數(shù)= 3.58×10。該旋轉(zhuǎn)引起的非定常流動在本文采用的旋轉(zhuǎn)坐標系下可視為準定常流動，因此，計算可按準定常進行。采用網(wǎng)格和無網(wǎng)格混合算法，計算區(qū)域離散整體用了235379 個網(wǎng)格單元和物面附近1200192 無 網(wǎng)格點，如圖11 所示。從圖中可見3 維Caradonna 旋翼模型和典型截面無網(wǎng)格和網(wǎng)格離散情況。為了準確捕捉黏性流動特性，法向第1 層無網(wǎng)格布點與物面距離為1×10，遠場落在離物體15倍的展向翼型弦長處。

圖11 旋翼表面及典型截面無網(wǎng)格和網(wǎng)格離散

在不同展向位置=0.50、0.68、0.80、0.96 的表面壓強系數(shù)計算值和試驗值比較如圖12 所示。從圖中可見，除了在個別點處與試驗值有所偏差外，計算結(jié)果與試驗值總體吻合，分布趨勢一致。

圖12 旋翼不同展向位置表面壓強系數(shù)

4 結(jié)論

（1）發(fā)展的混合算法是直接基于3 維旋轉(zhuǎn)坐標系下的N-S控制方程，相較于傳統(tǒng)的基于固定坐標系下的動網(wǎng)格算法，網(wǎng)格和局部無網(wǎng)格點云隨物體旋轉(zhuǎn)而無須重構(gòu)，旋轉(zhuǎn)角速度等重要的特征物理量也無須插值運算，有效避免了因重構(gòu)、插值等操作引起的誤差；

（2）混合算法繞流信息的跨區(qū)傳遞可利用無網(wǎng)格區(qū)點云取點的靈活性，借助于本文引入的無網(wǎng)格和網(wǎng)格對偶點實現(xiàn)。算例圖示等值線能跨區(qū)光滑過渡，證明本文發(fā)展的對偶點調(diào)整選取方法切實可行；

（3）發(fā)展的算法成功結(jié)合了LU-SGS 時間推進，并通過了模擬發(fā)動機槳葉旋轉(zhuǎn)運動的3 維懸停旋翼等典型復(fù)雜非定常繞流算例的考核運算，數(shù)值模擬結(jié)果與文獻計算或試驗結(jié)果吻合，展示出算法在模擬發(fā)動機槳葉等旋轉(zhuǎn)部件繞流問題方面具有重要的工程應(yīng)用前景。