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

鄧永輝，鄧永紅

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

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

鄧永輝，鄧永紅

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

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

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

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

1 模型的建立

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié)論

基于Fourier級數(shù)展開的拉氏變換反演算法，如Grump算法，最大的缺點就是離散誤差和截斷誤差依賴于自由參數(shù)，即通過選擇合適的稀有參數(shù)使離散誤差變得任意小，但同時截斷誤差又變得無窮大，反之亦然.Honig-Hirdes新算法［7］通過同步使用減少離散誤差的方法和加速Fourier級數(shù)收斂的方法以及近似計算最優(yōu)自由參數(shù)的方法克服了Grump等算法之不足.但Stehfest算法由于僅僅使用拉氏變換參數(shù)S的實部［4］，所以與使用復(fù)數(shù)的Honig-Hirdes算法相比，在拉氏變換與Galerkin有限元結(jié)合時，會喪失很多優(yōu)點，因為基于一個實數(shù)S的變換后的濃度剖面，不再是一個光滑的震蕩函數(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］蔣曉蓉.油藏數(shù)值模擬基礎(chǔ)［M］.成都:成都理工大學(xué)出版社，1998.

［5］羅煥炎，陳雨孫.地下水運(yùn)動的數(shù)值模擬［M］.北京:中國建筑工業(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ān)項目(2001BA60ZB－09－1)

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