夏艷華

（安徽理工大學(xué) 土木學(xué)院，安徽 淮南232001）

1 引言

地下水滲流模擬是評(píng)估人類各種工程活動(dòng)如地下水利用、石油開采、大壩建設(shè)等的重要手段［1］。由于地層形成經(jīng)歷了復(fù)雜的地質(zhì)演變，地層的各種屬性如滲透系數(shù)表現(xiàn)出多尺度異質(zhì)性，這使得地下水滲流模擬變得極為復(fù)雜。如果運(yùn)用標(biāo)準(zhǔn)的有限元計(jì)算，為了獲取足夠的精度，網(wǎng)格剖分需在細(xì)尺度上進(jìn)行。這將帶來如下兩個(gè)方面的主要困難：一是由于網(wǎng)格剖分在細(xì)尺度上進(jìn)行，將需要更多的計(jì)算資源，當(dāng)研究區(qū)域較大時(shí)，以現(xiàn)有的物質(zhì)基礎(chǔ)，計(jì)算將無法完成；另一是需要知道整個(gè)研究區(qū)域上計(jì)算參數(shù)的變化。然而，事實(shí)上在巖土工程領(lǐng)域，人們只能獲取局部的、有限的、稀疏的數(shù)據(jù)。

目前，研究隨機(jī)滲流場(chǎng)的主要方法是基于概率理論的 Monte Carlo法、矩方程法和攝動(dòng)法［2-5］，其中Monte Carlo法是最為有效的方法。最近有一種稱之為廣義多項(xiàng)式混沌的隨機(jī)計(jì)算方法成為應(yīng)用最為廣泛的方法［6］，此法將隨機(jī)場(chǎng)表達(dá)為輸入隨機(jī)變量的正交多項(xiàng)式，為了得到更好的收斂性，可以根據(jù)輸入的隨機(jī)變量選用不同類型的正交基。

但是這些隨機(jī)計(jì)算方法沒有考慮問題的多尺度性，當(dāng)遇到在多尺度上不均勻的問題時(shí)，這些方法往往顯得無能為力，因此需采用新的手段來模擬此類多尺度問題。目前研究者已提出了多種多尺度計(jì)算方法，在這些方法中，最常用的有多尺度有限元法（Ms-FEM）［7-8］、變分多尺度法（VMS）［9］和異質(zhì)多尺度法（HMM）［10］。

HMM是解決多尺度微分方程的一般方法。其計(jì)算理念由兩部分組成：一是宏觀尺度上求解器的選取，即在宏觀尺度上采用什么計(jì)算手段，目前通常采用有限元法；二是通過求解局部細(xì)觀尺度問題估算所需的宏觀諸如單元?jiǎng)偠染仃嚒⒘魍康任粗?。此?jì)算方法能有效地解決具有高震蕩系數(shù)的橢圓或雙曲線問題。但是多尺度計(jì)算方法同樣面臨著需要詳盡知道細(xì)觀局部求解區(qū)間的計(jì)算參數(shù)如滲透系數(shù)的問題［11］。事實(shí)上，我們只能獲取有限的、稀疏的、局部的數(shù)據(jù)，科學(xué)有效地描述這些數(shù)據(jù)往往需要借助隨機(jī)模型。

綜上所述，將隨機(jī)計(jì)算方法與HMM結(jié)合是解決異質(zhì)多尺度地層滲流問題的必經(jīng)途徑。因此提出結(jié)合有限元異質(zhì)多尺度方法和基于廣義多項(xiàng)式混沌的隨機(jī)配點(diǎn)法的SHMFE法模擬異質(zhì)隨機(jī)多尺度滲透模型。在此方法中，將對(duì)數(shù)滲透系數(shù)Y＝lnKε在細(xì)尺度上用KL分解展開，然后采用稀疏格網(wǎng)配點(diǎn)法離散隨機(jī)變量空間使隨機(jī)模型變?yōu)榇_定性問題［12］，最后采用FEHM法求解此確定性問題。并通過計(jì)算實(shí)例對(duì)SHMFE有效性加以證明。

2 多尺度滲流模型

2.1 微觀模型及均質(zhì)化

2.1.1 微觀模型

這里所研究的滲透模型其滲透系數(shù)取決于細(xì)尺度地層結(jié)構(gòu)且在空間上是不均勻的。在研究區(qū)域D上，滲透模型在細(xì)尺度ε上的微分方程如下：

式中：kε為細(xì)尺度上滲透系數(shù)張量，假定該張量為正定對(duì)稱張量；hε為細(xì)尺度上的水頭；hD為Dirichlet邊界?ΓD上給定水頭；gN為Neumann邊界?ΓN上給定流通量。

問題（1）可用標(biāo)準(zhǔn)的有限元進(jìn)行求解，為了達(dá)到理想的精度，運(yùn)用標(biāo)準(zhǔn)有限元求解時(shí)，需在模型最小尺度上進(jìn)行網(wǎng)格劃分。這意味著模擬需要大量的計(jì)算資源，當(dāng)最小尺度ε遠(yuǎn)遠(yuǎn)小于研究區(qū)域尺度時(shí)，對(duì)當(dāng)前的計(jì)算資源而言求解問題（1）是一個(gè)不可能完成的任務(wù)。這也是發(fā)展均質(zhì)化理論的主要?jiǎng)訖C(jī)。

2.1.2 均質(zhì)化模型

假定多尺度滲流模型是尺度分離的，根據(jù)均質(zhì)化理論［13］問題（1）的解hε弱收斂于問題（2）的解h0：

式中：k0為均質(zhì)化滲透系數(shù)張量，假定該張量為正定對(duì)稱張量。

問題（2）的解h0稱為均質(zhì)化解，它不依賴于細(xì)尺度ε。通常均質(zhì)化滲透系數(shù)張量k0是未知的，需要通過細(xì)尺度模型求解。

2.2 有限元異質(zhì)多尺度法

均質(zhì)化問題（2）可以通過有限元異質(zhì)多尺度法進(jìn)行求解。將問題（2）在宏觀尺度上進(jìn)行有限元近似，由于問題（2）的滲透系數(shù)k0未知，因此需通過微觀尺度問題（1）上的已知數(shù)據(jù)來求解。具體實(shí)現(xiàn)是通過在宏觀單元的每個(gè)積分點(diǎn)上求解微觀模型來計(jì)算宏觀單元的剛度矩陣。

問題（2）的弱解形式如下：

式（3）中的h0滿足Dirichlet邊界條件；v在Dirichlet邊界上為0。

式（3）可以用傳統(tǒng)的有限元法進(jìn)行離散，假定TH為研究區(qū)域D 的網(wǎng)格剖分，K為離散域的有限元，問題（2）弱解形式有限元法進(jìn)行離散的剛度矩陣通過如下的數(shù)值積分進(jìn)行計(jì)算：

式中：｛xl｝為有限元K高斯積分點(diǎn)；｛ωl｝為有限元K積分權(quán)重；｜K｜為有限元K的體積。

由于k0未知，式（4）中（▽v·k0▽h0）（xl）是未知的，該未知量通過求解積分點(diǎn)處局部微觀模型來近似得到，局部微觀模型如下：

式中：Iδ（xl）為以積分點(diǎn)xl為中心邊長為δ的立方體（三維）；水頭h（x）為宏觀水頭h在積分點(diǎn)xl處泰勒展開的線性近似。

微觀模型（5）的邊界條件是Dirichlet邊界條件，當(dāng)然微觀模型的邊界條件也可以是其他形式的邊界條件，如Neumann邊界條件或周期邊界條件等。

微觀模型（5）也可以通過有限元求解，這樣求解過程由兩個(gè)嵌套的有限元求解器組成。在求解微觀模型（5）的基礎(chǔ)上，（▽v·k0▽h0）（xl）通過下式計(jì)算得到：

將式（6）代入式（4）得：

3 隨機(jī)多尺度滲流模型

3.1 滲透系數(shù)隨機(jī)場(chǎng)

地層滲透系數(shù)在空間上本質(zhì)是異質(zhì)的、多尺度的，由于不可能獲取空間滲透系數(shù)完全的、詳細(xì)的信息，滲透系數(shù)需要用隨機(jī)理論來描述。一般認(rèn)為滲透系數(shù)的對(duì)數(shù)Y（x）＝ln（kε）為一弱二階平穩(wěn)隨機(jī)場(chǎng)，其均值在空間上為常數(shù)；協(xié)方差函數(shù)僅依賴于空間兩點(diǎn)之間的相對(duì)距離［14-15］。為了使用數(shù)值方法求解隨機(jī)滲流場(chǎng)問題，需對(duì)參數(shù)隨機(jī)場(chǎng)（如滲透系數(shù)）進(jìn)行離散化。這里通過Karhunen-Loμeve展開將參數(shù)隨機(jī)場(chǎng)Y（x）離散化［16-17］，并近似截取有限項(xiàng)：

式中μY（x）為參數(shù)隨機(jī)場(chǎng)Y（x）的均值；λi為隨機(jī)場(chǎng)Y（x）協(xié)方差函數(shù)的特征值；ψi為與特征值相對(duì)應(yīng)的正交特征函數(shù)，Yi為互不相關(guān)的隨機(jī)變量，如果Yi是高斯隨機(jī)變量，則他們是相互獨(dú)立的。

3.2 隨機(jī)多尺度滲流模型

如前所述，求解輸入?yún)?shù)為隨機(jī)變量的微分系統(tǒng)的傳統(tǒng)數(shù)值計(jì)算方法主要包括Monte Carlo法、矩方程法、攝動(dòng)法以及最近廣泛使用的廣義多項(xiàng)式混沌法。為了解決異質(zhì)地層隨機(jī)多尺度滲流問題，需在一個(gè)計(jì)算框架內(nèi)同時(shí)考慮隨機(jī)性和多尺度性，因此將廣義多項(xiàng)式混沌法與有限元異質(zhì)多尺度法結(jié)合起來是解決隨機(jī)多尺度滲流問題的一個(gè)有效途徑。

如果在計(jì)算中僅僅考慮滲透系數(shù)隨機(jī)場(chǎng)（也可以考慮其他類型的隨機(jī)場(chǎng)，如邊界條件隨機(jī)性、地層幾何形態(tài)隨機(jī)性等），考慮隨機(jī)輸入后的問題（1）變?yōu)槿缦码S機(jī)微分方程：

式中：Ω為隨機(jī)樣本空間。

根據(jù)展開式 （8），輸入隨機(jī)參數(shù)共有n個(gè)互不相關(guān)的隨機(jī)變量，令Y＝（Y1，Y2，…，Yn）為n維隨機(jī)變量，Ψ為隨機(jī)變量Y的樣本空間，ρ為Y的聯(lián)合分布密度，問題（9）的弱解形式為：

式中：y∈Ψ。

多種數(shù)值技術(shù)可以用來求解問題 （10）。在這些數(shù)值方法中，如同Monte Carlo方法一樣，配點(diǎn)法將問題（10）的概率空間和物理空間實(shí)現(xiàn)解耦使之變成確定性問題。這樣可以利用任何已有的數(shù)值方法來解決此確定性問題，可以有效利用已有的程序代碼。

問題（10）是建立在細(xì)尺度上的，這將面臨當(dāng)細(xì)尺度遠(yuǎn)小于研究區(qū)域D時(shí)無法求解的問題，因此需借助多尺度計(jì)算方法。這里將問題（10）的物理空間用有限元異質(zhì)多尺度法進(jìn)行計(jì)算。假定物理空間是尺度分離的，類比（7），問題 （10）的剛度矩陣用下面的數(shù)值積分公式來近似：

式中：Yi是n維隨機(jī)變量的Y一個(gè)實(shí)現(xiàn)。

從（11）中可以得知，滲透系數(shù)隨機(jī)場(chǎng)只需要在細(xì)尺度上局部區(qū)域進(jìn)行展開，隨機(jī)配點(diǎn)法的每一個(gè)實(shí)現(xiàn)實(shí)質(zhì)上是用來求解宏觀單元積分點(diǎn)處局部微觀模型問題。由于滲透系數(shù)隨機(jī)場(chǎng)只需要在微觀局部展開，式（8）只需截取少數(shù)幾項(xiàng)就可以達(dá)到所需精度。這一點(diǎn)對(duì)于隨機(jī)配點(diǎn)法的計(jì)算極為重要，因?yàn)椋?）中的n越大，隨機(jī)變量空間的維數(shù)越高，所需的配置點(diǎn)就越多，這將需要大量的計(jì)算資源。因此SHMFE能有效地緩解當(dāng)相關(guān)長度過小導(dǎo)致的概率空間維數(shù)過高的問題。SHMFE計(jì)算步驟概括如下：

（1）將滲透系數(shù)隨機(jī)場(chǎng)按式（8）展開，得到n維隨機(jī)變量Y，在此基礎(chǔ)上通過配點(diǎn)法生成Nc個(gè)概率空間積分點(diǎn)；

（2）在概率空間每個(gè)積分點(diǎn)處通過有限元異質(zhì)多尺度法求解確定性問題；

（3）按下面的公式計(jì)算隨機(jī)滲流場(chǎng)的統(tǒng)計(jì)量平均值和方差：

式中：ωi為積分權(quán)重。

4 算例

4.1 計(jì)算模型

算例選取的計(jì)算模型是一個(gè)具有不均勻多尺度滲透系數(shù)隨機(jī)場(chǎng)的二維穩(wěn)定滲流模型，用以說明SHMFE的計(jì)算過程及有效性。該計(jì)算模型來源于文獻(xiàn) ［17］，并對(duì)文獻(xiàn)中的模型加了適當(dāng)?shù)男薷摹?/p>

計(jì)算模型的物理區(qū)域?yàn)椋?1.5，0）×（-0.4，0.8）的矩形區(qū)域。模型物理區(qū)域左右兩側(cè)是定水頭邊界，總水頭分別為hD＝1（x＝-1.5）和hD＝0（x＝0）。上下側(cè)邊為零流量邊界。

滲透系數(shù)隨機(jī)場(chǎng)按式（8）近似離散。如果在宏觀尺度上進(jìn)行展開，對(duì)于協(xié)方差函數(shù)是離散指數(shù)型的隨機(jī)場(chǎng)，當(dāng)相關(guān)長度為0.1時(shí)，由于特征值衰減緩慢，式（8）需要截取100項(xiàng)才能達(dá)到規(guī)定的精度（截取項(xiàng)的特征值之和占總特征值之和的90%）。如果采用SHMFE方法，由于只需在局部微觀模型上展開，因而特征值衰減迅速，式（8）只需要截取少數(shù)幾項(xiàng)就能達(dá)到規(guī)定的精度（假定微觀模型的尺度與相關(guān)尺度相等，則只需4項(xiàng)即可）。因此，相對(duì)于隨機(jī)SFEM方法，SHMFE能有效地克服概率空間維過高的問題。SHMFE的另一個(gè)優(yōu)點(diǎn)是，計(jì)算僅需局部微觀模型區(qū)域統(tǒng)計(jì)數(shù)據(jù)，這與工程實(shí)踐一致，因?yàn)樵诠こ讨?，能獲取的只有局部數(shù)據(jù)，因此易獲得相對(duì)可靠的微觀區(qū)域統(tǒng)計(jì)模型。

算例僅為了說明SHMFE的有效性和可行性，而非實(shí)際工程，因此不具體考慮微觀尺度ε、細(xì)尺度下代表體積單元（即局部微觀模型尺度）以及統(tǒng)計(jì)模型相關(guān)尺度之間的關(guān)系。不失一般性，計(jì)算假定相關(guān)尺度各向同性且認(rèn)為微觀尺度ε、代表體積單元及相關(guān)尺度相等，并假定Y在局部微觀模型區(qū)域按式（13）展開［18］：

式中：ε為微觀尺度。

在算例中，考慮相關(guān)長度取1，0.1，0.05和0.01等4種計(jì)算方案。在這4種計(jì)算方案中，假定（13）中的uY（x）只與宏觀尺度有關(guān)，并令其等于sin（4＊x1＊x2）＋1。這4種方案均用SHMFE進(jìn)行模擬。為了檢驗(yàn)SHMFE可靠性，方案1除了采用SHMFE模擬外，還采用SFEM進(jìn)行模擬，并將2種計(jì)算方法的結(jié)果進(jìn)行對(duì)比。

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

圖1、圖2是方案1（即相關(guān)長度η＝1）采用兩種計(jì)算方法計(jì)算所得的水頭隨機(jī)場(chǎng)平均值和方差對(duì)比圖，圖1為平均值對(duì)比，圖2為方差對(duì)比，圖中虛線和實(shí)線分別為SFEM和SHMFE的計(jì)算結(jié)果。從圖中可知，兩種方法計(jì)算所得的結(jié)果是一致的，只在局部有細(xì)微的區(qū)別，這說明SHMFE方法可以獲得全尺度的SFEM法計(jì)算結(jié)果。由于4種方案隨機(jī)輸入?yún)?shù)（參見（13）式）均值相同，所以4種方案計(jì)算所得水頭隨機(jī)場(chǎng)的均值均相同（見圖1）。

圖1 相關(guān)尺度為1時(shí)水頭隨機(jī)場(chǎng)平均值比較（虛線和實(shí)線分別為SFEM及SHMFE的計(jì)算結(jié)果）

圖2 相關(guān)尺度為1時(shí)水頭隨機(jī)場(chǎng)方差比較（虛線和實(shí)線分別為SFEM及SHMFE的計(jì)算結(jié)果）

圖3 相關(guān)尺度為0.1時(shí)水頭隨機(jī)場(chǎng)方差

圖4 相關(guān)尺度為0.05時(shí)水頭隨機(jī)場(chǎng)方差

圖2到圖5是4種計(jì)算方案所得的水頭隨機(jī)場(chǎng)方差等值線圖。計(jì)算結(jié)果表明，水頭隨機(jī)場(chǎng)攝動(dòng)在研究的物理區(qū)域內(nèi)是不均勻的，攝動(dòng)主要集中在物理區(qū)域上側(cè)，并且受相關(guān)尺度大小影響顯著。當(dāng)相關(guān)尺度與區(qū)域尺度同數(shù)量級(jí)時(shí)，水頭隨機(jī)場(chǎng)攝動(dòng)集中在物理區(qū)域上側(cè)中間，并從中間向兩側(cè)遞減（見圖2）；隨著相關(guān)尺度減?。ㄏ嚓P(guān)尺度比物理區(qū)域尺度小1個(gè)數(shù)量級(jí)），攝動(dòng)集中區(qū)迅速向兩側(cè)擴(kuò)散，出現(xiàn)多個(gè)集中區(qū)（見圖3）；隨著相關(guān)尺度進(jìn)一步減?。ㄏ嚓P(guān)尺度比物理區(qū)域尺度小2個(gè)數(shù)量級(jí)），在物理區(qū)域的上側(cè)形成左右2個(gè)攝動(dòng)集中區(qū)，且逐漸趨于穩(wěn)定（見圖4、圖5），這說明相關(guān)尺度小于2個(gè)數(shù)量級(jí)后對(duì)計(jì)算結(jié)果基本無影響，事實(shí)上，此時(shí)（13）式中的隨機(jī)部分在宏觀尺度上可視為均質(zhì)的，即滿足尺度分離原則，水頭隨機(jī)場(chǎng)攝動(dòng)變化是微觀代表體積單元異質(zhì)所導(dǎo)致。

圖5 相關(guān)尺度為0.01時(shí)水頭隨機(jī)場(chǎng)方差

從上述算例可知，SHMFE只需在宏觀尺度上進(jìn)行局部微觀尺度計(jì)算就能很好地捕獲到微觀尺度異質(zhì)性對(duì)宏觀物理現(xiàn)象的影響，相對(duì)于SFEM而言具有更優(yōu)的計(jì)算效率。

5 結(jié)語

（1）為了模擬異質(zhì)多尺度地層穩(wěn)定滲流場(chǎng)，提出了一種結(jié)合有限元異質(zhì)多尺度法與隨機(jī)配點(diǎn)法的隨機(jī)異質(zhì)多尺度有限元法。算例表明，相對(duì)于SFEM而言，SHMFE能用更少的計(jì)算資源有效地捕獲異質(zhì)滲透系數(shù)隨機(jī)場(chǎng)對(duì)滲流隨機(jī)場(chǎng)的影響。該方法結(jié)合了HMM和SFEM的優(yōu)點(diǎn)，同時(shí)能有效地克服SFEM概率空間維過高的問題；另一方面，計(jì)算只需局部微觀模型數(shù)據(jù)，這與實(shí)際工程獲得的數(shù)據(jù)結(jié)果一致。

（2）本文沒有討論微觀尺度ε、代表體積單元和相關(guān)尺度的關(guān)系對(duì)計(jì)算結(jié)果的影響，也沒有討論尺度分離問題對(duì)計(jì)算結(jié)果的影響，這些問題需進(jìn)一步研究。

［1］Eaton TT.On the importance of geological heterogeneity for flow simulation［J］.Sedimentary Geology，2006，184（3-4）：187-201.

［2］P.E.Kloeden，E.Platen.Numerical Solution of Stochastic Differential Equations ［M］. Berlin： Springer-Verlag，1999.

［3］I.Karatzas，S.E.Shreve.Brownian Motion and Stochastic Calculus［M］.New York：Springer-Verlag，1988.

［4］G.S.Fishman.Monte Carlo：Concepts，algorithms，and applications［M］.New York：Springer-Verlag，1996.

［5］B.Oksendal.Stochastic Differential Equations：An Introduction with Applications ［M］. Berlin： Springer-Verlag，1998

［6］D.Xiu，G.E.Karniadakis.The Wiener-Askey polynomial chaos for stochastic differential equations［J］.SIAM J.Sci.Comput.，2002，24（2）：619-644

［7］Y.Efendiev，T.Y.Hou.Multiscale Finite Element：Theory and Applications［M］.New York：Springer，2009.

［8］T.Y.Hou，X.-H.Wu，Z.Cai.Convergence of a multiscale finite element method for elliptic problems with rapidly oscillating coefficients［J］.Mathematics of Computation 1999，68（227）：913-943

［9］T.J.R.Hughes，G.R.Feijo，L.Mazzei，J.-B.Quincy.The variational multiscale method—aparadigm for computational mechanics［J］.Computer Methods in Applied Mechanics and Engineering，1998，166（1-2）：3-24.

［10］E.Weinan，B.Engquist，X.Li，W.Ren，E.Vanden-Eijnden.Heterogeneous multiscale methods：a review［J］.Communications in Computational Physics，2007，2（3）：367-450.

［11］A.Abdulle，A.Nonnenmacher.A short and versatile finite element multiscale code for homogenization problems［J］.Computer Methods in Applied Mechanics and Engineering，2009，198（37-40）：2839-2859.

［12］F.Nobile，R.Tempone，C.G.Webster.A sparse grid stochastic collocation method for partial differential equations with random input data［J］.SIAM Journal on Numerical Analysis，2008，46（5）：2309-2345.

［13］D.Cioranescu，P.Donato.An Introduction to Homogenization［M］.Oxford：Oxford University Press，1999.

［14］李少龍，張家發(fā).非均質(zhì)土體中滲流場(chǎng)的非平穩(wěn)隨機(jī)模擬［J］.長江科學(xué)院院報(bào)，2009，10：49-53

［15］D.Xiu，J.S.Hesthaven.High-order collocation methods for differential equations with random inputs［J］.SIAM Journal on Scientific Computing，2005，27（3）：1118-1139.

［16］R.Ghanem，P.D.Spanos，Stochastic Finite Elements：A Spectral Approach ［M］. New York：Springer-Verlag，1991.

［17］Z.Chen，T.Y.Hou，A mixed multiscale finite element method for elliptic problems with oscillating coefficients［J］.Mathematics of Computation，2003，72（242）：541-576.