基于致動線方法的5 MW海上風力機氣動彈性分析

冷 峻，郜志騰，2，*，鄭小波，2，李 曄，2

(1. 上海交通大學 船舶海洋與建筑工程學院，上海 200240；2. 上海交通大學 多功能拖曳水池實驗室，上海 200240)

0 引言

工業革命以來，煤炭、石油和天然氣等傳統燃料面臨著資源日趨耗減的問題，大規模推動可再生能源技術部署被視為主要解決辦法之一，而風能在2020年所有可再生資源中增幅最大，產能達到2750億千瓦時[1-2]。目前所使用的風能主要來自陸上。海上風電自20世紀90年代起才不斷成熟，并逐漸向大規模商業化發展[3]。然而，海上風電的引入給風能企業帶來了成本挑戰，海上風電的設計原理和相關策略尚不清楚[4]。

風能未來的發展取決于風能技術的不斷更新和成本的不斷下降[5]。為了降低風能成本，風力機產業最通用的措施之一，就是設計出尺寸更大、效率更高的風力機。通過增大葉片長度，使得掃風面積增大，可以捕獲更多能量，在其他成本保持不變的同時，可以大幅降低單位能量成本[6-7]。葉片延長后，為了降低整體重量并降低制造成本，大型風力機葉片通常采用楊氏模量和剛度較低的復合材料，造成葉片柔性顯著增加，同時由于葉片幾何形狀也較為復雜，因此葉片氣動載荷和結構動力學響應之間存在強烈的氣動彈性耦合。同時，大撓度改變了有效掃風面積，進而導致功率輸出、力矩和推力也受到影響。此外，由于非定常氣動載荷與風力機結構之間的相互作用，容易發生葉片振動和顫振，嚴重時甚至可能影響風力機的使用壽命。

氣動彈性預測是葉片設計中的重要步驟，非定常氣動載荷和葉片變形之間復雜的相互作用對風力機性能和安全分析至關重要。相比于單獨的結構或氣動設計方法，結構動力學和空氣動力學的聯合分析可提供更全面和系統的數據，這將有助于提高葉片設計性能，同時保障整個風力機系統的安全可靠性，從而使風能發電更具競爭力。在風力機的結構動力學和空氣動力學研究過程中，考慮到模擬復雜的真實海洋狀態的困難性以及實驗的高成本[8-9]，數值模擬應用最為廣泛。

氣動彈性的數值分析已逐漸成為風電行業的重要研究內容，一些大型風能研究組織，如美國可再生能源國家實驗室（NREL）、Ris?DTU等，已經開發了多種氣動彈性研究工具[10-11]。這些工具大多采用葉素動量理論來模擬空氣動力學，而結構動力學建模則采用剛體多體動力學、模態法和傳統的線性有限元法。此外，來自其他學科領域（如航空航天、機械工程、土木工程）的許多研究人員也提出了不同的模型用于解決風力機氣動彈性問題。Wang等[12]基于ANSYS-FLUENT軟件和殼單元，建立了風力機單向耦合模型，將CFD建模計算出的氣動載荷作為邊界條件映射到有限元模型，而葉片撓度未映射回氣動模型。Li等[13]對5 MW浮式風力機，將結構梁單元和精細CFD模型耦合，并采用顯式風湍流模型進行求解。Mo等[14]將一維梁單元和葉素動量理論耦合，模擬了5 MW風力機葉輪，結果表明葉片振動和變形對氣動載荷有顯著影響，與穩態氣動模型相比，動態失速會導致葉片氣動載荷的劇烈波動，從而對疲勞設計產生很大影響。Meng等[15]結合致動線模型和Euler-Bernoulli梁理論，提出了一種彈性致動線模型，對風力機進行氣動彈性模擬，分析葉片振動速度對攻角的影響。Ma等[16]討論了旋轉剛化效應、應力加強效應、風力機葉片的非線性彎曲對風力機變形的影響，研究結果表明風力機的變形將影響遠尾流區的速度分布。Li等[17]結合變分漸近梁截面分析方法和Euler-Bernoulli梁模型，分析了彎扭耦合特性對風力機葉片動力學性能的影響。Qian等[18]分別采用線性模態疊加法和非線性幾何精確梁方法，計算了百米量級葉片在穩態來流（風）和湍流來流（風）條件下的動態響應，研究表明幾何非線性效應對于大型風力機氣動彈性分析不可忽略。

總的來說，學者們對風力機的氣動彈性耦合研究已經進行了許多嘗試，但仍然存在一些不容忽視的問題和困難。首先，三維葉片建模的復雜性和巨大的計算量仍然是一大挑戰。風力機數值模擬的目標之一就是減少建模和計算時間，同時提供足夠準確的預測，但傳統的結構動力分析工具無法滿足全面新設計的需求。同時，旋轉運動引起的大位移和復合材料彈性耦合引起的幾何非線性是不可忽略的，特別是對于大載荷下的柔性結構而言。另外，大多數現有的氣動彈性模型是基于小葉片撓度假設的傳統線性模型，但這些模型不再適用于經常經歷大撓度變形的長柔性葉片設計。

因此，針對以上大型風力機發展中出現的問題，本文基于柔性多體動力學、致動線方法和大渦模擬，建立了一種新型的雙向流固耦合模型，在避免了三維葉片建模的復雜性和巨大的計算量的同時，能夠準確預測結構載荷及動力學響應，評估風力機的功率和負載，并對風力機尾跡和周圍流場進行更詳細的模擬和分析。

1 計算方法

1.1 致動線方法

目前，計算流體力學模擬風力機及其周圍流場的方法主要有葉素動量法、致動盤法、致動線法和三維高保真模擬等。其中，致動線法憑借著計算量小、可以更好地解析葉片周圍流動細節等優勢，在風力機領域得到廣泛應用。致動線模型可用于精確再現風力機葉片的渦流系統，并適用于簡單的結構化網格，大量的研究證明了致動線法是預測風力機葉片非定常氣動載荷的有效方法。

在致動線方法中，葉片被簡化成一維的線，并被劃分成有限個致動點，如圖1所示。在獲取每個致動點上的風速Urel和攻角αrel后，通過查找相應的翼型數據表，獲得對應的升力系數CL和阻力系數CD，并計算每個致動點的升力L和 阻力D[19]：

圖1 致動線高斯投影原理Fig. 1 Projection of Gaussian function applied on the actuator line

其中，ρ為流場的密度，c為弦長，w為致動點所在截面寬度。

然后，根據高斯投影函數將升阻力分配到致動點周圍相應球形區域內流體網格中心的體積力f上：

其中，ε為高斯函數的投影半徑，r為網格中心到投影致動點的距離，F為升阻力的合力。

最后，將結果添加到Navier-Stokes方程中進行求解：

其中，ui為i方 向的流體速度，t為 時間，ρ為流體密度，xi為 流體網格坐標，p為 壓強，fi為對應方向的體積力。

1.2 柔性多體動力學方程

柔性多體系統由多個彈性和剛性組件組成，其中柔性體的變形不可以忽略。系統的動力學方程如下[20]：

其中，x(t)是 位移向量，M、C、K分別是質量矩陣、阻尼矩陣和剛度矩陣，x˙ 是x關 于時間的導數，F是外部力向量。

為了方便求解，將動力學方程等價為如下一階微分方程組：

其中，β為系統動量，?為完整約束，ψ為非完整約束，λ?和 λψ分別是對應約束的拉格朗日乘子，下標/x為 對x求導。

對方程組化簡并省略高階項，得到如下方程組：

其中，hb0為 常系數，I為單位矩陣，R為殘差向量。

在線性化后，采用隱式線性多步積分法對動力學方程進行迭代求解，詳細求解過程可參閱文獻[20]。此外，在可再生能源領域，相比于傳統梁模型，幾何精確梁模型[21-22]在分析細長柔性結構的幾何非線性方面具有優越性[23]，被逐漸應用到結構研究中[24-25]，本文將葉片簡化為幾何精確梁模型，每個梁單元具有3個節點，每個節點具有6個自由度，可以用于模擬葉片的大范圍旋轉運動和幾何非線性變形。

1.3 耦合機制

本文采用分區式耦合方法，將流體和結構求解器分別表示為算子F和S，求解器的輸入和輸出值分別對應結構運動學變量s如位移和速度，以及動力學變量f如升阻力，變量之間的關系如下：

為了準確及時地交換和映射葉片結構和流場之間的耦合信息，雙方求解器通過本地套接字在每個耦合時間步長交換數據。如圖2所示，在每個耦合時間步長中，結構求解器向流體求解器發送預測的運動學數據，即位移和速度，并從流體求解器接收相關的力數據，通過線性插值，投影到結構單元上。同時，流體求解器接收映射過來的運動學數據，通過致動線方法求解升力和阻力，并將這兩種體積力投影到每個致動點周圍的流體網格上。本文數值計算基于柔性多體動力學程序MBDyn[20]和計算流體力學軟件OpenFOAM[26]實現。

圖2 流固耦合計算流程圖Fig. 2 Flow chart of the fluid-structure interaction simulation

2 數值模型

2.1 NREL-5 MW風力機參數

本文選取的風力機模型為NREL提出的5 MW基準風力機，該風力機符合全球風力機產業尺度趨勢，被廣泛應用于數值模擬領域[27]。其參數如表1所示。

表1 NREL-5 MW基準風力機的主要參數Table 1 Key properties of the NREL-5 MW baseline wind turbine

2.2 數值模型

本文所有計算均在笛卡爾坐標系下進行，計算域為15D（長）×10D（寬）×10D（深），其中D為轉子直徑。總體布局如圖3所示。順流線方向，風力機上游為5D，風力機下游為10D，計算域足夠長、利于尾流充分發展；風力機輪轂和橫向邊界之間的距離為5D，可以認為橫向邊界對風力機幾乎沒有影響。來流入口邊界條件設為速度11.4 m/s的定常均勻來流，出口邊界條件為速度梯度為0。底部、頂部和橫向邊界條件設置為滑移邊界。如圖3所示，六面體網格用于離散計算區域，風力機附近的區域在三個方向上依次加密，最大網格尺寸為20 m，最小網格尺寸為1.25 m，網格總數為6891528。

圖3 流場計算域網格劃分Fig. 3 Refinement of the computational mesh

為了使預測出的體積力沿葉片分布足夠光滑，每個葉片被劃分為90個致動點。為了避免數值振蕩，并確保網格可以投影上體積力，ε應不小于最小網格尺寸的2倍[28]，因此選擇 ε=3。湍流模型采用標準的Smagorinsky模型，并使用大渦模擬方法模擬風力機尾跡。考慮到使用笛卡爾網格的致動線模型時，每個時間步長內葉尖通過的距離不應超過一倍網格長度，因此本算例流體求解器的時間步長設為0.005 s。在結構模型中，每個葉片被簡化為15個幾何精確梁單元，時間步長為0.001 s。算例共采用112個計算節點，總計算時長為93 h，計算所得數據均基于初始狀態60圈后達到穩定狀態的20圈內的平均值。

3 結果分析

3.1 模型驗證

為了驗證結構模型在靜力學和瞬態分析中的準確性，將其與商業軟件ANSYS和開源幾何精確梁求解器GEBT進行了比較。驗證模型為長60 m的懸臂梁，在梁的自由端施加100 kN的恒定荷載。選擇參考文獻[27]中的單個葉片橫截面材料特性：葉片截面分布質量為294.734 kg/m；翼型截面弦長方向剛度為1102.38 × 106N·m2；寬度方向截面剛度為3447.14 ×106N·m2；扭轉剛度為144.47 × 106N·m2；軸向拉伸剛度為1632.70×106N。梁在本文模型、GEBT和ANSYS中均劃分為10個單元，時間步長為0.0001 s。靜態和瞬態葉尖位移分別如圖4（a）和圖4（b）所示。可以看出，本文采用的結構模型所計算的結果與GEBT和ANSYS的結果吻合良好。

圖4 不同方法計算的葉尖位移Fig. 4 Tip displacement calculated with different methods

此外，將額定風速下5 MW風力機的本文氣動彈性計算結果與NREL官方報告中FAST程序的計算結果[27]進行了比較，結果列于表2。從數據可以看出：葉片結構變形、功率和推力均與FAST結果吻合較好；由于葉片的變形和振動，轉子的掃風面積減小，轉子功率和推力有所降低。

表2 與FAST的結果比較Table 2 Result comparison with FAST

3.2 尾跡分析

除了驗證和計算風輪的變形和基本性能，本文也分析了葉片結構響應以及氣動彈性對風力機尾跡的影響。圖5給出了風力機三個葉片在流向上的變形沿徑向分布情況，可以看出，葉尖處變形最大。圖6給出了葉尖在流向上的振動響應。風力機的三個葉片在5.35～5.65 m范圍內振動，其誘導的三個葉尖渦脫落，形成螺旋形尾跡，向下游移動，并在1D距離后略微破裂，如圖7所示。

圖5 葉片在流向上的變形沿徑向分布情況Fig. 5 Radial distribution of the blade deflection along the streamwise direction

圖6 葉尖沿流向振動的時間歷程Fig. 6 Time history of the tip vibration along the streamwise direction

圖8 給出了剛性葉片和柔性葉片算例對應的Q沿流向剖面的云圖。相比于圖8（a）剛性葉片模擬結果，圖8（b）中考慮葉片柔性后，葉片的變形影響了致動線的位置和體積力的分布，進而推遲了葉尖渦生成的原始位置。

圖8 Q沿流向切片圖Fig. 8 Q slice of the rotor along the streamwise direction

圖9 進一步分別提取了剛性葉片和柔性葉片模擬的葉尖渦位置，其中z用轉子半徑R參照。可以看出，葉片對葉尖渦的影響可持續到下游2D距離，并且在相同的軸向位置處，變形對葉尖渦渦強度的影響總是大于對根部渦附近的影響。

圖9 剛性葉片和柔性葉片的葉尖渦位置比較Fig. 9 Tip vortex locations compared between the rigid blade and the flexible blade

4 結論

本文基于柔性多體動力學、致動線方法和大渦模擬，建立了一種用于風力機氣動彈性和尾流分析的新型流固雙向耦合數值模型。與傳統致動線模型相比，本文氣動彈性模型基于幾何精確梁理論，進行了葉片自由度捕捉和重新建模，取代了傳統的致動線模型單一剛體旋轉運動描述，并保留了建模簡單便捷的優點。

經過驗證，本文流固雙向耦合模型在轉子功率、推力、結構變形等參數的計算上與其他方法所得結果相吻合。除此之外，本文模型通過和大渦模擬結合取代了葉素動量法，可以準確捕捉風力機的尾跡結構，包括葉尖渦和葉根渦。研究表明，葉片在風載荷作用下有顯著的振動響應，同時，葉片的變形影響了致動點和體積力的分布，最終影響風力機的性能和尾跡特性。

綜上所述，葉片的柔性在風力機氣動彈性設計中不可忽略。本文方法可以通過簡單的建模，避免三維建模的復雜性，并準確預測風力機性能和尾流特性，與傳統方法相比，更適用于現代兆瓦級復合材料風力機流固耦合模擬和尾流分析。

由于致動線法在來流速度和投影函數的選取上有一定局限性，需要在未來的工作中進行改善，使得葉片體積力分布更貼近真實葉片形狀。更具體的結構參數響應研究以及幾何非線性對葉片氣動的影響研究，也將在接下來的工作中展開。