基于BESO算法的高樁承臺式水平軸風力機支撐結(jié)構(gòu)優(yōu)化

占玲玉, 何文君, 周 岱,2, 韓兆龍,2, 朱宏博, 張 凱, 涂佳黃

(1. 上海交通大學(xué) 船舶海洋與建筑工程學(xué)院,上海 200240; 2. 上海交通大學(xué) 海洋工程國家重點實驗室,上海 200240; 3. 湘潭大學(xué) 土木工程與力學(xué)學(xué)院,湖南 湘潭 411105)

風力機大型化可有效提高風力發(fā)電量和風能利用率,是海上風力發(fā)電的世界發(fā)展趨勢[1].研發(fā)建造包括10、15和20 MW在內(nèi)的大型風力機是我國海洋綠色能源開發(fā)利用的戰(zhàn)略重點.然而隨著風力機轉(zhuǎn)子、塔筒尺寸的不斷增加[2],大型風力機結(jié)構(gòu)系統(tǒng)柔性增加,抗風能力明顯降低[3],風力機風致動力響應(yīng)更加復(fù)雜,系統(tǒng)風致動力災(zāi)變風險更加凸顯.

研究并優(yōu)化大型水平軸風力機支撐結(jié)構(gòu)對保證系統(tǒng)運行安全十分重要.風力機大型化導(dǎo)致塔筒頂部大型轉(zhuǎn)子和發(fā)電機質(zhì)量增加、支撐塔筒高度增大.常用的圓錐塔筒無法滿足風力機結(jié)構(gòu)的強度和剛度要求.迄今,學(xué)者們對大型水平軸風力機支撐結(jié)構(gòu)優(yōu)化已有一定研究.Gentils等[4]基于遺傳算法,考慮多準則約束,對美國國家可再生能源實驗室(NREL)提出的5 MW單樁式風力機(NREL 5 MW)塔架和基礎(chǔ)進行優(yōu)化分析并計算風力機支撐結(jié)構(gòu)優(yōu)化尺寸參數(shù);Oest等[5]基于疲勞損傷荷載、準靜態(tài)分析和動態(tài)分析,對海上風力機支撐結(jié)構(gòu)進行梯度優(yōu)化設(shè)計;葛旭等[6]結(jié)合有限元分析與進化策略,對海上三樁單立柱風力機支撐結(jié)構(gòu)主尺度參數(shù)進行多參數(shù)同步優(yōu)化設(shè)計;Tian等[7]對導(dǎo)管架式海上風力機的導(dǎo)管架基礎(chǔ)進行拓撲優(yōu)化,結(jié)合尺寸和形狀優(yōu)化實現(xiàn)結(jié)構(gòu)重分布,得到質(zhì)量更輕的導(dǎo)管架結(jié)構(gòu).上述研究僅為針對特定風力機結(jié)構(gòu)形式的參數(shù)優(yōu)化,未改變其支撐結(jié)構(gòu)形式,雖對風力機結(jié)構(gòu)穩(wěn)定性有一定提高,但設(shè)計方法具有局限性.相較于參數(shù)優(yōu)化,拓撲優(yōu)化方法具有更多的設(shè)計自由度,是優(yōu)化大型水平軸支撐結(jié)構(gòu)的新途徑.

拓撲優(yōu)化方法應(yīng)用于包括結(jié)構(gòu)優(yōu)化在內(nèi)的廣泛領(lǐng)域,其雙向漸進結(jié)構(gòu)優(yōu)化(BESO)算法優(yōu)化效果較好.王章駿等[8]采用BESO算法對吸能管進行拓撲優(yōu)化,提高了帶隔板薄壁方管的耐撞性能, 實現(xiàn)了吸能管輕量化的設(shè)計目標;劉昆等[9]基于BESO算法對船體梁腹板開孔進行設(shè)計優(yōu)化,提出適用于船體結(jié)構(gòu)孔型的BESO算法;BESO還在羅馬小體育宮[10]、卡塔爾國家會議中心[11]等建筑工程設(shè)計中得到應(yīng)用.目前,大型水平軸風力機塔筒支撐結(jié)構(gòu)拓撲優(yōu)化還鮮有研究.

基礎(chǔ)是風力機系統(tǒng)的關(guān)鍵結(jié)構(gòu)之一,高樁承臺基礎(chǔ)的結(jié)構(gòu)安全性能高,適合我國海床地質(zhì)條件,施工工藝成熟,在我國近海風電場建設(shè)中多有應(yīng)用[12].但不考慮樁-土相互作用會導(dǎo)致風力機位移響應(yīng)被嚴重低估[13-15],是結(jié)構(gòu)安全設(shè)計的重要隱患.

考慮樁土相互作用的風力機支撐結(jié)構(gòu)性能,開展支撐結(jié)構(gòu)優(yōu)化分析.基于反比例型刪除率的BESO算法,提出大型高樁承臺式水平軸風力機支撐結(jié)構(gòu)拓撲的優(yōu)化算法,運用有限元軟件ANSYS開展大型水平軸風力機支撐結(jié)構(gòu)尋優(yōu)分析,獲得一種風力機優(yōu)化支撐結(jié)構(gòu),顯著降低了風力機風致動力響應(yīng).

1 風力機結(jié)構(gòu)建模與分析

1.1 風力機支撐結(jié)構(gòu)建模

以應(yīng)用較廣的NREL 5 MW[16]風力機為背景,進行考慮樁土作用的近海坐底式高樁承臺大型水平軸風力機塔筒結(jié)構(gòu)優(yōu)化研究.高樁承臺式水平軸風力機初始模型如圖1(a)所示,其支撐系統(tǒng)包括上部塔筒結(jié)構(gòu)、承臺和群樁.水平軸風力機塔筒為圓筒式支撐結(jié)構(gòu),如圖1(b)所示,塔筒底部外徑D=6 m,塔筒高度H0=90 m[16].塔筒選用Q345鋼材[17],密度ρ0=7.85 g/m3,彈性模量E0=211 N/m2,泊松比ν0=0.3,風力機額定風速為11.4 m/s,具體參數(shù)如表1所示.

表1 5 MW水平軸風力機塔架參數(shù)Tab.1 5 MW horizontal axis wind turbine parameters

圖1 水平軸風力機塔筒結(jié)構(gòu)初始模型Fig.1 Horizontal axis wind turbine initial model

參照上海東海大橋海上大型風力發(fā)電工程項目采用高樁承臺式基礎(chǔ)[18],混凝土承臺直徑D1=20 m,承臺厚度S=5 m,設(shè)置4根直徑d1=3 m、壁厚t=20 mm的全直鋼管基樁,樁基長度H1=52 m,入土深度h=40 m.承臺材料選用C45標號混凝土,彈性模量E=33.5 GPa,密度為2 500 kg/m3,泊松比ν=0.2[19],基礎(chǔ)模型如圖2(a)所示.基于Jung等[20]和Chen等[21]確定土體模型和力學(xué)參數(shù):土體浮重度γ′、彈性模量E′、泊松比ν′、內(nèi)摩擦角φ′和不排水抗剪強度su、黏土泊松比νu、50%的最大有效應(yīng)力下的軸向應(yīng)變ε50、土層深度z、膨脹角ψ′、黏聚力c′、地基反力模量k,如圖2(b)所示.

圖2 基礎(chǔ)與土體模型Fig.2 Foundation geometric model and soil profiles

1.2 風力機支撐結(jié)構(gòu)風荷載

對高樁承臺大型水平軸風力機結(jié)構(gòu)系統(tǒng),基于計算流體動力學(xué)(CFD)方法并運用Star-CCM+計算軟件[22]數(shù)值模擬獲得風力機葉片氣動荷載,風力機風輪直徑為d,選取長方體風場計算域,計算域長度、寬度和高度分別為20d、8d和8d,如圖3(a)和3(b)所示.計算域劃分共 1 700 萬個網(wǎng)格.計算域中,風力機葉片圓柱體旋轉(zhuǎn)域在豎平面的直徑為1.5d,旋轉(zhuǎn)域高度即順風向尺度為0.75d;旋轉(zhuǎn)域之外的計算域等其他部分為靜止域.對旋轉(zhuǎn)域和靜止域設(shè)置交界面,采用滑移網(wǎng)格處理風輪運行旋轉(zhuǎn)問題,如圖3(c)和3(d)所示.葉片表面采用棱柱層網(wǎng)格,延伸因子為1.3,總厚度為0.08 m,無量綱化壁面距離為y+<1.

圖3 風場計算域和邊界條件 Fig.3 Wind field computation domain and boundary condition

對大型風力機支撐結(jié)構(gòu)采用錐筒結(jié)構(gòu)[23],塔筒高度為90 m.將塔筒簡化為受彎曲變形、軸向壓縮和扭轉(zhuǎn)變形的懸臂梁.風力機運行時,需考慮塔筒沿高度的風壓變化.沿塔筒高度l處的風荷載[24]為

(1)

式中:Cd為風阻力系數(shù),取值1.2;ρ為空氣密度;Dl為塔筒高度l處的塔筒截面直徑;V(l)為沿塔架高度l處的瞬時風速.

1.3 樁土相互作用

美國石油協(xié)會(API)規(guī)范推薦的土體p-y曲線,即非線性彈簧模型廣泛應(yīng)用于海洋工程樁基設(shè)計.引入用于分析風力機支撐結(jié)構(gòu)的樁土相互作用[25],該模型考慮了土體非線性、土層變化和土體種類等因素的影響,軟土中各層p-y曲線為

(2)

式中:p為深度X處單位樁長的極限水平土抗力標準值;yc為實驗室中對未擾動工試樣做不排水壓縮試驗時,其應(yīng)力達到最大應(yīng)力一半時的變形;y為泥面以下深度X處樁的側(cè)向水平變形;XR為泥面以下到土抗力減少區(qū)域底部的深度;pu為深度X處單位樁長的極限水平土抗力標準值.

對于砂土情況,p-y曲線為

(3)

式中:K為地基反力初始模量;A為計入靜力荷載和循環(huán)荷載條件的參數(shù).

基于API規(guī)范的p-y曲線法,計算循環(huán)荷載作用情況下的樁側(cè)土抗力,使用ANSYS軟件的Combine39單元模擬p-y曲線,該單元是具有非線性功能的單向拉壓單元,常用于樁土相互作用模擬.在搭建有限元模型時,根據(jù)計算得到的土體p-y曲線參數(shù),每間隔1 m設(shè)置一個Combine39土彈簧單元.在樁土相互作用模擬過程中,為簡化計算未考慮樁土循環(huán)弱化作用及樁端設(shè)置豎向簡支約束.

1.4 風力機支撐結(jié)構(gòu)系統(tǒng)振動響應(yīng)

采用 Newmark-β法[26]中的平均常加速度法對風力機支撐結(jié)構(gòu)系統(tǒng)風致振動進行瞬態(tài)分析,其中控制Newmark-β法精度和穩(wěn)定性的參數(shù)取值為γ=0.5,β=0.25.采用瑞利阻尼模擬風力機結(jié)構(gòu)阻尼,剛度阻尼系數(shù)取值為0.03[27].風力機結(jié)構(gòu)系統(tǒng)動力激勵下的振動方程[28]為

FI(t)+FC(t)+FK(t)

(4)

2 大型水平軸風力機支撐拓撲優(yōu)化

2.1 BESO算法的數(shù)學(xué)模型

針對大型水平軸風力機支撐結(jié)構(gòu)優(yōu)化問題,建立以結(jié)構(gòu)靜剛度最大為目標、以結(jié)構(gòu)體積為約束條件的數(shù)學(xué)力學(xué)模型[29],表示為

(5)

(6)

式中:CT為結(jié)構(gòu)平均柔順度;V*為目標體積比;vi為第i個單元的體積;n為迭代次數(shù);xi為設(shè)計變量;xmin和1分別代表某一單元在設(shè)計域內(nèi)缺失或存在.

2.2 單元靈敏度過濾及更新

為避免產(chǎn)生奇異單元剛度矩陣,算法以xmin代替0表示空單元,在該軟刪單元的BESO算法[30]中所有單元的材料彈性模量可根據(jù)材料插值模型設(shè)置為

(7)

式中:q為懲罰因子.

單元靈敏度表示增刪單元對結(jié)構(gòu)平均柔順度的影響程度,單元靈敏度大小為單元靈敏度的相對排序,表示為

(8)

式中:Ki和ui分別為單元剛度矩陣和單元位移向量.

為抑制優(yōu)化過程中出現(xiàn)實體單元與空洞單元交替出現(xiàn)的棋盤格等數(shù)值不穩(wěn)定現(xiàn)象[31],通過權(quán)重系數(shù)引入節(jié)點靈敏度節(jié)點概念:

(9)

在迭代過程中,根據(jù)節(jié)點靈敏度對單元靈敏度進行過濾并更新單元靈敏度,更新過程為

(10)

(11)

2.3 基于BESO算法的結(jié)構(gòu)拓撲優(yōu)化

鑒于水平軸風力機結(jié)構(gòu)系統(tǒng)對稱性,將支撐結(jié)構(gòu)三維優(yōu)化問題轉(zhuǎn)化為二維平面優(yōu)化,以風力機塔筒支撐形式和支撐點位置為優(yōu)化對象.將風力機塔筒主軸高度方向和塔筒半徑方向圍成的平面作為初始設(shè)計域,采用ANSYS軟件中的4節(jié)點Plane182單元對計算模型進行網(wǎng)格劃分.設(shè)置底部為固定端約束,選取最不利氣動荷載作為外部風荷載.采用反比例型刪除率優(yōu)化模型[32]對風力機支撐結(jié)構(gòu)進行拓撲優(yōu)化:

(12)

式中:ω為結(jié)構(gòu)材料的刪除率;ωmax為給定的最大刪除率,取0.057;ωmin為給定的最小刪除率,取0.012;Vi為優(yōu)化過程中每一步迭代過程的體積比;V*值取0.4.

在迭代過程中,根據(jù)刪除率對體積比進行更新

Vk+1=Vk(1±ω),k=1, 2, 3, …

(13)

式中:Vk為第k次迭代的體積比:Vk+1為第k+1次迭代的體積比.

當?shù)Y(jié)果滿足根據(jù)體積約束和收斂準則[31]時,終止迭代.為提高迭代效率,設(shè)置收斂準則,表示當前迭代步之前10次迭代中平均柔順度的相對變化量,當目標函數(shù)變化隨著迭代變化足夠小時,則判定迭代結(jié)束:

(14)

式中:Ck為第k次迭代的平均柔順度;τ為收斂因子,取值為0.001.

根據(jù)上述計算過程實現(xiàn)基于BESO算法的水平軸風力機支撐結(jié)構(gòu)優(yōu)化,流程如圖4所示.

圖4 拓撲優(yōu)化流程圖Fig.4 Flow chart of topology optimization

3 優(yōu)化前后模型受力特性

3.1 基于BESO算法的支撐結(jié)構(gòu)優(yōu)化

運用BESO雙向漸進結(jié)構(gòu)優(yōu)化算法,對5 MW大型水平軸風力機塔筒結(jié)構(gòu)(高度為90 m)進行拓撲優(yōu)化,如圖5所示.為防止風力機塔筒頂部發(fā)生過大位移,通常對風力機塔筒受彎矩作用較大的底部區(qū)域進行加固,所選取的底部拓撲區(qū)域為寬度a=10 m、長度b=50 m所圍成的平面,以此作為初始設(shè)計域Ω,采用4節(jié)點Plane182單元進行網(wǎng)格劃分.在設(shè)計域區(qū)域施加外部荷載,包括風力機機艙、轉(zhuǎn)子質(zhì)量等縱向荷載Fz,風力機葉片運動產(chǎn)生的推力、沿塔筒高度分布的梯度風荷載等在設(shè)計域頂部產(chǎn)生的荷載Fxy和彎矩Mxy.底部邊界條件為固定端約束,拓撲優(yōu)化的結(jié)構(gòu)目標體積分數(shù)為40%.

圖5 塔架結(jié)構(gòu)底部受力及拓撲設(shè)計區(qū)域Fig.5 Tower structure force diagram and topology design area

分析拓撲過程中體積比即單元刪除量變化和平均柔順度相對迭代次數(shù)的變化歷程,如圖6所示.拓撲迭代初期,拓撲單元較多,單元刪除量較大,拓撲迭代速率較快.隨著迭代進行,單元體積比逐漸減小,迭代20次后,結(jié)構(gòu)單元體積比開始趨近于目標體積約束條件,拓撲形態(tài)開始穩(wěn)定,最終得到反比例型刪除率的優(yōu)化模型柔順度量值為1 557.1 N·m.

圖6 迭代過程結(jié)構(gòu)平均柔順度和體積分數(shù)變化Fig.6 Average flexibility and volume fraction of optimized structure

對初始設(shè)計域,通過上述二維拓撲找形得到帶有側(cè)向支撐的塔筒形狀如圖7(a)所示,其側(cè)向支撐高度為48 m、寬度為8 m.基于支撐結(jié)構(gòu)的對稱性,考慮基礎(chǔ)為四樁布局,將二維拓撲結(jié)果擴展到三維空間,獲得的四支撐結(jié)構(gòu)優(yōu)化形式如圖7(b)所示.根據(jù)拓撲優(yōu)化結(jié)果,構(gòu)建高樁承臺式水平軸風力機支撐結(jié)構(gòu)模型.

圖7 塔筒支撐結(jié)構(gòu)拓撲優(yōu)化Fig.7 Topology optimization of tower support structures

3.2 優(yōu)化后風力機支撐結(jié)構(gòu)特性

基于CFD數(shù)值方法計算優(yōu)化后風力機支撐結(jié)構(gòu)的風致水平推力和風力機功率.為考慮最不利荷載工況,在進行CFD數(shù)值計算時,設(shè)速度入口風速等于切出風速為25 m/s,與邊界元法(BEM)、致動盤(AD)[33]及有限體積法(FVM)[34]中的結(jié)果進行對比分析,如圖8所示.可知,本文風力機功率、風力機旋轉(zhuǎn)時葉片對塔筒的水平推力(T=750 kN)均與相關(guān)文獻結(jié)果接近.

圖8 本文結(jié)果及其分析對比Fig.8 Power comparison between CFD model and reference

分別計算優(yōu)化后支撐結(jié)構(gòu)和原塔筒結(jié)構(gòu)的靜力特性,注意到風力機塔筒經(jīng)受周期性脈動風荷載作用,可將風力機旋轉(zhuǎn)多圈并穩(wěn)定后葉片對塔筒的水平推力用于計算最不利工況下的塔筒位移, 運用式(1)計算沿風機塔筒高度的風荷載.風力機支撐結(jié)構(gòu)順風向(x向)、橫風向(y向)和豎向(z向)最大動位移均發(fā)生于塔筒頂部,其中順風向位移明顯大于橫風向位移.針對5 MW大型水平軸風力機,分別計算初始支撐結(jié)構(gòu)和優(yōu)化后支撐結(jié)構(gòu)在額定風速下的塔頂位移.初始結(jié)構(gòu)塔頂順風向位移為550 mm,豎向位移為-2.5 mm,而優(yōu)化后結(jié)構(gòu)塔頂?shù)捻橈L向位移為250 mm,豎向位移為-1.7 mm.優(yōu)化后支撐結(jié)構(gòu)的位移明顯減少,順風向和豎向位移相對初始支撐結(jié)構(gòu)減少了54.5%和32%,如圖9所示.

圖9 初始結(jié)構(gòu)與優(yōu)化結(jié)構(gòu)位移響應(yīng)云圖Fig.9 Comparison of displacement of initial structure and optimized structure

3.3 樁土作用優(yōu)化后支撐結(jié)構(gòu)風致動力響應(yīng)分析

考慮樁土作用,分析5 MW風力機初始和優(yōu)化后結(jié)構(gòu)在順風向和豎向的動力響應(yīng).設(shè)來流風速為25 m/s,風力機葉片轉(zhuǎn)速為0.664 rad/s,風力機運動周期為9.463 s.風力機旋轉(zhuǎn)6～10圈后氣動性能基本穩(wěn)定,相應(yīng)的風力機旋轉(zhuǎn)時長為50～200 s.初始和優(yōu)化后結(jié)構(gòu)的塔筒頂端順風向位移(Ux)時程曲線如圖10所示.計算發(fā)現(xiàn)在t=68 s時刻,風力機結(jié)構(gòu)振幅較大.因此分析該時刻優(yōu)化前后支撐結(jié)構(gòu)的順風向和豎向位移,得到初始和優(yōu)化后塔頂順風向最大位移分別為 1 130 和300 mm,減少73.45%;豎向位移分別為 -49～-26 mm和-18～-8 mm,最大豎向位移減少了63.27%.由此可見,在保持與初始結(jié)構(gòu)體積相同情況下,優(yōu)化后風力機支撐結(jié)構(gòu)位移明顯減少,結(jié)構(gòu)剛度大于初始結(jié)構(gòu).

圖10 初始結(jié)構(gòu)和優(yōu)化結(jié)構(gòu)的塔頂順風向位移時程曲線Fig.10 Comparison of displacement time history curve of two models at the top of tower

考慮樁土作用的初始和優(yōu)化后結(jié)構(gòu)在50～200 s時段順風向的速度(vx)響應(yīng)時程和加速度(ax)響應(yīng)時程如圖11所示.可見優(yōu)化后結(jié)構(gòu)的速度和加速度響應(yīng)小于初始結(jié)構(gòu),風振特性得到顯著改善,在風力機旋轉(zhuǎn)多圈如100 s進入穩(wěn)態(tài)后改善效果更明顯,如塔頂加速度由初始結(jié)構(gòu)的5.02 m/s2減小為優(yōu)化后結(jié)構(gòu)的3.8 m/s2,減幅達24%.考慮樁土作用的支撐結(jié)構(gòu)振動效應(yīng)明顯減小,具有較好的減振效果.

圖11 初始和優(yōu)化結(jié)構(gòu)的塔頂順風向速度時程和加速度時程Fig.11 Comparison of velocity and acceleration time history curve at the top of tower

4 結(jié)論

針對大型高樁承臺式水平軸風力機支撐結(jié)構(gòu)系統(tǒng),運用CFD方法計算風致效應(yīng),考慮樁土作用效應(yīng),運用基于反比例函數(shù)的雙向漸進結(jié)構(gòu)優(yōu)化算法對風力機塔筒支撐結(jié)構(gòu)進行拓撲優(yōu)化,開展風力機結(jié)構(gòu)系統(tǒng)的靜動力響應(yīng)特性分析,得到結(jié)論如下:

(1) 運用基于反比例型刪除率的BESO算法并考慮樁土作用影響,開展水平軸風力機塔筒支撐結(jié)構(gòu)拓撲優(yōu)化,得到風力機支撐結(jié)構(gòu)的新形式.

(2) 優(yōu)化后支撐結(jié)構(gòu)風致效應(yīng)明顯減小,在保持結(jié)構(gòu)體積不變情況下,優(yōu)化后支撐結(jié)構(gòu)整體剛度明顯提高,塔頂順風向和豎向靜力位移分別減小54.5%和32%,動力分析下塔頂在順風向和豎向最大位移分別減少73.45%和63.27%.

運用BESO算法進行水平軸風力機支撐結(jié)構(gòu)的二維平面優(yōu)化,后續(xù)將對風力機支撐結(jié)構(gòu)開展三維優(yōu)化研究.