曹 騁 ，陳紅全 ，張加樂(lè) ，高煜堃

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

0 引言

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

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

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

1 控制方程

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

2.1 無(wú)網(wǎng)格算法

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

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

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

2.1.2 空間離散

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

其中

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

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

2.1.3 時(shí)間離散

對(duì)式（1）進(jìn)行空間離散后得到的半離散方程為

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

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

式中：為對(duì)下一物理時(shí)間步守恒變量值W的逼近，即當(dāng)虛擬時(shí)間層→時(shí)，有→W；()和()為點(diǎn)云中心點(diǎn)的衛(wèi)星點(diǎn)集分割出來(lái)的2個(gè)子集。

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

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

2.2 邊界條件

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

遠(yuǎn)場(chǎng)則基于帶源項(xiàng)的特征分析無(wú)反射邊界條件確定。

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

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

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

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

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

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

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

3 結(jié)果分析

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

3.1 NACA 0012翼型俯仰振蕩繞流

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

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

圖3 NACA 0012翼型周?chē)臒o(wú)網(wǎng)格區(qū)和網(wǎng)格區(qū)

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

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

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

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

圖7 不同時(shí)刻流場(chǎng)的壓強(qiáng)系數(shù)等值線(xiàn)

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

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

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

圖8 圓柱周?chē)臒o(wú)網(wǎng)格區(qū)和網(wǎng)格區(qū)

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

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

圖9 各時(shí)刻流線(xiàn)與試驗(yàn)流場(chǎng)比較

圖10 背風(fēng)區(qū)y = 0 上x(chóng) 方向速度

3.3 懸停旋翼黏性繞流

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

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

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

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

4 結(jié)論

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

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

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