基于動力學模態分解的大跨度平屋蓋風壓場研究

馮 帥，謝壯寧

(華南理工大學亞熱帶建筑科學國家重點實驗室，廣東，廣州 510641)

多次風災調查[1]結果顯示，大跨度屋蓋結構[2?3]通常具有質量較輕、柔性較大、自振頻率較低的特點，對風荷載十分敏感，其圍護結構在風荷載作用下頻頻發生破壞。大跨屋蓋結構形式多樣，風壓特性受其影響顯著，屋蓋表面脈動風壓分布可以反映建筑物在湍流邊界層下的風荷載特性[4]。一般來說，壓力場是一個隨機的、復雜的高維動態系統，很難直接理解其本質特征。而脈動風壓中隱藏的時空模式與相干結構和氣動機制密切相關，這些壓力模式對識別隨機變量場的基本物理機制和動態演化性質具有重要意義[5?6]。

為了提高流場動力學的分析效率，研究者提出了基于特征提取技術的模態分解方法—流場降階模型(ROM)，其本質是尋找一組低維的子空間(即流動模態或相干結構)，將高維、復雜非定常流場表示為子空間在低維坐標系上的疊加，從而在低維空間中描述流場演化。這種技術需要高維、大規模的流場數據作為樣本，可以直觀地展示出流場流動隨時間和空間的演化規律，對于流場的機理分析有重要意義。當前普遍采用的方法包括本征正交分解(POD)和動力學模態分解(DMD)2類。

POD技術將流場分解成若干空間正交模態，按照各個模態的能量(即特征值)大小進行排序，從而選擇出隨機場主要模態。這種方法作為一種數據壓縮和特征值萃取工具被廣泛應用于低矮建筑[7]、高層建筑[8]、大跨度建筑[9]的脈動風荷載及風致響應的研究中。盡管POD技術的應用范圍很廣，但由于協方差矩陣的存在，POD分析僅局限于變量的二階特征，并且其無法直接識別單頻動態相干結構解釋隨機系統的時間演化特征[10]。

DMD方法是一種從實驗測量或數值模擬流場中提取動力學信息的數據驅動算法，能夠用于分析流場的主要特征，或建立低階的流場動力學模型，該方法由Schmid基于動態系統Koopman分析提出[11]。其本質是將流動演化看做線性動力學過程，通過對整個過程的流場快照進行特征分析，得到表征流場信息的低階模態及其對應的特征值。DMD方法的最大特點在于分解得到的模態具有單一的頻率和增長率。另外，DMD可以直接通過各個模態的特征值表征流動演化過程，因此不需要額外建立控制方程。這種同時得到模態特征和動力學信息的特點，使DMD方法相比于POD而言，可以更好地說明模態在時間和空間上的物理意義。研究人員已經證明了DMD在解釋流場相干結構和動力學特性方面的有效性。例如，寇家慶等[12]針對DMD在流體力學研究的應用問題，綜述了DMD算法自提出以來的一系列改進以及對不同流動現象的應用，并通過典型測試算例說明DMD的應用過程，討論了DMD的研究現狀及未來發展方向。張弛等[13]利用DMD方法研究其主導脈動模態，提取出了相關模態的空間形態和脈動幅值，證明燃燒組織方式的改變會對火焰脈動的形態和規律產生影響。Luo等[6]探討了鈍體上隨機壓力場的動力特性，并將其應用于湍流邊界層中的棱柱體上的壓力場，并介紹了一種利用有限的風洞試驗數據進行氣動特性分析的算子理論方法[14]。Li等[15]通過數值模擬將DMD技術應用于雷諾數為22 000的方形棱柱繞流的原型風工程問題，以考察DMD在壓力流場分析中的準確性和實用價值。DMD方法對于模態排序沒有一個明確的規定，許多準則已經被提出用于選擇主導的DMD模態，包括α-準則[11]、I-準則[16]等。雖然DMD技術發展迅速且已被廣泛應用于許多領域[17?19]，但是在結構抗風領域的應用仍然較為少見。

本文應用POD和DMD方法對大跨度平度屋結構的隨機風壓場進行模態提取和機理分析。首先，簡述了兩種模態分解方法的理論及其不同特性；然后，通過風洞試驗獲取大跨度平屋蓋模型的隨機風壓場；最后，將兩種方法應用于風壓場中提取模態進行對比分析，并通過選取兩種方法的不同模態數對脈動風壓場進行重構比較，總結DMD和POD的獨特特征，對比兩種方法在研究屋蓋表面的隨機壓力場的內在差異。

1 POD方法和DMD方法

1.1 POD方法

任意時刻的流場xi(如速度、壓力、密度等)可以表示成平均值xˉ和脈動值x′i的疊加。POD是從一組二階統計意義上的空間數據中尋找最優正交基。要得到POD基，應求解脈動量的協方差矩陣C的特征值和特征向量。

式中，λj和uj分別為第j階模態的特征值和特征向量(也稱為POD基)。各階模態對應的模態系數aj為：

通過POD基與模態系數乘積的疊加可以表示脈動分量：

1.2 標準DMD方法及改進

通過試驗或數值模擬得到N個時刻快照，即[x1,x2,···,xN]，其中第i個時刻的快照表示為列向量xi，任意兩個快照之間的時間間隔均為Δt。基于Koopman算子理論假設[20]，流場xi+1與xi存在線性映射A，即：

式中，A為高維流場的系統矩陣，能夠反映系統的動態特征。構建快照矩陣X=[x1,x2,···,xN?1]和Y=[x2,x3,···,xN]，代入式(4)可得：

對于矩陣X，可提供一個矩陣A? 代替高維矩陣A，且這兩個矩陣相似。為尋求相似變換的正交子空間，可通過對X做奇異值分解得到。

式中：Σ為對角矩陣；U和V為酉矩陣，滿足UHU=I和VHV=I。矩陣的計算過程可視作Frobenius范數的最小化問題，結合式(5)～式(7)，可將A近似為：

為了方便分析DMD模態，定義特征值的對數映射定義為特征譜：

則第j個DMD模態對應的增長率gj和頻率fj分別為實部Re(λj)和虛部Im(λj)/2π。上述過程稱為“標準DMD方法”。

然而，只有當特征值的個數與系統的維數相符時，用標準DMD方法得到的動態系統的線性近似才是合理的。此外，標準DMD不適用于高度非線性的系統，特別是對于時間維度往往大于空間維度的試驗工況(例如風壓測點布置有限而采樣時間較長的大跨度屋蓋風洞試驗)。針對該不足之處，Le Clainche等[21]提出的高階動力學模態分解(HODMD)和Williams等[14]提出的精確動力學模態分解(EDMD)，二者的相似之處是通過引入嵌入維數(Taken's延遲嵌入定理[22])改進標準DMD方法(后文稱DMD方法)。

本文引入嵌入維數增廣大跨度平度屋隨機風壓場的試驗數據。將原始數據[x1,x2,···,xN]通過嵌入維數d進行擴展，即將延時坐標數據元素進行平移作為新的數據從而形成Hankel矩陣H：式中，hj為第j個時刻的Hankel矩陣列向量。嵌入維數d可以通過DMD模態特征譜的模是否等于或接近1確定(即DMD模態特征譜在復平面中是否位于或接近單位圓)。則時移矩陣對H1和H2分別為H1=[h1h2··· hj?1]，H2=[h2h3··· hj]。

數據以Hankel矩陣類型排列，通過將隨機風壓序列的延遲嵌入與DMD技術相結合，增強了隱藏在原始時間序列數據中的主導動力模態。最后將H1和H2代入標準DMD方法進行流場的動力分析。

對于DMD模態，本文采用Kou等[16]提出的模態排序準則，排序的模態不僅考慮了振幅大小，還包含模態的時間演變因素，綜合考慮模態對整個流場的貢獻，表示如下：

式中：Ij為第j階模態排序的參數；i為第i時刻；|| ||F為Frobenius范數；αj為第j個模態的振幅，代表了該模態對初始快照x1的貢獻，表示為：

原始快照可以近似地由DMD模態的線性組合表示：

式中：k為時間序號；xk為kΔt時刻的快照；j為第j階模態；αj?1為第j階DMD模態k時刻的模態系數，由于，特征值可以反映模態頻率是單頻，因此，模態系數的曲線是單頻衰減或增長曲線。根據流場的主導模態，通過式(14)可進行原脈動流場的重構。

1.3 POD與DMD的不同特性

DMD和POD都是具有無方程的數據驅動方法。POD根據能量從動力系統中捕捉到占優勢的空間特征，流場中的大部分能量可以用少數高階的模態來表示。由于協方差的關系，POD提取模式僅僅基于二階統計量。因此，POD缺乏捕捉精細流動結構和高階統計特征的能力。此外，POD的模態通常對應于多個頻率的耦合，而沒有實際的物理意義。并且由于不存在相鄰矩之間變量的動態關系的假設，所以POD方法不能確定所提取的模式的時間演化和穩定性。相反，DMD可以彌補POD的上述不足，對不同的流動現象具有更廣泛的適用性。DMD基于Koopman分析提供了復雜非線性動力學的線性近似，因此它可以捕捉時空相干結構特征，而不需要理解控制系統的物理方程。同時，DMD方法可以有效地提取相應頻率的相干結構模態，與POD側重于能量的方法不同，DMD提供了一種提取模式的時空演化和穩定性特征的新方法。因此，DMD在復雜動態系統的穩定性分析中具有優勢。

2 風洞試驗及數據

2.1 風洞試驗

試驗在華南理工大學5 m量級氣邊界層風洞(SCUT-1)中進行。流場地貌按照《建筑結構荷載規范》(GB 5009?2012)[23]中規定的B類地貌模擬，圖1(a)、圖1(b)分別表示風洞試驗的平均風速剖面V/Vref(Vref表示參考高度處平均風速)、湍流度剖面Iu和參考高度處的脈動風速功率譜，風洞試驗中屋蓋參考高度26.7 cm(即圖1(a)中zref)處的平均風速和湍流強度分別為9.1 m/s和0.12。采集的風壓以模型屋蓋高度處風壓作為無量綱化的參考風壓。風壓系數的表達式如下：

圖1 風洞試驗的流場參數Fig. 1 Flow field parameters of wind tunnel test

式中：Cpi(t)和pi(t)分別為測點i處的風壓系數序列和風壓序列；pa和p0分別為參考高度處的總壓和靜壓。平面屋蓋剛性測壓模型尺寸為200 cm×133.3 cm×26.7 cm(長×寬×高)，屋面共布置467個測點，試驗的幾何、風速、時間縮尺比分別為1/150、1/5、1/30。模型測點的布置原則為邊角區域加密，中間區域布置較疏，圖2為試驗模型和測點布置圖。試驗采樣頻率為300 Hz。

圖2 平屋面試驗模型和測點布置Fig. 2 Test model of flat roof and layout of tapping location

2.2 數據充分性分析

由于湍流結構和屋蓋結構之間的相互作用是非線性、混沌性的，因此在實際分析隨機風壓信號時，無法將信號中的動力特征信息充分提取。故本文通過引入嵌入維數進行動力系統的相空間重構，其基本思想是根據系統中任意分量的演化都是由與之相互作用的其他分量所決定的，相關分量的信息隱含在已知分量的發展過程中，為了重構一個等價的狀態空間，只需將已知分量在固定的時間延遲點上的測量數據作為新維處理，則可將原系統的許多性質保存，初步確定原系統的真實信息。圖3表示不同嵌入維數的DMD模態特征值分布，橫縱坐標分別表示復模態特征值的實部和虛部，顏色條表示模態排序參數Ij的數值大小。當特征點位于單位圓上或接近與單位圓周圍，則特征值是穩定或中性穩定的[24]。從圖3可以看出，特征值的數量隨著所取嵌入維數的增加而增加，因為較大的嵌入本質上擴展了動態系統的自由度和模態。由圖3(a)所示，當嵌入維數為1(即不考慮增廣數據)時，DMD特征值分散在單位圓內，表明標準DMD算法在獲取帶有噪聲的高維非線性數據的動力學方面需要更多的信息，因為理想情況下特征值應該位于或接近單位圓。從圖3(b)、圖3(c)、圖3(d)可得，隨著嵌入維數逐漸增大，特征值逐漸向單位圓靠近。當嵌入維數為8時，增廣數據矩陣維數為3736×3500，空間維度(467×8=3736)大于時間維度，絕大多數特征值落在單位圓上，嵌入過程挖掘了數據集中隱藏的模糊動態特征，滿足充分性要求。

圖3 不同嵌入維數的DMD模態特征值分布Fig. 3 Eigenvalue distribution of DMD modes with different embedding dimensions

3 試驗結果分析

3.1 風壓分布

大跨度平度屋表面的風壓可分為平均風壓和脈動風壓，圖4為0°和45°工況下的平均風壓系數pe和均方差風壓系數(反映風壓脈動特性)pe分布，從圖中發現，屋蓋表面主要以風吸力為主，每個工況的pe和pe分布輪廓相似。從圖4(a)、圖4(b)可以看出，當風垂直(0°)吹向屋蓋迎風前緣時，剪切層在迎風前緣處分離，形成破壞性的柱狀渦，故在分離區平均風吸力較大且脈動風壓波動較劇烈，pe和pe最大值出現在屋蓋角部，分別為?1.45和0.46。在柱狀渦的誘導下，隨著距離屋面迎風位置的增大，屋面風吸力先增大后逐漸減小并趨于平穩。由圖4(c)、圖4(d)所示，當來流與建筑迎風前緣成一定角度(45°)時會在屋蓋邊角區域產生錐形渦，錐形渦除具有旋渦截面的徑向和切向流速外，還具有一個沿迎風前緣的速度分量，在此作用下，旋渦內耗散的渦量得以不斷平衡和補充。因此，相比于柱狀渦，錐形渦更加穩定、持續，其作用下的風吸力也更強勁，pe和pe最大值同樣出現在屋蓋角部，分別為?3.48和0.70。

圖4 平面屋蓋風壓系數分布Fig. 4 Distribution of wind pressure coefficient on flat roof

3.2 模態分析

圖5給出了0°和45°風向角下大跨度平屋蓋脈動風壓場POD分解的前3階空間模態。POD(或DMD)模態僅表示屋蓋表面各點的相對壓力值，符號(正負號)意義不大。由圖5(a)～圖5(c)可以看出，0°風向角下POD分解的第1階模態在迎風角部區域脈動風壓較大，與圖4(b)中均方差風壓系數在角部區域分布相似，并且第2階、3階模態的極大值區域也出現在迎風前緣的角部區域，說明當0°風向角下POD分解模態的大部分脈動風壓能量包含在屋蓋的邊角區域，這與屋蓋的迎風前緣由于分離泡的分離和再附運動出現強勁風吸力的現象相吻合。從圖5(d)～圖5(f)可以得出，45°風向角下第1階和2階POD模態關于對角線互為一對反對稱模態，第3階模態關于對角線對稱，其模態分布在屋蓋對角線兩側形成兩個瓣狀，與圖4(d)中均方差風壓系數分布形狀相似，表現出明顯的錐形渦特征。

圖5 前3階POD模態Fig. 5 The first three-order POD modes

圖6為0°和45°風向角下大跨度平屋蓋脈動風壓場DMD分解的第1階～3階空間模態。本研究中DMD分解的第0階模態為頻率為0的靜態模態，反映了均勻流場特性的模態分布。為方便與POD模態進行對比，本文中取反映脈動風壓特性的前三階DMD模態進行分析。由圖6(a)～圖6(c)可以看出，在柱狀渦的作用下，0°風向角下DMD分解的3階模態能很好的捕捉到屋蓋邊角區域風壓的脈動特點，模態頻率分別為2.12 Hz、2.16 Hz和5.52 Hz，表明柱狀渦內部存在低頻運動，這將引起剪切層的拍打運動，進而形成屋蓋表面的風壓脈動[25]。增長率分別為?0.14、?0.16和?0.04，反映線性發展過程中各階模態的增長趨勢，如圖7(a)所示。前兩階模態在迎風前緣的模態分布相似，說明第2階模態是第1階模態的漂移模態，體現了在動力學發展過程中脈動風壓場隨時間的演化。由圖6(d)～圖6(f)可見，45°風向角下前三階模態均反映了錐形渦的特征，模態分布形狀與POD模態相似，模態頻率分別為0.08 Hz、0.32 Hz和0.94 Hz，第2階模態分布是由于渦軸的低頻搖擺運動使得屋蓋表面關于對角線對稱的兩點出現反相位吸力脈動。增長率分別為?1.12、?0.14和0.16，第3階模態的增長率大于0表示該模態能量在逐漸增大，這是由于錐形渦中沿迎風前緣的速度分量，在不斷平衡和補充旋渦內耗散的渦量，如圖7(b)所示。

圖6 前3階DMD模態Fig. 6 The first three-order DMD modes

圖7 前3階DMD模態系數隨時間快照的演化Fig. 7 The evolution of the first three-order DMD mode coefficients with time snapshot

通過比較發現，POD模態和DMD模態均能捕捉到屋蓋迎風前緣處破壞性旋渦的脈動特征。但從模態分布的數值結果對比可得，POD模態結果大于DMD模態結果一個數量級。這是由于兩種方法的模態時間演化機理不相同。如圖7和圖8所示，DMD方法提取的模態系數在固定頻率下具有穩定振幅或衰減的簡諧振蕩行為，能夠反映流場的時間特征，而POD的每個模態包含多個頻率的信息[26]，時間演化表現為隨機信號，這在一定程度上讓POD的脈動模態成為多個頻率段脈動的耦合。

圖8 前3階POD模態系數隨時間快照的演化Fig. 8 The evolution of the first three-order POD mode coefficients with time snapshot

3.3 脈動風壓場的重構

目前，風壓場重構的研究方法有POD方法[7]、徑向基[27]等。本文基于DMD方法對大跨度平屋蓋屋面風壓場進行重構，并與POD重構結果進行比較。

圖9所示POD和DMD模態歸一化累積能量隨模態數的變化示意圖。x軸為模態數的累積比，由累積模態與總模態的比值來定義，y軸為模態能量的累積比。由圖可得，兩種模態分解方法在45°工況達到累積模態能量的80%均比0°工況速率快，說明相比于柱狀渦，錐形渦作用下的脈動風壓波動更劇烈；而兩種方法的差別是，POD方法在0°和45°的累積模態能量達到80%所對應的模態比例分別為7.5%和5.6%，而DMD方法在0°和45°的累積模態能量達到80%所對應的模態比例分別為15.7%和13.6%。這可能與兩種模態分解方法的數據構造和模態選取方法不同，由于DMD方法引入嵌入維數構造維數更大的數據矩陣和在挑選主要模態上除了考慮振幅還考慮了時間演變的影響。直觀上在選取相同模態比例進行風壓場的重構時POD方法的表現應該優于DMD方法，但在下面的脈動風壓場重構中不盡其然，如圖10和圖11所示。

圖9 模態累積能量隨模態數的變化Fig. 9 The variation of mode cumulative energy with the number of modes

圖10和圖11分別表示0°和45°風向角下兩種模態分解方法在第600個快照(即時刻t=2 s)的脈動風壓場的重構。這里需要注意的是，DMD方法重構選取的模態數比例與POD方法相同(分別為7.5%、5.6%)，并且在重構脈動風壓場時將不包含零頻率的靜態模態。由圖可見，兩種方法的重構結果在描述原始脈動風壓場的迎風區域局部表現均較好，但DMD結果在整體輪廓上與原始風壓場更為契合，特別是在風壓場的再附和尾流區域。如前文所述，由于協方差的關系，POD提取模式僅僅基于二階統計量，造成POD缺乏捕捉精細流動結構和高階統計特征的能力。0°和45°風向角下DMD方法重構風壓場的模態頻率分別在0 Hz～22.8 Hz和0 Hz～18.9 Hz，說明低頻模態包含大部分脈動風壓能量。對比DMD和POD方法，DMD方法不僅可以從隨機風壓場中提取出脈動風壓場的主要結構，直接得到模態及對應的頻率，且可以判斷其穩定性，因此DMD方法可以同時得到流場在空間和時間上的主要特征，在揭示隨機風壓場流動機理時，DMD方法更具優勢。

圖10 0°風向角下脈動風壓場的重構Fig. 10 Reconstruction of fluctuating wind pressure field at 0°

圖11 45°風向角下脈動風壓場的重構Fig. 11 Reconstruction of fluctuating wind pressure field at 45°

4 結論

綜上所述，通過將POD和DMD模態分解方法應用在大跨度平屋蓋模型的隨機風壓場的分析中，對比兩種方法各自特征，得到如下結論：

(1) 通過將嵌入維數與DMD方法結合進行隨機風壓場動力系統的相空間重構，能夠挖掘數據集中隱藏的模糊動態特征，使分解得到的DMD模態更加中性穩定。

(2) 雖然POD與DMD算法迥異，但兩種方法分解的模態都能夠捕捉到大跨度平度屋迎風前緣處的破壞性旋渦的脈動特征。POD模態分布數值大于DMD模態的結果。這是由于DMD方法分解的是單頻模態，而POD的每個模態包含多個頻率的信息，這在一定程度上讓POD的脈動模態成為多個頻率段脈動的耦合，造成POD模態數值大于DMD模態數值。

(3) 相同比例的POD模態包含能量大于DMD模態，但當使用相同比例的模態進行脈動風壓場重構時，DMD結果比POD結果更能描述和契合原始脈動壓力場的局部特征，這是由于DMD是直接對壓力場進行重建，而POD主要是重建能量場。DMD方法分解的低頻模態包含大部分脈動風壓能量，解釋了脈動風壓場的主導頻率。因此在揭示隨機風壓場流動機理和特征，DMD方法更具優勢。