基于鍵相插值法的葉片振動(dòng)測(cè)量研究*

王維民， 任三群， 陳立芳， 邵化金

(北京化工大學(xué)高端機(jī)械裝備健康監(jiān)控與治愈化北京市重點(diǎn)實(shí)驗(yàn)室 北京，100029)

?

基于鍵相插值法的葉片振動(dòng)測(cè)量研究*

王維民， 任三群， 陳立芳， 邵化金

(北京化工大學(xué)高端機(jī)械裝備健康監(jiān)控與治愈化北京市重點(diǎn)實(shí)驗(yàn)室 北京，100029)

研究了轉(zhuǎn)速波動(dòng)對(duì)葉尖定時(shí)測(cè)量精度的影響及解決方法，以非接觸式葉尖定時(shí)測(cè)振法為基礎(chǔ)，引入鍵相插值法解決如升(降)速、轉(zhuǎn)速波動(dòng)等情況下的單鍵相渦輪機(jī)葉片振動(dòng)測(cè)量誤差較大的問題，同時(shí)對(duì)恒速運(yùn)轉(zhuǎn)下也可自適應(yīng)測(cè)量。首先通過建立仿真模型比較在不同波動(dòng)頻率下鍵相插值法與單鍵相法的結(jié)果，結(jié)果表明在波動(dòng)周期較大或波動(dòng)的轉(zhuǎn)速幅值較大等情況下，使用鍵相插值法計(jì)算的結(jié)果，其精確度明顯提高。最后分別通過模擬葉輪實(shí)驗(yàn)臺(tái)和真實(shí)葉輪實(shí)驗(yàn)臺(tái)進(jìn)行了實(shí)驗(yàn)驗(yàn)證。結(jié)果表明，筆者所提出的鍵相插值法可以有效地提高轉(zhuǎn)速波動(dòng)過程中的葉片振動(dòng)位移測(cè)量，對(duì)解決工程實(shí)際問題具有較大的參考價(jià)值。

葉尖定時(shí)； 葉片振動(dòng)； 轉(zhuǎn)速波動(dòng)； 單鍵相； 鍵相插值

引 言

在葉片振動(dòng)測(cè)量領(lǐng)域，葉尖定時(shí)法(blade tip timing,簡(jiǎn)稱BTT)是一種非常重要的非接觸式測(cè)量手段，相較于應(yīng)變片法，該方法由于可對(duì)整機(jī)葉片進(jìn)行測(cè)量，不會(huì)帶來氣流擾動(dòng)影響等優(yōu)點(diǎn)而受到眾多透平企業(yè)和研究人員的重視[1-2]。近年來，葉尖定時(shí)法主要用在葉片振動(dòng)測(cè)試[3]、非接觸式應(yīng)力測(cè)量[4]，裂紋故障的預(yù)警[5]、失諧葉盤的振型分析[6]以及葉盤系統(tǒng)的阻尼評(píng)價(jià)[7]等領(lǐng)域。

以上應(yīng)用對(duì)葉片振動(dòng)測(cè)量精度有較高的要求，如何精確測(cè)量葉片的振動(dòng)幅值成為關(guān)鍵。葉尖定時(shí)測(cè)量法基本原理為通過計(jì)算旋轉(zhuǎn)葉片經(jīng)過傳感器的理論到達(dá)時(shí)間和實(shí)際到達(dá)時(shí)間的差值，將該時(shí)間差值與對(duì)應(yīng)葉片的線速度相乘并經(jīng)過相關(guān)算法分析進(jìn)而獲得葉片振動(dòng)幅值、相位、頻率等信息[8-9]。設(shè)置鍵相的目的是計(jì)算旋轉(zhuǎn)轉(zhuǎn)速，及以鍵相為參考計(jì)算葉尖理論到達(dá)時(shí)間[10]。然而在實(shí)際的測(cè)量中，葉尖的理論到達(dá)時(shí)間受轉(zhuǎn)速波動(dòng)的影響，而實(shí)際達(dá)到時(shí)間又受到葉片加工誤差、軸向位移、徑向振動(dòng)等因素的影響。文獻(xiàn)[9]提出了轉(zhuǎn)速波動(dòng)的影響，但并未就具體的解決方案進(jìn)行深入的研究。文獻(xiàn)[1]研究了葉尖幾何形狀制造誤差多傳感器不對(duì)中的影響并提供了專利技術(shù)，文獻(xiàn)[12]提出了不依賴鍵相傳感器的葉片振動(dòng)測(cè)量方法，但該方法對(duì)葉片加工的誤差信息有較高的要求。

針對(duì)在計(jì)時(shí)過程中以單鍵相作為參考會(huì)導(dǎo)致計(jì)算實(shí)際到達(dá)的時(shí)間誤差增大的情況，筆者引入鍵相插值方法來計(jì)算實(shí)際到達(dá)時(shí)間，通過理論推導(dǎo)，建模仿真驗(yàn)證該方法的可行性，并通過兩個(gè)實(shí)驗(yàn)臺(tái)，開展實(shí)驗(yàn)研究。研究結(jié)果表明，鍵相插值法可以有效地提高轉(zhuǎn)速波動(dòng)過程中的葉片振動(dòng)位移測(cè)量，可為后續(xù)的葉片振動(dòng)參數(shù)分析提供較好的數(shù)據(jù)基礎(chǔ)。

1 葉片振動(dòng)位移分析方法

1.1 葉尖計(jì)時(shí)法簡(jiǎn)介(單鍵相法)

使用單鍵相時(shí)的葉尖計(jì)時(shí)測(cè)量原理示意圖如圖1所示。圖中：B1，B2，…分別表示葉片數(shù)；tkn，tkn+1為健相時(shí)間戳。

(1)

(2)

故

(3)

其中：n為轉(zhuǎn)數(shù)；b為葉片編號(hào)，1,2,3，…，B；B為葉片數(shù)；tKn,tKn+1為鍵相時(shí)間戳；vn為n轉(zhuǎn)時(shí)的線速度；R為葉片半徑；xDC為恒偏量或稱直流分量，可通過低轉(zhuǎn)速標(biāo)定或理論計(jì)算獲得；xbn為葉片b在n轉(zhuǎn)時(shí)的振動(dòng)位移。

當(dāng)運(yùn)行轉(zhuǎn)數(shù)恒定時(shí)，即Tn保持不變，上式對(duì)各葉片振動(dòng)位移的計(jì)算是準(zhǔn)確的。但是，當(dāng)葉片在旋轉(zhuǎn)的過程中，由于各種原因?qū)е罗D(zhuǎn)速不穩(wěn)定則會(huì)導(dǎo)致在計(jì)算葉片振動(dòng)位移時(shí)引入較大的計(jì)算誤差。故此引入鍵相插值法來降低由于單鍵相對(duì)振動(dòng)位移計(jì)算的誤差影響。

1.2 鍵相插值法理論推導(dǎo)

本方法主要針對(duì)單鍵相測(cè)振系統(tǒng)對(duì)葉片振動(dòng)位移計(jì)算方法作了改進(jìn)。在tKn,tKn+1之間通過線性插值的方法虛擬出B-1個(gè)參考鍵相點(diǎn)，由于這些參考鍵相點(diǎn)在單鍵相測(cè)振系統(tǒng)中是不存在的，稱為鍵相插值點(diǎn)。與第1個(gè)參考鍵相點(diǎn)tKn共同組成本轉(zhuǎn)內(nèi)B個(gè)葉片相對(duì)應(yīng)的參考鍵相。鍵相插值的原理示意圖見圖1，圖中分析對(duì)比在轉(zhuǎn)速波動(dòng)情況下的振動(dòng)位移的表示方法。

鍵相插值點(diǎn)可由下式確定

(4)

其中：B為葉片總數(shù)；tKbn為參考鍵相，即葉片b在n轉(zhuǎn)時(shí)的參考鍵相。

根據(jù)式(3)可得，鍵相插值法的葉片振動(dòng)位移為

(5)

化簡(jiǎn)得

(6)

故可得鍵相插值情況下，葉片振動(dòng)位移計(jì)算公式為

(7)

2 仿真模型建立與結(jié)果驗(yàn)證

2.1 仿真模型介紹

為了通過仿真方法驗(yàn)證鍵相插值法的可行性，需要得到與實(shí)際測(cè)量過程相類似的一系列時(shí)間點(diǎn)來做分析。一般葉片在旋轉(zhuǎn)過程中以近于諧振動(dòng)形式振動(dòng)的[13]，振動(dòng)位移方程可簡(jiǎn)化為

(8)

其中：A為振動(dòng)幅值；ω為振動(dòng)角頻率；φ0為初始相位。

仿真中為了簡(jiǎn)化復(fù)雜度，只考慮一階振動(dòng)，忽略高階振動(dòng)。將高速旋轉(zhuǎn)葉片的振動(dòng)過程簡(jiǎn)化為葉片旋轉(zhuǎn)和葉片振動(dòng)，單位統(tǒng)一用弧度表示，這也簡(jiǎn)化了葉片振動(dòng)方程的表述。其中主弧度方程表示葉片的旋轉(zhuǎn)過程(旋轉(zhuǎn)方程)

(9)

(10)

其中：N為仿真點(diǎn)數(shù)，與時(shí)間有關(guān)的量；Ω為在點(diǎn)i時(shí)的轉(zhuǎn)速，單位r/min；t為時(shí)間，本仿真是按1s為108步長(zhǎng)計(jì)算。

振動(dòng)弧度方程表示葉片在旋轉(zhuǎn)過程中的振動(dòng)，其中，f表示葉片的振動(dòng)頻率。由式(8)得

(11)

由式(9)和式(11)得出，葉片在N時(shí)(或t=N/108)的振動(dòng)累積弧度為

(12)

為了仿真葉片經(jīng)過傳感器時(shí)被傳感器捕捉到實(shí)際到達(dá)傳感器的時(shí)間。仿真中每一轉(zhuǎn)內(nèi)，使葉片的弧度累積到一定值時(shí)的第一個(gè)點(diǎn)即為實(shí)際到達(dá)傳感器的時(shí)間。為了方便計(jì)算，根據(jù)旋轉(zhuǎn)可得，檢測(cè)方程為

(13)

鍵相傳感器的檢測(cè)方式與上式相同，可減少編程計(jì)算的復(fù)雜度。

為了與實(shí)驗(yàn)相參照，假設(shè)葉片個(gè)數(shù)B=16，取鍵相與葉片1之間的夾角為0，各葉片等間距分布，葉片振動(dòng)初始相位為0。仿真的葉片編號(hào)為1, 4, 6, 8, 10, 12, 14, 15, 16，其對(duì)應(yīng)的葉片振動(dòng)頻率為1 806 Hz；葉片編號(hào)2的振動(dòng)頻率為1 807 Hz[14]。取A/R=0.004。則葉片1, 4, 6, 8, 10, 12, 14, 15, 16的振動(dòng)累積弧度方程可表示為

(14)

葉片2的振動(dòng)累積弧度方程為

(15)

轉(zhuǎn)速波動(dòng)方程為

(16)

其中：An為波動(dòng)轉(zhuǎn)速的波動(dòng)幅值；fn為波動(dòng)轉(zhuǎn)速的波動(dòng)頻率。

2.2 仿真及結(jié)果分析

在上述仿真方程的基礎(chǔ)上，基于4 800r/min做了以下6種情況的仿真實(shí)驗(yàn)(為了簡(jiǎn)化將初始相位角設(shè)為0)?；鶞?zhǔn)誤差是仿真實(shí)驗(yàn)中葉片經(jīng)過傳感器時(shí)的振動(dòng)位移獲得的振動(dòng)幅值與理論振動(dòng)幅值之比。

表1 轉(zhuǎn)速波動(dòng)的幅值和頻率

通過以上仿真結(jié)果可以得出，鍵相插值法較單鍵相法在轉(zhuǎn)速波動(dòng)過大，且波動(dòng)頻率不大時(shí)改善顯著。如圖2，3所示，當(dāng)轉(zhuǎn)速波動(dòng)幅值為5 r/min，頻率為11 Hz時(shí)基準(zhǔn)誤差與鍵相插值法誤差相近，而單鍵相法偏離較大，葉片編號(hào)10以后的數(shù)據(jù)誤差較大。同樣，當(dāng)波動(dòng)幅值增大到10和20 r/min時(shí)，鍵相插值法仍然具有較高的精度，而單鍵相法誤差逐步增大。當(dāng)波動(dòng)頻率為55 Hz時(shí)，在低波動(dòng)幅值下(5 r/min),鍵相插值法仍然具有較高的精度，但是當(dāng)波動(dòng)幅值增大到10 r/min甚至更高時(shí)(20 r/min)，鍵相插值法的改善效果會(huì)有所降低，但單鍵相法誤差更大，此時(shí)需要通過其他方法來解決由于轉(zhuǎn)速波動(dòng)過大而導(dǎo)致的葉片振動(dòng)測(cè)量誤差較大的問題。由此可見，在一定程度上引入鍵相插值會(huì)提高葉片振動(dòng)的測(cè)量精度，且不會(huì)帶來附加影響。

圖2 轉(zhuǎn)速波動(dòng)頻率為11 Hz時(shí)各波動(dòng)幅值下的誤差對(duì)比Fig.2 Comparing vibration amplitude error under different fluctuation amplitude at the speed fluctuation of 11 Hz

圖3 轉(zhuǎn)速波動(dòng)頻率為55 Hz時(shí)各波動(dòng)幅值下的誤差對(duì)比Fig.3 Comparing vibration amplitude error under different fluctuation amplitude at the speed fluctuation of 55 Hz

2.3 實(shí)驗(yàn)臺(tái)1對(duì)鍵相插值法的驗(yàn)證

實(shí)驗(yàn)臺(tái)由模擬渦輪葉片的柔性直葉片(葉片直徑100 mm，葉片個(gè)數(shù)為16)、cDAQ-9174采集設(shè)備(帶NI9402采集卡)、一套Bently 3300 XL 8 mm電渦流傳感器、兩套光纖傳感器、光電探測(cè)器、激光源、變頻器、無刷直流電機(jī)、驅(qū)動(dòng)器及其他輔件組成，如圖4所示。

圖4 直葉片振動(dòng)監(jiān)測(cè)實(shí)驗(yàn)臺(tái)及葉片F(xiàn)ig.4 Blade vibration monitoring apparatus

在上述實(shí)驗(yàn)臺(tái)的基礎(chǔ)上，進(jìn)行了4 200 r/min和2 400 r/min兩種轉(zhuǎn)速的測(cè)試。4 200 r/min時(shí)，取葉片編號(hào)1的振動(dòng)位移放大圖，如圖5所示。

圖5 單鍵相法測(cè)得葉片編號(hào)1的振動(dòng)位移Fig.5 Blade No.1 vibration displacement by KPIM

圖6 4 200 r/min時(shí)單鍵相法與鍵相插值法結(jié)果對(duì)比Fig.6 Comparing the results of single phase method and KPIM at the speed of 4 200 r/min

圖6為單鍵相法和鍵相插值法求得的葉片振動(dòng)幅值之間的對(duì)比。通過分析轉(zhuǎn)速情況可知，轉(zhuǎn)速波動(dòng)幅值小于6 r/min，波動(dòng)頻率超過30 Hz(且有噪聲干擾)，兩種方法測(cè)出的葉片振動(dòng)幅值接近，與上述仿真結(jié)論較為接近。運(yùn)行轉(zhuǎn)速為2 400 r/min時(shí)， 兩種方法獲得的葉片振動(dòng)幅值大小差別幾乎為0，如圖7所示。分析轉(zhuǎn)速發(fā)現(xiàn)，轉(zhuǎn)速波動(dòng)幅值小于2 r/min，且波動(dòng)頻率在20～30 Hz(有噪聲干擾)。以上兩種轉(zhuǎn)速實(shí)驗(yàn)與仿真結(jié)果得出的結(jié)論一致。

圖7 2 400 r/min下單鍵相法與鍵相插值法結(jié)果對(duì)比Fig.7 Comparing the results of single phase method and KPIM at the speed of 2 400 r/min

2.4 實(shí)驗(yàn)臺(tái)2介紹與實(shí)驗(yàn)結(jié)果

實(shí)驗(yàn)臺(tái)2中采用扭轉(zhuǎn)葉片，如圖8所示。葉盤外圓直徑137 mm，葉片個(gè)數(shù)為32，整個(gè)渦輪葉片采用鋁合金數(shù)控銑加工而成。

圖8 實(shí)驗(yàn)臺(tái)與扭葉片模型Fig.8 The test bench and twisted blade model

在上述實(shí)驗(yàn)臺(tái)的基礎(chǔ)上為了探究鍵相插值法的可行性，進(jìn)行了4 800 r/min的運(yùn)轉(zhuǎn)實(shí)驗(yàn)。運(yùn)用單鍵相法測(cè)得葉片編號(hào)1的振動(dòng)位移，如圖9所示。

通過圖9可以看到兩個(gè)現(xiàn)象:一是葉片偏向一側(cè)振動(dòng)位移較大，而另一側(cè)很?。欢请S著運(yùn)轉(zhuǎn)時(shí)間的增加，葉片的平衡位置會(huì)向上漂移。在對(duì)6 000,7 200,8 400,9 600 r/min這4個(gè)轉(zhuǎn)速的實(shí)驗(yàn)中，也都出現(xiàn)上述兩個(gè)現(xiàn)象。

現(xiàn)象一的發(fā)生是由于扭葉片本身結(jié)構(gòu)原因所導(dǎo)致的。葉片彎曲振動(dòng)可分為切向振動(dòng)和軸向振動(dòng)，分別沿最小主慣性軸和最大主慣性軸振動(dòng)，扭葉片的最小主慣性軸與軸線不平行，且左右的抗彎截面系數(shù)不同會(huì)導(dǎo)致葉片出現(xiàn)上述振動(dòng)現(xiàn)象。對(duì)現(xiàn)象二進(jìn)行了長(zhǎng)時(shí)間的運(yùn)轉(zhuǎn)實(shí)驗(yàn)與監(jiān)測(cè)，發(fā)現(xiàn)在運(yùn)轉(zhuǎn)超過4 min后該漂移現(xiàn)象趨于平緩，10 min之后平衡位置幾乎不變。這種現(xiàn)象發(fā)生原因是由于在初始運(yùn)轉(zhuǎn)的過程中葉片對(duì)熱膨脹效應(yīng)較為敏感，導(dǎo)致葉頂發(fā)生軸向位移，進(jìn)而表現(xiàn)為葉頂?shù)竭_(dá)傳感器的時(shí)間發(fā)生改變，隨著運(yùn)轉(zhuǎn)時(shí)間的增長(zhǎng)，熱膨脹效應(yīng)達(dá)到極限，進(jìn)而平衡位置不再發(fā)生改變。

圖9 葉片編號(hào)1振動(dòng)位移圖Fig.9 Blade No.1 vibration displacement diagram

相較于直葉片，扭葉片振動(dòng)形式更為復(fù)雜，導(dǎo)致鍵相插值法的使用受到較大限制，但扭葉片應(yīng)用更為普遍，需要繼續(xù)研究真實(shí)葉片振動(dòng)及測(cè)試方法。

3 結(jié) 論

1) 完成了鍵相插值法的理論推導(dǎo)，并建立了仿真數(shù)學(xué)模型。進(jìn)行了不同波動(dòng)轉(zhuǎn)速下的仿真，仿真的結(jié)果表明在對(duì)轉(zhuǎn)速波動(dòng)較大且波動(dòng)頻率不大時(shí)會(huì)對(duì)測(cè)量結(jié)果帶來很大的改善，同時(shí)對(duì)波動(dòng)頻率較高或恒轉(zhuǎn)速情況下可以自適應(yīng)。

2) 通過直葉片實(shí)驗(yàn)臺(tái)進(jìn)行了葉片振動(dòng)位移測(cè)量實(shí)驗(yàn)，在運(yùn)行轉(zhuǎn)速為2 400 r/min和4 200 r/min時(shí)，兩種方法測(cè)得的振動(dòng)幅值的大小較為接近，通過對(duì)轉(zhuǎn)速進(jìn)行分析可知由于轉(zhuǎn)速波動(dòng)很小導(dǎo)致，這也驗(yàn)證了仿真的結(jié)果。

3) 為了與工程應(yīng)用中使用的真實(shí)葉片更為接近，使用了扭葉片進(jìn)行相關(guān)實(shí)驗(yàn)，在監(jiān)測(cè)葉片振動(dòng)位移的過程中，出現(xiàn)了葉片單側(cè)振動(dòng)位移較大，另一側(cè)振動(dòng)位移很小的現(xiàn)象，同時(shí)隨著運(yùn)轉(zhuǎn)時(shí)間的增加平衡位置有漂動(dòng)。上述現(xiàn)象的出現(xiàn)，雖然不能為驗(yàn)證鍵相插值法提供依據(jù)，但為進(jìn)一步解決工程實(shí)際問題提供了更多的思考方向。

[1] Kaszynski A A, Brown J M. Accurate blade tip timing limits through geometry mistuning modeling[C]∥ASME Turbo Expo 2015: Turbine Technical Conference and Exposition. Montréal, Canada:American Society of Mechanical Engineers, 2015.

[2] Knappett D, Garcia J. Blade tip timing and strain gauge correlation on compressor blades[J]. Proceedings of the Institution of Mechanical Engineers, Part G: Journal of Aerospace Engineering, 2008, 222(4):497-506.

[3] Gallego-Garrido J, Dimitriadis G, Wright J R. A class of methods for the analysis of blade tip timing data from bladed assemblies undergoing simultaneous resonances—part i: theoretical development[J]. International Journal of Rotating Machinery, 2007, 2007:11.

[4] Pickering T M. Methods for validation of a turbomachinery rotor blade tip timing system[D]. Blacksburg: Virginia Polytechnic Institute and State University, 2014.

[5] Hu Zheng, Lin Jun, Chen Zhongsheng, et al. A non-uniformly under-sampled blade tip-timing signal reconstruction method for blade vibration monitoring[J]. Sensors, 2015, 15(2):2419-2437.

[6] Salhi B, Lardiès J, Berthillier M, et al. Modal parameter identification of mistuned bladed disks using tip timing data[J]. Journal of Sound & Vibration, 2008, 314(3):885-906.

[7] Abdelrhman A M, Hee L M, Leong M S, et al. Condition monitoring of blade in turbomachinery: a review[J]. Advances in Mechanical Engineering, 2014, 2014:10.

[8] 王萍. 葉尖定時(shí)方法在國外航空發(fā)動(dòng)機(jī)葉片振動(dòng)測(cè)量中的應(yīng)用綜述[J]. 航空科學(xué)技術(shù), 2013,6(6):5-9.

Wang Ping. Application review of blade tip timing method for foreign aero-engine blade vibration measurement[J]. Aeronautical Science & Technology, 2013,6(6):5-9. (in Chinese)

[9] Carrington I B. Development of blade tip timing data analysis techniques[D]. Manchester: University of Manchester, 2002.

[10]鄧勇, 盛鑫志, 薛秀生. 光纖葉片振幅監(jiān)測(cè)系統(tǒng)應(yīng)用試驗(yàn)研究[J]. 航空動(dòng)力學(xué)報(bào), 1992,2(2):156-158.

Deng Yong, Sheng Xinzhi, Xue Xiusheng. A fibre-optic amplitude monitor for blade vibration[J]. Journal of Aerospace Power,1992,2(2):156-158. (in Chinese)

[11]Twerdochlib M. Method of matching sensors in a multi-probe turbine blade vibration monitor:US,US7987725[P], 2011-08-02.

[12]Guo Haotian, Duan Fajie, Zhang Jilong. Blade resonance parameter identification based on tip-timing method without the once-per revolution sensor[J]. Mechanical Systems and Signal Processing, 2016, 66(2016): 625-639.

[13]薩勃洛斯基ИE. 渦輪機(jī)葉片振動(dòng)的非接觸測(cè)量[M].吳士祥,鄭叔琛,譯.北京:國防工業(yè)出版社,1986:40-43.[14]張玉貴. 煙氣輪機(jī)葉片振動(dòng)的非接觸式在線監(jiān)測(cè)關(guān)鍵技術(shù)研究[D]. 天津:天津大學(xué), 2008.

*國家自然科學(xué)基金面上資助項(xiàng)目(51275028);國家自然科學(xué)基金重點(diǎn)資助項(xiàng)目(51135001);中央高校基本科研業(yè)務(wù)費(fèi)專項(xiàng)資金資助項(xiàng)目(JD1508)

2015-04-17；

2015-09-20

10.16450/j.cnki.issn.1004-6801.2017.02.025

V216.2; TH45

王維民，男，1978年4月生，博士、教授、博士生導(dǎo)師。主要研究方向?yàn)樾D(zhuǎn)機(jī)械動(dòng)力學(xué)分析、狀態(tài)監(jiān)測(cè)及故障診斷。曾發(fā)表《An identification method for damping ratio in rotor system》(《Mechanical Systems and Signal Processing》2016,Vol.68-69)等論文。 E-mail：wwmbuct@163.com