孫芳錦，祝東涵，張大明

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

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

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

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

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

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

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

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

(1)

(2)

流體的邊界條件為

(3)

(4)

(5)

(6)

則對(duì)于t∈[0,T]，未知量對(duì)應(yīng)的基向量空間可以寫作，

(7)

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

(8)

為推導(dǎo)流體降階模型，需要對(duì)流體方程進(jìn)行時(shí)間離散，本文采用一階Euler時(shí)間離散[13]，則有

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

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

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

且對(duì)于時(shí)間系數(shù)，該函數(shù)應(yīng)滿足，

則最小殘差應(yīng)滿足如下關(guān)系，

將上述降階模型與風(fēng)力機(jī)葉片有限元模型進(jìn)行耦合計(jì)算，便可以計(jì)算出風(fēng)力機(jī)葉片周圍風(fēng)場(chǎng)分布、風(fēng)力機(jī)葉片的氣動(dòng)性能等參數(shù)。

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

圖1 降階模型計(jì)算流程圖

2 算例分析

本文采用本文降階模型對(duì)NREL V大型風(fēng)力機(jī)葉片[14]進(jìn)行了流固耦合計(jì)算，NREL V風(fēng)力機(jī)是一個(gè)數(shù)值風(fēng)力機(jī)模型，參考了多種商用風(fēng)力機(jī)原型機(jī)，美國(guó)國(guó)家可再 生能源實(shí)驗(yàn)室曾對(duì)其流固耦合作用進(jìn)行過(guò)試驗(yàn)研究。風(fēng)力機(jī)為上風(fēng)向3 葉片風(fēng)力機(jī)，采用變速變槳運(yùn)行控制方式。風(fēng)力機(jī)葉片的參數(shù)如表1所示。

表1 風(fēng)力機(jī)葉片參數(shù)[14]

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

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

表2 速度場(chǎng)和壓力場(chǎng)的特征值和能量比

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

從圖2和圖3可以看出，采用本文降階模型和經(jīng)典Galerkin方法計(jì)算得到的速度場(chǎng)和壓力場(chǎng)POD模態(tài)時(shí)間系數(shù)變化趨勢(shì)一致，且比較接近，說(shuō)明本文提出的最小殘差降階模型在計(jì)算中是正確的。

(a)一階模態(tài)

(a)一階模態(tài)

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

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

(a)軸向力

(a)揮舞方向變形

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

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

圖6 風(fēng)力機(jī)葉尖變形隨時(shí)間變化對(duì)比

圖7 風(fēng)力機(jī)葉尖變形隨方位角變化對(duì)比

為分析本文方法的計(jì)算效率和計(jì)算收斂情況，表3對(duì)比了不同時(shí)間步采用不同模態(tài)數(shù)(前10個(gè)模態(tài))時(shí)計(jì)算容差收斂的情況。

從表3中可以看出，本文的降階模型方法收斂速度較快，使用較少的時(shí)間步就能達(dá)到較小容差標(biāo)準(zhǔn)；且采用較少模態(tài)就可以達(dá)到較高計(jì)算精度，可以看到，當(dāng)采用6個(gè)模態(tài)就可將計(jì)算容差平均減少到10-3以下，說(shuō)明本文降階模型的計(jì)算效率高，且計(jì)算精度隨時(shí)間步的增加而增加，反映了本文的降階模型具有較高精度。

表3 不同時(shí)間步時(shí)采用不同模態(tài)數(shù)時(shí)計(jì)算容差收斂情況

為了進(jìn)一步說(shuō)明本文降階模型的計(jì)算效率，本文對(duì)比了采用降階模型和全階模型時(shí)的相關(guān)計(jì)算參數(shù)對(duì)比，如表4所示，這里主要對(duì)比了在不同計(jì)算收斂容差時(shí)，計(jì)算葉片和葉尖變形時(shí)的平均耦合迭代次數(shù)和平均耗費(fèi)計(jì)時(shí)。

表4 采用不同方法計(jì)算葉片變形的平均耦合迭代次數(shù)和計(jì)算計(jì)時(shí)

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

3 結(jié) 論

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

(1)采用本文降階模型和經(jīng)典Galerkin方法計(jì)算得到的速度場(chǎng)和壓力場(chǎng)POD模態(tài)時(shí)間系數(shù)變化趨勢(shì)一致，說(shuō)明本文提出的最小殘差降階模型在計(jì)算中是正確的。

(2)本文降階模型對(duì)風(fēng)力機(jī)葉片的流固耦合計(jì)算的受力和變形與已有文獻(xiàn)符合良好，證明了本文降階模型的準(zhǔn)確性。同時(shí)結(jié)果表明考慮流固耦合作用時(shí)，風(fēng)力機(jī)葉片的變形呈非線性分布，流固耦合作用對(duì)葉片變形有重大影響。

(3)隨著時(shí)間推進(jìn)，風(fēng)力機(jī)葉尖揮舞運(yùn)動(dòng)方向的變形要大于擺振運(yùn)動(dòng)方向的變形，葉尖揮舞運(yùn)動(dòng)方向的變形隨時(shí)間增加而增大，在約30 s達(dá)到最大。擺振運(yùn)動(dòng)方向的變形隨時(shí)間增加而減小。在設(shè)計(jì)時(shí)應(yīng)對(duì)葉尖揮舞運(yùn)動(dòng)方向的變形加以重視。隨方位角變化，風(fēng)力機(jī)葉尖揮舞運(yùn)動(dòng)方向的變形和擺振運(yùn)動(dòng)方向的變形基本都在90°和270°方位角時(shí)達(dá)到最大值。

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