郁嘉,趙柏陽,卜雪琴,林貴平,李志茂
(1.北京航空航天大學(xué)航空科學(xué)與工程學(xué)院,北京100191;2.航空科學(xué)技術(shù)國家實驗室,北京100191; 3.中航工業(yè)第一飛機設(shè)計研究院,陜西西安710089; 4.上海飛機設(shè)計研究院環(huán)控氧氣部,上海201203)
某型飛機發(fā)動機短艙熱氣防冰系統(tǒng)性能數(shù)值模擬
郁嘉1,2,趙柏陽3,卜雪琴1,*,林貴平1,李志茂4
(1.北京航空航天大學(xué)航空科學(xué)與工程學(xué)院,北京100191;2.航空科學(xué)技術(shù)國家實驗室,北京100191; 3.中航工業(yè)第一飛機設(shè)計研究院,陜西西安710089; 4.上海飛機設(shè)計研究院環(huán)控氧氣部,上海201203)
使用三維內(nèi)外強固傳熱耦合方法計算校核發(fā)動機短艙熱氣防冰系統(tǒng)的性能,并分析發(fā)動機進氣流量對蒙皮表面溫度的影響。內(nèi)、外部表面?zhèn)鳠嵯禂?shù)計算均采用純?nèi)S的CFD方法,在內(nèi)、外部網(wǎng)格數(shù)據(jù)交互時使用了距離加權(quán)反比插值法。通過計算獲得發(fā)動機短艙的局部水收集系數(shù)、蒙皮表面溫度的分布情況、各處溢流水量,并由此判定此防冰系統(tǒng)性能是否達到要求。分析表明此發(fā)動機短艙熱氣防冰系統(tǒng)符合防冰性能要求;當(dāng)發(fā)動機進氣流量增大時,蒙皮表面溫度下降,且溢流水量增加。
發(fā)動機;熱氣防冰系統(tǒng);表面溫度;發(fā)動機進氣流量
飛機結(jié)冰是危害飛行安全、導(dǎo)致飛行事故的重要因素,機翼和發(fā)動機短艙結(jié)冰對飛行安全的危害尤其嚴重。飛機結(jié)冰不僅會直接破壞飛機的氣動外形,導(dǎo)致升力下降。對于發(fā)動機而言,結(jié)冰還可引起發(fā)動機喘振,甚至冰塊脫落,損壞發(fā)動機葉片,嚴重時甚至?xí)?dǎo)致發(fā)動機熄火??梢妼Πl(fā)動機短艙進行防除冰設(shè)計的意義非常重大。目前飛機機翼與發(fā)動機多采用熱氣防冰系統(tǒng),需要從發(fā)動機引氣,但從發(fā)動機引氣的同時也會給發(fā)動機帶來一定的代償損失。因此必須在有限的發(fā)動機熱氣供給情況下保證防冰性能。
對飛機結(jié)冰的研究通常主要有三種:結(jié)冰飛行試驗、結(jié)冰風(fēng)洞試驗、數(shù)值模擬。其中前兩種試驗方式需耗費較大的人力物力,且由于條件限制,無法覆蓋所有的結(jié)冰條件。近年來,隨著計算機硬件成本的不斷降低,數(shù)值模擬方式在經(jīng)濟成本上較其他兩者有了巨大的優(yōu)勢,所以數(shù)值模擬方法逐漸成為了研究飛機結(jié)冰/防冰問題的重要手段[1-6]。目前國內(nèi)展開的防冰性能研究多是針對飛機機翼[7-10],針對發(fā)動機短艙的較少[11-13]。本文針對某飛機發(fā)動機短艙在給定條件下的防冰性能進行數(shù)值分析計算,采用了歐拉法計算水滴撞擊特性,利用距離反比插值法進行內(nèi)外流場數(shù)據(jù)交互,得到了短艙表面溫度和溢流水結(jié)果,并分析了發(fā)動機進氣流量對表面溫度的影響。
1.1 三維防冰表面溫度計算方法
影響熱氣防冰表面溫度的因素有外部空氣流動與換熱、腔內(nèi)部熱氣流動與換熱、蒙皮導(dǎo)熱,因此表面溫度是多場耦合的結(jié)果。計算防冰表面溫度時需要考慮多個計算區(qū)域:外部空氣-水滴兩相流場區(qū)域、發(fā)動機短艙蒙皮外表面上的水膜、蒙皮、防冰腔內(nèi)部,計算時需要對這4個計算區(qū)域進行耦合,最終得到收斂后的發(fā)動機短艙表面溫度。具體流程包括:計算外部整機空氣流場、在空氣流場的基礎(chǔ)上計算水滴場、計算防冰腔內(nèi)部熱氣流場、對蒙皮內(nèi)外熱流進行耦合計算[8]。
1.2 外部流場及水滴場計算
選取半整機模型劃分網(wǎng)格,采用了混合型網(wǎng)格。大部分區(qū)域采用非結(jié)構(gòu)化六面體網(wǎng)格繪制,僅在發(fā)動機附近流場區(qū)域使用結(jié)構(gòu)化網(wǎng)格繪制并對唇口附近進行加密,以保證短艙附近計算結(jié)果的準(zhǔn)確性與精確度。然后將劃分好的網(wǎng)格導(dǎo)入Fluent軟件中進行外部流場的計算。
空氣流場計算時,由于飛行速度較大且發(fā)動機周圍流場復(fù)雜,計算時采用Fluent內(nèi)置的三維粘性可壓流求解器。計算時使用如下的邊界條件:外部流場邊界使用遠場邊界條件,并將發(fā)動機入口段某截面處設(shè)為壓力出口,用匹配流量的方式模擬發(fā)動機的抽吸效應(yīng)對流場的影響。發(fā)動機外壁面采用恒壁溫邊界條件。計算收斂后得到空氣速度和對壁面的對流換熱熱流量。
圖1 計算網(wǎng)格Fig.1 Computational grid
水滴場使用Fluent中的用戶自定義標(biāo)量UDS (User Defined Scalars)進行計算。分別設(shè)置4個UDS:水滴x、y、z方向的速度、水滴相對容積分數(shù),并建立水滴場的連續(xù)性和動量方程,使用Fluent自帶求解器進行求解。計算時使用如下的邊界條件:流場入口處認為水滴速度與空氣速度一致并且水滴相對容積分數(shù)為1,認為水滴可以自由通過內(nèi)、外流場出口;當(dāng)水滴速度指向壁面時,認為水滴可以自由通過,反之指向流場時,說明此處為水滴遮蔽區(qū),局部水收集系數(shù)β=0。
使用上述方法計算得到水滴場相應(yīng)參數(shù)后,由定義可知歐拉法下β的計算公式為:
式中,u為水滴速度方向;n為發(fā)動機短艙壁面的單位法向矢量;V∞為自由流水滴速度;αn為水滴的相對容積分數(shù),如下式:
式中,α∞是自由流中水滴的容積分數(shù)。
1.3 防冰腔內(nèi)部流場計算
選取計算的是發(fā)動機短艙笛形管射流防冰系統(tǒng)。為計算防冰表面溫度分布,需要先對防冰腔內(nèi)部的熱氣流動情況進行計算。使用的計算模型為整個發(fā)動機短艙環(huán)狀防冰腔。射流孔均勻分布于整個笛型管上,熱氣從射流孔流出加熱蒙皮,然后在腔內(nèi)流動,并最終從排氣孔流出至外部環(huán)境。網(wǎng)格分為蒙皮導(dǎo)熱區(qū)域及內(nèi)部流動區(qū)域,內(nèi)部流場網(wǎng)格總數(shù)965萬,其中固體區(qū)域77萬,流體區(qū)域888萬,全部使用結(jié)構(gòu)化網(wǎng)格繪制。計算時使用如下的邊界條件:笛形管射流孔出口為計算區(qū)域的入口,設(shè)為壓力入口條件;排氣孔為計算區(qū)域的出口,設(shè)為壓力出口條件;防冰腔外壁面假設(shè)為均勻溫度293.15K,作為表面溫度迭代計算時的初始溫度。
圖2 發(fā)動機短艙外部及內(nèi)部構(gòu)造Fig.2 External and internal structure of engine nacelle
1.4 內(nèi)外流場數(shù)據(jù)交互
由于是通過外部網(wǎng)格計算來獲取短艙表面的局部水收集系數(shù)、熱流、剪切力等數(shù)據(jù),而計算防冰腔表面溫度時則是使用內(nèi)部網(wǎng)格計算,所以需要將外部網(wǎng)格計算時獲得的結(jié)果數(shù)據(jù)導(dǎo)入至內(nèi)部網(wǎng)格相應(yīng)的網(wǎng)格點上。使用距離加權(quán)反比插值法[14]來進行插值計算,計算公式如下:
1.5 蒙皮傳熱耦合計算
耦合計算蒙皮傳熱過程時需同時考慮蒙皮導(dǎo)熱及內(nèi)部流動換熱兩個方面,本文使用強固耦合方法[8]對其進行迭代求解,具體求解過程如下:
1)使用CFD方法計算得到蒙皮外表面換熱系數(shù);
2)對蒙皮外表面加載溫度邊界條件,使用UDF計算蒙皮表面水蒸發(fā)量、溢流水量及熱載荷;
3)引入松弛系數(shù)ω1,將松弛處理后的熱載荷q·″k作為邊界條件加載至蒙皮外表面,使用Fluent計算蒙皮導(dǎo)熱及內(nèi)部流動換熱,并得到蒙皮外表面溫度;
4)引入松弛系數(shù)ω2,對此溫度進行松弛處理,重新作為邊界條件加載至蒙皮外表面;
5)重復(fù)步驟1)~步驟4),直至熱載荷以及表面溫度計算結(jié)果收斂。
本文計算發(fā)動機短艙防冰性能情況,不考慮結(jié)冰過程,此時蒙皮熱載荷主要由如下4項構(gòu)成:
1)蒙皮與外部氣流之間的對流換熱;
2)蒙皮表面水蒸發(fā)所帶走的熱量;
3)外部流場中的水撞擊至蒙皮表面所帶來的能量;
4)水流入/流出計算區(qū)域所帶來/帶走的能量;
蒙皮總熱載荷計算公式如下:
式中:mimp、mevap、min、mout分別為撞擊水質(zhì)量流量、蒸發(fā)水質(zhì)量流量、流進控制體的溢流水質(zhì)量流量、流出控制體的溢流水質(zhì)量流量,cp,w為水的比熱容,ilv為汽化潛熱,Tref為參考溫度273.15 K,Trec為恢復(fù)溫度,定義如下:
式中:rc為恢復(fù)系數(shù),Te為附面層外邊界處溫度。
在計算表面換熱系數(shù)時采用CFD方法,并且在內(nèi)、外網(wǎng)格數(shù)據(jù)交互時,使用距離加權(quán)反比插值算法來保證數(shù)據(jù)準(zhǔn)確性。由于此方法在計算時是一個純?nèi)S過程,相比以前采用附面層積分方法計算n個二維截面(外部流場的)的對流換熱系數(shù),然后插值到三維內(nèi)部流場計算網(wǎng)格的外表面上,本文的方法更適合于三維問題。
2.1 計算狀態(tài)說明
使用上述發(fā)動機短艙熱氣防冰系統(tǒng)作為計算對象,進行三維內(nèi)外耦合計算,獲得了局部水收集系數(shù)、蒙皮外表面溫度、表面溢流水等結(jié)果。
算例已知飛行氣象條件:飛行馬赫數(shù)0.36,液態(tài)水含量2.3 g/m3,攻角7°,環(huán)境溫度266 K,壓力75 265Pa,水滴直徑20 μm。防冰系統(tǒng)引氣狀態(tài):熱氣噴口總壓211 313Pa,總溫493K。
2.2 計算結(jié)果
以下各曲線圖均以蒙皮外表面處的弧長作為x軸,相應(yīng)的計算結(jié)果值作為y軸。s=0為幾何駐點處,s為負數(shù)代表發(fā)動機短艙內(nèi)側(cè)表面,反之s為正數(shù)代表發(fā)動機短艙外側(cè)表面。為更好地對計算結(jié)果進行展示,截取了5個切面分別進行對比分析,取法如圖3所示。
插值計算結(jié)果中的局部水收集系數(shù)β值云圖如圖3所示。從圖3中可知,撞擊水量最大的區(qū)域仍位于幾何駐點附近。
圖3 發(fā)動機短艙表面局部水收集系數(shù)云圖Fig.3 Contour of the local water collection coefficient on engine nacelle
圖4、圖5為內(nèi)外耦合計算收斂后的蒙皮外表面溫度結(jié)果。從圖中可以看出:
1)由于笛形管射流孔均朝向短艙內(nèi)側(cè)表面,從圖中可以明顯看出,發(fā)動機短艙內(nèi)側(cè)表面溫度較外側(cè)高。
2)笛形管射流孔射流正對位置處的蒙皮表面明顯偏高,但遠離此區(qū)域后快速衰減,并達到周圍平均溫度。
圖4 發(fā)動機短艙表面溫度云圖Fig.4 Temperature contour of engine nacelle surface
圖5 發(fā)動機短艙表面溫度曲線圖Fig.5 Temperature curves of engine nacelle surface
3)0°、90°、225°、270°截面位于一列單孔和一列雙孔之間,90°和270°截面更靠近單孔,0°和225°截面更靠近雙孔。所以0°、90°、225°、270°截面對應(yīng)的溫度分布都出現(xiàn)了三峰,且90°和270°截面中間峰更顯著,而0°和225°截面對應(yīng)的溫度分布兩側(cè)峰較中間峰顯著。180°截面位于兩列雙孔之間,所以此截面處的溫度分布呈現(xiàn)了雙峰分布形式。
4)發(fā)動機短艙整體溫度基本高于0℃,作為重點關(guān)注區(qū)域的內(nèi)側(cè)表面溫度較高,且高于10℃,可滿足防冰系統(tǒng)需求。
圖6為內(nèi)外耦合計算收斂后的表面溢流水量結(jié)果。從圖中可以看出:
1)溢流水從駐點位置開始向短艙內(nèi)外兩側(cè)流動,內(nèi)外表面溢流水均未在防冰區(qū)內(nèi)全部蒸發(fā),會流向防冰區(qū)后部且可能形成冰脊;內(nèi)表面溢流水量明顯少于外表面。
2)笛形管噴口正對的短艙周向位置,沿其弦向方向溢流水明顯減少甚至全部蒸發(fā)。
3)在270°截面附近溢流水分布異常,可能是此處短艙與機身的連接件影響了附近流場,從而導(dǎo)致空氣與壁面之間剪切力方向發(fā)生改變,使得溢流水流動方向改變。
圖6 發(fā)動機短艙表面溢流水流量云圖(單位:g/s)Fig.6 Contour of the runback water mass flux on engine nacelle surface(unit:g/s)
液態(tài)水含量、環(huán)境溫度、飛行攻角以及飛行速度等對短艙防冰表面溫度的影響,與它們對飛機機翼的影響類似,本文不做贅述。這里主要考查發(fā)動機進氣流量對短艙表面溫度的影響。選取了如表1狀態(tài)點進行計算和對比分析,所有狀態(tài)下防冰系統(tǒng)熱氣噴口總壓205 214Pa,防冰系統(tǒng)熱氣噴口總溫493K。
圖7 發(fā)動機二維截面處的局部水收集系數(shù)分布Fig.7 Distribution of local water collection coefficient of engine 2D cross section
圖8 發(fā)動機二維截面處的表面溫度分布Fig.8 Distribution of surface tem perature of engine 2D cross section
發(fā)動機進氣流量變化時,首先會影響空氣流場,進而在此基礎(chǔ)上影響水滴流場。進氣道流量越大,發(fā)動機亦能更多地吸取周圍的空氣,從而使駐點向唇口外側(cè)(+s方向)延伸[15]。從圖7可見,發(fā)動機進氣流量增大時,駐點向發(fā)動機外側(cè)表面偏移。從0°到180°截面,局部水收集系數(shù)數(shù)值相比小流量時的增大量越來越小,甚至在180°截面時局部水收集系數(shù)反而減小。這是兩個影響因素作用的結(jié)果:一方面,進氣流量增大會導(dǎo)致更多的水滴撞擊;另一方面180°截面處的發(fā)動機短艙比較厚,且撞擊駐點離幾何駐點(s=0)最遠,使得局部水收集系數(shù)減小。兩者綜合作用導(dǎo)致在180°截面附近,進氣流量更大時局部水收集系數(shù)卻減小。圖8是不同發(fā)動機引氣流量下的溫度結(jié)果。隨著發(fā)動機進氣流量增大,更多氣流流入發(fā)動機進氣道,發(fā)動機內(nèi)側(cè)壁面對流換熱系數(shù)增大,而外側(cè)壁面對流換熱系數(shù)卻有所降低,使得發(fā)動機短艙表面溫度大部分區(qū)域降低,但在外側(cè)表面部分位置的溫度反而略微升高。另外從圖9中所示的發(fā)動機短艙表面溢流水量圖中,也可以明顯地看出駐點的偏移。溢流水在駐點附近會出現(xiàn)谷值,然后沿氣流方向流動;隨著撞擊水的加入,溢流水量逐漸增多直到撞擊極限附近達到最大;然后隨著蒸發(fā)的作用,溢流水量又逐漸減小。90°截面處溢流水分布曲線出現(xiàn)跳躍,原因是溢流水流動方向和此截面不在一個平面。由于撞擊水量增大,表面溫度降低,因此發(fā)動機表面溢流水總體呈上升趨勢。狀態(tài)2流量比較大,發(fā)動機內(nèi)部抽吸空氣能力增強,因此朝發(fā)動機內(nèi)部的剪切力較大,導(dǎo)致某些截面處內(nèi)表面的溢流水量較大。
圖9 發(fā)動機二維截面處的溢流水量分布Fig.9 Distribution of runback water mass flux of engine 2D cross section
對三維發(fā)動機短艙的防冰性能進行了計算和研究。利用CFD方法計算表面換熱系數(shù),以及距離加權(quán)反比插值法對內(nèi)外表面數(shù)據(jù)進行插值,更適用于三維防冰問題的計算。當(dāng)其他飛行氣象條件固定,只考慮發(fā)動機進氣流量的影響時,隨著進氣流量的增大,發(fā)動機短艙蒙皮表面溢流水量增加,蒙皮表面溫度有明顯的下降趨勢。本文的研究方法可為短艙防冰系統(tǒng)性能評估提供幫助。
[1]Messinger B L.Equilibrium temperature of an unheated icing surface as a function of air speed[J].Journal of the Aeronautical Sciences,1953,20(1):29-421.
[2]Wright W B.User manual for the NASA Glenn ice accretion code LEWICE version 2.2.2[R].NASA CR 2002-211793.
[3]Saeed F,Morency F,Paraschivoiu I.Numerical simulation of a hotair anti-icing system[R].AIAA 2000-0630.
[4]Honsek R,Habashi W G.FENSAP-ICE:Eulerian modeling of droplet impingement in the SLD regime of aircraft icing[R].AIAA2006-465.
[5]Bu X Q,Lin G P,Yu J,et al.Numerical simulation of an airfoil electrothermal anti-icing system[J].Proc.IMechE Part G:J.Aerospace Engineering,2013,227(10):1608-1622.DOI:10.1177/0954410 012463525.
[6]Shen X B,Lin G P,Yu J,et al.Three-dimensional numerical simulation of ice accretion at the engine inlet[J].Journal of Aircraft,2013,50(2):635-642.DOI:10.2514/1.C031992.
[7]Chang S N,Yuan M M,Huo X H,et al.Investigations of the bleed air anti-icing system for an aircraft wing[J].Journal of Aerospace Power,2008,23(6):1141-1145.(in Chinese)常士楠,袁美名,霍西恒,等.某型飛機機翼防冰系統(tǒng)計算分析[J].航空動力學(xué)報,2008,23(6):1141-1145.
[8]Bu X Q,Lin G P,Yu J.Three dimensional conjugate heat transfer simulation for the surface temperature of wing hot-air anti-icing system[J].Journal of Aerospace Power.2009,24(11):2495-2500.(in Chinese)卜雪琴,林貴平,郁嘉.三維內(nèi)外熱耦合計算熱氣防冰系統(tǒng)表面溫度[J].航空動力學(xué)報,2009,24(11):2495-2500.
[9]Pan X Y.Performance analysis of micro-ejector anti-icing tunnel[D].Nanjing University of Aeronautics and Astronautics,2006.(in Chinese)潘旭云.某型飛機防冰系統(tǒng)機翼防冰腔性能研究[D].南京航空航天大學(xué),2006.
[10]Huo X H,Wang D W,Li G P,et al.Research of wing anti-ice system performance validation for civil aircraft[J].Civil Aircraft Design&Research,2013,111(4):13-16.(in Chinese)霍西恒,王大偉,李革萍,等.某型飛機機翼防冰系統(tǒng)性能驗證研究[J].民用飛機設(shè)計與研究,2013,111(4):13-16.
[11]Chang S N,Ai S X.Design and calculation for the anti-icing system of an aircraft engine let[J].Journal of Beijing University of Aeronautics and Astronautics,2007,33(6):649-652.(in Chinese)常士楠,艾素霄.飛機發(fā)動機進氣道防冰系統(tǒng)的設(shè)計計算[J].北京航空航天大學(xué)學(xué)報,2007,33(6):649-652.
[12]Zhu Y F,F(xiàn)ang Y F,F(xiàn)eng W C.Design and calculation of aircraft nacelle anti-icing system[J].Journal of Aerospace Power,2012,27(6):1326-1331.(in Chinese)朱永峰,方玉峰,封文春.某型飛機發(fā)動機短艙防冰系統(tǒng)設(shè)計計算[J].航空動力學(xué)報,2012,27(6):1326-1331.
[13]Chen J S.Performance analysis of anti-icer in airplane engine inlet[D].Nanjing University of Aeronautics and Astronautics,2008.(in Chinese)陳景松.發(fā)動機進氣道前緣防冰腔性能研究[D].南京航空航天大學(xué),2008.
[14]Jin G D,Liu Y C,Niu W J.Comparison between inverse distance weighting method and Kriging[J].Journal of Changchun University of Technology,2003,24(3):53-57.(in Chinese)靳國棟,劉衍聰,牛文杰.距離加權(quán)反比插值法和克里金插值法的比較[J].長春工業(yè)大學(xué)學(xué)報,2003,24(3):53-57.
[15]Shen X B,Lin G P,Yang S H.Analysis on three dimensional water droplets impingement characteristics of engine inlet[J].Journal of Beijing University of Aeronautics and Astronautics,2011,37(1): 1-5.(in Chinese)申曉斌,林貴平,楊勝華.三維發(fā)動機進氣道水滴撞擊特性分析[J].北京航空航天大學(xué)學(xué)報,2011,37(1):1-5.
Numerical simulation of the performance of an engine nacelle hot-air anti-icing system
Yu Jia1,2,Zhao Boyang3,Bu Xueqin1,*,Lin Guiping1,Li Zhimao4
(1.School of Aeronautics Science and Engineering,Beihang University,Beijing 100191,China; 2.National Laboratory for Aeronautics and Astronautics,Beijing 100191,China; 3.AVIC The First Aircraft Institute,Xi’an 710089,China; 4.Environment Control System Department,Shanghai Aircraft Design and Research Institute of COMAC,Shanghai 201203,China)
A three-dimensional numerical simulation method with internal-external tight coupled heat transfer is developed to analysis the performance of an engine nacelle hot-air anti-icing system.Effect of engine inlet air mass flow rate on the skin temperature is also analyzed.A three-dimensional Computational Fluid Dynamic(CFD)method is applied to attain the heat transfer coefficient on both side of skin,and an inverse distance weighting algorithm is presented for the interpolation of data exchange between internal and external meshes.The local water collection coefficient on the engine nacelle,the distribution of the surface temperature and the runback water mass flow rate on the surface are obtained to determine.The simulation results indicate that this anti-icing system meets the requirement,the skin temperature decreases and the runback water mass flow rate increases as the engine inlet air mass flow rate rising.
engine;hot-air anti-icing system;surface temperature;engine inlet air mass flux
V211.3;V244.1+5
A
10.7638/kqdlxxb-2015.0212
0258-1825(2016)03-0302-06
2015-12-18;
2016-01-27
郁嘉(1979-),男,上海,講師,研究方向:飛行力學(xué)與飛行安全.E-mail:yujia@buaa.edu.cn
卜雪琴*(1982-),女,江西人,研究方向:飛機防除冰.E-mail:buxueqin@buaa.edu.cn
郁嘉,趙柏陽,卜雪琴,等.某型飛機發(fā)動機短艙熱氣防冰系統(tǒng)性能數(shù)值模擬[J].空氣動力學(xué)學(xué)報,2016,34(3):302-307.
10.7638/kqdlxxb-2015.0212 Yu J,Zhao B Y,Bu X Q,et al.Numerical simulation of the performance of an engine nacelle hot-air anti-icing system[J].Acta Aerodynamica Sinica,2016,34(3):302-307.