電熱除冰的熱力耦合特性及其對冰層的影響研究

肖春華，林貴平，桂業偉，肖京平

（1.中國空氣動力研究與發展中心 空氣動力學國家重點實驗室，四川 綿陽 621000；2.北京航空航天大學 航空科學與工程學院，北京 100083）

0 引 言

電熱除冰過程中，飛機蒙皮表面所結的冰層一方面受到氣動力載荷作用，另一方面又受到熱載荷的作用。在外部氣動力約束下，加熱除了融化冰層外，還將在冰層內部產生熱應力，這一熱力耦合特性是電熱除冰過程中的重要特征。

以往的電熱除冰計算研究主要圍繞除冰結構的傳熱、冰層融化和冰層的力學特性問題，目前針對電熱除冰過程中熱力耦合特性的研究還比較少。研究人員［1－6］陸續對電熱除冰問題進行了建模和計算，研究了多層電熱除冰結構的傳熱特性和冰層融化特性，但沒有考慮加熱對冰層力學特性的影響。另外，研究人員［7－10］對氣動載荷作用下的冰層內部應力分布進行了計算研究，研究發現低速條件下氣動載荷對冰層脫離的貢獻不大［9－10］，但研究未考慮加熱對冰層力學特性的影響，更沒有研究熱力耦合特性及其對冰層的破壞影響。因此，很有必要對電熱除冰過程中的熱力耦合特性及其對冰層的影響進行研究，以預測和掌握冰層的破壞規律以及對冰脫離的貢獻。

在電加熱條件下，耦合外部氣動載荷的作用，采用有限元法研究了電熱除冰的熱力耦合特性及其對冰層的影響。對加熱條件下冰層內部應力的分布特征進行了系統研究，比較了加熱和未加熱條件下冰層內部應力的差別，研究了熱力耦合特性及其對冰層破壞的影響，這對于電熱除冰理論和除冰技術的發展都具有現實的意義。

1 數值方法

1.1 計算模型

帶冰翼型為NACA0012翼型，表面所結冰型是明冰，外形由文獻［11］獲得。翼型的前緣局部放大后，得到明冰外形的清晰型線，見圖1。冰型外表面受到外部氣動載荷作用，內表面受到蒙皮的附著約束作用，當氣動載荷導致的冰層應力超過冰層和蒙皮間的抗剪強度時，冰層才可能脫離蒙皮表面。冰層的剝離主要考慮界面應力對冰層和蒙皮間的剪切破壞作用［12］，而冰層的破裂主要依據內部應力對冰層的壓縮／拉伸破壞作用，文中的抗壓強度取1.0MPa。界面應力分布的橫坐標s是沿冰層界面從上表面根部到下表面根部的曲線段長度，加熱界面是從弦長s1（0.018m）到s2（0.034m）處，參見圖2。

結冰計算條件：來流速度67m／s、壓強101300Pa、翼型弦長0.5334m、迎角4°、液態水含量1.0g／m3、水滴平均直徑20μm、結冰時間6min、結冰環境溫度－6℃。除冰計算條件：來流速度170m／s、環境溫度－6℃、加熱功率104～4×104W／m2。

圖1 翼型前緣明冰局部放大圖Fig.1 Local blowup of glaze ice on airfoil front

圖2 加熱界面熱載荷分布圖Fig.2 Distribution of heat flux at interface

1.2 計算方法

1.2.1 傳熱問題

（1）控制方程

采用瞬態熱傳導方程［13－15］

ρ、C、κ分別是冰的密度、比熱和導熱系數，T是溫度。

（2）定解條件

初始條件：取冰的初始溫度為T0。

邊界條件：參考國外文獻［1－6］，在電熱除冰過程中，冰層和外流場的交界面一般采用對流換熱條件，不考慮繼續結冰的情況。其中，是表面單位法向矢量，h是對流換熱系數，Ta是環境溫度，下標interface代表交界面。

冰層和蒙皮的界面采用恒熱流條件，加熱界面從弦長s1到s2處的熱流密度為Qconst，其它不加熱界面采用絕熱條件。

1.2.2 應力問題

（1）控制方程

采用平面彈性力學方程［9］，

平衡方程：

其中，fx、fy分別是x、y方向的單位體積力分量，σx、σy、τxy分別是x、y方向的正應力以及剪切應力。

幾何方程——應變－位移關系：

其中εx、εy、γxy分別是x、y方向的正應變及剪應變，u、v是x、y方向位移。

物理方程——應力－應變關系：

其中，σ是應力向量，ε是應變向量，D是彈性矩陣。

平面應力：E0＝E，v0＝v

（2）定解條件

幾何邊界條件：

其中，u、v是邊界上的位移，、是已知的位移。假定機翼是不發生變形的剛性體，冰層外表面是自由位移條件，冰層和蒙皮間是緊密粘附的無相對位移約束條件，當加熱界面的冰層融化變成水，由于水被包圍在冰層內，也采用無相對位移約束條件。

力的邊界條件：

其中，Fx和Fy是邊界上單位面積的力和是已知的面積力，這里指流場壓力，通過求解低速粘流的時均N－S方程組獲得［12］。直接將壓力載荷加到冰層外表面每個單元的界面中心，方向沿表面單元的法向。加熱界面在融化后只受壓不受剪切作用。

對于熱力耦合問題，應力計算所需的溫度載荷由前面的傳熱計算獲得。

（3）求解方法

為適應冰層不規則的外部形狀，應力計算采用有限元法和非結構網格技術，計算單元采用三節點平面應力三角形單元，線性代數方程組采用超松弛迭代法進行求解。冰的物性參數［14］見表1。

表1 冰型的物性參數表Table1 Material properties of the glaze ice

2 計算結果和分析

2.1 界面應力的分布

圖3～4是不同熱流密度下界面法向、切向應力的分布，加熱時間10s。由圖可知，隨著熱流密度的增大，應力變化越來越劇烈。在加熱界面附近，熱載荷產生的冰層應力非常大，甚至超過了抗拉／壓強度，相比之下，氣動載荷導致的冰層應力非常?。▓D3、4中的水平實線），幾乎可忽略。加熱界面的冰層受熱最嚴重，冰層受熱產生膨脹，膨脹受到阻礙產生較大的壓應力，所以界面法向應力在加熱界面為負，表現為壓應力。加熱界面兩側受熱較少，法向應力為正，表現為拉應力。加熱10s時，加熱界面的冰層還未融化，此時切向應力在加熱界面為正，表現為往上的剪切作用，在加熱界面兩側為負，表現為往下的剪切作用；而當加熱界面的冰層融化后，切向應力等于零。單位時間冰層受熱越厲害，造成加熱面及附近區域的溫差越大，冰塊受熱膨脹程度越大，導致該位置的壓應力越來越大，應力絕對值隨熱流密度的增大明顯增大，法向應力在加熱面附近區域均出現最大應力峰值。

圖3 不同熱流密度下界面法向應力分布Fig.3 Interface normal stress distribution under different heat fluxes

圖4 不同熱流密度下界面切向應力分布Fig.4 Interface shear stress distribution under different heat fluxes

圖中實曲線是熱流密度為零也即不加熱條件下的界面法向應力和切向應力的分布。可以看出，不加熱條件下的界面應力非常小，而加熱條件下的界面應力非常大，甚至超過了抗壓強度，造成冰層的破壞。所以，在表面冰層融化前，電加熱條件下，耦合氣動載荷，熱力耦合很容易造成冰層的破裂。加熱區及附近的冰與翼型界面應力大大超過氣動載荷的約束，但由于其它部位的界面應力還沒有超過抗壓強度，冰與物面間的粘附力仍然存在，所以在一定的加熱條件下冰沒有發生脫落，而只是發生了局部破裂。

2.2 應力隨熱流密度的變化

圖5是冰層第三主應力最大值隨熱流密度的變化，加熱時間分別為5、10、15和20s。由圖可知，相同時間下，冰層第三主應力最大絕對值隨熱流密度的增加逐漸增大，表現為壓應力不斷增大。因為熱流密度越大，相同時間內冰層的熱膨脹就越嚴重，由于外部氣動載荷和界面的約束，膨脹越顯著則產生的壓應力就越大，導致第三主應力最大值隨熱流密度的增加而增大。相同熱流密度下，第三主應力的絕對值隨加熱時間不斷增大。因為相同熱流密度下，加熱時間越長，單位時間內傳遞到冰層的熱量就越多，冰層受熱就越嚴重，膨脹也就越顯著，所以壓應力就越大。圖5中多數冰層的最大壓應力超過了抗壓強度，有的甚至達到4.0MPa，根據破壞條件，預測冰層內部將發生破裂現象。如果冰層具有均勻的物性參數，那么產生最大壓應力的部位將率先破裂。當冰層內部發生破裂后，應力得到釋放，如果要繼續產生破裂，需要進一步加熱造成更大的溫差，從而產生更大的應力。冰層內部其它部位的應力還未超過抗拉／壓強度，所以冰層產生的是局部的破裂而非完全的破碎。

圖5 第三主應力最大值隨熱流變化Fig.5 The third principal stress maximum vs.heat flux

圖6是熱流密度20kW／m2時冰層內部第三主應力分布云圖，加熱時間5s。由圖可知，在靠近加熱界面的區域第三主應力的絕對值最大。隨著加熱時間的延長，靠近加熱界面區域的第三主應力的絕對值越來越大，表現出壓應力的不斷增大。圖7是加熱5s時冰層內的溫度分布，由圖可知，此時冰層內部及界面還未開始融化。

圖6 第三主應力分布云圖（5s）Fig.6 Contour of the third principal stress（5s）

圖7 冰層內部溫度分布云圖（5s）Fig.7 Contour of temperature in ice（5s）

3 熱力耦合對冰層影響的原理性實驗

前面的計算表明，在電熱除冰過程中，冰層內部將產生很大的應力，局部區域甚至超過了抗壓強度，將導致冰層的破裂，為驗證結果特進行此原理性實驗。實驗在氣動中心0.3m×0.2m小型結冰風洞中進行（圖8），模型采用直徑30mm圓柱形電加熱裝置，共分3層，從外到內依次是硬鋁層、中間層和加熱層，模型橫截面示意圖見圖9。

圖8 結冰風洞第二試驗段照片Fig.8 Photo of primary test segment of icing wind tunnel

圖9 模型橫截面示意圖Fig.9 Sketch of model cross section

圖10是裂紋出現的時間及相應的表面溫度。結冰條件：風速30m／s、結冰控制溫度－8℃、結冰時間5min。除冰條件：采用連續性加熱模式、電壓100V（圖10（a））和200V（圖10（b））。由圖可知，熱力耦合將造成冰層的破裂，冰層出現裂紋時表面溫度均低于融點，表明裂紋的出現發生在表面冰層融化之前。裂紋出現的時間越早，相應的表面溫度就越低。

圖10 裂紋出現的時間及相應的表面溫度Fig.10 Time of crack and temperature of surface

圖11是加熱前后的冰層比較，包括加熱前、加熱30s后和加熱120s后的冰層。結冰條件：來流風速30m／s、結冰控制溫度為－8℃。除冰條件：采用連續性加熱模式、電壓200V。觀察裂紋出現前后的冰層可以發現，表面冰層還沒有融化時，裂紋就已經出現了。開啟電熱除冰裝置加熱30s后，冰層中下部出現了一條大約3～5cm長的裂紋，裂紋表面的方向和圓柱軸向呈很小的角度，此時的表面溫度顯示表面冰層還未融化。裂紋出現在表面冰層融化之前，先出現裂紋后發生融化，裂紋是由電熱除冰過程中熱力耦合產生的冰層內部應力造成的。加熱120s后，裂紋有輕微的擴大。隨著加熱的繼續，冰層有時會在主裂紋旁擴展出二次裂紋（圖12）。

圖11 加熱前后的冰層比較Fig.11 Comparison of ice before and after heating

圖12 出現主裂紋和二次裂紋的冰層Fig.12 Ice of main crack and two cracks

4 結 論

首先采用有限元法，在電加熱條件下，耦合外部氣動載荷的作用，對冰層的熱力耦合特性進行了計算研究，最后通過原理性實驗對熱力耦合特性對冰層的影響進行了驗證，研究獲得如下結論：

（1）熱力耦合產生的冰層應力非常大，局部區域明顯超出了抗拉／壓強度。法向應力在加熱界面為負，表現為壓應力，而在兩側界面為正，表現為拉應力。隨著熱流密度的增加，界面應力變化越來越劇烈，界面應力絕對值明顯增大；

（2）隨著熱流密度的增大，冰層受熱膨脹的程度越來越大，冰層內部第三主應力的最大絕對值不斷增大，表現為不斷增大的壓應力；

（3）冰層最大壓應力明顯超過了抗壓強度，冰層局部區域將發生破裂，裂紋的出現發生在表面冰層融化之前。冰層發生破裂后，應力得到了釋放，隨著加熱的繼續，冰層有時會在主裂紋旁擴展出二次裂紋。冰層其它區域的應力未超過抗拉／壓強度，故只產生局部的破裂而非完全的破碎。

［1］ BALIGA G.Numerical simulation of one－dimensional heat transfer in composite bodies with phase change［D］.To－ledo Ohio：University of Toledo，1980.

［2］ CHAO D F.Numerical simulation of two－dimensional heat transfer in composite bodies with application to deicing of aircraft components［D］.Toledo Ohio：University of Toledo，1983.

［3］ LEFFEL K L.A numerical and experimental investigation of electrothermal aircraft deicing［D］.Toledo Ohio：University of Toledo，1986.

［4］ MASIULANIEC K C.A numerical simulation of the full two－dimensional electrothermal de－icer pad［D］.Toledo Ohio：University of Toledo，1987.

［5］ WRIGHT W B，KEITH T G，DeWITT K J.Transient two－dimensional heat transfer through a composite body with application to deicing of aircraft components［R］.AIAA－88－0358，1988.

［6］ YASLIK A D，DEWITT K J，KEITH T G.Further developments in three－dimensional numerical simulation of electrothermal deicing systems［R］.AIAA－92－0528，1992.

［7］ SCAVUZZO R J，CHU M L，KELLACKEY C J.Structural analysis and properties of impact ices accreted on aircraft structures［R］.NASA－CR－198473，1996.

［8］ SCAVUZZO R J，CHU M L，OLSEN W A.Structural properties of impact ices accreted on aircraft structures［R］.NASA－CR－179580，1987.

［9］ SCAVUZZO R J，CHU M L，ANANTHASWAMY V.Influence of aerodynamic forces in ice shedding［J］.Journal of Aircraft，1994，31（3）：526－530.

［10］ 肖春華，桂業偉，易賢，等.結冰翼型表面明冰的壓力分布和應力計算［J］.航空動力學報，2009.

［11］ 易賢.飛機積冰的數值計算與積冰試驗相似準則研究［D］.四川綿陽：中國空氣動力研究與發展中心，2007.

［12］ LYNCH D A III，LUDWICZAK D R.Shear strength analysis of the aluminum／ice adhesive bond［R］.NASACR－107162，1996.

［13］ HUANG J R.Numerical simulation of an electrothermally de－iced aircraft surface using the finite element method［D］.Toledo Ohio：University of Toledo，1993.

［14］ 王勖成，邵敏.有限單元法基本原理和數值方法［M］.北京：清華大學出版社，1997.

［15］ 竹內洋一郎著.熱應力［M］.郭廷瑋，李安定譯.北京：科學出版社，1977.