基于DEM-FEM 耦合方法的海上風機結構冰激振動分析1)

楊冬寶 高俊松 劉建平 宋 礎 季順迎,2)

?(大連理工大學工業裝備結構分析國家重點實驗室，遼寧大連 116023)

?(中國三峽新能源(集團)股份有限公司，北京 100053)

??(上海勘測設計研究院有限公司，上海 200434)

引言

海上風能是最具有發展前景的可再生資源之一.我國漫長的海岸線和豐富的海風資源，為發展海上風電提供了有力的條件.雖然我國海上風電發展起步晚，但是自2018 年以來，我國海上風電新增裝機量已經超過英國和德國，成為全球最大的海上風電發展市場.

傳統的海上風機設計主要考慮空氣動力、水利環境和抗震設計等外界因素[1-2]，但是對于我國渤海及黃海北部寒區海域的風機設計，海冰也是無法忽略的重要因素之一.雖然在一些寒區風機抗冰設計中已經采用極端靜冰載荷校核風機結構抗冰性能，但是大量現場觀測和理論研究表明，由海冰引起的冰激振動對海洋結構和設備的危害遠遠超過極端靜載荷下結構整體安全問題[3].此外，現場測量和模型實驗表明直立或錐體基礎的海洋平臺結構都會發生強烈的冰激振動，相比于錐體結構的隨機振動，直立結構的穩態振動對于結構的危害更為嚴重[4-5].

海上風機結構主要分為適用于淺水海域的基礎固定式風機和適用于較深海域的浮式風機[6].在各種風機結構中，單樁式風機是目前國內外應用最為廣泛、研究最為集中的風機結構.對于風機與海冰的相互作用，其研究工作最早開始于20 世紀80年代丹麥、芬蘭和挪威等北歐國家[7].通過一系列單樁風機基礎與海冰相互作用的模型試驗探討了不同風機基礎結構模型的動力學特性和冰載荷特點[8-9].近些年來，隨著數值模擬技術的快速發展，Wang[10]采用Karna 冰力譜模型分析了冰載荷與風載荷對于風機疲勞的影響，得到冰載荷對于風機結構的振動響應和疲勞損傷的影響大于風載荷所造成的影響的結論；Shi 等[11]將半經驗海冰--結構作用模型與HAWC2 程序相耦合研究了具有抗冰錐的單樁式風機在不同冰速、冰厚工況下風機結構的動力響應.此外，美國能源部國家可再生能源實驗室(NREL)的FAST 風機軟件也引入了基于M??tt?nen-Blenkar 海冰模型IceDyn 模塊計算風機的冰載荷[12].相比于國外，國內對于風機研究主要通過數值模擬和室內模型試驗進行研究[13].天津大學在低溫冰水池對于不同直徑單樁風機基礎開展模型試驗，重點驗證各種冰載荷規范對于單樁風機適用性[7].一些國內學者基于風力譜和冰力函數也對各種固定式風機進行振動時域響應分析[14-16].

基于經典的海冰力學模型或冰力函數研究海冰與風機結構的相互作用更多側重于結構的響應，而對海冰破壞方式給予了大量的假設.近年來，基于非連續力學的離散元方法被廣泛地應用于海冰數值模擬中來有效模擬海冰破碎現象[17-18]，其中采用粘結--破碎特性的球體離散單元可以合理準確地模擬海冰與復雜海洋平臺結構、船體結構等相互作用過程[19-20].風機結構在海冰激勵下復雜的動力學響應計算可以采用有限元方法計算，因此基于計算參數傳遞的DEM-FEM 界面耦合模型[21]可計算單樁式風機結構的冰激振動.

為此，本文采用DEM-FEM 耦合方法建立海冰與單樁式風機的耦合模型，探討不同冰速、冰厚下單樁式風機的冰載荷及風機輪轂和下部樁基的冰激振動響應特征.

1 單樁式風機冰激振動的數值方法

在對風機整體結構的冰激振動特性模擬中，風機采用有限元方法進行動力學分析，海冰采用離散元方法模擬，兩者通過計算參數傳遞實現耦合分析.

1.1 海冰離散元方法

本文主要研究平整冰與風機相互作用過程，其中將平整冰模型離散為具有粘結--破碎特性的等直徑球體單元.為合理地模擬實際海域內海冰的邊界條件，對DEM 海冰模型的兩側邊界的球形單元給予y和z方向上固定的位移約束，后側給予恒定的流速，如圖1(a)所示.顆粒之間采用平行粘結模型[22]傳遞單元間的作用力，如圖1(b)所示.平行粘結模型之間的最大拉應力和最大剪應力可以通過經典梁的拉伸、扭轉和彎曲理論計算，即

圖1 海冰離散元方法Fig.1 DEM of sea ice

式中，Fn和Fs為粘結顆粒之間的法向和切向力；Mn和Ms為粘結顆粒的法向和切向力矩；A，R，J和I分別為粘結顆粒之間的粘結面積、等效粘結半徑、極慣性矩和慣性矩.

海冰離散元方法通過單元間的粘結失效模擬海冰內部裂紋的生成和擴展，其中采用拉剪分區斷裂準則對顆粒粘結失效進行判斷[23]，如圖2 所示.當兩個粘結顆粒之間的最大拉應力σmax或剪切應力τmax達到材料拉伸破壞強度σt或剪切破壞強度τt時，顆粒間的粘結失效海冰發生破壞.這里，粘結顆粒間的拉伸破壞強度σt和剪切破壞強度τt表示為

圖2 拉剪分區斷裂準則Fig.2 Fracture criterion with tensile and shear partition

海冰離散元參數σb和μb的選取直接影響海冰宏觀物理力學性質.通過DEM 模擬海冰的單軸壓縮試驗和三點彎曲試驗，可以得到海冰宏觀強度與顆粒直徑D、粘結強度σb、顆粒間摩擦系數μb等離散元微觀參數的關系，其可寫作[24-25]

式中，σf和σc為海冰的彎曲強度和壓縮強度；尺寸L可視為海冰的厚度，則L/D為平整冰中單元層數.

1.2 單樁式風機結構的有限元模型

建立帶有抗冰錐體的單樁式風機的有限元模型，該模型主要由樁基、抗冰錐、塔筒和輪轂葉片4 部分組成，如圖3(a)所示.為保證單樁式風機結構幾何形狀和整體結構振動響應的真實性，對風機的樁基和上部輪轂葉片作了適當的簡化.風機樁基和塔筒等主體部分由梁單元構造；葉片由三角形單元構造，葉片單元僅為保證結構模型幾何完整性不做動力學計算.風機抗冰錐體采用三角形平板殼單元進行有限元分析計算，其結構模型如圖3(b)所示.綜上所述，本文采用的單樁式風機模型主要由梁--殼組合單元構造，共有1229 個單元節點和1933個單元，風機各部分詳細信息列于表1 中.

圖3 單樁式風機有限元模型Fig.3 FEM of monopile-type wind turbine

表1 風機結構有限元模型的主要參數Table 1 Main parameters of finit element model Information of OWTs

為重點考慮風機在海冰作用下的振動特性，忽略空氣動力載荷和水動力載荷，并且采用6 倍樁徑法建立風機樁--土系統約束.此外風機上部的輪轂葉片的原型質量為213.7 t，簡化為一個集中質量單元.表2 列出簡化后模型的前5 階自振頻率.

表2 單樁式風機振動模態Table 2 Vibration mode of monopile-type wind turbine

采用Newmark 隱式積分求解算法對于風機動力學方程求解，結構的動力學方程為

式中，參數α 和β 通過結構的模態阻尼比和固有頻率計算得到，即α=0.171，β=0.014 6.

1.3 DEM-FEM 耦合模型

本文所采用的DEM-FEM 耦合算法是基于原耦合異步方法的基礎上[21].通過CUDA C++語言建立CPU-GPU 協同處理并行計算架構的DEM-FEM界面耦合模型.基于界面參數傳遞的DEM-FEM 耦合算法在有限元和離散元計算中需要二者進行參數信息傳遞達到強耦合的效果.在有限元計算中，需要實時獲得離散單元與結構之間的接觸力、接觸單元編號和局部接觸點等信息；同時在離散元計算中，又需要更新結構的坐標和單元速度信息，如圖4 所示.但是，GPU 設備端的DEM 并行計算所得到載荷信息結果是隨機分配在各個線程中，使得界面之間參數傳遞效率低下.本文采用CUDA 的Thrust 函數中inclusive scan 和sortbykey 函數對界面之間的參數進行高效的基數排序[21]，實現離散元與有限元之間高效的參數傳遞.

圖4 基于GPU 并行的DEM-FEM 耦合算法框架Fig.4 Schematic of a GPU-based parallel algorithm of coupled DEM-FEM

2 單樁式風機冰載荷分析及驗證

2.1 風機冰載荷分析

海洋風機的冰載荷和冰激振動與海冰的厚度和速度有著密切關系[12].針對50 年一遇海冰厚度為0.374 m，平均冰速為0.3 m/s[15,26]，本文依照冰速和冰厚兩個變量劃分20 種工況，具體參數列于表3.此外，在流體動力學方面，考慮水對海冰的拖曳力Fd和浮力Fb，表示為

式中，ρw為海水密度；為單元i浸入水體積；Cd為流對冰的拖曳系數；和vw為海冰和水的速度；z為向上單位矢量.

表3 數值模擬中的工況劃分Table 3 Working condition of numerical simulation

風機冰激振動的相關計算參數列于表4.當冰速v=0.2 m/s,冰厚h=0.2 m 時，海冰與風機結構相互作用過程如圖5 所示.海冰主要與風機錐體上錐作用發生彎曲破壞，破碎冰排發生攀爬與堆積現象，由此引發結構的冰激振動.通過DEM-FEM 耦合方法模擬平整海冰與單樁式風機結構相互作用，可以得到風機承受的冰載荷時程，如圖6 所示.該圖表明風機冰載荷呈現很強的隨機性與周期性.這是由于海冰與風機錐體作用發生彎曲破壞時，破碎冰排與錐體結構會產生周期性的攀爬、堆積和清除，并且冰排攀爬高度和堆積體積表現隨機性.由于冰載荷隨機性的存在，很難與現場實測數據或者經驗公式直接對比，因此本文采取脈沖統計方法[27]，統計數據峰值的均值作為冰載荷時程的冰載荷表征量.

表4 DEM-FEM 模擬中計算參數Table 4 Calculation parameters of DEM-FEM simulation

圖5 DEM-FEM 模擬海冰與單樁式風機相互作用Fig.5 Interaction between ice and monopile-type wind turbine of DEM-FEM simulation

圖6 風機冰載荷時程曲線Fig.6 Time history curve of ice load on the wind turbine

2.2 風機冰載荷驗證

為驗證模擬結果的準確性，對不同工況下的冰載荷峰值的均值與ISO-19906 標準[28]、IEC6100-3 規范[29]對比.對于帶有錐體的單樁式風機結構,兩個規范分別給出典型錐體結構的靜冰力公式.在IEC6100-3 規范中，采用三維窄結構的的Ralstion 公式，其海冰與正錐體相互作用時冰載荷為[29]

式中，A1，A2，A3和A4為無量綱參數，與海冰與錐體的摩擦系數μ和錐角α 有關；D為水線處錐體直徑；DT為錐頂直徑，即塔筒直徑.

ISO-19906 標準作為世界石油天然氣工業領域中寒區海洋結構物設計和評估的行業標準，其基于海冰塑性理論給出平整冰與錐體結構的水平方向冰載荷公式[28]

表5 列出ISO-19906 標準和IEC6100-3 規范中經驗公式中相關計算參數的具體取值.

表5 ISO-19906 標準和IEC6100-3 規范錐體結構冰載荷計算參數Table 5 Computational parameters of ice load on conical structure based on codes ISO-19906 and IEC6100-3

圖7 DEM-FEM 計算冰載荷峰值與經驗公式對比Fig.7 Comparison between the peck value of ice load in the DEM-FEM simulation and the empirical formulas

圖7 所示為DEM-FEM 耦合方法計算的冰載荷峰值的均值與ISO-19906 和IEC6100-3 經驗公式對比曲線，其中數值模擬統計的冰載荷為單一冰厚下不同冰速(0.1 ～0.5 m/s)冰載荷統計表征量的平均值，其具體數值列于表6.三者計算的冰載荷值與冰厚的平方成線性的關系.但是，兩個規范公式和DEM 模擬計算結果存在一定的差異.ISO-19906計算的冰載荷結果最大，數值模擬計算的冰載荷最小，并且模擬得到冰載荷由于海冰冰速變化存在著差異.這主要因為DEM-FEM 耦合模型中海冰和結構的宏細觀參數選取與經驗公式存在一定差異.此外，經驗公式更多強調其使用的普遍性，因此其表述結果也相對保守.

表6 DEM-FEM計算方法與經驗公式冰載荷計算結果對比Table 6 Calculation results of ice load in the DEM-FEM simulation and the empirical formulas

3 單樁式風機結構的冰激振動響應

為保證風力渦輪機具備最佳的空氣動力特性，通常要求結構的動力響應水平被控制在一個很小的水平范圍內.因此本文重點關注風電整體結構中基礎頂端和塔筒頂端的動力響應能力和動態響應特性.

3.1 單樁式風機結構的位移響應

圖8 為冰速0.2 m/s，冰厚0.2 m 工況下塔筒頂端和基礎頂端位移響應時程曲線.從圖8 可以發現，兩者的振動位移時程曲線存在明顯的差異：塔筒頂端振動周期T=3.461 s，基礎頂點振動周期T=0.708 s.此外，由于冰載荷表現為隨機性，使其風機兩部分在每個周期的振動幅值也存在明顯差別，塔筒頂部的振動幅值范圍大于基礎頂端的振動幅值.

圖8 單樁式風機位移響應Fig.8 Displacement response of monopile-type wind turbine

圖9 風機結構振動周期Fig.9 Vibration period of wind turbine

圖9 風機結構振動周期(續)Fig.9 Vibration period of wind turbine(continued)

為進一步說明風機塔筒頂端和基礎頂端振動周期的差異性，圖9(a)為風機塔筒頂端位移振動周期與結構一階自振周期的對比；圖9(b)為風機基礎頂端位移振動周期與結構三階自振周期的對比.結果表明不同的風機部分位移振動特性呈現出差異性：塔筒頂部的振動形式主要以結構的一階自振周期振動；基礎頂端的振動則以結構的三階自振周期振動.對于安裝有抗冰錐體的單樁式風機結構，其各部分在海冰作用下表現出不同動態響應行為，其行為與結構的基礎頻率密切相關.

為定性描述塔筒頂端和基礎頂端在各個工況下的振動幅值，這里引入結構振幅極值A，表示為

式中，di為第i個位移響應.圖10 為風機不同部位結構振幅極值A與冰厚平方的變化關系.對于不同的海冰工況，當冰速一定時，風機各部分結構的振動幅度與冰厚的平方成正相關.

由上述分析可知，塔筒頂端與基礎頂端在同一工況下表現出不同的動力敏感特性，因此定義塔筒頂端--基礎頂端振幅比γ 來描述動力敏感性的差異[15]，γ 表示為

圖10 振幅極值與冰厚的關系Fig.10 Relationship between maximum amplitude and ice thickness

式中，Atower為塔筒頂端的振幅極值；Afoundation為基礎頂端的振幅極值.圖11 為各個工況下γ 取值的分布直方圖.可以看出各個工況下，塔筒頂端的振幅與基礎頂端的振幅比存在差異，但是變化范圍很小(γ 平均值為3.273，方差為0.173).因此，在平整冰與風機錐體基礎作用引起的隨機振動過程中，由海冰輸入的能量主要由風機塔筒上部輪轂、葉片所吸收.通過γ 的分布可以看出，在隨機振動模式下，海冰輸入能量在風機上下兩部分分配變化很小.

圖11 不同工況下的風機平均振幅比γFig.11 Average amplitude ratio γ under different working conditions

3.2 單樁式風機結構振動加速度響應

圖12 為冰速0.2 m/s，冰厚0.2 m下塔筒頂端和基礎頂端振動加速度響應時程曲線.風機兩部分結構振動均表現出隨機性，符合抗冰錐體結構的振動特點[4,30].

圖12 單樁式風機振動響應Fig.12 Vibration response of monopile-type wind turbine

圖13 為塔筒頂端和基礎頂端振動加速度自功率譜密度曲線.風機兩部分結構的頻譜分析結果存在很大差異：塔筒頂端的加速度功率譜主要包含3 個頻率成分，即0.293 Hz,1.465 Hz 和3.516 Hz，分別對應表2 中風機結構的一、三和五階自振基頻；基礎頂端的加速度功率譜則只存在一個頻率成分，即1.465 Hz，對應結構的三階自振頻率.此外，風機兩部分振動加速度的頻域結果顯示結構振動的主頻主要以結構的三階自振頻率成分為主.

圖13 風機振動加速度頻譜分析Fig.13 PSD of the ice-inducted vibration acceleration of wind turbine

依據單樁式導管架結構現場原型監測結果，平整冰作用下結構振動加速度響應與冰厚平方和冰速乘積呈線性關系[31]，即

式中，β 為相關系數；amax為振動加速度峰值.圖14 為塔頂頂端和基礎頂端振動加速度峰值與冰厚、冰速的關系，風機兩部分振動加速度響應峰值均與冰厚平方和冰速乘積呈線性關系，兩者的相關系數R2分別為0.935,0.854.此外，通過對比兩條擬合曲線可知，對于任一工況，風機基礎頂端的振動加速度峰值大于塔筒頂端振動加速度峰值，這說明結構各部分振動響應存在著差異，距離海冰與結構相互作用位置越近，其振動響應越劇烈，更容易引起結構疲勞損傷，在后期的風機結構的冰激疲勞分析應重點關注.

圖14 振動加速度峰值與冰速、冰厚的關系Fig.14 Relation between the ice-inducted vibration and ice velocity,ice thickness

通過對風機結構的位移周期和加速度譜分析，在海冰與風機錐體基礎結構作用過程中，風機結構的樁基部分以結構的一階自振頻率發生振動.由于海冰斷裂長度和攀爬高度的隨機性，使得海冰破壞周期沒有恒定的值，風機結構產生冰激穩態振動的條件[32]很難出現.因此，海冰破碎作為一種外部激勵使得樁基在一階自振頻率持續振動，其振動形式表現為隨機振動形式.此外，由于風機結構獨有的特點：下部為大剛度樁基和上部為高柔度塔筒，使其動力特征表現為主從式結構特性[15].這種主從式結構特性使得在復雜的冰載荷作用下，風機塔筒(子結構)和樁基(主結構)表現為不同的響應行為，風機不同部位振動周期和加速度譜兩者出現差異.

4 結論

本文基于DEM-FEM 耦合模型，采用具有粘結破碎的球形離散單元構造海冰，采用梁單元和三角形板殼單元構造海上單樁式風機結構，模擬風機與平整冰相互作用過程，研究不同冰速、冰厚工況下，風機結構冰載荷特點和結構振動響應能力和行為.研究結果表明，具有抗冰錐體的單樁式風機所承受冰載荷表現一定的隨機性與周期性，并且通過與ISO-19906 標準、IEC6100-3 規范中公式對比，驗證了本文DEM-FEM 耦合模型計算單樁式風機冰載荷正確性，并得到結構冰載荷與冰厚平方成線性關系；在分析對比不同工況下風機基礎頂端和塔筒頂端的位移響應中，通過對兩部分振動周期的統計分析發現，塔筒頂部振動周期為結構的一階自振周期，而基礎頂端的振動周期為結構的三階自振周期.此外，通過引入振幅比，定性地描述結構的動力敏感特性，即塔筒頂端的振幅極值遠大于基礎頂端振幅極值，兩者比值在各個工況中基本不變.這表明由海冰輸入的冰力激振能量大部分由風機上部渦輪機吸收耗散，并且能量占比變化很小.在分析對比不同工況下風機基礎頂端和塔筒頂端的振動加速度響應中，得到兩部分的振動加速度峰值均與冰厚平方和冰速乘積呈線性關系，并且通過對加速度進行譜分析得到風機振動各部分的振動主頻為結構的三階自振頻率.此外，由于風機結構的動力特征表現為主從式結構特性，風機不同部位的振動響應行存在很大的差異.綜上，基于DEM-FEM耦合模型可為冰區單樁式風機的冰載荷、結構動力行為評估和結構冰激疲勞分析提供有益的數值計算方法，具有很強的工程應用背景.