大型風力機系統運行模態分析研究

鐘燦堂， 李德源， 汪顯能， 莫文威， 劉　雄

(1.廣東工業大學 機電工程學院，廣州　510006; 2.汕頭大學 能源研究所，廣東 汕頭　515063)

?

大型風力機系統運行模態分析研究

鐘燦堂1， 李德源1， 汪顯能1， 莫文威1， 劉雄2

(1.廣東工業大學 機電工程學院，廣州510006; 2.汕頭大學 能源研究所，廣東 汕頭515063)

摘要：針對大型風力機在風輪靜止、變速轉動下振動模態及變化特點，研究彈性變形、慣性及陀螺效應引起的系統各階模態變化及對系統氣彈穩定性影響。通過研究現有線性特征值分析方法，考慮大型風力機非線性特性及風輪轉動所致系統時變特性，基于多體系統動力學理論及混合多體系統HMBS(Hybrid Multi-body Systems)建模方法，結合動力學分析軟件ADAMS，分析靜止狀態整機系統線性特征值問題；考慮構件彈性變形及風輪旋轉，用剛性積分方法對系統非線性控制方程進行數值求解，通過傅里葉譜分析方法實現風輪旋轉下系統運轉模態識別，并討論、分析系統前十階模態變化及影響因素。研究結果可作為風力機系統氣彈穩定性判據，為避免共振、提高系統運行效率等提供有效的解決手段及分析方法。

關鍵詞：風力機；超級單元；模態分析；氣彈穩定性

現代風力機柔性構件如葉片、主軸及塔架在運行過程中產生的較大變形對其氣動性能、氣彈穩定性影響日益突出[1-2]。風力機靜止、運行中的模態(包括頻率、振型、阻尼比)及隨風輪轉速、風速變化會直接影響風力機的氣彈穩定性，如葉片相對風輪平面的振動方向對氣動阻尼有重要影響[3]。風輪旋轉及機艙偏航、俯仰運動產生的慣性力、陀螺力等對整個風力機的低階頻率及振形等均影響較大，會改變風力機的氣彈特性。因此準確模擬柔性風力機靜止與風輪轉動中各階模態并理清變化規律，在風力機氣彈穩定研究領域備受關注[4]。

由于風輪轉動及葉片、主軸、塔架等較大彈性變形，使柔性風力機動力學方程呈非線性及質量、剛度分布等結構參數時變性，導致模態分析困難。分析柔性風力機模態分析方法可歸結為兩類[5]：①用一組已知的力函數激勵風力機系統，基于動態響應信號，用系統模態識別方法提取風力機模態信息。Petersen等[6]用風力機氣彈分析軟件HAWC通過在塔頂施加頻率線性變化的正弦偏航及俯仰力矩，分析力矩響應功率譜，獲得某600 kW風力機頻率隨風輪轉速的變化規律及葉片擺振與風輪回旋的相互耦合關系，并通過實測驗證數值分析結果[7]。實測中一般激勵方法無法激起大型風力機振動，Tcherniak等[8-9]提出直接通過運行風力機響應信號提取系統模態方法，但輸入風況的隨機性及非各態歷經性質，其在理論及實際操作尚待進一步研究。此類方法對非線性風力機系統較有效，但僅能捕捉前若干階模態；②直接特征值分析，此需在某個平衡狀態對非線性方程進行線性化建立特征方程，但方程中質量、阻尼及剛度矩陣仍為時變的，一般含風輪轉速周期項，使特征方程建立困難，Bir等[10-11]采用多葉片坐標變換(Coleman 變換)方法消除系統微分方程中的周期項，實現風力機系統模態與氣彈穩定性分析。此類方法能較快速、準確確定各階模態，但對多自由度系統獲得顯式的線性特征方程存在一定困難。多葉片坐標變換僅對完全相同葉片、對稱分布風輪(三葉片風輪)有效。

柔性風力機主要由變形較小的輪轂、機艙及變形較大的葉片、主軸、塔架等構成。通過引入超級單元將柔性構件離散為以運動副及力元相聯的剛體系統，用剛、柔混合多體系統HMBS (Hybrid Multibody Systems)建立系統動力學模型[12]。以5 MW近海風力機[13]為研究對象，結合動力學分析軟件ADAMS建立風力機系統虛擬樣機模型進行動力學方程求解，并分析靜止狀態下風力機系統固有頻率、振型、葉片振動與主軸彎曲變形引起風輪偏航及俯仰運動相互耦合及旋轉風力機動力學方程的時變特性、周期系數處理等。通過在主軸施加角位移驅動提取特定位置響應(如力矩信號)，分析各階頻率、振型隨轉速變化規律及對慣性力、陀螺力影響。研究結果可對準確分析大型風力機模態、確定氣彈失穩邊界、改進風力機設計等均具有重要指導意義。

1風力機混合多體系統模型

借助機械系統動力學分析軟件ADAMS建立柔性風力機剛-柔混合多體模型。其中對柔性構件引入超級單元模型進行離散處理[14]，并據彈性梁彎曲、扭轉理論獲得各超級單元中力元彈簧剛度系數。NREL-5MW近海風力機中連接剛性、柔性子系統彈簧與阻尼器具體數值見文獻[15]。

1.1風力機系統離散

將變形較小的輪轂、機艙處理為剛體。為描述主軸彎曲與扭轉變形，將主軸離散為含5個自由度的超級單元。由于輪轂與主軸第1個剛體無相對運動，將兩者固結。主軸第1剛體與機艙通過球鉸連接模擬輪轂與主軸相對機艙的轉動；高速軸產生的反力矩轉移到主軸第4個剛體，該剛體與機艙通過旋轉鉸連接；3個葉片分別安裝在輪轂上，葉片與輪轂通過轉動鉸連接模擬葉片的變漿運動；機艙與塔架通過球鉸連接，球鉸3個自由度分別模擬機艙相對塔架偏航、俯仰及翻滾運動；塔架底部通過高強度螺栓與支撐基礎剛性連接。由此，風力機系統離散為含56個剛體、95個自由度的剛-柔混合多體系統。在ADAMS中建立整機裝配模型及慣性坐標系，見圖1。

圖1　風力機模型及傳動系統模型Fig.1 Turbine and shaft connection

1.2風力機系統動力學方程與求解

在ADAMS中用某剛體B的質心笛卡爾坐標及反映剛體方位的歐拉角作為廣義坐標，即

q=[x,y,z,ψ,θ,φ]T

(1)

采用不獨立的廣義坐標，系統動力學方程雖數量較大，但卻是高度稀疏耦合的微分代數方程[16]。

考慮約束方程，ADAMS利用帶乘子的第一類拉格朗日方程能量形式，得

(2)

式中：T為系統廣義坐標表達的系統動能；qj為廣義坐標；Qj為廣義力。

式(2)最右端項涉及約束方程及拉格朗日乘子表達式在廣義坐標qj方向的約束反力。該式建立的系統運動微分方程通常為關于N個廣義坐標qj的二階非線性微分-代數方程組，結合初始條件及數值積分法求出方程組的解，此數值解可代表系統的運動規律。對風力機系統而言，各部分運動會出現同時帶快、慢變分量情形，如機艙偏航、俯仰運動及風輪轉速等運動量與葉片揮舞、擺振等振動量變化速度相差較大，通常稱此微分方程為剛性方程，在ADAMS中用BDF(Backwards Differentiation Formulation)剛性積分法，具體求解見文獻[16]。

針對風力機系統氣彈穩定性問題，主要為系統各階模態及變化規律，問題歸結為如何求出系統各階模態的頻率、振型與阻尼比。據導出的系統非線性微分方程推導特征方程較困難，因風力機柔性構件存在振動及彈性變形，各廣義坐標、速度間存在耦合，且風輪轉動時系統是時變的，各葉片在不同時刻處不同位置，即系統質量與剛度分布亦是時變的，故求系統模態較困難。在某特殊情況下如對稱葉片風輪(風輪各葉片假設幾何形狀、力學特性相同，且對稱安裝)的動力學微分方程有其特殊性(周期系統)。小變形假設下可對此非線性微分方程進行線性化并解除各自由度間耦合，引入特定變換方法可獲得時變風力機系統特征方程并進行求解。

具有B個葉片的風力機系統，廣義坐標[17]可表示為

x=[x1,1…x1,Nbx2,1…x2,Nb…xB,1…xB,Nbxs,1…xs,Ns]T

(3)

式中：x第1個下標為葉片數；下標s為系統支撐構件(機艙、主軸、塔架等)在慣性坐標系下坐標；N=BNb+Ns為總變量數，Ns為支撐構件廣義坐標數，Nb為單葉片廣義坐標數。

式(3)中葉片廣義坐標一般用qi,n表達，即

x=[q1,1…q1,Nbq2,1…q2,Nb…qB,1…qB,Nbxs,1…xs,Ns]T

(4)

離散的風力機系統中各剛體動能Ti包括隨剛體質心移動、轉動的動能，即

(5)

式中：mi，Ii分別為某剛體質量及繞質心轉動慣量；vi，ωi分別為某剛體質心速度、角速度，可通過各廣義坐標及對時間導數表達。

廣義力定義為

(6)

將組集系統各剛體動能、廣義力代入式(2)，可得風力機系統動力學微分方程。具體見文獻[18]。

設式(4)中各變量相對未變形狀態為小量，即xi?1，并忽略二階以上高價小量，不引入氣動載荷，則由式(2)得系統線性化運動方程[17]為

(7)

式中：M，C，K分別為質量、阻尼、剛度矩陣，且均為時變矩陣，含T=2π/Ω的周期項，Ω為風輪轉速。

由于矩陣的時變性，即使獲得式(7)的線性化微分方程亦導不出特征方程。考慮靜止條件下，相對于靜平衡位置的運動微分方程式(7)中各矩陣為時不變矩陣，則可進行特征值分析獲得系統各級模態，其中葉片、風力機系統靜止狀態的模態分析即采用此方法，在ADAMS軟件中求出。

風力機在某確定轉速下運行即式(7)中含轉速Ω的周期項可用Lyapunov-Floquet進行分析，但較復雜，應用不便[18]。在直升機、風力機氣彈分析領域更多采用Coleman變換方法消除系統微分方程中周期項，實現風力機系統模態與氣彈穩定性分析。

多葉片變換實質即將在旋轉坐標系下定義的葉片運動與支撐構件運動轉換到慣性坐標下描述，消除控制方程中周期項，使式(7)成為時不變方程，從而導出模態分析的特征方程。將式(4)中葉片模態坐標變換為慣性系下多葉片坐標，即

qi,n=a0,n+a1,ncosΨi+b1,nsinΨi

(8)

由式(8)知，3個多葉片坐標a0,n，a1,n，b1,n分別替換3個葉片模態坐標q1,n，q2,n，q3,n，描述在慣性坐標系中風輪葉片運動。若設葉片模態n為揮舞彎曲模態，則a0,n為所有葉片同時揮舞振動，而a1,n，b1,n則分別描述俯仰、偏航運動。

狀態變量x多葉片變換為

x=B(t)z

(9)

式中：z為狀態向量，含多葉片坐標，即

(10)

轉換矩陣B為

(11)

(12)

若轉速Ω恒定，則B為常數陣，即

(13)

(14)

(MBR2+CBR+KB)z=0

(15)

式中：MB≡μBTMB，CB≡μBTCB，KB≡μBTKB為變換后的系統矩陣。

式(15)中如果風輪對稱，矩陣則為常數陣，即多葉片坐標下式(15)已消除周期項，變成時不變方程，則可定義風力機某轉速的特征值問題。

2靜止狀態下葉片與系統模態分析

靜止狀態指風力機系統處于靜平衡位置，風輪無旋轉且不考慮氣動載荷。采用線性化模型對葉片及系統進行特征值分析。

風力機系統前十階模態對應頻率見表1。表中NREL的FAST及ADAMS分析結果差異因FAST用16個自由度整機模型，而ADAMS用438個自由度整機多體系統模型，本文用ADAMS軟件分析，較NREL分析結果偏小，主要因風力機多體模型建立與NERL有所不同，尤其對主軸、葉片及塔筒等柔性構件處理引入超級單元概念。此外，NREL只用1個自由度分別描述主軸相對機艙運動、機艙相對塔架運動，而本文則用球鉸3個自由度分別描述兩種運動。

表1　風力機系統靜止時前十階頻率

各階模態對應振型見圖2。由圖2看出，第1、2階分別為塔架前后與側向彎曲模態，前后彎曲模態頻率略低于側向彎曲，此因該模態含部分風輪俯仰運動，風輪、機艙質量加大塔架該模態慣性所致。風輪、主軸的柔性對兩階模態影響較小。第3階模態為傳動鏈扭轉模態，即輪轂與發電機間相對扭轉(圖2(a))。該模態耦合三葉片同相擺振變形，其頻率值受發電機控制方式影響，恒速控制與異步發電機采用固定扭轉邊界條件建模，該模態與塔架側向彎曲相耦合，使該模態頻率低于塔架的彎曲頻率；而恒扭矩或功率控制型則采用自由扭轉邊界條件，模態頻率較塔架彎曲模態頻率高(表1)。

第4、5階為第1階偏航/俯仰模態對，與葉片揮舞模態相耦合，為兩階反對稱模態，即在偏航模態下振動時葉片1不動，葉片2、3的揮舞方向相反，此為葉片揮舞與機艙偏航運動耦合所致；而俯仰振動時葉片2、3的揮舞方向相同，但與葉片1相反，此為葉片揮舞與機艙俯仰運動耦合結果。第6階模態是第1階對稱揮舞模態，振動時3個葉片同相揮舞，葉片揮舞與塔架前后彎曲耦合，但兩種振動反相(圖2(b))。此耦合造成該模態頻率(0.688 5 Hz)略高于葉片1階揮舞頻率(0.674 5 Hz)。

第7、8階模態與葉片擺振模態相耦合，即兩階反對稱模態，且非常接近，約等于第1階葉片擺振模態頻率(圖2(c))。7階模態中葉片1不動，葉片2、3擺振方向相反但振幅相同；8階中1與2、3擺振方向相反，振幅為兩葉片的2倍。葉片擺振慣性力對風輪主軸力矩相互平衡，即各葉片擺振慣性力對主軸取矩之和等于零，該兩階模態與主軸扭轉無耦合；但各葉片擺振慣性力合力不為零，7階模態慣性力合力沿豎直方向，稱1階垂直擺振；8階模態慣性力合力沿水平方向，稱1階水平擺振，并與主軸彎曲變形耦合。文獻[7]研究表明，葉片的氣動阻尼與葉片振動方向相關，揮舞方向較擺振方向大，兩階模態與系統氣彈穩定性關聯較大，更易引起氣彈失穩。設計風力機葉片時可適當調整葉片各截面質量、剛度及扭轉中心等結構參數分布，改善葉片氣彈特性。

第9、10階模態為2階偏航/俯仰模態，機艙偏航與俯仰運動與第2階葉片揮舞模態耦合(圖2(d)) (文獻[18]針對某600 kW風力機模態分析，兩階模態仍與葉片第1階揮舞模態耦合，僅風輪偏航、俯仰與機艙偏航與俯仰反向)。

圖2　系統前十階模態振型Fig.2 The first ten order modes of the tuebine system

3運行狀態下葉片與系統模態分析

現代風力機柔性構件較大變形及運動對其動力特性影響不可忽略，因控制方程為非線性方程，要獲得系統模態需對系統施加特殊激勵，通過求解非線性動力學方程獲得系統時域響應，分析時-頻獲得系統各階頻率。分析時不考慮氣動載荷，對轉動軸最右邊與機艙用旋轉副連接的剛體(圖1)施加角位移驅動函數，使風輪從靜止狀態穩步旋轉到轉速ω1產生強迫振動。角位移驅動函數為

s=ωt=0, t=0)ω1t220, (0

(16)

式中：s為角位移；ω為瞬時角速度；ω1為輸入角速度；t為時間(0～120 s)，求解步長0.01 s。

在ADAMS所建模型上通過求解非線性動力方程，獲得轉速n=3 r/min、n=6 r/min、n=9 r/min、n=12 r/min的時域響應，選某特殊位置響應信號進行頻譜分析，可得各階模態對應頻率。選時域信號時，1、2階選塔架近地面的萬向節彎曲力矩響應；第3階選轉動軸中間轉動副力矩響應；第4～10階先從塔頂與機艙連接的球鉸力矩獲取俯仰、偏航力矩響應(因輪轂連接3葉片旋轉時模態會耦合到塔頂的球鉸副力矩中)，再從葉片近輪轂連接的萬向節揮舞方向力矩響應識別出4、5、6、9、10階固有頻率，從擺振方向力矩響應識別出7、8階固有頻率。

塔頂處俯仰力矩在風輪額定轉速n=12 r/min下時域與頻率響應見圖3。

圖3　塔頂俯仰力矩響應及頻譜圖Fig.3 The tilting moment in top of the tower and the relative spectrogram

風力機不同轉速的前十階頻率與對應振型說明見表2，各階頻率隨轉速變化曲線即Campell圖見圖4。由表2、圖4看出，前二階模態頻率幾乎不隨轉速變化，說明風輪運動與塔架彎曲變形耦合較小。

表2　風力機旋轉的前十階頻率

圖4　系統前十階頻率隨轉速變化曲線Fig.4 The Campbell diagram of the turbine system

風輪旋轉狀態下第4、5階模態為葉片第1階揮舞模態與風輪向后或向前回旋模態耦合，而第7、8階模態則為葉片第1階擺振模態與風輪向后回旋、向前回旋模態耦合。此兩對頻率隨轉速變化明顯，變化規律接近ωk+Ω及ωk-Ω，ωk為葉片揮舞或擺振頻率。此因風輪旋轉時機艙有偏航、俯仰偏轉運動，即風輪、機艙動坐標系作旋轉運動，出現哥氏加速度即風輪旋轉的陀螺效應引起反對稱模態頻率變化；而向后、向前回旋模態可用線性分析結果定性說明。葉片模態坐標[17]為

(17)

應特別關注該陀螺效應引起各階頻率變化對系統氣彈穩定性影響。文獻[7]研究表明，葉片氣動阻尼與振動方向相關，揮舞方向氣動阻尼較擺振方向大，風輪旋轉時向前回旋的氣動阻尼較向后回旋大。若系統某階模態氣動阻尼加結構阻尼為負，則系統發生氣彈失穩。由圖4知，第3、6階模態頻率在轉速3 r/min時發生交叉，此時揮舞與擺振耦合較嚴重，尤其額定轉速時第5、7階固有頻率較接近，因在揮舞振動中混入更多擺振，第8、9階亦呈接近趨勢。此為分析特定風力機系統氣彈失穩特性提供了，即設計風力機系統時在保證系統氣動性能基礎上可有意提高某階模態揮舞分量或前旋分量，使其氣動阻尼更大，保證系統不發生氣彈失穩。

第3階主軸扭轉模態與風輪對稱擺振耦合、第6階模態為對稱風輪揮舞模態，均無陀螺效應，頻率隨轉速略有增加，此由葉片轉動引起離心剛化效應所致。由于非線性模型中考慮葉片較大變形與振動耦合，使剛化效應出現的同時伴隨弱化，故兩頻率隨轉速增加不大，尤其擺振頻率高于揮舞頻率，弱化效應更大，頻率增加較揮舞小。

4結論

針對大型風力機系統模態分析方法及模態耦合對系統氣彈穩定性影響進行研究。結合動力學分析軟件ADAMS將風力機系統離散為剛-柔混合模型，實現風力機葉片及整機靜態線性特征值分析及風輪轉動下的非線性模態分析，結論如下：

(1)葉片各階模態相互耦合，葉片振動受系統其它構件如塔架、機艙、主軸動力特性影響較大，分析系統氣彈振動尤其失速誘導的葉片振動時須從整機動特性出發。

(2)風輪旋轉時出現葉片揮舞或擺振模態與風輪向前、后回旋模態耦合，風輪旋轉所致陀螺效應引起各階反對稱模態頻率隨風輪轉速分離，旋轉所致剛化效應因伴隨弱化效應，使對稱模態頻率隨轉速變化不大。

(3)系統各部變形相互耦合及慣性載荷與陀螺載荷會影響系統各階模態氣動阻尼，從而強烈影響系統氣彈穩定性。本文方法可為運行風力機氣彈穩定性判據及新風力機設計提供參考。

參 考 文 獻

[1] Riziotis V A, Voutsinas S G, Politis E S,et al. Aeroelastic stability of wind turbines: the problem, the methods and the issues[J]. Wind Energy, 2004(7):373-392.

[2] Meng Fan-zhong. Aero-elastic stability analysis for large-scale wind turbines[D]. Delft: Delft University of Technology, 2011.

[3] Kalles?e B S, Hansen M H. Some effects of large blade deflections on aeroelastic stability[C]//47th AIAA Aerospace Sciences Meeting.Florida, USA, 2009.

[4] Hansen M H. Aeroelastic instability problems for wind turbines[J]. Wind Energy, 2007 (10):551-557.

[5] Bir G, Jonkman J. Aeroelastic instabilities of large offshore and onshore wind turbines[C]//The EAWE 2007 Torque from Wind Conference Technical University of Denmark. Lyngby, Denmark, 2007.

[6] Petersen J T. Thomsen K, Madsen H A. Local blade whirl and global rotor whirl interaction[R]. Technical Report Ris?-R-1067, Ris? National Laboratary, Roskilde, 1998.

[7] Thomsen K, Petersen J T, Nim E, et al. A method for determinatoin of damping for edge blade vibrations[J]. Wind Energy, 2001(3):233-246.

[8] Tcherniak D, Chauhan S, Rossetti I, et al. Output only modal analysis on operating wind turbines: application to simulated data[C]//European Wind Energy Conferece.Warsaw, Poland, 2010.

[9] Allen M S, Sracic M W, Chauhan S, et al. Output only modal analysis of linear time periodic systems with application to wind turbine simulation data[J]. Mechanical Systems and Signal Processing, 2011(25):1174-1191.

[10] Bir G. Multiblade coordinate transformation and its application to wind turbine analysis[C]//The 2008 ASME Wind Energy Symposium. Nevada, USA, 2008.

[11] Hansen M H. Aeroelastic stability analysis of wind turbines using an eigenvalue approach[J]. Wind Energy, 2004(7):133-143.

[12] Zhao X Y Peter M, Wu J Y. A new multibody modelling methodology for wind turbine structures using a cardanic joint beam element[J]. Renewable Energy, 2007(32): 532-546.

[13] Jonkman J, Butterfield S, Musial W, et al. Definition of a 5-MW reference wind turbine for offshore system development[R]. National Renewable Energy Laboratory, 2009.

[14] Holierhoek J G. Aeroelasticity of large wind turbines[D]. Delft: Delft University of Technology, 2008.

[15] 李德源,莫文威,嚴修紅,等. 基于多體模型的水平軸風力機氣彈耦合分析[J].機械工程學報,2014,50(12):140-150.

LI De-yuan,MO Wen-wei,YAN Xiu-hong,et al.Aeroelastic analysis of horizontal axis wind turbine based on multi-body model[J].Journal of Mechanical Engineering,2014,50(12):140-150.

[16] MSC Software, Adams 2013/help[C].2013 .

[17] Hansen M H. Improved modal dynamics of wind turbines to avoid stall-induced vibrations[J]. Wind Energy, 2003(6): 179-195.

[18] Skjoldan P F, Hansen M H. On the similarity of the Coleman and Lyapunov-Floquet transformations for modal analysis of bladed rotor structures[J]. Journal of Sound and Vibration, 2009(327):424-439.

Operational modal analysis of a large wind turbine

ZHONGCan-tang1,LIDe-yuan1,WANGXian-neng1,MOWen-wei1,LIUXiong2

(1. School of Electromechanical Engineering, Guangdong University of Technology, Guangzhou 510006, China;2. Institute of Energy Science, Shantou University，Shantou 515063，China)

Abstract:Aiming at the vibration modes and its variation characteristics of a large wind turbine under static status and various rotational speeds, the variations of each mode caused by elastic deformation, inertia and gyroscopic effect and their effects on aeroelastic stability of the wind turbine system were studied. Considering the large wind turbine system’s non-linear characteristics and time-varying characteristics caused by the rotational wind rotor, the linear eigenvalue problems of the wind turbine were analyzed by using the dynamics software ADAMS, based on the dynamics theory of multi-body system and the modeling methodology of hybrid multi-body system (HMBS). A stiff integral method was then employed to solve the non-linear control equations of the system numerically considering the elastic deformation of flexible components and the rotation of the wind rotor. The operational modal identification was performed through Fourier spectral analysis. The variations of the system’s first ten modes and their influence factors were discussed and analyzed in detail. The results can be used as aeroelastic stability criteria, which offer an effective solution and analysis method to avoid resonance and improve operating efficiency of the wind turbine system.

Key words:wind turbine; superelement; modal analysis; aeroelastic stability

中圖分類號:TK83

文獻標志碼：A

DOI:10.13465/j.cnki.jvs.2016.06.022

通信作者李德源 男，博士后，教授，1965年生

收稿日期：2014-12-17修改稿收到日期：2015-03-20

基金項目：國家自然科學基金(51276043)；教育部高等學校博士點科研基金項目(20124402110005)；教育部科學技術研究重點項目(212130)資助

第一作者 鐘燦堂 男，碩士生，1988年生

E-mail：lidey@gdut.edu.cn