金銘君，李 強(qiáng)

(西北工業(yè)大學(xué) 燃燒、熱結(jié)構(gòu)與內(nèi)流場(chǎng)重點(diǎn)實(shí)驗(yàn)室，西安 710072)

沖擊作用下推進(jìn)劑變形的流固耦合分析方法

金銘君，李 強(qiáng)

(西北工業(yè)大學(xué) 燃燒、熱結(jié)構(gòu)與內(nèi)流場(chǎng)重點(diǎn)實(shí)驗(yàn)室，西安 710072)

固體火箭發(fā)動(dòng)機(jī)的點(diǎn)火過(guò)程是一個(gè)復(fù)雜多變的理化過(guò)程，具有時(shí)間短、升溫、升壓梯度大等特點(diǎn)。針對(duì)固體火箭發(fā)動(dòng)機(jī)點(diǎn)火過(guò)程中的裝藥結(jié)構(gòu)完整性問(wèn)題，文中建立了一套用于分析沖擊作用下固體推進(jìn)劑變形現(xiàn)象的仿真模型。采用RANS和ALE方法，分別對(duì)流體域和固體域進(jìn)行求解，以兩場(chǎng)獨(dú)立交叉耦合迭代的模式實(shí)現(xiàn)了仿真過(guò)程。以一個(gè)推進(jìn)劑冷流沖擊實(shí)驗(yàn)作為算例，對(duì)仿真模型進(jìn)行了驗(yàn)證，計(jì)算值與測(cè)量值間誤差不超過(guò)10%，仿真模型計(jì)算可靠，具有向固體火箭發(fā)動(dòng)機(jī)實(shí)際點(diǎn)火過(guò)程拓展的價(jià)值。

固體火箭發(fā)動(dòng)機(jī)；點(diǎn)火過(guò)程；推進(jìn)劑變形；流固耦合

0 引言

固體火箭發(fā)動(dòng)機(jī)的點(diǎn)火過(guò)程是一個(gè)短暫的快速升壓過(guò)程。其中，高溫高壓燃?xì)鈱?duì)固體推進(jìn)劑裝藥的沖擊作用，會(huì)導(dǎo)致裝藥結(jié)構(gòu)的變形、破損和移位，甚至導(dǎo)致發(fā)動(dòng)機(jī)失效爆炸。因此，為確保整個(gè)發(fā)動(dòng)機(jī)的工作性能及可靠性、安全性，對(duì)發(fā)動(dòng)機(jī)點(diǎn)火瞬態(tài)過(guò)程開(kāi)展裝藥結(jié)構(gòu)完整性及瞬態(tài)流暢變化研究是非常有必要的。一直以來(lái)，固體火箭發(fā)動(dòng)機(jī)的研制和定型主要以試驗(yàn)為主，周期長(zhǎng)、耗資大，加之燃燒室內(nèi)工作環(huán)境惡劣，點(diǎn)火過(guò)程又極為短暫，且復(fù)雜多變，對(duì)其進(jìn)行實(shí)驗(yàn)研究是極為困難的。近年來(lái)，隨著各種高精度、高效率計(jì)算方法的提出和完善，以及高速計(jì)算機(jī)的快速發(fā)展，計(jì)算機(jī)數(shù)值仿真技術(shù)憑借其高效、經(jīng)濟(jì)、安全等優(yōu)勢(shì)，在火箭發(fā)動(dòng)機(jī)的研制中，占據(jù)了愈發(fā)重要的地位。

國(guó)內(nèi)外大量研究表明，升壓梯度是破壞裝藥結(jié)構(gòu)完整性的重要因素[1]。在點(diǎn)火階段，過(guò)高的壓強(qiáng)峰或升壓速率會(huì)對(duì)推進(jìn)劑藥柱產(chǎn)生沖擊，造成藥柱變形或移位，而藥柱變形或移位又使得流動(dòng)區(qū)域發(fā)生了變化。對(duì)于這樣一個(gè)典型的流固交互問(wèn)題，采用流固耦合法對(duì)其分析，無(wú)疑是最佳選擇。1982年，Tante G L將流體和結(jié)構(gòu)之間的耦合關(guān)系作為動(dòng)力學(xué)模型，驗(yàn)證了燃燒室和藥柱的相對(duì)運(yùn)動(dòng)，建立了用于分析SRM點(diǎn)火爆炸原因的二維非定常模型。1991年，Titan IV SRMU地面試車失敗，掀起了針對(duì)點(diǎn)火沖擊影響的流固耦合方法研究熱潮。1991年，Johnston采用二維流場(chǎng)耦合結(jié)構(gòu)變形模型，計(jì)算了內(nèi)壓作用下藥柱變形及其對(duì)流場(chǎng)的影響，給出了較為安全的藥柱臨界模數(shù)[2]。此后，許多研究小組開(kāi)展了此方面的研究工作，并得到了有價(jià)值的結(jié)論。國(guó)內(nèi)報(bào)道的主要有海軍航空工程學(xué)院飛行器工程系通過(guò)MPCCI連接FLUENT和ABAQUS軟件，對(duì)SRM快速升壓過(guò)程的研究[1]；以及南京理工大學(xué)機(jī)械工程學(xué)院利用CFX和ANSYS軟件，對(duì)固體推進(jìn)劑裂紋內(nèi)點(diǎn)火過(guò)程的流固耦合研究等[3]。

本文著眼于具有良好適應(yīng)性的高精度SRM點(diǎn)火過(guò)程流固耦合方法研究，提出了一種適用的流固耦合計(jì)算分析模型。由于實(shí)際發(fā)動(dòng)機(jī)點(diǎn)火過(guò)程的推進(jìn)劑應(yīng)變和變形數(shù)據(jù)無(wú)法通過(guò)實(shí)驗(yàn)獲得，因此本文以一個(gè)推進(jìn)劑冷流沖擊實(shí)驗(yàn)作為算例，對(duì)計(jì)算模型進(jìn)行了校驗(yàn)，通過(guò)對(duì)比計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果中測(cè)試段內(nèi)的壓強(qiáng)曲線和推進(jìn)劑變形最大位移的大小，來(lái)對(duì)仿真模型的準(zhǔn)確性做出評(píng)價(jià)。

1 計(jì)算模型

1.1 流動(dòng)計(jì)算模型

針對(duì)氣流對(duì)推進(jìn)劑沖擊過(guò)程的流動(dòng)特性，對(duì)流動(dòng)模型做如下假設(shè)：

(1)計(jì)算中所涉及流體為有粘可壓縮牛頓流體，且可視為各向同性的連續(xù)介質(zhì)；

(2)冷流沖擊過(guò)程中流體溫度脈動(dòng)較小，忽略由溫度脈動(dòng)帶來(lái)的粘性系數(shù)的脈動(dòng)；

(3)流場(chǎng)運(yùn)算中忽略重力影響。

將密度ρ、速度ui、壓力p及溫度T經(jīng)上式處理后，代入有粘可壓縮N-S方程組，即可得到可壓縮雷諾控制方程組，進(jìn)行必要簡(jiǎn)化后，有[8]

(1)

(2)

(3)

1.2 湍流模型

本文采用渦粘性方法中的一方程模型：Spalart-Allmaras模型(S-A模型)來(lái)求解雷諾平均方程。

對(duì)于一方程模型的一般形式，有

(4)

式中SP為生成項(xiàng)；SD為耗散項(xiàng)；D為擴(kuò)散項(xiàng)。

(5)

其中

這里，d為到物面的最近距離。對(duì)于式中其他參變量的關(guān)系，有

g=r+Cw2(r6-r)

上述各式中，常系數(shù)取值如表1所示。

表1 S-A模型中常系數(shù)取值

1.3 推進(jìn)劑結(jié)構(gòu)力學(xué)計(jì)算模型

對(duì)推進(jìn)劑結(jié)構(gòu)模型做如下假設(shè)：

(1)認(rèn)為沖擊時(shí)間短暫，固體推進(jìn)劑的粘彈特性不能完全體現(xiàn)，將其近似為大變形彈性模型；

(2)將推進(jìn)劑材料視為均質(zhì)的各向同性Kirchhoff材料，忽略重力影響。

對(duì)于上述假設(shè)，材料本構(gòu)關(guān)系[9-10]滿足：

Sij=Cijk1Ek1=λEkkδij+2μEij

(6)

其中，Sij為基爾霍夫應(yīng)力；Eij為格林應(yīng)變；Cijk1為彈性模量的四階張量；λ和μ為拉梅(Lamé)常數(shù)，可表示為

λ=υE/(1+υ)(1-2υ)

μ=E/2(1+υ)

1.4 非線性有限元方法

采用ALE(Arbitrary Lagrangian-Eulerian)有限元方法，對(duì)結(jié)構(gòu)域進(jìn)行離散。在ALE方法中，由于質(zhì)量守恒被強(qiáng)加作為偏微分方程，所以需要建立一種弱形式[10]。將ALE域劃分為單元，定義所有非獨(dú)立變量為單元坐標(biāo)的函數(shù)，對(duì)于單元e，其ALE坐標(biāo)為

χ(ξe)=Φe(ξe)=χ1N1(ξe)

其中，ξe為單元e的坐標(biāo)。對(duì)于網(wǎng)格運(yùn)動(dòng)，有

其中，x1(t)為節(jié)點(diǎn)的運(yùn)動(dòng)。對(duì)于網(wǎng)格速度，有

對(duì)于有限元矩陣方程，有連續(xù)方程：

其中,Mρ、Lρ、Kρ分別為容量、轉(zhuǎn)換和散度矩陣：

對(duì)于動(dòng)量方程，有：

其中M和L分別是廣義質(zhì)量和傳遞矩陣，對(duì)應(yīng)于在參考構(gòu)型描述下的速度f(wàn)int和fext分別為內(nèi)力和外力向量：

1.5 耦合計(jì)算

目前，流固耦合問(wèn)題的求解方式有兩類：兩場(chǎng)交叉迭代和直接全場(chǎng)同時(shí)求解。后者計(jì)算量巨大，且具體實(shí)現(xiàn)十分復(fù)雜，故本文采用兩場(chǎng)分開(kāi)進(jìn)行交叉迭代計(jì)算，各自單獨(dú)使用獨(dú)立解算器和網(wǎng)格進(jìn)行解算的方式。在兩場(chǎng)數(shù)據(jù)交換的處理上，一般是依靠在兩套非匹配網(wǎng)格間進(jìn)行雙向插值來(lái)完成的，其精度受兩插值網(wǎng)格間重合度等諸多因素影響。本文采用在交界面處通過(guò)插值重構(gòu)一層虛網(wǎng)格進(jìn)行數(shù)據(jù)傳遞和相關(guān)變量解算的耦合方法，數(shù)據(jù)通過(guò)最小二乘法在兩場(chǎng)和虛擬網(wǎng)格間進(jìn)行傳遞，減少了來(lái)自雙向插值帶來(lái)的誤差；交界面節(jié)點(diǎn)力可在虛網(wǎng)格上直接進(jìn)行計(jì)算，節(jié)省了計(jì)算時(shí)間和計(jì)算量。設(shè)上角標(biāo)nf為流體解算器中的時(shí)間步，ns則為結(jié)構(gòu)解算器中的時(shí)間步；下標(biāo)f代表流體解算器中的變量，s代表結(jié)構(gòu)解算器中的變量，f-s則代表交界面處的變量。

通常固體膨脹波速cd比流體中波速c大，出于對(duì)算法穩(wěn)定性的考慮，固體解算器中的時(shí)間步長(zhǎng)ts通常小于流體解算器中的時(shí)間步長(zhǎng)tf。因此，兩場(chǎng)求解器中分別取其最佳時(shí)間步長(zhǎng)進(jìn)行交叉迭代，其具體流程和數(shù)據(jù)傳遞如圖1和圖2所示。

2 實(shí)驗(yàn)

2.1 實(shí)驗(yàn)方法

該實(shí)驗(yàn)采用一套冷流來(lái)流模擬裝置來(lái)模擬發(fā)動(dòng)機(jī)點(diǎn)火過(guò)程中壓強(qiáng)迅速升高的燃?xì)饬?，?duì)縮比的推進(jìn)劑進(jìn)行沖擊實(shí)驗(yàn)。實(shí)驗(yàn)開(kāi)始前，先在高速攝影儀驅(qū)動(dòng)程序界面上對(duì)推進(jìn)劑試件取一系列距離已知的標(biāo)定點(diǎn)；實(shí)驗(yàn)時(shí)，通過(guò)高速攝影儀來(lái)獲得實(shí)驗(yàn)過(guò)程的影像資料，并在測(cè)試段首尾兩端配置壓力傳感器來(lái)對(duì)內(nèi)腔壓強(qiáng)進(jìn)行監(jiān)測(cè)；實(shí)驗(yàn)結(jié)束后，對(duì)影像資料中標(biāo)定點(diǎn)間的距離變化情況進(jìn)行分析，獲得推進(jìn)劑形變、形變速度隨時(shí)間變化的曲線[11-13]。

2.2 實(shí)驗(yàn)裝置

實(shí)驗(yàn)裝置主要由3部分組成：氮?dú)夂銐汗庀到y(tǒng)，推進(jìn)劑變形測(cè)量段及數(shù)據(jù)采集系統(tǒng)。氮?dú)夂銐汗庀到y(tǒng)由4個(gè)氮?dú)鈿馄亢头€(wěn)壓裝置組成，可提供2～4.5 MPa的沖擊壓強(qiáng)。實(shí)驗(yàn)測(cè)量段內(nèi)腔尺寸為240 mm×200 mm×140 mm，內(nèi)有固定槽，可保證在測(cè)量過(guò)程中試件沒(méi)有移動(dòng)；前端為30 mm厚的有機(jī)玻璃，在保證測(cè)量段密閉性的同時(shí)，又兼顧了圖像采集系統(tǒng)對(duì)視野清晰度的要求；實(shí)驗(yàn)采用兩路測(cè)壓，測(cè)壓孔配置在測(cè)量段后蓋板上，測(cè)量點(diǎn)P1臨近氣流入口處，測(cè)量點(diǎn)P2臨近氣流出口位置。數(shù)據(jù)采集系統(tǒng)分為圖像采集系統(tǒng)和壓強(qiáng)測(cè)量系統(tǒng)兩部分。圖像采集部分采用美國(guó)Phanton高速攝影儀(V4.3)，所采用的傳感器分辨率為800×600，采樣率為300幀；壓強(qiáng)測(cè)量部分采用Dewetron并行數(shù)據(jù)采集系統(tǒng)。

2.3 實(shí)驗(yàn)方案

實(shí)驗(yàn)對(duì)不同沖擊壓強(qiáng)下推進(jìn)劑的變形情況進(jìn)行了測(cè)量，為了減小實(shí)驗(yàn)誤差，對(duì)每個(gè)工況進(jìn)行了3次實(shí)驗(yàn)，最終結(jié)果為3次實(shí)驗(yàn)結(jié)果的平均值。本文采用前文所述計(jì)算模型，對(duì)其中沖擊壓強(qiáng)為3.3、3.6、4.1 MPa的工況開(kāi)展了數(shù)值計(jì)算工作。

3 結(jié)果對(duì)比

3.1 實(shí)驗(yàn)結(jié)果

圖4給出了在沖擊壓強(qiáng)為3.3、3.6、4.1 MPa的工況下，氣流入口和出口處壓強(qiáng)隨時(shí)間變化曲線。其中，曲線p1代表氣流入口處壓強(qiáng)，p2代表氣流出口處壓強(qiáng)。

在各個(gè)工況下，沖擊壓強(qiáng)達(dá)到最大值需要約0.12 s，與實(shí)際發(fā)動(dòng)機(jī)點(diǎn)火過(guò)程的壓強(qiáng)建立時(shí)間大致相當(dāng)，可認(rèn)為該冷流推進(jìn)劑沖擊實(shí)驗(yàn)具備一定對(duì)真實(shí)發(fā)動(dòng)機(jī)點(diǎn)火過(guò)程的參照性。

在各個(gè)工況下，氣流入口與出口處存在明顯的壓強(qiáng)差，分析認(rèn)為這是造成推進(jìn)劑變形的主要原因之一。推進(jìn)劑形變的形狀與氮?dú)鉀_擊的位置有關(guān)，在該實(shí)驗(yàn)中，沖擊位置位于推進(jìn)劑的右側(cè)，推進(jìn)劑受沖擊后向左凹陷，位移最大處為氣流沖擊位置的中央處，如圖5所示。各工況下入-出口壓強(qiáng)差、推進(jìn)劑最大位移及形變速度如表2所示。

沖擊壓強(qiáng)/MPa前后壓差/MPa推進(jìn)劑形變大小/mm推進(jìn)劑形變速率/(mm/s)3．30．1610．70846．9363．60．2360．81753．1304．10．3291．19259．792

3.2 計(jì)算結(jié)果

參照上述實(shí)驗(yàn)，依據(jù)其設(shè)計(jì)參數(shù)配置物理模型和邊界條件，開(kāi)展了對(duì)推進(jìn)劑冷流沖擊的數(shù)值計(jì)算工作。流動(dòng)求解器采用了有限體積方法，空間對(duì)流項(xiàng)的離散采用了HLLC格式，粘性項(xiàng)差分采用WENO格式，時(shí)間項(xiàng)上采用四階Runge-Kutta積分方法。在計(jì)算中，參照實(shí)驗(yàn)測(cè)量到的推進(jìn)劑位移最大位置，設(shè)置了位移監(jiān)測(cè)點(diǎn)。由于實(shí)驗(yàn)結(jié)果分析中只給出了氮?dú)鉀_擊造成的推進(jìn)劑變形的最大位移結(jié)果，而未能給出具體位移情況隨時(shí)間變化的過(guò)程，考慮到氮?dú)獾臎_擊作用只是一個(gè)脈沖過(guò)程，因此在數(shù)值計(jì)算中，在氣流入口位置給定一個(gè)單脈沖過(guò)程的沖擊壓力作為初始邊界條件。計(jì)算域分為16個(gè)分區(qū)，流體計(jì)算網(wǎng)格總數(shù)約為28萬(wàn)，其計(jì)算結(jié)果如下：

首先，由表3可知，各工況下最大位移的計(jì)算值和測(cè)量值間誤差不超過(guò)10%，可認(rèn)為計(jì)算結(jié)果較為準(zhǔn)確。其次，在圖6中給出的推進(jìn)劑最大位移隨時(shí)間變化曲線可看到，在沖擊過(guò)后最初的0.02 s內(nèi)，3個(gè)工況下推進(jìn)劑均沒(méi)有發(fā)生明顯的形變。此后，推進(jìn)劑開(kāi)始變形，且形變量逐漸變大。對(duì)比圖4可知，在沖擊最初的0.02 s內(nèi)，雖然氮?dú)鈿庠春蛯?shí)驗(yàn)器內(nèi)的絕對(duì)壓差很大，但實(shí)驗(yàn)器內(nèi)氣流入口和出口間的相對(duì)壓差卻較小，因此推進(jìn)劑的變形很不明顯；而在0.02 s之后，隨氣流入口和出口間相對(duì)壓差的逐漸增大，推進(jìn)劑變形逐漸明顯，且在壓差達(dá)到最大值的時(shí)刻，推進(jìn)劑的位移變化幾乎同時(shí)達(dá)到最大值。最后，在沖擊過(guò)后的降壓階段，由于計(jì)算邊界無(wú)法準(zhǔn)確模擬氮?dú)鈮簭?qiáng)衰減的過(guò)程，可認(rèn)為此后的計(jì)算結(jié)果沒(méi)有意義。

沖擊壓強(qiáng)/MPa推進(jìn)劑最大位移/mm計(jì)算值實(shí)驗(yàn)值3．30．7210．7083．60．8920．8174．11．0861．192

由于實(shí)驗(yàn)條件所限，實(shí)驗(yàn)結(jié)果并沒(méi)有給出推進(jìn)劑變形大小隨沖擊壓強(qiáng)漸變的具體過(guò)程，因此在計(jì)算結(jié)果中，只有推進(jìn)劑形變的最大值和實(shí)驗(yàn)結(jié)果具有比照作用，但其得到的過(guò)程值可認(rèn)為是具有一定預(yù)驗(yàn)價(jià)值的。通過(guò)對(duì)比計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果中3個(gè)工況下的推進(jìn)劑變形最大位移以及測(cè)試段內(nèi)升壓過(guò)程的壓強(qiáng)曲線，驗(yàn)證了仿真模型的可行性和準(zhǔn)確性。

4 結(jié)束語(yǔ)

(1)為了對(duì)氣流沖擊作用下推進(jìn)劑的變形現(xiàn)象開(kāi)展數(shù)值仿真研究，本文提出了一套基于RANS和ALE方法的流固耦合仿真模型。

(2)通過(guò)在交界面構(gòu)造虛擬網(wǎng)格的方法，實(shí)現(xiàn)了兩場(chǎng)間的數(shù)據(jù)交互，既減少了來(lái)自雙向插值帶來(lái)的誤差，又節(jié)省了計(jì)算時(shí)間和計(jì)算量。

(3)以推進(jìn)劑冷流沖擊實(shí)驗(yàn)作為算例，對(duì)仿真模型進(jìn)行了校驗(yàn)，計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果誤差不超過(guò)10%，驗(yàn)證了仿真模型的可靠性和準(zhǔn)確性。

[1] 于勝春, 趙汝巖, 許濤, 等.固體火箭發(fā)動(dòng)機(jī)快速升壓過(guò)程的流固耦合分析[J].固體火箭技術(shù), 2008, 31(3):232-235.

[2] Johnston W A. Flow-structural analysis of the Ariane-5 solid rocket motor during ignition transient[R]. AIAA 2011-4056.

[3] 韓波, 周長(zhǎng)省, 陳雄.固體推進(jìn)劑裂紋內(nèi)點(diǎn)火過(guò)程流固耦合數(shù)值仿真[J].固體火箭技術(shù), 2011, 34(2):180-183.

[4] 閻超.計(jì)算流體力學(xué)方法及應(yīng)用[M].北京: 北京航空航天大學(xué)出版社, 2006.

[5] Mark Salita. Modern SRM ignition transient modeling(part 1): introduction and physical models[R]. AIAA 2001-3443.

[6] 王新月, 胡春波, 張堃元, 等.氣體動(dòng)力學(xué)基礎(chǔ)[M].西安: 西北工業(yè)大學(xué)出版社, 2006.

[7] David C Wilcox. Turbulence modeling for CFD[M]. California: DCW Industries, Inc., 1993.

[8] 陳矛章.粘性流體動(dòng)力學(xué)基礎(chǔ)[M].北京: 高等教育出版社, 1993.

[9] Toro E F. Riemann solvers and numerical methods for fluid dynamics: a practical introduction[M]. Berlin:Springer-Verlag Berlin and Heidelberg GmbH & Co. K, 1999.

[10] Ted Belytschko, Wing Kam Liu, Brian Moran. Nonlinear finite elements for continua and structures[M]. London: John Wiley & Sons Ltd., 1998.

[11] Soares D Jr, Rodrigues G G, Gonc alves K A. An efficient multi-time-step implicit-explicit method to analyze solid-fluid coupled systems discretized by unconditionally stable time-domain finite element procedures[J]. Computers and Structures, 2010, 88(5-6):387-398.

[12] Luciano Garelli, Rodrigo R Paz, Mario A Storti. Fluid-structure interaction study of the start-up of a rocket engine nozzle[J]. Computers & Fluids, 2010, 39(7):1208-1218.

[13] Guerraiche Bilal. Invstigation on mechanical behavior of composite propellant during ignition process[D].西安: 西北工業(yè)大學(xué), 2013.

(編輯：崔賢彬)

A solid-fluid interaction model for propellant deformation under impact condition

JIN Ming-jun，LI Qiang

(Science and Technology on Combustion, Internal Flow and Thermal-Structure Laboratory,Northwestern Polytechnical University, Xi'an 710072, China)

The ignition of SRM is a complicated and changeful physicochemical process, which has the features of extremely short duration and particularly high temperature & pressure gradient. An applicable solid-fluid interaction simulation model was established to analyze the phenomenon of propellant deformation under impact effect for the purpose of dealing with grain structure integrity problem during SRM ignition transient. RANS and ALE method were used independently in fluid and solid sub-domain. The interaction between two sub-domains was accomplished by the transmission of interface forces. The model was validated with a cold flow test. The results of simulation agree well with test data. It shows that the model has a relatively high accuracy and could be used to investigate the ignition process in a practical SRM.

solid rocket motor；ignition transient；propellant deformation；solid-fluid interaction

2016-02-25；

2016-10-27。

金銘君(1988—)，男，碩士生，從事固體火箭發(fā)動(dòng)機(jī)點(diǎn)火過(guò)程流固耦合研究。E-mail:jackking413@163.com

V512

A

1006-2793(2017)02-0158-06

10.7673/j.issn.1006-2793.2017.02.005