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

(1.浙江華東工程咨詢有限公司，浙江 杭州 310014；2.水資源與水電工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室(武漢大學(xué))，湖北 武漢 430072；3.武漢大學(xué)水工巖石力學(xué)教育部重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430072；4.武漢大學(xué)海綿城市建設(shè)水系統(tǒng)湖北省重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430072)

0 引 言

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

本文根據(jù)海上風(fēng)機(jī)的體形特征，采用線性濾波法中的自回歸(AR)模型，利用MATLAB平臺(tái)編制了海上風(fēng)機(jī)脈動(dòng)風(fēng)速時(shí)程模擬程序，并基于江蘇黃海某風(fēng)機(jī)工程項(xiàng)目，進(jìn)行實(shí)際案例計(jì)算。同時(shí)，采用美國(guó)可再生能源實(shí)驗(yàn)室(National Renewable Resource Laboratory， NREL)開發(fā)研制的TurbSim[12]軟件進(jìn)行結(jié)果對(duì)比。

1 脈動(dòng)風(fēng)基本特征

(1)

式中，x、y、z分別代表空間中任一點(diǎn)的三維坐標(biāo)；t為某一時(shí)刻。其中平均風(fēng)認(rèn)為是周期在10 min以上的風(fēng)速幅值不變部分，并且平均風(fēng)風(fēng)速沿高程方向的幅值變化規(guī)律可用指數(shù)函數(shù)表示，即

(2)

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

表1 常用脈動(dòng)風(fēng)速功率譜

圖1 Kaimal風(fēng)速譜

2 基于AR模型的脈動(dòng)風(fēng)速時(shí)程模擬

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

2.1 AR模型基本方程

采用AR模型對(duì)空間風(fēng)場(chǎng)進(jìn)行模擬時(shí)，需要考慮風(fēng)的空間相關(guān)性，因此將式(1)擴(kuò)展至空間M點(diǎn)形成基本方程，可表示為

(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點(diǎn)坐標(biāo)，i=1，2，…，M；V(X，Y，Z，t)、V(X，Y，Z，t-kΔt)分別為空間各點(diǎn)在t時(shí)刻及(t-kΔt)時(shí)刻的脈動(dòng)風(fēng)速時(shí)程向量；p為AR模型階數(shù)，可根據(jù)AIC準(zhǔn)則(最小信息準(zhǔn)則)選取，通常取p=4；Δt為模擬風(fēng)速時(shí)程的時(shí)間步長(zhǎng)；Ψk為AR模型自回歸系數(shù)矩陣，為M×M階方陣；L為M×M階下三角矩陣；n(t)是M維均值為0、方差為1相互獨(dú)立的白噪聲向量。為簡(jiǎn)化表述，將V(X，Y，Z，t)簡(jiǎn)寫為V(t)。

根據(jù)自相關(guān)函數(shù)的性質(zhì)，對(duì)式(4)兩側(cè)同乘VT(t-jΔt)，并做數(shù)學(xué)期望計(jì)算，得到

(8)

進(jìn)一步簡(jiǎn)化為矩陣形式(正則方程)，即

(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階相關(guān)函數(shù)Toeplitz矩陣，其表達(dá)式為

R=

(12)

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

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

2.2 求解自回歸系數(shù)矩陣

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

(14)

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

(15)

式中，Cohij(f)為空間兩點(diǎn)間相干函數(shù)[22]，用以表示風(fēng)速時(shí)程的空間相關(guān)特性，其三維表達(dá)式為

Cohij(f)=

(16)

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

2.3 求解隨機(jī)過程

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

(17)

3 工程實(shí)例模擬

依據(jù)前文所述基本求解思路，本文利用MATLAB數(shù)學(xué)分析平臺(tái)，編制基于AR模型的模擬程序，以對(duì)風(fēng)機(jī)塔筒、葉輪等空間多點(diǎn)脈動(dòng)風(fēng)速時(shí)程進(jìn)行模擬。

江蘇黃海如東海域某4 MW近海風(fēng)機(jī)塔筒高度為81.25 m，葉輪直徑為146 m。計(jì)算模型簡(jiǎn)化為平面結(jié)構(gòu)，采用笛卡爾直角坐標(biāo)系，原點(diǎn)位于風(fēng)機(jī)塔筒與海平面交界處，水平方向?yàn)閄軸，Y軸以豎直向上為正。塔筒部位共模擬8點(diǎn)，葉片各模擬3點(diǎn)，各模擬點(diǎn)基本位置如圖2所示。本文目標(biāo)譜采用Kaimal風(fēng)速譜，模擬主要需要場(chǎng)地參數(shù)、頻域模擬參數(shù)以及AR模型參數(shù)等。根據(jù)項(xiàng)目資料確定場(chǎng)地類型為A類，10 m處多年平均最大風(fēng)速為8 m/s；頻域起始模擬頻率定為0.001 Hz，截?cái)囝l率為10 Hz，共劃分為1 000份；采用4階AR模型，模擬總時(shí)長(zhǎng)為100 s，時(shí)間步長(zhǎng)為0.01 s。

圖2 海上風(fēng)機(jī)簡(jiǎn)化計(jì)算模型(單位：m)

同時(shí)，本文采用TurbSim程序進(jìn)行對(duì)比分析，該軟件常用于生成海上風(fēng)機(jī)結(jié)構(gòu)的風(fēng)速入流條件，支持自定義氣象邊界和模型條件，并可輸出不同文件類型全域風(fēng)速時(shí)程數(shù)據(jù)，以輔助各類結(jié)構(gòu)分析軟件計(jì)算，例如OpenFAST、Bladed等。

限于篇幅限制，本文僅選取節(jié)點(diǎn)8、11號(hào)進(jìn)行比較。如圖3、4所示，分別為各點(diǎn)風(fēng)速時(shí)程及功率譜密度函數(shù)模擬結(jié)果對(duì)比。

圖3 風(fēng)機(jī)部分控制點(diǎn)風(fēng)速時(shí)程

圖4 風(fēng)機(jī)部分控制點(diǎn)功率譜

由圖3可以看出，本文采用AR法得到的控制點(diǎn)整體風(fēng)速波動(dòng)大于TurbSim程序模擬結(jié)果，這是因?yàn)锳R模型中無法定義湍流強(qiáng)度，而在TurbSim程序中可通過IECturbc (IEC turbulence characteristic)進(jìn)行定義。在頻域中，兩者風(fēng)速功率譜模擬結(jié)果基本一致，并且與Kaimal目標(biāo)譜擬合效果良好，達(dá)到預(yù)期目標(biāo)。

同時(shí)，需要注意的是，TurbSim程序模擬風(fēng)場(chǎng)時(shí)，需要預(yù)先設(shè)定空間網(wǎng)格，如圖5所示。對(duì)于無法通過該網(wǎng)格確定的空間點(diǎn)，例如本文中葉片上控制點(diǎn)12～17，則無法直接獲取，只能選擇進(jìn)一步提高空間點(diǎn)分布密度以獲取相關(guān)結(jié)果，這無疑提高了計(jì)算成本。而本文對(duì)空間控制點(diǎn)沒有預(yù)設(shè)要求，例如控制點(diǎn)12脈動(dòng)風(fēng)速時(shí)程如圖6所示。

圖5 TurbSim程序預(yù)定義網(wǎng)格

圖6 控制點(diǎn)12號(hào)脈動(dòng)風(fēng)速時(shí)程

4 結(jié) 論

本文基于線性濾波法中自回歸模型及Kaimal脈動(dòng)風(fēng)速功率譜，編制了相關(guān)空間域脈動(dòng)風(fēng)速時(shí)程模擬軟件，并根據(jù)江蘇黃海某工程實(shí)際確定基本參數(shù)，對(duì)海上風(fēng)機(jī)結(jié)構(gòu)脈動(dòng)風(fēng)場(chǎng)進(jìn)行模擬計(jì)算，最后與TurbSim程序相對(duì)比。計(jì)算結(jié)果表明：

(1)對(duì)海上風(fēng)機(jī)的脈動(dòng)風(fēng)速時(shí)程模擬時(shí)，采用AR模型僅需通過少量地域氣象資料便可快速生成多組空間控制點(diǎn)脈動(dòng)風(fēng)速時(shí)程，同時(shí)計(jì)算過程迅速，操作簡(jiǎn)單。

(2)本文模擬結(jié)果雖然無法考慮湍流強(qiáng)度，但其功率譜密度函數(shù)與TurbSim程序結(jié)果及目標(biāo)Kaimal功率譜譜擬合效果較好，整體模擬精確。同時(shí)，由于本文程序?qū)τ?jì)算控制點(diǎn)無任何要求，因此可廣泛適用于各類風(fēng)機(jī)結(jié)構(gòu)。

(3)以當(dāng)前自回歸模型模擬海域脈動(dòng)風(fēng)速時(shí)，本文對(duì)風(fēng)機(jī)風(fēng)場(chǎng)進(jìn)行了簡(jiǎn)化，未考慮塔影效應(yīng)等影響，用于工程實(shí)際中時(shí)，仍需要輔以相關(guān)測(cè)風(fēng)塔實(shí)測(cè)信息或氣象記錄。