裂隙介質(zhì)溶質(zhì)運(yùn)移模型的數(shù)值模擬*

宋證遠(yuǎn), 李煥榮

(1.重慶工商大學(xué) 體育學(xué)院，重慶 400067；2.重慶工商大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，重慶 400067)

裂隙介質(zhì)溶質(zhì)運(yùn)移模型的數(shù)值模擬*

宋證遠(yuǎn)1, 李煥榮2

(1.重慶工商大學(xué) 體育學(xué)院，重慶 400067；2.重慶工商大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，重慶 400067)

考慮裂隙介質(zhì)中溶質(zhì)運(yùn)移問(wèn)題的數(shù)學(xué)模型，建立離散的廣義差分格式，給出不同的裂隙介質(zhì)阻滯系數(shù)下的數(shù)值算例，并將數(shù)值解與解析解進(jìn)行比較；比較結(jié)果表明：利用所建廣義差分格式求解裂隙介質(zhì)中的溶質(zhì)運(yùn)移問(wèn)題是可靠的，且該格式具有穩(wěn)定性和可實(shí)用性，可以用來(lái)數(shù)值模擬更加復(fù)雜的裂隙介質(zhì)中溶質(zhì)運(yùn)移的物理過(guò)程.

裂隙介質(zhì)；溶質(zhì)運(yùn)移問(wèn)題；數(shù)值模擬；廣義差分法

裂隙介質(zhì)中溶質(zhì)運(yùn)移的研究是從核工業(yè)發(fā)展起來(lái)后，為了核廢料的地質(zhì)儲(chǔ)存而發(fā)展起來(lái)的.隨著人類活動(dòng)不斷向地下深入，各種污染對(duì)地下深部裂隙水的威脅也逐漸加重，因此許多國(guó)家的學(xué)者開展了裂隙中溶質(zhì)運(yùn)移規(guī)律的探討和研究[1-5]．描述裂隙中溶質(zhì)運(yùn)移的控制方程為對(duì)流-擴(kuò)散方程，可以用此方程進(jìn)行數(shù)值模擬和預(yù)測(cè)，但最終要?dú)w結(jié)為對(duì)該方程進(jìn)行數(shù)值求解．在實(shí)際應(yīng)用中，數(shù)值計(jì)算方法主要以有限差分和有限元法為主[6-8]．但是用通常的有限元方法或有限差分法求出的解可能發(fā)生振蕩，即數(shù)值彌散，不符合物理要求[9]．采用廣義差分法[10-12]，既提高了計(jì)算精度，減少了計(jì)算量，又能有效地克服數(shù)值彌散．本文采用此方法來(lái)求解裂隙介質(zhì)中溶質(zhì)運(yùn)移方程，給出了不同裂隙介質(zhì)阻滯系數(shù)下的數(shù)值計(jì)算結(jié)果，并與文獻(xiàn)[5]中的試驗(yàn)測(cè)量值及其由解析-優(yōu)化法確定的解析解的計(jì)算結(jié)果進(jìn)行了比較，以探討其實(shí)用性．

1 裂隙介質(zhì)溶質(zhì)運(yùn)移模型及廣義差分格式

考慮單裂隙中溶質(zhì)運(yùn)移一維模型的定解問(wèn)題[5]:求c(z，t)，使得對(duì)于任意的T>0，滿足

(1)

式(1)中，c(z，t)表示溶質(zhì)濃度，Rd為阻滯系數(shù)，DL為縱向彌散系數(shù)，u為水流速度，c1(z)和c0分別表示溶質(zhì)初始濃度和上邊界的溶質(zhì)濃度，z軸垂直向下.

?u,v∈L2(I).

顯然，dimUh=dimVh=l.

(3)

(4)

(5)

(6)

(7)

2 數(shù)值算例

(8)

i=1,2,…,l

其中，

(9)

當(dāng)i=2,3,…,l-1時(shí)，式(8)可寫為

(10)

當(dāng)i=l時(shí)，式(8)可寫為

(11)

整理式(9)—(11)可得：

i=1時(shí)，

(12)

i=2,3,…,l-1時(shí)，

(13)

i=l時(shí)，

(14)

將式(12)—(14)寫成矩陣形式：

(15)

其中，

(16)

因此，只要給定時(shí)間步長(zhǎng)τ，空間剖分步長(zhǎng)h以及初邊值和參數(shù)Rd,DL,u的值，就可以由式(15)(16)求出問(wèn)題(1)的廣義差分解，即溶質(zhì)濃度的數(shù)值解．

在文獻(xiàn)[5]的裂隙介質(zhì)溶質(zhì)運(yùn)移試驗(yàn)中，初始濃度c1(z，0)=0，混凝土試件長(zhǎng)L=80 cm，裂隙開度為0.8 mm，通過(guò)調(diào)節(jié)進(jìn)水端和出水端的壓力，應(yīng)用達(dá)西定律求出裂隙中的水流速度u，得到了常壓下首次試驗(yàn)A的測(cè)量值(如圖1(a))和多次試驗(yàn)后的試驗(yàn)B的測(cè)量值(如圖1(b))，并根據(jù)水動(dòng)力彌散理論，用解析-優(yōu)化方法求出溶質(zhì)運(yùn)移的參數(shù)DL= 2.902 cm2/min，Rd=4.488(實(shí)驗(yàn)A)和Rd=1.0(實(shí)驗(yàn)B)，從而定解問(wèn)題(1)的解析解可以表示為

(17)

圖1 常壓下不同流量時(shí)通過(guò)距離濃度源23 cm處溶質(zhì)濃度與注入溶質(zhì)濃度的比值隨時(shí)間的變化曲線Fig.1 The profile of solute concentration to the injected solute concentration at 23 cm under atmospheric pressure

其中，文獻(xiàn)[16]中已導(dǎo)出誤差函數(shù)erfc(x)的近似表達(dá)式：

(18)

從圖1對(duì)比曲線圖可以看出，用擬合參數(shù)計(jì)算結(jié)果和實(shí)測(cè)結(jié)果很接近，表明由解析-優(yōu)化法確定參數(shù)得到的解析解(17)是可靠的．因此本文就只給出用廣義差分法求出的數(shù)值解與由式(17)求得的解析解的比較，如圖2和圖3，其中(a)表示實(shí)驗(yàn)A，(b)表示實(shí)驗(yàn)B.

圖2 通過(guò)距離濃度源不同距離處的溶質(zhì)濃度隨時(shí)間的變化曲線Fig.2 The solution concentration profile with time at different distances

從圖2和圖3可以看出兩個(gè)試驗(yàn)A和B受裂隙介質(zhì)的影響，試驗(yàn)A中溶質(zhì)濃度變化較試驗(yàn)B要慢，變化曲線也較平緩，這是由于試驗(yàn)A是在單裂隙介質(zhì)中的首次試驗(yàn)，裂隙壁對(duì)溶液中離子的吸附作用強(qiáng)烈，用相同的溶液進(jìn)行多次試驗(yàn)后，這種吸附作用會(huì)逐漸減弱，直至不產(chǎn)生吸附作用，故試驗(yàn)B的溶質(zhì)濃度變化要快且劇烈．并且從圖3還可以看出，當(dāng)溶液入滲開始后，接近地表處的溶質(zhì)濃度迅速上升，之后濃度變化不大，并逐步接近注入的溶質(zhì)濃度．

同時(shí)從圖2和圖3所表示的廣義差分解(點(diǎn)線)和解析解(實(shí)線)的曲線圖可以看出，本文格式的計(jì)算結(jié)果與解析解吻合得很好，即利用廣義差分格式求解裂隙介質(zhì)中的溶質(zhì)運(yùn)移問(wèn)題是可靠的，本文格式具有穩(wěn)定性和可實(shí)用性．

3 結(jié) 論

用廣義差分法求解裂隙介質(zhì)中溶質(zhì)運(yùn)移方程，構(gòu)造了此類對(duì)流-擴(kuò)散方程的廣義差分格式，進(jìn)行了數(shù)值模擬，并與文獻(xiàn)[5]中求得的解析表達(dá)式的計(jì)算結(jié)果進(jìn)行了比較．從數(shù)值計(jì)算過(guò)程和比較結(jié)果來(lái)看，應(yīng)用廣義差分格式后，不僅具有較高的精度，計(jì)算量少，計(jì)算程序?qū)崿F(xiàn)簡(jiǎn)單，而且能有效地克制數(shù)值彌散現(xiàn)象，因此可以用來(lái)數(shù)值模擬更加復(fù)雜的裂隙介質(zhì)中溶質(zhì)運(yùn)移的物理過(guò)程．

[1] HOEKSTRA P.Moisture Movement in Soil under Temperature Gradients with Cold Side Temperature Below Freezing[J].Water Resources Researeh,1966(2):241-250

[2] LUNARDINI V J.Heat Transfer in Cold Climates[M].Van Nostrand-Reinhold,New York,1981

[3] 安維東.渠道凍結(jié)時(shí)熱質(zhì)遷移的數(shù)值模擬[J].冰川凍土,1985(1):35-46

AN W D.Numerical Simulation Analysis Masstransfer under Freezing[J].Journal of Glaciology and Geocryology,1985(1):35-46

[4] 徐學(xué)祖,王家澄,張立新.凍土物理學(xué)[M].北京:科學(xué)出版社,2001

XU X Z,WANG J C,ZHANG L X.Frozen Soil Physics[M].Beijing:Science Press,2001

[5] 孫增奎,王連俊,白明洲,等.青藏鐵路多年凍土路堤溫度場(chǎng)的有限元分析[J].巖石力學(xué)與工程學(xué)報(bào)，2004,23(20):3454-3459

SUN Z K,WANG L J,BAI M Z,et al.Finite Element Analysis on Temperature Field of Qing-Tibet Railway Embankment on Permafrost[J].Chinese Journal of Rock Mechanics and Engineering,2004,23(20):3454-3459

[6] 米隆,賴遠(yuǎn)明,吳紫汪,等.凍土鐵路路基溫度特性的有限元分析[J].鐵道學(xué)報(bào),2003,25(2):62-67

MI L,LAI Y M,WU Z J,et al.The Finite Element Analysis for Temperature Property of Railway Embankments in Cold Regions[J].Journal of the China Railway Society,2003,25(2):62-67

[7] WANG T,ZHOU G Q.Neumann Stochastic Finite Element Method for Calculating Temperature Field of Frozen Soil Based on Random Field Theory[J].Sciences in Cold and Arid Regions,2013(4):488-497

[8] WANG C,QIU Z P.Interval Finite Difference Method for Steady-State Temperature Field Prediction with Interval Parameters[J].Acta Mechanica Sinica,2014,30(2):161-166

[9] WANG C,QIU Z P.Fuzzy Finite Difference Method for Heat Conduction Analysis with Uncertain Parameters[J].Acta Mechanica Sinica,2014,30(3):383-390

[10] 李煥榮.土壤水流與溶質(zhì)耦合運(yùn)移問(wèn)題的混合元-迎風(fēng)廣義差分法研究[J].計(jì)算數(shù)學(xué),2013,35(1):1-10

Li H R.Mixed Finite Element and Upwind Generalized Difference Method for Coupled Transport Models of Unsaturated Soil Water Flow and Solute[J].Mathematica Numerica Sinica,2013,35(1):1-10

[11] 李煥榮.水流入滲問(wèn)題的數(shù)值模擬[J].重慶工商大學(xué)學(xué)報(bào)(自然科學(xué)版),2009,26(3):213-218

Li H R.Numerical Simulation for the Infiltration Problems of Water Flow[J].Journal of Chongqing Technology and Business University (Natural Science Edition),2009,26(3):213-218

[12] 宋琿.季節(jié)性凍土地區(qū)基穩(wěn)定性的數(shù)值分析[D].成都:西南交通大學(xué)碩士學(xué)位論文，2006

SONG H.Numerical Analysis for the Stability of Seasonal Frozen Subgrade[D].Chengdu:Southwest Jiaotong Univer-sity,2006

責(zé)任編輯：李翠薇

Numerical Simulation for Solute Transport Model in Fractured Media

SONG Zheng-yuan1, LI Huan-rong2

(1. School of Physical Education, Chongqing Technology and Business University, Chongqing 400067, China; 2. School of Mathematics and Statistics, Chongqing Technology and Business University, Chongqing 400067, China)

This paper considers the mathematical model for solute transport in fractured media, sets up the discrete generalized difference scheme, gives numerical examples under different fractured media retardation coefficients and compares the numerical solution with analytic solution.Comparison results show that the solution of established generalized difference scheme to solute transport in fractured media is reliable, that this scheme is stable and practical and can be used to numerically simulate the physical process for the solute transport in more complex fractured media.

fractured media; solute transport; numerical simulation; generalized difference method

10.16055/j.issn.1672-058X.2017.0003.001

2016-11-17；

2016-12-21. * 基金項(xiàng)目：國(guó)家自然科學(xué)基金(11101453，11201508);重慶市自然科學(xué)基金(2013JCYJA20015)；重慶市教委科技項(xiàng)目(KJ1400602).

宋證遠(yuǎn)(1975-)，男，山東臨沂人，碩士，從事體育運(yùn)動(dòng)模型計(jì)算研究．

O241.82

A

1672-058X(2017)03-0001-06