多尺度有限元法在地下水模拟中的应用(1)

　2004年7月

水 利 学 报SH UI LI X UE BAO 第7期

收稿日期:2003212215

基金项目:国家自然科学基金(40172082);国家自然科学重点基金(40335045);博士点基金(0284002)

作者简介:薛禹群(1931-),男,江苏无锡人,科学院院士,主要研究方向为地下水模拟。

文章编号:055929350(2004)0720007207多尺度有限元法在地下水模拟中的应用

薛禹群1,2,叶淑君1,谢春红3,张云1

(11南京大学地球科学系,江苏南京　210093;21南京大学污染控制与资源化研究国家重点实验室,江苏南京　210093;

31南京大学数学系,江苏南京　210093)

摘要:本文详细介绍了多尺度有限元法的基本原理,并将其应用于非均质多孔介质中的流动问题,对水文地质参数按函数连续变化、渐变和突变3种非均质多孔介质中的二维地下水稳定流、非稳定流分别用多尺度有限元法和传统有限元法进行了计算。计算结果的对比表明,多尺度有限元法比传统有限元法有效,既节省计算量又有较高的精度。

关键词:多尺度有限元法;非均质;多孔介质;地下水数值模拟

中图分类号:P641文献标识码:A

地下水含水系统大多是非均质的,应用传统有限元法或有限差分法解这类非均质介质中的水流问题,需要非常精细的剖分才能获得较为精确的解。若计算区面积很大,精细剖分产生的结点数过多,往往会超出普通计算机的容量,或者使计算时间过长,所以,人们一直致力于寻求某种方法[1～9],既可以减少剖分的单元数又可以保证计算精度,多尺度有限元法MsFE M (Multiscale Finite E lement Method )就可以很好实现这一目标。求解椭圆型问题的多尺度有限元方法由H ou 等[1,2]提出,它无需在小尺度上精确求解就能正确抓住解的大尺度特征,它通过基函数满足局部微分算子来实现。多尺度有限元法是数学家提出用来求解椭圆型问题的新的数值方法,在水文地质领域的应用国内外基本上未见报道。作者应用后发现,这种方法不仅对椭圆型问题(稳定流问题)有效,对抛物型问题(即非稳定流问题)也适用,在解决非均质多孔介质中的水流问题时具有既节省计算量又保证计算精度等优点。本文首先介绍多尺度有限元法的基本原理,然后对参数连续变化、参数渐变和参数突变3种非均质条件下的二维地下水稳定流、非稳定流数学模型分别应用多尺度有限元法和传统有限元法求解,将两种方法的计算结果与解析解或采用精细剖分的有限元法的数值解进行比较。H ou 等[1]

只谈及MsFE M 在参数连续变化的简单二维稳定流问题中的应用,没有涉及参数渐变、参数突变以及非稳定流问题。1　多尺度有限元法基本原理

众所周知,应用传统有限元法进行地下水模拟时通常包括以下4个步骤:(1)将研究区离散成若干个单元,假定每个单元上介质是均质;(2)采用多项式插值来表示单元内的水头分布;(3)应用迦辽金或里兹等方法建立有限元方程,再将它们集合形成整个研究区的代数方程组;(4)解该代数方程组得到各结点的水头。其中第2步常采用的是线性多项式,即单元内水头呈线性分布,这是对单元内真实流场的一种近似。如果单元内介质为非均质,则误差将增大,因此传统有限元法要求单元内介质是均质的。若含水介质非均质性明显,则必须将单元剖分得很小,以便保证每个单元的参数是常数。由此不难看出,

应用传统有限元法模拟大区域高度非均质介质中的水流问题会有困难,原因在于采用的基函数不反映渗流的物理实质,不满足渗流微分方程。如果基函数本身满足渗流微分方程,就不必再要求单元内参数是常数,这种介质的非均质性就可以通过解渗流定解问题反映到基函数上。例如对于稳定流,可以使基

函数满足相应的椭圆方程,这就是多尺度有限元法(MsFE M )的基本思路。数学家[1,2]已证明求解系数高

度振荡的椭圆型边值问题,采用根据椭圆型边值问题构建基函数的MsFE M 收敛、高效和精确。

MsFE M 是针对多尺度问题提出的解决办法,如多孔介质的非均质性就包含了很多尺度,这种多尺度的非均质性通常用介质渗透性的多尺度波动来表示。对于一个具有很多尺度的地下水流问题,在所有小尺度上用有限元法求解,即使有现代的大型计算机,也是困难的。所有小尺度直接求解能提供物理过程在所有小尺度下的定量信息。但从工程应用角度,掌握多尺度系统的宏观性质就足够了。MsFE M 就是一种在大尺度上求解,而又能将小尺度特征反映到大尺度上的方法。对于大区域的模拟问题,若非均质性明显,MsFE M 非常适用,它能克服上述传统有限元法在这种条件下应用的困难。MsFE M 中每个单元的参数可以变化,参数变化对水头分布的影响反映到单元基函数中,因此较粗的单元剖分就能刻画出研究区参数的变化和流场分布,所以计算量要小得多。MsFE M 的关键是如何构造反映单元内非均质介质中渗流场分布特征的基函数,这也是它与传统有限元法的本质区别。

111　MsFEM 的基函数　MsFE M 的基函数是通过使其满足相应定解问题来构建的,对于地下水稳定流问题,则要求基函数满足相应的椭圆型边值问题。

假设求解一个二维稳定流问题,其对应的二阶渗流偏微分方程(椭圆型方程)为

-Δ・K (x ,y )ΔH =f (x ,y )∈D

(1)式中:K (x ,y )为渗透系数;H 为水头;f 为源汇项;D 为研究区。

图1　单元Δijk

示意将研究区D 剖分为三角形单元,以任一单元Δijk (图1)为例来说明基函数的构

建。要求单元Δijk 顶点i 处基函数

-Δ・K (x ,y )Δ

于0。j 、k 点处的基函数

可以了。

上述方法构建的基函数一般得不到解析解。通常将单元剖分为若干个子单元,用有限元法等数值法求基函数的数值解。尽管求基函数需要额外的计算量,但若传统有限元法单元大小和上述子单元大小相同,则它的计算量远大于多尺度有限元法[1]

。基函数的精确性对水流模型的最终结果影响甚小[2]。在通过上述方法得到单元各结点的基函数后,单元内任一点的水头H 和传统有限元法相同就可以表示为:H (x ,y )=H (x i ,y i )

对于二维非稳定流问题,H ou 等[1]未谈及,作者尝试同样用椭圆型边值问题来构建基函数,效果也很好。

112　MsFEM 基函数的边界条件　虽然水流模型的最终结果对基函数的精确性不敏感,但如果求解基函数的椭圆型定解问题时采用不同的边界条件对最终结果的精度影响则很大

[1,2],也就是说最终结果对

基函数的边界条件敏感。这里所谓基函数的边界条件主要有两种类型:一种是线性边界条件,基函数沿边界线性变化。例如

线性边界条件与传统有限元法中基函数在边界上的变化类似,基函数不反映边界上参数的变化,只是根据边界上各点的坐标关系,从边界的一端到另一端从1线性地变化为0。振荡边界条件则基函数