求解非均質(zhì)多孔介質(zhì)中隨機(jī)水流問(wèn)題的多尺度有限元降基方法

黃夢(mèng)杰, 賀新光

(1.湖南師范大學(xué) 資源與環(huán)境科學(xué)學(xué)院， 湖南 長(zhǎng)沙 410081；2.湖南師范大學(xué) 地理空間大數(shù)據(jù)挖掘與應(yīng)用湖南省重點(diǎn)實(shí)驗(yàn)室， 湖南 長(zhǎng)沙 410081)

1 研究背景

地下多孔介質(zhì)的水力特性通常是高度非均質(zhì)的，且在多種不同的空間尺度上發(fā)生變化。由于對(duì)介質(zhì)水力特性空間分布的不完全認(rèn)知和測(cè)量誤差，在實(shí)際應(yīng)用中，通常視水力特性參數(shù)的空間分布為空間隨機(jī)場(chǎng)，以致對(duì)多孔介質(zhì)中水流運(yùn)動(dòng)的模擬預(yù)測(cè)總是存在一定的不確定性。為了處理這種不確定性，有必要進(jìn)行地下多孔介質(zhì)中水流過(guò)程的隨機(jī)模擬[1-2]。然而，介質(zhì)水力特性的多尺度非均質(zhì)性與不確定性導(dǎo)致發(fā)展高效模擬非均質(zhì)多孔介質(zhì)中隨機(jī)多尺度水流模型的數(shù)值算法極具挑戰(zhàn)性[3]。如果應(yīng)用傳統(tǒng)的數(shù)值方法模擬這樣的隨機(jī)多尺度模型，其計(jì)算過(guò)程是非常耗時(shí)的，有時(shí)甚至是不可行的[4]。正如Xu[5]和Chung 等[6]所指出：以不確定性數(shù)據(jù)為特征的隨機(jī)多尺度問(wèn)題，目前依然是公開的、極具挑戰(zhàn)性的研究領(lǐng)域之一。因此，最近十幾年來(lái)，在計(jì)算科學(xué)與工程領(lǐng)域，發(fā)展模擬這類隨機(jī)多尺度數(shù)學(xué)物理模型的隨機(jī)多尺度數(shù)值方法的興趣一直穩(wěn)定地增長(zhǎng)[5-12]。

模擬地下多孔介質(zhì)中水流運(yùn)動(dòng)的挑戰(zhàn)之一是多孔介質(zhì)特性的空間多尺度性。這類問(wèn)題因分辨最細(xì)空間尺度信息所需的花費(fèi)，致使其在空間上直接的細(xì)尺度分辨模擬在計(jì)算上是非常復(fù)雜而耗時(shí)的，也是不切實(shí)際的。然而，直接忽略細(xì)尺度信息可能導(dǎo)致較大的誤差。為了切實(shí)有效地處理大尺度區(qū)域上的多尺度問(wèn)題，各種多尺度數(shù)值方法已被提出且受到了越來(lái)越多的重視，如：多尺度有限元方法[13](MsFEM)、變分多尺度方法[14]和非均質(zhì)多尺度方法[15]等。其中，MsFEM是多尺度數(shù)值方法的先驅(qū)之一，且許多別的多尺度方法共享了它的相似性[4,16]。一些研究已成功應(yīng)用MsFEM數(shù)值模擬地下非均質(zhì)多孔介質(zhì)中的水流問(wèn)題[3,16-19]。然而，當(dāng)MsFEM與隨機(jī)方法直接結(jié)合用于模擬隨機(jī)多尺度水流模型時(shí)，勢(shì)必導(dǎo)致需要計(jì)算大量的依賴于隨機(jī)參數(shù)樣本的多尺度有限元基函數(shù)，這是十分不利的[3]。因此，為了高效模擬隨機(jī)多尺度地下水流問(wèn)題，研究構(gòu)造獨(dú)立于隨機(jī)參數(shù)取樣的多尺度基函數(shù)是十分必要。

除了多尺度特性，模型輸入量的不確定性也為多孔介質(zhì)中的水流模擬帶來(lái)巨大的挑戰(zhàn)。為了盡量減少計(jì)算成本而又能維持良好的計(jì)算精度，近十幾年來(lái)，不少研究者基于Galerkin投影提出了一些求解隨機(jī)模型的降基方法(RBM)。RBM是一種有效的隨機(jī)模型降階方法，已成功應(yīng)用于各種研究與工程領(lǐng)域[20-24]。RBM依降基函數(shù)構(gòu)造方法的不同，可分為本征正交分解[25-26]( POD)和貪婪算法[20, 22](Greedy)。RBM的基本思想是：基于一組隨機(jī)抽取的參數(shù)樣本，構(gòu)造一組獨(dú)立于參數(shù)樣本的降基函數(shù)，然后對(duì)于任意新的隨機(jī)取樣，投影參數(shù)模型到降基函數(shù)空間而生成一個(gè)降階模型(ROM)。為加速在線計(jì)算，在降基方法中，離線-在線計(jì)算分離是強(qiáng)制性的。在離線階段，大量的基礎(chǔ)計(jì)算包括物理空間的有限元離散、降基函數(shù)的構(gòu)造和全局離散矩陣的組裝等。在在線階段，通過(guò)任意多次抽取參數(shù)樣本，求解參數(shù)化的降階模型，感興趣的輸出量可以被計(jì)算。然而，對(duì)于隨機(jī)多尺度模型，假如應(yīng)用傳統(tǒng)的有限元方法構(gòu)造全局的降基函數(shù)，則離線階段的計(jì)算因需分辨最細(xì)空間尺度信息而十分耗時(shí)。為降低離線階段的計(jì)算成本，最新的研究[4,27]建議：應(yīng)用MsFEM計(jì)算局部輔助函數(shù)，并進(jìn)一步應(yīng)用POD生成多尺度降基函數(shù)空間以構(gòu)造隨機(jī)多尺度問(wèn)題的降階模型。

本文遵循最新研究[4]為處理橢圓型隨機(jī)多尺度問(wèn)題的主要思想，通過(guò)有效耦合多尺度有限元和降基方法而發(fā)展一種高效求解非均質(zhì)多孔介質(zhì)中的隨機(jī)多尺度水流問(wèn)題(拋物型)的多尺度有限元降基方法(RMsBM)。而且，為了加速在線計(jì)算，不同于文獻(xiàn)[4]，本文采用一種矩陣離散的經(jīng)驗(yàn)插值方法[28](MDEIM)仿射分解參數(shù)依賴的代數(shù)系統(tǒng)以實(shí)現(xiàn)水流降階模型的離線-在線計(jì)算分離。Negri 等[28]指出：MDEIM方法以一種非浸潤(rùn)、高效而純代數(shù)的方式處理非仿射的參數(shù)化離散系統(tǒng)。最后，為表明所提出的RMsBM的性能與效率，本文執(zhí)行了幾個(gè)有代表性的數(shù)值試驗(yàn)。結(jié)果表明：本文為求解非均質(zhì)多孔介質(zhì)中的隨機(jī)多尺度飽和水流問(wèn)題提出了一種高效的數(shù)值模擬方法。

2 模型問(wèn)題

本節(jié)介紹非均質(zhì)隨機(jī)多孔介質(zhì)中的承壓水流控制方程及其全階模型和降階模型。

2.1 水流控制方程

令D?R2為二維有界物理域，(0,T]為時(shí)間域，D×(0,T]為笛卡兒積，μ∈Ω?Rp為p維隨機(jī)參數(shù)向量。非均質(zhì)隨機(jī)多孔介質(zhì)中的非穩(wěn)定飽和水流可用含參數(shù)的拋物型偏微分方程描述如下。

R(x,t) (x,t)∈D×(0,T]

(1)

初始條件：

ψ(x,t,μ)|t=0=ψ0(x)

(2)

邊界條件：

ψ(x,t,μ)|ΓD=ψΓ(x,t)

(3)

-K(x,μ)ψ(x,t,μ)·n|ΓN=q(x,t)

(4)

式中：x∈D為空間坐標(biāo)；ψ為水力水頭；Ss為比貯水系數(shù)；K為滲透系數(shù)；R為源/匯項(xiàng)，n為單位外法向量；ΓD和ΓN分別為區(qū)域的Dirichlet和Neumann邊界，且ΓD∩ΓN=?。

對(duì)于非均質(zhì)多孔介質(zhì)，K一般在隨機(jī)參數(shù)空間振蕩且在物理空間高度非均質(zhì)。然而，在任何情況下，Ss的變化幅度均小于lnK[29]，因此，本研究假設(shè)Ss為區(qū)域D內(nèi)的一個(gè)確定性函數(shù)。對(duì)于含參數(shù)的水流問(wèn)題(1)，解的統(tǒng)計(jì)性質(zhì)如水頭均值和方差是令人感興趣的。為此，可首先利用Monte Carlo或稀網(wǎng)格隨機(jī)配點(diǎn)方法對(duì)隨機(jī)參數(shù)空間Ω進(jìn)行抽樣，然后，對(duì)于任意給定的參數(shù)μ∈Ω，應(yīng)用有限元等確定性方法求解參數(shù)化的水流問(wèn)題(1)，并可最終獲得解的統(tǒng)計(jì)性質(zhì)。然而，對(duì)于一個(gè)給定的參數(shù)μ，為捕捉細(xì)尺度非均質(zhì)性，仍需要在物理空間執(zhí)行高分辨模擬。這樣，如果直接在一個(gè)細(xì)網(wǎng)格上對(duì)流動(dòng)問(wèn)題執(zhí)行全階模擬，則計(jì)算成本是非常高的。為了提高計(jì)算效率，本文發(fā)展一種基于多尺度有限元基函數(shù)的降基方法求解參數(shù)方程(1)。

2.2 全階模型

(5)

(6)

(7)

式中：矩陣CNf和ANf(μ)在第(i,j)處的元素分別為(CNf)ij和(ANf)ij；向量b的第i個(gè)元素為(bNf)i。本文稱方程式(7)為水流問(wèn)題(1)的全階模型。為了獲得全離散的解，可以采用Euler或其他方法離散半離散的方程式(7)。本文應(yīng)用常用的Crank-Nicolson(CN)方法構(gòu)造方程式(7)的全離散格式。關(guān)于CN全離散格式的更多細(xì)節(jié)，請(qǐng)參閱He等[3]或其他相關(guān)文獻(xiàn)。

2.3 降階模型

(8)

3 多尺度有限元降基方法

3.1 多尺度有限元基函數(shù)的構(gòu)造

(9)

(10)

式中：gi為定義在κ的邊界?κ上且與節(jié)點(diǎn)i相關(guān)的邊界條件；d為粗單元的節(jié)點(diǎn)總數(shù)。

3.2 多尺度降基函數(shù)的生成

Vsnap={φims(μn):i=1,2，…,Nc,n=1,2，…,Nopt}

(11)

(12)

(13)

3.3 離線-在線計(jì)算

(CNr)ij=(φirm,SSφjrm),(ANr)ij

(14)

降階模型(8)全離散后，可得一個(gè)Nr元的線性代數(shù)方程組。根據(jù)方程(14),方程(8)的降階矩陣和向量均涉及φirm(i=1,2，…,Nr)的內(nèi)積計(jì)算。因降基函數(shù)φirm屬于有限元空間Xh，故其可用Nf個(gè)標(biāo)準(zhǔn)的有限元基函數(shù)ξk(k=1,2，…,Nf)∈Xh表示，即：

(15)

令(Z)ki=zki,k=1,2，…,Nf,i=1,2，…,Nr。將方程(15)代入式(14)，可得：

CNr=ZTCNfZ,ANr(μ)=ZTANf(μ)Z,

bNr=ZTbNf

(16)

式中：矩陣Z∈RNf×Nr在第(k,i)處的元素為(Z)ki，CNf、ANf(μ)和bNf已由方程(7)所定義。因?yàn)榫仃嘮和CNf以及向量bNf均獨(dú)立于隨機(jī)參數(shù)μ，所以，降階矩陣CNr和向量bNr均只需在離線階段計(jì)算一次。然而，為了得到降階模型的解ΨNr(μ)，依賴于參數(shù)μ的全階矩陣ANf(μ)需重新計(jì)算而導(dǎo)致巨大的計(jì)算量。因此，ANr(μ)的有效估計(jì)對(duì)于實(shí)現(xiàn)降階模型的離線-在線計(jì)算分離是至關(guān)重要的。如果矩陣ANf(μ)可用依賴于μ的系數(shù)加權(quán)的常數(shù)矩陣的仿射組合表示，則加權(quán)和的每一項(xiàng)可離線投影到多尺度降基函數(shù)空間。假設(shè)ANf(μ)可仿射表示為：

(17)

則有：

(18)

4 數(shù)值模擬

本節(jié)執(zhí)行幾個(gè)數(shù)值算例以表明所提出的多尺度有限元降基方法求解非均質(zhì)多孔介質(zhì)中隨機(jī)飽和水流問(wèn)題的性能。為此，考慮非均質(zhì)含水層對(duì)水位突變的響應(yīng)，即：假如在時(shí)刻t=0時(shí)，突然將含水層邊界上的水位從15 m下降到12 m，并模擬一段時(shí)間內(nèi)水頭隨時(shí)間的變化。假設(shè)含水層的比貯水系數(shù)和厚度分別為8×10-4/m和10 m，則隨機(jī)水流問(wèn)題(1)的初始條件和邊界條件分別為：

ψ(x,0,μ)=15m inD

(19)

ψ(x,t,μ)=12m on(0,T]×?D

(20)

并假設(shè)滲透系數(shù)對(duì)數(shù)Y=ln(K(x,μ))是一個(gè)二階平穩(wěn)的高斯隨機(jī)場(chǎng)，其兩點(diǎn)指數(shù)協(xié)方差函數(shù)cov[Y]為：

cov[Y](x1,y1;x2,y2)=

(21)

式中： (xi,yi)為空間坐標(biāo)；σ2為隨機(jī)場(chǎng)Y的方差；lx和ly分別為x和t方向的相關(guān)長(zhǎng)度。任給一個(gè)協(xié)方差函數(shù)，其隨機(jī)場(chǎng)Y可由截?cái)嗟腒LE展式生成[3]：

(22)

式中：隨機(jī)向量μ:=(μ1,μ2,…,μN(yùn)k)∈RNk,μi(i=1,2，…,Nk)為獨(dú)立同分布的標(biāo)準(zhǔn)正態(tài)隨機(jī)變量，λi和bi(x)分別為協(xié)方差函數(shù)的特征值和特征函數(shù)。此時(shí)，水流方程(1)則是一個(gè)Nk維的隨機(jī)參數(shù)問(wèn)題。因輸入場(chǎng)K(x,μ)=exp(Y(x,μ))關(guān)于參數(shù)μ是非仿射的，為了實(shí)現(xiàn)水流問(wèn)題(1)的ROM的離線-在線計(jì)算分離，本研究應(yīng)用MDEIM對(duì)非仿射的參數(shù)矩陣執(zhí)行仿射分解。

為了不失一般性，本節(jié)所有的數(shù)值試驗(yàn)均考慮同一空間尺度的物理區(qū)域D。令區(qū)域D為1個(gè)200 m × 200 m的矩形。然后，1個(gè)一致的有限元網(wǎng)格將區(qū)域D劃分成Nx×Ny個(gè)子矩形，并稱之為Nx×Ny的網(wǎng)格。時(shí)間離散步長(zhǎng)固定為Δt=5 min。為表明多尺度降基方法的性能，本節(jié)在粗網(wǎng)格上應(yīng)用所提出的RMsBM和在細(xì)網(wǎng)格上應(yīng)用標(biāo)準(zhǔn)的FEM分別求解水流問(wèn)題(1)的ROM與FOM模型，并對(duì)所獲得的解進(jìn)行比較。下列的數(shù)值算例中，標(biāo)準(zhǔn)的FEM采用128×128的細(xì)網(wǎng)格求解水流問(wèn)題的FOM模型，并將其解作為參考解。為了評(píng)估粗網(wǎng)格尺寸對(duì)RMsBM的影響，下面的例子將考慮3種不同的粗網(wǎng)格，即8×8、16×16和32×32的粗網(wǎng)格。為了方便，將RMsBM方法獲得的ROM解稱之為RMsBM解。為估計(jì)RMsBM解的精度，本研究使用如下定義的相對(duì)L2誤差：

(23)

式中：Nf為細(xì)網(wǎng)格上的節(jié)點(diǎn)總數(shù);ψr和ψf分別為水流問(wèn)題的RMsBM解和參考解。

4.1 局部基函數(shù)數(shù)目的影響

本小節(jié)旨在討論局部多尺度基函數(shù)的數(shù)目對(duì)所提出的RMsBM性能的影響。假定隨機(jī)場(chǎng)Y的方差σ2=2.0，其幾何平均等于-6.125，且相關(guān)結(jié)構(gòu)是各向異性的，并具有相關(guān)長(zhǎng)度Lx=40 m和Ly=30 m。此數(shù)值算例中，使用KLE展式，取Nk=25項(xiàng)生成水力參數(shù)傳導(dǎo)場(chǎng)K(x,μ)=exp(Y(x,μ)。此時(shí)，水流方程(1)則是一個(gè)25維的隨機(jī)參數(shù)問(wèn)題，并在一個(gè)16×16的粗網(wǎng)格上求解該問(wèn)題的RMsBM解。為此，在RMsBM的離線階段，首先從隨機(jī)參數(shù)空間選取350個(gè)25維的參數(shù)樣本，并應(yīng)用MsFEM在粗網(wǎng)格上分別計(jì)算所選樣本的局部snapshot函數(shù)，然后，對(duì)所選樣本構(gòu)造的snapshot函數(shù)分別在其相應(yīng)的每一個(gè)局部領(lǐng)域執(zhí)行本征正交分解，最終獲得局部多尺度降基函數(shù)。在RMsBM的在線階段，則使用已構(gòu)造的局部多尺度降基函數(shù)求解任意給定的隨機(jī)參數(shù)樣本所確定的水流問(wèn)題的RMsBM解。

為了評(píng)估局部基函數(shù)數(shù)目對(duì)RMsBM性能的影響，本算例分別考慮粗網(wǎng)格的每個(gè)粗節(jié)點(diǎn)處均使用同一數(shù)目N=1，2，…，5的局部多尺度降基函數(shù)的情形。應(yīng)注意到：128×128細(xì)網(wǎng)格上的FOM離散方程組的自由度為Nf=16 641，然而，當(dāng)N=1,2,…,5時(shí)，16×16粗網(wǎng)格上的ROM離散方程組的自由度僅分別為Nr= 289、578、869、1 156和1 445。顯然，RMsBM的ROM可以提供比標(biāo)準(zhǔn)的FOM明顯更快的在線計(jì)算，且RMsBM的在線計(jì)算隨局部基函數(shù)數(shù)目的減少而加速。為了評(píng)估RMsBM解的近似精度，本例從參數(shù)空間中再隨機(jī)抽取2000個(gè)新的參數(shù)樣本，并使用RMsBM在離線階段已構(gòu)造的多尺度降基函數(shù)空間分別求解每一新參數(shù)樣本所確定的水流問(wèn)題(1)的RMsBM解。然后，分別計(jì)算這2000個(gè)樣本的RMsBM解和參考解的均值和方差。圖1展示了在不同時(shí)刻t=8、14和20 h，每個(gè)粗節(jié)點(diǎn)處使用不同數(shù)目局部多尺度降基函數(shù)的RMsBM解的均值和方差的相對(duì)L2誤差。從圖1可知，當(dāng)局部降基函數(shù)的數(shù)目從1開始增加時(shí)，3個(gè)不同時(shí)刻的水頭均值和方差的相對(duì)誤差均快速衰減，然后，當(dāng)局部基函數(shù)數(shù)目從3增至5時(shí)，所有的相對(duì)誤差都一直保持很小的值。此外，圖1也表明，每個(gè)粗網(wǎng)格節(jié)點(diǎn)處僅具有3個(gè)局部基函數(shù)的RMsBM可為參考解的均值和方差提供良好的近似逼近。說(shuō)明在該情況下，當(dāng)局部降基函數(shù)數(shù)目N=3時(shí)，RMsBM可以在近似精度和在線計(jì)算成本之間取得良好的平衡。具體地，對(duì)于一個(gè)參數(shù)樣本，為獲得其在時(shí)刻t=20 h的解，具有3個(gè)局部基函數(shù)的RMsBM和標(biāo)準(zhǔn)的FEM 的平均在線CPU時(shí)間分別是26和53 s。這表明：相較于標(biāo)準(zhǔn)的FEM，對(duì)于這種情況，具有3個(gè)局部基函數(shù)的RMsBM可節(jié)省大約51%的在線CPU時(shí)間。為了進(jìn)一步比較，圖2和3分別描繪了在時(shí)刻t=8 h，由具有1個(gè)和3個(gè)局部基函數(shù)的RMsBM和標(biāo)準(zhǔn)的FEM所獲得的水頭均值和方差的空間分布。從圖2和3可以看出：由具有不同數(shù)目的局部基函數(shù)的RMsBM獲得的水頭均值均能很好地匹配參考的水頭均值，且由具有3個(gè)局部基函數(shù)的RMsBM得到的水頭方差也與參考的方差具有良好的一致性。然而，由每個(gè)粗網(wǎng)格節(jié)點(diǎn)處僅具有1個(gè)局部基函數(shù)的RMsBM計(jì)算的水頭方差沒(méi)能很好地匹配參考的方差。這表明：使用適當(dāng)數(shù)目的局部多尺度降基函數(shù)對(duì)RMsBM方法的效率與精度非常重要。為了直觀展示不同參數(shù)樣本的個(gè)體誤差，圖4提供了前300個(gè)樣本在時(shí)刻t=8 h、每個(gè)粗節(jié)點(diǎn)處具有不同數(shù)目局部基函數(shù)的RMsBM解的相對(duì)L2誤差。從圖4可以觀察到，具有3個(gè)局部基函數(shù)的RMsBM為這300個(gè)樣本提供了最準(zhǔn)確的解，而僅具有1個(gè)局部基函數(shù)的RMsBM則給予了最不精確的解。此外，圖4也表明，相比于具有1個(gè)或2個(gè)局部基函數(shù)的RMsBM，具有3個(gè)局部基函數(shù)的RMsBM所獲得解的相對(duì)誤差對(duì)參數(shù)樣本的選擇是不敏感的。這些結(jié)果表明：在16×16粗網(wǎng)格上，每個(gè)節(jié)點(diǎn)處具有3個(gè)局部降基函數(shù)的RMsBM可以為128×128細(xì)網(wǎng)格上的參數(shù)水流問(wèn)題的FOM模型提供穩(wěn)定而準(zhǔn)確的ROM近似解。

4.2 粗網(wǎng)格尺寸的影響

本節(jié)使用8×8、16×16和32×32粗網(wǎng)格考察不同粗網(wǎng)格尺寸對(duì)RMsBM性能的影響。為此，假定滲透系數(shù)隨機(jī)場(chǎng)K(x,μ)的所有參數(shù)與4.1節(jié)相同。對(duì)不同的粗網(wǎng)格，RMsBM在其離線階段構(gòu)造多尺度降基函數(shù)之后，在其在線階段，分別求解隨機(jī)選取的2 000個(gè)參數(shù)樣本中每一樣本所確定的水流問(wèn)題(1)的ROM解。圖5展示了使用不同數(shù)目局部基函數(shù)的RMsBM在3個(gè)不同粗網(wǎng)格上在時(shí)刻t=8 h獲得的RMsBM解的均值和方差的相對(duì)誤差。從圖5中可觀察到誤差的兩個(gè)主要特征：一是對(duì)于1個(gè)給定的粗網(wǎng)格，RMsBM解的誤差隨局部基函數(shù)數(shù)目的增加而單調(diào)減少，特別是8×8和16×16粗網(wǎng)格，當(dāng)每節(jié)點(diǎn)處的局部基函數(shù)由1個(gè)增至3個(gè)時(shí)，RMsBM解的誤差顯著減少；二是隨著粗網(wǎng)格的細(xì)化，每個(gè)節(jié)點(diǎn)具有固定數(shù)目局部降基函數(shù)(1、2或3個(gè))的RMsBM解的近似精度普遍提高。然而，每個(gè)節(jié)點(diǎn)具有3個(gè)局部基函數(shù)時(shí)，RMsBM在32×32粗網(wǎng)格上并沒(méi)有提供比其在16×16粗網(wǎng)格上更高精度的解。這表明：在這種情況下，RMsBM存在一個(gè)最優(yōu)的粗網(wǎng)格(16×16)。此外，還觀察到：RMsBM在8×8粗網(wǎng)格上每節(jié)點(diǎn)具有6個(gè)局部基函數(shù)的解和其在16×16粗網(wǎng)格上具有3個(gè)局部基函數(shù)的解具有幾乎一樣的近似精度。

圖1 具有不同數(shù)目局部基函數(shù)的RMsBM解在3個(gè)不同時(shí)刻的均值和方差的相對(duì)L2誤差

圖2 參考解與具有1個(gè)和3個(gè)局部基函數(shù)的RMsBM解在時(shí)刻t=8 h的水頭均值對(duì)比

圖3 參考解與具有1個(gè)和3個(gè)局部基函數(shù)的RMsBM解在時(shí)刻t=8 h的水頭方差對(duì)比

圖4 具有不同數(shù)目局部基函數(shù)的RMsBM對(duì)前300個(gè)參數(shù)樣本在時(shí)刻t=8 h解的相對(duì)L2誤差對(duì)比

圖6給出了8×8粗網(wǎng)格上每節(jié)點(diǎn)具有3個(gè)和6個(gè)局部基函數(shù)的RMsBM解和參考解在時(shí)刻t=8 h的方差的對(duì)比。從圖6可看出：對(duì)于8×8粗網(wǎng)格，當(dāng)每節(jié)點(diǎn)處的局部基函數(shù)從3個(gè)增至6個(gè)時(shí)，RMsBM解的方差可以很好地逼近參考方差。然而，隨每節(jié)點(diǎn)局部基函數(shù)數(shù)目的增加，RMsBM的在線階段則需要更多的計(jì)算成本。顯然，8×8、16×16和32×32粗網(wǎng)格的每個(gè)粗單元分別包含16×16、8×8和4×4的細(xì)網(wǎng)格。從對(duì)RMsBM的描述可以看出：在離線階段構(gòu)造局部snapshot函數(shù)時(shí)，局部離散問(wèn)題的自由度隨著粗網(wǎng)格的細(xì)化而減少，但其局部問(wèn)題的總數(shù)卻隨之增加。因此，本文建議：在這種情況下，為了取得解的近似精度和計(jì)算成本之間的良好平衡，可使用16×16粗網(wǎng)格和每節(jié)點(diǎn)3個(gè)局部基函數(shù)的RMsBM有效地降階水流問(wèn)題的全階模型。

圖5 具有不同數(shù)目局部基函數(shù)的RMsBM在3個(gè)不同粗網(wǎng)格上解的均值和方差的相對(duì)誤差

圖6 參考解與8×8粗網(wǎng)格上具有3個(gè)和6個(gè)局部基函數(shù)的RMsBM解的方差對(duì)比

4.3 空間變異性的影響

本節(jié)探討滲透系數(shù)場(chǎng)的空間變異性(σ2)對(duì)所提出的RMsBM性能的影響。為此，考慮具有3個(gè)不同方差 1.5、2.0和2.5的隨機(jī)參數(shù)輸入場(chǎng)，并假設(shè)其余的相關(guān)參數(shù)與4.1節(jié)中的相同?；?.2節(jié)的討論，本節(jié)僅考慮RMsBM在16×16粗網(wǎng)格上執(zhí)行模型降階。同樣地，對(duì)每個(gè)輸入場(chǎng)隨機(jī)選取的2 000個(gè)樣本，分別求解每個(gè)樣本所確定的水流問(wèn)題(1)的RMsBM解和參考解。顯然，水流模型輸出水頭的空間變異性隨其輸入?yún)?shù)場(chǎng)σ2的增大而增強(qiáng)。當(dāng)σ2較大時(shí)，期望所構(gòu)造的多尺度降基函數(shù)也能在粗網(wǎng)格上有效捕捉輸入場(chǎng)的高空間變異性。為此，圖7展示了具有不同數(shù)目局部基函數(shù)的RMsBM在16×16粗網(wǎng)格上為3個(gè)不同隨機(jī)參數(shù)輸入場(chǎng)獲得的解在時(shí)刻t=8 h的均值和方差的相對(duì)誤差。圖7表明：盡管輸入場(chǎng)方差的不同會(huì)導(dǎo)致輸出水頭空間變異性的不同，但RMsBM對(duì)3個(gè)不同的σ2給出了幾乎相同的水頭均值和方差的相對(duì)誤差。這里，特別注釋：當(dāng)σ2=2.5時(shí)，隨機(jī)滲透系數(shù)場(chǎng)的許多實(shí)現(xiàn)達(dá)7個(gè)數(shù)量級(jí)的變化。

圖7 不同數(shù)目局部基函數(shù)的RMsBM對(duì)不同隨機(jī)場(chǎng)的解在時(shí)刻t=8 h的均值和方差的相對(duì)誤差

這表明：所提出的RMsBM對(duì)具有高空間變異性輸入場(chǎng)的隨機(jī)水流問(wèn)題的降階具有很強(qiáng)的穩(wěn)健性。

圖8對(duì)3個(gè)不同的空間變異，給出了參考解與使用3個(gè)局部基函數(shù)的RMsBM解在時(shí)刻t=8 h和截面y=100 m處的均值和方差的對(duì)比。從圖8可觀察到：具有3個(gè)局部基函數(shù)的RMsBM解的均值和方差均能很好地匹配所有這些情況的參考均值和方差。這表明：在適當(dāng)?shù)拇志W(wǎng)格上使用合適數(shù)目的局部基函數(shù)的RMsBM能有效地捕獲具有高空間變異性的細(xì)尺度異質(zhì)性對(duì)粗尺度解的影響。

圖8 不同隨機(jī)場(chǎng)的參考解與具有3個(gè)局部基函數(shù)的RMsBM解在時(shí)刻t=8 h和截面y=100 m處的均值和方差的對(duì)比

5 結(jié) 論

本文通過(guò)有效耦合多尺度有限元和降基方法，提出并討論了一種求解非均質(zhì)多孔介質(zhì)中的隨機(jī)參數(shù)飽和水流問(wèn)題的多尺度有限元降基方法。該方法的有效性在于構(gòu)造一組獨(dú)立于隨機(jī)參數(shù)取樣的多尺度有限元降基函數(shù)和生成一個(gè)降階多尺度模型，并應(yīng)用矩陣離散的經(jīng)驗(yàn)插值方法仿射分解非仿射的隨機(jī)參數(shù)問(wèn)題的離散系統(tǒng)以加速降階模型的在線計(jì)算。所提出的RMsBM允許其計(jì)算過(guò)程執(zhí)行一個(gè)離線-在線計(jì)算分離，即在RMsBM離線階段，構(gòu)造多尺度有限元降基函數(shù)，而在在線階段，則使用已構(gòu)造的多尺度降基函數(shù)求解任意參數(shù)樣本所確定的水流問(wèn)題的解。RMsBM的離線計(jì)算可能較為耗時(shí)，但其在線計(jì)算是快速有效的，而快速有效的在線計(jì)算對(duì)預(yù)測(cè)任意參數(shù)樣本所對(duì)應(yīng)的模型輸出量和評(píng)估隨機(jī)輸入量對(duì)輸出量的影響是十分重要的。數(shù)值算例表明：

(1)RMsBM通過(guò)選擇一個(gè)合適的粗網(wǎng)格和一個(gè)最優(yōu)數(shù)目的降基函數(shù)，可實(shí)現(xiàn)隨機(jī)參數(shù)水流問(wèn)題求解在近似計(jì)算精度和在線計(jì)算成本之間的良好平衡。

(2)RMsBM對(duì)具有高空間變異性輸入場(chǎng)的隨機(jī)水流問(wèn)題的降階模擬提供了一個(gè)穩(wěn)健而精確的替代模型。

總之，所提出的RMsBM無(wú)需在全局細(xì)網(wǎng)格上分辨所有的隨機(jī)小尺度信息，即可在適當(dāng)?shù)拇志W(wǎng)格上應(yīng)用合適數(shù)目的局部降基函數(shù)，有效捕獲隨機(jī)輸入場(chǎng)的細(xì)尺度非均質(zhì)性對(duì)水流模型輸出場(chǎng)的影響，從而大大降低了計(jì)算的復(fù)雜性，節(jié)省了計(jì)算成本。因此，方法可為大尺度區(qū)域地下水污染狀況的實(shí)時(shí)監(jiān)測(cè)分析和放射性廢物地下儲(chǔ)放的安全性評(píng)價(jià)提供一種有力的數(shù)值模擬工具。