孫露露，楊衛(wèi)波

(揚(yáng)州大學(xué)水利與能源動(dòng)力工程學(xué)院，江蘇揚(yáng)州 225127)

溫度梯度作用下非飽和土壤中溶質(zhì)遷移過程的數(shù)值模擬

孫露露，楊衛(wèi)波*

(揚(yáng)州大學(xué)水利與能源動(dòng)力工程學(xué)院，江蘇揚(yáng)州 225127)

根據(jù)非飽和土壤非等溫Richard方程、傳熱與溶質(zhì)傳遞方程，建立溶質(zhì)遷移模型，研究溫度梯度作用下入滲速度對(duì)土壤溫度變化及溶質(zhì)遷移規(guī)律的影響.數(shù)值模擬的結(jié)果表明:溶質(zhì)遷移過程中，土壤溫度與溶質(zhì)的物質(zhì)的量濃度均隨時(shí)間延續(xù)而逐漸升高，隨深度的加深而逐漸降低;入滲速度的增加可促進(jìn)土壤中熱量的擴(kuò)散和對(duì)流，增強(qiáng)溶質(zhì)的對(duì)流和彌散作用.

非飽和土壤;溫度梯度;滲流;溶質(zhì)遷移;數(shù)值模擬

近年來，溶質(zhì)在土壤中的遷移機(jī)制已成為當(dāng)前土壤環(huán)境控制領(lǐng)域的一個(gè)研究熱點(diǎn).非飽和土壤是由固、液、氣三相組成的多孔介質(zhì)，溶質(zhì)在沿非飽和土壤的水流運(yùn)動(dòng)過程中會(huì)導(dǎo)致土壤污染.Viotti等［1］建立了一維非飽和土壤中溶質(zhì)遷移方程，研究了不同輸入?yún)?shù)對(duì)溶質(zhì)濃度分布的影響;Kumar等［2］通過無網(wǎng)格Galerkin(element free galerkin，EFG)法分析了二維非飽和多孔介質(zhì)中溶質(zhì)遷移問題，并與有限元法進(jìn)行了算例比較;Javadi等［3］對(duì)非飽和土壤中水、氣和溶質(zhì)的遷移進(jìn)行了案例分析;Ghasemzadeh［4］研究了變形多孔介質(zhì)中熱量、水流和溶質(zhì)的遷移過程;Winiarski等［5］實(shí)驗(yàn)研究了非飽和土壤的非均質(zhì)性對(duì)多相流和溶質(zhì)遷移特性的影響;眭素剛等［6］開展了非飽和土壤彌散、滲流特性試驗(yàn)，并模擬分析了分層土壤中溶質(zhì)的濃度分布;左自波等［7］研究了不同降雨條件下非飽和帶孔隙水壓力、含水量變化及溶質(zhì)的遷移規(guī)律.鑒于溶質(zhì)類型的多樣性，本文以某一特定濃度的無機(jī)液體A為代表，在忽略吸附、衰減及化學(xué)性質(zhì)的條件下，根據(jù)多孔介質(zhì)傳熱、傳質(zhì)及滲流理論，建立非飽和土壤中溶質(zhì)遷移的數(shù)學(xué)模型，研究入滲速度對(duì)非飽和土壤溫度及溶質(zhì)的物質(zhì)的量濃度c(A)分布規(guī)律的影響，以及滲流作用下非飽和土壤的傳熱傳質(zhì)規(guī)律.

1 計(jì)算模型

1.1 物理模型

為模擬滲流影響下土壤溫度和溶質(zhì)濃度的時(shí)空變化規(guī)律，以非飽和土壤為研究對(duì)象，建立溶質(zhì)遷移模型，進(jìn)行滲流場(chǎng)、溫度場(chǎng)和溶質(zhì)傳遞場(chǎng)三場(chǎng)耦合的二維數(shù)值模擬，研究溫度及溶質(zhì)的物質(zhì)的量濃度的分布規(guī)律，模擬對(duì)象結(jié)構(gòu)如圖1所示.圖1中ABCD區(qū)域?yàn)榉秋柡屯寥?，模型尺寸為AB=1 m，AD=0.5 m;E(0.5，0.45)為溶質(zhì)源.溶質(zhì)遷移在豎直方向上分別取5個(gè)觀測(cè)點(diǎn)，1#(0.5，0.425)，2#(0.5，0.4)，3#(0.5，0.35)，4#(0.5，0.25)，5#(0.5，0.1).

圖1 非飽和土壤中溶質(zhì)遷移模型示意圖Fig.1 Schematic diagram of physicalmodelw ith solute transport in unsaturated soils

1.2 數(shù)學(xué)模型

由于非飽和土壤中滲流、傳熱及溶質(zhì)傳遞是一個(gè)多因素互相耦合的復(fù)雜過程，為便于模型建立與求解，現(xiàn)作如下簡化假設(shè):①土壤水不可壓縮;②土壤土質(zhì)均勻，各向同性;③不存在氣體、固體和其他補(bǔ)給及化學(xué)反應(yīng)來源;④水體的密度、導(dǎo)熱性和比熱容恒定;⑤忽略溫度變化對(duì)土壤熱物理性質(zhì)的影響.

1)滲流控制方程.非飽和土壤滲流滿足Darcy定律［8］:

式中v為滲流速度;k為滲透系數(shù);H為水力水頭.

基于質(zhì)量守恒原理，非飽和土壤二維滲流的非等溫Richard方程［9］為

式中θ為土壤體積含水率;τ為時(shí)間;kx，ky分別為土壤x與y方向的等溫滲透系數(shù);kTx，kTy分別為土壤x與y方向的非等溫滲透系數(shù);T為溫度.

非飽和土壤中體積含水率和土壤水吸力呈非線性關(guān)系，可由土壤物理學(xué)領(lǐng)域最為普遍使用的van Genuchten土—水特征曲線通用方程［10］

表示，式中θ為土壤體積含水率;θr為殘余體積含水率;θs為飽和體積含水率;ψ為土壤水吸力;n，m，α為模型參數(shù)，其中n為土—水特征曲線指數(shù)，與曲線的陡緩程度有關(guān)且滿足n=(1－m)－1.

將van Genuchten土—水特征曲線與土壤滲透系數(shù)Mualem模型［11］聯(lián)合，可得土壤滲透系數(shù)的解析表達(dá)式:

式中Se=(θ－θr)/(θs－θr)為有效飽和度;kr=k/ks為相對(duì)滲透系數(shù)，其中ks為飽和滲透系數(shù);m為模型參數(shù).

2)傳熱控制方程.當(dāng)土壤中的滲流為穩(wěn)定流時(shí)，熱傳輸方程［8］330為

式中cs為土壤比熱容;λ為土壤熱導(dǎo)率;cw為土壤中流體(水)的比熱容;u為滲透速度.

3)溶質(zhì)傳遞方程.基于溶質(zhì)的質(zhì)量守恒定律，可得非飽和土壤的二維溶質(zhì)傳遞方程［8］231:

4)定解條件的確定.

①滲流問題.對(duì)于非飽和土壤模型，其頂部CD為速度入口邊界條件;AB為常水頭邊界條件，水頭為0;其余邊界為無流動(dòng)邊界.非飽和土壤ABCD區(qū)域，水頭壓力初始值為負(fù)數(shù).

②傳熱問題.非飽和土壤模型的AB和CD為溫度邊界條件，分別為300 K和308 K，其余邊界為周期性熱條件邊界.非飽和土壤ABCD區(qū)域初始溫度為邊界條件下穩(wěn)定的溫度場(chǎng)，如圖2所示.

③溶質(zhì)傳遞問題.非飽和土壤模型中點(diǎn)E為溶質(zhì)源，通量邊界條件為 0.5 mol·(m2·s)－1，其余為無通量邊界條件.整個(gè)土壤模型區(qū)域水中溶質(zhì)的初始濃度為0.

圖2 非飽和土壤溶質(zhì)遷移前穩(wěn)態(tài)溫度場(chǎng)(K)Fig.2 Tem perature distribution before solute transport in unsaturated soils

1.3 計(jì)算參數(shù)及數(shù)值計(jì)算

采用Comsol Multiphysics多物理場(chǎng)耦合分析軟件中多孔介質(zhì)地下水流模塊、傳熱模塊及物質(zhì)傳遞模塊，進(jìn)行3個(gè)物理場(chǎng)的耦合模擬計(jì)算.非飽和土壤滲流計(jì)算參數(shù)由土壤水特性軟件Spaw得到:滲透系數(shù)k為3.5 ×10－6m·s－1，飽和含水率 θs為0.43，殘余含水率 θr為 0.08，土壤密度 ρs為 1 400 kg·m－3，關(guān)系常數(shù)α，n分別為1.74和1.38;傳熱計(jì)算參數(shù):土壤熱導(dǎo)率 λs為 2 W·(m·K)－1，土壤比熱容 cs為 2 000 J·(kg·K)－1，水的密度 ρw為1 000 kg·m－3，水的熱導(dǎo)率 λw為0.54W·(m·K)－1，水的比熱容 cw為4 200 J·(kg·K)－1;溶質(zhì)傳遞計(jì)算參數(shù):分子擴(kuò)散系數(shù)D為8.25×10－7m2·s－1，水平彌散度αr為0.001 m，垂直彌散度αz為0.002 m.

本文采用三角形網(wǎng)格對(duì)非飽和土壤進(jìn)行自由網(wǎng)格劃分，近溶質(zhì)源區(qū)域網(wǎng)格較密，并進(jìn)行了網(wǎng)格精度檢驗(yàn).由于基本假設(shè)中忽略了溫度變化對(duì)土壤熱物理性質(zhì)的影響，且水體的物性參數(shù)恒定，所以實(shí)際大氣溫度波動(dòng)對(duì)溶質(zhì)遷移的影響較小，可忽略不計(jì).恒定的環(huán)境溫度條件下，土壤具有穩(wěn)定且呈線性分布的溫度場(chǎng)，形成滲流條件下溶質(zhì)遷移模型的初始溫度梯度，見圖2.

2 計(jì)算結(jié)果及分析

2.1 不同深度處土壤溫度及溶質(zhì)的物質(zhì)的量濃度分布

對(duì)入滲速度為2×10－6m·s－1時(shí)非飽和土壤中溶質(zhì)遷移模型進(jìn)行數(shù)值計(jì)算，分別得到第25 h時(shí)刻非飽和土壤溫度及溶質(zhì)的物質(zhì)的量濃度分布云圖、不同深度處非飽和土壤溫度及溶質(zhì)的物質(zhì)的量濃度變化曲線，見圖3～6.

圖3 第25 h時(shí)刻土壤溫度(K)分布云圖Fig.3 Tem perature distribution cloud graph at 25 h

圖4 第25 h時(shí)刻溶質(zhì)的物質(zhì)的量濃度分布云圖Fig.4 Concentration distribution cloud graph at 25 h

由圖3可觀察得到，土壤溫度由上至下逐步降低，與土壤初始時(shí)的溫度相比有了較大的提升.這主要是由于入滲雨水溫度高于土壤本身的初始溫度，隨著雨水滲流進(jìn)入土壤，與土壤固體骨架及水分進(jìn)行導(dǎo)熱和對(duì)流換熱，從而導(dǎo)致溫度逐漸升高.由圖4可以觀察得到，靠近溶質(zhì)源區(qū)域溶質(zhì)的c(A)越高，遠(yuǎn)離溶質(zhì)源區(qū)域c(A)越低，且溶質(zhì)沿豎直方向的遷移較水平方向更大.這是由于溶質(zhì)在土壤中的遷移與溶質(zhì)彌散度、水流在土壤中的滲透速度相關(guān)，且豎直方向的溶質(zhì)彌散度及水流的滲透速度均大于水平方向.

圖5 不同深度處土壤溫度變化曲線Fig.5 Variations of soil tem perature at different depths

圖6 不同深度處溶質(zhì)的物質(zhì)的量濃度變化曲線Fig.6 Variations of concentration at different depths

由圖5可觀察得到，1#處土壤溫度最高，5#處最低，且整體隨著時(shí)間增加而升高.在0～5 h期間，土壤溫度有一定波動(dòng)，越深處波動(dòng)越明顯.這是由于滲流前期土壤逐漸飽和，下層土壤含水率增大導(dǎo)致土壤熱導(dǎo)率減小，這與上層土壤熱導(dǎo)率差異較大;滲流后期土壤基本飽和，豎直方向熱導(dǎo)率差異逐漸縮小.由圖6可觀察得到，溶質(zhì)的c(A)隨時(shí)間變化而不斷增大，越接近溶質(zhì)源區(qū)域溶質(zhì)的c(A)越高且增幅更明顯.這是由于越接近土壤地表，溶質(zhì)的c(A)變化受滲流影響越顯著，且隨著深度增加影響越小.

2.2 不同入滲速度時(shí)土壤溫度及溶質(zhì)的物質(zhì)的量濃度分布

為分析入滲速度對(duì)非飽和土壤溫度及溶質(zhì)的物質(zhì)的量濃度分布特性的影響，對(duì)入滲速度分別為1×10－6，2×10－6，3×10－6m·s－1時(shí)溶質(zhì)遷移模型進(jìn)行數(shù)值計(jì)算，分別得到不同入滲速度時(shí)非飽和土壤中1#處的溫度及溶質(zhì)的物質(zhì)的量濃度變化曲線，見圖7～8.

圖7 不同入滲速度時(shí)1#處土壤溫度變化曲線Fig.7 Variations of soil tem perature at 1#for different infiltration rates

圖8 不同入滲速度時(shí)1#處溶質(zhì)的物質(zhì)的量濃度變化曲線Fig.8 Variations of concentration at 1#for different infiltration rates

由圖7觀察得到，土壤溫度隨時(shí)間變化大致逐漸升高，且入滲速度越大土壤溫升越快;當(dāng)v=1×10－6m·s－1時(shí)，土壤溫度波動(dòng)明顯.這是因?yàn)闈B流初期，入滲速度越小土壤的飽和時(shí)間越長，導(dǎo)致上、下層土壤之間的熱導(dǎo)率差異縮減緩慢，從而使得土壤溫度產(chǎn)生波動(dòng).由圖8可觀察得到，對(duì)于同一個(gè)觀測(cè)點(diǎn)1#，在土壤滲流初期，滲流速度越小1#處溶質(zhì)的c(A)越大;滲流后期，入滲速度越大溶質(zhì)的c(A)越大.這主要是由于滲流初期土壤非飽和，溶質(zhì)的對(duì)流作用大于彌散作用，當(dāng)入滲速度加快時(shí)溶質(zhì)對(duì)流作用越明顯，所以溶質(zhì)遷移越快，溶質(zhì)的c(A)增加緩慢;而滲流后期土壤達(dá)到飽和，此時(shí)彌散作用起主導(dǎo)作用，入滲速度越大水動(dòng)力彌散系數(shù)越大，溶質(zhì)的c(A)增加越快.綜合圖7，8可以得到，入滲速度的增加可以促進(jìn)土壤中熱量的擴(kuò)散和對(duì)流，增強(qiáng)溶質(zhì)的對(duì)流作用和彌散作用.

3 結(jié)論

1)滲流初期，上、下層土壤含水率不同導(dǎo)致土壤熱導(dǎo)率不均勻，對(duì)土壤溫度變化具有一定影響.

2)越接近土壤地表，溶質(zhì)的物質(zhì)的量濃度變化受滲流影響越顯著，且隨著深度增加，影響越來越小.

3)入滲速度的增加可以促進(jìn)土壤中熱量的擴(kuò)散和對(duì)流，增強(qiáng)溶質(zhì)的對(duì)流作用和彌散作用.

［1］ VIOTTIP，PAPINIM P，STRACQUALURSIN，et al.Contaminant transport in an unsaturated soil:laboratory tests and numerical simulationmodel as procedure for parameters evaluation ［J］.Ecol Model，2005，182(2):131-148.

［2］ KUMAR R P，DODAGOUDAR G R.Meshfree analysis of two-dimensional contaminant transport through unsaturated porousmedia using EFGM［J］.Int JNumer Meth Biomed Engng，2010，26(12):1797-1816.

［3］ JAVADIA A，AL-NAJJAR M M，EVANSB.Numericalmodeling of contaminant transport through soils:case study［J］.JGeotech Geoenviron Eng，2008，134(2):214-230.

［4］ GHASEMZADEH H.Heat and contaminant transport in unsaturated soil［J］.Int JCiv Eng，2008，6(2):90-107.

［5］ WINIARSKIT，LASSABATERE L，ANGULO-JARAMILLO R，et al.Characterization of the heterogeneous flow and pollutant transfer in the unsaturated zone in the fluvio-glacial deposit［J］.Procedia Environ Sci，2013，19:955-964.

［6］ 眭素剛，徐世光，劉文連，等.某渣場(chǎng)污染物在非飽和巖土介質(zhì)中的遷移模擬研究［J］.中南大學(xué)學(xué)報(bào)(自然科學(xué)版)，2013，44(5):2173-2180.

［7］ 左自波，張璐璐，王建華.降雨條件下非飽和土中污染物遷移的數(shù)值模擬［J］.地下空間與工程學(xué)報(bào)，2011，7(S1):1347-1352.

［8］仵彥卿.多孔介質(zhì)滲流與污染物遷移數(shù)學(xué)模型［M］.北京:科學(xué)出版社，2011:188-339.

［9］王華軍，路俊超，楊賓，等.溫度梯度下土壤鹽分遷移過程的數(shù)值模擬［J］.太陽能學(xué)報(bào)，2014，35(6):1086-1091.

［10］ van GENUCHTEN M T.A closed form equation for predicting the hydraulic conductivity of unsaturated soils［J］.Soil Sci Soc Am J，1980，44(5):892-898.

［11］ MUALEM Y.A new model for predicting the hydraulic conductivity of unsaturated porousmedia［J］.Water Resour Res，1976，12(3):513-522.

Numerical simulation of solute transport in unsaturated soils under tem perature gradient

SUN Lulu，YANGWeibo*

(Sch of Hydr，Energy ＆ Power Engin，Yangzhou Univ，Yangzhou 225127，China)

A model for solute transport in unsaturated soils is established based on the coupled governing equations of non-isothermal flow，heat transfer and solute transport.Numerical analyses are performed to investigate the effects of infiltration rates on the soil temperature variations and solute transport rules under temperature gradient.The results prove that soil temperature and solute concentration of different depths are mainly increasing over time，and both gradually reduce along the depth direction.And furthermore，the increase of the infiltration rate could promote the convection diffusion of heat in the soil and strengthen the advection and dispersion effect of solute.

unsaturated soils;temperature gradient;seepage;solute transport;numerical simulation

TK 124

A

1007-824X(2015)03-0046-05

2015-03-27.* 聯(lián)系人，E-mail:yangwb2004@163.com.

中國科學(xué)院可再生能源重點(diǎn)實(shí)驗(yàn)室資助項(xiàng)目(y507k51001);江蘇省自然科學(xué)基金資助項(xiàng)目(BK20141278);熱流科學(xué)與工程教育部重點(diǎn)實(shí)驗(yàn)室(西安交通大學(xué))開放基金資助項(xiàng)目(KLTFSE2014KF05);廣西建筑新能源與節(jié)能重點(diǎn)實(shí)驗(yàn)室開放基金資助項(xiàng)目(桂林能15-J-22-3).

孫露露，楊衛(wèi)波.溫度梯度作用下非飽和土壤中溶質(zhì)遷移過程的數(shù)值模擬［J］.揚(yáng)州大學(xué)學(xué)報(bào)(自然科學(xué)版)，2015，18(3):46-50.

(責(zé)任編輯 賈慧鳴，青 禾)