求解非均质多孔介质中随机水流问题的多尺度有限元降基方法

黄梦杰, 贺新光

(1.湖南师范大学 资源与环境科学学院， 湖南 长沙 410081；2.湖南师范大学 地理空间大数据挖掘与应用湖南省重点实验室， 湖南 长沙 410081)

1 研究背景

地下多孔介质的水力特性通常是高度非均质的，且在多种不同的空间尺度上发生变化。由于对介质水力特性空间分布的不完全认知和测量误差，在实际应用中，通常视水力特性参数的空间分布为空间随机场，以致对多孔介质中水流运动的模拟预测总是存在一定的不确定性。为了处理这种不确定性，有必要进行地下多孔介质中水流过程的随机模拟[1-2]。然而，介质水力特性的多尺度非均质性与不确定性导致发展高效模拟非均质多孔介质中随机多尺度水流模型的数值算法极具挑战性[3]。如果应用传统的数值方法模拟这样的随机多尺度模型，其计算过程是非常耗时的，有时甚至是不可行的[4]。正如Xu[5]和Chung 等[6]所指出：以不确定性数据为特征的随机多尺度问题，目前依然是公开的、极具挑战性的研究领域之一。因此，最近十几年来，在计算科学与工程领域，发展模拟这类随机多尺度数学物理模型的随机多尺度数值方法的兴趣一直稳定地增长[5-12]。

模拟地下多孔介质中水流运动的挑战之一是多孔介质特性的空间多尺度性。这类问题因分辨最细空间尺度信息所需的花费，致使其在空间上直接的细尺度分辨模拟在计算上是非常复杂而耗时的，也是不切实际的。然而，直接忽略细尺度信息可能导致较大的误差。为了切实有效地处理大尺度区域上的多尺度问题，各种多尺度数值方法已被提出且受到了越来越多的重视，如：多尺度有限元方法[13](MsFEM)、变分多尺度方法[14]和非均质多尺度方法[15]等。其中，MsFEM是多尺度数值方法的先驱之一，且许多别的多尺度方法共享了它的相似性[4,16]。一些研究已成功应用MsFEM数值模拟地下非均质多孔介质中的水流问题[3,16-19]。然而，当MsFEM与随机方法直接结合用于模拟随机多尺度水流模型时，势必导致需要计算大量的依赖于随机参数样本的多尺度有限元基函数，这是十分不利的[3]。因此，为了高效模拟随机多尺度地下水流问题，研究构造独立于随机参数取样的多尺度基函数是十分必要。

除了多尺度特性，模型输入量的不确定性也为多孔介质中的水流模拟带来巨大的挑战。为了尽量减少计算成本而又能维持良好的计算精度，近十几年来，不少研究者基于Galerkin投影提出了一些求解随机模型的降基方法(RBM)。RBM是一种有效的随机模型降阶方法，已成功应用于各种研究与工程领域[20-24]。RBM依降基函数构造方法的不同，可分为本征正交分解[25-26]( POD)和贪婪算法[20, 22](Greedy)。RBM的基本思想是：基于一组随机抽取的参数样本，构造一组独立于参数样本的降基函数，然后对于任意新的随机取样，投影参数模型到降基函数空间而生成一个降阶模型(ROM)。为加速在线计算，在降基方法中，离线-在线计算分离是强制性的。在离线阶段，大量的基础计算包括物理空间的有限元离散、降基函数的构造和全局离散矩阵的组装等。在在线阶段，通过任意多次抽取参数样本，求解参数化的降阶模型，感兴趣的输出量可以被计算。然而，对于随机多尺度模型，假如应用传统的有限元方法构造全局的降基函数，则离线阶段的计算因需分辨最细空间尺度信息而十分耗时。为降低离线阶段的计算成本，最新的研究[4,27]建议：应用MsFEM计算局部辅助函数，并进一步应用POD生成多尺度降基函数空间以构造随机多尺度问题的降阶模型。

本文遵循最新研究[4]为处理椭圆型随机多尺度问题的主要思想，通过有效耦合多尺度有限元和降基方法而发展一种高效求解非均质多孔介质中的随机多尺度水流问题(抛物型)的多尺度有限元降基方法(RMsBM)。而且，为了加速在线计算，不同于文献[4]，本文采用一种矩阵离散的经验插值方法[28](MDEIM)仿射分解参数依赖的代数系统以实现水流降阶模型的离线-在线计算分离。Negri 等[28]指出：MDEIM方法以一种非浸润、高效而纯代数的方式处理非仿射的参数化离散系统。最后，为表明所提出的RMsBM的性能与效率，本文执行了几个有代表性的数值试验。结果表明：本文为求解非均质多孔介质中的随机多尺度饱和水流问题提出了一种高效的数值模拟方法。

2 模型问题

本节介绍非均质随机多孔介质中的承压水流控制方程及其全阶模型和降阶模型。

2.1 水流控制方程

令D⊂R2为二维有界物理域，(0,T]为时间域，D×(0,T]为笛卡儿积，μ∈Ω⊂Rp为p维随机参数向量。非均质随机多孔介质中的非稳定饱和水流可用含参数的抛物型偏微分方程描述如下。

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为空间坐标；ψ为水力水头；Ss为比贮水系数；K为渗透系数；R为源/汇项，n为单位外法向量；ΓD和ΓN分别为区域的Dirichlet和Neumann边界，且ΓD∩ΓN=Ø。

对于非均质多孔介质，K一般在随机参数空间振荡且在物理空间高度非均质。然而，在任何情况下，Ss的变化幅度均小于lnK[29]，因此，本研究假设Ss为区域D内的一个确定性函数。对于含参数的水流问题(1)，解的统计性质如水头均值和方差是令人感兴趣的。为此，可首先利用Monte Carlo或稀网格随机配点方法对随机参数空间Ω进行抽样，然后，对于任意给定的参数μ∈Ω，应用有限元等确定性方法求解参数化的水流问题(1)，并可最终获得解的统计性质。然而，对于一个给定的参数μ，为捕捉细尺度非均质性，仍需要在物理空间执行高分辨模拟。这样，如果直接在一个细网格上对流动问题执行全阶模拟，则计算成本是非常高的。为了提高计算效率，本文发展一种基于多尺度有限元基函数的降基方法求解参数方程(1)。

2.2 全阶模型

(5)

(6)

(7)

式中：矩阵CNf和ANf(μ)在第(i,j)处的元素分别为(CNf)ij和(ANf)ij；向量b的第i个元素为(bNf)i。本文称方程式(7)为水流问题(1)的全阶模型。为了获得全离散的解，可以采用Euler或其他方法离散半离散的方程式(7)。本文应用常用的Crank-Nicolson(CN)方法构造方程式(7)的全离散格式。关于CN全离散格式的更多细节，请参阅He等[3]或其他相关文献。

2.3 降阶模型

(8)

3 多尺度有限元降基方法

3.1 多尺度有限元基函数的构造

(9)

(10)

式中：gi为定义在κ的边界∂κ上且与节点i相关的边界条件；d为粗单元的节点总数。

3.2 多尺度降基函数的生成

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

(11)

(12)

(13)

3.3 离线-在线计算

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

(14)

降阶模型(8)全离散后，可得一个Nr元的线性代数方程组。根据方程(14),方程(8)的降阶矩阵和向量均涉及φirm(i=1,2，…,Nr)的内积计算。因降基函数φirm属于有限元空间Xh，故其可用Nf个标准的有限元基函数ξ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)所定义。因为矩阵Z和CNf以及向量bNf均独立于随机参数μ，所以，降阶矩阵CNr和向量bNr均只需在离线阶段计算一次。然而，为了得到降阶模型的解ΨNr(μ)，依赖于参数μ的全阶矩阵ANf(μ)需重新计算而导致巨大的计算量。因此，ANr(μ)的有效估计对于实现降阶模型的离线-在线计算分离是至关重要的。如果矩阵ANf(μ)可用依赖于μ的系数加权的常数矩阵的仿射组合表示，则加权和的每一项可离线投影到多尺度降基函数空间。假设ANf(μ)可仿射表示为：

(17)

则有：

(18)

4 数值模拟

本节执行几个数值算例以表明所提出的多尺度有限元降基方法求解非均质多孔介质中随机饱和水流问题的性能。为此，考虑非均质含水层对水位突变的响应，即：假如在时刻t=0时，突然将含水层边界上的水位从15 m下降到12 m，并模拟一段时间内水头随时间的变化。假设含水层的比贮水系数和厚度分别为8×10-4/m和10 m，则随机水流问题(1)的初始条件和边界条件分别为：

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

(19)

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

(20)

并假设渗透系数对数Y=ln(K(x,μ))是一个二阶平稳的高斯随机场，其两点指数协方差函数cov[Y]为：

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

(21)

式中： (xi,yi)为空间坐标；σ2为随机场Y的方差；lx和ly分别为x和t方向的相关长度。任给一个协方差函数，其随机场Y可由截断的KLE展式生成[3]：

(22)

式中：随机向量μ:=(μ1,μ2,…,μNk)∈RNk,μi(i=1,2，…,Nk)为独立同分布的标准正态随机变量，λi和bi(x)分别为协方差函数的特征值和特征函数。此时，水流方程(1)则是一个Nk维的随机参数问题。因输入场K(x,μ)=exp(Y(x,μ))关于参数μ是非仿射的，为了实现水流问题(1)的ROM的离线-在线计算分离，本研究应用MDEIM对非仿射的参数矩阵执行仿射分解。

为了不失一般性，本节所有的数值试验均考虑同一空间尺度的物理区域D。令区域D为1个200 m × 200 m的矩形。然后，1个一致的有限元网格将区域D划分成Nx×Ny个子矩形，并称之为Nx×Ny的网格。时间离散步长固定为Δt=5 min。为表明多尺度降基方法的性能，本节在粗网格上应用所提出的RMsBM和在细网格上应用标准的FEM分别求解水流问题(1)的ROM与FOM模型，并对所获得的解进行比较。下列的数值算例中，标准的FEM采用128×128的细网格求解水流问题的FOM模型，并将其解作为参考解。为了评估粗网格尺寸对RMsBM的影响，下面的例子将考虑3种不同的粗网格，即8×8、16×16和32×32的粗网格。为了方便，将RMsBM方法获得的ROM解称之为RMsBM解。为估计RMsBM解的精度，本研究使用如下定义的相对L2误差：

(23)

式中：Nf为细网格上的节点总数;ψr和ψf分别为水流问题的RMsBM解和参考解。

4.1 局部基函数数目的影响

本小节旨在讨论局部多尺度基函数的数目对所提出的RMsBM性能的影响。假定随机场Y的方差σ2=2.0，其几何平均等于-6.125，且相关结构是各向异性的，并具有相关长度Lx=40 m和Ly=30 m。此数值算例中，使用KLE展式，取Nk=25项生成水力参数传导场K(x,μ)=exp(Y(x,μ)。此时，水流方程(1)则是一个25维的随机参数问题，并在一个16×16的粗网格上求解该问题的RMsBM解。为此，在RMsBM的离线阶段，首先从随机参数空间选取350个25维的参数样本，并应用MsFEM在粗网格上分别计算所选样本的局部snapshot函数，然后，对所选样本构造的snapshot函数分别在其相应的每一个局部领域执行本征正交分解，最终获得局部多尺度降基函数。在RMsBM的在线阶段，则使用已构造的局部多尺度降基函数求解任意给定的随机参数样本所确定的水流问题的RMsBM解。

为了评估局部基函数数目对RMsBM性能的影响，本算例分别考虑粗网格的每个粗节点处均使用同一数目N=1，2，…，5的局部多尺度降基函数的情形。应注意到：128×128细网格上的FOM离散方程组的自由度为Nf=16 641，然而，当N=1,2,…,5时，16×16粗网格上的ROM离散方程组的自由度仅分别为Nr= 289、578、869、1 156和1 445。显然，RMsBM的ROM可以提供比标准的FOM明显更快的在线计算，且RMsBM的在线计算随局部基函数数目的减少而加速。为了评估RMsBM解的近似精度，本例从参数空间中再随机抽取2000个新的参数样本，并使用RMsBM在离线阶段已构造的多尺度降基函数空间分别求解每一新参数样本所确定的水流问题(1)的RMsBM解。然后，分别计算这2000个样本的RMsBM解和参考解的均值和方差。图1展示了在不同时刻t=8、14和20 h，每个粗节点处使用不同数目局部多尺度降基函数的RMsBM解的均值和方差的相对L2误差。从图1可知，当局部降基函数的数目从1开始增加时，3个不同时刻的水头均值和方差的相对误差均快速衰减，然后，当局部基函数数目从3增至5时，所有的相对误差都一直保持很小的值。此外，图1也表明，每个粗网格节点处仅具有3个局部基函数的RMsBM可为参考解的均值和方差提供良好的近似逼近。说明在该情况下，当局部降基函数数目N=3时，RMsBM可以在近似精度和在线计算成本之间取得良好的平衡。具体地，对于一个参数样本，为获得其在时刻t=20 h的解，具有3个局部基函数的RMsBM和标准的FEM 的平均在线CPU时间分别是26和53 s。这表明：相较于标准的FEM，对于这种情况，具有3个局部基函数的RMsBM可节省大约51%的在线CPU时间。为了进一步比较，图2和3分别描绘了在时刻t=8 h，由具有1个和3个局部基函数的RMsBM和标准的FEM所获得的水头均值和方差的空间分布。从图2和3可以看出：由具有不同数目的局部基函数的RMsBM获得的水头均值均能很好地匹配参考的水头均值，且由具有3个局部基函数的RMsBM得到的水头方差也与参考的方差具有良好的一致性。然而，由每个粗网格节点处仅具有1个局部基函数的RMsBM计算的水头方差没能很好地匹配参考的方差。这表明：使用适当数目的局部多尺度降基函数对RMsBM方法的效率与精度非常重要。为了直观展示不同参数样本的个体误差，图4提供了前300个样本在时刻t=8 h、每个粗节点处具有不同数目局部基函数的RMsBM解的相对L2误差。从图4可以观察到，具有3个局部基函数的RMsBM为这300个样本提供了最准确的解，而仅具有1个局部基函数的RMsBM则给予了最不精确的解。此外，图4也表明，相比于具有1个或2个局部基函数的RMsBM，具有3个局部基函数的RMsBM所获得解的相对误差对参数样本的选择是不敏感的。这些结果表明：在16×16粗网格上，每个节点处具有3个局部降基函数的RMsBM可以为128×128细网格上的参数水流问题的FOM模型提供稳定而准确的ROM近似解。

4.2 粗网格尺寸的影响

本节使用8×8、16×16和32×32粗网格考察不同粗网格尺寸对RMsBM性能的影响。为此，假定渗透系数随机场K(x,μ)的所有参数与4.1节相同。对不同的粗网格，RMsBM在其离线阶段构造多尺度降基函数之后，在其在线阶段，分别求解随机选取的2 000个参数样本中每一样本所确定的水流问题(1)的ROM解。图5展示了使用不同数目局部基函数的RMsBM在3个不同粗网格上在时刻t=8 h获得的RMsBM解的均值和方差的相对误差。从图5中可观察到误差的两个主要特征：一是对于1个给定的粗网格，RMsBM解的误差随局部基函数数目的增加而单调减少，特别是8×8和16×16粗网格，当每节点处的局部基函数由1个增至3个时，RMsBM解的误差显著减少；二是随着粗网格的细化，每个节点具有固定数目局部降基函数(1、2或3个)的RMsBM解的近似精度普遍提高。然而，每个节点具有3个局部基函数时，RMsBM在32×32粗网格上并没有提供比其在16×16粗网格上更高精度的解。这表明：在这种情况下，RMsBM存在一个最优的粗网格(16×16)。此外，还观察到：RMsBM在8×8粗网格上每节点具有6个局部基函数的解和其在16×16粗网格上具有3个局部基函数的解具有几乎一样的近似精度。

图1 具有不同数目局部基函数的RMsBM解在3个不同时刻的均值和方差的相对L2误差

图2 参考解与具有1个和3个局部基函数的RMsBM解在时刻t=8 h的水头均值对比

图3 参考解与具有1个和3个局部基函数的RMsBM解在时刻t=8 h的水头方差对比

图4 具有不同数目局部基函数的RMsBM对前300个参数样本在时刻t=8 h解的相对L2误差对比

图6给出了8×8粗网格上每节点具有3个和6个局部基函数的RMsBM解和参考解在时刻t=8 h的方差的对比。从图6可看出：对于8×8粗网格，当每节点处的局部基函数从3个增至6个时，RMsBM解的方差可以很好地逼近参考方差。然而，随每节点局部基函数数目的增加，RMsBM的在线阶段则需要更多的计算成本。显然，8×8、16×16和32×32粗网格的每个粗单元分别包含16×16、8×8和4×4的细网格。从对RMsBM的描述可以看出：在离线阶段构造局部snapshot函数时，局部离散问题的自由度随着粗网格的细化而减少，但其局部问题的总数却随之增加。因此，本文建议：在这种情况下，为了取得解的近似精度和计算成本之间的良好平衡，可使用16×16粗网格和每节点3个局部基函数的RMsBM有效地降阶水流问题的全阶模型。

图5 具有不同数目局部基函数的RMsBM在3个不同粗网格上解的均值和方差的相对误差

图6 参考解与8×8粗网格上具有3个和6个局部基函数的RMsBM解的方差对比

4.3 空间变异性的影响

本节探讨渗透系数场的空间变异性(σ2)对所提出的RMsBM性能的影响。为此，考虑具有3个不同方差 1.5、2.0和2.5的随机参数输入场，并假设其余的相关参数与4.1节中的相同。基于4.2节的讨论，本节仅考虑RMsBM在16×16粗网格上执行模型降阶。同样地，对每个输入场随机选取的2 000个样本，分别求解每个样本所确定的水流问题(1)的RMsBM解和参考解。显然，水流模型输出水头的空间变异性随其输入参数场σ2的增大而增强。当σ2较大时，期望所构造的多尺度降基函数也能在粗网格上有效捕捉输入场的高空间变异性。为此，图7展示了具有不同数目局部基函数的RMsBM在16×16粗网格上为3个不同随机参数输入场获得的解在时刻t=8 h的均值和方差的相对误差。图7表明：尽管输入场方差的不同会导致输出水头空间变异性的不同，但RMsBM对3个不同的σ2给出了几乎相同的水头均值和方差的相对误差。这里，特别注释：当σ2=2.5时，随机渗透系数场的许多实现达7个数量级的变化。

图7 不同数目局部基函数的RMsBM对不同随机场的解在时刻t=8 h的均值和方差的相对误差

这表明：所提出的RMsBM对具有高空间变异性输入场的随机水流问题的降阶具有很强的稳健性。

图8对3个不同的空间变异，给出了参考解与使用3个局部基函数的RMsBM解在时刻t=8 h和截面y=100 m处的均值和方差的对比。从图8可观察到：具有3个局部基函数的RMsBM解的均值和方差均能很好地匹配所有这些情况的参考均值和方差。这表明：在适当的粗网格上使用合适数目的局部基函数的RMsBM能有效地捕获具有高空间变异性的细尺度异质性对粗尺度解的影响。

图8 不同随机场的参考解与具有3个局部基函数的RMsBM解在时刻t=8 h和截面y=100 m处的均值和方差的对比

5 结 论

本文通过有效耦合多尺度有限元和降基方法，提出并讨论了一种求解非均质多孔介质中的随机参数饱和水流问题的多尺度有限元降基方法。该方法的有效性在于构造一组独立于随机参数取样的多尺度有限元降基函数和生成一个降阶多尺度模型，并应用矩阵离散的经验插值方法仿射分解非仿射的随机参数问题的离散系统以加速降阶模型的在线计算。所提出的RMsBM允许其计算过程执行一个离线-在线计算分离，即在RMsBM离线阶段，构造多尺度有限元降基函数，而在在线阶段，则使用已构造的多尺度降基函数求解任意参数样本所确定的水流问题的解。RMsBM的离线计算可能较为耗时，但其在线计算是快速有效的，而快速有效的在线计算对预测任意参数样本所对应的模型输出量和评估随机输入量对输出量的影响是十分重要的。数值算例表明：

(1)RMsBM通过选择一个合适的粗网格和一个最优数目的降基函数，可实现随机参数水流问题求解在近似计算精度和在线计算成本之间的良好平衡。

(2)RMsBM对具有高空间变异性输入场的随机水流问题的降阶模拟提供了一个稳健而精确的替代模型。

总之，所提出的RMsBM无需在全局细网格上分辨所有的随机小尺度信息，即可在适当的粗网格上应用合适数目的局部降基函数，有效捕获随机输入场的细尺度非均质性对水流模型输出场的影响，从而大大降低了计算的复杂性，节省了计算成本。因此，方法可为大尺度区域地下水污染状况的实时监测分析和放射性废物地下储放的安全性评价提供一种有力的数值模拟工具。