風扇葉片飛失顯式動力學仿真網格與時間步長研究

徐 雪 ，李宏新 ，馮國全 ，3

（1.中國航發沈陽發動機研究所，沈陽 110015；2.中國航空發動機集團有限公司，北京 100097；3.遼寧省航空發動機沖擊力學重點試驗室，沈陽 110015）

0 引言

風扇葉片飛失（Fan Blade Out，FBO）后發動機動力學行為及其影響的設計技術是民用航空發動機安全性領域最關鍵的技術壁壘之一，是各適航體系都要求的重要試驗[1-3]。近年來，隨著顯式動力學計算軟件逐步成熟，憑借其對撞擊、碰摩和斷裂等典型高度瞬態、高度非線性物理過程的有效模擬，成為FBO 設計和仿真工作的重要方法[4]。

在該領域，Nicolas 等[5]開展了FBO 過程仿真與整機試驗測試結果的對比分析，提出風扇轉子和靜子的模型簡化方法；Robin[6]進行了整機FBO仿真分析中風扇轉子葉尖嚴重碰摩的瞬態動力學行為計算，描述了葉尖碰摩的微觀物理過程；Sinha[7-8]系統地發展了考慮碰摩的轉子動力學行為計算理論，并給出了典型算例的仿真分析結果以及與整機試驗數據的對比；馬輝等[9]開展了旋轉葉片-機匣系統碰摩動力學研究，提出轉子盤、葉片和機匣碰摩系統的動力學簡化分析方法；劉璐璐等[10]開展了縮比的FBO 整機模擬試驗與仿真研究，獲得安裝結構的外傳振動特點。由于顯式動力學的算法特點，其網格尺度和時間步長設置在分析中具有關鍵作用，不僅影響計算精度和資源的需求，還直接影響計算的收斂性和物理過程時序，因此網格無關性和時間步長無關性研究對于瞬態顯式動力學分析結果具有決定性影響，但是在以往的研究中對于這方面的分析較少。

本文基于商用有限元顯式動力學計算軟件，開展了FBO 過程中風扇葉片的網格與時間步長無關性研究。

1 顯式動力學網格與時間步長的理論關系

1.1 計算方法和主要影響因素

顯式動力學分析的基本動力學方程[11]為

式中：[M](t)為系統質量矩陣；[C+CT+G](t)為阻尼矩陣，其中C為系統阻尼，CT為接觸阻尼，G為陀螺力矩；[K+KT](t)為剛度矩陣，其中K為系統剛度，KT為接觸剛度；{u}(t)為系統自由度向量；{F+FT}(t)為力向量，其中F為系統外力和慣性力，FT為接觸相互作用力；?t為時間增量步長；β、γ為Newmark 法中加速度矢量和速度矢量的時間步參數占比控制常數。

顯式動力學仿真計算所需的計算資源和總運算時間取決于每個增量時間步中方程解算時間以及總的計算步數。方程解算時間一方面取決于矩陣和力向量的數據量大小（即計算模型網格的總自由度數），另一方面取決于上述矩陣中有多少數據需要參與到碰摩、損傷及穿透等復雜的迭代計算過程，例如采用罰函數法進行局部碰摩行為的迭代計算。這2 種因素都與結構的具體設計特點和物理過程有關，需要根據仿真問題進行具體分析。

總的計算步數由仿真目標時長和時間步長?t決定，尤其是?t與網格的尺度有緊密聯系，而且與計算過程的精度和穩定程度息息相關，也需要根據仿真問題進行具體分析。β、γ只能間接地在小范圍內調節時域積分計算的速度和效率。

由上述分析可知，針對給定的仿真結構和物理過程，能夠直接影響計算速度和結果精度的只有網格劃分和?t2 種因素，而且二者之間又有緊密聯系。同時由上述計算方法可知，在顯式計算中沒有前后時間步的迭代收斂計算，因此前序時間步中形成的計算偏差會隨時間步的不斷增加逐步積累放大。

1.2 網格與時間步長的理論關系

在?t的選擇上，由于瞬態顯式動力學分析所關注的物理過程包括碰撞（時長約10-2s 數量級）、沖擊（時長約10-3s 數量級）和爆炸（時長約10-4s 數量級）[13]等，由于其時間歷程很短，可能與計算網格中應力波傳遞時間尺度相近（如圖1 所示），因此必須考慮應力波的影響。

圖1 物理過程結構

根據顯式動力學計算理論[14]，積分步長必須小于臨界時間步長，否則會使1 個計算時間步長內應力波傳播通過多個網格，使結點參數變化劇烈，從而使計算過程不穩定。臨界時間步長取決于計算網格中最小尺度單元的特征長度

材料中縱波的傳播速度為

式中：?tcr為臨界時間步長；Lch為網格特征長度；c為材料中縱波的傳播速度，即材料中的聲速；E為材料壓縮彈性模量；μ為材料泊松比；ρ為材料密度。

該特征長度在不同的網格類型和不同的計算程序中有多種設置方法，以單個六面體實體單元為例，Lch通常有最短邊長、最短對角線以及其他計算方法

式中：lx、ly、lz分別為x、y、z方向的單元邊長。

以六面體單元最短邊長法計算臨界時間步長為例：如果材料 為TC4 鈦合金[15]，在100 ℃時，E=110 GPa，μ=0.34，ρ=4440 kg/m3，算得聲速為6175 m/s，網格最短邊長為10 mm，則其開展瞬態顯式動力學分析的臨界時間步長為1.62×10-6s。

上述研究的臨界時間步長是從計算穩定性角度給出的時間步長的理論上限，小于這個理論步長計算過程也可能順利完成，但是其計算結果可能不符合真實的物理過程，因此還需要根據實際物理過程進一步計算分析。

2 FBO的物理過程及有限元建模

2.1 FBO的物理過程

FBO事件主要有以下幾個典型的瞬態物理過程：

（1）飛失的風扇葉片與機匣和后續葉片的碰撞損傷過程；

（2）飛失后風扇轉子不平衡量顯著增大，低壓轉子在不平衡離心力作用下變形破壞過程；

（3）結構載荷不斷增大，觸發結構保險破壞失效[4]過程；

（4）剩余低壓轉子穩定在風車轉速持續運轉過程[16]。

前3 個物理過程不是孤立的，而是綜合體現在FBO 后的很短時間內（轉子旋轉幾周，總時長在10-2～10-1s 數量級）不平衡轉子與機匣的碰摩、變形和快速減速過程，主要約束條件是適航要求中葉片包容和大振動載荷外傳等條款[1-3]；第4個物理過程實際上也包含剩余葉片與機匣碰摩的物理過程，但是由于持續時間較長，一般采用隱式動力學方法，主要約束條件是適航要求中轉子持續旋轉條款[1-3]。

由于第3個物理過程與結構保險的具體結構形式息息相關，而且在仿真應用中多采用帶有門限值的接觸條件方法來模擬，并且第4 個物理過程關注的時間范圍和采用的分析方法不在本文討論的范圍內，所以本文中只研究前2個物理過程。

2.2 模型簡化、網格劃分及邊界條件設置

根據新一代雙轉子大涵道比渦扇發動機結構設計特點[17]和具體結構參數[18-19]，通過一系列由簡到繁的模型簡化試算分析，逐漸積累對于具體結構設計特點及物理過程的認識，建立了顯式動力學分析簡化模型，進行網格劃分和參數設置，如圖2 所示并見表1。

圖2 仿真模型及網格

表1 風扇葉片參數及仿真目標參數

2.2.1 模型簡化

為了聚焦葉片與風扇機匣（包括易磨耗涂層）的碰摩過程，排除前序葉片對于后續葉片的影響，采用僅保留1 片葉片的風扇轉子，但是葉片的材料、根/尖半徑、不同截面的弦長和厚度以及葉片的扭轉角等主要參數與真實發動機的基本一致；為了準確模擬低壓轉子的動力學行為，低壓轉子的結構形式、尺度、相互連接的拓撲關系以及材料與真實發動機的基本一致；為了準確模擬葉尖與機匣碰摩所產生的接觸剛度和阻尼，以及易磨耗涂層的剝落消耗過程，風扇機匣和易磨耗涂層的結構和材料與真實發動機的基本一致；為了降低仿真結構的復雜性，低壓轉子的3 主軸承都采用彈簧單元來模擬。

2.2.2 網格劃分

為了更加準確地控制網格的尺度，并更好地模擬碰摩損傷過程，全部模型采用六面體網格；在FBO 物理過程中，風扇葉尖與風扇機匣的碰摩是接觸速度最高、破壞作用最明顯、對于結構系統影響也最大的細節過程，因此選擇風扇葉片和易磨耗涂層2 個部件（碰摩對）作為主要關注部件，其網格尺度根據研究需要采用不同的尺度；其余部件網格出于簡化計算的需要，采用較大的尺度，而且各計算項目保持一致。

2.2.3 邊界條件設置

低壓轉子設置轉動邊界條件模擬轉子轉動行為；風扇機匣和轉子軸承的靜子端采用固定邊界條件來約束整個轉靜子模型系統；轉/靜子之間事先不限定是否接觸，由軟件根據計算結果自行判斷是否接觸，接觸后以恒定的摩擦系數方法確定接觸面的邊界條件。

2.3 測量參數及測點布置

針對前2 個物理過程，選擇3 個有代表性的輸出參數，并根據模型的實際情況布置參數測點，如圖3所示。

圖3 葉尖碰摩的測試點位置

為了反映FBO 狀態下，離心力與葉片碰摩產生的附加接觸剛度和阻尼的共同影響，在風扇盤前緣軸心位置設置了軸心軌跡測點，輸出形式為以時間為X軸的空間曲線，衡量標準為不同網格尺度和時間步長的計算項目得到的空間曲線的一致性；為了反映葉片碰摩過程中葉片的動力學響應和損傷情況，在葉片中部沿葉高方向等距均布了5 個應力響應測點，應力取值采用米塞斯等效應力，形式為相同測試位置的應力時序變化，衡量標準為不同網格尺度和時間步長的計算項目得到相同位置應力隨時間變化的趨勢和量值的一致性；為了反映轉子通過軸承外傳的振動載荷，在1 支點軸承水平位置設置了振動加速度測點，形式為同一測試位置的加速度時序變化，衡量標準為不同網格尺度和時間步長的計算項目得到加速度時序的變化趨勢和量值的一致性。

3 基于DOE 方法的參數研究方案及總體完成情況

采用網格尺度和時間步長為輸入參數，通過大范圍的初步試算，縮小參數范圍，最終確定每個參數的4個取值，通過排列組合形成16種不同計算參數的仿真計算項目表。在相同的硬件和軟件平臺下開展了上述計算項目表中所有的仿真計算，完成情況、總計算耗時及系統節點數量見表2。

表2 不同參數的仿真計算項目及計算情況

4 仿真結果分析

4.1 網格無關性檢驗

在仿真分析中，時間步長為3×10-7s 和2×10-7s2 組不同網格尺度的仿真計算全部完成，分組對比分析如下。

4.1.1 網格尺度對轉子軸心軌跡的影響

以時間步長同為3×10-7s的4個仿真計算為例，如圖4 所示。從圖中可見，網格最粗的計算結果A1-B3 線在第0.02 s 左右就與其他3 條曲線出現明顯差別，而其他3條曲線在第0.05～0.06 s和第0.08～0.10 s也出現了較大偏差，從仿真視頻來看，在第0.05～0.06 s 和第0.08～0.10 s 都出現了葉片與機匣的劇烈碰摩。

圖4 時間步長為3×10-7 s的3維軌跡

以時間步長同為3×10-7s 的4 個仿真計算為例開展同樣分析的結論基本一致，只是第0.08 s 之后的虛線變化時序與上述分析略有差別。

4.1.2 網格尺度對葉片應力的影響

以時間步長為2×10-7s 的4 個仿真計算結果為例，在不同的網格尺度下，葉身正中部測點的應力對比如圖5 所示。從圖中可見，網格尺度為40 mm 的計算結果曲線與其他3 條明顯不一致，隨著網格逐漸加密，葉身應力的差異逐漸縮小，盡管其總的趨勢基本一致，但是在第0.08 s 之后具體的應力大小變化的時間順序規律不一致。

圖5 時間步長為2×10-7 s時不同網格尺度計算的葉中位置應力

4.1.3 網格尺度對支點外傳振動的影響

以時間步長為2×10-7s 的4 個仿真計算結果為例，在不同的網格尺度下，1 支點外傳振動對比分析如圖6 所示。從圖中可見，在A-F區域范圍內，各條曲線的偏差較大且沒有規律，在其余時間范圍，幾條曲線符合得較好，但是隨著時間的增加曲線偏差變大。

圖6 在相同時間步長下不同網格尺度對于支點振動的影響

仿真視頻的結果如圖7 所示。從圖中可見，從A～F 的6 個時間范圍內都發生了風扇葉片與機匣的劇烈碰摩。

圖7 第0.034 s的仿真結果

總體來看，對于時間步長為2×10-7s 和3×10-7s 的2 組計算，網格尺度為20 mm 和15 mm 的計算精度已經十分接近，但是后者的計算資源需求卻是前者的近2倍。

4.2 時間步長無關性檢驗

在上述仿真分析中，只有網格尺度為40 mm 的1組不同時間步長的仿真完成了全部計算，分組對比分析如下。

4.2.1 時間步長對于轉子軸心軌跡的影響

對于給定的網格尺度，第0.06 s 之前不同時間步長的軸心軌跡符合得都很好，第0.06 s 之后時間步長最大的A1-B1 線和A1-B2 線逐漸出現偏差，第0.09 s之后幾條曲線偏差逐漸增大，且時間步長越大偏差增大越快，如圖8所示。

圖8 網格尺度為40 mm的軸心3維軌跡

4.2.2 時間步長對葉片應力的影響

對于2 個給定的網格，第0.35 s 之前的不同時間步長的計算曲線符合得很好，之后不同時間步長計算的偏差明顯，但是總體的趨勢和量級基本一致，如圖9所示。

圖9 不同時間步長下的葉片中部同一節點的應力

4.2.3 時間步長對支點外傳振動的影響

以網格尺度為40 mm 的4 個仿真計算結果為例，在不同的時間步長下，1 支點外傳振動對比分析如圖10 所示。從圖中可見，在A～F 范圍內，各條曲線的偏差較大且沒有規律，在其余時間范圍內，幾條曲線符合得較好，但是隨著時間的增加曲線偏差變大。

圖10 相同時間步長情況下不同網格尺度對于支點振動的影響

在A～F 的6 個時間范圍與前述網格無關性分析中的時間范圍基本一致，因此分析各條曲線偏差較大且沒有規律的計算結果也是因風扇葉片與機匣的劇烈碰摩造成的。

總體來看，對于網格尺度40mm 的計算，時間步長為3×10-7s 和2×10-7s 的計算符合性已經比較好，但是后者的計算資源需求卻是前者的近1.5倍。

5 結論

（1）顯式動力學的計算效率和結果精度與網格和時間步長緊密相關，所以在任何工程分析之前都應該針對具體的物理過程和結構設計特點，開展網格和時間步長的研究，以確定分析過程中的關鍵參數取值；

（2）通過理論計算確定的臨界時間步長可以作為計算過程穩定性的判據，但是不能作為計算結果精確性的判據，就本文中的FBO 結構特點和物理過程而言，工程可用的時間步長約為理論計算值的1/10左右；

（3）對于某一特定的物理過程和目標時間，參與碰摩、損傷、沖擊等物理過程的構件存在著網格尺度和時間步長的臨界值，進一步提升網格密度和降低時間步長（提升時間分辨率），對于計算精度的提升作用有限，但是卻要耗費大量的計算資源，不參與上述物理過程的構件可以采用較粗的網格以提高計算速度；

（4）碰摩過程對于FBO的顯式動力學仿真中參與碰摩構件的應力和振動載荷結果精度有著重要影響；

（5）由計算方法決定，顯式計算的偏差會隨時間積累逐步放大。

6 展望

本文開展的仿真分析基于某商用顯式動力學軟件平臺，由于不同軟件平臺的計算方法基本一致，但是在軟件設置上和具體參數處理上有一定區別，因此上述結論在其他軟件平臺上的規律、趨勢大體一致，但是量值還需要針對具體的結構特點、物理過程和軟件設置進行專門的分析來確定。

受條件所限，本文對于計算結構偏差僅按照仿真結果的一致性來評判，在真實的工程應用中應該開展若干部件級試驗驗證，通過試驗結果進一步作為仿真過程的評判依據來校準仿真設置，可以提供更加貼近工程實踐的仿真分析結果。