王淵博， 李 春,2， 繆維跑

(1.上海理工大學 能源與動力工程學院，上海 200093；2.上海理工大學 能源與動力工程學院，上海市動力工程多相流動與傳熱重點實驗室，上海 200093)

基于致動線方法的風力機組風場布置優(yōu)化

王淵博1， 李 春1,2， 繆維跑1

(1.上海理工大學 能源與動力工程學院，上海 200093；2.上海理工大學 能源與動力工程學院，上海市動力工程多相流動與傳熱重點實驗室，上海 200093)

基于致動線方法，利用OpenFOAM對5種不同間距串列布置的具有3臺5 MW風力機的風場進行了模擬.通過調(diào)整風力機組間的位置，對比不同風場布置時的總輸出功率，并結合輪轂高度處流場分布及尾跡速度輪廓線，分析尾跡對風場下游風力機的影響.結果表明：串列布置時，處于尾跡中的風力機功率下降可達82.72%；改變風場中風力機組相對位置可使風場布置得到優(yōu)化，雖個別風力機性能有所下降，但風場總輸出功率最大可提高11.8%.

風場布置； 總輸出功率； 致動線； 風力機； OpenFOAM

風能在世界各國的能源供給中發(fā)揮著十分重要的作用[1-3].風場上游風力機產(chǎn)生的尾跡會導致下游速度虧損及湍流強度增加[4-5]，不僅大幅降低了下游風力機的輸出功率，浪費風場風和土地資源，更為重要的是還會導致葉片疲勞載荷顯著提高，降低風力機壽命周期，甚至導致葉片損壞事故[6].因此，優(yōu)化風場布置以減小風力機組間尾跡的影響顯得尤為重要[7-8].

目前，獲取風力機氣動性能的主要方法有葉素動量理論(BEM)方法、渦尾跡法和求解N-S方程的計算流體力學(CFD)方法.在上述方法中，BEM方法盡管計算效率高，但不能計算風力機尾跡流場，尚無法衡量風場上游風力機尾跡對下游風力機氣動性能的影響[9].與BEM方法相比，渦尾跡法有相對高的精度，且可以預測風力機尾跡，但隨風力機尺寸的增大，風力機繞流場渦結構跨越多個尺度，渦尾跡法對大尺度空間內(nèi)小渦發(fā)展的預測能力較低[10].因此，風力機氣動性能研究逐漸轉向更加精確的CFD方法[11]，即采用相關數(shù)值模擬方法求解N-S方程以獲取流場變量.鐵庚等[12]利用Fluent軟件對FFA-W3-211、FFA-W3-301和NA-CA63-215 3種翼型進行數(shù)值模擬，并與實驗數(shù)據(jù)進行對比，驗證了二維黏性不可壓數(shù)值模擬的可靠性.楊賀[13]利用Fluent軟件，基于10 kW水平軸風力機進行了葉片外形設計和風輪氣動性能數(shù)值模擬，最終得到葉片各葉素參數(shù)和風輪轉矩，結果驗證了設計的合理性.張昇龍[14]利用Fluent軟件對某商用1.5 MW水平軸風力機在3種不同風速條件下的風輪尾跡和氣動性能進行數(shù)值模擬，將計算結果與相似實驗測量值進行比較，證明了三維黏性不可壓模型的正確性.賈彥等[15]利用WAsP軟件及風資源數(shù)據(jù)進行風場模擬計算，建立了由多臺風力機組成的風場，比較不同風力機組布局下風場的年凈發(fā)電量，并給出了風場最優(yōu)風力機組布局.上述研究中，鐵庚等[12-14]僅對風力機葉片或不超過2 MW的單臺風力機進行CFD數(shù)值模擬，這是因為CFD方法需建立大量近壁面網(wǎng)格以求解邊界層流動，模擬多臺兆瓦級風力機流場需耗費巨大的計算資源，計算成本難以承受.賈彥等[15]雖對多臺風力機風場進行模擬，但WAsP軟件基于線性流動控制方程及半經(jīng)驗的尾跡模型，且假設流動不發(fā)生分離，計算結果與真實風場存在一定誤差，同時無法獲得風力機周圍流場信息[16].

Sorensen等[17]提出致動線方法，結合了BEM方法計算速度快的優(yōu)點與CFD方法計算精度高的優(yōu)點，通過較少的計算資源可獲得很好的流場信息，尤其適用于風力機組三維流動模擬.筆者基于致動線方法，通過開源CFD軟件OpenFOAM模擬了5種不同間距串列布置的3臺5 MW風力機繞流流場，研究風電機組間距對風電場整體性能的影響，通過分析風力機流場特征和尾跡速度分布，揭示上游尾跡干擾下游風力機吸收風能的流動機理，為未來大型風電場優(yōu)化布局提供參考.

1 致動線方法

1.1 致動線方法控制方程

致動線方法控制方程在笛卡爾坐標系下的張量形式為：

(1)

式中：下標i、j使用了愛因斯坦求和約定；t為時間；xi、xj為笛卡爾坐標系的坐標軸；vi、vj為速度沿坐標軸方向的分量；ρ為密度；p為壓強；μ為動力黏度；源項fi為體積力，表征葉片對流體的反作用力，為致動線方法控制方程特有，傳統(tǒng)CFD方法控制方程忽略此項.

1.2 體積力

1.2.1 體積力計算

根據(jù)牛頓第三定律，可利用流場對葉片的作用力求得葉片對流場的體積力.由于致動線方法中無實體壁面，需通過其他模型獲取葉片的作用力以達到與實體葉片相同的效果.以BEM方法作為附加模型，在每個時間步讀取流場信息，根據(jù)翼型氣動特性計算葉片的作用力，計算公式如下：

(2)

(3)

(4)

式中：W為當?shù)匾硇偷南鄬λ俣?；c為當?shù)匾硇拖议L；Cτ為切向力系數(shù)；τ為切向單位向量；Cn為法向力系數(shù)；n為法向單位向量；Fs為當?shù)匾碚箚挝婚L度所受到的作用力.

當?shù)匾碚箚挝婚L度葉片對流場的體積力為:

(5)

1.2.2 體積力投射

葉片對流場的體積力需投射至致動線附近的計算域網(wǎng)格，通過致動線方法控制方程式(1)求解下一時間步的流場.體積力投射規(guī)律對模擬結果的影響極大，當投射范圍較小時，數(shù)值計算將產(chǎn)生嚴重振蕩，甚至無法收斂；當投射范圍過大時，投射區(qū)域網(wǎng)格所屬體積力過小，流場狀態(tài)改變微弱，尤其是葉尖渦和葉根渦被過分抹平，與真實情況不相符[18].

為解決上述問題，將體積力以三維高斯分布向流場投射[17]，表達式如下：

(6)

式中：r為體積力中心至投射點的距離；ε為分布因子，控制體積力分布的集中程度.

Troldborg[19]對致動線方法中體積力的分布形式進行深入研究后指出，若采用高斯分布向流場投射體積力，分布因子ε=2Δs時的效果與真實流場最接近，其中Δs為致動線分布點間距.此外，為提高計算效率，同時兼顧精確度，需確定投射邊界，邊界上投射的體積力大小為最大體積力的5%，邊界之外不再進行體積力投射[20].

2 構造致動線求解器

PisoFoam是OpenFOAM自帶的標準求解器，可以求解湍流及層流流動，基于瞬態(tài)Piso算法編寫，相比Simple算法，該算法計算效率更高[21].圖1為瞬態(tài)Piso算法流程圖.

圖1 瞬態(tài)Piso算法流程圖

基于PisoFoam求解器構造新求解器PisoFoamTurbine，專用于基于致動線方法對風力機的模擬.與標準求解器PisoFoam相比，新求解器PisoFoamTurbine建立了“turbines”類，該類包括定義體積力數(shù)組、求解體積力大小和投射體積力等一系列操作.計算過程中，通過該類求解相應范圍內(nèi)體積力的大小，改變動量方程，并在式(1)源項中添加體積力.

PisoFoamTurbine求解器中動量方程變?yōu)橐韵滦问剑?/p>

fvVectorMatrix UEqn

(

fvm::ddt(U)

+ fvm::div(phi, U)

+ turbulence->divDevReff(U)

- turbines.force()

);

其中,turbines.force()函數(shù)是用來實現(xiàn)流場中體積力計算及完成相應空間分布的操作，其余計算均與PisoFoam一致.

3 計算模型與計算域網(wǎng)格劃分

3.1 計算模型

計算模型選用美國可再生能源實驗室(NREL)設計的5 MW風力機，主要參數(shù)如表1所示.

表1 NREL 5 MW風力機主要參數(shù)

3.2 風力機組布置

以D表示風輪直徑大小，在27D×6D×6D的模擬風場中串列布置3臺NREL 5 MW風力機，分別為WT1、WT2和WT3，如圖2所示.WT1距入口3D；WT3與WT1的距離為固定值14D；WT2與WT1的距離L為變量，取值如表2所示.由于尾跡長度小于5D時湍流強度急劇增強[22-23]，風力機疲勞載荷增大、壽命縮短，故L最小取5D以保證WT2不在強尾跡區(qū)，最大取9D以保證WT3不在強尾跡區(qū).

圖2 模擬風場中風力機組的布置

算例編號12345L5D6D7D8D9D

3.3 數(shù)值計算模型

Sarmast等[24]詳細研究了分布因子及網(wǎng)格尺度對致動線方法計算結果的影響：當分布因子確定時，優(yōu)化網(wǎng)格可加速計算結果收斂并提高準確性.因此，對模擬風場中風輪所在水平圓柱區(qū)域的網(wǎng)格進行加密，并在網(wǎng)格加密區(qū)與粗網(wǎng)格區(qū)之間添加2層網(wǎng)格過渡帶,以避免數(shù)值計算振蕩，如圖3所示.網(wǎng)格尺度從內(nèi)向外逐層加倍，其中加密區(qū)網(wǎng)格尺度為3.94 m×3.94 m×3.94 m，與Fleming等[25]模擬風場的網(wǎng)格尺度相近.

(a) 主視圖

(b) 三維透視圖

模擬風場的邊界條件類型如圖4所示.進口設為速度入口，方向垂直于進口平面，大小為11.4 m/s；出口設為壓力出口，大小為標準大氣壓；地面和頂面設為滑移邊界條件；兩個側面設為周期邊界條件，即計算域一個側面網(wǎng)格節(jié)點與另一個側面網(wǎng)格節(jié)點一一對應，緊鄰周期邊界一側計算域外的“鏡像單元”信息由緊鄰另一側的計算域內(nèi)的單元提供.

圖4 模擬風場的邊界條件類型

數(shù)值計算采用大渦模擬(LES)方法，流動變量通過網(wǎng)格過濾操作分為網(wǎng)格尺度(GS)和亞格子尺度(SGS)2部分.引入各向異性的殘余應力張量模型封閉亞格子模型，利用Smagorinsky渦粘假設參數(shù)化，將模型分為2部分：線性渦粘模型和基于混合長假設的亞尺度渦粘模型[26].

4 結果與分析

通過致動線方法在每個時間步求出風力機功率并實時輸出.當風從計算域入口發(fā)展到出口之后再計算200 s，風力機的輸出功率即取此時間段內(nèi)的平均值.圖5為不同L處各風力機輸出功率的數(shù)值模擬結果.

圖5 風力機的輸出功率

從圖5可以看出，當L分別為5D、6D、7D、8D和9D時，風場中3臺風力機總輸出功率分別為7.719 1 MW、7.477 4 MW、7.321 6 MW、7.014 8 MW和6.808 2 MW，即隨L的增大總輸出功率逐漸減小.通過最小二乘法擬合總輸出功率與L的關系，可得到式(7)：

P=-0.228 5l+8.867 3

(7)

式中：P為風場總輸出功率，MW；l為L中的數(shù)值部分.

隨著L的變化，WT2和WT3的輸出功率變化明顯.理論上，下游風場中的風力機要達到與上游風力機相同的輸出功率，則需與上游風場風力機保持足夠大的距離[27].WT1的輸出功率隨著L的增大稍微有所增大，證明下游風力機對WT1也有輕微影響，且隨著L增大此影響逐漸減弱.WT1輸出功率變化不大的原因是其處于風場上游且入流風速不變.

WT2和WT3的輸出功率比見表3.輸出功率比定義為WT2和WT3的輸出功率分別與WT1輸出功率的比值.風力機尾跡的湍流強度具有非線性特征[28]，所以WT2和WT3的輸出功率比與L呈現(xiàn)非線性關系.由表3可見，WT2的輸出功率沒有超過WT1輸出功率的21%，此現(xiàn)象在順排風場中較為常見[29-30].當下游風力機位于上游風力機尾跡區(qū)域時，尾跡的速度虧損和湍流強度增加造成下游風力機入流風速減小且極不穩(wěn)定，導致下游風力機輸出功率嚴重下降.因為WT3所處的尾跡由WT2和WT1的尾跡疊加而成，這種疊加尾跡的不穩(wěn)定性比非疊加尾跡高，可促進其與周圍流體的混合和能量交換，加速尾跡恢復，使尾跡速度虧損得以快速補償，湍流強度大幅減小，為下游風力機提供質(zhì)量相對較高的入流風，所以WT3的輸出功率整體上比WT2的高，這也證明了所研究內(nèi)容的準確性及可靠性.

表3 WT2和WT3的輸出功率比

通過最小二乘法分別擬合WT2和WT3輸出功率比與L的關系，可得到如式(8)和式(9)所示的關系式.

R2=0.008 7l+0.130 1

(8)

R3=-0.055 4l+0.662 2

(9)

式中：R2為WT2的輸出功率比；R3為WT3的輸出功率比.

借助速度云圖可更直觀地看出L對風力機組的影響，圖6為不同L處風力機輪轂高度處的速度云圖.從圖6可以看出，隨著L的增大，WT2的低速尾跡區(qū)域逐漸增大.當L=5D和L=6D時，WT2與WT3間距較大，加之WT1尾跡的疊加，WT2尾跡與周圍流體發(fā)生了較為充分的混合和能量交換，因此接近WT3處的尾跡速度得到一定提升，WT3入流速度較大，其輸出功率較WT2的高.當L=7D時，WT2的低速尾跡發(fā)展至WT3風輪平面，導致WT3入流速度較低，是WT3輸出功率降低的主要原因.隨著L的進一步增大，當L=9D時，WT2與WT3間距較小，WT2尾跡尚未與周圍流體發(fā)生足夠的混合和能量交換，尾跡速度虧損嚴重，極大地降低了WT3的輸出功率；與此同時WT2與WT1間距較大，WT2輸出功率受WT1尾跡影響較小，WT2輸出功率已高于WT3輸出功率.上述分析進一步說明了文中計算的可行性與實用性.

圖6 輪轂高度處速度云圖

此外，從圖6可以看出，WT3尾跡比WT2尾跡和WT1尾跡更為復雜，非穩(wěn)定流動明顯,因此與周圍流體的混合和能量交換速度更充分，使得尾跡得到快速恢復.

從圖5可知，風場總功率的變化主要由WT3引起，因此有必要進一步分析WT3的尾跡.圖7為不同L處WT3各風輪直徑輪轂高度處的速度輪廓線.從圖7可以看出，L取值不同時，WT3各風輪直徑輪轂高度處的速度輪廓線變化較大.L=9D時，WT3與WT2間距最小，WT2的尾跡還未充分恢復，導致WT3的入流風速非常低，此時速度輪廓線凹的幅度最大，WT3輸出功率相應最小.L=5D時，速度輪廓線相對比較平緩，說明WT3所處尾跡速度虧損較小.此時,WT3距WT2較遠，WT2的尾跡有足夠的時間和空間與周圍流體進行混合和能量交換，尾跡恢復充分.L分別為7D、8D和9D時，速度輪廓線逐步向上抬升，WT3的入流風速和輸出功率相應逐漸增大.由此可以看出所得結論的統(tǒng)一性和正確性.此外，從圖7還可以看出，尾跡附近區(qū)域風速略高于入流風速，這是風場上游風力機造成的風加速現(xiàn)象.

圖7 速度輪廓線

5 結 論

(1) 風力機尾跡處會出現(xiàn)速度虧損和湍流強度增加的現(xiàn)象，處于尾跡中的風力機，其入流風速相應減小且極不穩(wěn)定，輸出功率大幅降低.

(2) 疊加尾跡的非穩(wěn)定流動更加明顯，與周圍流體混合得更加充分，恢復更加完全.

(3) 適當調(diào)整風場中風力機的相對位置，可以優(yōu)化風場布置，從而大幅提高風場總輸出功率.文中模擬風場WT2與WT1的距離為L=5D時,雖然WT2的輸出功率較WT1減小82.72%，但風場總輸出功率提高了11.8%.

[1] 李少華, 匡青峰, 吳殿文, 等. 1.2 MW風力機整機流場的數(shù)值模擬[J]. 動力工程學報, 2011,31(7): 551-556.

LI Shaohua, KUANG Qingfeng, WU Dianwen, et al. Numerical simulation on flow field of a 1.2 MW wind turbine[J]. Journal of Chinese Society of Power Engineering, 2011, 31(7): 551-557.

[2] SAHU B K, HILOIDHARI M, BARUAH D C. Global trend in wind power with special focus on the top five wind power producing countries[J]. Renewable and Sustainable Energy Reviews, 2013, 19(1): 348-359.

[3] CHENG Shikun, LI Zifu, MANG H P, et al. A review of prefabricated biogas digesters in China[J]. Renewable and Sustainable Energy Reviews, 2013, 28: 738-748.

[4] 韓中合, 李引, 季劍. 考慮風剪切的1.3 MW風力機整機三維定常流動數(shù)值研究[J]. 動力工程學報, 2011, 31(10): 779-783.

HAN Zhonghe, LI Yin, JI Jian. Numerical study on 3D steady flow of a 1.3 MW wind turbine considering wind shear factor[J]. Journal of Chinese Society of Power Engineering, 2011, 31(10): 779-783.

[5] 黎作武, 賀德馨. 風能工程中流體力學問題的研究現(xiàn)狀與進展[J]. 力學進展, 2013, 43(5): 472-525.

LI Zuowu, HE Dexin. Reviews of fluid dynamics researches in wind energy engineering[J]. Advances in Mechanics, 2013, 43(5): 472-525.

[6] 王強. 水平軸風力機三維空氣動力學計算模型研究[D]. 北京: 中國科學院大學, 2014.

[7] OZTURK U A, NORMAN B A. Heuristic methods for wind energy conversion system positioning[J]. Electric Power Systems Research, 2004, 70(3): 179-185.

[8] INDASI V S, LYNCH M, MCGANN B, et al. WAsP model performance verification using lidar data[J]. International Journal of Energy and Environmental Engineering, 2016, 7(1): 105-113.

[9] SHENKAR R. Design and optimization of planar and nonplanar wind turbine blades using vortex methods[D]. Denmark: Technical University of Denmark, 2010.

[10] 夏巍巍. 水平軸風力機尾流結構及其對風電場布局影響的研究[D]. 蘭州: 蘭州理工大學, 2014.

[11] MIAO Weipao, LI Chun, JUN Yang, et al. Numerical investigation of wake control strategies for maximizing the power generation of wind farm[J]. Journal of Solar Energy Engineering, 2016, 138(3): 034501.

[12] 鐵庚, 祁文軍. 水平軸風力機葉片翼型的氣動特性數(shù)值模擬[J]. 流體機械, 2013, 41(3): 29-33.

TIE Geng, QI Wenjun. Numerical simulation on the aerodynamic performance of aerofoil of horizontal wind turbine blade[J]. Fluid Machinery, 2013, 41(3): 29-33.

[13] 楊賀. 水平軸風力機風輪氣動性能研究[D]. 呼和浩特: 內(nèi)蒙古工業(yè)大學, 2014.

[14] 張昇龍. 水平軸風力機尾流場及其氣動性能的研究[D]. 蘭州: 蘭州理工大學, 2014.

[15] 賈彥, 劉璇, 李華, 等. 考慮尾流效應對風電場機組布局的影響分析[J]. 可再生能源, 2014, 32(4): 429-435.

JIA Yan, LIU Xuan, LI Hua, et al. Analysis of wind farm units layout considering wake effect[J]. Renewable Energy Resources, 2014, 32(4): 429-435.

[16] 田琳琳, 趙寧, 武從海, 等. 復雜地形風電場的機組布局優(yōu)化[J]. 南京航空航天大學學報, 2013, 45(4): 503-509.

TIAN Linlin, ZHAO Ning, WU Conghai, et al. Optimal positioning of wind turbines on complex wind farm[J]. Journal of Nanjing University of Aeronautics & Astronautics, 2013, 45(4): 503-509.

[17] SORENSEN J N, SHEN Wenzhong. Numerical modeling of wind turbine wakes[J]. Journal of Fluids Engineering, 2002, 124(2): 393-399.

[18] 劉磊. 風力機葉片非定常氣動特性的研究[D]. 北京: 中國科學院大學, 2012.

[19] TROLDBORG N. Actuator line modeling of wind turbine wakes[D]. Denmark: Technical University of Denmark, 2008.

[20] JIN Wenjie. Numerical simulation of wind turbine wakes based on actuator line method in NEK5000[D]. Stockholm, Sweden: Royal Institute of Technology, 2013.

[21] 李東岳. 淺談Piso，Simple算法[EB/OL]. (2016-03-12) [2016-05-10]. http://www.dyfluid.Com/PISOSIMPLE.html.

[22] MO J O, CHOUDHRY A, ARJOMANDI M, et al. Large eddy simulation of the wind turbine wake characteristics in the numerical wind tunnel model[J]. Journal of Wind Engineering and Industrial Aerodynamics, 2013, 112: 11-24.

[23] WILSON J M, DAVIS C J, VENAYAGAMOORTHY S K, et al. Comparisons of horizontal-axis wind turbine wake interaction models[J]. Journal of Solar Energy Engineering, 2015, 137(3): 031001.

[24] SARMAST S, SEGALINI A, MIKKELSEN R F, et al. Comparison of the near-wake between actuator-line simulations and a simplified vortex model of a horizontal-axis wind turbine[J]. Wind Energy, 2016, 19(3): 471-481.

[25] FLEMING P, GEBRAAD P M O, LEE S, et al. Simulation comparison of wake mitigation control strategies for a two-turbine case[J]. Wind Energy, 2014, 18(2): 2135-3143.

[26] 卞鳳嬌, 徐宇, 王強, 等. 基于OpenFOAM的風力機致動線模型研究[J]. 工程熱物理學報, 2016, 37(1): 72-75.

BIAN Fengjiao, XU Yu, WANG Qiang, et al. Numerical study of actuator line model of wind turbine based on OpenFOAM platform[J]. Journal of Engineering Thermophysics, 2016, 37(1): 72-75.

[27] CHOWDHURY S, ZHANG Jie, MESSAC A, et al. Unrestricted wind farm layout optimization (UWFLO): investigating key factors influencing the maximum power generation[J]. Renewable Energy, 2012, 38(1): 16-30.

[28] SHAKOOR R, HASSAN M Y, RAHEEM A, et al. Wake effect modeling: a review of wind farm layout optimization using Jensen's model[J]. Renewable and Sustainable Energy Reviews, 2016, 58: 1048-1059.

[29] WU Yuting, PORTé-AGEL F. Simulation of turbulent flow inside and above wind farms: model validation and layout effects[J]. Boundary-Layer Meteorology, 2012, 146(2): 181-205.

[30] NEWMAN J, LEBRON J, MENEVEAU C, et al. Streamwise development of the wind turbine boundary layer over a model wind turbine array[J]. Physics of Fluids, 2013, 25(8): 231-234.

Layout Optimization of the Wind Farm in a Wind Power Station Based on Actuator Line Method

WANGYuanbo1,LIChun1,2,MIAOWeipao1

(1. School of Energy and Power Engineering, University of Shanghai for Science and Technology, Shanghai 200093, China; 2. Shanghai Key Laboratory of Multiphase Flow and Heat Transfer in Power Engineering, School of Energy and Power Engineering, University of Shanghai for Science and Technology, Shanghai 200093, China)

Based on actuator line method, numerical simulations were conducted on the wind farm with three 5 MW wind turbines in five different aligned arrangements using OpenFOAM. By comparing the total power output, the velocity fields and velocity profiles at the hub height of wind turbines for different wind farm layouts, the effects of upstream turbine wake on the downstream wind turbine were investigated. Results show that in an aligned arrangement, the power output may be reduced by 82.72% for the wind turbine lying in the wake. Changing the relative position of wind turbines reasonably could help to optimize the wind farm layout, when the total power output may be increased by a maximum of 11.8%, even if the output of an individual wind turbine declines somewhat.

wind farm layout; total power output; actuator line; wind turbine; OpenFOAM

2016-05-26

2016-08-04

國家自然科學基金資助項目(51176129，51676131)；上海市科學技術委員會資助項目(13DZ2260900)

王淵博(1991-)，男，陜西澄城人，碩士研究生，研究方向為風力發(fā)電. 李 春(通信作者)，男，教授，博導，電話(Tel.)：15216702797；E-mail：lichun_usst@163.com．

1674-7607(2017)05-0418-07

TK83

A 學科分類號：480.60