溶质运移模型的有限元数值解

邓永辉，邓永红

(湖南财政经济学院，湖南 长沙 410205)

溶质运移模型的有限元数值解

邓永辉，邓永红

(湖南财政经济学院，湖南 长沙 410205)

将混合拉普拉斯变换有限单元法引入到求解首采区卤水动态二维模型的溶质运移问题中，能够有限地消除在求解对流占优的地下水溶质运移问题时产生的数值扩散和过量的现象，具有一步到位、局部求解的优点，最后将该方法应用到具有空间一阶导数项的对流弥散方程，以检验该方法的数值有效性和求解溶质运移模型的能力.

对流扩散方程;拉普拉斯变换;混合法

察尔汗盐湖首采区含水介质为石盐层，高浓度卤水溶液(矿化度300～450 g·L－1)在这种易溶岩介质中的溶质运移涉及到一个突出的问题是卤水与介质间的物相转化即固液转化问题［1－2］.而对于在开采过程中卤水与易溶介质之间的固液转化问题，由于涉及到矿床地质、水文地质、物理化学、地球化学、多孔介质水动力弥散理论及水盐体系相图原理等多方面的知识，目前国内外尚无可借鉴的经验和理论依据.

1 模型的建立

1.1 对流扩散方程在建立目标模型之前先了解一下一维单向的流场的对流扩散方程［1］.

对于一维的对流扩散模型，可以用质量守衡方程求出微分方程解等方法来确定.周期一个长为Δx圆柱形内物质迁移表示为

其中，Dl为扩散系数(l表示扩散方向)，U为真实渗透速度，C为液体浓度.

1.2 二维对流扩散方程从二维流的渗流场中割离出一个微分单元体，该单元的宽度为Δx，长度为Δy，厚度为M.

水沿x轴方向从左面流入单元体;沿y轴方向从前流入.经时段Δt后分别经Δx，Δy从对面流出.

考察渗流引起的浓度变化:

a)顺x轴方向由渗流引起物质迁移即渗流扩散模型，可用质量守衡方程求出微分方程解等方法来确定.在长为Δx的单元体内物质迁移可表示为

b)顺y轴方向由渗流引起物质迁移即渗流扩散模型，可用质量守衡方程求出微分方程解等方法来确定.在长为Δy的单元体内物质迁移可表示为

c)输入输出可以包括弥(扩)散和对流引起的现象.式(3)和(4)被截面积和Δt去除后，使后二项系数为零，即除以ΔxΔyΔt分别可得到以下方程

1.3 汇源补给项和固液转化项计算区含水层垂直向交换量包括大气降水补给量、晶间卤水蒸发量、渠系采卤(回渗)量和下伏含水层的越流补给量.

汇源补给量引起的浓度变化=CQ－C.

溶质运移方程中固液转化是人们长期探索的问题之一，笔者对此不作太多的研究，假设固液转化系数f是常数，也就是说假设由于固液转化带来的浓度变化MμfC为常量.一般来讲，研究固液转化的方法主要是非平衡化学法，假定地下水系统中有几种不同的物理﹑化学和生物化学作用过程，用平衡化学法判断这些作用是否平衡，用反应动力学描述固液转化速率.但在目前，非化学平衡法还处于探索阶段，尤其对高浓度卤水的计算，还没有一种较为成熟的方法.

1.4 最终数学模型通过上述几种简单的模型的推导和分析，再结合首采区能引起浓度变化的各种因素，如对流、扩散、排泄、补给等等，可以导出最终目标模型

其中，V1，V2为渗透速度(L·s－1);D11，D22为弥散系数(L2·s－1).式(10) 只表明浓度随时间的变化的规律，要求出微分方程的解还需要一些定解条件，求出在特定条件下浓度的值，在计算区内边界条件如下

2 用HLTFEM解目标方程

将HLTFEM求解思路引入到察尔汗盐湖采卤方案中溶质运移的计算，尽管HLTFEM严格受限于稳定流场线性溶质传输，但是在察尔汗盐湖首采区流速和传输参数可以是空间变量的函数，求解区域可以是不规则的，允许边界条件是时间变量的一般函数，这就使得这种新的没有时间步长、定点求解的计算方法仍适于溶质运移问题的模拟.

在求解溶质运移方程时，由于该问题的复杂性，因此，文中忽略固液转化作用，暂不考虑汇源项，则溶质运移方程可写为［3］

设基本函数

其中，(xj，yj)是第j个结点的坐标.

将区域Ω剖分成三角形网，三角形的顶点取为结点.设任一三角形单元(△)的3个顶点的结点号码为i，j，k，坐标分别为(xi，yi)，(xj，yj)和(xk，yk)，规定和这3 个结点相联系的基函数在单元(△)中的值为［6］

首先，形成［A］和F在该单元中的部分，其次所有单元叠加形成整体的［A］和F，并结合边界条件，式(17)就建立起所需要的方程组［A］C+F=0，由式(24)知，系数矩阵是高度稀疏非对称复值矩阵.为了节省计算机内存，使用压缩存贮的技巧，将方程非零系数按最大带宽存入二维数组中，然后根据各计算结点的相邻结点编号和相邻结点的个数，采用高斯消去法求解此二维数组，即可获得象空间的浓度分布.

根据离散的有限元方程组的解，对拉氏变换后结点的浓度Cj进行反演.若以L－1表示拉氏反演，则式(18)可化成

其中，Cj(t)是在结点j处时间域的浓度.

本文采用Honig和Hirdes提出的基于Fourier级数展开的拉氏反演新算法进行数值反演［7］，其反演公式为

3 结论

基于Fourier级数展开的拉氏变换反演算法，如Grump算法，最大的缺点就是离散误差和截断误差依赖于自由参数，即通过选择合适的稀有参数使离散误差变得任意小，但同时截断误差又变得无穷大，反之亦然.Honig-Hirdes新算法［7］通过同步使用减少离散误差的方法和加速Fourier级数收敛的方法以及近似计算最优自由参数的方法克服了Grump等算法之不足.但Stehfest算法由于仅仅使用拉氏变换参数S的实部［4］，所以与使用复数的Honig-Hirdes算法相比，在拉氏变换与Galerkin有限元结合时，会丧失很多优点，因为基于一个实数S的变换后的浓度剖面，不再是一个光滑的震荡函数而与时间域中的浓度剖面的特性相似.

［1］孙纳正.地下水污染——数学模型和数值方法［M］.北京:北京地质出版社，1989.

［2］ROACHE P J.Computational fluid dynamics［M］.Hermosa:Albuquereque，1976:446 －447.

［3］王文科.地下水有限分析数值模拟的理论与方法［M］.西安:陕西科学技术出版社，1996:102－126.

［4］蒋晓蓉.油藏数值模拟基础［M］.成都:成都理工大学出版社，1998.

［5］罗焕炎，陈雨孙.地下水运动的数值模拟［M］.北京:中国建筑工业出版社，2001.

［6］DURBIN F.Numerical inversion of Laplace transform:an efficient improvement to Durbner and Abare’s method［J］.Comp.J.，1993，17:371 －376.

Finite Element Numerical Solution for a Solute Transport Model

DENG Yong-hui，DENG Yong-hong

(Hunan University of Finance and Economics，Changsha 410205，China)

In the paper，the hybrid laplace transform finite element method was used for solving solute transport problems of dynamic two-dimensional model of brine water in the first exploitation area，which can limitedly eliminate numerical diffusion and overdo phenomena from the solving convection dominated underwater solute transport problems，whose advantages were one step and local solving.And the method was used for convectivediffusion equation which have first derivative to test the numerical effectiveness and the capability to solve solute transport model.

convective-diffusion equation;laplace transform;hybrid method

O 35 < class="emphasis_bold">文献标志码:A

A

1004－1729(2011)01－0025－04

2010－10－11

国家攻关项目(2001BA60ZB－09－1)

邓永辉(1979－)，女，湖南常德人，湖南财政经济学院讲师，硕士.