基于侵入混沌多項式法的隨機多孔介質(zhì)內(nèi)順磁性流體熱磁對流不確定度量化1)

姜昌偉 謝云峰 石 爾 劉代飛 李 杰 胡章茂

(長沙理工大學(xué)能源與動力工程學(xué)院,長沙 410114)

引言

重力與磁場力相互作用引起的熱磁對流現(xiàn)象廣泛應(yīng)用于提高晶體生產(chǎn)質(zhì)量、驅(qū)動無機械元件及研究氧濃度傳感器等[1-2].數(shù)值計算是研究熱磁對流驅(qū)動作用機理的重要手段之一,研究人員已開展了許多磁致順磁性流體熱磁對流方面的數(shù)值計算研究[3].He 等[4]提出了一種精確的分離譜元方法系統(tǒng)研究了不同磁感應(yīng)強度與加熱部件位置下矩形方腔內(nèi)空氣熱磁對流的變化規(guī)律.Kaneda 等[5]采用數(shù)值研究方法揭示了磁場力影響方腔內(nèi)順磁性流體熱磁對流的作用機理.Zeng 等[6]和姜昌偉等[7]數(shù)值研究了梯度磁場作用下多孔介質(zhì)方腔內(nèi)順磁性或逆磁性流體的傳熱特性.Zhang 等[8]應(yīng)用Lattice-Boltzmann 方法研究了四極磁場作用下多孔介質(zhì)方腔內(nèi)順磁性流體的熱磁對流與熵產(chǎn).

目前,梯度磁場下多孔介質(zhì)內(nèi)流體熱磁對流的數(shù)值計算普通采用確定性數(shù)學(xué)與幾何結(jié)構(gòu)、邊界條件與物性參數(shù)等.然而現(xiàn)實中的熱磁對流問題中數(shù)理模型、物性參數(shù)與邊界條件等存在著大量的不確定性,如果忽略這些不確定性因素將不能準(zhǔn)確反映真實的熱磁對流過程,因此理論上講必須考慮這些不確定性因素對熱磁對流控制方程及其數(shù)值模擬結(jié)果的影響.梯度磁場下多孔介質(zhì)內(nèi)流體熱磁對流過程中,由于多孔介質(zhì)具有空間變異性,造成多孔介質(zhì)孔隙空間分布呈現(xiàn)出一定的隨機性,確定性熱磁對流控制方程難以描述多孔介質(zhì)孔隙率不確定性對熱磁對流的影響,因此需要發(fā)展出一種把孔隙率作為隨機輸入的熱磁對流不確定性量化分析方法.

不確定性量化研究已成為近年來國際上熱門的課題[9].傳統(tǒng)的不確定性量化方法有蒙特卡羅(Monte-Carlo)法和攝動(Perturbation)法等.葉坤等[10]基于不確定性與全局靈敏度分析方法分析了高超聲速舵面氣動熱的不確定度.Wang 和Qiu[11]基于攝動有限元法,提出了模糊攝動有限元法,用以預(yù)測材料性質(zhì)、外部載荷與邊界條件等具有模糊參數(shù)的不確定性溫度場.雖然蒙特卡羅法需要消耗大量的計算資源,但由于該方法可以應(yīng)用現(xiàn)有的求解器進行求解,因此其廣泛應(yīng)用于計算流體力學(xué)不確定度量化,同時其他計算流體力學(xué)不確定性分析方法均采用該方法作為基準(zhǔn)方法;而攝動法一般應(yīng)用于隨機參數(shù)變異系數(shù)比較小的場合.

近年來,為解決傳統(tǒng)不確定性分析方法高昂的計算成本,研究者先后提出隨機配點法、隨機有限元法、譜方法及多項式混沌展開法等多種數(shù)值計算方法,其中多項式混沌展開法應(yīng)用最為廣泛.多項式混沌展開法最早由Wiener[12]提出,早期在解決湍流問題中得到了較多應(yīng)用.根據(jù)與求解器的耦合方式不同,多項式混沌展開法分為非侵入式(non-intrusive)不確定性量化方法和侵入式(intrusive)不確定性量化方法等.

非侵入式方法基于概率方法,通過在隨機空間內(nèi)采樣獲得若干配置點,然后在各配置點上調(diào)用現(xiàn)有的求解器運行求解,并基于統(tǒng)計方法獲得相關(guān)參數(shù)的統(tǒng)計特征.非侵入式方法包括概率配點法、隨機配置法、稀疏網(wǎng)格配點法等,該方法有利于采用經(jīng)過完好確認(rèn)的CFD 軟件,避免了修改控制方程給工業(yè)應(yīng)用引入新誤差的風(fēng)險.目前,非侵入式方法廣泛應(yīng)用于流體流動與傳熱、燃燒與化學(xué)反應(yīng)等領(lǐng)域的不確定性量化[13-17].Narayan 和Zhou[18]與Guo 等[19]提出了非結(jié)構(gòu)化網(wǎng)格隨機配置與最小二乘法多項式近似兩種不確定性量化方法,并從理論上給出了這兩種方法的穩(wěn)定性與精度的結(jié)果,顯示出其在不確定性量化中的巨大潛力.基于非侵入式方法,黃明等開展了渦輪動葉凹槽狀葉頂氣熱性能、高超聲速吸氣推進飛行器進氣道性能、各向異性陶瓷基復(fù)合材料渦輪葉片熱響應(yīng)等的不確定性量化分析,數(shù)值模擬與實驗數(shù)據(jù)吻合一致驗證了非侵入式方法的有效性,為工程中的高維不確定性量化問題提供了很好的解決思路[20-24].Rajabi 等[25]對多孔介質(zhì)方腔內(nèi)雙擴散自然對流基準(zhǔn)問題進行了詳細(xì)的不確定性傳播分析及全局靈敏度分析以確定模型參數(shù)不確定性對腔體內(nèi)流動、傳熱與傳質(zhì)過程的影響.

侵入式方法在隨機維數(shù)較低的情況下具有一定的精度優(yōu)勢.但該方法需要求解一個復(fù)雜的控制方程組,同時需要重新開發(fā)新的求解器,這增加了其在計算流體力學(xué)應(yīng)用中的難度,并且隨機維數(shù)一旦增加,這種優(yōu)勢立即消除.Xiu 和Karniadakis[26]提出了廣義多項式混沌展開法,并給出了每個不確定性參數(shù)分布最佳正交基函數(shù)的證明.王曉東和康順[27]等將譜隨機不確定性法與有限差分法相結(jié)合,得到不確定性場控制方程的隨機解耦方程組,發(fā)展出用來隨機自然對流模擬的不確定性量化方法.Li 等[28]基于廣義多項式混沌展開和Galerkin 投影將隨機對流擴散方程轉(zhuǎn)化為一組耦合的確定性方程組,研究了兩類隨機輸入分布數(shù)值格式對不確定性量化結(jié)果的影響.Em-Alrani 等[29]和Ma 和Zabaras[30]提出了一種穩(wěn)定化隨機有限元方法來求解隨機不可壓縮Navier-Stokes 方程與隨機多孔介質(zhì)內(nèi)流體自然對流.Chakraborty 和Chowdhury[31]提出了一種稱為基于廣義方差分析(ANOVA)的Galerkin 投影新算法,并成功應(yīng)用于求解隨機Navier-Stokes 方程與隨機穩(wěn)態(tài)擴散問題.由于侵入式方法需要修改控制方程及重新編寫程序代碼,這給侵入式方法的工程應(yīng)用帶來了困難,因此近年來侵入式方法的相關(guān)文獻很少.

本文提出了一種基于侵入式多項式混沌展開法的隨機多孔介質(zhì)內(nèi)順磁性流體熱磁對流不確定性量化方法,構(gòu)建出熱磁對流不確定性量化的隨機數(shù)理模型與有限元程序框架,對隨機多孔介質(zhì)內(nèi)順磁性流體熱磁對流進行不確定性量化.

1 物理模型與控制方程

1.1 物理模型

圖1 給出熱磁對流物理模型和坐標(biāo)系統(tǒng),4 塊磁體置于充滿順磁性流體的多孔介質(zhì)方腔四周.多孔介質(zhì)方腔左右側(cè)壁面分別等溫加熱與冷卻,其他壁面絕熱.方腔無量綱尺寸、磁體無量綱尺寸、方腔與磁體無量綱距離分別為1,5/6,1/8.

圖1 熱磁對流物理模型Fig.1 Physical model of thermomagnetic convection

1.2 控制方程

順磁性流體熱磁對流控制方程包括質(zhì)量守恒方程、多孔介質(zhì)動量守恒方程和能量守恒方程,其無量綱控制方程如下.

質(zhì)量守恒方程

動量守恒方程

能量守恒方程

熱磁對流無量綱控制方程中磁場分布采用標(biāo)量磁位法計算,通過求解Maxwell 方程組獲得磁感應(yīng)強度分布B

其中,B為無量綱磁感應(yīng)強度,(Bx,By);H為無量綱磁場強度,(Hx,Hy);U為無量綱速度,(U,V);P為無量綱壓力;θ為無量綱溫度;τ為無量綱時間;ε為孔隙率;eg為重力方程的單位矢量;Pr為普朗特數(shù);Ra為瑞利數(shù);Da為達西數(shù);γ為磁場力數(shù).

2 侵入式不確定性量化方法

2.1 輸出隨機參數(shù)Karhunen-Loeve (KL)展開

考慮多孔介質(zhì)孔隙率ε具有不確定性,以ε(X,ω)代表輸入隨機參數(shù)孔隙率,其中X∈D(參數(shù)空間),ω ∈Ω(概率空間).輸入隨機參數(shù)可以分解為ε(X,ω)=εˉ(X)+ε′(X,ω),其中 εˉ(X) 表示孔隙率均值,ε′(X,ω)為孔隙率擾動.通常采用協(xié)方差函數(shù)Rhh(X1,X2)=〈ε′(X1,ω),ε′(X2,ω)〉來描述輸入隨機參數(shù)隨機場的空間性質(zhì),協(xié)方差函數(shù)可分解為[32]

式中 λi和 φi分別為特征值和特征函數(shù).特征值和特征函數(shù)可通過求解Fredholm 方程的分析解得到

由此,ε(X,ω) 可表示為

式中 ξi(ω) 為0 均值單位方差的正交高斯隨機變量.

將KL 展開截斷成有限項,截取前m項為KL 展開最終表述形式為

對于二維隨機場,協(xié)方差函數(shù)可寫為Rhh(X1,X2)=σ2exp(-|X1-X2|/η-|Y1-Y2|/η),其中σ2與η 分別為方差及相關(guān)長度.

2.2 輸出響應(yīng)量多項式混沌展開

輸出響應(yīng)量分布特征未知,可以采用多項式混沌展開描述輸出響應(yīng)量不確定性,其緊湊格式可表示為

實際應(yīng)用中多項式混沌擴展通常截斷成有限項.因此,方程(10)可簡化為

擴展項總數(shù)(P+1)由隨機維數(shù)m與多項式展開最高階數(shù)n決定

因此控制方程中的速度、溫度及壓力通過多項式混沌展開可表示為

2.3 輸入隨機參數(shù)非線性函數(shù)隨機譜展開

孔隙率 ε(X,ω) 基于KL 展開寫成

控制方程(2) 包含 ε(X,ω) 的非線性函數(shù),因此首先需要應(yīng)用拉丁超立方抽樣蒙特卡羅法將非線性函數(shù)表示為多項式基

2.4 基本譜分解技術(shù)解耦無量綱控制方程

將方程式(13)～ 式(15)代入控制方程式(1)～式(3),可得方程式(16)～ 式(18)

然后對方程式(16)～ 式(18)進行Galerkin 映射,并基于多項式的正交性可得方程

采用同樣的方法可以得到速度與溫度邊界條件:(1)速度邊界

(2)溫度邊界

附錄給出了m=1,n=2 時解耦控制方程式(19)～式(21)的多項式混沌展開式.

2.5 不確定性量化程序流程

基于侵入式多項式混沌展開法的熱磁對流不確定性量化程序流程如圖2 所示,該算法的計算流程如下:

圖2 不確定性量化流程圖Fig.2 Flowchart of uncertainty quantification

(1)應(yīng)用KL 展開對已知協(xié)方差函數(shù)的輸入隨機參數(shù)進行展開,應(yīng)用拉丁超立方抽樣蒙特卡羅法對多孔介質(zhì)孔隙率非線性項進行多項式混沌展開;

(2)采用多項式混沌對輸出響應(yīng)量進行展開,基于譜分解技術(shù)將熱磁對流控制方程解耦成一組確定性控制方程;

(3)應(yīng)用磁場控制方程求解熱磁對流控制方程中的磁場力源項 ?B2;

(4)通過插值方法對輸入隨機參數(shù)特征值與特征函數(shù)、磁場分布進行插值,應(yīng)用有限元求解解耦熱磁對流控制方程獲得輸出響應(yīng)量多項式混沌展開系數(shù);

(5)利用所得多項式混沌展開系數(shù),評估輸出響應(yīng)量的統(tǒng)計特性、輸出響應(yīng)量的混沌效應(yīng)和輸出響應(yīng)量的概率密度分布函數(shù)(PDF)與累積分布函數(shù)(CDF).

3 結(jié)果分析與討論

隨機多孔介質(zhì)內(nèi)順磁性流體熱磁對流不確定性量化研究中固定Pr=1,Ra=1 × 105,Da=1 × 10-3,γ=25,(X)=0.6,η=1.0,σ=0.1.圖3 與圖4 分別給出了輸入隨機參數(shù)KL 展開的特征值與特征函數(shù).由圖3 可知,特征值隨著隨機維數(shù)增加而快速衰減,取前6 項特征值和特征函數(shù)即滿足要求,因此本研究中取m=6.表1 給出了不同多項式混沌展開階數(shù)下熱壁面平均Nusselt 數(shù)的平均值和標(biāo)準(zhǔn)偏差,當(dāng)n=2 時可以保證不確定性量化精度.因此,P+1=(m+n)!/(m!n!)=(6+2)!/(6!2!)=28,共需耦合求解28組控制方程,即112 個控制方程.解耦熱磁對流控制方程采用有限元求解,網(wǎng)格劃分?jǐn)?shù)目為60 × 60.

圖3 孔隙率KL 展開中協(xié)方差核特征值Fig.3 Eigenvalues of covariance kernel used in the KL expansion of porosity

圖4 孔隙率KL 展開中協(xié)方差核特征函數(shù)Fig.4 Eigenfunction of covariance kennel used in the KL expansion of porosity

圖4 孔隙率KL 展開中協(xié)方差核特征函數(shù)(續(xù))Fig.4 Eigenfunction of covariance kennel used in the KL expansion of porosity (continued)

表1 多項式混沌展開階數(shù)對平均Nusselt 數(shù)均值與標(biāo)準(zhǔn)偏差的影響Table 1 Influence of order of polynomial chaos expansion on mean value and standard deviation of average Nusselt number

3.1 輸出響應(yīng)量均值與標(biāo)準(zhǔn)偏差

圖5 與圖6 分別給出了基于侵入式多項式混沌展開法與蒙特卡羅法的輸出響應(yīng)量均值與標(biāo)準(zhǔn)偏差,兩種方法下的輸出響應(yīng)量均值與標(biāo)準(zhǔn)偏差非常吻合,表明侵入式多項式混沌展開法能夠有效地模擬多孔介質(zhì)孔隙率不確定性在順磁性流體熱磁對流的演化與傳播.從圖5 可看出,順磁性流體受到磁浮升力與重力浮升力共同作用形成相互對持的局面,方腔內(nèi)順磁性流體流動呈現(xiàn)分層流動現(xiàn)象,方腔上部形成一個逆時針方向流動的漩渦而下部形成一個順時針方向流動的漩渦.溫度場分布規(guī)律表明,熱壁面等溫線稠密區(qū)出現(xiàn)上下部而冷壁面等溫線稠密區(qū)出現(xiàn)在中部.從圖6 可以看出,由于孔隙率不確定性通過動量方程演化進而影響能量方程,因此其對溫度場的影響明顯小于對速度場的影響.由于腔體壁面的黏滯效應(yīng),流函數(shù)最大標(biāo)準(zhǔn)偏差區(qū)域出現(xiàn)在方腔左側(cè)中部,而方腔壁面附近流函數(shù)標(biāo)準(zhǔn)偏差較小,溫度場最大標(biāo)準(zhǔn)偏差區(qū)域出現(xiàn)方腔的左右側(cè)中部.

圖5 不同不確定性量化方法下流函數(shù)與溫度場均值Fig.5 Mean value of stream function and temperature field with different uncertainty quantification methods

圖6 不同不確定性量化方法下流函數(shù)與溫度場標(biāo)準(zhǔn)偏差Fig.6 Standard deviation of stream function and temperature field with different uncertainty quantification methods

圖7 是基于侵入式多項式混沌展開法與蒙特卡羅法的熱壁面局部Nusselt 數(shù)分布均值與標(biāo)準(zhǔn)偏差比較,可以看出兩者吻合良好.從熱壁面局部Nusselt 數(shù)均值變化規(guī)律來看,順磁性流體在重力與磁場力的相互作用下,由于方腔左側(cè)壁面下部區(qū)域磁場力與重力相互協(xié)同,因此該區(qū)域傳熱能力增強;而上部區(qū)域,磁場力與重力相互對抗,傳熱能力減弱.熱壁面局部Nusselt 數(shù)標(biāo)準(zhǔn)偏差與均值具有類似的變化規(guī)律.圖8 給出了兩種方法下熱壁面平均Nusselt 數(shù)概率密度分布函數(shù)與累積分布函數(shù),兩者吻合良好.

圖7 不同不確定性量化方法下熱壁面局部Nusselt 數(shù)分布Fig.7 Distribution of local Nusselt number at the hot wall with different uncertainty quantification methods

圖8 不同不確定性量化方法下熱壁面平均Nusselt 數(shù)概率密度函數(shù)和累積分布函數(shù)Fig.8 Probability density function and cumulative distribution function of average Nusselt number at hot wall with different uncertainty quantification methods

3.2 輸出響應(yīng)量混沌效應(yīng)

圖9 給出了一階模式下的流函數(shù),第1 模式流函數(shù)與流函數(shù)平均值具有相似的輪廓,而流函數(shù)第2～ 6 個模式呈現(xiàn)出混沌效應(yīng).圖10 給出了一階模式下的溫度場,一階模式下的溫度場在溫度場不確定性中占主導(dǎo)地位,第1 模式下的溫度場與溫度場標(biāo)準(zhǔn)偏差相類似,并且呈現(xiàn)出與溫度場平均值明顯的差異,而溫度場第2～ 6 個模式具有一定的混沌效應(yīng).從控制方程可知,多孔介質(zhì)孔隙率不確定性主要通過動量方程影響著方腔內(nèi)順磁性流體流動,而順磁性流體流動通過能量方程演化導(dǎo)致溫度場波動,即溫度場不確定性主要由速度場控制,這就造成孔隙率不確定性對溫度場的影響小于對速度場的影響.最后,應(yīng)該注意的是流函數(shù)與溫度場的平均值至少比一階模式下的流函數(shù)與溫度場高一個數(shù)量級.圖11 給出了二階模式下的部分流函數(shù)與溫度場.二階模式下的流函數(shù)與溫度場呈現(xiàn)出顯著的混沌效應(yīng),一階模式下的流函數(shù)與溫度場數(shù)值遠(yuǎn)遠(yuǎn)高于二階模式下的流函數(shù)與溫度場數(shù)值,反映出隨機輸出響應(yīng)量多項式混沌展開具有快速收斂的特性.

圖9 流函數(shù)一階模式Fig.9 First-order modes of stream function

圖9 流函數(shù)一階模式(續(xù))Fig.9 First-order modes of stream function (continued)

圖10 溫度場一階模式Fig.10 First-order modes of temperature field

圖11 流函數(shù)與溫度場二階模式Fig.11 Second-order modes of stream function and temperature field

4 結(jié)論

提出了一種基于侵入式多項式混沌展開的隨機多孔介質(zhì)內(nèi)順磁性流體熱磁對流不確定性量化方法,研究了孔隙率不確定性在熱磁對流控制方程中的傳播規(guī)律與演化特征,主要結(jié)論如下.

(1) 基于侵入式多項式混沌展開實現(xiàn)了隨機多孔介質(zhì)內(nèi)順磁流性流體熱磁對流的不確定性量化.通過譜分解將熱磁對流隨機控制方程轉(zhuǎn)化為一組解耦確定性控制方程并編制了相應(yīng)的計算程序.該方法為計算流體力學(xué)不確定性量化問題的求解提供了一條有效的途徑.

(2) 量化了隨機多孔介質(zhì)內(nèi)順磁性流體熱磁對流的傳播規(guī)律與統(tǒng)計特征.輸入隨機參數(shù)通過隨機熱磁對流控制方程傳播與演化,進而影響著順磁性流體熱磁對流,順磁性流體熱磁對流呈現(xiàn)出顯著的混沌效應(yīng).輸出響應(yīng)量呈現(xiàn)出快速收斂的特性,一階模式下流函數(shù)與溫度場數(shù)值至少比相應(yīng)的平均值低一個數(shù)量級,而二階模式下流函數(shù)與溫度場數(shù)值遠(yuǎn)小于一階模式的數(shù)值.

(3) 比較了基于侵入式多項式混沌展開法與蒙特卡羅法兩種方法下隨機多孔介質(zhì)內(nèi)順磁性流體熱磁對流問題的統(tǒng)計特征,統(tǒng)計特征非常吻合.計算效率方面,與蒙特卡羅法相比,侵入式多項式混沌展開法在相同的精度下具有較小的計算量.

附錄

隨機維數(shù)m=1,多項式混沌展開最高階數(shù)n=2 時控制方程式(19)～ 式(21)的多項式混沌展開式如下

模式0 (方程式(A1)～ 式(A3))

模式1 (方程式(A4)～ 式(A6))

模式2 (方程式(A7)～式(A10))