基于時域氣聲伴隨算法的直升機(jī)轉(zhuǎn)子葉片形狀優(yōu)化分析

何 磊，魏海軍

(1.安徽農(nóng)業(yè)大學(xué) 經(jīng)濟(jì)技術(shù)學(xué)院，安徽 合肥 230011；2.上海海事大學(xué) 商船學(xué)院，上海 201306)

直升機(jī)機(jī)動靈活，可以在低空飛行，戰(zhàn)略打擊，空中救援，航拍偵查，被廣泛應(yīng)用于軍事領(lǐng)域以及民用航空。旋翼作為直升機(jī)的核心動力部件，研究提升其性能對直升機(jī)的動力性能、飛行性能、飛行穩(wěn)定性有著重要作用。

伴隨方程已被廣泛用于計算流體力學(xué)的氣動外形優(yōu)化問題，是一個非常強(qiáng)大的工具，該方法將流場控制方程作為等式約束引入目標(biāo)函數(shù)，將物體邊界作為控制函數(shù)，創(chuàng)新性地把優(yōu)化設(shè)計問題轉(zhuǎn)化為控制問題，這為復(fù)雜氣動外形優(yōu)化設(shè)計提供了可能，因?yàn)樗试S以基本參數(shù)輸入，來計算一組給定輸入的目標(biāo)函數(shù)的靈敏度導(dǎo)數(shù)。單次梯度信息求解計算量只相當(dāng)于兩倍的流場計算量，與設(shè)計變量的數(shù)目無關(guān)。基于伴隨的穩(wěn)態(tài)問題形狀優(yōu)化被廣泛應(yīng)用于復(fù)雜氣動外形優(yōu)化設(shè)計的二維和三維問題[1]。Mavriplis使用非定常伴隨方法，開展了空氣動力學(xué)中三維問題的應(yīng)用研究[2]。Nielsen等人采用離散型伴隨方法，使用ANSYS CFX開展直升機(jī)旋翼的外形氣動優(yōu)化設(shè)計的研究[3]。李亞等研究了在噪聲計算中采用了FW-H方程、結(jié)合旋轉(zhuǎn)偶極子輻射、渦聲方程三種方法，實(shí)現(xiàn)了大渦模擬預(yù)報螺旋槳輻射噪聲的三種聲學(xué)方法[4]。南京航空航天大學(xué)的吳琪建立了一套適合于先進(jìn)旋翼氣動外形綜合優(yōu)化設(shè)計分析的黏性伴隨方法和代理模型方法相結(jié)合的混合設(shè)計方法，實(shí)現(xiàn)了基于黏性伴隨方法的旋翼先進(jìn)氣動外形優(yōu)化設(shè)計分析[5]。伴隨方法在多學(xué)科設(shè)計優(yōu)化問題的應(yīng)用主要集中在固定翼和旋轉(zhuǎn)翼飛機(jī)的流固耦合上[6]。本文將三維多學(xué)科伴隨方法推廣到旋翼機(jī)的降噪問題。由于對旋翼機(jī)噪聲的要求越來越嚴(yán)格，旋翼氣動聲學(xué)在直升機(jī)的設(shè)計過程中起著主導(dǎo)作用。通常，轉(zhuǎn)子噪聲已通過選擇轉(zhuǎn)子葉片配置、翼型、平面形狀和葉尖形狀來解決。雖然這些方法主要依靠風(fēng)洞實(shí)驗(yàn)和飛機(jī)設(shè)計工程師的經(jīng)驗(yàn)，但很少有人嘗試在降噪背景下設(shè)計優(yōu)化翼型的形狀。在之前開展的項(xiàng)目研究中，我們推導(dǎo)了二維翼型的耦合伴隨氣動聲學(xué)問題，并證明了非定常氣動聲學(xué)優(yōu)化的可行性，這里我們研究了一種三維非定常伴隨方法，使用基于梯度優(yōu)化用于直升機(jī)轉(zhuǎn)子氣動聲學(xué)應(yīng)用的方法，以便通過應(yīng)用最佳形狀修正，在不產(chǎn)生顯著性能損失的情況下使得直升機(jī)轉(zhuǎn)子噪聲最小化。

1 時域氣聲伴隨算法的計算模型的構(gòu)建

1.1 基于CFD方法的直升機(jī)轉(zhuǎn)子葉片流場計算模型

關(guān)于流體的數(shù)值模擬分析公式，本文中使用的是基本三維非結(jié)構(gòu)網(wǎng)格雷諾平均Navier-Stokes流體數(shù)值模擬方法。ANSYS CFX已被廣泛用于驗(yàn)證穩(wěn)態(tài)和時變流，并包含離散切線和伴隨靈敏度能力，這已被證明用于穩(wěn)態(tài)和時變流問題的優(yōu)化。因此，本文將僅對這些公式進(jìn)行簡要描述，其他細(xì)節(jié)可在以前的參考文獻(xiàn)中獲得[7]。流體的數(shù)值模擬基于Navier-Stokes方程的保守形式，可寫成：

(1)

對于移動網(wǎng)格問題，這些問題以任意拉格朗日-歐拉形式表示為：

(2)

使用以頂點(diǎn)為中心的中值雙控制體積公式，該公式為二階精度，其中閉合控制體積周圍的無粘通量積分S離散為：

(3)

(4)

以簡化的形式展示了不同時間水平面上U和x的函數(shù)關(guān)系，如下所示：

Rn(Un,Un-1,Un-2,xn,xn-1,xn-2)=0

(5)

在每個時間步長n，隱式殘差相對于未知解向量Un線性化，并使用牛頓法求解如下：

(6)

其中，向量U表示所有時間步的流量值，并且該方程中的每個(塊)行對應(yīng)于方程(5)中給出的特定時間步的解。方程(6)表示所有時間步的聯(lián)立解，在實(shí)踐中，由于每個新的時間步取決于前兩個時間水平，因此通過使用前向塊替換(即時間上的前向積分)的牛頓法求解。

表1 直升機(jī)葉片轉(zhuǎn)子參數(shù)

1.2 推導(dǎo)氣動靈敏度分析公式

基本敏感性分析的實(shí)施遵循參考文獻(xiàn)中開發(fā)的策略[7]。考慮一個任意的目標(biāo)函數(shù)L，該函數(shù)使用非定常流動解集U和非定常網(wǎng)格解集x進(jìn)行計算，表達(dá)式為：

L=L(U,x)

(7)

假設(shè)狀態(tài)變量(即U，x)依賴于一些輸入設(shè)計參數(shù)D，目標(biāo)函數(shù)L對設(shè)計輸入的總靈敏度可以表示為設(shè)計輸入的狀態(tài)靈敏度向量和狀態(tài)變量的目標(biāo)靈敏度向量之間的內(nèi)積，如下所示：

(8)

如前所述的非線性流動殘差算子和線性彈性網(wǎng)格殘差算子提供了約束，這些約束可以在整個空間和時間域上以一般形式表示為：

(9)

當(dāng)對設(shè)計輸入進(jìn)行線性化時，推導(dǎo)得到：

(10)

這些構(gòu)成了正向靈敏度或切線靈敏度方程。然后，網(wǎng)格和流量靈敏度矢量可以代入方程(8)，以獲得目標(biāo)相對于設(shè)計變量D的完全靈敏度。

對于每個設(shè)計參數(shù)D，正向靈敏度方法需要方程(10)的新解。另一方面，伴隨方法可以以近似獨(dú)立于設(shè)計變量數(shù)目的代價獲得任意數(shù)目的設(shè)計輸入D的靈敏度。伴隨問題可通過將方程(10)預(yù)乘大耦合矩陣的逆，并將靈敏度的結(jié)果表達(dá)式代入方程(8)，并將伴隨變量定義為系統(tǒng)的解來獲得：

(11)

其中ΛU和Λx分別是流動和網(wǎng)格伴隨變量。最終目標(biāo)靈敏度可通過以下方式獲得：

(12)

回顧方程(11)和(12)適用于整個時域，方程(13)表示正向時間積分，而方程(12)中的回代過程導(dǎo)致時間上的反向積分，從最后一個物理時間步，然后進(jìn)入初始時間步。

1.3 前飛旋翼槳距公式

由于本文的目標(biāo)是優(yōu)化配平前飛狀態(tài)的旋翼，因此必須將旋翼槳葉槳距角作為控制輸入，以便旋翼能夠獲得所需的推力和側(cè)向力矩[8]。俯仰角可以表示為平均俯仰角和幾個諧波分量的組合：

(13)

其中ψ是轉(zhuǎn)子方位角。在本文中，我們只使用一個諧波，因此i=1，設(shè)計變量D=[θ0，θci，θsi]的選擇應(yīng)使轉(zhuǎn)子得到調(diào)整，其中θ0是總參數(shù)，θci，θsi是兩個周期螺距參數(shù)。為了將旋翼槳距驅(qū)動納入計算流體動力學(xué)分析中，需要一個葉片表面槳距驅(qū)動方程(剩余形式)：

Sθ=xSθ+[Tθ(θ(D))]x=0

(14)

和網(wǎng)格運(yùn)動方程

G′(xθ,xsθ)=0

(15)

其中，xsθ表示俯仰后的葉片表面網(wǎng)格坐標(biāo)，xθ表示俯仰和方位旋轉(zhuǎn)后網(wǎng)格運(yùn)動的內(nèi)部網(wǎng)格坐標(biāo)，xp表示俯仰和方位旋轉(zhuǎn)后的表面和內(nèi)部網(wǎng)格坐標(biāo)，而U表示的是非定常流動解。

2 轉(zhuǎn)子葉片噪聲優(yōu)化分析

在本文中建立的模型是前飛狀態(tài)的 NACA0012 四葉旋翼，其自由流馬赫數(shù)為0.095，翼尖馬赫數(shù)為0.638，雷諾數(shù)為2×106。在此條件下，非線性流動效應(yīng)可以忽略不計，因此可以使用機(jī)身和機(jī)翼組成的不可滲透的集成表面。聲學(xué)表面與ANSYS CFX轉(zhuǎn)子表面網(wǎng)格重合，表面上的每個節(jié)點(diǎn)都是一個聲源，通過 FW-H 集成過程在觀察者位置產(chǎn)生聲壓[9-10]。將聲學(xué)模塊的線性化與流動求解器的線性化相結(jié)合，可以計算聲學(xué)目標(biāo)函數(shù)相對于在基于梯度的形狀優(yōu)化過程中使用的設(shè)計變量的完整向量的伴隨靈敏度。

2.1 聲學(xué)分析公式：FW-H方程

FW-H方程可以用微分形式表示為：

(16)

其中，

其中，公式(17)表示流體氣動彈性非穩(wěn)定力關(guān)系函數(shù)，公式(18)表示流體質(zhì)量變化函數(shù)。f表示所定義的求解積分面，是時間和位移的函數(shù)。ρ表示空氣總密度，p表示總壓力。ui為流體速度，而vi為轉(zhuǎn)子葉面速度，c0為自由流體產(chǎn)生的聲速。ni是聲學(xué)表面的單位法線。然后通過方程(22)給出積分解，忽略四階無窮小量。

其中y是觀察者位置，t是觀察者時間，r是源x和觀察者y之間的距離，M是表面馬赫數(shù)。方程(22)是基于時域氣聲伴隨算法的FW-H方程，p′T(y,t)代表由厚度產(chǎn)生的聲壓在葉面面壁邊界上的積分，p′L(y,t)代表由載荷產(chǎn)生的聲壓在葉面面壁邊界上的積分，其中含1/r的項(xiàng)為遠(yuǎn)場噪聲，含1/r2的項(xiàng)為近場噪聲。

表2 聲學(xué)目標(biāo)函數(shù)相對于直升機(jī)轉(zhuǎn)子葉片的觀察者位置(R是轉(zhuǎn)子葉片半徑)

2.2 基于時域氣聲伴隨算法的FW-H靈敏度公式

計算耦合氣動聲學(xué)問題的靈敏度需要將聲學(xué)模型相對于設(shè)計變量進(jìn)行線性化[11]。如前一節(jié)所述FW-H積分方程的靈敏度公式與其他參考文獻(xiàn)中詳述的CFD數(shù)值模擬曲線的線性化非常相似。首先，聲碼的正向線性化是通過離散FW-H積分方程的精確微分設(shè)計的。然后轉(zhuǎn)置正向線性化，并從相反的方向推導(dǎo)求解得到關(guān)于設(shè)計變量完整向量的伴隨靈敏度公式。針對聲碼的復(fù)雜微分形式使用正向靈敏度進(jìn)行驗(yàn)證。最后，通過伴隨靈敏度公式來驗(yàn)證切線的雙重性問題。

聲學(xué)目標(biāo)函數(shù)的切線靈敏度公式，是計算關(guān)于單個設(shè)計變量的目標(biāo)函數(shù)的正向線性化靈敏度。本文將觀察者位置和時間、聲壓關(guān)系表示為：

p′(y,t,D)=F(U(D),x(D))

(23)

其中D是設(shè)計變量的向量表示，F(xiàn)(U(D)，x(D))表示以數(shù)值方式評估方程 (23) 所需的所有離散運(yùn)算。在這個表達(dá)式中，U(D)和x(D)是時間積分過程中每個時間步的流動和網(wǎng)格解。聲壓切線靈敏度時間方程可表示為：

(24)

所以推導(dǎo)得到定義的聲學(xué)目標(biāo)函數(shù)的切線線性化方程表示為：

(25)

伴隨靈敏度公式，伴隨靈敏度可以通過轉(zhuǎn)換正切靈敏度公式得到。轉(zhuǎn)換方程(25)得到

(26)

3 轉(zhuǎn)子葉片算例驗(yàn)證

相應(yīng)的轉(zhuǎn)子轉(zhuǎn)速為1100轉(zhuǎn)每分鐘(旋翼軸向進(jìn)給比μ=0.15)。ANSYS CFX數(shù)值模擬由2.32×106個節(jié)點(diǎn)所組成的網(wǎng)格，進(jìn)行了兩次旋翼葉片旋轉(zhuǎn)試驗(yàn)，網(wǎng)格由2度時間步長、每個時間步長6個流固耦合迭代和位于距轉(zhuǎn)子輪轂兩個半徑的靜止平面內(nèi)觀測器組成，方位角ψ= 180度。進(jìn)行了優(yōu)化設(shè)計計算。如圖1所示的結(jié)果，在這兩種情況下，對比發(fā)現(xiàn)本文所建立的FW-H模型實(shí)現(xiàn)與實(shí)驗(yàn)值之間一致性非常好[11]。表明本文計算方法和建?？梢杂行У仡A(yù)測前飛狀態(tài)的旋翼轉(zhuǎn)子葉片的噪聲。

圖1 某型直升機(jī)葉片轉(zhuǎn)子結(jié)構(gòu)圖

基于時域氣聲伴隨算法的直升機(jī)轉(zhuǎn)子葉片F(xiàn)W-H積分方程的靈敏度優(yōu)化收斂結(jié)果如圖2所示，經(jīng)過12次非線性迭代和20次設(shè)計循環(huán)之后，優(yōu)化問題的可行性和最優(yōu)性伴隨靈敏度提高了3個數(shù)量級。對直升機(jī)四葉NACA0012旋翼進(jìn)行配平之后，其推力和側(cè)向力矩關(guān)于時間歷程的數(shù)值變化進(jìn)一步表明了本文計算方法和建?？梢杂行У仡A(yù)測前飛狀態(tài)的旋翼轉(zhuǎn)子葉片的噪聲。圖3是直升機(jī)轉(zhuǎn)子旋翼葉片的厚度噪聲、載荷和總聲壓隨時間的變化曲線，厚度噪聲僅受到主要作用于加載壓力的優(yōu)化的輕微影響，從而得到最小化的總聲壓。噪音最小化是通過改變厚的翼型實(shí)現(xiàn)的，特別是在內(nèi)側(cè)部分。由于轉(zhuǎn)子扭矩不包括在優(yōu)化范圍內(nèi)，因此盡管代表了當(dāng)前數(shù)值模擬的優(yōu)化方法可實(shí)現(xiàn)最大的降噪量，但這里實(shí)現(xiàn)的降噪是以犧牲直升機(jī)旋翼轉(zhuǎn)子性能為代價的，如圖4所示。最后得到了氣動聲學(xué)優(yōu)化直升機(jī)轉(zhuǎn)子旋翼葉片旋翼的優(yōu)化翼型形狀，如圖5所示。

圖2 直升機(jī)轉(zhuǎn)子旋翼計算噪聲聲壓與實(shí)驗(yàn)值對比曲線

圖3 直升機(jī)轉(zhuǎn)子旋翼配平非線性迭代函數(shù)優(yōu)化的可行解和最優(yōu)解對比曲線

圖4 直升機(jī)轉(zhuǎn)子旋翼葉片的厚度噪聲、載荷和總聲壓隨時間的變化曲線

圖5 氣動聲學(xué)優(yōu)化直升機(jī)轉(zhuǎn)子旋翼葉片旋翼的優(yōu)化翼型形狀

4 結(jié)論

本文建立了一種直升機(jī)氣動聲學(xué)時變問題的靈敏度分析公式。采用ANSYS CFX數(shù)值模擬方法與FW-H聲學(xué)模塊耦合的混合氣動聲學(xué)方法，將直升機(jī)轉(zhuǎn)子噪聲信號傳播給遠(yuǎn)場觀測者。得到以下結(jié)論：

(1)推導(dǎo)了與氣動聲學(xué)分析問題類似的切向靈敏度公式和伴隨靈敏度公式。首次驗(yàn)證了全耦合氣動聲學(xué)分析公式能有效地預(yù)測四葉NACA0012旋翼在前飛時的噪聲。在驗(yàn)證成功的基礎(chǔ)上，采用伴隨公式對旋翼噪聲信號進(jìn)行優(yōu)化。

(2)前飛狀態(tài)的直升機(jī)旋翼旋轉(zhuǎn)總噪聲、厚度噪聲、載荷和總聲壓優(yōu)化后數(shù)值有明顯降低的趨勢，厚度噪聲僅受到主載荷壓力的優(yōu)化輕微影響，從而得到最小化的總聲壓。噪音最小化是通過改變厚的翼型實(shí)現(xiàn)的，特別是在內(nèi)側(cè)部分。

(3)本方法對其他汽車噪聲降噪處理、航空固定翼機(jī)型機(jī)翼氣動聲學(xué)優(yōu)化都有較強(qiáng)的參考意義。