風力機葉片流固耦合計算的降階模型研究

孫芳錦，祝東涵，張大明

(1.廣西巖土力學與工程重點實驗室，廣西 桂林 541004；2.桂林理工大學 土木與建筑工程學院，廣西 桂林 541004；3.遼寧工程技術大學 土木工程學院，遼寧 阜新 123000；4.桂林理工大學 信息科學與工程學院，廣西 桂林 541006)

風能具有蘊藏量豐富、可再生、無污染等特性，成為新能源重要的發展方向。風力發電作為風能利用的最主要形式，近年來得到快速發展。風力機葉片是風力機的發電來源，其性能直接影響發電機的效率。目前風力機葉片柔性大，有較強的幾何非線性，在運行工況下，如陣風作用等都會引起風輪葉片發生變形和振動，而振動產生的變形又會改變空氣流場分布，產生了風荷載與風力機葉片的流固耦合作用。流固耦合作用會引起葉片的氣彈變形，可能導致葉片振動過大，或振動失穩，甚至導致疲勞破壞[1]。因此建立準確方法研究風荷載與風機葉片的流固耦合作用，并分析葉片氣彈振動的變化規律，對于為風力機葉片的合理設計、保證風力機的安全高效運行具有重要意義。

隨著計算機軟硬件技術的迅速提高，數值模擬方法已發展成為研究風荷載與風力機葉片流固耦合作用的重要工具。由于風力機葉片外形復雜、非線性強，其流固耦合計算考慮的自由度多，通常會耗費大量機時。降階模型是研究流固耦合問題的重要手段之一[2-3]，降階模型屬于經驗譜方法，假定流場包含一系列特征值和特征向量，其最重要的特性可以通過降階模型獲得。建立降階模型可以有效減小模型尺寸，節省計算機時和計算資源，因此采用降階模型進行風力機葉片的流固耦合計算，對于準確進行風力機葉片設計有著重要的意義。流固耦合計算方法可以分為強耦合法和弱耦合法，它們指的是空氣流體控制方程和結構控制方程間的耦合關系，弱耦合法是在每一時間步內先對流體控制方程和結構控制方程分別獨立求解，然后將作用在流體域模型上的氣動力荷載通過流體和結構的交界面傳遞給結構域模型，用來預測結構的位移，然后再將結構的位移作為新的荷載傳遞給流體域，此過程如此反復，直到結果收斂于指定值。強耦合方法是指在每一時間步內對流體控制方程和結構控制方程同時聯立求解，即在每一時間步內流體域和結構域互相傳遞氣動力和結構位移的過程是同時進行的。對于風力機這種大型柔性結構的流固耦合計算，弱耦合法更加實用，本文即針對弱耦合法，建立風力機葉片流固耦合計算的降階模型。

建立降階模型的第一步是進行流場分解，分解方法包括正交分解法(POD)，Vorterra級數方法和諧波平衡法等。其中正交分解法(POD)是比較常用的方法，它是從時域或頻域計算中提取結構的動力信息[4-6]，為了從試驗或數值模擬中提取流場特征向量，Sirovich[7]提出了“快照”的概念構建降階模型。目前國內外對降階模型在風力機葉片流固耦合計算中的研究還非常有限，降階代理模型曾被用于預測失速條件下風力機的非穩態氣動力荷載[8]，本征正交分解的降階模型曾用于復雜流體動力系統的低維建模[9]，文獻[10]提出了風力機葉片非定常氣動力的降階模型方法，建立了葉段振動狀態下的非定常氣動力模型,用以模擬葉片變形與氣流耦合作用下的附加非定常氣動力,給出多葉段模型的氣動彈性響應歷程，文獻[11]基于POD方法，結合邊界有限元法和降階模型對水平軸風力機的氣動特性進行了研究。

針對風力機葉片流固耦合計算復雜，耗時量大的局限性，本文采用降階模型對風力機葉片進行流固耦合計算。針對風力機葉片設計參數多，參數空間維數大的特點，采用本征正交分解(POD)方法建立流固耦合計算中空氣流體方程的降階模型，利用POD方法將流體控制方程的速度和壓力等未知量進行基模態分解并用一小部分空間模態表示，對流體方程進行時間離散，再在上述縮減空間上采用Galerkin投影方法獲得流體控制方程的最小殘差投影，構造流體方程的最小殘差降階模型。再將流體降階模型與風力機葉片進行流固耦合計算，達到減少計算時間，提高計算效率的目的。

1 基于本征正交分解的最小殘差降階模型

本文研究的是弱耦合問題，這里將采用本征正交分解(POD)建立空氣流體方程的降階模型。其主要思路是利用POD方法將流體控制方程的速度和壓力等未知量進行基模態分解，用一小部分空間模態表示，再對各基向量進行投影重構，構造最小殘差的降階模型，達到減少計算時間，提高計算效率的目的。

流體控制方程由為不可壓黏性Navie-Stokes方程(N-S方程)，即連續方程和動量方程，

(1)

(2)

流體的邊界條件為

(3)

(4)

(5)

(6)

則對于t∈[0,T]，未知量對應的基向量空間可以寫作，

(7)

在POD的基向量空間上將流體方程的未知量可以寫作，

(8)

為推導流體降階模型，需要對流體方程進行時間離散，本文采用一階Euler時間離散[13]，則有

在上述縮減空間上采用Galerkin投影方法獲得流體控制方程的最小殘差投影，進而獲得降階模型，那么對于?t∈[0,T],n=1,2,…,N，在時間步n時，流體方程的降階模型可以寫作，

其中ζn+1={an+1,bn+1}T

那么上述函數的變分形式可以寫作，

且對于時間系數，該函數應滿足，

則最小殘差應滿足如下關系，

將上述降階模型與風力機葉片有限元模型進行耦合計算，便可以計算出風力機葉片周圍風場分布、風力機葉片的氣動性能等參數。

圖1給出了本文降階模型的流程和主要步驟。

圖1 降階模型計算流程圖

2 算例分析

本文采用本文降階模型對NREL V大型風力機葉片[14]進行了流固耦合計算，NREL V風力機是一個數值風力機模型，參考了多種商用風力機原型機，美國國家可再 生能源實驗室曾對其流固耦合作用進行過試驗研究。風力機為上風向3 葉片風力機，采用變速變槳運行控制方式。風力機葉片的參數如表1所示。

表1 風力機葉片參數[14]

首先采用POD快照方法，組成流場狀態矩陣，利用奇異值分解,可得未知量對應的基向量空間。這里基于FLUENT平臺，按照上述NREL大型風力機葉片的初始條件，計算流體速度場和壓力場的POD基，得到t=0 s到20 s的速度場和相應壓力場，共200組即為快照。再應用本征正交分解(POD)，得到基向量及對應的特征值，得到的每一個特征值可以用來刻畫所對應的模態所捕捉到的能量，也就是說特征值的大小可以反應該基向量在整個系統中所占的權重。表2分別給出了速度場和壓力場的模態特征值和能量累計比例。

從表2數據可以看出，速度場的前15階和壓力場的前14階特征向量能夠捕獲系統約99%的能量，因此選擇速度場的前15階和壓力場的前14階特征向量作為降階模型的POD正交基，作為原模型的替代。

表2 速度場和壓力場的特征值和能量比

為證明本文提出的最小殘差降階模型的正確性，這里分別采用本文方法和經典Galerkin法計算了時間離散系數，這里分別給出了速度場和壓力場前4階POD模態的時間系數對比，如圖2和圖3所示。

從圖2和圖3可以看出，采用本文降階模型和經典Galerkin方法計算得到的速度場和壓力場POD模態時間系數變化趨勢一致，且比較接近，說明本文提出的最小殘差降階模型在計算中是正確的。

(a)一階模態

(a)一階模態

采用本文提出的降階模型計算了風力機葉片流固耦合計算中的受力和變形(風速為11.4 m/s)，分別如圖4和圖5所示，并將計算結果與已有文獻[15]進行了對比，以驗證本文降階模型在風力機葉片流固耦合計算中的準確性。

從圖4可以看出，采用本文降階模型得到的葉片軸向力和切向力與已有文獻[11]變化趨勢一致，且結果比較接近，可以看出，軸向力在葉片的中部靠后的位置達到最大值，切向力在葉片的中后部達到最大值。觀察圖5，可以看出，降階模型計算得到的葉片在揮舞方向和擺振方向的變形與已有文獻結果變化趨勢一致，看出葉片的主要變形區域發生在葉片的中上部分，且越接近中上部，變形越大。因此，由流固耦合作用引起的風力機葉片大變形必須在大型風力機葉片設計中引起足夠重視，并在必要時限制相應的風力機葉片變形。

(a)軸向力

(a)揮舞方向變形

同時本文還采用提出的降階模型計算了葉尖的變形情況(風速為11.4 m/s)，迎角為30°，圖6給出了葉尖隨時間變化的變形情況，圖7給出了葉尖隨方位角變化的變形情況(0°方位角是指葉片垂直向上)，并與CFD計算結果進行了對比。

觀察圖6和圖7可以看出，本文降階模型計算的葉尖隨時間和方位角變化的變形與CFD計算結果變化趨勢一致，說明了本文降階模型的正確性，以及在風力機葉片流固耦合計算中的有效性。從圖6看出，隨著時間推進，風力機葉尖揮舞運動方向的變形要大于擺振運動方向的變形，葉尖揮舞運動方向的變形隨時間增加而增大，在約30 s達到最大時，變化較小，變形基本保持不變。擺振運動方向的變形隨時間增加而減小，在約15 s達到最小值后變形基本保持不變，因此在設計時應對葉尖揮舞運動方向的變形加以重視。從圖7可以看出，隨方位角變化，風力機葉尖揮舞運動方向的變形和擺振運動方向的變形變化趨勢比較接近，基本都是在90°和270°方位角時達到最大值。

圖6 風力機葉尖變形隨時間變化對比

圖7 風力機葉尖變形隨方位角變化對比

為分析本文方法的計算效率和計算收斂情況，表3對比了不同時間步采用不同模態數(前10個模態)時計算容差收斂的情況。

從表3中可以看出，本文的降階模型方法收斂速度較快，使用較少的時間步就能達到較小容差標準；且采用較少模態就可以達到較高計算精度，可以看到，當采用6個模態就可將計算容差平均減少到10-3以下，說明本文降階模型的計算效率高，且計算精度隨時間步的增加而增加，反映了本文的降階模型具有較高精度。

表3 不同時間步時采用不同模態數時計算容差收斂情況

為了進一步說明本文降階模型的計算效率，本文對比了采用降階模型和全階模型時的相關計算參數對比，如表4所示，這里主要對比了在不同計算收斂容差時，計算葉片和葉尖變形時的平均耦合迭代次數和平均耗費計時。

表4 采用不同方法計算葉片變形的平均耦合迭代次數和計算計時

分析表4中的數據可以看出，達到同樣的計算容差時，采用本文的降階模型的平均迭代次數和平均耗時都要比全階模型少，其中平均耦合迭代次數比全階模型平均減少約42%，平均耗費機時比全階模型平均減少了約40%，主要原因在于構造最小殘差的降階模型只是用一小部分空間模態表示各參數，在縮減空間上獲得流體控制方程的最小殘差投影，因此達到了降低計算成本的目的，說明了本文降階模型在達到計算準確性的同時還提高了計算效率，是正確可靠的。

3 結 論

本文提出了基于POD方法，建立了風力機葉片的最小殘差投影降階模型，將該降階模型應用于風力機葉片的流固耦合分析中，取得了較好效果，得到的主要結論有：

(1)采用本文降階模型和經典Galerkin方法計算得到的速度場和壓力場POD模態時間系數變化趨勢一致，說明本文提出的最小殘差降階模型在計算中是正確的。

(2)本文降階模型對風力機葉片的流固耦合計算的受力和變形與已有文獻符合良好，證明了本文降階模型的準確性。同時結果表明考慮流固耦合作用時，風力機葉片的變形呈非線性分布，流固耦合作用對葉片變形有重大影響。

(3)隨著時間推進，風力機葉尖揮舞運動方向的變形要大于擺振運動方向的變形，葉尖揮舞運動方向的變形隨時間增加而增大，在約30 s達到最大。擺振運動方向的變形隨時間增加而減小。在設計時應對葉尖揮舞運動方向的變形加以重視。隨方位角變化，風力機葉尖揮舞運動方向的變形和擺振運動方向的變形基本都在90°和270°方位角時達到最大值。

(4)本文的降階模型方法收斂速度較快，使用較少的時間步就能達到較小容差，且采用較少模態就可以達到較高計算精度。且與全階模型相比，采用本文的降階模型進行流固耦合計算時，平均耦合迭代次數比全階模型平均減少約42%，平均耗費機時比全階模型平均減少了約40%，大大提高了計算效率。