肖春華，林貴平，桂業(yè)偉，肖京平

（1.中國空氣動(dòng)力研究與發(fā)展中心 空氣動(dòng)力學(xué)國家重點(diǎn)實(shí)驗(yàn)室，四川 綿陽 621000；2.北京航空航天大學(xué) 航空科學(xué)與工程學(xué)院，北京 100083）

0 引 言

電熱除冰過程中，飛機(jī)蒙皮表面所結(jié)的冰層一方面受到氣動(dòng)力載荷作用，另一方面又受到熱載荷的作用。在外部氣動(dòng)力約束下，加熱除了融化冰層外，還將在冰層內(nèi)部產(chǎn)生熱應(yīng)力，這一熱力耦合特性是電熱除冰過程中的重要特征。

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

在電加熱條件下，耦合外部氣動(dòng)載荷的作用，采用有限元法研究了電熱除冰的熱力耦合特性及其對(duì)冰層的影響。對(duì)加熱條件下冰層內(nèi)部應(yīng)力的分布特征進(jìn)行了系統(tǒng)研究，比較了加熱和未加熱條件下冰層內(nèi)部應(yīng)力的差別，研究了熱力耦合特性及其對(duì)冰層破壞的影響，這對(duì)于電熱除冰理論和除冰技術(shù)的發(fā)展都具有現(xiàn)實(shí)的意義。

1 數(shù)值方法

1.1 計(jì)算模型

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

結(jié)冰計(jì)算條件：來流速度67m／s、壓強(qiáng)101300Pa、翼型弦長(zhǎng)0.5334m、迎角4°、液態(tài)水含量1.0g／m3、水滴平均直徑20μm、結(jié)冰時(shí)間6min、結(jié)冰環(huán)境溫度－6℃。除冰計(jì)算條件：來流速度170m／s、環(huán)境溫度－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 計(jì)算方法

1.2.1 傳熱問題

（1）控制方程

采用瞬態(tài)熱傳導(dǎo)方程［13－15］

ρ、C、κ分別是冰的密度、比熱和導(dǎo)熱系數(shù)，T是溫度。

（2）定解條件

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

邊界條件：參考國外文獻(xiàn)［1－6］，在電熱除冰過程中，冰層和外流場(chǎng)的交界面一般采用對(duì)流換熱條件，不考慮繼續(xù)結(jié)冰的情況。其中，是表面單位法向矢量，h是對(duì)流換熱系數(shù)，Ta是環(huán)境溫度，下標(biāo)interface代表交界面。

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

1.2.2 應(yīng)力問題

（1）控制方程

采用平面彈性力學(xué)方程［9］，

平衡方程：

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

幾何方程——應(yīng)變－位移關(guān)系：

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

物理方程——應(yīng)力－應(yīng)變關(guān)系：

其中，σ是應(yīng)力向量，ε是應(yīng)變向量，D是彈性矩陣。

平面應(yīng)力：E0＝E，v0＝v

（2）定解條件

幾何邊界條件：

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

力的邊界條件：

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

對(duì)于熱力耦合問題，應(yīng)力計(jì)算所需的溫度載荷由前面的傳熱計(jì)算獲得。

（3）求解方法

為適應(yīng)冰層不規(guī)則的外部形狀，應(yīng)力計(jì)算采用有限元法和非結(jié)構(gòu)網(wǎng)格技術(shù)，計(jì)算單元采用三節(jié)點(diǎn)平面應(yīng)力三角形單元，線性代數(shù)方程組采用超松弛迭代法進(jìn)行求解。冰的物性參數(shù)［14］見表1。

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

2 計(jì)算結(jié)果和分析

2.1 界面應(yīng)力的分布

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

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

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

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

2.2 應(yīng)力隨熱流密度的變化

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

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

圖6是熱流密度20kW／m2時(shí)冰層內(nèi)部第三主應(yīng)力分布云圖，加熱時(shí)間5s。由圖可知，在靠近加熱界面的區(qū)域第三主應(yīng)力的絕對(duì)值最大。隨著加熱時(shí)間的延長(zhǎng)，靠近加熱界面區(qū)域的第三主應(yīng)力的絕對(duì)值越來越大，表現(xiàn)出壓應(yīng)力的不斷增大。圖7是加熱5s時(shí)冰層內(nèi)的溫度分布，由圖可知，此時(shí)冰層內(nèi)部及界面還未開始融化。

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

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

3 熱力耦合對(duì)冰層影響的原理性實(shí)驗(yàn)

前面的計(jì)算表明，在電熱除冰過程中，冰層內(nèi)部將產(chǎn)生很大的應(yīng)力，局部區(qū)域甚至超過了抗壓強(qiáng)度，將導(dǎo)致冰層的破裂，為驗(yàn)證結(jié)果特進(jìn)行此原理性實(shí)驗(yàn)。實(shí)驗(yàn)在氣動(dòng)中心0.3m×0.2m小型結(jié)冰風(fēng)洞中進(jìn)行（圖8），模型采用直徑30mm圓柱形電加熱裝置，共分3層，從外到內(nèi)依次是硬鋁層、中間層和加熱層，模型橫截面示意圖見圖9。

圖8 結(jié)冰風(fēng)洞第二試驗(yàn)段照片F(xiàn)ig.8 Photo of primary test segment of icing wind tunnel

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

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

圖10 裂紋出現(xiàn)的時(shí)間及相應(yīng)的表面溫度Fig.10 Time of crack and temperature of surface

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

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

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

4 結(jié) 論

首先采用有限元法，在電加熱條件下，耦合外部氣動(dòng)載荷的作用，對(duì)冰層的熱力耦合特性進(jìn)行了計(jì)算研究，最后通過原理性實(shí)驗(yàn)對(duì)熱力耦合特性對(duì)冰層的影響進(jìn)行了驗(yàn)證，研究獲得如下結(jié)論：

（1）熱力耦合產(chǎn)生的冰層應(yīng)力非常大，局部區(qū)域明顯超出了抗拉／壓強(qiáng)度。法向應(yīng)力在加熱界面為負(fù)，表現(xiàn)為壓應(yīng)力，而在兩側(cè)界面為正，表現(xiàn)為拉應(yīng)力。隨著熱流密度的增加，界面應(yīng)力變化越來越劇烈，界面應(yīng)力絕對(duì)值明顯增大；

（2）隨著熱流密度的增大，冰層受熱膨脹的程度越來越大，冰層內(nèi)部第三主應(yīng)力的最大絕對(duì)值不斷增大，表現(xiàn)為不斷增大的壓應(yīng)力；

（3）冰層最大壓應(yīng)力明顯超過了抗壓強(qiáng)度，冰層局部區(qū)域?qū)l(fā)生破裂，裂紋的出現(xiàn)發(fā)生在表面冰層融化之前。冰層發(fā)生破裂后，應(yīng)力得到了釋放，隨著加熱的繼續(xù)，冰層有時(shí)會(huì)在主裂紋旁擴(kuò)展出二次裂紋。冰層其它區(qū)域的應(yīng)力未超過抗拉／壓強(qiáng)度，故只產(chǎn)生局部的破裂而非完全的破碎。

［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］ 肖春華，桂業(yè)偉，易賢，等.結(jié)冰翼型表面明冰的壓力分布和應(yīng)力計(jì)算［J］.航空動(dòng)力學(xué)報(bào)，2009.

［11］ 易賢.飛機(jī)積冰的數(shù)值計(jì)算與積冰試驗(yàn)相似準(zhǔn)則研究［D］.四川綿陽：中國空氣動(dòng)力研究與發(fā)展中心，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］ 王勖成，邵敏.有限單元法基本原理和數(shù)值方法［M］.北京：清華大學(xué)出版社，1997.

［15］ 竹內(nèi)洋一郎著.熱應(yīng)力［M］.郭廷瑋，李安定譯.北京：科學(xué)出版社，1977.