占玲玉, 何文君, 周 岱,2, 韓兆龍,2, 朱宏博, 張 凱, 涂佳黃
(1. 上海交通大學(xué) 船舶海洋與建筑工程學(xué)院,上海 200240; 2. 上海交通大學(xué) 海洋工程國家重點實驗室,上海 200240; 3. 湘潭大學(xué) 土木工程與力學(xué)學(xué)院,湖南 湘潭 411105)
風(fēng)力機大型化可有效提高風(fēng)力發(fā)電量和風(fēng)能利用率,是海上風(fēng)力發(fā)電的世界發(fā)展趨勢[1].研發(fā)建造包括10、15和20 MW在內(nèi)的大型風(fēng)力機是我國海洋綠色能源開發(fā)利用的戰(zhàn)略重點.然而隨著風(fēng)力機轉(zhuǎn)子、塔筒尺寸的不斷增加[2],大型風(fēng)力機結(jié)構(gòu)系統(tǒng)柔性增加,抗風(fēng)能力明顯降低[3],風(fēng)力機風(fēng)致動力響應(yīng)更加復(fù)雜,系統(tǒng)風(fēng)致動力災(zāi)變風(fēng)險更加凸顯.
研究并優(yōu)化大型水平軸風(fēng)力機支撐結(jié)構(gòu)對保證系統(tǒng)運行安全十分重要.風(fēng)力機大型化導(dǎo)致塔筒頂部大型轉(zhuǎn)子和發(fā)電機質(zhì)量增加、支撐塔筒高度增大.常用的圓錐塔筒無法滿足風(fēng)力機結(jié)構(gòu)的強度和剛度要求.迄今,學(xué)者們對大型水平軸風(fēng)力機支撐結(jié)構(gòu)優(yōu)化已有一定研究.Gentils等[4]基于遺傳算法,考慮多準(zhǔn)則約束,對美國國家可再生能源實驗室(NREL)提出的5 MW單樁式風(fēng)力機(NREL 5 MW)塔架和基礎(chǔ)進行優(yōu)化分析并計算風(fēng)力機支撐結(jié)構(gòu)優(yōu)化尺寸參數(shù);Oest等[5]基于疲勞損傷荷載、準(zhǔn)靜態(tài)分析和動態(tài)分析,對海上風(fēng)力機支撐結(jié)構(gòu)進行梯度優(yōu)化設(shè)計;葛旭等[6]結(jié)合有限元分析與進化策略,對海上三樁單立柱風(fēng)力機支撐結(jié)構(gòu)主尺度參數(shù)進行多參數(shù)同步優(yōu)化設(shè)計;Tian等[7]對導(dǎo)管架式海上風(fēng)力機的導(dǎo)管架基礎(chǔ)進行拓撲優(yōu)化,結(jié)合尺寸和形狀優(yōu)化實現(xiàn)結(jié)構(gòu)重分布,得到質(zhì)量更輕的導(dǎo)管架結(jié)構(gòu).上述研究僅為針對特定風(fēng)力機結(jié)構(gòu)形式的參數(shù)優(yōu)化,未改變其支撐結(jié)構(gòu)形式,雖對風(fēng)力機結(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è)計目標(biāo);劉昆等[9]基于BESO算法對船體梁腹板開孔進行設(shè)計優(yōu)化,提出適用于船體結(jié)構(gòu)孔型的BESO算法;BESO還在羅馬小體育宮[10]、卡塔爾國家會議中心[11]等建筑工程設(shè)計中得到應(yīng)用.目前,大型水平軸風(fēng)力機塔筒支撐結(jié)構(gòu)拓撲優(yōu)化還鮮有研究.
基礎(chǔ)是風(fēng)力機系統(tǒng)的關(guān)鍵結(jié)構(gòu)之一,高樁承臺基礎(chǔ)的結(jié)構(gòu)安全性能高,適合我國海床地質(zhì)條件,施工工藝成熟,在我國近海風(fēng)電場建設(shè)中多有應(yīng)用[12].但不考慮樁-土相互作用會導(dǎo)致風(fēng)力機位移響應(yīng)被嚴(yán)重低估[13-15],是結(jié)構(gòu)安全設(shè)計的重要隱患.
考慮樁土相互作用的風(fēng)力機支撐結(jié)構(gòu)性能,開展支撐結(jié)構(gòu)優(yōu)化分析.基于反比例型刪除率的BESO算法,提出大型高樁承臺式水平軸風(fēng)力機支撐結(jié)構(gòu)拓撲的優(yōu)化算法,運用有限元軟件ANSYS開展大型水平軸風(fēng)力機支撐結(jié)構(gòu)尋優(yōu)分析,獲得一種風(fēng)力機優(yōu)化支撐結(jié)構(gòu),顯著降低了風(fēng)力機風(fēng)致動力響應(yīng).
以應(yīng)用較廣的NREL 5 MW[16]風(fēng)力機為背景,進行考慮樁土作用的近海坐底式高樁承臺大型水平軸風(fēng)力機塔筒結(jié)構(gòu)優(yōu)化研究.高樁承臺式水平軸風(fēng)力機初始模型如圖1(a)所示,其支撐系統(tǒng)包括上部塔筒結(jié)構(gòu)、承臺和群樁.水平軸風(fēng)力機塔筒為圓筒式支撐結(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,風(fēng)力機額定風(fēng)速為11.4 m/s,具體參數(shù)如表1所示.
表1 5 MW水平軸風(fēng)力機塔架參數(shù)Tab.1 5 MW horizontal axis wind turbine parameters
圖1 水平軸風(fēng)力機塔筒結(jié)構(gòu)初始模型Fig.1 Horizontal axis wind turbine initial model
參照上海東海大橋海上大型風(fēng)力發(fā)電工程項目采用高樁承臺式基礎(chǔ)[18],混凝土承臺直徑D1=20 m,承臺厚度S=5 m,設(shè)置4根直徑d1=3 m、壁厚t=20 mm的全直鋼管基樁,樁基長度H1=52 m,入土深度h=40 m.承臺材料選用C45標(biāo)號混凝土,彈性模量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
對高樁承臺大型水平軸風(fēng)力機結(jié)構(gòu)系統(tǒng),基于計算流體動力學(xué)(CFD)方法并運用Star-CCM+計算軟件[22]數(shù)值模擬獲得風(fēng)力機葉片氣動荷載,風(fēng)力機風(fēng)輪直徑為d,選取長方體風(fēng)場計算域,計算域長度、寬度和高度分別為20d、8d和8d,如圖3(a)和3(b)所示.計算域劃分共 1 700 萬個網(wǎng)格.計算域中,風(fēng)力機葉片圓柱體旋轉(zhuǎn)域在豎平面的直徑為1.5d,旋轉(zhuǎn)域高度即順風(fēng)向尺度為0.75d;旋轉(zhuǎn)域之外的計算域等其他部分為靜止域.對旋轉(zhuǎn)域和靜止域設(shè)置交界面,采用滑移網(wǎng)格處理風(fēng)輪運行旋轉(zhuǎn)問題,如圖3(c)和3(d)所示.葉片表面采用棱柱層網(wǎng)格,延伸因子為1.3,總厚度為0.08 m,無量綱化壁面距離為y+<1.
圖3 風(fēng)場計算域和邊界條件 Fig.3 Wind field computation domain and boundary condition
對大型風(fēng)力機支撐結(jié)構(gòu)采用錐筒結(jié)構(gòu)[23],塔筒高度為90 m.將塔筒簡化為受彎曲變形、軸向壓縮和扭轉(zhuǎn)變形的懸臂梁.風(fēng)力機運行時,需考慮塔筒沿高度的風(fēng)壓變化.沿塔筒高度l處的風(fēng)荷載[24]為
(1)
式中:Cd為風(fēng)阻力系數(shù),取值1.2;ρ為空氣密度;Dl為塔筒高度l處的塔筒截面直徑;V(l)為沿塔架高度l處的瞬時風(fēng)速.
美國石油協(xié)會(API)規(guī)范推薦的土體p-y曲線,即非線性彈簧模型廣泛應(yīng)用于海洋工程樁基設(shè)計.引入用于分析風(fēng)力機支撐結(jié)構(gòu)的樁土相互作用[25],該模型考慮了土體非線性、土層變化和土體種類等因素的影響,軟土中各層p-y曲線為
(2)
式中:p為深度X處單位樁長的極限水平土抗力標(biāo)準(zhǔn)值;yc為實驗室中對未擾動工試樣做不排水壓縮試驗時,其應(yīng)力達到最大應(yīng)力一半時的變形;y為泥面以下深度X處樁的側(cè)向水平變形;XR為泥面以下到土抗力減少區(qū)域底部的深度;pu為深度X處單位樁長的極限水平土抗力標(biāo)準(zhǔn)值.
對于砂土情況,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è)置豎向簡支約束.
采用 Newmark-β法[26]中的平均常加速度法對風(fēng)力機支撐結(jié)構(gòu)系統(tǒng)風(fēng)致振動進行瞬態(tài)分析,其中控制Newmark-β法精度和穩(wěn)定性的參數(shù)取值為γ=0.5,β=0.25.采用瑞利阻尼模擬風(fēng)力機結(jié)構(gòu)阻尼,剛度阻尼系數(shù)取值為0.03[27].風(fēng)力機結(jié)構(gòu)系統(tǒng)動力激勵下的振動方程[28]為
FI(t)+FC(t)+FK(t)
(4)
針對大型水平軸風(fēng)力機支撐結(jié)構(gòu)優(yōu)化問題,建立以結(jié)構(gòu)靜剛度最大為目標(biāo)、以結(jié)構(gòu)體積為約束條件的數(shù)學(xué)力學(xué)模型[29],表示為
(5)
(6)
式中:CT為結(jié)構(gòu)平均柔順度;V*為目標(biāo)體積比;vi為第i個單元的體積;n為迭代次數(shù);xi為設(shè)計變量;xmin和1分別代表某一單元在設(shè)計域內(nèi)缺失或存在.
為避免產(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)
鑒于水平軸風(fēng)力機結(jié)構(gòu)系統(tǒng)對稱性,將支撐結(jié)構(gòu)三維優(yōu)化問題轉(zhuǎn)化為二維平面優(yōu)化,以風(fēng)力機塔筒支撐形式和支撐點位置為優(yōu)化對象.將風(fēng)力機塔筒主軸高度方向和塔筒半徑方向圍成的平面作為初始設(shè)計域,采用ANSYS軟件中的4節(jié)點Plane182單元對計算模型進行網(wǎng)格劃分.設(shè)置底部為固定端約束,選取最不利氣動荷載作為外部風(fēng)荷載.采用反比例型刪除率優(yōu)化模型[32]對風(fēng)力機支撐結(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次迭代的體積比.
當(dāng)?shù)Y(jié)果滿足根據(jù)體積約束和收斂準(zhǔn)則[31]時,終止迭代.為提高迭代效率,設(shè)置收斂準(zhǔn)則,表示當(dāng)前迭代步之前10次迭代中平均柔順度的相對變化量,當(dāng)目標(biāo)函數(shù)變化隨著迭代變化足夠小時,則判定迭代結(jié)束:
(14)
式中:Ck為第k次迭代的平均柔順度;τ為收斂因子,取值為0.001.
根據(jù)上述計算過程實現(xiàn)基于BESO算法的水平軸風(fēng)力機支撐結(jié)構(gòu)優(yōu)化,流程如圖4所示.
圖4 拓撲優(yōu)化流程圖Fig.4 Flow chart of topology optimization
運用BESO雙向漸進結(jié)構(gòu)優(yōu)化算法,對5 MW大型水平軸風(fēng)力機塔筒結(jié)構(gòu)(高度為90 m)進行拓撲優(yōu)化,如圖5所示.為防止風(fēng)力機塔筒頂部發(fā)生過大位移,通常對風(fēng)力機塔筒受彎矩作用較大的底部區(qū)域進行加固,所選取的底部拓撲區(qū)域為寬度a=10 m、長度b=50 m所圍成的平面,以此作為初始設(shè)計域Ω,采用4節(jié)點Plane182單元進行網(wǎng)格劃分.在設(shè)計域區(qū)域施加外部荷載,包括風(fēng)力機機艙、轉(zhuǎn)子質(zhì)量等縱向荷載Fz,風(fēng)力機葉片運動產(chǎn)生的推力、沿塔筒高度分布的梯度風(fēng)荷載等在設(shè)計域頂部產(chǎn)生的荷載Fxy和彎矩Mxy.底部邊界條件為固定端約束,拓撲優(yōu)化的結(jié)構(gòu)目標(biāo)體積分?jǐn)?shù)為40%.
圖5 塔架結(jié)構(gòu)底部受力及拓撲設(shè)計區(qū)域Fig.5 Tower structure force diagram and topology design area
分析拓撲過程中體積比即單元刪除量變化和平均柔順度相對迭代次數(shù)的變化歷程,如圖6所示.拓撲迭代初期,拓撲單元較多,單元刪除量較大,拓撲迭代速率較快.隨著迭代進行,單元體積比逐漸減小,迭代20次后,結(jié)構(gòu)單元體積比開始趨近于目標(biāo)體積約束條件,拓撲形態(tài)開始穩(wěn)定,最終得到反比例型刪除率的優(yōu)化模型柔順度量值為1 557.1 N·m.
圖6 迭代過程結(jié)構(gòu)平均柔順度和體積分?jǐn)?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)建高樁承臺式水平軸風(fēng)力機支撐結(jié)構(gòu)模型.
圖7 塔筒支撐結(jié)構(gòu)拓撲優(yōu)化Fig.7 Topology optimization of tower support structures
基于CFD數(shù)值方法計算優(yōu)化后風(fēng)力機支撐結(jié)構(gòu)的風(fēng)致水平推力和風(fēng)力機功率.為考慮最不利荷載工況,在進行CFD數(shù)值計算時,設(shè)速度入口風(fēng)速等于切出風(fēng)速為25 m/s,與邊界元法(BEM)、致動盤(AD)[33]及有限體積法(FVM)[34]中的結(jié)果進行對比分析,如圖8所示.可知,本文風(fēng)力機功率、風(fēng)力機旋轉(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)的靜力特性,注意到風(fēng)力機塔筒經(jīng)受周期性脈動風(fēng)荷載作用,可將風(fēng)力機旋轉(zhuǎn)多圈并穩(wěn)定后葉片對塔筒的水平推力用于計算最不利工況下的塔筒位移, 運用式(1)計算沿風(fēng)機塔筒高度的風(fēng)荷載.風(fēng)力機支撐結(jié)構(gòu)順風(fēng)向(x向)、橫風(fēng)向(y向)和豎向(z向)最大動位移均發(fā)生于塔筒頂部,其中順風(fēng)向位移明顯大于橫風(fēng)向位移.針對5 MW大型水平軸風(fēng)力機,分別計算初始支撐結(jié)構(gòu)和優(yōu)化后支撐結(jié)構(gòu)在額定風(fēng)速下的塔頂位移.初始結(jié)構(gòu)塔頂順風(fēng)向位移為550 mm,豎向位移為-2.5 mm,而優(yōu)化后結(jié)構(gòu)塔頂?shù)捻橈L(fēng)向位移為250 mm,豎向位移為-1.7 mm.優(yōu)化后支撐結(jié)構(gòu)的位移明顯減少,順風(fēng)向和豎向位移相對初始支撐結(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
考慮樁土作用,分析5 MW風(fēng)力機初始和優(yōu)化后結(jié)構(gòu)在順風(fēng)向和豎向的動力響應(yīng).設(shè)來流風(fēng)速為25 m/s,風(fēng)力機葉片轉(zhuǎn)速為0.664 rad/s,風(fēng)力機運動周期為9.463 s.風(fēng)力機旋轉(zhuǎn)6~10圈后氣動性能基本穩(wěn)定,相應(yīng)的風(fēng)力機旋轉(zhuǎn)時長為50~200 s.初始和優(yōu)化后結(jié)構(gòu)的塔筒頂端順風(fēng)向位移(Ux)時程曲線如圖10所示.計算發(fā)現(xiàn)在t=68 s時刻,風(fēng)力機結(jié)構(gòu)振幅較大.因此分析該時刻優(yōu)化前后支撐結(jié)構(gòu)的順風(fēng)向和豎向位移,得到初始和優(yōu)化后塔頂順風(fēng)向最大位移分別為 1 130 和300 mm,減少73.45%;豎向位移分別為 -49~-26 mm和-18~-8 mm,最大豎向位移減少了63.27%.由此可見,在保持與初始結(jié)構(gòu)體積相同情況下,優(yōu)化后風(fēng)力機支撐結(jié)構(gòu)位移明顯減少,結(jié)構(gòu)剛度大于初始結(jié)構(gòu).
圖10 初始結(jié)構(gòu)和優(yōu)化結(jié)構(gòu)的塔頂順風(fēng)向位移時程曲線Fig.10 Comparison of displacement time history curve of two models at the top of tower
考慮樁土作用的初始和優(yōu)化后結(jié)構(gòu)在50~200 s時段順風(fēng)向的速度(vx)響應(yīng)時程和加速度(ax)響應(yīng)時程如圖11所示.可見優(yōu)化后結(jié)構(gòu)的速度和加速度響應(yīng)小于初始結(jié)構(gòu),風(fēng)振特性得到顯著改善,在風(fēng)力機旋轉(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)的塔頂順風(fēng)向速度時程和加速度時程Fig.11 Comparison of velocity and acceleration time history curve at the top of tower
針對大型高樁承臺式水平軸風(fēng)力機支撐結(jié)構(gòu)系統(tǒng),運用CFD方法計算風(fēng)致效應(yīng),考慮樁土作用效應(yīng),運用基于反比例函數(shù)的雙向漸進結(jié)構(gòu)優(yōu)化算法對風(fēng)力機塔筒支撐結(jié)構(gòu)進行拓撲優(yōu)化,開展風(fēng)力機結(jié)構(gòu)系統(tǒng)的靜動力響應(yīng)特性分析,得到結(jié)論如下:
(1) 運用基于反比例型刪除率的BESO算法并考慮樁土作用影響,開展水平軸風(fēng)力機塔筒支撐結(jié)構(gòu)拓撲優(yōu)化,得到風(fēng)力機支撐結(jié)構(gòu)的新形式.
(2) 優(yōu)化后支撐結(jié)構(gòu)風(fēng)致效應(yīng)明顯減小,在保持結(jié)構(gòu)體積不變情況下,優(yōu)化后支撐結(jié)構(gòu)整體剛度明顯提高,塔頂順風(fēng)向和豎向靜力位移分別減小54.5%和32%,動力分析下塔頂在順風(fēng)向和豎向最大位移分別減少73.45%和63.27%.
運用BESO算法進行水平軸風(fēng)力機支撐結(jié)構(gòu)的二維平面優(yōu)化,后續(xù)將對風(fēng)力機支撐結(jié)構(gòu)開展三維優(yōu)化研究.