基于POD算法的數據后處理應用程序

2019-01-07 08:32:22張青山陳偉民董國祥

上海船舶運輸科學研究所學報 2018年4期

張青山，陳偉民，董國祥

(上海船舶運輸科學研究所航運技術與安全國家重點實驗室，上海 200135)

0 引言

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

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

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

1 本征正交分解技術

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

POD的基本數學思想為：假設空間上有m個離散采樣點x1,x2,…，xm，單次采樣即可同時得到這m個離散點的值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)

具體相關矩陣為

(2)

對時間相關矩陣的特征值和特征向量進行求解，有

CA=λA

(3)

式(3)中：特征值λ為非負實數；A為對應的特征向量。由此可構造特征函數

Φ=U·A

(4)

(5)

2 應用程序包的構成

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

圖1 POD數據后處理程序結構框圖

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

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

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

3個模塊的相關代碼都采用C++語言編制，模塊之間的調用和相關的數據轉換等都在程序內部實現，無需額外生成外部文件進行頻繁的輸入/輸出操作。

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

圖2 POD數據后處理主體界面

3 POD數據后處理程序的驗證

為驗證POD數據后處理程序的可靠性，以經典的自由來流圓柱繞流為例，將POD數據后處理程序應用到8 000幅PIV測量得到的瞬態速度場中。

3.1 試驗布置

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

3.2 POD流場分解

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

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

圖3 PIV測試區域基本示意

b) 第2階模態

c) 第3階模態

d) 第4階模態

圖5從圓柱尾流中提取的前4階POD模態

圖6為前4階POD模態對應的模態系數。從圖6中可看出：前2階模態系數呈現周期性變化規律，同時在時間上有一定的錯位，看似存在一定的聯系；后2階模態系數同樣有一定的周期性，這是因為其模態與卡門渦街有一定的相似之處，但周期性相對前2階模態系數較差。

a) 第1階模態系數

b) 第2階模態系數

c) 第3階模態系數

3.3 POD流場重構

由POD分解可知，圓柱尾流中包含復雜的流動，既有大尺度的卡門渦街結構，又有大量的小渦和散渦。為更加清晰地識別流場中的主導結構，基于POD數據后處理軟件對流場進行重構，重構結果見圖7。

a) 原始速度場U

c) 重構速度場U

d) 重構速度場V

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

3.4 相位平均

由上述分析可知，圓柱繞流尾流中卡門渦街呈現出周期性脫落的特征。由此確定，POD模態系數與旋渦脫落過程存在著一定的聯系。OUDHEUSDEN等[18]在研究自由來流情況下方柱分離再附流動的流動特性時發現POD分解方法中前2階模態系數與旋渦脫落過程之間存在特殊的關系，同時提出一種根據前2階模態的系數確定流場中相關結構不同時刻相位的方法。因此，在進行相位分析是需有2個模態的系數，從而確定旋渦脫落過程中各時刻的相位角。

圖8為卡門渦街脫落對應模態系數的相關分布。從圖8中可看出，卡門渦街模態系數呈圓形分布。設定(0,0)為原點，即可將整個圖像分成N等分的扇形，此時有N個相位角。本文選取8個相位角，并對流場進行相位平均。

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