柯世堂， 余 瑋， 王同光

(1.南京航空航天大學 土木工程系，南京 210016；2.南京航空航天大學 江蘇省風力機設計高技術研究重點實驗室，南京 210016)

基于大渦模擬考慮葉片停機位置大型風力機風振響應分析

柯世堂1,2， 余 瑋1,2， 王同光2

(1.南京航空航天大學 土木工程系，南京 210016；2.南京航空航天大學 江蘇省風力機設計高技術研究重點實驗室，南京 210016)

為研究停機狀態(tài)下不同葉片位置對大型風力機塔架-葉片體系風振響應的影響，以某3 MW大型水平軸三葉片風力機為研究對象?；诖鬁u模擬(LES)方法對葉片八個不同停機位置下的風力機體系流場和氣動力性能進行了數(shù)值模擬，并通過與國內(nèi)外實測數(shù)據(jù)對比驗證了該方法的有效性；結合有限元法對考慮葉片不同停機位置的風力機塔架-葉片耦合模型進行了動力特性和風振響應時域分析。主要結論為：不同葉片停機位置對風力機塔架繞流特性和氣動力分布影響顯著，塔架中上部風振響應受葉片不同停機位置的影響最大，尤其是迎風面和背風面的脈動位移均方差響應較大，相應最大值出現(xiàn)在工況八下塔頂350°位置；最大塔底彎矩出現(xiàn)在工況二的環(huán)向330°位置；工況五下三葉片葉尖順風向位移響應峰值>2.5 m；研究表明在進行大型風力機抗風設計時應考慮不同葉片停機位置的影響。

大型風力機體系；停機狀態(tài)；大渦模擬；風振響應；氣動力分布

大型風力機屬于典型的風敏感結構，其塔架-葉片耦合結構具有風場復雜、阻尼小、自振頻率低且分布密集等特性，風荷載往往成為其設計的主要控制荷載[1-2]。由于作用在風力機的外部載荷主要是由風況條件決定，高風速條件下風力機將處于停機狀態(tài)，葉片不同停機位置將引起風力機體系表面流場和氣動力分布的差異從而導致風振響應不同。因此，對考慮葉片不同停機位置下大型風力機體系風振響應特性進行研究可為設計人員進行高風速下風力機抗風設計提供參考。

目前，針對大型風力機的風致動力響應分析，王振宇等[3]基于Shinozuka理論，采用諧波疊加法分別模擬良態(tài)風和臺風時程，對良態(tài)風和臺風作用下的風力機進行動力時程分析，研究發(fā)現(xiàn)臺風風向突然偏轉90°將會使塔架發(fā)生屈服破壞；陳小波等[4]分別分析了屋頂集風效應和離心剛化效應對風力機塔架風振反應的影響，結果表明屋頂集風效應增強了塔架風振反應而離心剛化效應減小了塔架風振反應；KWON等[5]基于風力機塔架-葉片耦合模型提出了一種新的陣風因子計算方法，并通過算例分析表明陣風荷載效應對于大型風力機塔架體系的局部等效靜風荷載的分布的影響顯著；MURTAGH等[6]采用了時域強迫振動分析方法對風力機耦合體系進行了順風向響應分析，結果表明不考慮塔架葉片相互影響會導致塔頂?shù)娘L致響應減??；董輝等[7]基于FLUENT軟件對2 MW水平軸風力機在額定工況下遭受暴雨時的風場和風驅雨進行分析，研究了風雨共同作用下風力機風雨荷載的大小和分布規(guī)律，結果表明雨滴對塔架的沖擊作用不可忽略；金志昊等[8]對風力機葉片采用時域和頻域方法進行了順風向風致振動計算，分析了暴雨停機狀態(tài)下葉片極限位移最不利情況。已有研究大多是針對葉片固定工況下大型風力機氣動力數(shù)值模擬和風振響應的分析，鮮有系統(tǒng)考慮葉片不同停機位置下風力機塔架-葉片體系的風振響應特性研究。

鑒于此，本文以南京航空航天大學自主研發(fā)的某3 MW水平軸三葉片風力機為研究對象，首先基于大渦模(Large Eddy Simulation, LES)擬法對停機狀態(tài)下風力機體系8種計算工況(由葉片旋轉全過程狀態(tài)下和塔架的相對位置確定)進行了數(shù)值模擬，獲取了塔架和葉片表面典型測點的氣動力時程作為后續(xù)風振響應分析的輸入?yún)?shù)；然后結合有限元法對考慮葉片不同停機位置的風力機塔架-葉片耦合模型進行了動力特性和風振響應時域分析，最終對比研究了不同停機位置對于大型風力機體系流場特性、氣動力分布和風振響應的影響，相關結論對大型風力機抗風的精細化設計具有參考意義。

1 動力特性分析

以某3 MW風力發(fā)電機組為研究對象，其中塔架高85 m，頂部直徑2.0 m，底部直徑2.5 m，塔架為通長變厚度結構，頂壁厚30 mm，底壁厚60 mm；機艙長12.0 m，寬4.0 m，高4.0 m；風輪切入風速3.5 m/s，額定風速12.5 m/s，切出風速25 m/s，風輪傾角5°，風輪轉速范圍為9～19 r/min，偏航速度為0.5度/秒；各葉片之間成120°夾角，沿周向均勻分布，葉片長度為44.5 m。

根據(jù)葉片與塔架的相對位置，并考慮到三葉片體系旋轉過程中存在的周期性，以葉片與豎直方向夾角為0°作為初始狀態(tài)，依次順時針旋轉15°設置8個計算工況?；贏NSYS平臺分別建立不同工況下風力機塔架-葉片一體化有限元模型，其中塔架和葉片采用shell63單元，機艙采用beam188梁單元模擬，忽略其內(nèi)部結構，圓形筏基基礎的單元類型為solid65，尺寸為直徑24 m，高度2 m，基礎底端固結，地基與基礎的作用采用combin14彈簧單元模擬。風力機各部件通過多點約束單元耦合形成塔架-葉片一體化有限元模型[9-10]。依據(jù)效率與精度均衡的原則，整機模型共劃分3 473個單元。不同計算工況下有限元模型見圖1。

圖1 不同計算工況有限元模型示意圖

表1給出了采用分塊Lanczos方法計算提取的風力機典型階的固有頻率，可以看出葉片和塔架的耦合效應明顯降低了風力機系統(tǒng)的基頻，不同計算工況下風力機系統(tǒng)基頻均很低，僅為0.209 Hz左右，且各模態(tài)之間頻率間隔較小。葉片不同停機位置主要影響風力機體系的低階固有頻率，各工況低階頻率出現(xiàn)較小的差別，高階頻率一致?？梢娡C狀態(tài)下葉片不同位置對風力機體系的頻率和振型影響微弱，低階振型主要以葉片帶動機艙及塔架進行前后揮舞和左右擺動，高階模態(tài)出現(xiàn)塔架及葉片本身的結構變形和失穩(wěn)形態(tài)。

表1 不同計算工況下風力機典型階固有頻率列表

2 數(shù)值模擬

2.1 計算域及網(wǎng)格劃分

為保證流動能夠充分發(fā)展，計算流域取12D×5D×5D(流向x×展向y×豎向z，D為風輪直徑)。風力機置于距離計算域入口3D，從而保證尾流的充分流動。由于葉片表面扭曲復雜故采用混合網(wǎng)格離散形式，將整個計算域分為內(nèi)外兩個部分：核心區(qū)域采用四面體網(wǎng)格，并對風力機周圍局部網(wǎng)格進行加密，外圍區(qū)域采用高質量六面體結構網(wǎng)格，網(wǎng)格總數(shù)為795萬，計算域及網(wǎng)格具體劃分見圖2。

2.2 邊界條件及計算方法

計算域進口采用速度入口，按照B類地貌設置風速剖面以及湍流強度剖面，其中湍流強度定義式為

(1)

式中：地面粗糙指數(shù)為0.15，10 m參考高度處的基本風速為25 m/s。

對應該風力機的切出風速，并通過用戶自定義函數(shù)實現(xiàn)上述入流邊界條件與FLUENT的連接(見圖3)；計算域出口設置為壓力出口，相對壓力為0；計算域

地面以及風力機表面采用無滑移壁面，計算域兩側面和頂面采用對稱邊界條件。

數(shù)值計算采用3D單精度、分離式求解器，由于風力機所處流場屬于非定常且擾流情況復雜，基于大渦模擬法能夠對風力機復雜的流場進行更好地模擬[11]。亞格子尺度選用Smagorinsky-Lilly模型，壓力速度耦合方程組求解采用SIMPLEC格式，壓力項離散采用standard格式，動力離散采用BCD(Bounded Central Differencing)格式，瞬態(tài)方程采用二階隱式，控制方程的計算殘差設置為1×10-6，以模型高度和模型高度處的平均風速計算時間步長定為0.001 s。

2.3 結果分析及有效性驗證

圖4和圖5分別給出了不同工況下風力機體系典型截面的渦量分布以及速度流線圖。由圖4和圖5可知，葉片和塔架之間的相互干涉作用將引起風力機體系大范圍的渦量增值區(qū)域，且葉片不同停機位置引起的渦量分布差異明顯，葉片與塔架相對位置越接近時塔架背后形成的尾渦區(qū)域越明顯；同時葉片與塔架的相對位置不同將引起葉片引起的尾渦不同，塔架的存在破壞了葉片尾渦，而葉片尾渦也將改變塔架渦的形狀；來流由葉片流經(jīng)塔架時，塔架周圍形成了明顯的速度增大區(qū)，并且隨著葉片與塔架之間的位置靠近程度增加，塔架周圍流體加速越顯著。

(a) 網(wǎng)格區(qū)域劃分

(b) x-y平面網(wǎng)格分布

(c) x-z平面網(wǎng)格分布

圖3 平均風速剖面及湍流強度分布Fig.3 Vertical profiles of wind speed and turbulence intensity

風力機體系表面共設置126個測點，其中塔架沿子午向布置8層測點，每層環(huán)向12個測點，各葉片沿翼展方向布置10個測點[12-13]，基于以上大渦模擬法獲得風力機體系表面風壓系數(shù)時程曲線，以作為風荷載時程輸入?yún)?shù)。限于篇幅，圖6僅給出了塔架中上部以及葉片典型測點的壓力系數(shù)時程曲線，塔架迎風面呈正壓分布，塔架兩側及背風面呈現(xiàn)出負壓分布且兩側的負壓明顯大于塔架背風面。由于葉片的阻擋作用使得來流出現(xiàn)撞擊、分離、附著以及漩渦脫落等特征湍流的情況，同時塔架與葉片之間的相互干擾作用增大使得塔架表面壓力系數(shù)的均方差較大。

圖4 不同工況風力機塔架典型截面渦量分布圖Fig.4 Vorticity distribution of the wind turbine tower under different working conditions

圖5 不同工況風力機塔架典型截面速度流線圖Fig.5 Velocity contour distribution with streamlines of the wind turbine tower under different working conditions

圖6 風力機體系典型測點壓力系數(shù)時程曲線Fig.6 The history curve of pressure coefficient of wind turbine typical points

圖7給出了風力機塔架下部受葉片干擾影響較小的典型截面沿環(huán)向分布的平均壓力系數(shù)和脈動風壓系數(shù)，并與規(guī)范以及國內(nèi)外實測分布曲線[14-15]進行對比，由圖7可知，大渦模擬與規(guī)范的平均風壓系數(shù)沿環(huán)向分布規(guī)律基本一致，僅在背風區(qū)數(shù)值略小于規(guī)范值；脈動風壓分布曲線分布在國內(nèi)外實測曲線中間，且沿環(huán)向分布規(guī)律較為接近，考慮到脈動風壓分布與實測塔所處的地形、來流湍流和周邊干擾密切相關，本文大渦模擬獲得的脈動風壓分布趨勢和數(shù)值均在已有實測結果包絡中，對比結果也驗證了數(shù)值模擬的有效性，可用于后續(xù)的風振響應時域分析。

(a) 平均風壓

(b) 脈動風壓

3 風振響應分析

風力機塔架-葉片耦合模型風振響應時域計算基于ANSYS軟件平臺，采用大渦模擬得到的風壓系數(shù)作為風荷載時程輸入?yún)?shù)，采用隱式方法Newmark和HHT來求解瞬態(tài)動力學平衡方程，其中各模態(tài)阻尼均為0.02，積分時間步長取為0.5 s，加載時間步數(shù)為1 024 步。

3.1 塔架響應

圖8給出了不同工況下塔架徑向位移響應均值的分布圖，其中負值表示徑向朝內(nèi)，正值表示徑向朝外。由圖8可知，不同工況下塔架徑向位移均值分布規(guī)律較一致，葉片不同位置對塔架位移主要影響塔架中上部。徑向位移隨著塔架高度的增加而逐漸增大，最大正負位移均出現(xiàn)在塔頂處350°和170°位置，以工況六位移最大為0.193 m，工況二位移最小為0.133 m；不同計算工況下塔架位移極值在迎風面和背風面的分布范圍不同，極值區(qū)域隨著葉片與塔架位置相對干涉作用增大而增大，工況六極值區(qū)域最大在塔架50 m以上，工況二極值區(qū)域最小在塔架70 m以上。

圖9給出了不同計算工況下塔架位移脈動響應的均方差分布圖，由圖9可知，① 塔架位移脈動響應的均方差隨著塔架高度增加而逐漸增加；② 塔架中上部葉片與塔架之間的相互干涉作用顯著，由于葉片和塔架背風面產(chǎn)生渦旋的脫落，同一高度環(huán)向分別在塔架迎風面和背風面范圍出現(xiàn)兩個峰值區(qū)域；③ 脈動響應均方差極值區(qū)域最小出現(xiàn)在工況三(>80 m)，最大在工況六和工況八(>55 m)；④ 位移響應均方差最大均出現(xiàn)在塔頂350°位置，工況八數(shù)值最大為0.292，工況三最小為0.134。

圖8 不同計算工況塔架徑向位移響應均值Fig.8 The radial displacement of mean response of tower under different calculate working conditions

為對比不同葉片位置引起塔架底部彎矩值的差異，以工況一作為初始狀態(tài)，將各工況下塔底內(nèi)力響應值與之作差，圖10給出了不同計算工況下塔底內(nèi)力響應的特征值。由圖10可知：① 塔底徑向彎矩均值在0°、60°、90°以及330°位置出現(xiàn)差異，其中工況五下較明顯，工況二0°位置數(shù)值相差最大達到256.93 N·m；② 塔底環(huán)向彎矩均值以工況二和工況五的差異較明顯，數(shù)值相差最大位于工況二330°；③ 塔底彎矩脈動響應均方差以工況六和工況八下0°和330°位置脈動增加更顯著。

圖9 不同計算工況塔架徑向響應均方差Fig.9 RMS of radial displacement response of tower under different calculate working conditions

圖10 不同計算工況塔架底部內(nèi)力響應特征值Fig.10 Characteristic force responses of tower bottom under different calculate working conditions

3.2 葉片響應

表2給出了不同計算工況下葉根彎矩與剪力的特征值，由表可以發(fā)現(xiàn)葉片不同停機位置引起了葉根內(nèi)力特征值的顯著差異，剪力與彎矩均值均隨著塔架與葉片相對位置的靠近呈現(xiàn)減小的趨勢，以工況一內(nèi)力均值最大，工況五內(nèi)力均值最小；工況二和工況八內(nèi)力脈動響應的均方差較大。

定義工況一中與豎直方向夾角為0°的葉片編號為1，順時針旋轉編號葉片2、葉片3。表3給出了不同計算工況下風力機各葉片葉尖順風向位移響應的峰值，其中峰值因子取為2.5。

表2 不同計算工況葉根內(nèi)力特征值列表

由表3可知：① 工況五各葉片極大負值均>-2.50 m，工況三各葉片的極值位移較小；② 其余工況下均出現(xiàn)單個葉片的極大負值>-2.00 m，其中最大負值出現(xiàn)在工況八葉片3葉尖處；③ 葉片2不同停機位置下會對塔架產(chǎn)生遮擋效應，以工況五下葉片遮擋效應最顯著，極大負值位移達到-2.530 m，此時葉片與塔架之間相互干擾作用最大。因此，葉片不同位置導致葉尖位移出現(xiàn)明顯差異，為保證停機狀態(tài)下葉片與塔架之間的安全性，應考慮最不利工況下葉尖的位移。

表3 不同計算工況葉尖順風向位移峰值列表

4 結 論

基于大渦模擬和有限元法，對葉片不同停機位置下風力機塔架-葉片體系的流場作用、氣動力分布和風振響應特性進行了系統(tǒng)研究，主要結論如下：

(1)停機狀態(tài)下葉片不同位置對風力機體系的頻率和振型影響微弱，低階振型主要以葉片帶動機艙及塔架進行前后揮舞和左右擺動，高階模態(tài)出現(xiàn)塔架及葉片本身的結構變形和失穩(wěn)形態(tài)。

(2)葉片和塔架之間的相互干涉作用將引起風力機體系大范圍的渦量增值區(qū)域，葉片與塔架相對位置越接近時塔架背后形成的尾渦區(qū)域越明顯。

(3)葉片不同停機位置顯著影響塔架中上部的位移響應，其中工況六的位移響應均值最大，各工況下位移最大值均出現(xiàn)在塔頂170°和350°處，脈動響應均方差在迎風面和背風面出現(xiàn)峰值區(qū)域，最大值出現(xiàn)在工況八塔頂350°處。

(4)塔底彎矩響應特征值在環(huán)向0°和330°處差異顯著，其中均值以工況二和工況五最明顯，脈動響應均方差以工況六和工況八差異最顯著。

(5)葉根剪力與彎矩均值隨著塔架與葉片相對位置的靠近呈現(xiàn)減小的趨勢，工況五中三葉片葉尖位移峰值均>-2.50 m，葉尖最大位移峰值出現(xiàn)在工況八。

[ 1 ] 柯世堂，王同光. 偏航狀態(tài)下風力機塔架-葉片耦合結構氣彈響應分析[J]. 振動與沖擊, 2015, 34(18): 33-38. KE Shitang， WANG Tongguang. Aero-elastic responses analysis of wind turbine tower-blade coupled model in yaw condition[J]. Journal of Vibration and Shock, 2015, 34(18): 33-38.

[ 2 ] KE S T, GE Y J, WANG T G, et al. Wind field simulation and wind-induced responses of large wind turbine tower-blade coupled structure[J]. The Structural Design of Tall and Special Buildings, 2015, 24(8): 571-590.

[ 3 ] 王振宇, 張彪, 趙艷, 等. 臺風作用下風力機塔架振動響應研究[J]. 太陽能學報, 2013, 34(8): 1434-1442. WANG Zhenyu, ZHANG Biao, ZHAO Yan， et al. Dynamic response of wind tuebine under typhoon[J]. Acta Energiae Solaris Sinica， 2013, 34(8): 1434-1442.

[ 4 ] 陳小波, 李靜, 陳健云. 考慮葉片離心剛化效應的風力機塔架風振反應分析[J]. 工程力學, 2010, 27(1): 240-245. CHEN Xiaobo, LI Jing, CHEN Jianyu. Wind-include response analysis of wind turbine tower considering centrifugal stiffening effect of blade[J]. Engineering Mechanics， 2010, 27(1): 240-245.

[ 5 ] KWON D K, KAREEM A, BUTLER K. Gust-front loading effects on wind turbine tower systems[J]. Journal of Wind Engineering & Industrial Aerodynamics, 2012, 104(3):109-115.

[ 6 ] MURTAGH P J, BASU B, BRODERICK B M. Along-wind response of a wind turbine tower with blade coupling subjected to rotationally sampled wind loading[J]. Engineering Structures, 2005, 27(8):1209-1219.

[ 7 ] 董輝, 高乾豐, 鄧宗偉, 等. 大型風力機風雨荷載特性數(shù)值研究[J]. 振動與沖擊, 2015, 34(15): 17-22. DONG Hui, GAO Qianfeng, DENG Zongwei, et al. Numerical simulations for wind and rain loads of large-scale wind turbines[J]. Journal of Vibration and Shock, 2015, 34(15): 17-22.

[ 8 ] 金志昊, 范宣華, 蘇先樾,等. 風力機葉片順風向風致振動研究[J]. 南京航空航天大學學報, 2011, 43(5): 677-681. JIN Zhihao, FAN Xuanhua, SU Xianyue, et al. Research on blade wind-included vibration in wind direction[J]. Journal of Nanjing University of Aeronautics & Astronautics， 2011, 43(5): 677-681.

[ 9 ] 柯世堂, 曹九發(fā), 王瓏, 等. 風力機塔架-葉片耦合模型風致響應時域分析[J]. 湖南大學學報(自然科學版), 2014, 41(4): 87-93.KE Shitang, CAO Jiufa, WANG Long， et al. Time-domain analysis of the wind-included responses of the coupled model of wind turbine tower-blade coupled system[J]. Journal of Hunan University(Natural Science)， 2014, 41(4): 87-93.

[10] 柯世堂, 王同光, 陳少林, 等. 大型風力機全機風振響應和等效靜力風荷載[J]. 浙江大學學報(工學版), 2014, 48(4): 686-692. KE Shitang, WANG Tongguang,CHEN Shaolin, et al. Wind-included responses and equivalent static wind load of large wind turbine system[J]. Journal of Zhejiang University(Engineering Science), 2014, 48(4): 686-692.

[11] LUO Kun, ZHANG Sanxia , GAO Zhiying, et al. Large-eddy simulation and wind-tunnel measurement of aerodynamics and aeroacoustics of a horizontal-axis wind turbine[J]. Renewable Energy, 2014, 77: 351-362.

[12] 彭文春, 鄧宗偉, 高乾豐, 等. 風機塔筒流固耦合分析與受力監(jiān)測研究[J]. 工程力學, 2015, 32(7): 136-142. PENG Wenchun, DENG Zongwei, GAO Qianfeng, et al. Fluid-solid interaction analysis and stess monitoring research of wind turbine tower[J]. Engineering Mechanics, 2015, 32(7): 136-142.

[13] 黃帥, 宋波, 賀文山, 等. 考慮脈動風影響的風電塔風致動力響應數(shù)值分析與現(xiàn)場監(jiān)測比較[J]. 土木工程學報, 2012(增刊1): 102-106. HUANG Shuai, SONG Bo, HE Wenshan, et al. Numerical analysis of dynamic response and in-site monitoring of wind power tower considering the influence of the fluctuating wind load[J]. China Civil Engineering Journal, 2012(Sup1): 102-106.

[14] 李曉娜,陸煜,劉慶寬,等. 圓形截面煙囪風致干擾效應試驗研究[J]. 工程力學, 2015(增刊1): 159-162. LI Xiaona, LU Yu, LIU Qingkuan, et al. Experimental study on wind-included interference effects of circular section chimneys[J]. Engineering Mechanics, 2015(Sup1): 159-162.

[15] NISHIMURA H, TANIIKE Y. Aerodynamic characteristics of fluctuating forces on a circular cylinder[J]. Journal of Wind Engineering & Industrial Aerodynamics, 2001, 89(1): 713-723.

Wind-included vibration response analysis for a large wind turbine blade-tower systembased on large eddy simulation

KE Shitang1，2, YU Wei1，2, WANG Tongguang2

(1.Department of Civil Engineering, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China；2.Key Laboratory of Jiangsu Province High Tech Design of Wind Turbine, Nanjing University ofAeronautics and Astronautics, Nanjing 210016, China)

In order to study the effects of different blade positions on wind-included vibration responses of a large wind turbine tower-blade system under a stopped status, a 3 MW wind turbine with 3 blades developed by Nanjing University of Aeronautics and Astronautics was studied. First of all, the effects of eight different blade positions on the flow field and aerodynamic performances of the wind turbine system under stopped status were simulated with the large eddy simulation(LES). In order to verify the effectiveness of the LES method, the simulation results were compared with the actual test data of home and abroad. Then the finite element method was used to analyge dynamic characteristics and wind indaced vibration responses of the wind turbine tower-blade coupled model with different blade positions under a stopped status. The results showed that different blade positions have obvious effects on the flow characteristis of the wind turbine tower frame and the aerodynamic force distribution;the effects of different blade positions on the wind-induced vibration responses of the upper-middle part of the tower are the maximum;the mean square deviation between fluctuating displacements of windward side and leeward side is larger,its maximum value occurs at the tower top 350°position under the operation condition of 8;the maximum bending moment of the towerbottom appears at the annular 330°postion under the operation condition of 2;the peak displacements of 3 blade tips along the wind reach more than 2.5 m under the operation condition of 5.The study indicated that the effects of different blade positions under a stopped status should be considered for the anti-wind design of large wind turbine systems.

large wind turbine system；stopped status；large eddy simulation；wind-included vibration response；aerodynamic force distribution

國家973計劃項目(2014CB046200)；中國博士后科學基金(2015T80551)

2015-10-28 修改稿收到日期：2016-02-25

柯世堂 男，博士，副教授，1982年生

TK83

A

10.13465/j.cnki.jvs.2017.07.014