Spar型海上浮式風力機系統(tǒng)運動耦合計算方法

尚景宏,趙玉娜，張亮,胡長洪，丁曉朦

(1.哈爾濱工程大學 深海工程技術(shù)研究中心，黑龍江 哈爾濱 150001；2.中海石油(中國)有限公司 天津分公司，天津 300452;3.日本九州大學 應用力學研究所，福岡縣 春日 816-8580)

?

Spar型海上浮式風力機系統(tǒng)運動耦合計算方法

尚景宏1,2,趙玉娜1，張亮1,胡長洪3，丁曉朦1

(1.哈爾濱工程大學 深海工程技術(shù)研究中心，黑龍江 哈爾濱 150001；2.中海石油(中國)有限公司 天津分公司，天津 300452;3.日本九州大學 應用力學研究所，福岡縣 春日 816-8580)

針對海上浮式風力機系統(tǒng)(FOWTs)與風浪流相互作用問題，建立了時域非線性耦合動力學方程和流體動力數(shù)學模型，采用Fortran語言編程、時域頻域轉(zhuǎn)換和龍格-庫塔法迭代求解方法，分析定常風、規(guī)則波和隨機風浪環(huán)境下，Spar型FOWTs除艏搖之外的5自由度運動響應特性。結(jié)果表明，額定定常風速時，風力機系統(tǒng)5自由度運動響應的平均值和峰峰值都受影響，而其他定常風速不影響縱向運動響應的峰峰值；規(guī)則波對5自由度運動響應的平均值影響很小，但峰峰值明顯增大；在隨機風浪下，F(xiàn)OWTs的縱蕩(搖)和升沉運動分別不同程度地出現(xiàn)明顯低頻性和波頻性，而且與風速相關(guān)。分析方法及其結(jié)果為海上FOWTs的運動性能設計提供參考。

Spar型海上浮式風力機系統(tǒng)；時域非線性耦合；運動響應特性；規(guī)則波；隨機風浪

近年來，能源開發(fā)與環(huán)境保護的需求強力推動著風能科學與技術(shù)的研究，產(chǎn)業(yè)規(guī)模和新技術(shù)不斷發(fā)展，2009～2014年全球風電裝機容量的年均增長率高達18.47%[1]。

海上風能因分布廣泛、發(fā)電穩(wěn)定等優(yōu)勢而受到重視，海上風場建設和裝機容量年均增長率大于風電總裝機增長率，發(fā)展十分迅速[1-3]。海上風力機按照支撐結(jié)構(gòu)可分為固定式和浮式兩種，其中浮式風力機(floating offshore wind turbine,FOWT)可作業(yè)于深水海域，是海上風力機的主要發(fā)展方向[3-5]。根據(jù)FOWT基礎的形式，可進一步分為Spar式、張力腿式以及半潛式等。其中，Spar平臺是浮式風力機支撐結(jié)構(gòu)中應用較廣泛的形式，其結(jié)構(gòu)設計可使整個系統(tǒng)重心降低[6]，在各種工況下[7-8]仍能保持良好穩(wěn)性；另外，Spar式風力機在風浪環(huán)境下載荷與運動較小，經(jīng)濟性比張力腿式有很大的優(yōu)勢[7]。目前對于Spar型FOWTs的研究已經(jīng)取得了顯著成果，在運動響應、氣動性能等方面發(fā)展了一些理論模型，在線性頻域內(nèi)進行系統(tǒng)性能預報[7-8]，并且開展相關(guān)的試驗技術(shù)研究。2009年，世界上第一座漂浮式風力機Hywind在挪威附近的北海正式啟用。然而，F(xiàn)OWTs在風浪中的運動是耦合的、非線性的。就Spar式平臺本身已發(fā)展了非線性理論預報方法[9-11]。Hanson等[16-18]把HAWC2和SIMO/RIFLEX集成在一起，計算了FOWT平臺在各種風浪工況下的性能，并與模型試驗結(jié)果對比。Shim等用WAMIT計算NREL張力腿式風力機系統(tǒng)的水動力性能，用有限元方法分析錨泊系統(tǒng)載荷和強度，對風力機、浮體和錨泊系統(tǒng)耦合運動進行頻域、時域分析。但對于平臺上安裝大型風力機的FOWTs，在風浪流載荷下的耦合運動研究主要限于線性簡化模型[12]，其耦合特性規(guī)律和機理并不明確。本文對Spar型海上FOWTs建立時域非線性運動模型，應用時域耦合以及時頻轉(zhuǎn)換的方法研究定常風、規(guī)則波對系統(tǒng)運動響應的影響以及該系統(tǒng)在隨機風浪下的運動響應特性。

1　計算方法

考慮Spar型FOWTs(含風力機組、浮式平臺及其錨系)在風浪流作用下的運動。假設平臺為剛體，平臺瞬時姿態(tài)由平臺歐拉角Θ=(α,β,γ)T、重心瞬時位置X0=(x,y,z)T、重心運動線速度U0=(ux,uy,uz)T及角速度Ω=(ωx,ωy,ωz)T來描述。矢量形式的系統(tǒng)非線性動力學方程[10]組寫為

(1)

式中：M、I分別為系統(tǒng)的質(zhì)量和轉(zhuǎn)動慣量矩陣；F、N分別為作用于重心處的合外力和力矩，表示為

(2)

式中：下標分別表示作用于平臺重心處的靜水回復力、波浪力、流力、風力(含風輪)和系泊力及其它們的力矩。式(1)中B是平臺繞x、y、z軸的瞬時角位移和角速度間的關(guān)系矩陣，由平臺歐拉角表示為[10]

(3)

關(guān)于風力機氣動載荷數(shù)值計算模型，分別基于葉素理論、廣義動態(tài)尾流法[14]計算葉元誘導速度，基于B-L(Beddoes-Leishman)動態(tài)失速理論[15]修正葉元氣動力系數(shù)，在此基礎上計算風輪氣動載荷。廣義動態(tài)尾流法模擬偏航和動態(tài)入流效應精度高，而B-L動態(tài)失速理論可修正動態(tài)入流引起的動態(tài)失速，以進一步提高風力機氣動載荷的計算精度，是較適合FOWT復雜氣動載荷的數(shù)值模型?？紤]到Spar平臺的結(jié)構(gòu)特征，作用于平臺的水動力載荷由細長體理論計算，系泊線載荷采用準靜態(tài)懸鏈線法計算。此外，考慮了風力機的變槳控制，其控制模式與Hywind OC3相同。

在方程組(1)中，各物理量是時間的函數(shù)，各方程的參數(shù)之間相互關(guān)聯(lián)，而且待求物理量存在乘積項，因此該方程組為非線性的耦合運動方程，需要在時域中迭代求解。為方便數(shù)值計算，將方程(1)轉(zhuǎn)化為

(4)

方程(4)為二階常微分方程。方程的求解采用四階龍格-庫塔法[12]，取時間步長Δt，在時域中步進迭代，直至收斂至設定精度?；谏鲜鼋o出的FOWTs方程組、流體動力數(shù)值模型及其計算方法，編制了Fortran時域耦合程序。數(shù)值計算流程如圖1所示。

圖1　數(shù)值計算流程Fig.1　Flowchart of the numerical model

2　耦合程序驗證

美國NREL設計的OC3-Hywind[13]浮式風力機系統(tǒng)的支撐平臺屬典型的Spar型。將該FOWTs作為計算模型，驗證本文時域耦合程序的正確性。取工況1(環(huán)境參數(shù)見表1)運動響應計算結(jié)果，與文獻[14]所給數(shù)據(jù)進行對比分析，如圖2所示。

表1　工況1環(huán)境參數(shù)

文獻[14]中用HAWC2、Accoina SESAM等5種不同計算工具進行模擬分析，不同計算工具得出的計算結(jié)果有一定的差異。這里取文獻中計算結(jié)果峰值最大值“上限”、谷值最小值“下限”以及一組有代表性的曲線“文獻值”與本文計算結(jié)果“計算值”對比。

圖2　工況1 FOWTs響應時歷曲線Fig.2　Response history of the FOWTs in case I

從圖2可以看出，本文耦合模型計算出的縱蕩、升沉、艏搖運動響應都在文獻計算值的上下限范圍內(nèi)，且運動響應曲線趨勢與文獻中有代表性的運動響應曲線變化趨勢一致；在風輪氣動響應方面，本文風輪轉(zhuǎn)速計算值峰值比文獻風輪轉(zhuǎn)速峰值上限大0.5 r·m-1，使風輪功率峰值也超出文獻峰值。總體來說，雖然本文數(shù)值模型計算的結(jié)果與文獻值略有差異，但得出的風力機系統(tǒng)各項性能曲線與文獻中性能曲線走勢相一致，在FOWTs耦合運動分析和選型設計上具有應用和參考價值。

3　定常風對FOWTs運動響應的影響

在工作工況下，風力機的氣動推力會使FOWTs產(chǎn)生縱蕩和縱搖運動，風力機的扭矩會使系統(tǒng)產(chǎn)生橫搖運動；同時，在橫搖和縱搖運動的作用下，上部風力機會偏航和俯仰，產(chǎn)生偏航力矩和俯仰力矩，從而引起系統(tǒng)艏搖，因此，定常風會對FOWTs帶來多方面影響。選取3種工況(表2)對定常風作用下FOWTs的水/氣動力響應進行計算和分析。

表2　定常風工況環(huán)境參數(shù)

除艏搖的5自由度運動響應平均值及峰峰值列于圖3，風輪推力、扭矩和功率的平均值和峰峰值列于圖4。在該FOWTs設計時，為避免艏搖運動過大影響風力機性能，加入了艏搖彈簧抑制艏搖運動，致系統(tǒng)艏搖運動很小(約10-3度數(shù)量級)，因此，不對該系統(tǒng)的艏搖運動作詳細分析。

對于峰峰值，從圖3(b)、(d)、(f)看出，在8 m/s(較低)和16 m/s(較高)風速下，系統(tǒng)的縱向(縱蕩、升沉和縱搖)運動均可平衡到某個運行點，運動響應的峰峰值幾乎為零；而在額定風速11.2 m/s下，系統(tǒng)的縱向運動呈周期性變化，峰峰值分別為6.77 m、2.06 m、14°。此外，11.2 m/s風速下系統(tǒng)橫向(橫蕩和橫搖)運動響應峰峰值也明顯比8 m/s及16 m/s風速下的峰峰值大很多(圖3(h)、(j))。這是因為在額定風速附近風力機推力變化較為劇烈(圖4(b))，從風力機的穩(wěn)態(tài)推力曲線來說，在低于額定風速時推力隨風速的增加而增大，而在高于額定風速時推力隨風速的增加而減小，因此，在額定風速附近時，風力機系統(tǒng)的整體響應很難平衡在某一點，而是呈現(xiàn)出隨時間周期性變化的特點。

對于平均值，風力機系統(tǒng)在11.2 m/s風速下縱向運動平均值較大(圖3(a)、(c)、(e))。這主要歸因于風力機的氣動推力(圖4(a))，從風力機穩(wěn)態(tài)推力曲線也可進行解釋，風力機在額定風速11.2 m/s附近的推力明顯大于8 m/s和16 m/s風速附近的推力，推力大導致系統(tǒng)的縱向運動響應值變大。另外，橫蕩和橫搖的平均值隨著風速的增加而不斷增大(圖3(g)、(i))，這歸因于風力機的氣動扭矩(圖4(c))，且從穩(wěn)態(tài)扭矩曲線來講，低于額定風速時扭矩隨風速的增加而增大，高于額定風速時變化較小，因此，風力機的扭矩平均值在8 m/s附近時最小、16 m/s附近時最大，從而導致在16 m/s風速時系統(tǒng)的橫蕩和橫搖運動均值最大。

圖3　工況1～3 FOWTs運動響應均值和峰峰值對比Fig.3　Comparison of motion in case 1～3

圖4　工況1～3風輪氣動載荷隨風速變化Fig.4　Aerodynamic loading of the turbine rotor in case 1～3

4　規(guī)則波對FOWTs運動響應的影響

從FOWTs風輪的推力、扭矩以及功率穩(wěn)態(tài)曲線可知，在低風速、額定風速和高風速等不同風速區(qū)域，風力機氣動性能有明顯的差別。著重分析同一規(guī)則波中不同風速下風力機氣動性能以及系統(tǒng)運動性能的影響，三種計算工況環(huán)境參數(shù)見表3。

表3　工況4～6環(huán)境參數(shù)

工況4～6與工況1～3下風力機系統(tǒng)的5自由度運動響應結(jié)果見圖5，風力機氣動載荷結(jié)果見圖6。從圖6看出，規(guī)則波對FOWT平臺的縱蕩、縱搖、橫蕩、橫搖運動響應平均值的影響較小，這是因為這4個自由度的運動響應主要歸因于風輪的氣動推力和扭矩(圖6(a)、(b))。

從圖5、圖6看出，在定常風和規(guī)則波作用下，浮式風力機系統(tǒng)的各自由度運動響應、葉輪氣動載荷響應峰峰值較無浪工況都變大。

圖5　工況4～6FOWTs運動響應均值和峰峰值對比Fig.5　Comparison of FOWT’s motion in case 4～6

圖6　工況4～6風輪氣動載荷隨風速變化Fig.6　Aerodynamic loading of the turbine rotor in case 4～6

圖7和圖8分別給出了FOWTs風輪偏航角、俯仰角和偏航力矩、俯仰力矩在各風速下的平均值和峰峰值計算結(jié)果。

從圖7(a)、(b)看出，在8 m/s和16 m/s風速下，風輪的偏航角平均值幾乎為零，偏航角的峰峰值在0.1°～0.2°范圍內(nèi)，可認為風力機風輪不偏航；而在11.2 m/s風速下，偏航角平均值為-0.02°、峰峰值為1°左右，雖然比其他兩個風速下的平均值和峰峰值大很多，但該角度仍很小，也認為不偏航。

從圖7(c)、(d)看出，在8、16 m/s風速下，風輪的俯仰角平均值基本在3°左右，而11.2 m/s風速下為4.5°，約為前者的1.5倍；11.2 m/s風速下風輪俯仰角的峰峰值高達15°，而其它兩個風速下均低于2°。

綜上可見，該FOWTs風輪的俯仰明顯比偏航嚴重很多，這從圖8所示的力矩進一步說明。在各風速下風力機偏航力矩平均值的絕對值均小于4 kN·m，峰峰值不超過50 kN·m；而俯仰力矩平均值最大可達180 kN·m，峰峰值高達500 kN·m。

風力機俯仰會明顯影響風輪的氣動性能，同時俯仰角峰峰值較大，還會明顯影響風力機的疲勞性能，因此，應該針對風輪的俯仰性能，對FOWTs進行仔細的優(yōu)化設計。

圖7　工況4～6FOWTs風輪偏航和俯仰角均值與峰峰值Fig.7　Yaw and tilt angle of the FOWTs in case 4～6

圖8　工況4～6FOWTs風輪偏航和俯仰力矩均值與峰峰值Fig.8　Yaw and tilt moment of the FOWTs in case 4～6

5　湍流風不規(guī)則波下FOWTs運動響應特性分析

從3、4節(jié)可知，在低于和高于額定風速區(qū)域，定常風和規(guī)則波對Spar型FOWTs運動響應特性的影響規(guī)律基本一致。因此，對湍流風、不規(guī)則波下FOWTs運動響應特性的分析，只取額定風速11.2 m/s(工況7)和高風速區(qū)域的代表性風速16.0 m/s(工況8)，參數(shù)見表4。

表4　工況參數(shù)

圖9為工況7、8的歸一化后的湍流風譜和波浪譜?？梢钥闯?，本文湍流風數(shù)值模型模擬出的風譜能與目標風譜基本吻合(圖9(a)、(b))，不規(guī)則波浪數(shù)值模型模擬出的波浪譜也很好地與目標波浪譜吻合(圖9(c))，表明了本文湍流風以及不規(guī)則波數(shù)值模型的正確性。

圖9　工況7、8環(huán)境風浪譜Fig.9　Wind and wave spectra for case 7、8

圖10　工況7、8FOWTs縱蕩(搖)、升沉運動響應譜Fig.10　Surge, pitch and heave spectra for case 7、8

圖10為工況7、工況8下Spar式FOWTs縱蕩、縱搖及升沉運動歸一化后的響應譜。從圖10(a)、(b)看出，在額定及高風速下FOWTs縱蕩運動均呈現(xiàn)出明顯的低頻特性，這說明湍流風的低頻特性對FOWTs的縱蕩運動影響明顯。從圖10(c)、(d)看出，在額定和高風速區(qū)域FOWTs的縱搖運動響應譜特性有明顯不同，在額定風速處，系統(tǒng)的縱搖運動呈現(xiàn)出明顯的低頻性，而在高風速區(qū)域，系統(tǒng)的縱搖運動呈現(xiàn)出明顯的波頻特性。從圖10(e)、(f)看出，與額定風速處相比，高風速區(qū)域FOWTs升沉運動響應波頻特性明顯增強。這主要是因為在額定風速附近，風輪的氣動推力變化較為劇烈，對FOWTs運動響應影響較大，而湍流風呈現(xiàn)出顯著的低頻特性，故致使額定風速處的運動響應低頻特性要比高風速區(qū)域的運動響應低頻特性明顯。

6　結(jié)論

1)考慮浮式風力機系統(tǒng)在風浪流聯(lián)合作用下的瞬時受力與運動，建立了FOWTs的時域非線性耦合動力學方程組、各部件的流體動力數(shù)學模型、以及在時域中迭代求解耦合方程的方法，編寫了Fortran程序源代碼。針對OC3-Hywind系統(tǒng)，與已有結(jié)果對比驗證了數(shù)學模型和程序代碼的有效性，能夠應用于典型Spar式FOWTs耦合性能分析。

2)定常風影響FOWTs運動響應的程度與風速范圍有關(guān)。在額定風速處，定常風對5自由度運動響應的平均值和峰峰值都有顯著影響。而在低于和高于額定風速的范圍，僅對5自由度運動響應的均值和橫蕩(搖)運動的峰峰值產(chǎn)生一定影響。

3)規(guī)則波對FOWTs5自由度運動響應的平均值影響很小，但峰峰值明顯增大。風浪環(huán)境對風輪推力、扭矩和功率的平均值和峰峰值的影響特性分別與其對縱蕩(搖)、升沉的影響特性一致；從穩(wěn)態(tài)推力曲線的不同區(qū)域?qū)Ρ龋L浪環(huán)境在額定風速處對FOWTs運動響應和風輪氣動性能影響最大；此外，在定常風、規(guī)則波作用下，風輪俯仰運動劇烈，影響風力機氣動及疲勞性能，需針對俯仰性能對FOWTs進行優(yōu)化設計。

4)隨機風浪環(huán)境對FOWTs的縱蕩運動響應呈現(xiàn)出明顯的低頻性；系統(tǒng)的縱搖運動響應在額定風速處低頻性占主導地位，而在高于額定風速區(qū)波頻性占主導地位；在額定及高風速區(qū)，系統(tǒng)升沉運動響應的低頻性和波頻性都較為明顯，在高風速區(qū)的波頻性比額定風速處明顯增強。

[1]徐濤. 2014年世界風能產(chǎn)業(yè)概況[J]. 風能產(chǎn)業(yè), 2015(10): 22-28.

[2]SAHIN A D. A review of research and development of wind energy in Turkey[J]. Clean-soil, air, water, 2008, 36(9): 734-742.

[3]TWIDELL J, GAUDIOSI G. Offshore wind power[M]. (s.l.): Multi-science Publishing Co. Ltd, 2009: 2-14.

[4]張亮, 吳海濤. 海上浮式風力機技術(shù)現(xiàn)狀及難點[J]. 風能產(chǎn)業(yè), 2013(7): 22-28.

[5]WANG Yi, DUAN Menglan, SHANG Jinghong. Application of an abandoned jacket for an offshore structure base of wind turbine in Bohai heavy ice conditions[C]//Proceedings of the 19th International Offshore and Polar Engineering Conference. Osaka, Japan, 2009: 384-389.

[6]WITHEE J E. Fully coupled dynamic analysis of a floating wind turbine system[D]. Cambridge, MA, USA: Massachusetts Institute of Technology, 2004.

[7]KARIMIRAD M. Dynamic response of floating wind turbines[J]. Scientia Iranica, 2010, 17(2): 146-156.

[8]張亮, 趙玉娜, 馬勇, 等. Spar式風力機平臺設計及水動力影響因素研究[J]. 哈爾濱工程大學學報, 2015, 36(1): 19-23.

ZHANG Liang, ZHAO Yu’na, MA Yong, et al. Design of a Spar-type wind turbine platform and analysis of the hydrodynamic influencing factors[J]. Journal of Harbin engineering university, 2015, 36(1): 19-23.

[9]MA Q W, PATEL M H. On the non-linear forces acting on a floating spar platform in ocean waves [J]. Applied ocean research, 2001, 23(1): 29-40.

[10]MA Q W. Numerical simulation of nonlinear interaction between structures and steep waves [D]. London: University of London, 1998.

[11]CHITRAPU A S, SAHA S, SALPEKAR V Y. Time-domain simulation of spar platform response in random waves and current[C]//Proceedings of the 17th International Conference on Offshore Mechanics and Arctic Engineering. Lisbon Portugal, 1998.

[12]SHI Wei, PARK H C, CHUNG C W, et al. Time domain and frequency domain characterization of floating offshore wind turbine[C]//The Proceedings of the 9th ISOPE Pacific/Asia Offshore Mechanics Symposium, PACOMS-2010. Busan, Korea, 2010: 91-98.

[13]JONKMAN J. Definition of the floating system for phase IV of OC3[R]. Technical Report NREL/TP-500-47535. National Renewable Energy Laboratory, 2010.

[14]JONKMAN J, MUSIAL W. Offshore code comparison collaboration (OC3) for IEA task 23 offshore wind technology and deployment[R]. (s.l.):Technical Report NREL/TP-500-48191. National Renewable Energy Laboratory, 2010.

[15]HANSEN A C, BUTTERFIELD C P. Aerodynamics of horizontal-axis wind turbines[J]. Annual review of fluid mechanics, 1993, 25: 115-149.

[16]NIELSEN F G, HANSON T D, SKAARE B. Integrated dynamic analysis of floating offshore wind turbines[C]//Proceedings of 25th International Conference on Offshore Mechanics and Arctic Engineering. Hamburg, Germany,2006.

[17]LARSEN T J, HANSON T D. A method to avoid negative damped low frequent tower vibrations for a floating, pitch controlled wind turbine [J]. Journal of physics: conference series, 2007, 75(1): 012073.

[18]SKAARE B, HANSON T D, NIELSON F G. Importance of control strategies on fatigue life of floating wind turbines[C]//Proceedings of the 26th International Conference on Offshore Mechanics and Arctic Engineering. San Diego, CA, United States, 2007.

本文引用格式：

尚景宏,趙玉娜，張亮，等. Spar型海上浮式風力機系統(tǒng)運動耦合計算方法[J]. 哈爾濱工程大學學報， 2016, 37(9): 1163-1171.

SHANG Jinghong, ZHAO Yuna, ZHANG Liang，et al. Coupled method for predicting motions of Spar-type offshore floating wind turbine systems[J]. Journal of Harbin Engineering University, 2016, 37(9): 1163-1171.

Coupled method for predicting motions of Spar-type offshore floating wind turbine systems

SHANG Jinghong1,2, ZHAO Yuna1, ZHANG Liang1, HU Changhong3, DING Xiaomeng1

(1.Deepwater Engineering Research Center， Harbin Engineering University, Harbin 150001, China; 2.CNOOC China Ltd.Tianjin 300452, China；3. Research Institute for Applied Mechanics, Kyushu University, Fukuoka 816-8580, Japan)

In this study, we established non-linear time-domain coupled dynamic equations and aero/hydro-dynamical models to study the interaction between wind waves and floating offshore turbine systems (FOWTs). We used Fortran code, time-frequency domain transformation, and the Runge-Kutta iteration method to solve nonlinear equations. We analyzed the 5DOF motion response characteristics of Spar FOWTs, except yawing, under steady wind wave, regular wave, and random wind-wave conditions. The results show that constant wind affects the average and peak-to-peak values of surge, pitch, heave, sway, and roll at the rated wind speeds, while it does not affect the peak-to-peak values of sway at other wind speeds. In addition, regular waves have only marginal effect on the average values of surge, pitch, heave, sway, and roll, but enlarge their peak-to-peak values. Under random wind and wave conditions, we characterized the system's sway and pitch motions for low and high wave frequencies , depending on the wind speed. The results of this study provide

for the design and hydrodynamic analysis of offshore floating wind turbine systems.

Spar-type offshore floating wind turbine system; non-linear time-domain coupled method; characteristics of motion response; regular wave; random wind and wave condition

2015-10-28.

時間：2016-09-07.

國家自然科學基金項目(11572094)；國家教育部博士點基金項目(P012213003)；高等學校學科創(chuàng)新引智計劃“111工程”(B07019).

尚景宏(1978-), 男,博士研究生；

趙玉娜,E-mail: zhangliang@hrbeu.edu.cn.

10.11990/jheu.201510073

TN911.7

A

1006-7043(2016)09-1163-09

張亮(1959-), 男,教授，博士生導師；

趙玉娜(1991-),女,博士研究生.

網(wǎng)絡出版地址：http://www.cnki.net/kcms/detail/23.1390.u.20160918.1547.004.html