李迺璐,穆安樂,Balas M J
(1.揚州大學(xué)水利與能源動力工程學(xué)院,江蘇揚州225127;2.西安理工大學(xué)機(jī)械與精密儀器工程學(xué)院,西安710048 3.美國安柏瑞德航空航天大學(xué)航空學(xué)院,美國32114)
基于B-L氣動模型的旋轉(zhuǎn)水平風(fēng)機(jī)葉片經(jīng)典顫振穩(wěn)定性分析
李迺璐1,穆安樂2,Balas M J3
(1.揚州大學(xué)水利與能源動力工程學(xué)院,江蘇揚州225127;2.西安理工大學(xué)機(jī)械與精密儀器工程學(xué)院,西安710048 3.美國安柏瑞德航空航天大學(xué)航空學(xué)院,美國32114)
研究旋轉(zhuǎn)水平軸風(fēng)力機(jī)葉片周期時變氣彈系統(tǒng)的經(jīng)典顫振穩(wěn)定性特性。葉片結(jié)構(gòu)采用具有揮舞和扭轉(zhuǎn)耦合的典型界面振動模型,引入Beddoes-Leishman氣動模型為旋轉(zhuǎn)葉片提供低攻角處周期時變的非定常氣動力;為了研究顫振邊界,利用標(biāo)量風(fēng)速和揮舞/扭轉(zhuǎn)固有頻率比對所建立的旋轉(zhuǎn)葉片氣彈模型進(jìn)行變型。在此基礎(chǔ)上,通過時域響應(yīng)曲線分析旋轉(zhuǎn)葉片揮舞自由度和扭轉(zhuǎn)自由度氣彈穩(wěn)定性特性,分析了標(biāo)量速度,葉片剛度,揮舞/扭轉(zhuǎn)固有頻率比和結(jié)構(gòu)阻尼的影響,揭示了旋轉(zhuǎn)葉片經(jīng)典顫振邊界的變化規(guī)律且其準(zhǔn)確性得到驗證。
旋轉(zhuǎn)風(fēng)機(jī)葉片;經(jīng)典顫振;氣彈模型;穩(wěn)定性
近年來世界各國對綠色能源的需求導(dǎo)致大型風(fēng)機(jī)的裝機(jī)容量不斷上升。隨著機(jī)型的增大,風(fēng)力機(jī)葉片展向長度也隨之加大。目前兆瓦級的葉片尺寸都>30 m,因此在慣性力、彈性力和復(fù)雜氣動負(fù)載力耦合作用下,葉片會出現(xiàn)顫振現(xiàn)象。經(jīng)典顫振問題是研究大型風(fēng)力機(jī)葉片安全穩(wěn)定運行的一個重要問題。水平風(fēng)力機(jī)葉片的經(jīng)典顫振發(fā)生在葉片運行低攻角狀態(tài)的勢流中,為氣流流動基本附著無明顯分離情況下,風(fēng)機(jī)葉片扭轉(zhuǎn)自由度和揮舞自由度產(chǎn)生的自激振蕩[1]。由于顫振分析研究涉及葉片結(jié)構(gòu),葉片周圍復(fù)雜非定常氣動力以及葉片流固耦合等一系列復(fù)雜困難的問題,大部分經(jīng)典顫振研究都集中于葉片定常氣動力的顫振分析研究[2-4]。目前引入非定常氣動力的葉片經(jīng)典顫振研究也局限于葉片氣彈系統(tǒng)在靜態(tài)點的顫振特性研究。任勇生等[5]針對非時變?nèi)~片氣彈系統(tǒng),采用特征法畫出根軌跡曲線來分析葉片顫振穩(wěn)定性,Kallesoe[6]也采用特征值和特征向量來分析靜態(tài)攻角下的葉片氣彈特性。而目前針對旋轉(zhuǎn)風(fēng)力機(jī)葉片的研究又大多集中于復(fù)雜葉片結(jié)構(gòu)方程式的求解和線性化[7-9],或是基于風(fēng)洞實驗的葉片氣動力特性。
因此本文基于標(biāo)量化的揮舞/扭轉(zhuǎn)振動界面,引入Beddoes-Leishman非定常氣動模型為旋轉(zhuǎn)葉片在低攻角處提供周期時變非定常勢流氣動力,研究旋轉(zhuǎn)風(fēng)機(jī)葉片動態(tài)時變氣彈系統(tǒng)的經(jīng)典顫振穩(wěn)定性特性。該研究分析了標(biāo)量風(fēng)速,揮舞扭轉(zhuǎn)固有頻率比,葉片剛度和結(jié)構(gòu)阻尼對旋轉(zhuǎn)葉片氣彈穩(wěn)定性的影響,揭示了旋轉(zhuǎn)動態(tài)葉片經(jīng)典顫振邊界的變化規(guī)律。
1.1 Beddoes-Leishman氣動模型
Beddoes-Lesihman(B-L)模型不僅可以描述葉片的氣動失速特性,還可以描述葉片在勢流中的氣動特性。在本文中此模型用來計算旋轉(zhuǎn)葉片在勢流中的周期時變非定常氣動力。B-L氣動力模型方程如下[10]:
式中:x1,x2,x3,x4為四個氣動模型狀態(tài)量,前兩個方程式(x1,x2)為動態(tài)勢流特性,后兩個方程式(x3,x4)為動態(tài)氣流分離特性。αE為有效攻角,fst為靜態(tài)氣流分離點。α0為靜態(tài)升力系數(shù)為零時的攻角,b1,b2,A1,A2分別為時間遲延和幅值常量。Tu,Tp,Tf為時間常量,由此為基礎(chǔ)非定常氣動力系數(shù)可表示為:
1.2 周期時變勢流氣動力模擬
在對NACA0012葉片進(jìn)行風(fēng)洞實驗得到靜態(tài)氣動系數(shù)實驗數(shù)據(jù)的基礎(chǔ)上,定義勢流中低攻角為幅值變化[2o,8o]的正弦曲線,利用Beddoes-Lesihman模型計算出旋轉(zhuǎn)葉片周期時變的非定常氣動力。本文利用Matlab/Simulink進(jìn)行仿真和模擬(見圖1)。圖1(a)顯示了動態(tài)氣動升力的時域響應(yīng),可以發(fā)現(xiàn)隨著攻角的周期性變化,氣動負(fù)載升力也為周期性變化曲線;圖1 (b)顯示了動態(tài)氣動升力與攻角的關(guān)系。圖1中虛線為靜態(tài)氣動升力風(fēng)洞實驗數(shù)據(jù),實線為動態(tài)非定常氣動升力,攻角周期性變化范圍在線性區(qū)域內(nèi),表示葉片旋轉(zhuǎn)持續(xù)工作在勢流中,此時完全沒有分離流的影響。計算出的動態(tài)升力系數(shù)偏離靜態(tài)值較小,基本呈線性變化,符合勢流中的旋轉(zhuǎn)葉片氣動特性并可用于經(jīng)典顫振研究。
圖1 旋轉(zhuǎn)葉片非定常氣動升力Fig.1 Unsteady aerodynamic lift of the rotating blade
圖1所示旋轉(zhuǎn)葉片的振動運動包括揮舞和扭轉(zhuǎn),分別由揮舞彎曲位移h和扭轉(zhuǎn)偏移角度θ來表示。b為葉片截面的半弦長,重心與彈性軸的距離為xθb,葉片的振動受制于彎曲力Qh和扭轉(zhuǎn)力Qθ(見圖2)。
圖2 葉片截面結(jié)構(gòu)模型Fig.2 The structuralmodel of the blade section
兩自由度揮舞-扭轉(zhuǎn)耦合的結(jié)構(gòu)模型運行方程可表示為[12]:
式中:m為葉片截面的質(zhì)量,Iθ為關(guān)于彈性軸的質(zhì)量矩,Ch,Cθ為結(jié)構(gòu)阻尼系數(shù),Kh,Kθ為揮舞,扭轉(zhuǎn)彈簧常量。彎曲力Qh和扭轉(zhuǎn)力Qθ表示為:
式中:ρ為空氣密度,V為相對風(fēng)速,Cl,θ,Cm,θ為旋轉(zhuǎn)葉片的氣動升力和氣動力矩系數(shù)。定義無綱量揮舞彎曲位移為h=h/b,標(biāo)量風(fēng)速V*=V/(bωθ),Ch=2ζhωhm,Cθ=2ζθωθIθ,Kh=ωh2m,Kθ=ωθ2Iθ,其中ζh,ζθ為揮舞和扭轉(zhuǎn)的臨界阻尼系數(shù),ωh,ωθ為解耦的旋轉(zhuǎn)葉片揮舞和扭轉(zhuǎn)固有頻率,定義rθ2=Iθ/mb2,結(jié)構(gòu)式(3),式(4)可表示為:
式中:M為馬赫數(shù),μ=m/πρb2,ωh/ωθ為旋轉(zhuǎn)葉片揮舞/扭轉(zhuǎn)固有頻率比,Cl,θ,Cm,θ為B-L模型計算的旋轉(zhuǎn)葉片非定常氣動力。由于非定常氣動力為周期時變函數(shù),此旋轉(zhuǎn)葉片的氣彈方程為周期時變系統(tǒng)。
數(shù)值計算的主要參數(shù)基于旋轉(zhuǎn)葉片翼型NACA0012[12],m=12.387 kgm-1,b=0.135 m,xθ= 0.246 6,ρ=1.225 kg/m3,M=0.3。根據(jù)建立的旋轉(zhuǎn)葉片氣彈系統(tǒng)方程,本文研究標(biāo)量速度,結(jié)構(gòu)阻尼,葉片剛度,固有頻率比等對旋轉(zhuǎn)葉片經(jīng)典顫振穩(wěn)定性特性的影響,以及對顫振邊界的特性進(jìn)行分析。
3.1 顫振速度判定
首先不考慮結(jié)構(gòu)阻尼的情況下,研究在不同的無綱量風(fēng)速V*下葉片截面揮舞和扭轉(zhuǎn)自由度的氣彈穩(wěn)定性特性。給定參數(shù)ωh/ωθ=0.2,無綱量旋轉(zhuǎn)葉片扭轉(zhuǎn)固有頻率ωθ=1,ζh=ζθ=0。通過不斷增加V*值,直到揮舞響應(yīng)和扭轉(zhuǎn)響應(yīng)中至少有一個開始出現(xiàn)不穩(wěn)定震蕩即判斷對應(yīng)的V*值為顫振速度。圖3顯示了給定不同V*值的旋轉(zhuǎn)葉片截面氣彈系統(tǒng)的時域響應(yīng)。結(jié)果顯示V*在[1,7]區(qū)域內(nèi)旋轉(zhuǎn)葉片的扭轉(zhuǎn)自由度是處于臨界穩(wěn)定的,揮舞顫振發(fā)生在V*=7.1,此時扭轉(zhuǎn)響應(yīng)曲線開始出現(xiàn)不穩(wěn)定震蕩;而此時旋轉(zhuǎn)葉片的揮舞自由度并沒有出現(xiàn)明顯的不穩(wěn)定震蕩現(xiàn)象。這說明扭轉(zhuǎn)振動相對揮舞振動隨風(fēng)速的變化更敏感,此模型參數(shù)下的旋轉(zhuǎn)葉片經(jīng)典顫振速度為V*=7.1。當(dāng)V*=9大于顫振速度時,揮舞和扭轉(zhuǎn)響應(yīng)曲線的不穩(wěn)定震蕩更加嚴(yán)重。Bruining[13]中采用半經(jīng)驗?zāi)P蜑槿~片提供氣動負(fù)載,與此結(jié)構(gòu)模型參數(shù)類似,處于低攻角的靜態(tài)葉片在V*=4時保持穩(wěn)定,在V*=8時氣彈不穩(wěn)定。因此,判定的顫振速度V*值是符合葉片經(jīng)典顫振氣彈特性的,其準(zhǔn)確性得到驗證。
圖3 揮舞/扭轉(zhuǎn)顫振時域響應(yīng)曲線Fig.3 Time-domain responses of flapwise and torsional deflection
3.2 葉片剛度影響
考慮葉片剛度的影響,保持標(biāo)量風(fēng)速V*=9和揮舞扭轉(zhuǎn)固有頻率比ωh/ωθ=0.2不變,增加旋轉(zhuǎn)葉片揮舞固有頻率和扭轉(zhuǎn)固有頻率。圖4顯示了不同葉片剛度下的旋轉(zhuǎn)葉片顫振響應(yīng)曲線。圖中虛線表示參數(shù)ωθ=1,ωh=0.2(剛度較小)對應(yīng)的響應(yīng)曲線,實線表示參數(shù)ωθ=2,ωh=0.4(剛度較大)對應(yīng)的響應(yīng)曲線。由圖4可見,當(dāng)固有頻率增大時,揮舞自由度的不穩(wěn)定震蕩變?yōu)榉€(wěn)定震蕩,扭轉(zhuǎn)自由度的不穩(wěn)定震蕩幅值明顯減小,提高了旋轉(zhuǎn)葉片的穩(wěn)定性;隨著固有頻率的增加,揮舞震蕩幅值相較扭轉(zhuǎn)震蕩幅值減少的更明顯。由此可知,剛度較大的旋轉(zhuǎn)葉片穩(wěn)定性更好并且揮舞顫振可以得到更為有效的抑制。
圖4 不同葉片剛度下的顫振響應(yīng)曲線Fig.4 Vibration responses based on different blade stiffness
3.3 葉片阻尼影響
保持其他參數(shù)不變,不同葉片阻尼因素對應(yīng)的旋轉(zhuǎn)葉片揮舞/扭轉(zhuǎn)彎曲響應(yīng)曲線(見圖5)。隨著結(jié)構(gòu)阻尼的增加,顫振響應(yīng)可以從臨界穩(wěn)定達(dá)到漸進(jìn)穩(wěn)定,當(dāng)結(jié)構(gòu)阻尼進(jìn)一步增加時,顫振幅值也隨著減小,同時達(dá)到穩(wěn)定的時間也大大縮短。由此可見,結(jié)構(gòu)阻尼可以明顯改善旋轉(zhuǎn)葉片經(jīng)典顫振的氣彈穩(wěn)定性。
3.4 揮舞/扭轉(zhuǎn)固有頻率比影響
考慮葉片揮舞和扭轉(zhuǎn)固有頻率比的影響,保持其他參數(shù)不變,增加旋轉(zhuǎn)葉片的揮舞和扭轉(zhuǎn)固有頻率之比,通過不同固有頻率比下的系統(tǒng)時域響應(yīng)曲線分析來確定其顫振速度。圖6顯示了固有頻率比ωh/ωθ在[0.5,1.5]區(qū)域內(nèi)旋轉(zhuǎn)葉片的揮舞和扭轉(zhuǎn)響應(yīng)曲線。由圖可見,隨著ωh/ωθ的增加,揮舞自由度和扭轉(zhuǎn)自由度由臨界穩(wěn)定變?yōu)椴环€(wěn)定震蕩最后又回到穩(wěn)定震蕩,說明ωh/ωθ的變化對于揮舞自由度和扭轉(zhuǎn)自由度都有較大的影響。同時扭轉(zhuǎn)自由度在ωh/ωθωh/ωθ值為[1,1.5]區(qū)域內(nèi)時域響應(yīng)由發(fā)散震蕩變?yōu)樗p震蕩,說明在這個區(qū)間內(nèi)固有頻率比的變化對扭轉(zhuǎn)自由度有著非常大的影響。結(jié)果表明,揮舞/扭轉(zhuǎn)固有頻率比的變化對于揮舞/扭轉(zhuǎn)自由度穩(wěn)定性影響都比較大,同時扭轉(zhuǎn)自由度更為敏感。
圖5 不同葉片阻尼下的顫振響應(yīng)曲線Fig.5 Vibration responses based on different blade damping
圖6 不同固有頻率比下的顫振響應(yīng)Fig.6 Vibration responses based on different ratio between flapwise and torsional natural frequencies
3.5 顫振邊界特性
為了研究不同固有頻率比下旋轉(zhuǎn)葉片的穩(wěn)定性及旋轉(zhuǎn)葉片經(jīng)典顫振邊界的特性,利用時域曲線穩(wěn)定性分析ωh/ωθ=0.2,0.5,1.25,1.5,4.0時對應(yīng)的旋轉(zhuǎn)葉片的顫振臨界速度V,臨界葉片剛度ωθ和臨界結(jié)構(gòu)阻尼ζh,ζθ。顫振臨界速度的計算分析基于參數(shù)ωθ=1,ζh=ζθ=0;臨界剛度的計算分析基于參數(shù)V*=2,ζh= ζθ=0;臨界阻尼的計算分析基于參數(shù)ωθ=1,V*=8。由表1可知,在固有頻率比較低時,如[0.2,1.25]區(qū)域內(nèi),顫振速度隨著固有頻率比的增加而減小;而在固有頻率比較高時,如[1.5,4]區(qū)域內(nèi),顫振速度隨著固有頻率比的增加而增加。在ωh/ωθ值為1.5時,具有最小的顫振速度為V=2.2。因此,旋轉(zhuǎn)葉片經(jīng)典顫振邊界隨著固有頻率比的增加會先降低再升高,并且會在中間某個位置達(dá)到最小值。說明旋轉(zhuǎn)葉片在揮舞和扭轉(zhuǎn)固有頻率比值較低或較高時的氣彈穩(wěn)定性較好,隨著低固有頻率比的增加和高固有頻率比的減小,旋轉(zhuǎn)葉片的氣彈穩(wěn)定性都會降低。此結(jié)論與Theodorsen[14]穩(wěn)定邊界的變化趨勢一致,說明了所建模型的準(zhǔn)確性。同時在風(fēng)速較小并不考慮結(jié)構(gòu)阻尼影響的情況下分析葉片臨界剛度,當(dāng)固有頻率比<1(扭轉(zhuǎn)固有頻率較高)時,隨著扭轉(zhuǎn)固有頻率的減小,需要較高的剛度來保證葉片顫振扥穩(wěn)定性;當(dāng)固有頻率比>1(揮舞固有頻率較高)時,較高的揮舞固有頻率就可以保證穩(wěn)定性而不需要較高的葉片剛度。最后,考慮風(fēng)速較大時固定葉片剛度下的臨界結(jié)構(gòu)阻尼,阻尼邊界受固有頻率比的影響和剛度邊界類似,但是阻尼邊界值普遍較小,說明結(jié)構(gòu)阻尼有很好的顫振抑制作用,特別當(dāng)揮舞固有頻率明顯高于扭轉(zhuǎn)固有頻率時,結(jié)構(gòu)阻尼取任意值都是穩(wěn)定的。
表1 旋轉(zhuǎn)葉片顫振邊界Tab.1 Flutter speed of the rotating w ind turbine blade
本文基于標(biāo)量化的葉片截面結(jié)構(gòu)模型和B-L氣動模型,針對周期時變的旋轉(zhuǎn)水平軸風(fēng)力機(jī)葉片氣彈系統(tǒng)進(jìn)行了經(jīng)典顫振穩(wěn)定性分析研究。利用時域曲線分析顫振速度,其準(zhǔn)確性進(jìn)行了驗證。結(jié)果表明,剛度較大的旋轉(zhuǎn)葉片有較好的氣彈穩(wěn)定性;結(jié)構(gòu)阻尼可以有效的抑制旋轉(zhuǎn)葉片顫振的發(fā)生;揮舞扭轉(zhuǎn)固有頻率比對扭轉(zhuǎn)自由度和揮舞自由度的穩(wěn)定性影響都比較大。顫振邊界的研究表明,旋轉(zhuǎn)葉片在固有頻率比較低和較高的時候有良好的氣彈穩(wěn)定性,在中間值域內(nèi)氣彈穩(wěn)定性會變差;較高的扭轉(zhuǎn)固有頻率仍然需要一定的葉片剛度和結(jié)構(gòu)阻尼來保證旋轉(zhuǎn)葉片的顫振穩(wěn)定性,而具有較高揮舞固有頻率的葉片則對剛度和結(jié)構(gòu)阻尼要求較低。
[1]孫麗萍,王昊,丁嬌嬌.風(fēng)力發(fā)電機(jī)葉片的氣動彈性及顫振研究綜述[J].液壓與氣動,2012,(10):1-5.
SUN Li-ping,WANG Hao,DING Jiao-jiao.A summery on the aeroelasticity and flutter ofwind turbine blade[J].Journalof Chinese Hydraulics and Pneumatics,2012,(10):1-5.
[2]Reddy K K,Chen J,Behal A.Multi-input/multi-output adaptive output feedback control design for aeroelastic vibration suppression[J].Journal of Guidance,Control,and Dynamics,2007,30(4):1040-1048.
[3]Kukreja S.Non-linear system identification for aeroelastic systems with application to experimental data[M].NASA TM-214641,2008.
[4]McEver M A,Ardelean E V,Cole D G.Active control and closed-loop identification of flutter instability in typical section airfoil[J].Journal of Guidance,Control,and Dynamics,2007,30(3):733-740.
[5]任勇生,張明輝.水平軸風(fēng)力機(jī)葉片的彎扭耦合氣彈穩(wěn)定性研究[J].振動與沖擊,2010,29(7):196-200.
REN Yong-sheng,ZHANG Ming-h(huán)ui.Aeroelastic stability of a horizontal axis wind turbine blade with bending-torsion coupled[J].Journal of Vibration and Shock,2010,29(7): 196-200.
[6]Kallesoe B.S.,A low-order model for analysing effects of blade fatigue load control[J].Wind Energy,2006,(9): 421-436.
[7]Park JH,Park H Y,Jeong S Y,et al.Linear vibration analysis of rotating wind-turbine blade[J].Current Applied Physics,2010,10(2):332-334.
[8]莫文威,李德源,夏鴻建,等.水平軸風(fēng)力機(jī)柔性葉片多體動力學(xué)建模與動力特性分析[J].振動與沖擊,2013,32 (22):99-105.
MOWen-wei,LIDe-yuan,XIA Hong-jian,et al.Multibody dynamic modeling and dynamic characteristics analysis of flexible blades for a horizontal axiswind turbine[J].Journal of Vibration and Shock,2013,32(22):99-105.
[9]王丹,陳予恕,曹慶杰,等.兩自由度彎曲耦合渦輪機(jī)械葉片動力學(xué)分析[J].振動與沖擊,2014,33(7):199-204.
WANG Dan,CHEN Yu-shu,CAO Qing-jie,et al.Dynamic analysis of a 2-DOF turbomachine blade with coupling of bending and torsion[J].Journal of Vibration and Shock,2014,33(7):199-204.
[10]Hansen M H,Gaunaa M,Madsen H A.A beddoes-leishman type dynamic stallmodel in state-space and indicial formulations[R].Riso National Laboratory,Denmark,2004.
[11]Mahajan A J,Kaza K R V.Semi-empiricalmodel for prediction of unsteady forces on an airfoil with application to flutter[J].Journal of Fluids and Structures,1993(7):87-103.
[12]Galvanetto U,Peiro J.An assessment of some effects of the nonsmoothness of the leishman-beddoes dynamic stall model on the nonlinear dynamics of a typical aerofoil section[J].Journal of Fluids and Structures,2008(24):151-163.
[13]Bruining A,Timmer W A.Airfoil characteristics of rotating wind turbine blades[J].Journal of Wind Engineering and Industrial Aerodynamics,1992,39(3):35-39.
[14]Thoedorsen T,Garrick IE.Mechanism of flutter,a theoretical and experimental investigation of the flutter problem[R].NACA Report685,1940.
Classical flutter stability of rotating horizontal w ind turbine blades based on B-L aeroelasticmodel
LINai-lu1,MU An-le2,BalasM J3
(1.School of Hydraulic,Energy and Power Engineering,Yangzhou University,Yangzhou 225127,China;
2.School of Mechanical and Precision Instrument Engineering,Xi'an Technology University,Xi'an 200215,China; 3.Aerospace Engineering,Embry-Riddle Aeronautical University,F(xiàn)L 32114,USA)
Classical flutter stability of a rotating horizontal axiswind turbine blade as a periodic,time-varying and aeroelastic system was studied.The structuralmodel of the bladewas a classic vibration system with 2-DOF(flapwise and torsional blade oscillation),and periodic time-varying aerodynamic loads were offered with Beddoes-Leishman dynamic model at the low angle of attack.The normalized structural model was used to analyze the effects of normalized wind velocity,blade stiffness,ratio between flapwise and torsional natural frequencies and blade damping on the aeroelastic stability of the rotating blade.The flutter boundary was found to reveal the characteristics of classic flutter of the rotating blade.
rotating wind turbine blade;classic flutter;Beddoes-Leishman;stability
V211.47;TK8
A
10.13465/j.cnki.jvs.2015.23.030
江蘇省高校自然科學(xué)研究面上項目(14KJB480006);教育部留學(xué)回國人員科研啟動基金資助項目;揚州大學(xué)科技創(chuàng)新培育基金(2014CXJ028);國家自然科學(xué)基金(51075326);美國能源部項目(DESC0001261)
2014-07-25修改稿收到日期:2014-10-17
李迺璐女,博士,副教授,1985年生