佟勝喜, 李東輝, 趙慶賀, 高 峰, 范彥銘
(1.沈陽飛機設計研究所, 沈陽 110035; 2.遼寧通用航空研究院, 沈陽 110136; 3.沈陽航空航天大學航空宇航學院, 沈陽 110136)
為了降低航空運輸帶來的碳排放,采用新能源代替燃油的綠色航空成為航空界的研究熱點之一。波音公司、空中客車公司,以及美國航空航天局(NASA)、德國航空航天研究院(DLR)等研究機構均開展了關于電動飛機技術的相關研究。中國在電動飛機領域也開展了相關的研發工作。中國商飛公司等單位共同發起研制了“靈雀H”氫燃料電混合動力飛機,中國航空研究院聯合荷蘭宇航院積極開展電動飛機領域的國際合作,遼寧通用航空研究院開展了系列電動飛機的研制工作,其中RX1E雙座電動飛機、RX1E-A增程型雙座電動飛機已經取得適航證并進入市場,目前正在進行四座電動飛機和水上飛機的研制工作。
電動飛機使用電能作為能量來源,主要有太陽能飛機、燃料電池飛機和蓄電池飛機等。現階段電動飛機主要以鋰離子電池組作為能源儲備裝置,受鋰離子電池水平的限制,電動飛機較同等量級的燃油飛機的航時、航程等都偏低。為了增加電動飛機的航時,需要盡可能地提高飛機的巡航升阻比等氣動指標[1-2]。因此,電動飛機與常規飛機相比,在氣動布局方面需采用更大的展弦比。
飛機是一個彈性體,在氣動載荷的作用下,機翼發生明顯變形,當氣動載荷為定常載荷時,結構的彈性變形是一個緩慢的過程,氣動力和運動與時間無關,靜氣動彈性問題[3-7]就是研究飛行器在定常氣動載荷作用下的變形問題。氣動力和結構變形的相互耦合,氣動彈性變形影響著飛行條件下的載荷重新分布、升力分布、阻力、操縱效率、飛機的配平以及靜穩定性和操縱性等。
靜氣動彈性問題假設結構的彈性變形是一個十分緩慢的過程,如果機翼的結構設計不當,例如機翼的剛心與壓力中心距離過大,當飛行器在較大的速壓下飛行時,作用在機翼上的扭矩過大直至超過機翼的強度極限,導致機翼結構破壞。
當飛行速度比較低時,彈性變形的影響很小。隨著飛行速度的增加,彈性變形的影響也越來越嚴大,甚至會使機翼變得不穩定,或使操縱面失效或反效。萬志強等[8]通過使用基于偶極子網格法及線性氣動力影響系數矩陣的線性靜氣動彈性分析方法對飛機進行了不同飛行高度、不同縱向過載下的靜氣動彈性響應特性分析。鄧立東等[9]考慮了氣動力的非線性問題,開展了飛機非線性飛行載荷計算方法研究。謝長川等[10]介紹了大展弦比大柔性飛機的氣動彈性研究成果,著重提出了具有結構幾何非線性的飛機的氣動彈性和飛行動力學問題。嚴德等[11]開展了基于試驗氣動力的縱向機動飛行載荷分析,氣動彈性對大展弦比機翼升力的展向分布的影響是飛機結構設計者考慮的一個重要因素。楊希明等[12]闡述了中外在氣動彈性風洞試驗開展的主要研究工作,指出了氣動彈性風洞試驗在飛行器研制中的重要意義。
遼寧通用航空研究院設計的RX1E雙座電動飛機采用了較大的展弦比,其展弦比達到17.52,大展弦比飛機機翼的彈性變形相對較大,因此其受到氣動力作用后的彈性變形也會更加明顯。
現通過面元法[13-14]和計算流體力學/計算固體力學(computational fluid dynamics/computational structural dynamics,CFD/CSD)耦合方法得到銳翔電動飛機機翼在風洞試驗條件下的彈性變形,并分析彈性變形對縱向氣動特性的影響。
開展雙座電動飛機風洞試驗時,試驗模型為 1∶6 全金屬模型,翼展2.417 m,模型結構如圖1所示。模型主要連接件使用合金鋼,機身、機翼等部分使用鋁合金加工。試驗風速為70 m/s。對機翼部分單獨建模并進行計算,使用CFD與機翼結構耦合計算機翼的變形,對比分析機翼彈性變形前后氣動力變化和載荷變化情況。

圖1 機翼風洞試驗模型Fig.1 The wind tunnel test model
相比于CFD,面元法可以快速高效地獲得氣動載荷,高階面元法可實現對復雜外形飛行器三維氣動力建模,相比于平板面元法能更好地模擬真實飛行器表面的流動情況。
首先采用Patran[15]劃分機翼有限元網格,圖2給出了機翼有限元網格以及采用獲得的彈性變形。
為了縮短時域方法計算周期,減少存儲量,結構計算采用簡化的線性模態疊加方法,通過在物面施加近似邊界條件,反復迭代求解得到結構體的實際變形。結構動力方程為

(1)

氣動彈性分析需要對氣動模型與結構模型進行耦合,使氣動數據與結構數據相互傳遞。而氣動彈性問題一般需要兩個變換:①將結構模型的位移轉換到氣動模型的位移;②將氣動模型的氣動力轉換成結構模型中的等效力。氣動模型與結構模型耦合示意圖如圖3所示。

圖2 機翼機構網格及位移變形Fig.2 Wing structure mesh and displacement deformation

圖3 結構網格和氣動網格數據傳遞Fig.3 Structural grid and aerodynamic grid data transfer
位移變換需要一個樣條矩陣,公式為
uk=Gkgug
(2)
式(2)中:Gkg為樣條矩陣;ug為結構模型中的位移;uk為氣動模型中的位移。通過Gkg將結構模型中的位移轉換成氣動模型中的位移。
氣動力等效通過虛功原理完成,空氣動力Fk及作用在結構網格點上的等效值Fg再虛位移上做的功相同,公式為

(3)
式(3)中:δuk與δug分別表示氣動力及結構點的虛位移。
將式(2)代入到式(3)中得到

(4)
由于虛位移的任意性,可得到氣動力Fk與Fg的關系為

(5)
為對比驗證面元法獲得的機翼彈性變形,采用了CFD/CSD流固耦合方法[16-20]進行計算。耦合迭代的過程中需要在氣動模型和結構模型間進行數據交換,即將氣動力插值到結構網格上,而把結構位移插值到氣動網格上。另外,結構位移使得機翼物面變形,氣動網格也要隨之變形,因此,還需要在每一步迭代中實現網格變形。
采用基于非結構背景網格的動網格方法來實現CFD網格的快速更新。該方法結合了彈簧網格方法與代數背景網格方法的優點,變形能力強、計算效率高且魯棒性好。
氣動網格如圖4所示,分別給出了機翼網格的整體視圖及局部放大圖。物面附近為六面體網格,過渡區域為金字塔網格,網格規模為5.32×106。

圖4 機翼氣動網格Fig.4 The wing grid partial view
采用有限體積法求解RANS方程,三維可壓縮RANS方程可以寫為

(6)
式(6)中:w是守恒變量;F、FV分別為無黏和黏性通量。

(7)

空間離散采用二階中心格式,時間推進采用LUSGS隱式時間推進。湍流模型采用kω-SST兩方程湍流模型。
在中國航空工業空氣動力研究院FL-8風洞中進行全機縮比模型的風洞試驗,圖5給出了風洞試驗模型在風洞中安裝的照片。使用翼身組合體與光機身的試驗結果差值得到機翼的氣動力試驗結果。試驗中對支架干擾、洞壁干擾進行了修正。

圖5 FL-8風洞中的試驗模型Fig.5 Test model installed in FL-8 wind tunnel
靜氣動彈性問題研究飛行器彈性變形對定常氣動載荷分布的影響以及由氣動載荷所產生的靜變形的穩定性。結構彈性變形將引起氣動載荷的重新分布,氣動載荷的重新分布不僅導致內部結構載荷和應力的重新分布,也會改變空氣動力穩定性導數。
數值模擬的來流速度為70 m/s,質心位于根弦25%位置處,剛性軸為Z軸。計算時與風洞試驗相同的條件進行計算。

圖6 機翼變形前后的氣動特性Fig.6 Aerodynamic characteristics of wings before and after deformation
縱向氣動力數據對比結果如表1所示,表中CL0為零迎角升力系數,CLα為升力線斜率,CD為阻力系數,α為迎角,Cm為俯仰力矩系數,dCm/dCL為俯仰力矩系數Cm對升力系數CL的導數,表征縱向靜安定裕度。圖6給出了面元法、CFD與風洞試驗結果的對比曲線,參考點均位于25%平均氣動弦長。通過數據和曲線可以看出,對于升力系數,CFD結果的升力線斜率略小于試驗結果,面元法在線性段的值較試驗結果整體約小0.04。受彈性變形影響,升力系數增加,阻力系數減小,相同迎角下的升阻比提高,相同升力系數下的升阻比基本沒有變化。升阻比的CFD結果較試驗偏小的原因主要是CFD的阻力系數偏大約0.005。
面元法與CFD計算的縱向靜安定裕度遠小于風洞試驗結果。與剛性模型相比,彈性變形后的縱向靜安定裕度顯著提高,迎角2°時dCm/dCL值由剛性的-0.046變化至彈性的-0.061,且彈性模型的靜安定裕度隨升力系數的提高逐漸增大,到6°迎角時dCm/dCL≈ -0.08,剛性模型變化很小。彈性模型的靜安定裕度隨升力系數的提高逐漸增大的趨勢與風洞試驗結果相同。

表1 縱向氣動力特性對比
根據CFD結果,分析彈性模型變形后彎矩、扭矩和剪切應力的變化情況。圖7給出了不同迎角下的彎矩、扭矩及剪力,其中GM表示剛性模型,TM表示彈性模型。通過曲線可以看出,變形前后的彎矩及剪力變化不大,扭矩變化較為明顯。

圖7 機翼變形前后的氣動載荷Fig.7 Aerodynamic loads before and after wing deformation
為獲得電動飛機的彈性變形,采用面元法及CFD/CSD耦合方法實現對機翼靜氣動彈性的數值模擬,并開展了靜氣動彈性風洞試驗,結論如下。
(1)采用面元法獲得的外載荷,在氣動力的線性段與CFD/CSD趨勢一致,升力系數約小0.05,計算結果表明該方法能滿足工程實際中計算結果和計算高效的要求。
(2)受彈性變形影響,升力系數增加,阻力系數減小,相同升力系數下的升阻比基本沒有變化。彈性變形對俯仰力矩系數影響顯著,變形后的縱向靜安定裕度顯著提高,dCm/dCL由剛性的-0.046變化至彈性的-0.061,且彈性模型的靜安定裕度隨升力系數的提高逐漸增大,到6°迎角時dCm/dCL≈-0.08,剛性模型變化很小。大展弦比機翼在翼稍附近出現了明顯的彈性變形。
(3)考慮彈性變形影響后的CFD/CSD縱向靜安定性計算結果與風洞試驗表現出了相同的隨迎角增加趨勢。
(4)機翼扭矩受彈性變形影響較大,彎矩和剪力幾乎沒有變化。