裂隙介质广泛存在于地层结构中, 其内部裂隙结构复杂、分布随机、裂隙长度同基质材料差异很大, 具有强烈非均质和明显多尺度特征。为了解决这类问题, 各种多尺度、升尺度算法^{[1, 2]}被相继提出, 其中多尺度有限元法(MsFEM)^{[3, 4]}的应用最为广泛。近两年, 裂隙介质渗流的多尺度分析成为力学界的研究热点之一。Zhang等^[5]提出了基于离散裂隙模型的多尺度有限元法; Hajibeygi等^[6]提出了基于分层断裂模型的迭代多尺度有限体积法。由于离散裂隙模型和分层断裂模型的裂隙特征参数获取困难, 对于裂隙较多的地层计算量巨大, 为此, 本文采用双重介质模型求解裂隙介质的渗流问题, 并将多尺度有限元法引入到该模型中。

传统多尺度有限元法求解实际问题时多尺度基函数的构造和组装程序编写复杂, 运算耗费大量时间。为了解决这一问题, 本文对传统多尺度有限元法进行了改进, 将应变光滑技术引入到多尺度有限元法中, 提出了光滑多尺度有限元法(SMsFEM), 并编写了计算程序, 通过求解裂隙岩体渗流问题验证了该算法的正确性和高效性。

双重介质模型最初由前苏联学者Barenblatt^[7]于1960年提出, 是目前在水文地质领域应用最广的裂隙介质渗流模型。该模型将裂隙和孔隙系统地处理为两个平行的连续体, 两者之间存在水交换, 其值取决于两者的水头差。控制方程如下:

式中:下标f表示裂隙; m表示孔隙; Ω 为求解域; h为压力水头; c为贮水系数; k为渗透系数, 一般是随机非均质的; α _w 为一阶质量传递系数, w为孔隙度; Γ ₁为已知水头边界; 为边界上的水头值; Γ ₂ 为已知流量边界; q为边界上的流量值; n_x、n_y 为边界外法线的方向余弦。

在双重介质模型中, 裂隙和孔隙的渗透系数均存在强烈的非均质性, 并且裂隙和孔隙内流体的流动存在耦合关系。本文分别构建裂隙水头和孔隙水头的多尺度基函数, 使其满足各自的椭圆微分方程。采用初始时刻的全局细尺度解确定多尺度基函数的边界条件, 使基函数在捕捉局部非均质性的同时可以捕捉到裂隙和孔隙内流体流动的耦合信息。并采用应变光滑技术构造宏观矩阵, 进一步提高算法的精度和效率。

选取1.5 m× 1.5 m的裂隙岩体作为研究对象进行分析。首先对模型进行空间离散, 将研究区域划分为6× 6的粗网格和30× 30的细网格, 如图1所示。在每个细网格内, 材料参数是均质的, 以便使用常规有限元法在细网格上进行求解。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 模型粗细网格的剖分Fig.1 Fine and coarse scale meshes

为了体现裂隙和孔隙介质的非均质性, 将细网格单元的渗透率参数取为不同值, 图2为裂隙介质的渗透率分布图。孔隙的渗透率取为坐标 的函数(k_m=0.1x), 其他参数的选取如下:w_f=0.1; w_m=0.9; α _w=0.33 d^-1m^-1; c_f=0.01 m^-1; c_m=0.1 m^-1; Δ t=0.005 d; θ =2/3。初始水头值取为2.0 m, 左右边界水头值分别为2.5、2.0 m, 上下边界为不透水边界。取稳态条件下的细尺度有限元解来确定多尺度基函数的边界条件。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 kf分布图Fig.2 Distribution of kf

分别采用光滑多尺度有限元(SMsFEM)、多尺度有限元(MsFEM)和在细网格上的常规有限元(FEM-F)三种算法求解裂隙介质渗流问题, 对计算的结果和耗时进行对比分析。

选取0.05 d的裂隙和孔隙水头的SMsFEM和FEM-F的结果绘制等值线图进行比较, 如图3所示。从图中可以看出:采用常规有限元在细网格上的计算结果与光滑多尺度有限元在粗网格上的计算结果基本一致。可见, 采用有限全局信息的光滑多尺度有限元法可以有效地捕捉材料的微观信息, 在较粗的网格上得到较精确的结果, 验证了SMsFEM的正确性。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 水头等值线图Fig.3 Contour graphs of hydraulic head

图4为SMsFEM、MsFEM和FEM-F三种算法在t=0.05 d和t=0.1 d, 1~7点处(如图1所示)求解的裂隙及孔隙水头值对比曲线。从图中可以看出:SMsFEM和MsFEM的值与FEM-F的值吻合得很好, 其中SMsFEM的结果更接近FEM-F的值。主要是由于SMsFEM算法在组装宏观矩阵时将面积分转化为线积分, 克服了传统有限元计算刚度过硬的缺点, 由此提高了计算结果的精度。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 SMsFEM、MsFEM和FEM-F求解的水头值对比曲线Fig.4 Hydraulic head of SMsFEM、MsFEM and FEM-F

用SMsFEM、MsFEM两种算法计算宏观矩阵时计算耗时分别为14.3512 s和19.1325 s, 采用SMsFEM的计算时间比MsFEM快25%。对于网格数目众多的实际工程问题, 此种优势会更加明显。可见, 将应变光滑技术引入到多尺度有限元中可以有效提高计算效率。

表1为SMsFEM、FEM两种算法的计算耗时对比。通过对表中数据进行处理可以得到, 在t=0.5 d时SMsFEM的计算耗时仅为FEM-F的17.5%; 在t=2.5 d时SMsFEM的计算耗时仅为FEM-F的7.4%; 在t=5 d时SMsFEM的计算耗时仅为FEM-F的6.1%。可见, 采用SMsFEM算法在粗网格上计算非均质裂隙介质的渗流问题大大节省了计算时间, 并且随着运算时间的增加, SMsFEM同FEM-F计算耗时的比值越来越小, SMsFEM算法的优势越来越明显。

表1

Table 1

表1(Table 1)

表1 SMsFEM和FEM计算耗时对比 Table 1 Comparison of CPU time between SMsFEM and FEM t=0.5 d t=2.5 d t=5 d SMsFEM FEM SMsFEM FEM SMsFEM FEM CPU 时间/s 24.31 138.76 52.12 702.90 86.14 1404.28

表1 SMsFEM和FEM计算耗时对比 Table 1 Comparison of CPU time between SMsFEM and FEM

(1)采用光滑多尺度有限元法求解的结果同有限元法在细网格上求解的结果基本一致, 并且大大节省了计算的时间。

(2)本文提出的光滑多尺度有限元法比常规的多尺度有限元法计算耗时更少, 精度更高。

(3)利用有限全局信息构造的多尺度基函数可以更好地捕捉材料的非均质性及裂隙和孔隙介质内流体流动的耦合关系。

本文提出的光滑多尺度有限元法可以解决几乎所有常规多尺度有限元法可以解决的问题。比如:多孔介质内单相流和多相流的流动问题、复合材料的热传导问题等, 可见该方法具有广阔的应用前景。