海水入侵问题的间断有限元方法

【摘要】： 海水入侵是在许多沿海国家和地区都普遍存在的问题,它是指海水侵入地下淡水蓄水层的复杂运动,一般是由自然水环境的变化和社会经济发展所造成的。由于淡水的密度小于海水的密度,所以在地下蓄水层中,淡水总是存在于海水的上方,并且在咸淡水界面上达到一种动态的压力平衡。沿海国家和地区为了满足人类生活和经济发展的需要,往往会过度地开采地下水(淡水),这就会引发海水入侵问题。当淡水的开采量超过了它的补给量,海岸带附近的淡水水头(水位)就会随之下降,进而导致淡水在咸淡水界面处的产生压力减小,压力平衡遭到破坏。在压力作用下,海水会侵染淡水层,以达到新的平衡。当海水侵染到淡水开采井时,开采出的水就不再适合饮用和生产,这将给沿海地区的人类生活和社会经济发展带来严重的危害。目前,全世界有几十个国家和地区的几百个地方都发现了海水入侵问题,因此,这一问题受到了国际社会的广泛关注,许多国家都在积极开展海水入侵问题的研究和治理工作。 近些年来,对海水入侵问题的数学模型及数值解法的研究已经取得了一些成果。基于咸淡水界面为两种流体混合的过渡带的假设,Henry[51]得到了稳定流在匀质介质、简化边界条件情况下的解析解。1970年,Pinder[67]最早提出了海水入侵过渡带模型,并提出了Henry模型的有限元数值解法。Segol等人[77,78]发展了地下水位和浓度相互依赖的剖面二维有限元模型。Huyakorn等[54]提出了依赖浓度的地下水流方程和溶质运移方程,研究了三维情形的有限元模型。但是,上述海水入侵的研究工作都是在一些具体的假设条件得到的,并不能真实的反映海水入侵过程。90年代初,薛禹群等[85,86]研制我国第一个三维有限元海水入侵模型。他们考虑了速度随盐分变化的情形,但是没有给出对防止海水入侵的各类工程的数值预测；参见[47]及其引文。自上世纪末,袁益让,梁栋和芮洪兴教授深入研究了海水入侵数学模型及其强对流占优和非线性的特点,提出了处理海水入侵模型的特征有限元[88]、迎风分数步有限差分[89,90,91,92]等方法,对海水入侵防治工程实施后咸淡水变化运移的真实过程进行了定量模拟。类似于文献[39,38,44]中处理油藏问题的方法,这些方法克服了普通差分或有限元方法容易产生数值振荡的缺点,但也存在缺点：迎风方法引入的额外的扩散会影响结果的精度,而特征方法不便于处理非线性问题[45]。另外,对于出现间断解的情形,上述方法不能很好的解决。本文在上述研究成果的基础上,充分考虑了海水入侵模型中相互耦合的两个方程的特性,选用混合元和间断有限元相结合的方法来处理这一问题,并相应地提出半离散和全离散格式。 1973年,Rccd和Hill[[69]在处理中子输运方程——不依赖时间的线性双曲型方程时,提出了第一个间断有限元方法。自那时起,求解双曲型方程的间断伽略金方法得到了广泛的关注和积极的发展。间断的思想与其他数值技巧相结合,发展出了许多不同的数值方法,如Johnson[55]等人将迎风的思想融入到间断伽略金方法中,用来处理双曲型方程；Cockburn等人将高精度有限差分和有限体积方法的技术(如迎风、数值流通量、斜率限制器等,参见[57])与伽略金方法相结合,提出了Runge-Kutta间断伽略金方法[24,23,22,20,25],这也是双曲型方程间断伽略金方法的主要进展[21]。为了能够处理椭圆型和对流扩散型方程,Cockburn等人又提出了局部间断伽略金方法[26,15]。同样在70年代,求解椭圆型和抛物型方程的间断伽略金有限元方法(也被称为内罚方法)也被独立地提出来,并得到了深入的研究和发展[37,7,84,2]。近几年来,许多新型的内罚方法相继被提出和研究,包括对称内罚伽略金方法[2],Baumann-Oden方法[9,64],非对称内罚伽略金方法[74,72,75],以及不完全内罚伽略金方法[34]等；更多细节可参见综述性文献[3,14,4,70]。 由于与连续的伽略金有限元方法相比,间断伽略金方法具有局部守恒,易于实现hp自适应,可以更好的逼近间断解等优点[40],所以它们被广泛的应用于许多实际问题。尤其以在计算流体力学领域的应用最为广泛,涉及不可压缩[80,81]和可压缩[16,30,31]混溶驱动问题,Navier-Stokes问题[8,48,49,50,65],无粘性可压流体[35],两相流[43,41,42],气体动力学[17],半导体等器件的数值模拟[17,60],多孔介质流体[82,73]和多孔弹性介质流体[66]等。另外,它们还被应用于弹性力学[61,71]等其他物理问题中。注意到上述实际问题的模型大多是对流扩散方程,所以处理此类问题的数值格式往往是通过将两类间断伽略金方法相结合得到的,即采用内罚方法处理扩散项,而用双曲型间断伽略金方法处理对流项。不仅如此,在处理对流扩散问题,特别是对流占优的情形时,往往在间断伽略金方法中引入一些有限差分和有限元方法的技巧,如迎风[55,79,52,81,82,31,70],斜率限制器[70]和数值流通量[35,36]等,以取得更好的逼近效果。 混合元方法自出现以来,就得到了广泛和深入地研究[11,68,62,63],主要用于求解椭圆型和抛物型偏微分方程[56,13,12]。混合元方法首先通过引入一个新的变量(一般是未知函数导数的某种线性组合),将最初的方程化成混合形式,然后对混合变分形式使用有限元方法进行求解。这样可以将未知函数和新引入的变量同时求出,并且对这两个函数的逼近都可以达到最优阶。基于这一优点,早在上世纪八十年代,混合元方法就被用于与其他方法相结合,共同求解耦合的偏微分方程系统,如多孔介质混溶驱动模型[38,44]。近几年来,随着间断伽略金方法的发展,这种新型的方法也经常与混合元联系在一起使用,来处理一些复杂的耦合系统,例如不可压缩[80,81]和可压缩[16,30,31]混溶驱动问题,多孔弹性流体问题[66]等。由于海水入侵的数学模型是一个由一个抛物方程和一个对流扩散方程组成的耦合系统,本论文也采用类似于上述的想法来处理这一问题,即利用混合元方法求解抛物方程,用间断伽略金方法求解对流扩散方程。 在芮洪兴教授的精心指导下,作者完成博士论文的写作。在本论文中,作者充分研究了袁益让等人在文献[88,89,90,91,92]中考虑的海水入侵数学模型,提出了几种处理该模型问题的数值格式。本文所采用的数值方法是混合元方法和间断伽略金有限元方法,又由于此海水入侵的数学模型是耦合的偏微分方程系统,所以本文中还用到了迭代方法。应用这些方法,作者提出了连续时间格式和隐式离散时间格式。对这些格式,我们在理论上分析了格式的收敛性,得到了先验的误差估计。为了提高计算的速度和效率,作者改进了全隐离散时间格式,提出了一个半隐的离散格式,并在最后给出了并行求解此格式的区域分解方法。本文共分为四个章节,下面简要介绍一下各章的主要内容。 第一章,我们做一些准备工作。首先,给出了在文献[88,89,90,91,92]中提出并研究的海水入侵的数学模型,该模型是由一个抛物方程和一个对流扩散方程组成的耦合系统。对于此系统,我们给出了两种不同的边值条件,并在后面几章分别给出相应的处理格式。接下来,给出了一些定义和记号,包括Sobolev空间和分片Sobolev空间及它们的范数、空间剖分及其属性、标量和矢量跳跃和平均的定义等。最后,以引理的形式给出了在误差估计中必需的定理和不等式,这些引理在许多的书籍和文献中都曾被提到或证明。 第二章,对于前一章给出的数学模型和第一个边值条件,我们提出一个连续时间格式。混合元方法被用来处理数学模型中的抛物方程,而间断伽略金有限元方法(对称内罚方法)被用来处理系统的对流扩散方程,其中对流项只是单纯做了内积处理。在给出了上述方程的变分形式的基础上,我们提出了迭代的连续时间格式。对于该格式,得到了L2范数和H1范数下的最优误差估计。这一部分内容主要出自文章[58],该文章已被Journal of Systems Science and Complexity(SCI)接收。 接着,针对上述连续时间格式,提出了一些时间离散方法,涉及显式和隐式方法以及一些高阶方法。对于这些方法,我们简要地研究了它们的稳定性和收敛性,并讨论了其优缺点。 第三章,我们主要考虑附加了第二个边值条件的海水入侵模型。模型中的抛物方程仍然采用混合有限元方法来求解。而对于对流扩散方程,为了更好的逼近方程中的对流项,我们使用迎风间断伽略金有限元方法来处理。另外,通过引入一个参数λ和θ,我们在同一个表达式中考虑了两种不同的内罚方法：对称内罚和不对称内罚方法以及两种时间离散方法：向后Euler方法和Crank-Nicolson方法。在此基础上,我们得到了全隐式的迭代格式。理论上的hp先验误差估计和设计的数值实验结果都证明此格式的有效性和可靠性。这些结果主要出自文章[59]。 众所周知,非线性问题的全隐格式的求解效率很低,而显格式是条件稳定的,时间步长会受到限制。为了改变这种两难处境,我们提出了一种半隐格式,即对于线性的扩散项采用隐式处理,而非线性的对流项采用显式处理。理论分析和数值实验都表明此格式的计算精度与全隐格式相差无几。 第四章,我们讨论了一种求解抛物方程的基于间断有限元方法的区域分解算法。由于上一章中的半隐格式可以看成是对抛物方程的全隐格式,所以这一章的并行方法可以有效的提高半隐格式的计算速度。我们以一个标准的抛物方程为例,提出了区域分解算法,并且考虑它的稳定性和收敛性,最后给出了数值算例以验证我们的结论。这一章的结果出自文章[76]。