旋轉帽罩表面積冰甩脫位置的數值模擬

傅　亮，陳　勇，王安正

（上海交通大學機械與動力工程學院，上海200240）

旋轉帽罩表面積冰甩脫位置的數值模擬

傅亮，陳勇，王安正

（上海交通大學機械與動力工程學院，上海200240）

為了更好地理解旋轉帽罩表面積冰的甩脫機理和過程，采用水杯法在萬能試驗機上分別測得積冰/鋁板界面黏附強度值和積冰內部拉伸強度值。參考試驗測得的這2個數據，結合內聚力模型（CZM）和擴展有限元法（X-FEM），在有限元軟件A N SY S中分別模擬得到積冰與旋轉帽罩之間的界面分離區域和冰形內部裂紋路徑及其擴展過程，對旋轉帽罩表面積冰甩脫的位置和體積作出預測。預測結果與帽罩表面積冰甩脫試驗結果作比較，證明了CZM和X-FEM耦合的數值模擬方法可行。為實際工程應用中的熱氣或者電防冰系統的布置提供參考。

積冰甩脫；黏附強度；拉伸強度；內聚力模型；擴展有限元法；旋轉帽罩；航空發動機

0　引言

當飛機在結冰條件下飛行時，發動機的旋轉帽罩會因為大氣中過冷水滴的撞擊、積聚而結冰。帽罩積冰會影響發動機進氣流場，降低氣動性能，并且甩脫的積冰可能會打壞風扇葉片甚至壓氣機部件，影響飛行安全。因此有必要對帽罩積冰在甩脫的具體位置及體積進行預測，為實際工程應用中的熱氣防冰或電防冰系統的布置提供參考。

積冰/界面黏附強度與很多因素有關，比如基底材料的表面特性、積冰的形成條件、空氣溫度、凍結速度等。國內外學者采用了多種試驗方法測量積冰/金屬界面黏附強度［1］，但不同學者提供的結果相差較大。Loughborough［2］發現冰/銅和冰/鋁界面黏附強度分別為0.85MPa和1.52MPa；Druez［3］發現冰/鋁界面黏附強度在0.075～0.12MPa間隨表面粗糙度的增大而增大；Raraty和Tabor［4］試驗發現，在-6℃和-15℃條件下，冰/拋光不銹鋼界面黏附強度分別為0.03MPa和0.16MPa；丁金波等［5］也通過試驗對不同粗糙度金屬表面的積冰黏附性能做了研究，黏附強度數值在0.02～0.13MPa之間；李清英等［6］利用電脈沖除冰系統來減小積冰的黏附強度；W.Dong等［7］對旋轉帽罩進行熱分析，利用熱氣來減小積冰的黏附強度。

冰內部物理屬性也受很多因素的影響，比如環境溫度、形成方式、雜質比例等，不同試驗測得的冰內拉伸強度差別也很大。Xian和Scavuzzo［8］發現當測試溫度從-23.3℃升高到-3.9℃的過程中，積冰拉伸強度從1.172MPa減小到0.827MPa；Reich［9］試驗測得霜冰和明冰的拉伸強度不一樣，霜冰的拉伸強度是0.12MPa，明冰的拉伸強度為0.62～2.3MPa。

本文利用水杯法分別測量得到積冰/鋁板界面黏附強度值和積冰內部拉伸強度值，然后在有限元軟件ANSYS中，利用內聚力模型定義冰形與帽罩表面接觸單元，模擬得到冰形與旋轉帽罩之間的界面分離區域；利用擴展有限元法模擬得到積冰內部裂紋路徑及其擴展過程，最終預測旋轉帽罩表面積冰甩脫的位置和體積，并與之前學者的試驗結果［10］進行比較。

1　試驗方法

理論上只有當界面分離和冰內部裂紋同時出現時，旋轉帽罩表面的積冰才會發生甩脫，如圖1紫色方框部分所示。積冰與帽罩表面界面分離的難易是由積冰/界面黏附強度決定的，而冰內部裂紋形成的難易是由積冰內部拉伸強度決定的。

圖1　旋轉帽罩表面積冰甩脫

1.1積冰/鋁板界面黏附強度測試

采用水杯法測量積冰/鋁板界面黏附強度如圖2所示。其中1為數據采集計算機，內置測試軟件為Material Test4.1；2為UTM6104萬能拉伸試驗機；3為量程為500 N的拉力傳感器；4為倒置的水杯，杯口直徑為20mm；5為鋁板；6為拉伸機固定底座。測量前，先通過水杯底面的小孔給杯中注入3 mL蒸餾水，約占水杯容積的2/3。然后將水杯和鋁板整體放入冰柜里，冰柜溫度設置為-15℃，凍結24 h后取出，用夾具將鋁板固定在拉伸機底座上，水杯端通過1段柔性連接固定在上方夾具上。設置拉伸機向上位移速度參數為5mm/min，最大位移量為20mm，啟動拉伸機，直到水杯中冰塊與鋁板界面完全分離后停止，輸出拉力與時間的關系圖像，得到過程中最大拉力值Fmax和水杯冰塊本身的重力值Gclamp，則界面黏附強度σn為

式中：S=3.14×10-4m2，為界面接觸面積，即水杯底面積；σn為積冰/鋁板界面黏附強度，為了減小試驗誤差，多次重復試驗取平均值，最終算得0.07MPa。

圖2　積冰/鋁板界面黏附強度測試裝置

圖3　冰拉伸強度試驗裝置及結果

1.2積冰內部拉伸強度測試

同樣采用水杯法測量積冰內部拉伸強度如圖3所示。將2個完全相同的水杯對接好，并注入適量水，放置冰柜中凍結24 h后，2個水杯通過其中的1塊整冰牢固地黏在一起，將2個水杯上下固定在拉伸試驗機上，拉伸機參數設置與第1.1節的相同，啟動拉伸機，直到2個水杯完全分離后停止，輸出拉力與時間的關系（圖3），得到過程中最大拉伸應力F'max和夾具本身的重力G'clamp，則利用式（1）計算得到積冰內部拉伸強度σt=0.38MPa。

2　數值方法

2.1內聚力模型

內聚力模型是由Alfano和Crisfield［11］首先提出的，Chen Y等［12］利用內聚力模型研究發動機風扇葉片表面積冰的甩脫問題。內聚力模型通過ANSYS軟件中的接觸單元來實現，如圖4所示。OB顯示了線彈性加載階段，BD為線性軟化階段。在B點法向/切向接觸應力達到最大值，OB的斜率代表了材料的脆性或黏性。一旦達到材料的開裂強度（B點），其損傷會隨循環載荷的增大而累積，而在起裂后的加載/卸載過程中，其斜率與界面損傷程度di有關，界面應力σt與開裂位移δi的關系式為

式中：i=n，t，m，分別對應純法向型、純剪切型和混合模式裂紋。界面損傷程度di定義為

圖4　內聚力模型

當△i≤1時di=0，而△i≥1時0＜di≤1，其余參數表示為

三角形OBD的面積為界面斷裂過程中釋放的能量即斷裂能Gc，一般通過試驗獲得，其表達式為

積冰/帽罩界面分離既不是單純的法向分離，又不是單純的切向滑移，在這種混合模式下，斷裂能滿足

其中

利用ANSYS中的內聚力模型對積冰/帽罩界面破壞區域進行模擬，模擬過程中最重要的2個參數是界面黏附強度σmax和界面斷裂能Gc，數值模擬中的斷裂能Gc=0.2 N/mm。

2.2擴展有限元法

擴展有限元法由Belytschko和Black［13］首先提出，方修君等［14］利用擴展有限元法對3點彎梁的開裂過程進行了模擬。應用擴展有限元法時，不需要實時地對裂紋尖端區域重新進行網格劃分，因為與傳統有限元法相比，擴展有限元法修正了有限元的位移函數，改進了對不連續邊界的描述方法，使其獨立于網格劃分，擴展有限元法的位移函數為

傳統有限元部分

式中：N（ix）為節點i對應的傳統有限元法形函數；ui為傳統有限元法節點i位移矢量；αi為節點i的擴展自由度變量；F（x）為裂紋尖端彈性漸進函數；j為裂紋尖端節點自由度變量；H（x）為階躍函數

式中：x*為考察點x在裂紋面上的投影點；n為裂紋面的單位外法向。

利用X-FEM模擬冰形內部具體的裂紋路徑，設置正確的裂紋初始單元十分重要，所以要先對冰形內部作應力分析，找到最大應力位置，參考試驗值，設置極限拉伸強度為0.38MPa。

2.3數值模擬流程

分別利用ANSYS中自帶的內聚力模型和擴展有限元方法對積冰/帽罩金屬表面的界面分離及積冰內部的裂紋擴展進行模擬，分析流程如圖5所示。由于缺少冰風洞等試驗條件，帽罩、冰形幾何模型的選擇與文獻［10］的相符。王健等通過試驗手段分析旋轉帽罩表面冰生長/脫落過程，得到如圖6（a）的幾何模型，本文利用GetData提取出帽罩及最厚的冰形曲線，在ANSYS里完成幾何建模，幾何模型中旋轉帽罩的錐底直徑為100mm，錐角為84°。胡婭萍等［15］發現帽罩錐角的不同對積冰影響很大。積冰的物性參數如下，密度為897kg/m3，彈性模量為9000MPa，泊松比為0.325，拉伸強度設為0.38MPa；帽罩部分的物性參數如下，密度為2700kg/m3，彈性模量70 GPa，泊松比為0.3，積冰在帽罩表面的黏附強度設為0.07MPa。

圖5　數值模擬分析流程

有限元模型如圖6所示。其中圖6（b）為帽罩和冰形2維有限元網格，紫色為帽罩部分，有2556個四邊形網格，積冰部分有568個四邊形網格。模擬過程中進行過網格無關性檢驗，發現加密網格對結果影響不大。

圖6　積冰/帽罩幾何和有限元模型

3　數值模擬結果分析

3.1界面分離區域

采用ANSYS內置的內聚力模型模擬界面分離區域，因為幾何模型的對稱性，在冰形和旋轉帽罩表面之間建立接觸對，采用2維3節點接觸單元CONTA172，帽罩表面設為目標面，冰形界面設為接觸面，法向和切向剛度均采用默認值，接觸類型設為初始接觸類型。約束和載荷條件與文獻［10］的一致，在帽罩底面所有節點設置軸向約束，帽罩100%轉速設為1591.74 rad/s，計算求解。積冰與旋轉帽罩接觸應力分布如圖7所示，根據最大應力強度準則，當界面應力等于0.07MPa時，假設該部分界面發生分離，分離界面如圖8紅色區域所示。當帽罩轉速遞增時，損傷系數最大的位置始終不變，距積冰底面17.5mm，如圖8所示。

圖7　積冰/帽罩界面接觸應力分布

圖8　不同轉速下界面損傷系數di分布

3.2冰形內部應力

在模擬積冰內裂紋路徑之前，要先對積冰作整體應力分析，找到冰形中拉伸應力最大的位置，以便準確設置裂紋初始單元。計算過程中的接觸、約束條件、邊界條件設置參照第3.1節，積冰在特定轉速（70%、80%、90%、100%）下內部拉伸應力分布如圖9所示。從圖中可見，隨著轉速的增加，最大應力位置保持不變（恰好距離積冰最底部也是17.5mm，與圖8界面損傷系數最大的位置相一致），該位置是冰內最容易發生斷裂的位置，應該在該區域設置裂紋初始單元。

圖9　不同轉速冰形內部拉伸應力分布

3.3冰形內部裂紋擴展

擴展有限元法模擬冰形內部具體的裂紋位置，以最大拉應力理論作為積冰內部斷裂準則，極限拉伸強度為0.38MPa。由于冰形周向對稱，所以可以將整體冰形簡化為2維進行處理，總時間設為1，最小時間步長設為0.0001。裂紋擴展路徑如圖10所示。從圖中可見，裂紋近似平行于軸線方向，裂紋尖端的拉伸應力為0.38MPa，與極限拉伸強度相一致。裂紋擴展過程4個時刻裂紋附近的位移如圖11所示，裂紋擴展順序由a至d。

3.4模擬結果驗證

王健等［10］利用高速攝影機拍攝的帽罩旋轉過程中積冰甩脫瞬間如圖12（a）所示。從圖中可見，試驗中距離積冰底部5～20mm的區域的積冰發生甩脫。本文積冰甩脫位置的模擬結果如圖12（b）所示，其中紫色線框部分是積冰預測甩脫區域，與積冰底部的距離在3～17.5mm之間，積冰甩脫的位置都是在靠近積冰底部的后半段，說明本文采用的內聚力模型和擴展有限元法耦合的數值模擬方法可行。

圖10　裂紋完全擴展拉伸應力

圖11　裂紋擴展過程

圖12　試驗結果與模擬結果的對比

應力/MPa .189515 42232.3 84464.4 126696 168929 211161 253393 295625 337857 380089

4　結論

本文利用水杯法分別測量積冰/鋁板界面黏附強度值σn和積冰內部拉伸強度值σt；并采用數值模擬方法，利用內聚力模型對界面接觸單元進行定義，對界面分離區域進行預測，再利用擴展有限元法對積冰內部應力位置進行裂紋初始單元設置，進而有效模擬出積冰內部裂紋擴展路徑，最終確定積冰甩脫位置，主要結論如下：

（1）內聚力模型與擴展有限元法的耦合能有效地模擬預測得到帽罩表面積冰甩脫的位置，也為其他部件表面積冰甩脫的分析提供新的思路。

（2）當試驗環境溫度在-15℃時，由水杯法測得的積冰/鋁板界面黏附強度為0.07MPa，積冰內部拉伸強度為0.38MPa。

［1］Javan-Mashmool M，Volat C，Farzaneh M.A new method for measuring ice adhesion strength at anice-substrate interface［J］.Hydrological Processes，2006，20（4）：645-655.

［2］Loughborough D L.Reduction of the adhesion of ice to de-icer surfaces［J］.Journal of the Aeronautical Sciences，2012，13（3）：126-134.

［3］Druez J，Phan C L，Laforte J L.Adhesion of glaze and rime on aluminium electrical conductors［J］.Transactions-Canadian Society for Mechanical Engineering，1978，5（4）：215-220.

［4］Raraty L E，Tabor D.The adhesion and strength properties of ice［J］. Proceedings of the Royal Society a Mathematical Physical&Engineering Sciences，1958（6）：184-201.

［5］丁金波，董威.表面粗糙度對冰凍黏強度影響試驗研究［J］.航空發動機，2012，38（4）：42-46. DING Jingbo，DONG Wei.Experimental study of influence of surface roughness on ice adhesion［J］.Aeroengine，2012，38（4）：42-46.（in Chinese）

［6］李清英，朱春玲，白天.電脈沖除冰系統的除冰試驗與數值模擬［J］.航空動力學報，2012，27（2）：350-356. LI Qingying，ZHU Chunlin，BAI Tian.De-icing experiment and numerical simulation of the electro-impulse de-icing system［J］.Journal of Aerospace Power，2012，27（2）：350-356.（in Chinese）

［7］Dong W，Zhu J，Zheng M.Thermal analysis and testing of a cone with leading edge hot air anti-icing system［R］.AIAA-2014-0739.

［8］Xian X，Chu M L，Scavuzzo R J，et al.An experimental evaluation of the tensile strength of impact ice［J］.Journal of Materials Science Letters，1978，8（10）：119-130.

［9］Reich A D.Ice property/structure variations across the glaze/rime transition［R］.AIAA-1992-0296.

［10］王健，胡婭萍，吉洪湖，等.旋轉整流罩積冰生長與脫落過程的試驗［J］.航空動力學報，2014，29（6）：1352-1357. WANG Jian，HU Yaping，JI Honghu，et al.Experiment of ice accrection and shedding on rotating spinner［J］.Journal of Aerospace Power，2014，29（6）：1352-1357.（in Chinese）

［11］Alfano G，Crisfield M A.Finite element interface models for the delamination analysis of laminated composites：mechanical and computational issues［J］.International Journal for Numerical Methods in Engineering，2001，50：1701-1736.

［12］Chen Y，Dong W，Wang Z.Numerical simulation of ice shedding from a fan blade［R］.ASME 2015-GT-42265.

［13］Belytschko T，Black T.Elastic crack growth in finite elements with minimal remeshing［J］.International Journal for Numerical Methods in Engineering，1999，45（5）：601-620.

［14］方修君，金峰，王進廷.基于擴展有限元法的粘聚裂紋模型［J］.清華大學學報，2007，47（3）：344-347. FANG Xiujun，JIN Feng，WANG Jinting.Cohesive crack model based on extended finite element method［J］.Journal of Tsinghua University，2007，47（3）：344-347.（in Chinese）

［15］胡婭萍，吉洪湖，王健，等.錐角對旋轉整流罩積冰影響的模擬試驗［J］.航空動力學報，2014，29（3）：495-503. HU Yaping，JI Honghu，WANG Jian，et al.Experiment on effect of cone angle on ice accretion of rotating spinner［J］.Journal of Aerospace Power，2014，29（3）：495-503.（in Chinese）

（編輯：栗樞）

Numerical Simulation of Ice Shedding from Spinning Cone

FU Liang，CHEN Yong，WANG An-zheng
（School of Mechanical Engineering，Shanghai Jiaotong University，Shanghai 200240，China）

In order to understand the mechanism and process of ice shedding for spinning cone，the cup method was introduced to measure the ice/Aluminum interfacial adhesive strength and the ice internal tensile strength.The Cohesive Zone Model（CZM）and the eXtended Finite Element Method（X-FEM）were introduced to simulate the ice shedding process from a spinning cone.The CZM was applied to simulate the initiation and propagation of ice/spinning cone interface crack.The X-FEM was introduced to model crack growth inside the ice.The research method provides a reference for ice shedding analysis of other parts of aeroengine.The comparation between simulation and experiment results show that the method is feasible and could also be used to further study optimization of the electrothermal and hot-air deicing system.

ice shedding；interfacial strength；tension strength；Cohesive Zone Model；eXtended Finite Element Method；spinning cone；aeroengine

V 228.8

Adoi：10.13477/j.cnki.aeroengine.2016.05.012

2016-04-30基金項目：國家自然科學基金（51376122）資助

傅亮（1990），男，在讀碩士研究生，研究方向為航空發動機結冰力學分析；E-mail：984292889@qq.com。

引用格式：傅亮，陳勇，王安正.旋轉帽罩表面積冰甩脫位置的數值模擬［J］.航空發動機，2016，42（5）：70-75FULiang，CHENYong，WANGAnzheng. Numericalsimulationoficesheddingfromspinningcone［J］.Aeroengine，2016，42（5）：70-75.