劉瑞莉, 徐公勇, 宋淳宸, 邱順冬
(1. 西南交通大學,四川成都 610031; 2. 成都興城人居地產投資集團有限公司, 四川成都 610031; 3.建研科技股份有限公司,北京 100000)
隨著計算機技術的快速發展,蒙特卡洛模擬在工程實際上得到越來越廣泛的應用。與頻域方法相比,蒙特卡洛模擬方法在考慮結構非線性、系統隨機性以及其它相關問題等方面具有明顯的優勢。例如,其被廣泛應用于考慮結構和氣動力非線性的大跨橋梁抖振響應分析[1]。此外,蒙特卡洛模擬方法也常被用作檢驗頻域方法計算精度的標準。
目前,蒙特卡洛模擬方法主要有譜表示方法和線性濾波法。譜表示方法由于具有嚴密的理論推導以及容易使用等特點而得到廣泛的應用,而快速傅里葉變換(FFT)算法的引入使得該類方法的模擬效率得到顯著提高[2-4]。最初,多點平穩隨機過程的譜表示方法為單指標方法,模擬出的樣本不具備各態歷經性[5]。隨后,Deodatis[6]將其拓展到了多指標的譜表示方法,以實現各態歷經性的要求。為了進一步提高模擬效率,研究者從減少功率譜矩陣分解時間方面提出了多種改進的譜表示方法[7-9]。近期,Peng等[10]提出了一種基于隨機場的風場模擬方法,該方法不僅可避免功率譜矩陣分解,同時可利用二維FFT極大提高模擬效率。
線性濾波法又稱白噪聲濾波法,即將均值為零的高斯白噪聲通過設計的濾波器,從而獲得具有目標統計特性的隨機過程。該類方法包含自回歸(Autoregressive, AR)模型、移動平均 (Moving averages, MA)模型以及自回歸和移動平均模型(Autoregressive and moving-average, ARMA)模型[11-13]。與ARMA模型相比,AR模型的參數相對較少,因此更加便于使用,其已被廣泛應用于多點風速時程的模擬。Novak等[14]對基于AR模型的風速模擬誤差進行了詳細分析。舒新玲等[15]采用AR模型模擬了大跨屋蓋的隨機脈動風速時程,并提出VC與Matlab混合編程模擬的快速實現方法。李春祥等[16]基于AR模型模擬了一幢高度為200 m超高層建筑的風速時程,并驗證了其可行性。李明等[17]基于西堠門大橋風場模擬的算例,對AR模型的定階進行了討論,并認為階數過高與過低均將影響模擬精度,選擇合適的模型階數即可得到滿意的精度又可節省計算時間。
基于上述研究基礎,本文對譜表示方法與基于AR模型的風場模擬方法進行對比研究。首先,分別對譜表示方法和基于AR模型的模擬方法進行了介紹。隨后,以某橋梁風場模擬為算例,對兩種方法在模擬精度和效率方面進行了對比。最后,給出了相應的結論。
考慮一個均值為零的n變量平穩隨機過程X(t)=[X1(t),X2(t),...,Xn(t)]T,設其功率譜密度矩陣可表示如下:
(1)
式中:ω為頻率;Sjk(ω),j,k=1, 2, …,n;j≠k為互功率譜密度函數,一般為復數;當j=k時,Sjj(ω)為自功率譜密度函數。功率譜矩陣元素具有如下性質:
Sjj(ω)=Sjj(-ω),j=1, 2,…,n
(2)
j,k=1, 2, …,n;j≠k
(3)
其中*表示復數共軛。由于功率譜矩陣S(ω)為非負定的Hermitian矩陣,因此其可進一步分解如下:
S(ω)=H(ω)HT*(ω)
(4)
式中:H(ω)為下三角矩陣,其元素Hjk(ω)可表示為如下的極坐標形式:
Hjk(ω)=|Hjk(ω)|eiθjk(ω)
(5)
(6)
式中:Im和Re分別表示復數的虛部和實部;上述多點隨機過程被分解為n個完全不相干的子多點隨機過程的疊加,且每個子過程的元素之間完全相干。因此,多點隨機過程X(t)可采用下式模擬[6]:

j=1, 2,…,n
(7)
式中:Δω=ωu/N;ωu為截止頻率;N為頻率離散點數;ωml=(l-m/n)Δω,m=1, 2, ...,n;φml為在[0,2π]上服從均勻分布的隨機相位角。值得說明,當ωml退化為ωl=lΔω時,上述經典的雙指標模擬公式即退化為單指標模擬公式。
設上節的n變量平穩隨機過程X(t)可表示為如下的多點AR模型形式[9]:
(8)
式中p為模型階數;Φr為n×n的自回歸系數矩陣;L為n×n的系數矩陣;Y(t)=[Y1(t),Y2(t),...,Yn(t)]T是均值為0的高斯白噪聲隨機過程,且滿足下式:
E[Y(t)YT(s)]=Inδ(s-t)
(9)
式中δ(·)為狄拉克δ函數;In為n×n的單位矩陣。
將式兩邊同乘于X(t-qΔt),然后取數學期望可得:
q=1,2,…,p
(10)
式中RXX(-qΔt)=E[X(t)XT(t-qΔt)]。上式可進一步表示為如下的矩陣形式:
(11)

(12)
由式(11),可求得不同時刻的自回歸系數矩陣:
(13)
將式(8)兩邊同乘于YT(t)和XT(t),并求期望后可分別得到:
RXY(0)=L
(14)
(15)
將式(14)代入式(15),可得如下的表達式:
(16)
設L為下三角矩陣,其可以通過對等式右邊部分進行Cholesky分解得到。在求得Φr和L矩陣后,則可以通過該AR模型進行多點風速隨機過程的模擬。
為了對比譜表示方法和基于AR模型的模擬方法在模擬精度和效率上的差異,對泰州長江大橋主梁上的順風向脈動風速進行了模擬。該橋為連接泰州市和揚中市的三塔懸索橋,其跨徑布置為390 m、1 080 m、1 080 m及390 m,如圖1所示。設模擬點的間距為7.2 m,則共需要模擬301個點。

圖1 泰州長江大橋示意
設主梁模擬點的自功率譜相等,且均用如下的雙邊Kaimal譜描述:
(17)

此外,采用如下的Davenport相干函數描述不同模擬點之間的相干性:
(18)
式中:λ=7為衰減系數;ξ表示兩點之間的距離。
為了對比,譜表示方法考慮雙指標和單指標兩種工況,基于多點AR模型的模擬方法選取較為常用的4階和6階模型兩種工況[16]。
在精度對比時,截止頻率均取為4π。譜表示方法的頻率離散點數為1 024,時間步長為0.25 s。在基于AR模型的模擬方法中,頻率離散點數為2 048,時間步長為0.2 s。限于篇幅,本文僅給出基于4階AR模型模擬方法得出的風速時程樣本,如圖2所示。可以看出,模擬點1與2之間的風速時程接近,而模擬點76則與前兩者差別較大,其原因是模擬點76與點1之間的距離已超500 m,相關性較小。
圖3給出了不同計算工況下模擬時程的估計功率譜與理論功率譜的對比。可以看出,雙指標和單指標譜表示方法的模擬精度均高于AR模型的模擬精度。此外,雙指標譜表示方法的精度高于單指標譜表示方法,4階AR模型的模擬精度則高于6階AR模型。圖4給出了不同工況下的自/互相關函數對比圖。可以看出,雙指標譜表示方法的精度同樣明顯高于其它模擬方法。而單指標譜表示方法在相關函數的模擬精度上與基于AR模型的模擬方法差別不大。
在進行模擬效率對比時,所有工況的截止頻率和頻率離散點數均分別取為4π和1 024,模擬點數分別取32、64、128、256以及512五種工況。本算例的所有工況均在一臺處理器為Intel (R) Core (TM) i7-7500U,主頻為2.70GHz,內存為8G的筆記本電腦上進行,計算結果如表1所示。可以看出,當模擬點數為64時,基于4階AR模型的模擬方法效率為單指標譜表示方法的2倍,而為雙指標譜表示方法的80倍。因此,基于AR模型的模擬方法效率明顯高于譜表示方法,尤其是與雙指標譜表示方法相比。此外,4階和6階AR模型的模擬效率相差較小。在所需內存上,當模擬點數分別為512和128時,單指標和雙指標譜表示方法即出現內存不足的現象。由此可見,譜表示方法所需要的內存明顯大于基于AR模型的模擬方法。

(a)模擬點1

(b)模擬點2

(c)模擬點76

(a) 單指標譜表示方法

(b) 雙指標譜表示方法

(c) 4階AR模型

(d) 6階AR模型
本文首先對譜表示方法與基于AR模型的風場模擬方法分別進行了介紹,然后通過數值算例詳細對比了單指標譜表示方法、多指標譜表示方法、基于4階AR模型以及6階AR模型的風場模擬方法在模擬精度和效率上的差異,得出了以下結論:
(1)雙指標譜表示模擬方法的模擬精度明顯高于單指標譜表示方法以及基于AR模型的模擬方法。

(a)模擬點1

(b) 模擬點2

(c)模擬點1與2

(d) 模擬點1與76

模擬點數單指標譜表示雙指標譜表示4階AR模型6階AR模型 322.4924.020.881.07646.06237.812.953.2112822.41/11.9616.08256112.82/68.1080.49512//458.46491.73
備注:表中/表示電腦內存不足。
(2)單指標譜表示方法在功率譜的模擬精度上高于基于AR模型的模擬方法,而在相關函數模擬精度上,兩者之間差別不明顯。
(3)基于AR模型的模擬方法在模擬效率以及消耗內存上的性能均明顯優于譜表示方法。因此,當風場模擬點數很大時,可考慮使用基于AR模型的模擬方法。