汪江，陳喆，朱虹

(1.安徽省電力科學(xué)研究院，合肥市，230601;2.東南大學(xué)土木工程學(xué)院，南京市，210096)

0 引言

當(dāng)前輸電塔線設(shè)計(jì)時(shí)風(fēng)荷載一般按照靜力風(fēng)考慮，忽略脈動(dòng)風(fēng)的影響［1］。實(shí)際上塔線體系中各部分動(dòng)力耦合作用明顯［2］，忽略脈動(dòng)風(fēng)作用可能會(huì)導(dǎo)致計(jì)算結(jié)果存在較大誤差。但如采用現(xiàn)場(chǎng)實(shí)測(cè)和風(fēng)洞試驗(yàn)則所需周期長(zhǎng)、花費(fèi)高，因此利用計(jì)算機(jī)程序針對(duì)輸電塔結(jié)構(gòu)進(jìn)行脈動(dòng)風(fēng)速時(shí)程的精確模擬具有重要意義。

1 風(fēng)速時(shí)程模擬的AR法

AR法通過(guò)線性自回歸過(guò)濾器，將人工產(chǎn)生的均值為零的白噪聲隨機(jī)數(shù)，輸出為具有特定頻譜特性的隨機(jī)數(shù)序列［3］。

M個(gè)點(diǎn)空間相關(guān)脈動(dòng)風(fēng)速時(shí)程V(X，Y，Z，t)的AR模型可表示為

式中:Ψk為AR模型的自回歸系數(shù)矩陣;N(t)為M維獨(dú)立正態(tài)分布隨機(jī)過(guò)程向量。為簡(jiǎn)便起見，將V(X，Y，Z，t)簡(jiǎn)寫為 V(t)，則根據(jù)式(1)可得

對(duì)式(2)做數(shù)學(xué)期望計(jì)算，并結(jié)合自相關(guān)函數(shù)的定義和性質(zhì)，可得

式中:rV、rN分別為V(t)、N(t)的協(xié)方差矩陣。矩陣形式為

式(5)即為 AR模型的正則方程，式中，Ψ={I，Ψ1，…，Ψp}T，為(p+1)M ×M 階矩陣，I是 M 階單位矩陣;Op為pM×M階零矩陣;R為(p+1)M×(p+1)M階自相關(guān)Toeplitz矩陣

式中:rV(mΔt)是M×M階矩陣，即:

式中:m=1，…，p。

根據(jù)隨機(jī)振動(dòng)理論，功率譜密度函數(shù)與相關(guān)函數(shù)滿足Wiener-Khintchine公式，即:

式中:i，j=1，…，M。當(dāng) i≠j時(shí)，Sij(n)為空間第 i，j點(diǎn)的互功率譜密度函數(shù);當(dāng)i=j時(shí)，Sij(n)為第i點(diǎn)自功率譜密度函數(shù)。

通過(guò)式(7)、(8)求得式(6)后，代入式(5)，即可求解得到AR模型的系數(shù)矩陣Ψk和N(t)的協(xié)方差矩陣rN。

求得N(t)的協(xié)方差矩陣 rN后，N(t)按下式確定:

求出AR模型的系數(shù)矩陣Ψk和獨(dú)立隨機(jī)過(guò)程向量N(t)后，并假定初始時(shí)刻之前的風(fēng)速為0，即t≤0時(shí)，V(t)為零向量，則由式(1)可得

式中:jΔt=0，…，T;k=1，…，p。計(jì)算可以得到 M 個(gè)空間點(diǎn)的具有時(shí)間空間相關(guān)、時(shí)間步長(zhǎng)為Δt的離散脈動(dòng)風(fēng)速時(shí)程向量。

2 混合編程及模擬程序?qū)崿F(xiàn)

2.1 混合編程方法

Matlab語(yǔ)言執(zhí)行時(shí)程序逐行解釋，當(dāng)包含較多循環(huán)語(yǔ)句時(shí)存在執(zhí)行效率低的缺點(diǎn)［4］。而 Visual C++擁有強(qiáng)大的執(zhí)行效率，可與Matlab互補(bǔ)以便充分發(fā)揮各自的優(yōu)勢(shì)，適合大規(guī)模數(shù)值運(yùn)算。通過(guò)Matlab接口程序(application programming interface，API)，混合編程有多種實(shí)現(xiàn)方式［5］。

本文采用Visual C++編寫具有一定格式的程序(*.c函數(shù)的形式)，經(jīng)過(guò)編譯生成的Matlab動(dòng)態(tài)鏈接子程序(*.dll函數(shù)的形式)，即MEX文件，然后以子函數(shù)或子過(guò)程的形式提供給Matlab應(yīng)用程序，可以實(shí)現(xiàn)直接調(diào)用［6］。

2.2 程序?qū)崿F(xiàn)

Matlab強(qiáng)大的矩陣運(yùn)算功能及豐富的庫(kù)函數(shù)，可簡(jiǎn)化風(fēng)速模擬程序的編制，節(jié)省代碼并減小出錯(cuò)，提高程序開發(fā)效率。如偽隨機(jī)數(shù)的生成，通過(guò)調(diào)用rand函數(shù)，理論上其樣本空間為21492個(gè)［7］。其他如矩陣相乘、除，矩陣的Cholesky分解，快速傅里葉變換(fast fourier transform，F(xiàn)FT)運(yùn)算，風(fēng)速時(shí)程的功率譜估計(jì)及繪圖等均可以通過(guò)Matlab相應(yīng)的運(yùn)算格式和調(diào)用相應(yīng)的庫(kù)函數(shù)來(lái)實(shí)現(xiàn)，庫(kù)函數(shù)可根據(jù)問(wèn)題自動(dòng)采用最優(yōu)算法，保證了計(jì)算精度和模擬求解的穩(wěn)定性。

在編制采用AR法模擬風(fēng)速的程序時(shí)，為了保證數(shù)值積分運(yùn)算求解精度，須經(jīng)過(guò)大量循環(huán)迭代運(yùn)算。Romberg定積分算法可提高計(jì)算效率和求解精度。為測(cè)試編制好的 Visual C++語(yǔ)言 Mex源文件rombergwk.c的迭代計(jì)算是否正確，下面采用該程序計(jì)算積分

＞＞ rombergwk(1，1，1，2，10－4) ［注:前2 個(gè)入口參數(shù)對(duì)此積分不需要，故任意取為1。其中截?cái)嗾`差(計(jì)算精度)為 10－4］

與ln2=0.693147相比，該結(jié)果的精度較高，而且還可以通過(guò)最后1個(gè)入口參數(shù)(10－4)來(lái)控制精度，證明了該程序迭代計(jì)算的正確性。通過(guò)結(jié)合各類風(fēng)速譜函數(shù)編制程序時(shí)應(yīng)寫成Mex文件，在Matlab中編譯連接后方可實(shí)現(xiàn)其對(duì)Visual C++語(yǔ)言程序的直接調(diào)用。

風(fēng)速模擬程序運(yùn)行環(huán)境為Windows XP SP3系統(tǒng)，AMD Turion(tm)64 X2 TL －581.90 GHz處理器，1.00 GB內(nèi)存，Matlab 版本為7.0。采用混編技術(shù)計(jì)算用時(shí)1267.5 s，未采用混編技術(shù)計(jì)算用時(shí)2520.7 s，幾乎是采用混編技術(shù)用時(shí)的2倍。由此可看出混合編程后，運(yùn)算速度與執(zhí)行效率大為增加。

3 算例分析

自立式輸電塔屬于空間桿件體系，可采用空間桿單元和梁?jiǎn)卧⑤旊娝挠邢拊Ｐ停?］。鑒于輸電塔自重較輕，自重二階效應(yīng)不明顯，數(shù)值分析時(shí)可將輸電塔視為線性結(jié)構(gòu)進(jìn)行計(jì)算［9］。圖1為2E2-SZ2型220 kV輸電塔三維模型。

圖1 2E2-SZ2型220 kV輸電塔三維模型Fig.1 220 kV 2E2-SZ2 transmission tower model

塔架總高 58.7 m，節(jié)點(diǎn)總數(shù) 624，單元總數(shù)1602，其中迎風(fēng)面節(jié)點(diǎn)總數(shù)為161，迎風(fēng)面單元總數(shù)為325。

利用自編AR脈動(dòng)風(fēng)模擬程序，得到了考慮空間相關(guān)的多點(diǎn)水平脈動(dòng)風(fēng)模擬結(jié)果，模擬過(guò)程中，主要參數(shù)由表1給出。

表1 風(fēng)模擬程序主要參數(shù)Tab.1 Main parameters of simulation program

采用Welch法［10］得到功率譜。部分節(jié)點(diǎn)的水平風(fēng)速模擬時(shí)程曲線及功率譜如圖2所示。

對(duì)輸電塔結(jié)構(gòu)節(jié)點(diǎn)沿高度繪制功率譜后發(fā)現(xiàn)，高度位于35 m以下的節(jié)點(diǎn)，低頻部分與高頻部分模擬譜與目標(biāo)譜之間逼近較好，個(gè)別節(jié)點(diǎn)中部分模擬譜與目標(biāo)譜誤差略大;高度位于35 m以上部分高低頻部分模擬譜均與目標(biāo)譜有一定誤差，且隨著高度增加有擴(kuò)大的趨勢(shì)。事實(shí)上，若采用Bartlett［11］法繪制功率譜，低頻部分為平直型折線，高頻部分為近似以目標(biāo)譜為分界，呈上下震蕩型折線，低頻部分誤差更為明顯。因此綜合考慮為減小瞬態(tài)效應(yīng)，計(jì)算時(shí)可考慮忽略前面部分時(shí)程點(diǎn)數(shù)據(jù)。程序運(yùn)行期間快速穩(wěn)定，說(shuō)明Visual C++與Matlab混合編程切實(shí)有效，圖2中各點(diǎn)模擬功率譜與目標(biāo)譜總體吻合良好，證明所編制風(fēng)速時(shí)程程序的正確性。需要注意的是，模擬時(shí)間步長(zhǎng)不宜小于0.1 s，時(shí)間步長(zhǎng)過(guò)小時(shí)模擬精度大為降低，易產(chǎn)生錯(cuò)誤。

圖2 節(jié)點(diǎn)水平風(fēng)速模擬時(shí)程曲線及功率譜Fig.2 The simulation time-history curves of horizontal wind speeds and power spectrum

4 結(jié)語(yǔ)

(1)三維空間風(fēng)速時(shí)程的精確模擬是進(jìn)行輸電塔結(jié)構(gòu)非線性風(fēng)振響應(yīng)分析的基礎(chǔ)。本文針對(duì)輸電塔的三維空間分析模型，采用混合編程機(jī)制編制的風(fēng)速時(shí)程模擬程序具有良好的計(jì)算效率和穩(wěn)定性。

(2)不同高度處節(jié)點(diǎn)的脈動(dòng)風(fēng)速時(shí)程趨勢(shì)基本接近，但隨著高度的增加，脈動(dòng)風(fēng)速波動(dòng)范圍略有減小，表明輸電塔結(jié)構(gòu)體系的脈動(dòng)特性隨著高度的增加略有降低。

(3)由于Davenport譜紊流尺度沿高不變，對(duì)于一般構(gòu)筑物，如常見輸電塔、預(yù)應(yīng)力網(wǎng)殼結(jié)構(gòu)、橋梁等結(jié)構(gòu)均可用本風(fēng)速時(shí)程模擬程序模擬風(fēng)速，因而本程序有一定的通用性。

［1］DL/T 5092—1999P 110～500 kV架空送電線路設(shè)計(jì)技術(shù)規(guī)程［S］.北京:中國(guó)電力出版社，1999.

［2］ MomomuraY，MarukawaH，OkamuraT，etal. Full-scale measurements of wind-induced vibration of a transmission line system in a mountainous area［J］.Journal of Wind Engineering and Industrial Aerodynamics，1997，72(1):241-252.

［3］周云.結(jié)構(gòu)風(fēng)振控制的設(shè)計(jì)方法與應(yīng)用［M］.北京:科學(xué)出版社，2009.

［4］劉保柱，蘇彥華，張宏林.Matlab7.0從入門到精通［M］.北京:人民郵電出版社，2010.

［5］王素立.Matlab混合編程與工程應(yīng)用［M］.北京:清華大學(xué)出版社，2008.

［6］董振海.Matlab編譯程序和外部接口［M］.北京:國(guó)防工業(yè)出版社，2010.

［7］舒新琳，周岱.風(fēng)速時(shí)程AR模型及其快速實(shí)現(xiàn)［J］.空間結(jié)構(gòu)，2003，9(4):27-32.

［8］徐建國(guó)，石華軍，郭勇，等.大跨越輸電塔線體系有限元建模［J］.電力建設(shè)，2011，32(3):15-19.

［9］王肇民，Peil U.塔桅結(jié)構(gòu)［M］.上海:同濟(jì)大學(xué)出版社，1989:196-253.

［10］李志敏.預(yù)應(yīng)力網(wǎng)殼結(jié)構(gòu)的風(fēng)振響應(yīng)研究［D］.南京:東南大學(xué)，2007.

［11］鄒鯤，袁俊泉，龔享銥.Matlab6.x信號(hào)處理［M］.北京:清華大學(xué)出版社，2002.