基于OpenMP的旋翼氣動(dòng)噪聲并行計(jì)算

任明霞，楊愛明

(復(fù)旦大學(xué) 航空航天系，上海 200433)

相對(duì)于固定翼飛機(jī)，直升機(jī)因其獨(dú)特的起降方式，在軍用和民用領(lǐng)域得到了廣泛的應(yīng)用.但是直升機(jī)面臨的一個(gè)重要問題就是噪聲問題.在民用領(lǐng)域，該問題嚴(yán)重影響乘坐舒適性和地面環(huán)境；在軍用領(lǐng)域，該問題嚴(yán)重影響直升機(jī)的隱身性能.旋翼噪聲的精確數(shù)值預(yù)測(cè)方法作為低噪聲旋翼優(yōu)化設(shè)計(jì)的先導(dǎo)技術(shù)，無疑具有重要的工程應(yīng)用價(jià)值.

圖1 旋翼噪聲計(jì)算技術(shù)流程圖Fig.1 Technology flowchart of rotor noise calculation

如果觀測(cè)點(diǎn)距離直升機(jī)較近，直接利用計(jì)算流體動(dòng)力學(xué)(Computational Fluid Dynamics, CFD)方法計(jì)算觀測(cè)點(diǎn)處的噪聲是可行的.但是對(duì)于遠(yuǎn)場(chǎng)噪聲問題，觀測(cè)點(diǎn)與直升機(jī)距離較遠(yuǎn)，直接應(yīng)用CFD方法計(jì)算遠(yuǎn)場(chǎng)噪聲問題是不現(xiàn)實(shí)的，其主要原因在于數(shù)值耗散和數(shù)值色散的影響.目前主流的噪聲計(jì)算方法是一種將CFD與聲比擬相結(jié)合的方法.在近場(chǎng)，通過數(shù)值方法求解Euler方程或Navier-Stokes(N-S)方程，得到聲源控制面上的信息；在遠(yuǎn)場(chǎng)通過求解描述聲波傳播的波動(dòng)方程得到遠(yuǎn)場(chǎng)的壓力波動(dòng).Ffowcs Williams和Hawkings[1]于1969年在Lighthill聲比擬研究的基礎(chǔ)上，將N-S方程重新整理成非齊次波動(dòng)方程的形式，得到了著名的FW-H方程.該方程可以描述運(yùn)動(dòng)物體產(chǎn)生的遠(yuǎn)場(chǎng)噪聲問題，目前旋翼噪聲的數(shù)值預(yù)測(cè)大多采用該方程.由于FW-H方程是一個(gè)微分方程，直接求解比較困難，NASA的科學(xué)家Farassat[2]提出了一種積分方法，得到了FW-H方程的兩個(gè)求解公式，即F1和F1A.該積分方法甚至可以將聲源面取為可穿透的控制面.本文對(duì)旋翼遠(yuǎn)場(chǎng)噪聲的數(shù)值計(jì)算采用可穿透的控制面，通過F1A公式進(jìn)行求解.

根據(jù)Farassat的積分公式求解FW-H方程，需要先得到直升機(jī)旋翼的流場(chǎng)信息，圖1為旋翼噪聲計(jì)算的技術(shù)路線示意圖.在本文工作中，旋翼流場(chǎng)的數(shù)值計(jì)算采用格心格式的有限體積方法[3],時(shí)間推進(jìn)采用雙時(shí)間推進(jìn)方法[4-5]，子迭代由LU-SSOR方法[6]完成.旋翼渦流場(chǎng)的CFD計(jì)算中，由于隱式LU-SSOR方法在每個(gè)迭代過程中占據(jù)主要的計(jì)算時(shí)間，因此對(duì)LU-SSOR算法進(jìn)行并行計(jì)算無疑可以提高旋翼渦流場(chǎng)的計(jì)算效率，大大縮減程序運(yùn)行時(shí)間.

并行編程是使用程序語(yǔ)言顯式地進(jìn)行說明，從而實(shí)現(xiàn)將計(jì)算任務(wù)中不同部分分配給不同的CPU同時(shí)執(zhí)行.目前，工程技術(shù)人員常用的編程語(yǔ)言是消息傳遞接口[7](Message Passing Interface, MPI)和直接控制共享內(nèi)存式并行編程的應(yīng)用程序接口[8](Open Multi-Processing, OpenMP).MPI需要明確劃分?jǐn)?shù)據(jù)結(jié)構(gòu)并重構(gòu)源程序，編程困難且開發(fā)周期長(zhǎng)，而OpenMP主要針對(duì)循環(huán)進(jìn)行并行，能有效克服并行編程的可移植性和擴(kuò)散性能差的缺點(diǎn).基于此，本文選用OpenMP進(jìn)行源程序的并行優(yōu)化,針對(duì)隱式LU-SSOR算法存在較強(qiáng)數(shù)據(jù)依賴性的問題，發(fā)展了兩種并行策略.

1 數(shù)值方法

1.1 流場(chǎng)計(jì)算方法

在絕對(duì)坐標(biāo)系下，積分形式的RANS方程為：

(1)

式中：W為守恒變量；Hc和Hv為對(duì)流通量和粘性通量.表達(dá)式如下：

(2)

(3)

式(2)和(3)中：ρ，(u,v,w)，p，E，H，k，T和τ分別為流體密度，流體速度在3個(gè)坐標(biāo)方向上的分量，壓強(qiáng)，單位體積的總能，單位體積的總焓，熱傳導(dǎo)系數(shù)，溫度和粘性應(yīng)力；V和Vb分別為流體速度矢量和網(wǎng)格速度矢量.為了封閉方程，還需要引入如下的狀態(tài)方程：

p=(γ-1)ρ[E-0.5(u2+v2+w2)].

(4)

在計(jì)算旋翼的懸停流場(chǎng)時(shí)，由于在旋轉(zhuǎn)坐標(biāo)系下觀察流體的絕對(duì)運(yùn)動(dòng)，流動(dòng)是定常的.因此將控制方程在旋轉(zhuǎn)坐標(biāo)系下表達(dá)無疑是方便的.引入如下公式：

(5)

式中：

B=Bxi+Byj+Bzk=Br=Bxrir+Byrjr+Bzrkr,

(6)

(7)

(i,j,k)和(ir,jr,kr)分別為絕對(duì)坐標(biāo)系和旋轉(zhuǎn)坐標(biāo)系下的坐標(biāo)軸方向單位矢量；ω為旋轉(zhuǎn)坐標(biāo)系相對(duì)于絕對(duì)坐標(biāo)系的旋轉(zhuǎn)角速度矢量.將式(5)應(yīng)用于絕對(duì)坐標(biāo)系下的控制方程(1)中的3個(gè)動(dòng)量方程，則可以得到如下的旋轉(zhuǎn)坐標(biāo)系下的控制方程：

(8)

(9)

與絕對(duì)坐標(biāo)系下的控制方程(1)相比，方程(8)中多了一項(xiàng)源項(xiàng).方程(9)中各物理量的介紹在此省略.式(8)和(9)中，下標(biāo)r代表旋轉(zhuǎn)坐標(biāo)系.值得指出的是，式(8)和(9)中所有矢量均為絕對(duì)矢量，僅對(duì)其分量在旋轉(zhuǎn)坐標(biāo)系下表達(dá)而已.

旋翼流場(chǎng)的數(shù)值計(jì)算采用格心格式的有限體積方法.對(duì)于任意控制體Vi,j,k，應(yīng)用式(1)可得

(10)

(11)

式中：

(12)

1.2 時(shí)間推進(jìn)

用二階向后差分近似式(10)中的物理時(shí)間導(dǎo)數(shù)項(xiàng)可得

(13)

為了得到物理時(shí)間t上的時(shí)間精確解Wn+1，我們?cè)谑?13)中引入守恒變量對(duì)偽時(shí)間τ的導(dǎo)數(shù)項(xiàng)，則可以得到：

(14)

顯然，偽時(shí)間方向上的定常收斂解即為物理時(shí)間上的時(shí)間精確解.為了方便起見，我們將待求變量Wn+1重新記為W*，則式(14)可改寫為：

(15)

式(15)在偽時(shí)間方向上的時(shí)間推進(jìn)格式采用隱式LU-SSOR方法.首先采用一階向后差分近似偽時(shí)間導(dǎo)數(shù)項(xiàng)

(16)

然后將殘值R進(jìn)行線化處理(僅處理無粘部分)，即：

(17)

將式(17)代入式(16)，則可得到：

(18)

式中：A，B和C分別為3個(gè)方向上的Jacobian矩陣，即

交界面處的Jacobian矩陣按照如下的迎風(fēng)格式處理：

(19)

其他方向依此類推.式(19)中A±的定義為：

式中：ρA為A的譜半徑；σ為一個(gè)略大于1的常數(shù)；I為單位矩陣.上述定義可以保證A+的所有信息向正方向傳播，而A-的所有信息向負(fù)方向傳播.定義如下3個(gè)算子：

(20)

則式(18)可以寫成如下的算子形式：

(L+D+U)ΔW*=RHS,

(21)

RHS由下式定義：

(22)

又因?yàn)?/p>

(L+D+U)ΔW*≈(L+D)D-1(U+D)ΔW*,

(23)

則我們可以得到如下的方程：

(L+D)D-1(U+D)ΔW*=RHS.

(24)

該方程的求解可以分為兩步，即上掃和下掃過程：

更新 (W*)m+1=(W*)m+1+ΔW*.

值得指出的是，為了提高計(jì)算效率，一般將偽時(shí)間步長(zhǎng)取為無窮大，則子迭代過程近似為高效的Newton迭代.

2 OpenMP并行計(jì)算

由以上推導(dǎo)過程可知，LU-SSOR方法存在很強(qiáng)的數(shù)據(jù)依賴性，因此不能利用一般的OpenMP并行語(yǔ)句進(jìn)行并行處理.基于此，本文根據(jù)NAS并行基準(zhǔn)[9]發(fā)展了超平面(Hyperplane)和管道流(Pipelining)兩種并行策略.

2.1 超平面策略

由于隱式LU-SSOR存在很強(qiáng)的數(shù)據(jù)依賴性，因此不能簡(jiǎn)單地在I/J/K任意一個(gè)方向上做OpenMP并行計(jì)算.超平面算法的思想是構(gòu)造平面M=I+J+K，計(jì)算沿著構(gòu)造平面進(jìn)行，然后將任務(wù)在J方向上分配給不同的線程，實(shí)現(xiàn)并行計(jì)算.

以二維情況為例，如圖2，超平面方法將計(jì)算空間劃分為一個(gè)個(gè)的平面，M=I+J，然后在J方向上進(jìn)行并行計(jì)算.實(shí)心圓為流場(chǎng)信息已知的網(wǎng)格點(diǎn)，空心圓為待求解流場(chǎng)信息的網(wǎng)格點(diǎn).以點(diǎn)(Ia,Ja)為例，在用雙時(shí)間算法進(jìn)行推進(jìn)過程中，點(diǎn)(Ia,Ja)的流場(chǎng)值依賴于點(diǎn)(Ia-1,Ja)和點(diǎn)(Ia,Ja-1).超平面算法可以保證點(diǎn)(Ia-1,Ja)和點(diǎn)(Ia,Ja-1)的流場(chǎng)值先于點(diǎn)(Ia,Ja)更新，從而解決了數(shù)據(jù)依賴性問題.即對(duì)于任意一個(gè)超平面L，其上的網(wǎng)格點(diǎn)數(shù)值僅僅與上一個(gè)超平面(L-1)上的網(wǎng)格點(diǎn)有關(guān)，超平面內(nèi)的網(wǎng)格點(diǎn)之間沒有數(shù)據(jù)依賴性.

2.2 管道流策略

管道流方法的核心思想是在最外層方向上(比如K方向)建造一條管道，在該管道內(nèi)，每個(gè)線程相當(dāng)于一個(gè)流水線.然后將里層(比如J方向)的任務(wù)靜態(tài)地分解成許多子任務(wù)，并分配給相應(yīng)的線程.

圖2 超平面算法示意圖Fig.2 Hyperplane algorithm schematic

圖3 管道流算法示意圖Fig.3 Pipelining algorithm schematic

圖3給出了擁有n條線程(T0，T1，T2，…，Tn-1)或稱管道流線的管道流算法示意圖.在K方向上建造管道，使每個(gè)線程都將K的取值遍歷一遍.對(duì)于第N層(K=N)，J方向上的任務(wù)被分成M份，并依次分配給n個(gè)不同的線程.因?yàn)榇嬖跀?shù)據(jù)依賴，任務(wù)隊(duì)列中任務(wù)KNJM的執(zhí)行需要在任務(wù)KN-1JM和KNJM-1之后.可通過恰當(dāng)?shù)卦O(shè)置標(biāo)志數(shù)組實(shí)現(xiàn)相應(yīng)的線程同步，從而解決數(shù)據(jù)依賴問題.

3 算 例

3.1 并行測(cè)試

為了驗(yàn)證發(fā)展的針對(duì)隱式LU-SSOR方法的兩種并行策略，選取展弦比為13.71的UH-1H槳葉進(jìn)行測(cè)試.槳葉網(wǎng)格的大小為149(弦向)×41(法向)×81(展向)，背景網(wǎng)格的大小為241(周向)×81(縱向)×101(徑向).計(jì)算在Intel(R)Xeon(R)CPU E5-2650 v3 @2.30GHZ服務(wù)器上進(jìn)行，該服務(wù)器共有40個(gè)線程.

以一個(gè)掃描周期為例，利用超平面(Hyperplane)和管道流(Pipelining)兩種方法對(duì)原有串行程序進(jìn)行并行化處理，得到一個(gè)掃描周期內(nèi)的耗時(shí)t(單位： s)隨所用線程個(gè)數(shù)n(單位： 個(gè))變化的關(guān)系圖，如圖4(見第592頁(yè)).可以看出，當(dāng)線程個(gè)數(shù)達(dá)到18時(shí)，并行加速比達(dá)到3.13，說明本文發(fā)展的兩種并行策略使流場(chǎng)計(jì)算的效率得到大幅度提高,因此直升機(jī)旋翼噪聲計(jì)算的周期大大縮減.

3.2 Caradonna-Tung旋翼有升力懸停

圖4 超平面和管道流并行效果Fig.4 Parallel effects of hyperplane and pipelining

Caradonna-Tung旋翼[10]由兩片槳葉組成，槳葉形狀為矩形，展弦比為6，翼型為NACA0012，槳葉半徑為1.143m，弦長(zhǎng)19.05cm，翼型無扭轉(zhuǎn)，無尖削.本算例的計(jì)算狀態(tài)為： 槳尖馬赫數(shù)為0.439，總距角為8°.

圖5 背景網(wǎng)格(左)和槳葉網(wǎng)格(右)Fig.5 Background grid(left) and rotor mesh(right)

計(jì)算網(wǎng)格采用兩塊重疊網(wǎng)格，背景網(wǎng)格大小為221(周向)×121(縱向)×121(徑向)，槳葉網(wǎng)格大小為121(弦向)×81(法向)×81(展向)，如圖5所示.背景網(wǎng)格在周向均勻分布，在垂直方向于槳盤平面附近進(jìn)行加密，徑向在槳尖附近進(jìn)行加密.旋翼槳葉網(wǎng)格的生成采用了雙曲型網(wǎng)格生成方法[11]，只需要給定初始表面網(wǎng)格分布即可自動(dòng)生成體網(wǎng)格.為了提高重疊網(wǎng)格的計(jì)算效率，避免旋翼槳葉網(wǎng)格在計(jì)算過程中出現(xiàn)重疊，外邊界的大小沿展向從小到大分布.

本算例并行計(jì)算共耗時(shí)12.5h,而串行計(jì)算耗時(shí)39.9h,計(jì)算速度提高了3.2倍.圖6給出了計(jì)算得到的不同展向位置翼型壓力系數(shù)(CP)分布與試驗(yàn)值的對(duì)比，計(jì)算結(jié)果與試驗(yàn)值吻合很好，說明了本文發(fā)展的并行計(jì)算方法的準(zhǔn)確性.

圖6 不同展向位置翼型壓力系數(shù)分布Fig.6 Distribution of pressure coefficients at different span positions

對(duì)該算例進(jìn)行網(wǎng)格無關(guān)性分析，選取不同的網(wǎng)格數(shù)量級(jí)進(jìn)行計(jì)算，如表1.圖7給出了不同網(wǎng)格數(shù)量級(jí)計(jì)算得到的懸停拉力系數(shù)與NASA理論計(jì)算數(shù)值的對(duì)比，橫坐標(biāo)為網(wǎng)格數(shù)量n(萬(wàn))，縱坐標(biāo)為拉力系數(shù)CT.由圖7可知，該算例當(dāng)旋翼網(wǎng)格數(shù)量達(dá)到60萬(wàn)，背景網(wǎng)格數(shù)量達(dá)到200萬(wàn)時(shí)，計(jì)算收斂.

圖7 懸停拉力系數(shù)與理論值對(duì)比Fig.7 Comparison of hover thrust coefficient with theoretical value

表1 旋翼網(wǎng)格和背景網(wǎng)格分布

注： 括號(hào)內(nèi)的數(shù)值為網(wǎng)格的數(shù)量.

3.3 AH-1/OLS有升力前飛

AH-1/OLS模型旋翼[12-13]由兩片槳葉組成，槳葉平面形狀為矩形，槳葉帶8.2°的線性負(fù)扭轉(zhuǎn)，展弦比為9.22，靜止?fàn)顟B(tài)下，在距離旋轉(zhuǎn)中心75%展長(zhǎng)處扭轉(zhuǎn)角為0°，翼型為OLS翼型[14-15].槳葉在前飛過程中不僅有旋轉(zhuǎn)運(yùn)動(dòng)，還包括揮舞、變距、擺振等運(yùn)動(dòng).若定義旋翼槳葉的方位角為Ψ，變距角為θ，揮舞角為β，擺振角為δ，則槳葉的運(yùn)動(dòng)規(guī)律為：

ψ=ψ0+ωt,
β=β0+βlccos(ψ)+βlssin(ψ),
θ=θ0+θlccos(ψ)+θlssin(ψ),
δ=δ0+δlccos(ψ)+δlssin(ψ).

對(duì)AH-1/OLS旋翼的10014工況進(jìn)行遠(yuǎn)場(chǎng)距離相關(guān)性分析，對(duì)應(yīng)的計(jì)算狀態(tài)為：

Mtip=0.664,μ=0.164,αTPP=1°,Re=1.6×106,CT=0.0054.

變距和揮舞運(yùn)動(dòng)規(guī)律如下：

θ(t)=6.14°+0.9°cos(ψ)-1.39°sin(ψ),
β(t)=0.5°+1.0°cos(ψ).

計(jì)算所用背景網(wǎng)格為圓柱形網(wǎng)格，在計(jì)算過程中也隨槳葉做旋轉(zhuǎn)運(yùn)動(dòng).本文選取不同的遠(yuǎn)場(chǎng)距離進(jìn)行了計(jì)算比較，如圖8所示，r表示背景網(wǎng)格的半徑，d表示背景網(wǎng)格的縱向長(zhǎng)度.表2給出了4種計(jì)算狀態(tài)和不同遠(yuǎn)場(chǎng)距離的拉力系數(shù)與試驗(yàn)值的誤差，其中R為旋翼槳盤半徑.由表可知隨著遠(yuǎn)場(chǎng)距離的增大，計(jì)算收斂.

圖8 遠(yuǎn)場(chǎng)邊界示意圖Fig.8 Far-field boundary schematic

表2 遠(yuǎn)場(chǎng)距離選擇(a)和拉力系數(shù)的計(jì)算值與誤差(b)

為了將數(shù)值計(jì)算的噪聲與試驗(yàn)值作對(duì)比，選取該旋翼有噪聲試驗(yàn)數(shù)據(jù)的3017工況進(jìn)行噪聲計(jì)算.3017工況的計(jì)算狀態(tài)為：Mtip=0.663,μ=0.162,αTPP=1°,Re=1.6×106.變距運(yùn)動(dòng)規(guī)律為：

θ(t)=5.03°+1.98°cos(ψ)-2.04°sin(ψ).

選取旋翼控制面計(jì)算噪聲.圖9(見第594頁(yè))給出了旋翼控制面上的壓強(qiáng)分布.圖10(見第594頁(yè))給出了計(jì)算得到的噪聲與試驗(yàn)值的對(duì)比，橫坐標(biāo)為方位角ψ，縱坐標(biāo)為壓力脈動(dòng)p′，單位為Pa.計(jì)算結(jié)果與實(shí)驗(yàn)值吻合較好，說明本文發(fā)展的旋翼氣動(dòng)噪聲數(shù)值預(yù)測(cè)方法的可靠性.在該工況下，串行計(jì)算耗時(shí)112.94h，并行計(jì)算耗時(shí)37.8h，計(jì)算速度提高了2.99倍.

圖9 旋翼控制面上的壓強(qiáng)分布示意圖Fig.9 Diagram of pressure distribution on rotor control surface

圖10 計(jì)算噪聲與試驗(yàn)值對(duì)比Fig.10 Comparison of calculated noise with test data

4 結(jié) 語(yǔ)

本文基于CFD/FW-H方法對(duì)旋翼的氣動(dòng)噪聲進(jìn)行了數(shù)值預(yù)測(cè)并基于OpenMP對(duì)源程序進(jìn)行了并行優(yōu)化.針對(duì)隱式LU-SSOR算法存在數(shù)據(jù)依賴性難以實(shí)現(xiàn)并行的問題，發(fā)展了超平面(Hyperplane)和管道流(Pipelining)兩種并行策略，實(shí)現(xiàn)了LU-SSOR算法的OpenMP并行計(jì)算.算例顯示，當(dāng)并行計(jì)算的線程達(dá)到18時(shí)，兩種并行策略使得旋翼渦流場(chǎng)的計(jì)算速度提高3倍以上.因此，旋翼噪聲的計(jì)算周期大大縮短，說明本文發(fā)展的兩種并行策略具有重要的工程應(yīng)用價(jià)值.