基于POD算法的數(shù)據(jù)后處理應(yīng)用程序

張青山， 陳偉民， 董國祥

(上海船舶運(yùn)輸科學(xué)研究所 航運(yùn)技術(shù)與安全國家重點(diǎn)實(shí)驗(yàn)室， 上海 200135)

0 引 言

尾流、射流和分離流等流體在流動(dòng)過程中表現(xiàn)出的主要特征是具有不同空間尺度和時(shí)間尺度的相干結(jié)構(gòu)。流體在流動(dòng)過程中多主體結(jié)構(gòu)的共同存在使得識(shí)別每個(gè)相干結(jié)構(gòu)的時(shí)空特征具有一定的難度，這種識(shí)別在研究傳熱、傳質(zhì)和流動(dòng)噪聲等物理過程方面具有重要作用。近年來，大渦模擬(Large Eddy Simulation, LES)和直接數(shù)值模擬(Direct Numerical Simulation, DNS)的廣泛應(yīng)用及時(shí)間分辨粒子圖像測速儀(Time Resolved-Particle Image Velocimetry, TR-PIV)的迅速發(fā)展使得以合理的空間分辨率和時(shí)間分辨率來精確捕獲瞬時(shí)變化的流量成為可能，進(jìn)而使得從大量的數(shù)值和試驗(yàn)數(shù)據(jù)中提取顯著的相干結(jié)構(gòu)更加顯著、精確。

SIROVICH[1]和LUMLEY[2]的研究證實(shí)本征正交分解(Proper Orthogonal Decomposition, POD)在近幾十年被廣泛用于從流場中提取并確定相干結(jié)構(gòu)。LUMLEY[2]、NOACK等[3]和LIU等[4]的研究結(jié)果表明，POD算法通過對采用快照序列計(jì)算得到的時(shí)間和空間相關(guān)矩陣進(jìn)行對角化，能提取最有能量的結(jié)構(gòu)，按照其自身的能量和相關(guān)結(jié)構(gòu)的層次進(jìn)行排序。MEYER等[5]使用立體PIV獲得的橫流中的湍流和射流的瞬時(shí)速度場進(jìn)行POD分析，結(jié)果表明尾渦是主要的動(dòng)力學(xué)結(jié)構(gòu)，與噴流核心的強(qiáng)相互作用有關(guān)，而噴流剪切層旋渦是局部的，且不占優(yōu)勢。PROTHIN等[6]根據(jù)基于時(shí)間分辨立體粒子圖像測速儀(TR-SPIV)獲得的NACA0015翼型周圍的流動(dòng)，采用基于POD系數(shù)的相位平均來強(qiáng)調(diào)尾流的不穩(wěn)定特性。DERI等[7]通過層析PIV測量得到的流場的三維POD重建，量化圍繞薄平板埋藏在高雷諾數(shù)湍流中的三維相干渦旋結(jié)構(gòu)。文獻(xiàn)[8]～文獻(xiàn)[11]基于PIV測量得到的平板繞流的二維流場，利用POD分解技術(shù)提取平板兩側(cè)及下游存在的大尺度相干結(jié)構(gòu)，發(fā)現(xiàn)模態(tài)系數(shù)之間有一定的聯(lián)系，可利用模型系數(shù)對流場進(jìn)行相位平均分析。這些研究雖然能實(shí)現(xiàn)POD算法的功能，但需針對不同的數(shù)據(jù)結(jié)果作相應(yīng)的調(diào)整，對于實(shí)際工程應(yīng)用而言操作繁瑣，不能保證便捷、迅速地實(shí)現(xiàn)。

本文基于POD基礎(chǔ)理論開發(fā)POD的應(yīng)用程序，設(shè)計(jì)輸入模塊、POD核心模塊和輸出模塊，并采用C++語言編制核心代碼；此外，為方便用戶使用，采用QT設(shè)計(jì)用戶界面，對功能模塊和用戶界面進(jìn)行打包并形成應(yīng)用程序包。

1 本征正交分解技術(shù)

POD作為一種處理大批量數(shù)據(jù)的方法，廣泛應(yīng)用于天氣預(yù)報(bào)、圖像識(shí)別、信號分析、數(shù)據(jù)壓縮、隨機(jī)過程和海洋學(xué)等領(lǐng)域中。POD分解可從系統(tǒng)已知的試驗(yàn)解或數(shù)值解中得到一系列基函數(shù)(這些基函數(shù)常被稱作“模態(tài)”)，并保證這些基函數(shù)在最小二乘意義下是最優(yōu)的正交基函數(shù)。在湍流速度場研究中，POD分解能將相干結(jié)構(gòu)及其包含的能量對應(yīng)起來，從而使POD能用來識(shí)別湍流流場中各種能級的結(jié)構(gòu)。LUMLEY[12-14]將POD引入到湍流研究領(lǐng)域中，通過對空間速度相關(guān)函數(shù)進(jìn)行正交分解來識(shí)別流場中的相干結(jié)構(gòu)。SIROVICH[1]通過對HOLMES等[15]提出的直接POD方法進(jìn)行改進(jìn)，提出快照方法(Snapshot POD)。該方法利用時(shí)空轉(zhuǎn)換原理解決空間點(diǎn)數(shù)太多造成空間矩陣龐大而難以直接用POD求解的問題，使采用POD方法處理復(fù)雜流場成為可能。

POD的基本數(shù)學(xué)思想為：假設(shè)空間上有m個(gè)離散采樣點(diǎn)x1,x2,…，xm，單次采樣即可同時(shí)得到這m個(gè)離散點(diǎn)的值u(x1),u(x2),…，u(xm)；共采集N次，得到空間場(u1(x1),u1(x2),…，u1(xm))，(u2(x1),u2(x2),…，u2(xm)),…, (uN(x1),uN(x2),…，uN(xm))，以矩陣的形式表示為

(1)

具體相關(guān)矩陣為

(2)

對時(shí)間相關(guān)矩陣的特征值和特征向量進(jìn)行求解，有

CA=λA

(3)

式(3)中：特征值λ為非負(fù)實(shí)數(shù)；A為對應(yīng)的特征向量。由此可構(gòu)造特征函數(shù)

Φ=U·A

(4)

(5)

2 應(yīng)用程序包的構(gòu)成

POD數(shù)據(jù)后處理應(yīng)用程序包由數(shù)據(jù)輸入模塊、POD數(shù)據(jù)處理核心模塊和數(shù)據(jù)輸出模塊等3個(gè)功能模塊組成(見圖1)。

圖1 POD數(shù)據(jù)后處理程序結(jié)構(gòu)框圖

1) 數(shù)據(jù)輸入模塊的主要功能是從外部磁盤中讀取POD數(shù)據(jù)后處理所需的原始數(shù)據(jù)，包括文件列表管理(獲取指定目錄及子目錄中的相關(guān)信息)、文件格式篩選(根據(jù)用戶指定的格式提取相應(yīng)格式的文件)、文件結(jié)構(gòu)解析及讀取(解析單個(gè)文件中的數(shù)據(jù)結(jié)構(gòu)并讀取數(shù)據(jù))和數(shù)據(jù)轉(zhuǎn)換(將從文件中讀取的數(shù)據(jù)轉(zhuǎn)為POD數(shù)據(jù)處理核心模塊所需的類型)等4個(gè)子模塊。

2) POD數(shù)據(jù)處理核心模塊的主要功能是對數(shù)據(jù)進(jìn)行POD分析，包括POD分解(進(jìn)行POD分解計(jì)算，分別計(jì)算特征值、特征模態(tài)和模態(tài)系數(shù))、POD重構(gòu)(基于POD分解得到的特征值、模態(tài)和模態(tài)系數(shù)對流場進(jìn)行重構(gòu))和相位平均(基于重構(gòu)的流場進(jìn)行相位平均，同時(shí)計(jì)算流場瞬態(tài)相位角和相位流場信息等)等3個(gè)子模塊[4, 8, 10-11]。

3) 數(shù)據(jù)輸出模塊的主要功能是將POD分析得到的結(jié)構(gòu)導(dǎo)入/保存到外部磁盤中，包含結(jié)果列表管理(根據(jù)用戶選擇的POD功能模塊生成相應(yīng)的文件列表)、輸出格式管理(根據(jù)用戶指定的格式篩選相應(yīng)的數(shù)據(jù)轉(zhuǎn)換方式)、數(shù)據(jù)轉(zhuǎn)換(根據(jù)用戶指定的格式對POD分析結(jié)果進(jìn)行轉(zhuǎn)換)和輸出文件解析及保存(解析用戶指定的格式，同時(shí)將轉(zhuǎn)換后得到的POD分析結(jié)果保存到外部硬盤中)等4個(gè)子模塊。

3個(gè)模塊的相關(guān)代碼都采用C++語言編制，模塊之間的調(diào)用和相關(guān)的數(shù)據(jù)轉(zhuǎn)換等都在程序內(nèi)部實(shí)現(xiàn)，無需額外生成外部文件進(jìn)行頻繁的輸入/輸出操作。

為方便使用，采用QT語言設(shè)計(jì)用戶交界面，避免用戶誤輸出或誤操作導(dǎo)致程序崩潰，基本的界面設(shè)計(jì)見圖2。同時(shí)，為使用戶熟悉程序的使用，在界面中添加提示子模塊，在輸入?yún)?shù)時(shí)提示用戶，確保輸入的參數(shù)正確和程序運(yùn)行正常。此外，設(shè)計(jì)多線程管理功能，保證程序運(yùn)行時(shí)主界面不卡頓。

圖2 POD數(shù)據(jù)后處理主體界面

3 POD數(shù)據(jù)后處理程序的驗(yàn)證

為驗(yàn)證POD數(shù)據(jù)后處理程序的可靠性，以經(jīng)典的自由來流圓柱繞流為例，將POD數(shù)據(jù)后處理程序應(yīng)用到8 000幅PIV測量得到的瞬態(tài)速度場中。

3.1 試驗(yàn)布置

PIV試驗(yàn)在一個(gè)低速循環(huán)水槽[17]中進(jìn)行，試驗(yàn)段截面尺寸為100 mm×100 mm。直徑D=8 mm的圓柱體直接橫跨在試驗(yàn)段中，安裝在半深的高度上。圖3為PIV測試區(qū)域基本示意。自由來流速度為1 m/s。圓柱體的下游采用激光片光源照亮，采用高速相機(jī)拍攝，拍攝區(qū)域的大小為10D×8D。經(jīng)過PIV互相關(guān)算法計(jì)算之后即可得到速度場。

3.2 POD流場分解

在界面中設(shè)定好原始數(shù)據(jù)所在目錄、原始場數(shù)量、原始場類型、輸入格式和輸出格式等參數(shù)之后，開始進(jìn)行POD分解計(jì)算，得到特征值、特征模態(tài)和模態(tài)系數(shù)等參數(shù)。

圖4為前100階POD模態(tài)對應(yīng)的特征值，這些特征值代表對應(yīng)模態(tài)所包含的能量。從圖4中可看出，隨著模態(tài)階數(shù)的增加，特征值逐漸減小。前4階POD模態(tài)的含能較高，因此提取前4階POD模態(tài)并繪制速度云圖(見圖5)。從圖5中可看出，前2階POD模態(tài)為卡門渦街脫落，第3階和第4階POD模態(tài)與卡門渦街脫落有關(guān)。由此可知：特征值較大的模態(tài)對應(yīng)的模態(tài)結(jié)構(gòu)的空間尺度較大；特征值較小的模態(tài)對應(yīng)的模態(tài)結(jié)構(gòu)的空間尺度較小。因此，只需關(guān)注大尺度的流動(dòng)結(jié)構(gòu)即可。

圖3 PIV測試區(qū)域基本示意

b) 第2階模態(tài)

c) 第3階模態(tài)

d) 第4階模態(tài)

圖5從圓柱尾流中提取的前4階POD模態(tài)

圖6為前4階POD模態(tài)對應(yīng)的模態(tài)系數(shù)。從圖6中可看出：前2階模態(tài)系數(shù)呈現(xiàn)周期性變化規(guī)律，同時(shí)在時(shí)間上有一定的錯(cuò)位，看似存在一定的聯(lián)系；后2階模態(tài)系數(shù)同樣有一定的周期性，這是因?yàn)槠淠B(tài)與卡門渦街有一定的相似之處，但周期性相對前2階模態(tài)系數(shù)較差。

a) 第1階模態(tài)系數(shù)

b) 第2階模態(tài)系數(shù)

c) 第3階模態(tài)系數(shù)

3.3 POD流場重構(gòu)

由POD分解可知，圓柱尾流中包含復(fù)雜的流動(dòng)，既有大尺度的卡門渦街結(jié)構(gòu)，又有大量的小渦和散渦。為更加清晰地識(shí)別流場中的主導(dǎo)結(jié)構(gòu)，基于POD數(shù)據(jù)后處理軟件對流場進(jìn)行重構(gòu)，重構(gòu)結(jié)果見圖7。

a) 原始速度場U

c) 重構(gòu)速度場U

d) 重構(gòu)速度場V

通過比較發(fā)現(xiàn)：POD重構(gòu)保留了原始速度場的整體特性，圓柱繞流下游的卡門渦街脫落現(xiàn)象基本上與原始速度場一致；同時(shí)，重構(gòu)速度場剔除了尺度較小的旋渦結(jié)構(gòu)，使得大尺度相關(guān)結(jié)構(gòu)的空間結(jié)構(gòu)更加光滑、突顯。

3.4 相位平均

由上述分析可知，圓柱繞流尾流中卡門渦街呈現(xiàn)出周期性脫落的特征。由此確定，POD模態(tài)系數(shù)與旋渦脫落過程存在著一定的聯(lián)系。OUDHEUSDEN等[18]在研究自由來流情況下方柱分離再附流動(dòng)的流動(dòng)特性時(shí)發(fā)現(xiàn)POD分解方法中前2階模態(tài)系數(shù)與旋渦脫落過程之間存在特殊的關(guān)系，同時(shí)提出一種根據(jù)前2階模態(tài)的系數(shù)確定流場中相關(guān)結(jié)構(gòu)不同時(shí)刻相位的方法。因此，在進(jìn)行相位分析是需有2個(gè)模態(tài)的系數(shù)，從而確定旋渦脫落過程中各時(shí)刻的相位角。

圖8為卡門渦街脫落對應(yīng)模態(tài)系數(shù)的相關(guān)分布。從圖8中可看出，卡門渦街模態(tài)系數(shù)呈圓形分布。設(shè)定(0,0)為原點(diǎn)，即可將整個(gè)圖像分成N等分的扇形，此時(shí)有N個(gè)相位角。本文選取8個(gè)相位角，并對流場進(jìn)行相位平均。

圖9為相位平均過程中參與每個(gè)相位角平均的流場個(gè)數(shù)。從圖9中可看出，在8 000個(gè)瞬態(tài)場中，參與這8個(gè)相位角平均的流場的數(shù)量基本接近，為80～100個(gè)。

圖8 卡門渦街脫落對應(yīng)模態(tài)系數(shù)的相關(guān)分布

圖9 相位平均過程中參與每個(gè)相位角平均的流場個(gè)數(shù)

圖10為重構(gòu)流場的法向速度場，給出了一個(gè)卡門渦街脫落周期內(nèi)相同時(shí)間間隔的8個(gè)相位角的速度場，同時(shí)移除了平均場。從圖10中可清晰地看到圓柱繞流尾流中卡門渦街交替脫落的特性，在向下游輸運(yùn)的過程中先逐漸增長，然后逐漸耗散。

a) T/8

b) 2T/8

c) 3T/8

d) 4T/8

e) 5T/8

f) 6T/8

g) 7T/8

h) 8T/8

圖10重構(gòu)流場的法向速度場(移除時(shí)均平均速度場)

4 結(jié) 語

本文基于本征正交分解理論編制POD數(shù)據(jù)后處理程序，同時(shí)以經(jīng)典的圓柱繞流為例驗(yàn)證POD數(shù)據(jù)后處理程序的各項(xiàng)功能，得到以下結(jié)論：

1) POD數(shù)據(jù)后處理程序運(yùn)行穩(wěn)定，在POD分析計(jì)算過程中不卡頓；

2) 程序的各子模塊設(shè)定基本合理，能獨(dú)立完成相應(yīng)的功能，如輸入模塊只負(fù)責(zé)原始數(shù)據(jù)載入，POD核心模塊負(fù)責(zé)POD分解、POD重構(gòu)和相位平均，輸出模塊負(fù)責(zé)將POD分析結(jié)構(gòu)導(dǎo)出到外部硬盤中；

3) POD數(shù)據(jù)后處理程序能很好地實(shí)現(xiàn)POD數(shù)據(jù)分解、POD數(shù)據(jù)重構(gòu)和相位平均等功能；

4) 基于POD數(shù)據(jù)后處理程序，可獲取更加豐富的流場信息，是研究流場的非定常特征的新方法。