楊夢姚，毛璐璐，韓兆龍,b,c，周 岱,b,c，雷 航，曹 宇

(上海交通大學(xué) a. 船舶海洋與建筑工程學(xué)院；b. 海洋工程國家重點實驗室；c. 水動力學(xué)教育部重點實驗室， 上海 200240)

風(fēng)能是一種重要的能源.2018年底，我國風(fēng)電累計裝機容量2.1億千瓦，占全球風(fēng)電市場的三分之一以上[1].風(fēng)力發(fā)電機根據(jù)轉(zhuǎn)動軸與來風(fēng)方向的相對位置分為水平軸風(fēng)力機和垂直軸風(fēng)力機兩類.以H型風(fēng)機為代表的垂直軸風(fēng)力機具有產(chǎn)生噪聲更小、整體結(jié)構(gòu)更簡單、生命周期更長、發(fā)電效率更高等特點，具有廣闊的發(fā)展前景[2].H型垂直軸風(fēng)力機通常由2～6片與轉(zhuǎn)軸平行的葉片組成，其結(jié)構(gòu)特征使得必有葉片會處于迎風(fēng)位置，易產(chǎn)生離心力，主軸螺栓連接處易斷裂，葉片在風(fēng)荷載的作用下易產(chǎn)生振動問題，抗風(fēng)能力較弱[3].

圖1 垂直軸風(fēng)力機風(fēng)振與減振技術(shù)路線圖Fig.1 Flowchart of wind vibration and vibration reduction of VAWT

近年來，研究者對垂直軸風(fēng)機(VAWT)風(fēng)振響應(yīng)開展了一些研究[4-7].潘博等[8]提出研究垂直軸風(fēng)機必須考慮風(fēng)致振動及氣動力與葉片的共振.Campos等[9]提出結(jié)構(gòu)柔度問題在風(fēng)機生命周期評定中十分重要.對結(jié)構(gòu)耗能減振，設(shè)置阻尼器是較為常用的方法.劉保文等[10]對大跨度斜拉橋橋塔設(shè)置縱向阻尼器，有效減小了主塔在地震荷載作用下的內(nèi)力，使得內(nèi)力重新分布，改善了結(jié)構(gòu)受力性能.周云等[11]在某高層建筑上布置若干黏滯阻尼器，減震效果良好.目前，關(guān)于垂直軸風(fēng)力機減振的研究較少.本文針對H型垂直軸風(fēng)力機振動問題，在風(fēng)機結(jié)構(gòu)的不同位置上分別布置阻尼器，采用數(shù)值模擬的方法分析其耗能減振能力.

基于計算流體動力學(xué)，采用STAR-CCM+軟件建立垂直軸風(fēng)力機的空氣動力學(xué)模型，計算風(fēng)機旋轉(zhuǎn)時葉片所受氣動荷載.將氣動荷載在MATLAB軟件中進行擬合，獲得時序性風(fēng)機葉片荷載函數(shù)，并將其施加到基于ANSYS軟件建立的垂直軸風(fēng)力機的結(jié)構(gòu)模型上，從而進行動力學(xué)分析，得到風(fēng)機結(jié)構(gòu)的位移響應(yīng)時程.為抑制風(fēng)力機的振動，在風(fēng)機不同位置處施加基于COMBIN40單元而模擬的阻尼器，計算阻尼器耗能減振的能力.

1 風(fēng)力機力學(xué)模型

本文的技術(shù)路線如圖1所示.H型垂直軸風(fēng)力機的幾何模型如圖2所示(其中數(shù)字為節(jié)點編號)，其通常由2～6片葉片組成.目前2葉片或3葉片風(fēng)機在風(fēng)電市場上應(yīng)用較廣[1]，茲選取典型3葉片H型垂直軸風(fēng)力發(fā)電機為對象，研究風(fēng)機旋轉(zhuǎn)時風(fēng)荷載對風(fēng)機振動穩(wěn)定性的影響.垂直軸風(fēng)機多采用對稱翼型，故本文中風(fēng)機的葉片翼型截面選擇NACA-0021[12]，D型主梁處材料采用碳纖維，葉片蒙皮處材料采用玻璃纖維.主軸與支桿材料采用Q345鋼材[13]，彈性模量E=200 GPa，泊松比μ=0.17，密度ρ=7.8×103kg/m3.為減小計算量，選用某大型10 MW風(fēng)機的100∶1縮尺模型.風(fēng)機模型的幾何尺寸為：弦長C=0.265 m，支撐桿長D=0.8 m，葉片高度H=1.2 m，主軸和撐桿截面皆取圓截面，主軸直徑為0.2 m，撐桿直徑為0.1 m.

圖2 H型垂直軸風(fēng)力發(fā)電機的幾何模型Fig.2 Geometric model of a H-rotor type vertical axis wind turbine

2 分析方法

2.1 風(fēng)力機風(fēng)荷載模擬

由于風(fēng)機翼型截面的復(fù)雜性，故較難手動計算出風(fēng)機表面荷載的數(shù)值.計算流體動力學(xué)方法常應(yīng)用于解決流固耦合問題[14]，茲基于計算流體動力學(xué)理論開展H型垂直軸風(fēng)力機葉片表面風(fēng)荷載時程的數(shù)值模擬研究.采用雷諾平均(Reynolds Averaged Navier-Stokes, RANS)和標(biāo)準(zhǔn)k-ε湍流模型[15]作為數(shù)值模擬方法，運用STAR-CCM+軟件進行數(shù)值計算.雷諾平均納維-斯托克斯(N-S)方程為

(1)

式中：ui、uj分別為i、j方向的速度分量，i、j為坐標(biāo)軸方向；t為時間；xi、xj為坐標(biāo)；ρ為流體密度；p為壓力；ν為動力黏度系數(shù)；上標(biāo)“ ′ ”代表來流的脈動值.

由于本文的研究目的側(cè)重于風(fēng)機減振，所以對于風(fēng)機的氣動荷載采用一些簡化計算，具體如下：① 風(fēng)力機處于自然大氣邊界層中，因其葉片高度不大，故可不考慮風(fēng)速梯度的影響.② 考慮到在實際工程中，近地面的湍流強度較大.且當(dāng)風(fēng)力機轉(zhuǎn)速恒定時，來流風(fēng)速越低，風(fēng)機葉尖速越大，因風(fēng)輪轉(zhuǎn)動引起的湍流強度隨之增大.在這種高葉尖速比的條件下，湍流強度的大小會對風(fēng)力機的氣動力造成一定影響[16].為簡化模型，本文的計算條件參考相關(guān)風(fēng)洞實驗條件設(shè)置[16],其中來流風(fēng)速為8 m/s，湍流度為0.005.③ 因本文的研究對象是風(fēng)力機的縮尺模型(1∶100)，故本文的雷諾數(shù)是原模型的1/100，根據(jù)文獻[16]，雷諾數(shù)越大，風(fēng)機的空氣動力學(xué)性能越好，同時，由于葉片升力系數(shù)的增大，轉(zhuǎn)子轉(zhuǎn)矩也更大.因此本文所采用的縮尺模型的氣動載荷特性與原模型存在一定的偏差，模型的折算氣動載荷略小于實際載荷.④ 在CFD數(shù)值模擬過程中，為簡化計算，未考慮葉片結(jié)構(gòu)變形的影響.

圖3 風(fēng)場計算域與邊界條件Fig.3 Calculation domain and boundary conditions of wind field

為減少邊界條件的影響，在選取風(fēng)場計算域時，取其尺寸為31D×12D×5H[17]即24.8 m×9.6 m×6 m.將垂直軸風(fēng)機的中心置于計算域的前1/3處，以確保計算域風(fēng)場流動發(fā)展完全，如圖3所示.入口選用速度邊界條件，計算風(fēng)速取8 m/s；出口選用自由擴展邊界條件；風(fēng)輪區(qū)域流體定義為旋轉(zhuǎn)流體，風(fēng)輪定義為相對旋轉(zhuǎn)固體壁面；壁面選用滑移壁面.對風(fēng)場計算域進行網(wǎng)格劃分，使近葉片周圍的網(wǎng)格較密集，遠離葉片的網(wǎng)格較稀疏.由此共劃得棱柱形網(wǎng)格約1.89×106個，風(fēng)場計算域的阻塞率小于3%，無量綱化的壁面距離y+<1，滿足數(shù)值模擬對于網(wǎng)格的要求[18].

圖4 葉片氣動荷載時程曲線與函數(shù)擬合結(jié)果對比Fig.4 Comparison of blade aerodynamic load time history curve and function fitting results

本文計算得到葉片的轉(zhuǎn)矩荷載時程曲線結(jié)果如圖4所示.圖中：T0為轉(zhuǎn)矩，F(xiàn)T為切向力，F(xiàn)N為法向力.將其與Lei等[19]的轉(zhuǎn)矩時程曲線進行對比，發(fā)現(xiàn)具有相同的波形且荷載數(shù)值接近.

由于風(fēng)機的氣動荷載呈現(xiàn)周期性，所以可采用周期函數(shù)對其進行擬合.在MATLAB中，采用多階正弦函數(shù)對風(fēng)機的切向力、法向力及轉(zhuǎn)矩時程曲線進行擬合，獲得相應(yīng)的荷載函數(shù),便于在ANSYS中對結(jié)構(gòu)進行瞬態(tài)分析.以葉片的法向力為例，任一時間序列的荷載函數(shù)可近似模擬為

(2)

式中：fn(t)為t時間下的葉片法向力；ai、bi及ci分別為第i個正弦函數(shù)的各項系數(shù).

如圖4所示，法向力、切向力及轉(zhuǎn)矩的荷載函數(shù)曲線的擬合效果較好，幾乎與氣動荷載時程曲線保持一致.因此，采用多階正弦函數(shù)對氣動荷載數(shù)據(jù)進行擬合的方法是可行的.

2.2 結(jié)構(gòu)風(fēng)致振動響應(yīng)分析

運用ANSYS/Mechanical APDL 18.2軟件開展有限元數(shù)值模擬研究，采用Newmark-β方法進行瞬態(tài)分析.Newmark-β方法是一種隱式時間積分方案，具有二階時間精度[20]，此處選取γ=0.5，β=0.25(γ和β分別為控制數(shù)值分析精度和穩(wěn)定性的參數(shù))，即平均常加速度法.邊界約束條件為底端固接，在葉片表面施加氣動荷載，用四面體型三維實體單元對模型進行網(wǎng)格劃分，共得到5.234×104個網(wǎng)格單元.圖5為葉片表面所施加的荷載示意圖.

圖5 葉片施加荷載示意圖Fig.5 Schematic diagram of blade loading

風(fēng)機結(jié)構(gòu)在荷載激勵下的運動方程可以表示為

(3)

用有限元方法求解上述方程，可以得到各節(jié)點的風(fēng)致響應(yīng).

2.3 阻尼器模擬

在建筑物或空間結(jié)構(gòu)上安裝單個或多個阻尼器，可以抑制結(jié)構(gòu)的振動響應(yīng)[21-22].阻尼器會隨著風(fēng)機的振動而振動，其產(chǎn)生的慣性力作用可部分抵消風(fēng)荷載對風(fēng)機所產(chǎn)生的力作用，從而減小振動水平.為驗證阻尼器在垂直軸風(fēng)力機上的減振性能，在原結(jié)構(gòu)(模型1)的主軸與支桿連接處，即節(jié)點2和節(jié)點3處分別增加1個阻尼器，得到模型2和模型3，如圖6(a)、(b)所示.

圖6 加阻尼器后風(fēng)力機模型示意圖與彈簧-阻尼單元Fig.6 VAWT model with dampers and spring-damping unit

本文參考文獻[23]，采用多維減振阻尼器，它是6個自由度的并聯(lián)機構(gòu)，由法蘭盤、連接耳板及阻尼元件3部分組成.其阻尼原件主要是由6個黏彈性阻尼單元構(gòu)成的，在ANSYS軟件中，將其簡化為3個一維軸向COMBIN40 單元，提供三向平動自由度.3個一維轉(zhuǎn)動COMBIN40單元提供三向轉(zhuǎn)動自由度[23].COMBIN40單元是由相互平行的彈簧滑動器和阻尼器組成，并且串聯(lián)著一個間隙控制器的彈簧單元.在本文的計算中，將其簡化為相互平行的彈簧和阻尼器單元，如圖6(c)所示，其中k為彈簧剛度，c為阻尼系數(shù).

圖7 懸臂梁受動力荷載的時程響應(yīng)Fig.7 Time-history response of cantilever beam at dynamic load

阻尼元件的剛度和阻尼系數(shù)對于減振效果的影響非常明顯.具體體現(xiàn)為：當(dāng)阻尼元件剛度不變時，峰值位移在整體范圍內(nèi)隨著阻尼系數(shù)的增大而減小，當(dāng)剛度足夠大時，阻尼系數(shù)的改變對峰值位移的影響不大.當(dāng)阻尼系數(shù)不變時，峰值位移在整體范圍內(nèi)隨著剛度的增大而減小，當(dāng)阻尼系數(shù)足夠大時，剛度的改變對峰值位移的影響不大.基于上述規(guī)律并經(jīng)過反復(fù)試算，最終假定阻尼元件剛度為5×107N/m，阻尼系數(shù)為2.5×107N/(m·s).

3 數(shù)值算例驗證

3.1 動力作用下懸臂梁變形驗證

在結(jié)構(gòu)力學(xué)中，二維懸臂梁問題既是經(jīng)典問題也是基準(zhǔn)問題，被廣泛用于驗證或測試研究者提出的各種數(shù)值方法.為驗證本文方法的可靠性，考慮了在簡諧荷載作用下的懸臂梁變形問題.問題定義及荷載曲線如圖7(a)、(b)所示.懸臂梁的長度L=10 m，梁截面寬度b=1 m，高度h=1 m，彈性模量E=2.3×1010Pa，泊松比μ=0.27.在懸臂梁自由端施加簡諧荷載D(t)，荷載函數(shù)為D(t)=1 000sint，如圖7(b)所示.計算所得最大位移為0.177 mm，如圖7(c)所示，圖中Uy為懸臂梁末端y向位移時程，與文獻結(jié)果[24]對比，響應(yīng)幅值與響應(yīng)頻率與文獻結(jié)果相符.

3.2 時間步長選取

對于運用ANSYS分析瞬態(tài)動力學(xué)問題，選取不同的時間步長和時間步會對結(jié)果產(chǎn)生影響.時間步長越小，計算精度越高，但相應(yīng)會造成計算資源的浪費，因此在數(shù)值模擬前應(yīng)當(dāng)先選取合適的時間步長.茲對計算模型進行在3種時間步長下的模擬，即：0.004、0.002及0.001 s，分別得到對應(yīng)的結(jié)構(gòu)關(guān)于x軸的底部彎矩M0(見表1)：

(4)

式中:hi為結(jié)構(gòu)物的高度；P(x)為x高度處結(jié)構(gòu)表面壓力；A，Ai為受力面積.

表1 不同時間步長下底部彎矩比較

從表1看出，當(dāng)時間步長為0.002 s時，Mx為107.47 N·m，與時間步長取0.001 s時所得的結(jié)果相差僅0.7%，可看作近似相等.故時間步長為0.002 s時，既能保證所得結(jié)果的精確性，也能節(jié)省計算資源.茲選取時間步長為0.002 s.

4 數(shù)值結(jié)果與分析

4.1 風(fēng)荷載模擬

本文僅考慮來流風(fēng)速8 m/s，風(fēng)機轉(zhuǎn)子轉(zhuǎn)速17.92 rad/s[19]的工況下，垂直軸風(fēng)機所承受的氣動荷載.風(fēng)機旋轉(zhuǎn)1個周期所需的時間 0.350 6 s.圖8所示為第6個周期計算結(jié)束即t=2.103 6 s時其中1個葉片的壓力(p′)分布，圖中坐標(biāo)系為葉片的局部坐標(biāo)系.結(jié)果表明：較大的正壓力出現(xiàn)在葉片翼型前緣，與文獻結(jié)果相符[25].圖4所示為葉片上的瞬時風(fēng)壓隨時間的變化趨勢.可見，在前兩個周期內(nèi)，相鄰周期的壓力峰值差距較大，在風(fēng)機旋轉(zhuǎn)了2～3個周期之后，壓力數(shù)值趨于穩(wěn)定.

圖8 風(fēng)機葉片壓力云圖Fig.8 Pressure cloud diagram of fan blade

4.2 風(fēng)機結(jié)構(gòu)的風(fēng)致響應(yīng)分析

在運用ANSYS軟件計算垂直軸風(fēng)力機的振動響應(yīng)時，由于垂直軸風(fēng)力機旋轉(zhuǎn)過程中葉片的氣動荷載隨時間呈周期性變化，所以，將周期性變化的氣動荷載加載到風(fēng)機結(jié)構(gòu)上，即是把原來的風(fēng)載恒定，風(fēng)機轉(zhuǎn)動的問題轉(zhuǎn)化成了風(fēng)機不動，而風(fēng)載隨時間變化，從而計算風(fēng)機的振動響應(yīng).在風(fēng)機結(jié)構(gòu)響應(yīng)分析中，為簡化計算，未考慮結(jié)構(gòu)整體轉(zhuǎn)動的自由度.同時，本文中未考慮風(fēng)機轉(zhuǎn)動對風(fēng)振響應(yīng)的影響，風(fēng)機轉(zhuǎn)動可能會使風(fēng)機的風(fēng)振響應(yīng)更為復(fù)雜，在后續(xù)研究中將進一步展開分析.

分析發(fā)現(xiàn)t=10 s時，結(jié)構(gòu)的振動幅值較大，故選取t=10 s時刻下的總位移云圖分析，如圖9所示,
圖中Usum為風(fēng)力機結(jié)構(gòu)的總位移.圖9(a)所示為t=10 s時，原結(jié)構(gòu)(模型1)的總位移云圖.可見，模型1的最大位移發(fā)生在位于葉片頂端的節(jié)點16處，達到2.7 mm.即在該時刻下，風(fēng)機的危險點為節(jié)點16，由于風(fēng)機葉片的對稱性，可重點關(guān)注葉片頂端節(jié)點.圖9(b)所示為t=10 s時，在上部支撐桿處加了阻尼器后的模型2總位移圖，此時最大位移發(fā)生在位于葉片底端的節(jié)點9處，數(shù)值為1.5 mm，下降達44.4%.圖9(c)所示為t=10 s時，在下部支撐桿處加了阻尼器后模型3的總位移圖，此時最大位移仍舊發(fā)生在位于葉片頂端，數(shù)值為1.6 mm，下降達40.7%.

為更清晰地對比結(jié)構(gòu)的三向位移，用圖10表示t=10 s時3個模型在x、y及z向的位移云圖，其中，Ux、Uy及Uz分別對應(yīng)模型在x、y及z向的位移.由圖10對比x向位移數(shù)值，可見模型1的Ux分布在-2.7～1.0 mm，模型2的Ux分布在-0.5～0.8 mm，模型3的Ux分布在-0.9～0.9 mm，位移最大值全部發(fā)生在葉片頂端或底端.對比y向位移數(shù)值，可見模型1的Uy分布在-1.5～2.0 mm，模型2的Uy分布在-0.4～1.2 mm，模型3的Uy分布在-0.5～1.4 mm.對比z向位移數(shù)值，可見模型1的Uz分布在-0.13～0.25 mm，模型2的Uz分布在-0.04～0.04 mm，模型3的Uz分布在-0.05～0.05 mm.3個結(jié)構(gòu)的x向和y向的位移數(shù)值較大且處于同一數(shù)量級，z向位移較小.

圖9 不同垂直軸風(fēng)力機模型總位移云圖對比Fig.9 Comparison of total displacement contour map of different VAWT models

圖10 風(fēng)力機結(jié)構(gòu)三向位移云圖對比Fig.10 Comparison of three-way displacement contour maps of VAWT

由以上分析可知，原模型危險點出現(xiàn)在葉片頂端，故取節(jié)點8為對象，研究3個模型在風(fēng)機葉片頂端的風(fēng)致響應(yīng).又因x、y向位移數(shù)值較大，z向位移較小，故僅對比節(jié)點8的x、y向位移響應(yīng).圖11所示為3個模型節(jié)點8的風(fēng)致響應(yīng)對比.圖中Ux，8表示風(fēng)力機節(jié)點8在x方向的位移，Uy,8表示風(fēng)力機節(jié)點8在y方向的位移.可見，模型2和模型3在葉片頂端節(jié)點的位移抑制效果較為顯著.表2所示為3個模型的節(jié)點8的三向位移改善率統(tǒng)計數(shù)據(jù).可見，模型2在x方向的平均位移改善率約為44%，在y方向的平均位移改善率約為54%，在z方向的平均位移改善率約為5%.模型3在x方向的平均位移改善率約為40%，y方向的平均位移改善率約為53%，但在z方向未見明顯改善.

圖11 3個模型節(jié)點8的風(fēng)致響應(yīng)對比Fig.11 Comparison of wind-induced response of three models at node 8

表2 3個模型節(jié)點8的三向風(fēng)致響應(yīng)統(tǒng)計

5 結(jié)論

本文針對H型垂直軸風(fēng)力機振動問題，基于計算流體動力學(xué)方法，運用STAR-CCM+軟件數(shù)值模擬風(fēng)機旋轉(zhuǎn)時葉片所受氣動荷載，將氣動荷載在MATLAB軟件中實施函數(shù)擬合，獲得時序性風(fēng)機葉片荷載函數(shù)，再基于ANSYS軟件開展時程分析，得到風(fēng)機結(jié)構(gòu)的三向位移響應(yīng)時程.在風(fēng)機不同位置處施加基于COMBIN40單元而模擬的阻尼器，數(shù)值模擬計算其耗能減振能力.主要結(jié)論如下：

(1) H型垂直軸風(fēng)力機的風(fēng)致位移響應(yīng)曲線具有一定的周期性.t=10 s時，風(fēng)機最大位移出現(xiàn)在葉片的頂端，風(fēng)致x、y向位移處于同一數(shù)量級，z向位移較小.

(2) 在H型垂直軸風(fēng)力機的主軸與支桿連接處布置阻尼器可降低結(jié)構(gòu)位移響應(yīng)，總位移最大降幅達44%.

(3) 阻尼器布置的位置與結(jié)構(gòu)的位移改變密切相關(guān).t=10 s時，在近風(fēng)機葉片頂端連桿處布置阻尼器，結(jié)構(gòu)最大位移出現(xiàn)在風(fēng)機葉片底端，位移下降達44%；在近風(fēng)機葉片底端連桿處布置阻尼器，最大位移則依舊出現(xiàn)在了風(fēng)機葉片頂端，位移下降達40.7%.

本文成果可為垂直軸風(fēng)力機減振研究提供一定的技術(shù)參考.