錐導(dǎo)乘波體氣動(dòng)外形優(yōu)化與分析

孟希慧，張慶兵

(北京電子工程總體研究所，北京　100854)

?

錐導(dǎo)乘波體氣動(dòng)外形優(yōu)化與分析

孟?；?，張慶兵

(北京電子工程總體研究所，北京100854)

以高超聲速巡航飛行器為應(yīng)用背景,在Ma=6,H=30 km設(shè)計(jì)條件下，對(duì)錐導(dǎo)乘波體進(jìn)行氣動(dòng)外形優(yōu)化設(shè)計(jì)。首先以升阻比為優(yōu)化目標(biāo)，利用遺傳算法對(duì)錐導(dǎo)乘波體進(jìn)行氣動(dòng)力優(yōu)化；然后對(duì)基于氣動(dòng)力優(yōu)化得到的乘波體進(jìn)行前緣鈍化研究，詳細(xì)分析了乘波體前緣的3種鈍化半徑對(duì)其氣動(dòng)力與氣動(dòng)熱的影響。結(jié)果表明，采用遺傳算法對(duì)乘波體工程估算的氣動(dòng)力進(jìn)行優(yōu)化是可靠的。對(duì)乘波體進(jìn)行前緣鈍化可以有效降低最大熱流密度，但同時(shí)也會(huì)降低其升阻比。隨著鈍化半徑的增大，乘波體升阻比降低較為明顯，但對(duì)熱流密度的影響逐漸減弱，因此將乘波體應(yīng)用于高超聲速巡航飛行器時(shí)應(yīng)綜合考慮鈍化對(duì)其氣動(dòng)力和氣動(dòng)熱的影響，尋找最佳平衡點(diǎn)。

高超聲速；乘波體；優(yōu)化設(shè)計(jì)；鈍化；氣動(dòng)力；氣動(dòng)熱

0　引言

高超聲速巡航飛行器需要具備射程長(zhǎng)、機(jī)動(dòng)性強(qiáng)的能力，所以需要使用高升阻比的氣動(dòng)外形。傳統(tǒng)氣動(dòng)外形會(huì)隨著馬赫數(shù)的提高遇到升阻比屏障[1]，而乘波體是追求高升阻比、突破高超聲速飛行器“升阻比屏障”的一種有效手段。

國(guó)內(nèi)外針對(duì)乘波體的生成方法、乘波體的氣動(dòng)力優(yōu)化研究較多。在生成方法研究方面，1959年英國(guó)的Nonweiler[2]提出利用楔形流場(chǎng)構(gòu)造“∧”形乘波體，被認(rèn)為是最早的乘波體，但該構(gòu)型容積率較小，且存在較多不利于氣動(dòng)熱防護(hù)的銳邊緣，因此實(shí)用價(jià)值不高。1980年Rasmussen等[3]提出基于錐形流場(chǎng)生成的乘波體，與“Λ”形乘波體相比，錐導(dǎo)乘波體顯著提高了容積率和升阻比，但是這種構(gòu)型的發(fā)動(dòng)機(jī)進(jìn)口截面流動(dòng)存在橫向流動(dòng)，會(huì)對(duì)吸氣式發(fā)動(dòng)機(jī)的工作造成不利影響。Takashima 和 Lewis等[4]提出乘波體設(shè)計(jì)的楔錐法以及 Sobieczky[5]提出的吻切錐理論較好地解決了這個(gè)問(wèn)題。在氣動(dòng)力優(yōu)化方法研究方面，Rasmussen[6]等采用高超聲速小擾動(dòng)理論對(duì)錐導(dǎo)乘波體進(jìn)行了無(wú)黏氣動(dòng)力優(yōu)化。1987年，Boweutt等[7]將黏性效應(yīng)引入乘波體的優(yōu)化過(guò)程中，開(kāi)展了軸對(duì)稱體繞流的乘波體黏性優(yōu)化研究。而乘波體在高超聲速飛行條件下氣動(dòng)加熱相當(dāng)嚴(yán)重，可能造成乘波體結(jié)構(gòu)的破壞失效，所以針對(duì)乘波體進(jìn)行氣動(dòng)力熱方面的研究是有必要的。

本文針對(duì)Ma=6，H=30 km設(shè)計(jì)條件下的乘波體進(jìn)行氣動(dòng)力熱研究，分析了3種鈍化半徑對(duì)乘波體氣動(dòng)力熱的影響，得到了隨著鈍化半徑的增大，乘波體升阻比降低較為明顯，但對(duì)熱流密度的影響逐漸減弱的結(jié)論，研究結(jié)果對(duì)未來(lái)乘波體設(shè)計(jì)具有較強(qiáng)的工程應(yīng)用價(jià)值。

1　錐導(dǎo)乘波體的生成與氣動(dòng)力優(yōu)化

1.1錐導(dǎo)乘波體的生成

基于錐型流場(chǎng)理論，在設(shè)計(jì)條件下(如表1所示)生成錐導(dǎo)乘波體，生成原理如圖1所示，具體生成過(guò)程參照文獻(xiàn)[8]。

1.2氣動(dòng)力的工程估算

升阻比是衡量飛行器氣動(dòng)性能的一個(gè)重要性能參數(shù)，也是乘波體區(qū)別于其他構(gòu)型的重要性能指標(biāo)。

表1　來(lái)流條件

圖1　錐導(dǎo)乘波體生成原理Fig.1　Sketch of cone-derived waverider configuration design

利用Matlab軟件編程得到乘波體表面各面元的壓力系數(shù)和摩擦力系數(shù)，從而通過(guò)數(shù)值積分來(lái)求得乘波體的升力系數(shù)、阻力系數(shù)和升阻比，具體計(jì)算過(guò)程參照文獻(xiàn)[8-9]。

1.2.1無(wú)黏氣動(dòng)力計(jì)算

由于乘波體的上表面被設(shè)計(jì)為平行于來(lái)流方向的自由流面，壓強(qiáng)系數(shù)為0，所以不提供升力，只需考慮下表面的壓強(qiáng)系數(shù)cp即可。將下表面劃分為三角微元后，根據(jù)各點(diǎn)氣動(dòng)參數(shù)便可利用式(1)求得壓強(qiáng)系數(shù)cp，從而利用式(2)和(3)求得無(wú)黏升力系數(shù)cLP和波阻系數(shù)cDP。

(1)

(2)

(3)

式中：Ayn為各面元面積在Oxz平面上的投影;Azn為各面元面積在Oxy平面上的投影Sp為乘波體浸濕面積。

1.2.2有黏氣動(dòng)力計(jì)算

對(duì)于真實(shí)流動(dòng)需要考慮黏性，利用米多爾-斯馬特參考溫度法可以確定出各面元的摩擦力系數(shù)cf，從而利用式(7)和(8)求得由摩擦力引起的升力系數(shù)cLf和阻力系數(shù)cDf。

高超聲速層流到湍流的轉(zhuǎn)捩問(wèn)題還沒(méi)有成熟的解決方法，因此本文假定乘波體的全表面均為充分發(fā)展的湍流。

米多爾-斯馬特參考溫度法給出平板上的湍流局部摩擦系數(shù)為

(4)

(5)

(6)

從而可求得各微元由摩擦力引起的升力系數(shù)cLf和阻力系數(shù)cDf：

(7)

(8)

式中：Avyn為各面元面積在微元y方向速度上的投影；Avzn為各面元面積在微元z方向速度上的投影。

1.3錐導(dǎo)乘波體的氣動(dòng)力優(yōu)化

1.3.1優(yōu)化方法

遺傳算法是一種現(xiàn)代智能全局尋優(yōu)算法，其原理是效仿生物界中的“物競(jìng)天擇，適者生存”。遺傳算法可以反復(fù)修改對(duì)于個(gè)體解決方案的種群，在每一步中隨機(jī)地從當(dāng)前種群中選擇若干個(gè)體作為父代，使用它們產(chǎn)生子代，在連續(xù)若干代后，種群朝著優(yōu)化解的方向進(jìn)化。該算法不受搜索空間限制性條件(如可微、連續(xù)、單峰)的制約，不需要其他輔助性信息(如導(dǎo)數(shù))，這些特點(diǎn)使得遺傳算法已成功地應(yīng)用到難以用傳統(tǒng)方法求解的復(fù)雜優(yōu)化問(wèn)題之中[10-11]。

P=(x2,x3,x4,x5,y1,y2,y3,y4).

(9)

未優(yōu)化的乘波體前緣線在底面內(nèi)(Oxy平面)的投影如圖2所示。

圖2　未優(yōu)化的乘波體前緣線在底面內(nèi)的投影Fig.2　Projection of the leading edge of non-optimized waverider

為了減少數(shù)學(xué)模型的自由度，提高優(yōu)化速率，將5個(gè)點(diǎn)取為等間距，所以有

(10)

描述前緣線的數(shù)學(xué)模型簡(jiǎn)化為

P=(x5,y1,y2,y3,y4).

(11)

本文針對(duì)激波角為12°，在設(shè)計(jì)條件下，以升阻比最大作為優(yōu)化目標(biāo)，利用Matlab軟件中的遺傳算法工具箱對(duì)乘波體的外形進(jìn)行優(yōu)化。

1.3.2優(yōu)化結(jié)果

圖3給出Ma=6,H=30 km,β=12°條件下乘波體的優(yōu)化收斂曲線，可以看出結(jié)果在第51代基本收斂，升阻比達(dá)到最大值5.874 61，與第1代相比提升了約14%，最大升阻比對(duì)應(yīng)的最優(yōu)個(gè)體值為xbest=(0.830 9,0.416 9,0.466 0,0.570 5,0.670 2)，因此采用遺傳算法對(duì)乘波體工程估算的升阻比進(jìn)行優(yōu)化是可行的。

圖3　Ma=6,H=30 km,β=12°條件下乘波體的優(yōu)化收斂曲線Fig.3　Optimization history at Ma=6, H=30 km,β=12°for waverider

從優(yōu)化得到的前緣線出發(fā)，利用流線追蹤技術(shù)得到乘波體下表面的點(diǎn)云坐標(biāo)，將點(diǎn)云輸入Pro/E工程軟件進(jìn)行曲面擬合，得到下表面，而乘波體上表面為自由來(lái)流流面，便可獲得最優(yōu)乘波體的三維模型, 如圖4所示，其外形參數(shù)見(jiàn)表2。

圖4　最優(yōu)乘波體三維圖(Ma=6,H=30 km)Fig.4　View of optimum waverider(Ma=6,H=30 km)

表2Ma=6,H=30 km,β=12°時(shí)最優(yōu)乘波體的外形參數(shù)

Table 2　Geometry parameters of optimum waverider at Ma=6,H=30 km, β=12°

2　優(yōu)化后的乘波體氣動(dòng)力數(shù)值模擬

對(duì)于上述工程算法,采用數(shù)值計(jì)算方法對(duì)Ma=6,H=30 km,β=12°條件下最優(yōu)乘波體的氣動(dòng)力性能參數(shù)進(jìn)行驗(yàn)證，對(duì)比結(jié)果如表3所示。

表3　Ma=6,H=30 km,β=12°條件下最優(yōu)乘波體氣動(dòng)力性能參數(shù)

由上表可以看出，數(shù)值模擬的計(jì)算結(jié)果與優(yōu)化過(guò)程中采用的工程估算方法能夠較好吻合，3組氣動(dòng)力參數(shù)的最大誤差在5%之內(nèi)。

3　前緣鈍化對(duì)乘波體氣動(dòng)力熱影響分析

3.1計(jì)算模型的前緣鈍化

在用于計(jì)算氣動(dòng)力性能的乘波體模型中,其前緣是趨近于無(wú)限薄的，但在實(shí)際應(yīng)用中,無(wú)論是從加工的角度還是從熱防護(hù)的角度，乘波體前緣都應(yīng)是進(jìn)行一定鈍化的。

本文采用的鈍化方式為保持乘波體下表面不變，將上表面向上偏移一段距離(偏移距離為鈍化直徑)，之后將上下表面采用圓弧過(guò)渡。采用這樣的鈍化方法優(yōu)點(diǎn)在于可以盡量減少原有前緣線的位置變化，以減少其對(duì)于乘波體氣動(dòng)力性能的影響[12-13]，鈍化后的乘波體如圖5所示。

圖5　鈍化后的乘波體Fig.5　Waverider with blunt leading edge

3.2計(jì)算格式與網(wǎng)格

根據(jù)文獻(xiàn)[14-15]中的結(jié)論，Roe格式在熱流計(jì)算精確性方面具有優(yōu)勢(shì)；與格式效應(yīng)相比，網(wǎng)格對(duì)熱流影響更大，而且駐點(diǎn)附近切向網(wǎng)格尺度對(duì)熱流計(jì)算結(jié)果影響甚微，法向網(wǎng)格方面僅有物面第1層網(wǎng)格高度對(duì)熱流影響巨大。本文計(jì)算中壁面附近網(wǎng)格高度達(dá)到1e-5 m，如圖6所示。

圖6　乘波體頭部壁面局部網(wǎng)格Fig.6　Local mesh of the nose of waverider

3.3計(jì)算結(jié)果與分析

Vanmol在文獻(xiàn)[16]中提出對(duì)于主動(dòng)冷卻來(lái)說(shuō)，前緣半徑最少需要 1 cm。利用上述鈍化方法，對(duì)乘波體分別以1，1.5，2 cm 3種鈍化半徑值進(jìn)行鈍化，在設(shè)計(jì)條件下進(jìn)行數(shù)值計(jì)算，得到氣動(dòng)力參數(shù)見(jiàn)表4。

表4　3種鈍化半徑對(duì)應(yīng)的氣動(dòng)力參數(shù)

表4為利用CFD數(shù)值模擬得到的3種鈍化半徑后的氣動(dòng)力參數(shù)，可以看出，鈍化后乘波體的阻力系數(shù)有所上升，而升力系數(shù)幾乎不變，導(dǎo)致升阻比隨鈍化半徑增大而變小。但對(duì)于這里所討論的長(zhǎng)度為3.2 m的乘波體，考慮到防熱需求的前緣鈍化半徑并沒(méi)有從本質(zhì)上降低乘波體的氣動(dòng)力性能，仍具有較高升阻比。

鈍化半徑為1 cm時(shí)乘波體升阻比為4.721 8，鈍化半徑為1.5 cm時(shí)乘波體升阻比為4.548 2，與1 cm鈍化相比降低了3.68%；當(dāng)鈍化半徑為2 cm時(shí)乘波體升阻比為4.126 0，與1.5 cm鈍化相比降低了9.28%。從圖7可以看出，隨著鈍化半徑的增大，前緣鈍化對(duì)于升阻比的負(fù)面影響越來(lái)越強(qiáng)。

圖7　升阻比隨鈍化半徑變化趨勢(shì)Fig.7　L/D at different blunt radius

對(duì)比鈍化前后對(duì)稱面處壓力云圖如圖8，9所示，可以看出未鈍化時(shí)，激波基本附著于前緣上；鈍化半徑為2 cm時(shí)，前緣激波脫體，局部形成弓形激波，導(dǎo)致乘波體受到的波阻增加。將未鈍化和鈍化半徑為2 cm的波阻系數(shù)和摩阻系數(shù)結(jié)果列表，如表5所示，發(fā)現(xiàn)鈍化后阻力系數(shù)顯著增大主要是因?yàn)椴ㄗ柘禂?shù)顯著增大，而摩擦阻力系數(shù)基本不變。

表5　鈍化前后的壓差阻力系數(shù)與摩擦阻力系數(shù)

圖8　未鈍化時(shí)對(duì)稱面處壓力云圖Fig.8　Contours of static pressure with non-blunt leading edge

圖9　鈍化半徑為2 cm時(shí)對(duì)稱面處壓力云圖Fig.9　Pressure contours at symmetry plane with 2 cm blunt radius

采用Roe計(jì)算格式進(jìn)行氣動(dòng)熱計(jì)算，為了說(shuō)明利用CFD數(shù)值仿真方法計(jì)算氣動(dòng)熱的效果，本文結(jié)合半球模型進(jìn)行驗(yàn)證。將半球模型的CFD結(jié)果與精度較高的Kemp-Riddle駐點(diǎn)熱流計(jì)算公式[17](式(12))計(jì)算結(jié)果進(jìn)行對(duì)比，如表6所示。

(12)

表6　半球模型3種半徑對(duì)應(yīng)的駐點(diǎn)熱流

由表6可以看出，半球模型CFD的計(jì)算結(jié)果與Kemp-Riddle駐點(diǎn)熱流計(jì)算公式的結(jié)果相差較小，3種半徑對(duì)應(yīng)的結(jié)果最大相差10.4%，因此可以認(rèn)為利用CFD數(shù)值模擬計(jì)算熱流的結(jié)果是可靠的。利用CFD數(shù)值模擬方法計(jì)算乘波體外形的駐點(diǎn)熱流結(jié)果如表7所示，最大熱流密度隨半徑變化趨勢(shì)如圖10所示。

表7　3種鈍化半徑對(duì)應(yīng)的駐點(diǎn)熱流

鈍化半徑為1 cm時(shí)乘波體最大熱流密度為387.081 1 kW/m2，鈍化半徑為1.5 cm時(shí)乘波體最大熱流密度降至335.960 4 kW/m2，與1 cm鈍化相比降低了13.21%；當(dāng)鈍化半徑增大到2 cm時(shí)，乘波體最大熱流密度降至306.621 5 kW/m2，與1.5 cm鈍化相比降低了8.73%。從圖10可以看出，隨著鈍化半徑的增大，前緣鈍化對(duì)于最大熱流密度的影響越來(lái)越小。

圖10　最大熱流密度隨鈍化半徑變化趨勢(shì)Fig.10　qmax at different blunt radiuses

4　結(jié)論

本文在Ma=6,H=30 km設(shè)計(jì)條件下，針對(duì)錐導(dǎo)乘波體進(jìn)行氣動(dòng)外形優(yōu)化研究，詳細(xì)分析了3種前緣鈍化半徑對(duì)氣動(dòng)力熱的影響，主要研究結(jié)論如下：

(1) 采用遺傳算法對(duì)乘波體進(jìn)行氣動(dòng)力優(yōu)化，并將采用遺傳算法優(yōu)化后得到的工程估算結(jié)果與CFD數(shù)值模擬結(jié)果進(jìn)行對(duì)比，發(fā)現(xiàn)2種方法的計(jì)算結(jié)果能夠較好地吻合，表明采用遺傳算法對(duì)乘波體工程估算的氣動(dòng)力進(jìn)行優(yōu)化是可靠的。

(2) 對(duì)乘波體進(jìn)行前緣鈍化可以有效降低最大熱流密度，但同時(shí)也會(huì)降低其升阻比。

(3) 隨著乘波體前緣鈍化半徑的增大，鈍化對(duì)升阻比的影響越來(lái)越強(qiáng)，而對(duì)最大熱流密度的影響越來(lái)越弱。因此在進(jìn)行乘波體設(shè)計(jì)時(shí)，應(yīng)綜合考慮氣動(dòng)力和氣動(dòng)熱的設(shè)計(jì)指標(biāo)，選擇合理的前緣鈍化半徑，從而達(dá)到既能有效地降低熱流密度又能保證高升阻比的目的。

[1]KUCHEMANN D. The Aerodynamic Design of Aircraft[M].Oxford: Pergamon Press, 1978.

[2]NOWEILER T R F. Aerodynamic Problems of Manned Space Vehicles [J]. Journal of the Royal Aeronautical Society, 1963,67(9):521-528.

[3]RASMUSSEN M L. Waverider Configurations Derived from Inclined Circular and Elliptic Cones[J]. Journal of Spacecraft and Rocket，1980,17(6):537-545.

[4]TAKASHIMA N, LEWIS M J. Waverider Configurations Based on Nonaxisymmetric Flow Fields for Engine-AirFrame Integration[R]. AIAA-94-0380, 1994.

[5]SOBIECZKY H. Hypersonic Waverider Design from Given Shock Waves[C]∥ Proceedings of the First International Hypersonic Waverider Symposium. Maryland University of Maryland, 1990.

[6]KIM B S, RASMUSSEN M L, JISEHKE M C. Optimization of Waverider Configurations Generated from Axisymmetric Conical Flows[J]. Journal of Spacecraft and Rocket，1983,20(5):461-469.

[7]BOWCUTT K G, ANDERSON J D. Viscous Optimized Waveriders[R]. AIAA-87-0272，1987.

[8]楊海江.乘波體氣動(dòng)外形設(shè)計(jì)與計(jì)算[D].南京：南京航空航天大學(xué)，2008.

YANG Hai-jiang. Waverider Aerodynamic Configuration Design and Aerodynamic Performance Calculation [D]. Nanjing: Nanjing University of Aeronautics and Astronautics, 2008.

[9]小約翰 D 安德森.高超聲速和高溫氣體動(dòng)力學(xué)[M].楊永，李棟，譯.北京：航空工業(yè)出版社，2013.

John D Anderson Jr. Hypersonic and High-Temperature Gas Dynamics[M].YANG Yong,LI Dong,Translated.Beijing: Aviation Industry Press,2013.

[10]雷英杰，張善文.MATLAB遺傳算法工具箱及應(yīng)用[M].西安:西安電子科技大學(xué)出版社,2004:1-10.

LEI Ying-jie, ZHANG Shan-wen. The Genetic Algorithm Toolbox of MATLAB and Its Application [M].Xi′an:Xidian University Press, 2004:1-10.

[11]史峰,王輝,郁磊.MATLAB智能算法30個(gè)案例分析 [M].北京:北京航空航天大學(xué)出版社, 2011：1-2.

SHI Feng, WANG Hui, YU Lei. 30 Cases of MATLAB Intelligent Algorithm[M].Beijing: Beihang University Press, 2011：1-2.[12]梅東牧,武哲,李天.乘波飛行器的優(yōu)化設(shè)計(jì)和氣動(dòng)熱計(jì)算研究[J]. 航空計(jì)算技術(shù),2008,38(6):15-16.

MEI Dong-mu, WU Zhe, LI Tian. Research of the Waverider Optimization Design and Aerothermo Dynamics[J].Aeronautical Computing Technique, 2008,38(6):15-16.

[13]Wilson Santos. Bluntness Effects on Lift-to-Drag Ratio of Leading Edges for Hypersonic Waverider Configurations[C]∥The 18th AIAA/3AF International Space Planes and Hypersonic Systems and Technologies Conference, France, September 24-28, 2012.

[14]閆超,禹建軍,李君哲.熱流CFD計(jì)算中格式和網(wǎng)格效應(yīng)若干問(wèn)題研究[J].空氣動(dòng)力學(xué)報(bào),2006, 24(1): 125-130.

YAN Chao, YU Jian-jun, LI Jun-zhe. Scheme Effect and Grid Dependency in CFD Computations of Heat Transfer[J]. Acta Aerodynamica Sinica, 2006, 24(1): 125-130.

[15]張建偉.高超聲速再入體氣動(dòng)熱數(shù)值模擬研究[D].濟(jì)南：山東大學(xué)，2012.

ZHANG Jian-wei. Numberical Study of Aeroheating for Hypersonic Reentry Bodies[D]. Jinan:Shandong University,2012.

[16]VANMOL D. Heat Transfer Characteristics of Hypersonic Waveriders with Emphasis on the Leading Edge Effects[D].University of Maryland, College Park, 1991．

[17]姜貴慶，劉連元.高速氣流傳熱與燒蝕熱防護(hù)[M].北京：國(guó)防工業(yè)出版社，2003: 35-36.

JIANG Gui-qing, LIU Lian-yuan. Heat Transfer of Hypersonic Gas and Ablation Thermal Protection[M]. Beijing: National Defence Industry Press,2003: 35-36.

Optimization and Analysis of Cone-Derived Waverider

MENG Xi-hui, ZHANG Qing-bing

(Beijing Institute of Electronic System Engineering, Beijing 100854，China )

The optimization and analysis of cone-derived waverider is completed at the conditions ofMa=6,H=30 km. This design condition corresponds to typical hypersonic cruise conditions. The optimization of aerodynamics is completed with maximum lift over drag ratio as the first optimized objection. The CFD is used to analyze the influence of bluntness on the performance of waverider．The relations of aerodynamics and aerodynamic heating to blunt radius are obtained at three blunt radiuses. The results show that bluntness could reduce the maximum heat flux effectively, but it also reduces the aerodynamic performance. As the blunt radius increases, the aerodynamic performance reduced significantly, but its influence to the heat flux wears drops. When the aerodynamic shape of waverider is designed, the influence of bluntness to aerodynamics and aerodynamic heating should be considered synthetically and a compromise should be made between them.

hypersonic; waverider; optimized design; blunt; aerodynamics; aerodynamic heating

2015-08-18；

2015-09-20

有

孟希慧(1990-)，女，天津人。碩士生，研究方向?yàn)轱w行器設(shè)計(jì)。

通信地址：100854北京142信箱30分箱E-mail：15802952724@163.com

10.3969/j.issn.1009-086x.2016.04.005

TJ761.6; TJ760.1

A

1009-086X(2016)-04-0024-07