薛曉川,王志強

(北京航空航天大學儀器科學與光電工程學院，北京 100191)

磁懸浮儲能飛輪振動狀態(tài)監(jiān)測研究

薛曉川,王志強

(北京航空航天大學儀器科學與光電工程學院，北京 100191)

針對磁懸浮儲能飛輪系統(tǒng)的不平衡振動會影響轉(zhuǎn)子穩(wěn)定運轉(zhuǎn)和系統(tǒng)正常工作的問題，對其擾動原因進行研究。通過建立含有轉(zhuǎn)子不平衡量的廣義干擾方程和磁懸浮飛輪動力學方程，對不平衡振動中的同頻、倍頻及模態(tài)分量進行理論分析。設計并實現(xiàn)了一套基于LabVIEW-DAQ的磁懸浮飛輪振動狀態(tài)監(jiān)測系統(tǒng)，完成對儲能飛輪振動信號的實時采集、處理和顯示，實現(xiàn)了飛輪的振動狀態(tài)監(jiān)測,并對一臺300 kW儲能飛輪進行試驗驗證。針對儲能飛輪工作的強電環(huán)境引起的信號噪聲，對原始信號進行小波閾值消噪處理。對消噪后的振動信號進行頻域分析，獲取某一轉(zhuǎn)速下飛輪的振動頻譜圖和不同轉(zhuǎn)速下飛輪振動信號的三維瀑布圖。試驗結(jié)果表明，該飛輪振動頻譜中存在穩(wěn)定的飛輪旋轉(zhuǎn)頻率成分及其倍頻分量。該系統(tǒng)操作簡單，通用性強，具有一定的市場前景和應用價值。

儲能飛輪； 磁懸?。?不平衡振動； 頻率成分； 數(shù)據(jù)采集； 小波消噪； 頻域分析； 狀態(tài)監(jiān)測

0 引言

飛輪儲能技術作為一種新型儲能方式，以其高效、環(huán)保等優(yōu)勢，具有廣泛的發(fā)展前景，并已在電網(wǎng)系統(tǒng)、軌道交通等領域發(fā)揮了巨大作用[1]。其原理是通過飛輪升降速，實現(xiàn)電能的存儲和釋放。相比機械軸承，磁軸承具有無摩擦、體積小、長壽命等特點[2]，將磁軸承應用于儲能飛輪系統(tǒng)是近幾年的研究熱點。

磁懸浮儲能飛輪工作時常伴有不平衡振動，會影響其工作性能，嚴重時甚至會引起飛輪失穩(wěn)。對振動源進行識別和監(jiān)測，已成為飛輪研究的重要問題。文獻[3]、文獻[4]中的飛輪振動測試系統(tǒng)，能夠?qū)崿F(xiàn)小型反作用力飛輪的振動監(jiān)測，但不適用于大型儲能飛輪。文獻[5]采用廣義Lagrange-Maxwell方程，從機電耦合特性的角度，對磁懸浮飛輪系統(tǒng)振動進行了理論分析和仿真驗證，但缺乏試驗驗證。

本文針對磁懸浮飛輪的不平衡振動問題，通過建立廣義干擾和轉(zhuǎn)子動力學模型，對振動源中的同頻、倍頻和模態(tài)分量進行理論分析，設計并實現(xiàn)了一套基于LabVIEW-DAQ的飛輪振動監(jiān)測系統(tǒng)，以識別振動源中的各個分量并進行分析，為飛輪的故障診斷和狀態(tài)預測提供依據(jù)。

1 磁懸浮飛輪振動特性分析

儲能飛輪作為復雜的機電耦合系統(tǒng)，其振動原因較為復雜，但具有基本相同的表現(xiàn)形式。根據(jù)其頻域特征，可主要劃分為同頻分量、倍頻分量和模態(tài)分量?，F(xiàn)通過建立數(shù)學模型的方法，對引起磁懸浮飛輪轉(zhuǎn)子不平衡的各種分量進行分析。

1.1 飛輪振動中的同頻分量

磁懸浮飛輪正常工作時，振動成分中存在同頻分量，其頻率等于飛輪的轉(zhuǎn)動頻率，且振動幅值較大。通過建立廣義干擾和轉(zhuǎn)子動力學模型的方法，對其影響因素進行分析。轉(zhuǎn)子廣義干擾示意圖如圖1所示。以飛輪質(zhì)心O為坐標原點建立空間直角坐標系，Z軸為轉(zhuǎn)子轉(zhuǎn)軸，X平面、Y平面過原點且與Z軸互相垂直。X軸、Y軸方向隨轉(zhuǎn)子旋轉(zhuǎn)而不斷變化。

圖1 轉(zhuǎn)子廣義干擾示意圖

轉(zhuǎn)子不平衡量是引起飛輪不平衡振動的重要原因，這源于加工制造工藝的限制,可以通過轉(zhuǎn)子動平衡校正來減小，但無法徹底消除。其表現(xiàn)形式為飛輪工作時旋轉(zhuǎn)主軸和慣性主軸不重合。磁懸浮飛輪的工作原理是將飛輪轉(zhuǎn)子懸浮于轉(zhuǎn)子的幾何軸，因此當飛輪繞轉(zhuǎn)子幾何軸旋轉(zhuǎn)時，飛輪的慣性主軸就會繞幾何軸轉(zhuǎn)動[6]，產(chǎn)生慣性力和慣性力偶矩使轉(zhuǎn)子發(fā)生不平衡振動?，F(xiàn)從廣義干擾的角度，對飛輪的不平衡振動原因進行分析。

根據(jù)不平衡的定義可知，不平衡量存在時的一種表現(xiàn)形式為力和力矩的廣義干擾作用于飛輪[7]。ω(t)為飛輪的轉(zhuǎn)動角頻率；F(t)為所有引起飛輪振動的轉(zhuǎn)子慣性離心力的合力，即主矢；M(t)為慣性離心力矩的合力偶矩，即主矩；θ為初始相位差。F(t)和M(t)的方向隨轉(zhuǎn)子的轉(zhuǎn)動而改變，其頻率均等于飛輪的旋轉(zhuǎn)頻率。分別將F、M在X軸、Y軸上進行分解，得：

(1)

(2)

將主矢和主矩的共同作用看作一種廣義干擾，在X軸、Y軸上進行分解，分量為Ax(t)和Ay(t)。由于主矢和主矩的頻率均等于飛輪轉(zhuǎn)頻，故廣義干擾分量的頻率也與飛輪旋轉(zhuǎn)頻率相等。設其幅值分別為Ax0和Ay0，φ為干擾量相對于坐標系的初始相位，則有：

(3)

在不考慮磁軸承控制的前提下，單獨分析引入轉(zhuǎn)子不平衡激勵帶來的影響。此時可將飛輪簡化為一個單自由度的線性系統(tǒng)。廣義干擾Ax(t)、Ay(t)為系統(tǒng)的激勵源，可以看出在X方向和Y方向上，會產(chǎn)生穩(wěn)定的同頻分量，表現(xiàn)為飛輪以等同于轉(zhuǎn)頻的頻率振動。

此外，轉(zhuǎn)子的不平衡量還會激發(fā)磁軸承電流剛度力和位移剛度力中的同頻分量。在磁懸浮飛輪系統(tǒng)中，轉(zhuǎn)子的不平衡量會導致磁軸承中心平面對應點的偏移。磁懸浮飛輪動力學簡化圖如圖2所示。

圖2 磁懸浮飛輪動力學簡化圖

圖2中：外圓為徑向磁軸承，內(nèi)圓為飛輪轉(zhuǎn)子。

以磁軸承的幾何中心O′為坐標原點建立直角坐標系，X軸、Y軸的方向與位移傳感器的安裝方向一致。由于轉(zhuǎn)子不平衡量的存在，使得其幾何軸和慣性主軸不重合。轉(zhuǎn)子以ω(t)的轉(zhuǎn)速作逆時針旋轉(zhuǎn),O為轉(zhuǎn)子幾何中心，P為轉(zhuǎn)子質(zhì)心，ε為OP向量的模長。以O為坐標原點，建立轉(zhuǎn)子的轉(zhuǎn)動坐標系，對應坐標軸為M和N，且方向隨轉(zhuǎn)子轉(zhuǎn)動時刻變化。OP向量與M軸的夾角為θ′。t時刻M軸與X軸夾角為ω(t)t。X(t)和Y(t)為轉(zhuǎn)子的幾何中心相對于磁軸承幾何中心的位移分量，為磁軸承的控制量。X′(t)和Y′(t)為轉(zhuǎn)子的質(zhì)心相對于磁軸承幾何中心的位移分量。

由圖2中的幾何關系，得到如下方程組：

(4)

磁懸浮飛輪正常工作時的控制電流一般小于偏置電流，此時可得到磁懸浮飛輪徑向轉(zhuǎn)子動力學方程，如式(5)所示。

(5)

式中：Fx、Fy分別為磁軸承在X和Y方向的控制力；μ為真空磁導率；n為線圈匝數(shù)；S為磁鐵的橫截面積；σ為磁鐵磁極偏角；iy為偏置電流；ix、iy分別為X和Y方

向的控制電流；X0、Y0分別為轉(zhuǎn)子在平衡位置時X方向和Y方向的位移。

當轉(zhuǎn)子在平衡位置處的位移X=X0、Y=Y0時，可將方程組作線性化處理，分別對電流、位移求偏導，在X、Y趨于0的條件下作一階Taylor展開,得：

(6)

化簡得：

(7)

實際控制時，控制電流為轉(zhuǎn)子位移的PID函數(shù)，因此,不平衡引入的擾動會同時引入到位移項和電流項之中，代入式(6)可得:

(8)

由此可知，轉(zhuǎn)子的不平衡量在磁軸承控制的過程中引入同頻干擾。

1.2 飛輪振動中的倍頻分量

在對飛輪進行磁懸浮時，需根據(jù)位移傳感器的測量值對磁軸承的軸向和徑向位移進行控制。飛輪倍頻干擾主要由傳感器諧波噪聲引起，最終表現(xiàn)為對磁軸承控制電流項的影響。其數(shù)學表達式為：

(9)

式中：λx、λy分別為位移傳感器引入的噪聲干擾；λx,i、λy,i分別為各次諧波對應的干擾幅值；Ψx,i、Ψy,i分別為各次諧波干擾項所對應的初始相位。

將λx、λy代入式(8),可得：

(10)

由此可知，傳感器倍頻擾動會在磁軸承轉(zhuǎn)子的動力學方程中引入倍頻分量。

1.3 飛輪振動中的模態(tài)分量

除同頻、倍頻擾動外，模態(tài)分量也是引起飛輪不平衡振動、制約飛輪轉(zhuǎn)速的重要原因。這種由模態(tài)分量引起的飛輪擾動被稱為結(jié)構(gòu)模態(tài)振動，它取決于飛輪系統(tǒng)的固有頻率。在飛輪轉(zhuǎn)速范圍內(nèi)，存在轉(zhuǎn)子固有頻率和飛輪旋轉(zhuǎn)頻率相等的情況，即引發(fā)共振[8]。系統(tǒng)中的共振頻率由質(zhì)量、剛度和衰減系數(shù)決定[9]。因此，為了避免因不平衡激勵引發(fā)的共振現(xiàn)象，應使飛輪轉(zhuǎn)頻盡量遠離飛輪的固有頻率。飛輪轉(zhuǎn)子的各階模態(tài)固有頻率可由有限元分析軟件計算得到。

2 飛輪振動狀態(tài)監(jiān)測系統(tǒng)設計

為了獲取某一時刻磁懸浮飛輪的振動數(shù)據(jù)，本文設計了一套飛輪振動狀態(tài)監(jiān)測系統(tǒng)，能夠?qū)崿F(xiàn)振動信號的采集、處理、顯示和記錄。后期可通過對數(shù)據(jù)進行小波閾值去噪、快速傅里葉變換(fast Fourier transform,FFT)頻譜分析、三維瀑布圖分析的方法進行驗證，以識別飛輪不同狀態(tài)下的振動對應頻率，進而判斷飛輪的運行狀態(tài)。

總體設計思路為：通過固連在飛輪殼體上的多個壓電式加速度計，獲取飛輪的振動信號。振動信號通過恒流適配器，傳送到集成有數(shù)據(jù)采集卡的工控機中。使用編寫的上位機軟件對其進行控制，實現(xiàn)數(shù)據(jù)的采集，處理和記錄。系統(tǒng)總體設計框圖如圖3所示。

圖3 系統(tǒng)總體設計框圖

2.1 系統(tǒng)硬件

振動監(jiān)測系統(tǒng)的硬件部分主要由工控機、NI PCI-6255數(shù)據(jù)采集卡、CT5200恒流適配器和CT-1050LC加速度計等模塊實現(xiàn)。硬件部分設計框圖如圖4所示。

圖4 硬件部分設計框圖

數(shù)據(jù)采集卡選用NI公司的PCI-6255數(shù)據(jù)采集卡，采用擁有較高的分辨率和敏感度的高度精M系列模塊。傳感器安裝位置如圖5所示。

圖5 傳感器安裝位置示意圖

測量敏感元件選用CT1050系列壓電式加速度傳感器。該傳感器基于壓電效應原理，具有頻域?qū)?、靈敏度高、橫向靈敏度小等優(yōu)點。測試系統(tǒng)使用六個壓電式加速度傳感器與飛輪殼體固連。兩組互相正交、水平放置，測量飛輪徑向方向振動；兩個傳感器豎直放置，測量輪體軸向方向振動。

2.2 系統(tǒng)軟件

LabVIEW虛擬儀器具有靈活性好、開發(fā)費用低、周期短、人機交互性好的特點。本試驗系統(tǒng)軟件部分采用圖形化編程軟件LabVIEW平臺中的DAQ數(shù)據(jù)采集功能[10]。其流程是：先設定數(shù)據(jù)采集的工作狀態(tài)，如通道數(shù)、采樣率、數(shù)據(jù)記錄文件路徑等；再運行振動信號數(shù)據(jù)采集程序，并將各個通道振動信號的時域和頻域波形圖實時顯示在主面板上。軟件部分設計框圖如圖6所示。

圖6 軟件部分設計框圖

3 試驗驗證與數(shù)據(jù)分析

利用試驗室研制的一臺300 kW磁懸浮儲能飛輪進行試驗驗證。其主要技術指標為：額定功率為300 kW，額定轉(zhuǎn)速為30 000 r/min，整體質(zhì)量為855 kg，轉(zhuǎn)子直徑為314 mm，位移傳感器分辨率為2.6×10-3V/μm。

3.1 飛輪模態(tài)分析

首先分析振動模態(tài)分量的影響。在實際飛輪系統(tǒng)中，當飛輪工作頻率與其固有頻率相等時，會產(chǎn)生共振現(xiàn)象，使系統(tǒng)無法發(fā)揮最大效能。使用ANSYS軟件對飛輪轉(zhuǎn)子進行有限元分析，可得到X方向一階、二階固有頻率為697.6 Hz和1 844.4 Hz，對應振型分別為X方向一階、二階彎曲。由于儲能飛輪的額定轉(zhuǎn)速為30 000 r/min(轉(zhuǎn)頻為500 Hz)，低于轉(zhuǎn)子的一階彎曲模態(tài)頻率，因此可以將其視為剛性轉(zhuǎn)子，忽略擾動中模態(tài)分量的影響。

3.2 小波閾值去噪

對監(jiān)測系統(tǒng)獲取的飛輪振動信號進行數(shù)學處理。由于實測信號中不可避免地存在噪聲，為了保證試驗結(jié)果的準確性，可運用小波變換的自適應閾值方法對振動原始信號進行消噪處理，提取信號中的有用成分。具體實現(xiàn)手段是先對原始信號進行小波變換，得到小波變換系數(shù)計算閾值；舍去小于閾值的變換系數(shù)，對于大于閾值的變換系數(shù)，采用軟閾值方法進行處理；作小波逆變換實現(xiàn)信號的還原，實現(xiàn)圖像去噪的目的[11]。具體實現(xiàn)步驟如下。

①選取DB4和Sym4小波基對原始振動信號進行2層分解。

②計算相應子空間的均值、中值和均方根，確定最優(yōu)閾值。

③對大于閾值的系數(shù)作軟閾值處理。

④小波包重構(gòu)信號恢復。

振動信號時域波形圖如圖7所示。對比圖7(a)和圖7(b)可以發(fā)現(xiàn)，小波去噪效果明顯，驗證了該方法的有效性和正確性。

圖7 振動信號時域波形圖

3.3 FFT頻譜分析

頻域特性是數(shù)字信號的基本特性。從頻域角度對信號進行分析與處理，容易對信號的特性進行深入了解[12]。振動信號頻譜分析法是在通用信號處理理論基礎上發(fā)展起來的，它是診斷工作的依據(jù)[13]。傅里葉變換是一種積分變換，它建立了信號時域與頻域聯(lián)系的紐帶。

借助MATLAB軟件對采集到的小波去噪后的信號進行FFT處理，得到振動信號的頻域圖。通過同一轉(zhuǎn)速下某傳感器測量飛輪振動信號頻域圖，進行頻域分析，可以得出各個頻率下對應的振動幅值大小，識別對應的特征頻率。將對應的飛輪轉(zhuǎn)速、頻率、振幅繪制成三維瀑布圖，分析飛輪的振動幅值大小及其對應頻率隨轉(zhuǎn)速的變化規(guī)律。在對徑向方向傳感器測量值進行頻譜分析時，盡管X方向和Y方向的受力以及磁軸承的剛度及阻尼略有不同，但振動的本質(zhì)規(guī)律相同。因此取X方向某一傳感器的測量結(jié)果作FFT并繪制頻譜圖。飛輪轉(zhuǎn)速為4 500 r/min時的頻域波形圖如圖8所示。

不考慮信號頻譜中直流分量的影響，對磁懸浮飛輪的振動頻率進行分析。

圖8 振動頻譜圖(4 500 r/min)

由圖8可以明顯地看出，系統(tǒng)的主振頻率為75.1 Hz，基本等于飛輪在4 500 r/min狀態(tài)時的同頻振動頻率，誤差為0.13%，其振動幅值達到了0.023 07g(g為重力加速度)。其振動源是由轉(zhuǎn)子不平衡量引起的慣性離心力和力矩，以及磁軸承電流剛度力和位移剛度力中的同頻分量。這些因素的疊加會引起飛輪系統(tǒng)的同頻振動。

同時，從圖8還可以看出，系統(tǒng)中存在倍頻振動分量，即轉(zhuǎn)頻的高次諧波。其振動幅值小于同頻分量，以二倍頻、三倍頻較為明顯。該倍頻振動分量與磁軸承位移傳感器引入的諧波噪聲分量有關。信號在50 Hz左右處存在穩(wěn)定的分量，且幅值較大。該分量由工頻電信號干擾造成，實際分析時可略去。如需避免工頻干擾，可以采取使用抗電磁干擾的采集卡等措施，改善系統(tǒng)性能。

儲能飛輪徑向振動信號在某轉(zhuǎn)速范圍內(nèi)的三維轉(zhuǎn)子降速瀑布圖如圖9所示。

圖9 三維轉(zhuǎn)子降速瀑布圖

在磁懸浮飛輪自由降速的過程中，忽略直流分量和工頻干擾的影響，系統(tǒng)的振動主頻是飛輪的工作轉(zhuǎn)速頻率，且振動幅值最大。采用提高轉(zhuǎn)子動平衡試驗精度等方法減小此同頻振動量。除此之外，還可以明顯地看出，系統(tǒng)中存在轉(zhuǎn)頻的倍頻成分，以2次、3次諧波為主。通過降低位移傳感器諧波噪聲，可減小倍頻分量對系統(tǒng)的影響。同時，應保證系統(tǒng)的額定工作轉(zhuǎn)速遠離計算得到的轉(zhuǎn)子一階固有頻率，防止共振現(xiàn)象的發(fā)生。試驗現(xiàn)象與之前的理論分析相吻合，同時也表明飛輪振動狀態(tài)監(jiān)測系統(tǒng)能夠較好地完成實時采集、處理和監(jiān)測的預期功能。

4 結(jié)束語

本文以磁懸浮飛輪的不平衡振動為研究對象，分析了飛輪振動成分中的同頻、倍頻和模態(tài)分量。在廣義干擾方程和轉(zhuǎn)子動力學方程的基礎上，闡述了各種振動分量的原因和影響因素。設計并實現(xiàn)了一種飛輪振動狀態(tài)監(jiān)測系統(tǒng)，對飛輪的振動數(shù)據(jù)進行實時的采集、處理和顯示。后期對振動數(shù)據(jù)進行小波閾值消噪和頻譜分析。試驗結(jié)果表明，磁懸浮飛輪正常工作時，存在穩(wěn)定的同頻和倍頻振動分量。驗證了理論分析的合理性和本系統(tǒng)的可行性，為飛輪故障診斷和狀態(tài)預測提供依據(jù)。

[1] 陳廣宇,張華東,于寶鑫,等.飛輪儲能關鍵技術及其應用[J].電氣應用,2015,34(24):53-58.

[2] 潤章,朱熀秋,朱利東.磁懸浮飛輪軸向混合磁軸承參數(shù)設計及分析[J].機械設計與研究,2015,31(1):39-43.

[3] TANIWAKI S,OHKAMI Y.Experimental and numerical analysis of reaction wheel disturbances[J].JSME International Journal,2003,46(2):519-526.

[4] OH H S,CHEON D I.Precision measurements of reaction wheel disturbances with frequency compensation process[J].Journal of Mechanical Science and Technology,2005,19(1):136-143.

[5] 陳峻峰,劉昆,梁文杰,等.磁懸浮飛輪儲能系統(tǒng)機電耦合非線性動力學研究[J].動力學與控制學報,2013,11(3):225-234.

[6] 劉彬,房建成,劉剛,等.磁懸浮飛輪不平衡振動控制方法與試驗研究[J].機械工程學報,2010,46(12):188-194.

[7] 王全武,虎剛.飛輪振動頻譜特征的初步理論分析和驗證[J].空間控制技術與應用,2008,34(4):42-46.

[8] 臧孝華,王全武.振動測試技術在飛輪動平衡調(diào)整中的應用[J].新技術新工藝,2013(1):51-53.

[9] 王光軍,胡萬明.高速回轉(zhuǎn)機械振動產(chǎn)生的原因、測試標準與防治[J].氯堿工業(yè),2006(5):44-44.

[10]王平,楊濤,侯守全,等.LabVIEW中DAQ數(shù)據(jù)采集系統(tǒng)設計[J].自動化儀表,2015,36(7):31-33.

[11]朱希安.小波分析及其在數(shù)字圖像處理中的應用[M].北京：電子工業(yè)出版社,2012:88-95.

[12]劉小群,周云波.基于Matlab的DFT及FFT頻譜分析[J].山西電子技術,2010(4):48-49.

[13]梅海青,馮毅.基于振動信號頻譜分析的回轉(zhuǎn)機械故障診斷方法[J].流體機械,2003,31(1):29-31.

Study on the Vibration Condition Monitoring
for Magnetic Levitation Energy Storage Flywheel

XUE Xiaochuan,WANG Zhiqiang
(College of Instrument Science and Optoelectronic Engineering,Beihang University,Beijing 100191,China)

In magnetic levitation energy storage flywheel system,the unbalanced vibration may influence the stability of the rotor and the normal operation of system;aiming at these problems,the reasons of the disturbance are studied.Through establishing the generalized disturbance equation containing unbalance quantity of rotor,and the dynamics equation of magnetic levitation flywheel,the theoretical analysis of the same frequency,double frequency and modal components in the unbalance vibration is carried out.A vibration condition monitoring system based on LabVIEW-DAQ for the magnetic levitation flywheel is designed and implemented.With this system,real-time acquisition,processing and display of the flywheel vibration signal can be completed,and the vibration condition monitoring of flywheel is implemented,and the experimental verification for a certain 300 kW energy storage flywheel is conducted.Aiming at the signal noise caused by the strong electric working environment of flywheel,the original signals are processed with wavelet threshold denosing,and the frequency domain syslysis is accomplished for the vibration signal after denoising,to obtain the vibration spectrum of the flywheel at a certain speed; and the three dimensional waterfall diagram of flywheel vibration signal under different speeds.The experimental results show that in the vibration spectrum of the flywheel,both the stable flywheel rotating frequency component and the double frequency component exist.The system is simple and versatile,so it has good market prospect and application value.

Energy storage flywheel; Magnetic levitation; Unbalance vibration; Frequency component; Data acquisition; Wavelet manic elimination; Frequency domain analysis; Condition monitoring

薛曉川(1993—)，男，在讀碩士研究生，主要從事電機伺服控制技術、集成電路設計、數(shù)字信號處理等方向的研究。 E-mail：sleepingriver@163.com。 王志強(通信作者)，男，博士，講師，主要從事先進空間慣性執(zhí)行機構(gòu)驅(qū)動控制技術、交流永磁電機控制技術、電力電子技術的研究。E-mail：wangzhiqiang@buaa.edu.cn。

TH825;TP806+.1

A

10.16086/j.cnki.issn1000-0380.201707022

修改稿收到日期：2017-01-31