基于AR模型的海上風機脈動風速時程模擬

駱光杰，朱洪澤，郭 健，蘇 凱，3，4

(1.浙江華東工程咨詢有限公司，浙江 杭州 310014；2.水資源與水電工程科學國家重點實驗室(武漢大學)，湖北 武漢 430072；3.武漢大學水工巖石力學教育部重點實驗室，湖北 武漢 430072；4.武漢大學海綿城市建設水系統湖北省重點實驗室，湖北 武漢 430072)

0 引 言

近年來風電新能源行業不斷發展，近海風機由于其具備占地少、噪音小，風資源豐富等優勢[1]成為當前風電能源發展的首選。對于兆瓦級風機而言，由于結構高聳且柔性較大，同時，海上運行條件復雜多變，尤其以風荷載變化最為明顯，對結構安全提出了嚴峻考驗。無論采用葉素-動量定理[2]或簡化伯努利方程[3]計算風機荷載，還是通過指數律風輪廓線計算塔筒風荷載[4-5]，均需要首先獲取結構風速變化情況。重大工程建設中，常采用風洞試驗或現場實測等手段確定風荷載，但前者耗資巨大，并且結果不具普適性[6]，后者存在布設困難，并難以獲取整體風場分布情況?；跀抵的M方法，可利用較少數據快速獲取大量代表性風速時程曲線，在風機結構設計和強度校核過程中，具備突出的作用和意義。相關數值模擬方法主要有線性濾波法、諧波疊加法、逆傅里葉變換法、小波分析法等[7]，不同方法存在各自優缺點，但線性濾波法中自回歸模型(Auto-regressive model， AR)擁有計算量小和速度快的特點，近年來在相關領域研究中使用較多。溫鵬等[8]通過自回歸算法對陸上風機進行風速模擬，同時利用模型殘差平方和進行模型階數判定；劉文洋等[9]針對大跨度空間結構，設計MATLAB程序以實現風速記錄模擬與整理存儲；張文福等[10]針對細化的標量過程AR模型和向量過程AR模型，從計算精度及效率角度出發進行了討論分析；舒新玲等[11]通過VC和MATLAB混合編程實現空間風速場的模擬。

本文根據海上風機的體形特征，采用線性濾波法中的自回歸(AR)模型，利用MATLAB平臺編制了海上風機脈動風速時程模擬程序，并基于江蘇黃海某風機工程項目，進行實際案例計算。同時，采用美國可再生能源實驗室(National Renewable Resource Laboratory， NREL)開發研制的TurbSim[12]軟件進行結果對比。

1 脈動風基本特征

(1)

式中，x、y、z分別代表空間中任一點的三維坐標；t為某一時刻。其中平均風認為是周期在10 min以上的風速幅值不變部分，并且平均風風速沿高程方向的幅值變化規律可用指數函數表示，即

(2)

而脈動風是由風的不規則運動產生，在工程中常假定為平穩隨機過程，即在任意時刻的風速概率分布完全相同的隨機過程。因此常通過功率譜函數進行時域上的脈動風速時程模擬。目前，根據對強風觀測記錄處理方式的不同。存在兩類常用于海上風機結構設計的脈動風速功率譜(power spectral density， PSD)，包括Davenport譜[15]、Kaimal譜[16]、Von-Karman譜等[17-20]。常見功率譜如表1所示，其中Kaimal風速譜由于計算簡單，并且考慮了風速譜隨時間變化情況，美國、歐洲規范均采用該風速譜用于風荷載模擬等問題研究[4]，本文即采用該風速譜進行計算，風速譜模型示意如圖1所示。

表1 常用脈動風速功率譜

圖1 Kaimal風速譜

2 基于AR模型的脈動風速時程模擬

AR模型是將均值為0的白噪聲隨機序列通過濾波器，使其輸出具有指定譜特征的隨機過程[21]。由于不涉及其他變量，整體計算相對簡單，非常有利于脈動風速時程的模擬。

2.1 AR模型基本方程

采用AR模型對空間風場進行模擬時，需要考慮風的空間相關性，因此將式(1)擴展至空間M點形成基本方程，可表示為

(3)

其中，

X=[x1，x2，…，xM]T

(4)

Y=[y1，y2，…，yM]T

(5)

Z=[z1，z2，…，zM]T

(6)

N(t)=L·n(t)

(7)

式中，(xi，yi，zi)為空間第i點坐標，i=1，2，…，M；V(X，Y，Z，t)、V(X，Y，Z，t-kΔt)分別為空間各點在t時刻及(t-kΔt)時刻的脈動風速時程向量；p為AR模型階數，可根據AIC準則(最小信息準則)選取，通常取p=4；Δt為模擬風速時程的時間步長；Ψk為AR模型自回歸系數矩陣，為M×M階方陣；L為M×M階下三角矩陣；n(t)是M維均值為0、方差為1相互獨立的白噪聲向量。為簡化表述，將V(X，Y，Z，t)簡寫為V(t)。

根據自相關函數的性質，對式(4)兩側同乘VT(t-jΔt)，并做數學期望計算，得到

(8)

進一步簡化為矩陣形式(正則方程)，即

(9)

d=L·LT

(10)

Ψ=[I，Ψ1，…，ΨP]T

(11)

式中，Ψ為(p+1)M×M階矩陣；I為M階單位陣；O為PM×M階零矩陣；R為(p+1)M×(p+1)M階相關函數Toeplitz矩陣，其表達式為

R=

(12)

需要特別注意的是，R矩陣中任意元素R(kΔt)均是M階方陣，k=0，…，p，表示為

式中，矩陣中各元素Rij(kΔt)代表時間相差kΔt時，空間任意兩點i與j間的相關函數，其中i，j=1，…，M。

2.2 求解自回歸系數矩陣

對于自相關矩陣R，可由Wiener-Khintchine公式可知式(13)矩陣中各元素Rij(kΔt)可由功率譜函數表示為

(14)

式中，當i=j時，功率譜函數Sij(f)即為自功率譜密度函數；而i≠j時，則稱為互功率譜密度函數，并且可由兩點處自功率譜函數求得，計算公式為

(15)

式中，Cohij(f)為空間兩點間相干函數[22]，用以表示風速時程的空間相關特性，其三維表達式為

Cohij(f)=

(16)

通過選擇表1中脈動風速譜模型，可計算空間M點各點自功率譜密度函數，再由式(14)～(16)可得到相關系數矩陣R，求解正則方程可得到AR模型自回歸系數Ψ和矩陣d。

2.3 求解隨機過程

對求解得到的d矩陣進行Cholesky分解，可得下三角矩陣L，進而通過隨機白噪聲向量n(t)得到獨立隨機過程向量N(t)。通常為方便計算，假定初始時刻及之前時刻各點風速均為0，即：t≤0時，V(t)=0。在此基礎上，將自回歸系數Ψ和矩陣d代入基本方程(1)并可表示為

(17)

3 工程實例模擬

依據前文所述基本求解思路，本文利用MATLAB數學分析平臺，編制基于AR模型的模擬程序，以對風機塔筒、葉輪等空間多點脈動風速時程進行模擬。

江蘇黃海如東海域某4 MW近海風機塔筒高度為81.25 m，葉輪直徑為146 m。計算模型簡化為平面結構，采用笛卡爾直角坐標系，原點位于風機塔筒與海平面交界處，水平方向為X軸，Y軸以豎直向上為正。塔筒部位共模擬8點，葉片各模擬3點，各模擬點基本位置如圖2所示。本文目標譜采用Kaimal風速譜，模擬主要需要場地參數、頻域模擬參數以及AR模型參數等。根據項目資料確定場地類型為A類，10 m處多年平均最大風速為8 m/s；頻域起始模擬頻率定為0.001 Hz，截斷頻率為10 Hz，共劃分為1 000份；采用4階AR模型，模擬總時長為100 s，時間步長為0.01 s。

圖2 海上風機簡化計算模型(單位：m)

同時，本文采用TurbSim程序進行對比分析，該軟件常用于生成海上風機結構的風速入流條件，支持自定義氣象邊界和模型條件，并可輸出不同文件類型全域風速時程數據，以輔助各類結構分析軟件計算，例如OpenFAST、Bladed等。

限于篇幅限制，本文僅選取節點8、11號進行比較。如圖3、4所示，分別為各點風速時程及功率譜密度函數模擬結果對比。

圖3 風機部分控制點風速時程

圖4 風機部分控制點功率譜

由圖3可以看出，本文采用AR法得到的控制點整體風速波動大于TurbSim程序模擬結果，這是因為AR模型中無法定義湍流強度，而在TurbSim程序中可通過IECturbc (IEC turbulence characteristic)進行定義。在頻域中，兩者風速功率譜模擬結果基本一致，并且與Kaimal目標譜擬合效果良好，達到預期目標。

同時，需要注意的是，TurbSim程序模擬風場時，需要預先設定空間網格，如圖5所示。對于無法通過該網格確定的空間點，例如本文中葉片上控制點12～17，則無法直接獲取，只能選擇進一步提高空間點分布密度以獲取相關結果，這無疑提高了計算成本。而本文對空間控制點沒有預設要求，例如控制點12脈動風速時程如圖6所示。

圖5 TurbSim程序預定義網格

圖6 控制點12號脈動風速時程

4 結 論

本文基于線性濾波法中自回歸模型及Kaimal脈動風速功率譜，編制了相關空間域脈動風速時程模擬軟件，并根據江蘇黃海某工程實際確定基本參數，對海上風機結構脈動風場進行模擬計算，最后與TurbSim程序相對比。計算結果表明：

(1)對海上風機的脈動風速時程模擬時，采用AR模型僅需通過少量地域氣象資料便可快速生成多組空間控制點脈動風速時程，同時計算過程迅速，操作簡單。

(2)本文模擬結果雖然無法考慮湍流強度，但其功率譜密度函數與TurbSim程序結果及目標Kaimal功率譜譜擬合效果較好，整體模擬精確。同時，由于本文程序對計算控制點無任何要求，因此可廣泛適用于各類風機結構。

(3)以當前自回歸模型模擬海域脈動風速時，本文對風機風場進行了簡化，未考慮塔影效應等影響，用于工程實際中時，仍需要輔以相關測風塔實測信息或氣象記錄。