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

(上海交通大學(xué) a. 船舶海洋與建筑工程學(xué)院； b. 海洋工程國家重點(diǎn)實(shí)驗(yàn)室， 上海 200240)

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

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

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

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

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

1 數(shù)學(xué)模型

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

1.1 問題定義

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

Y=f(x,t,ξ)

(1)

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

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

(2)

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

為了求解實(shí)際問題，通常需要將式(2)進(jìn)行截?cái)?若將混沌多項(xiàng)式展開到p階，則截?cái)嗪蟮亩囗?xiàng)式形式如下：

(3)

式中：混沌多項(xiàng)式的系數(shù)共有Nc個(gè)，

(4)

1.2 混沌多項(xiàng)式

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

表1 不同分布對(duì)應(yīng)的混沌多項(xiàng)式

1.3 非侵入式混沌多項(xiàng)式方法

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

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

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

(5)

寫成矩陣與向量的形式為

Φc=y

(6)

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

c=(ΦTΦ)-1ΦTy

(7)

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

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

2 采樣方法

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

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

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

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

3 二維隨機(jī)阻曳流算例

3.1 求解實(shí)例

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

圖1 二維阻曳流模型

3.2 控制方程

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

(10)

x方向運(yùn)動(dòng)方程：

(11)

y方向運(yùn)動(dòng)方程：

(12)

(13)

式中：τ為切應(yīng)力，滿足牛頓流體本構(gòu)方程.

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

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

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

4 計(jì)算結(jié)果與討論

4.1 計(jì)算結(jié)果

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

應(yīng)用MC法進(jìn)行不確定度量化的主要步驟如下：

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

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

(3) 在這些采樣點(diǎn)上進(jìn)行CFD計(jì)算；

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

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

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

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

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

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

圖3 MC-RS方法12次模擬結(jié)果

圖4 MC-LHS方法12次模擬結(jié)果

圖5 MC-LHS方法48次模擬結(jié)果

圖6 MC-RS方法400次模擬結(jié)果

4.2 結(jié)果分析

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

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

5 結(jié)語

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

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