鄧永輝，鄧永紅

(湖南財(cái)政經(jīng)濟(jì)學(xué)院，湖南 長(zhǎng)沙 410205)

溶質(zhì)運(yùn)移模型的有限元數(shù)值解

鄧永輝，鄧永紅

(湖南財(cái)政經(jīng)濟(jì)學(xué)院，湖南 長(zhǎng)沙 410205)

將混合拉普拉斯變換有限單元法引入到求解首采區(qū)鹵水動(dòng)態(tài)二維模型的溶質(zhì)運(yùn)移問(wèn)題中，能夠有限地消除在求解對(duì)流占優(yōu)的地下水溶質(zhì)運(yùn)移問(wèn)題時(shí)產(chǎn)生的數(shù)值擴(kuò)散和過(guò)量的現(xiàn)象，具有一步到位、局部求解的優(yōu)點(diǎn)，最后將該方法應(yīng)用到具有空間一階導(dǎo)數(shù)項(xiàng)的對(duì)流彌散方程，以檢驗(yàn)該方法的數(shù)值有效性和求解溶質(zhì)運(yùn)移模型的能力.

對(duì)流擴(kuò)散方程;拉普拉斯變換;混合法

察爾汗鹽湖首采區(qū)含水介質(zhì)為石鹽層，高濃度鹵水溶液(礦化度300～450 g·L－1)在這種易溶巖介質(zhì)中的溶質(zhì)運(yùn)移涉及到一個(gè)突出的問(wèn)題是鹵水與介質(zhì)間的物相轉(zhuǎn)化即固液轉(zhuǎn)化問(wèn)題［1－2］.而對(duì)于在開(kāi)采過(guò)程中鹵水與易溶介質(zhì)之間的固液轉(zhuǎn)化問(wèn)題，由于涉及到礦床地質(zhì)、水文地質(zhì)、物理化學(xué)、地球化學(xué)、多孔介質(zhì)水動(dòng)力彌散理論及水鹽體系相圖原理等多方面的知識(shí)，目前國(guó)內(nèi)外尚無(wú)可借鑒的經(jīng)驗(yàn)和理論依據(jù).

1 模型的建立

1.1 對(duì)流擴(kuò)散方程在建立目標(biāo)模型之前先了解一下一維單向的流場(chǎng)的對(duì)流擴(kuò)散方程［1］.

對(duì)于一維的對(duì)流擴(kuò)散模型，可以用質(zhì)量守衡方程求出微分方程解等方法來(lái)確定.周期一個(gè)長(zhǎng)為Δx圓柱形內(nèi)物質(zhì)遷移表示為

其中，Dl為擴(kuò)散系數(shù)(l表示擴(kuò)散方向)，U為真實(shí)滲透速度，C為液體濃度.

1.2 二維對(duì)流擴(kuò)散方程從二維流的滲流場(chǎng)中割離出一個(gè)微分單元體，該單元的寬度為Δx，長(zhǎng)度為Δy，厚度為M.

水沿x軸方向從左面流入單元體;沿y軸方向從前流入.經(jīng)時(shí)段Δt后分別經(jīng)Δx，Δy從對(duì)面流出.

考察滲流引起的濃度變化:

a)順x軸方向由滲流引起物質(zhì)遷移即滲流擴(kuò)散模型，可用質(zhì)量守衡方程求出微分方程解等方法來(lái)確定.在長(zhǎng)為Δx的單元體內(nèi)物質(zhì)遷移可表示為

b)順y軸方向由滲流引起物質(zhì)遷移即滲流擴(kuò)散模型，可用質(zhì)量守衡方程求出微分方程解等方法來(lái)確定.在長(zhǎng)為Δy的單元體內(nèi)物質(zhì)遷移可表示為

c)輸入輸出可以包括彌(擴(kuò))散和對(duì)流引起的現(xiàn)象.式(3)和(4)被截面積和Δt去除后，使后二項(xiàng)系數(shù)為零，即除以ΔxΔyΔt分別可得到以下方程

1.3 匯源補(bǔ)給項(xiàng)和固液轉(zhuǎn)化項(xiàng)計(jì)算區(qū)含水層垂直向交換量包括大氣降水補(bǔ)給量、晶間鹵水蒸發(fā)量、渠系采鹵(回滲)量和下伏含水層的越流補(bǔ)給量.

匯源補(bǔ)給量引起的濃度變化=CQ－C.

溶質(zhì)運(yùn)移方程中固液轉(zhuǎn)化是人們長(zhǎng)期探索的問(wèn)題之一，筆者對(duì)此不作太多的研究，假設(shè)固液轉(zhuǎn)化系數(shù)f是常數(shù)，也就是說(shuō)假設(shè)由于固液轉(zhuǎn)化帶來(lái)的濃度變化MμfC為常量.一般來(lái)講，研究固液轉(zhuǎn)化的方法主要是非平衡化學(xué)法，假定地下水系統(tǒng)中有幾種不同的物理﹑化學(xué)和生物化學(xué)作用過(guò)程，用平衡化學(xué)法判斷這些作用是否平衡，用反應(yīng)動(dòng)力學(xué)描述固液轉(zhuǎn)化速率.但在目前，非化學(xué)平衡法還處于探索階段，尤其對(duì)高濃度鹵水的計(jì)算，還沒(méi)有一種較為成熟的方法.

1.4 最終數(shù)學(xué)模型通過(guò)上述幾種簡(jiǎn)單的模型的推導(dǎo)和分析，再結(jié)合首采區(qū)能引起濃度變化的各種因素，如對(duì)流、擴(kuò)散、排泄、補(bǔ)給等等，可以導(dǎo)出最終目標(biāo)模型

其中，V1，V2為滲透速度(L·s－1);D11，D22為彌散系數(shù)(L2·s－1).式(10) 只表明濃度隨時(shí)間的變化的規(guī)律，要求出微分方程的解還需要一些定解條件，求出在特定條件下濃度的值，在計(jì)算區(qū)內(nèi)邊界條件如下

2 用HLTFEM解目標(biāo)方程

將HLTFEM求解思路引入到察爾汗鹽湖采鹵方案中溶質(zhì)運(yùn)移的計(jì)算，盡管HLTFEM嚴(yán)格受限于穩(wěn)定流場(chǎng)線(xiàn)性溶質(zhì)傳輸，但是在察爾汗鹽湖首采區(qū)流速和傳輸參數(shù)可以是空間變量的函數(shù)，求解區(qū)域可以是不規(guī)則的，允許邊界條件是時(shí)間變量的一般函數(shù)，這就使得這種新的沒(méi)有時(shí)間步長(zhǎng)、定點(diǎn)求解的計(jì)算方法仍適于溶質(zhì)運(yùn)移問(wèn)題的模擬.

在求解溶質(zhì)運(yùn)移方程時(shí)，由于該問(wèn)題的復(fù)雜性，因此，文中忽略固液轉(zhuǎn)化作用，暫不考慮匯源項(xiàng)，則溶質(zhì)運(yùn)移方程可寫(xiě)為［3］

設(shè)基本函數(shù)

其中，(xj，yj)是第j個(gè)結(jié)點(diǎn)的坐標(biāo).

將區(qū)域Ω剖分成三角形網(wǎng)，三角形的頂點(diǎn)取為結(jié)點(diǎn).設(shè)任一三角形單元(△)的3個(gè)頂點(diǎn)的結(jié)點(diǎn)號(hào)碼為i，j，k，坐標(biāo)分別為(xi，yi)，(xj，yj)和(xk，yk)，規(guī)定和這3 個(gè)結(jié)點(diǎn)相聯(lián)系的基函數(shù)在單元(△)中的值為［6］

首先，形成［A］和F在該單元中的部分，其次所有單元疊加形成整體的［A］和F，并結(jié)合邊界條件，式(17)就建立起所需要的方程組［A］C+F=0，由式(24)知，系數(shù)矩陣是高度稀疏非對(duì)稱(chēng)復(fù)值矩陣.為了節(jié)省計(jì)算機(jī)內(nèi)存，使用壓縮存貯的技巧，將方程非零系數(shù)按最大帶寬存入二維數(shù)組中，然后根據(jù)各計(jì)算結(jié)點(diǎn)的相鄰結(jié)點(diǎn)編號(hào)和相鄰結(jié)點(diǎn)的個(gè)數(shù)，采用高斯消去法求解此二維數(shù)組，即可獲得象空間的濃度分布.

根據(jù)離散的有限元方程組的解，對(duì)拉氏變換后結(jié)點(diǎn)的濃度Cj進(jìn)行反演.若以L(fǎng)－1表示拉氏反演，則式(18)可化成

其中，Cj(t)是在結(jié)點(diǎn)j處時(shí)間域的濃度.

本文采用Honig和Hirdes提出的基于Fourier級(jí)數(shù)展開(kāi)的拉氏反演新算法進(jìn)行數(shù)值反演［7］，其反演公式為

3 結(jié)論

基于Fourier級(jí)數(shù)展開(kāi)的拉氏變換反演算法，如Grump算法，最大的缺點(diǎn)就是離散誤差和截?cái)嗾`差依賴(lài)于自由參數(shù)，即通過(guò)選擇合適的稀有參數(shù)使離散誤差變得任意小，但同時(shí)截?cái)嗾`差又變得無(wú)窮大，反之亦然.Honig-Hirdes新算法［7］通過(guò)同步使用減少離散誤差的方法和加速Fourier級(jí)數(shù)收斂的方法以及近似計(jì)算最優(yōu)自由參數(shù)的方法克服了Grump等算法之不足.但Stehfest算法由于僅僅使用拉氏變換參數(shù)S的實(shí)部［4］，所以與使用復(fù)數(shù)的Honig-Hirdes算法相比，在拉氏變換與Galerkin有限元結(jié)合時(shí)，會(huì)喪失很多優(yōu)點(diǎn)，因?yàn)榛谝粋€(gè)實(shí)數(shù)S的變換后的濃度剖面，不再是一個(gè)光滑的震蕩函數(shù)而與時(shí)間域中的濃度剖面的特性相似.

［1］孫納正.地下水污染——數(shù)學(xué)模型和數(shù)值方法［M］.北京:北京地質(zhì)出版社，1989.

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

［3］王文科.地下水有限分析數(shù)值模擬的理論與方法［M］.西安:陜西科學(xué)技術(shù)出版社，1996:102－126.

［4］蔣曉蓉.油藏?cái)?shù)值模擬基礎(chǔ)［M］.成都:成都理工大學(xué)出版社，1998.

［5］羅煥炎，陳雨孫.地下水運(yùn)動(dòng)的數(shù)值模擬［M］.北京:中國(guó)建筑工業(yè)出版社，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">文獻(xiàn)標(biāo)志碼:A

A

1004－1729(2011)01－0025－04

2010－10－11

國(guó)家攻關(guān)項(xiàng)目(2001BA60ZB－09－1)

鄧永輝(1979－)，女，湖南常德人，湖南財(cái)政經(jīng)濟(jì)學(xué)院講師，碩士.