基于非侵入式混沌多項式法的隨機阻曳流CFD模擬不確定度量化

夏 立， 鄒早建,b， 袁 帥， 曾智華

(上海交通大學 a. 船舶海洋與建筑工程學院； b. 海洋工程國家重點實驗室， 上海 200240)

隨著計算機科學技術的飛速發(fā)展，計算流體動力學(Computational Fluid Dynamics，CFD)技術在工業(yè)生產(chǎn)和科學研究中得到了廣泛的應用.而CFD的計算質(zhì)量和可信度一直備受關注.目前船舶CFD不確定度分析研究主要集中在驗證與確認(Verification and Validation, V&V).驗證(Verification)是對數(shù)值誤差進行分析，以保證數(shù)值方法的準確性；而確認(Validation)主要將數(shù)值計算結果和實驗數(shù)據(jù)進行對比，以評估CFD數(shù)值模擬和真實物理模型的相似程度.在V&V分析中，CFD模型和數(shù)值方法是確定的，流動參數(shù)、幾何模型和邊界條件等因素都是確定量，獲得的結果也是確定的量.然而在真實的物理環(huán)境中，許多物理量往往是不確定的，比如流動參數(shù)、幾何模型和邊界條件.這些不確定量的累積，可能會對系統(tǒng)的響應造成很大的影響.

實際上，CFD不確定度量化主要包含兩方面的內(nèi)容：其一是確定隨機變量的分布情況；其二是量化由隨機變量引起的不確定性對模型輸出參數(shù)的影響.本文重點關注第二方面，即假定隨機變量的分布情況是已知的，研究不確定性如何在流場中傳播.

不確定度量化的傳統(tǒng)方法可分為兩種：一是非統(tǒng)計方法，如局部敏感度分析、矩量法和攝動法，可以解決不確定性較小的問題或線性系統(tǒng)問題；二是統(tǒng)計方法，如蒙特卡洛(Monte Carlo，MC)法等.MC法雖可求解非線性系統(tǒng)的不確定度量化問題，但其本質(zhì)上是基于多次重復獨立的數(shù)值模擬.使用MC法進行不確定度量化研究，需要對n維空間中的隨機參數(shù)進行采樣，這些采樣點即為計算的工況，最后根據(jù)計算結果進行統(tǒng)計分析，得到響應量的統(tǒng)計特性.若要使響應量的期望概率水平達到10-n，至少需要生成10n+2個樣本點[1].故該方法計算效率低、成本高，在工程實際問題中，往往被當作最后的求解手段.

近年來，混沌多項式(Polynomial Chaos，PC)方法被用于不確定度量化研究.該方法源于Wiener[2]提出的齊次混沌方法，而后由Ghanem等[3]發(fā)展到固體力學領域.近年來，該方法也被應用到CFD領域[4-18].在船舶計算水動力學方面，不確定度量化研究相對很少.He等[19]使用非侵入式混沌多項式方法，將船體形狀、航速、波浪分別作為隨機的輸入變量，就Delft Catamaran船的水動力性能進行了不確定度量化分析.Diez等[20]使用代理模型、正交法和Karhunen-Loève展開法就Delft Catamaran船的阻力和下蹲進行了不確定度量化分析.Stern等[21]就船舶水動力學的不確定度量化問題及其驗證方法進行了綜述.

本文比較了驗證與確認方法和不確定度量化的區(qū)別，分析了不確定度量化的研究現(xiàn)狀；介紹了不確定度量化問題的定義以及用于解決該問題的非侵入式混沌多項式方法，并以隨機二維阻曳流為例分析和討論了計算結果.

1 數(shù)學模型

根據(jù)是否需要對求解器進行修改，PC方法分為侵入式混沌多項式(Intrusive Polynomial Chaos, IPC)法和非侵入式混沌多項式(Non-Intrusive Polynomial Chaos, NIPC)法.其中IPC法需要對求解器進行大量修改，其具體步驟為將控制方程的解展開成混沌多項式的形式，并將其代入控制方程、邊界條件、初始條件等，通過改寫后的方程求取混沌多項式的待定系數(shù).IPC法改寫求解器的工作量巨大，目前只應用于簡單的隨機CFD不確定度量化問題.NIPC法不需要改變CFD求解器，只需將其當作一個黑箱，隨機輸入變量和響應量之間的關系是通過滿足一定條件的正交基函數(shù)的線性組合來表示的，即通過混沌多項式展開式來表示.因此，當混沌多項式的待定系數(shù)確定后，就可以獲得響應量的全部概率分布信息，大大提高了計算效率.近年來，NIPC法由于其高效性，逐漸被應用于不確定度量化研究.本文采用NIPC法，以隨機二維阻曳流為例進行不確定度量化分析.

1.1 問題定義

不確定度量化問題可以定義為已知隨機變量ξ(輸入變量)的分布，確定響應Y(輸出變量)的分布情況，如下式：

Y=f(x,t,ξ)

(1)

式中：x為空間變量；t為時間變量；ξ為滿足一定分布的隨機變量組成的向量.

根據(jù)混沌多項式的理論，若函數(shù)空間是以內(nèi)積定義為L2范數(shù)的Hilbert空間，且空間中的輸入變量ξ是獨立的，則響應Y可以分成確定的和隨機的兩部分：

(2)

式中：ci為待定的系數(shù)，只和確定的空間變量及時間變量有關；d為隨機變量的個數(shù)；Φi(i=1,2,…)為混沌多項式，只和隨機變量ξ有關.

為了求解實際問題，通常需要將式(2)進行截斷.若將混沌多項式展開到p階，則截斷后的多項式形式如下：

(3)

式中：混沌多項式的系數(shù)共有Nc個，

(4)

1.2 混沌多項式

若隨機變量服從高斯分布，則最優(yōu)的多項式是Hermite多項式.Xiu和Karniadakis發(fā)展了基于Wiener-Askey方法的混沌多項式求解方法[4]，應用該方法可以求出多種服從典型分布的隨機變量對應的混沌多項式，如表1所示.

表1 不同分布對應的混沌多項式

1.3 非侵入式混沌多項式方法

NIPC法可分為基于正交的方法和基于回歸的方法兩種.基于正交的方法主要是利用混沌多項式的正交性，求取混沌多項式的系數(shù).但這種方法需要積分，而數(shù)值積分方法通常采用高斯積分法，所以采樣點一般選為高斯點.然而，該方法會隨著不確定量的維數(shù)增加而產(chǎn)生大量的采樣點，大大增加工作量，即發(fā)生所謂的維數(shù)災難.

本文采用的是基于回歸的NIPC法.該方法首先需要生成滿足分布的特定采樣點，然后在這些點上進行CFD計算，得到響應Y后，通過混沌多項式的系數(shù)求出所關心的量的分布情況.

超定方程可以表示成下式：

(5)

寫成矩陣與向量的形式為

Φc=y

(6)

待定系數(shù)可以通過最小二乘法求得：

c=(ΦTΦ)-1ΦTy

(7)

與基于正交的方法不同，該方法需要很好地配置采樣點，而不是選取高斯積分時的高斯點.采樣點的個數(shù)通常取2倍的混沌多項式的待定系數(shù)的個數(shù)[14]，即2Nc.

在求取混沌多項式的待定系數(shù)后，響應Y的均值與方差可以通過以下兩式得到：

2 采樣方法

本文采用隨機采樣(Random Sampling，RS)法和拉丁超立方采樣(Latin Hypercube Sampling，LHS)法進行采樣.這兩種算法均能生成符合特定分布的隨機數(shù).LHS法本質(zhì)上是一種受約束的分層隨機采樣方法，針對有眾多隨機變量的問題特別有效.LHS法在滿足隨機變量的分布情況的同時，可以大大減少采樣點.從有著d維不確定量的系統(tǒng)中抽取N個樣本的具體步驟如下：

(1) 將每一維度分成等概率的N個區(qū)間(例如，均勻分布時區(qū)間的長度相同)；

(2) 在上述各個區(qū)間中隨機抽取一個點，且滿足每個與軸垂直的超平面最多含有一個樣本.

可以看出，樣本總數(shù)N并不會隨著隨機變量的維數(shù)而增加.對于高維問題，LHS法可以大大減少采樣點，從而減小計算量.

3 二維隨機阻曳流算例

3.1 求解實例

如圖1所示，以下平板的某一點為原點建立坐標系，其中x方向與上平板的拖曳方向一致，y方向垂直于x軸向上. 兩平行平板間為牛頓流體，其黏度為μ=1 kPa·s.上下壁面均滿足壁面無滑移假設，即流體速度與壁面速度一致.上平板拖曳速度為u=0.01 m/s，v=0；下平板為固面，u=0，v=0.本文將左右兩邊的入口壓力p1和出口壓力p2處理為不確定量，且滿足均勻分布U(950 Pa, 1 050 Pa).采用速度-壓力有限元方法[22]對算例進行求解.

圖1 二維阻曳流模型

3.2 控制方程

該問題實際為求解在進出口壓力和拖曳復合作用下定常的二維流場分布.其連續(xù)性方程可以表示為

(10)

x方向運動方程：

(11)

y方向運動方程：

(12)

(13)

式中：τ為切應力，滿足牛頓流體本構方程.

3.3 網(wǎng)格收斂性分析

本文使用四邊形單元對計算域進行離散；速度單元階次比壓力單元高一次，以便提高計算精度.為確定合適的網(wǎng)格數(shù)量，首先進行網(wǎng)格收斂性分析，采用GCI(Grid Convergence Index)方法[23]進行網(wǎng)格收斂性驗證.網(wǎng)格收斂性分析必須在確定條件下進行.本文網(wǎng)格收斂性分析采用的工況為入口壓力p1=1 kPa，出口壓力p2=0.共采用了粗、中、細3套網(wǎng)格，網(wǎng)格數(shù)分別為9×9、14×14以及19×19.根據(jù)GCI方法，獲得細網(wǎng)格、中等密度網(wǎng)格和粗網(wǎng)格對應的數(shù)值計算結果φ1、φ2、φ3，據(jù)此可得到近似相對誤差、外推相對誤差和細網(wǎng)格收斂指數(shù).需要指出的是，當ε21=φ2-φ1和ε32=φ3-φ2都接近于0時，可以認為已經(jīng)得到了精確解.本文以出口邊界中點處的速度值作為網(wǎng)格收斂性分析的標準，求得的φ1、φ2、φ3分別為 10.624 9、10.625 0和 10.625 0 mm/s.可以認為這3套網(wǎng)格都得到了精確解.最終采用的網(wǎng)格離散情況如圖2所示.

圖2 計算域網(wǎng)格離散

4 計算結果與討論

4.1 計算結果

本文分別采用MC法和基于回歸的NIPC法對二維隨機阻曳流進行了不確定度量化分析.其中MC法選用了基于RS的MC法(MC-RS)和基于LHS的MC法(MC-LHS).本算例所關心的輸出量為下平板的壓力分布，即對下平板各點處的壓力p進行不確定度量化，得到其均值與方差.

應用MC法進行不確定度量化的主要步驟如下：

(1) 確定隨機變量的分布；

(2) 用蒙特卡洛法生成采樣點；

(3) 在這些采樣點上進行CFD計算；

(4) 統(tǒng)計這些計算結果，得到所關心的量的均值、方差等統(tǒng)計特性.

應用基于回歸的NIPC法進行不確定度量化的主要步驟，其前3步與應用MC法的相同，其后步驟如下：

(4) 根據(jù)隨機變量分布，確定混沌多項式的形式；

(5) 根據(jù)計算結果，擬合出混沌多項式的待定系數(shù)；

(6) 根據(jù)所求系數(shù)，得出所關心的量的均值、方差等統(tǒng)計特性.

計算結果如圖3～6所示.其中各幅圖的圖(a)展示了各個采樣點上的CFD計算結果；圖(b)展示了MC法的結果，即對圖(a)所展示的計算結果進行統(tǒng)計分析得到的統(tǒng)計特性；圖(c)展示了兩種方法(MC和NIPC)求得的統(tǒng)計特性的相對誤差η；圖(d)展示了NIPC法的結果，即對圖(a)所展示的計算結果按NIPC法進行后處理，在求取混沌多項式系數(shù)后，按式(8)和(9)得到統(tǒng)計特性.

圖3 MC-RS方法12次模擬結果

圖4 MC-LHS方法12次模擬結果

圖5 MC-LHS方法48次模擬結果

圖6 MC-RS方法400次模擬結果

4.2 結果分析

首先分析采用基于兩種采樣方法的MC法進行不確定度量化研究的結果.對比圖3(c)和圖4(c)可以發(fā)現(xiàn)，MC-RS和MC-LHS兩者得到的均值結果較好，但方差有較大誤差.從圖5(c)和圖6(c)可以看出，隨著模擬次數(shù)的增加，方差的誤差也隨之減小.

再比較采用MC法與NIPC法分別進行不確定度量化研究的結果.從圖3～6的圖(b)可以看出，使用MC法進行不確定度量化，輸出量的統(tǒng)計特性隨著模擬次數(shù)的增加而緩慢收斂；而從圖(d)可以看出，使用NIPC法進行不確定度量化，其結果收斂較快且并不依賴采樣方法.從圖3(d)、圖4(d)、圖5(d)與圖6(c)的比較可以發(fā)現(xiàn)，經(jīng)過12、48次模擬，并用NIPC法處理后所得到的下平板壓力的統(tǒng)計特性，與用MC法模擬400次后統(tǒng)計出的均值、方差相差無幾，說明NIPC法可以更好地對該隨機CFD問題進行不確定度量化，即在生成的采樣點較少時，NIPC法也能很好地進行不確定度量化.

5 結語

本文分別采用MC法和基于回歸的NIPC法研究了由于邊界條件引起的不確定性在平板間阻曳流流場中的傳播，驗證了NIPC法求解不確定度量化問題的有效性.結果表明：應用基于不同采樣方法的MC法進行不確定度量化的結果相差不大，但效果都遠不如NIPC法.

下一步將開展混沌多項式展開的階數(shù)對不確定性傳播的影響的研究，以及基于正交的NIPC法與基于回歸的NIPC法的比較.就實際的船舶水動力學問題而言，隨機因素也大量存在，如海洋環(huán)境中的隨機海流、隨機海浪，船體制造工藝引起的誤差，船體質(zhì)量分布偏差等，今后可將NIPC法應用于隨機的基于CFD的船舶水動力性能預報問題.