一維對流-擴(kuò)散試驗(yàn)各種邊界條件及其統(tǒng)一形式解析解

張文杰，趙 培，賈文強(qiáng)

（上海大學(xué) 土木工程系，上海 200072）

1 引 言

壓實(shí)黏土襯墊或黏土類防滲屏障廣泛應(yīng)用于填埋場及其他污染控制領(lǐng)域，除壓實(shí)土料的滲透系數(shù)外，擴(kuò)散和吸附性能對污染防控效果都有顯著影響。室內(nèi)土柱對流-擴(kuò)散試驗(yàn)常用來研究各種黏土類屏障材料中污染物的遷移規(guī)律，通過對試驗(yàn)數(shù)據(jù)的擬合分析，可以確定擴(kuò)散系數(shù)和阻滯因子。獲取用于擬合的試驗(yàn)數(shù)據(jù)有兩種方法：一種是截面分析，即經(jīng)過一定試驗(yàn)歷時(shí)后切片測試，獲取溶質(zhì)濃度沿深度變化曲線，該方法的優(yōu)點(diǎn)是土柱不必完全擊穿，因此，試驗(yàn)歷時(shí)相對短；另一種是累積濃度法[1]，獲得累積出流濃度隨時(shí)間變化曲線，該法操作簡單，但通常需要使出流濃度達(dá)到源濃度以獲得完整的擊穿曲線，因此，試驗(yàn)歷時(shí)要長得多。試驗(yàn)數(shù)據(jù)的擬合多采用一維對流-擴(kuò)散解析解[2-4]，然而不同邊界條件對應(yīng)不同的解析解，采用何種邊界條件關(guān)系到數(shù)據(jù)擬合結(jié)果的正確性或者準(zhǔn)確程度。

Shackelford[5-6]使用半無限邊界解析解擬合土柱試驗(yàn)數(shù)據(jù)，認(rèn)為土柱高度范圍內(nèi)溶質(zhì)濃度可以用無限遠(yuǎn)處濃度梯度為0的解析解來求；Barry等[7-8]在對上文的討論中提出不同意見，認(rèn)為如果自由出流邊界的溶液收集方法導(dǎo)致出流不連續(xù)，土柱底部應(yīng)采用零濃度梯度或零濃度邊界，且二者等價(jià)[9]；Rabideau等[10]認(rèn)為，零濃度梯度邊界不適合低滲流速度的情形，若存在垂直于擴(kuò)散方向的水流，則應(yīng)采用零濃度邊界；余闖[11]、謝海建[12]等建立了污染物一維擴(kuò)散解析解并進(jìn)行驗(yàn)證，其下邊界均采用零濃度邊界；詹良通等[13]擬合土柱離心試驗(yàn)數(shù)據(jù)時(shí)使用了半無限邊界的解析解。由此可見，國內(nèi)外學(xué)者在進(jìn)行土柱試驗(yàn)數(shù)據(jù)擬合或驗(yàn)證解析解時(shí)采用了不同的邊界條件。邊界條件到底該如何選擇？使用不同邊界條件解析解進(jìn)行數(shù)據(jù)擬合產(chǎn)生的誤差是否可以忽略等問題尚需進(jìn)一步研究。

本文探討了土柱試驗(yàn)通常涉及到的邊界條件及其適用條件，給出了統(tǒng)一形式的邊界條件，并求得統(tǒng)一形式的有限厚度土層對流-擴(kuò)散-吸附解析解，通過算例對比了各邊界條件解析解計(jì)算結(jié)果的差異，分析了因采用不同邊界條件進(jìn)行數(shù)據(jù)擬合產(chǎn)生的誤差。

2 土柱試驗(yàn)的邊界條件

一維土柱試驗(yàn)中，因?qū)α骱蛿U(kuò)散引起的溶質(zhì)通量可表達(dá)為[5]

式中：J為土孔隙流體中溶質(zhì)通量；θ為體積含水率，對于飽和土θ等于孔隙率n；v為孔隙水流速，通常認(rèn)為，v=V/θ，V為達(dá)西流速；cr為土孔隙流體中溶質(zhì)濃度；z為從土柱頂面開始算起的溶質(zhì)運(yùn)移距離（向下為正）；t為時(shí)間；D為水動(dòng)力彌散系數(shù)，是同時(shí)考慮分子擴(kuò)散和機(jī)械彌散的綜合參數(shù)。

考慮溶質(zhì)質(zhì)量守恒，由連續(xù)性條件知：

式中：Rd為考慮線性吸附的阻滯因子。

聯(lián)立式（1）、（2），即得到考慮對流、擴(kuò)散和吸附的土柱中溶質(zhì)一維運(yùn)移控制方程：

若忽略背景濃度，則初始條件為

為便于擬合分析，通常試驗(yàn)的土柱頂面采用溶液循環(huán)方法，或放置足夠體積的已知濃度溶液，以保證在試驗(yàn)過程中入流濃度不變或近似不變，因此，上邊界條件為

式中：0+為土柱上表面內(nèi)側(cè)的入流邊界。

底邊界條件較復(fù)雜，根據(jù)試驗(yàn)方法和底部出流條件，分為零濃度邊界、自由出流邊界、半無限邊界、柯西邊界4種。

2.1 零濃度邊界

當(dāng)土柱底部有與底面平行流動(dòng)的純水，流速和流量足夠時(shí)，流水將溶質(zhì)從土柱底面帶走的速度足夠快，或者土柱底面處的流動(dòng)不連續(xù)，此時(shí)底面邊界可設(shè)置為零濃度邊界：

式中：H為試樣厚度。

此時(shí)，底面孔隙中的溶質(zhì)能完全擴(kuò)散至邊界外，這一邊界條件也稱Dirichlet邊界，試驗(yàn)中土柱底部以純水循環(huán)可（接近）實(shí)現(xiàn)此邊界。在試驗(yàn)初期，溶質(zhì)遠(yuǎn)未遷移至底面附近，底面處濃度和濃度梯度均為 0，可用此邊界條件對應(yīng)的解析解擬合試驗(yàn)數(shù)據(jù)；而當(dāng)試驗(yàn)進(jìn)行足夠長時(shí)間后，此邊界條件對應(yīng)的底面處濃度梯度為（負(fù)）無窮大，而濃度仍為0，這與實(shí)際不符，因?yàn)橥林鶕舸┖蟮酌鎯?nèi)邊界會(huì)有少量溶質(zhì)積聚，特別是土柱內(nèi)溶質(zhì)遷移以對流為主導(dǎo)時(shí)，此時(shí)內(nèi)邊界濃度不為 0。因此，對應(yīng)的邊界條件應(yīng)表述為

式中：H+為土柱底面外邊界。

2.2 自由出流邊界

當(dāng)?shù)酌嫣巸H水壓力梯度不為 0，即擴(kuò)散作用相比對流來說可以忽略時(shí)，定義底部為自由出流邊界，此時(shí)出流溶質(zhì)的濃度梯度為0：

式中：cf為某時(shí)段內(nèi)按出流量計(jì)算的平均濃度，當(dāng)土柱內(nèi)溶質(zhì)遷移以對流為主導(dǎo)時(shí)，cf與cr接近，隨對流作用比例下降，cf與cr之間的差異增大，流量為0即純擴(kuò)散時(shí)cf無意義。用Péclet特征數(shù)Pe表示對流與擴(kuò)散在溶質(zhì)遷移中所占比例，Parlange等[14]認(rèn)為，當(dāng)Pe＞4時(shí)，近似有cf=cr，此時(shí)邊界條件為

這一邊界也稱為Neumann邊界。通常土柱試驗(yàn)在完全擊穿（以源濃度出流）之前，底面后出流的溶液濃度大于先前出流的溶液，因此，濃度梯度不為0，但若底邊界處孔隙中的溶液能被及時(shí)吸走，則底邊界處沒有濃度差，即不發(fā)生擴(kuò)散，此時(shí)可認(rèn)為濃度梯度為 0，故此邊界也被稱為完全吸收邊界，對應(yīng)一個(gè)理想狀態(tài)。通常粗顆粒土料的土柱試驗(yàn)較多采用這一邊界條件下的解析解進(jìn)行擬合，但試驗(yàn)時(shí)應(yīng)控制對流為主導(dǎo)，并注意底部濃度梯度，以接近這一邊界的上述物理本質(zhì)，而當(dāng)土料滲透系數(shù)很小，以擴(kuò)散為主導(dǎo)時(shí)，邊界處孔隙中的溶液不易排出，不宜采用此邊界。

2.3 半無限邊界

半無限邊界是指在無限深度處濃度梯度為0：

土柱試驗(yàn)中若土柱底部與收集到的出流溶液相接觸，且出流溶液的濃度測量采用下面兩種方式：一是出流達(dá)到一定體積后從底部移除并測量將之作為這段時(shí)間的平均出流濃度；二是持續(xù)測量累積濃度，在這兩種情況下，底邊界內(nèi)外濃度均不同，即濃度梯度不為 0。邊界條件式（10）對應(yīng)的解析解在z=H處的濃度梯度不為0，取某時(shí)刻t的濃度計(jì)算值可擬合濃度剖面，取z=H處的濃度計(jì)算值可擬合平均出流濃度，或經(jīng)變換后用于擬合累積濃度[1]。雖然式（10）定義無限遠(yuǎn)處的邊界，但擬合時(shí)流速v已知，此定義與對流無關(guān)，擬合造成的誤差僅來源于底邊界的擴(kuò)散速度，即式（10）確定的底邊界處濃度梯度值與試驗(yàn)操作對應(yīng)的真實(shí)濃度梯度之間的誤差。

2.4 柯西邊界

由上述可知，當(dāng)土柱底部被擊穿時(shí)，無論零濃度邊界（對應(yīng)濃度梯度為無窮大）還是自由出流邊界（對應(yīng)濃度梯度為 0）均基于理想假定，實(shí)際上土柱試驗(yàn)在擊穿后底邊界處濃度不為 0，而濃度梯度一般處于0與負(fù)無窮大之間，且出流濃度和濃度梯度并非固定值，這樣的邊界稱為混合邊界，即Cauchy邊界：

式中：μ、λ為描述邊界條件的參數(shù)。λ正比于擴(kuò)散系數(shù)，通常不為 0，這樣，當(dāng)接觸參數(shù)μ趨向無限大，式（11）退化成零濃度邊界，當(dāng)μ=0，式（11）退化成零濃度梯度邊界，μ取其他值時(shí)表示底邊界處濃度和濃度梯度均不為 0，且濃度梯度值隨濃度而變化。因此，式（11）可以看做各邊界條件的統(tǒng)一形式。式中的兩個(gè)參數(shù)μ、λ實(shí)際上可用一個(gè)值μ/λ表示（單位為1/m），稱為Cauchy邊界參數(shù)。

3 統(tǒng)一形式邊界條件下的解析解

由式（3）～（5）、（11）定義的數(shù)學(xué)模型無法直接求解，可通過變換式

Péclet數(shù)Pe常對有限域級數(shù)解的斂散性有較大影響，滿足收斂所需的項(xiàng)數(shù)m隨Pe增大而增大，依據(jù)此級數(shù)形式解析解編制簡單程序，其中m取20 000，可計(jì)算Pe=30對流為主導(dǎo)的問題。

4 不同邊界條件解析解對比及擬合誤差分析

式（21）對應(yīng)的邊界條件是前述各種邊界條件的統(tǒng)一形式，因此，可以方便地對比不同邊界條件下解析解的差別。算例中取土柱高度H=5 cm；孔隙水流速v= 3 × 1 0-6cm/s；水動(dòng)力彌散系數(shù)D=5×10-6cm2/s；忽略吸附作用，取Rd=1；試樣頂部為常濃度邊界，底部分別取上述各種邊界；模擬溶質(zhì)運(yùn)移時(shí)間分別為2、9 d。計(jì)算得到2、9 d后溶質(zhì)濃度分布剖面分別如圖1(a)、1(b)所示。

圖1 不同底部邊界條件解析解對比Fig.1 Analytical solutions for different boundary conditions

圖1(a)為溶質(zhì)運(yùn)移2 d后的濃度剖面線，可見土柱未被擊穿（擊穿指試樣底部出流達(dá)到某個(gè)濃度值或相對濃度值），圖中各剖面線重合，說明此時(shí)各種邊界條件下解析解結(jié)果相同。此時(shí)若使用解析解擬合試驗(yàn)數(shù)據(jù)，使用各邊界條件均可，不存在差別。但需要指出的是，為了獲得盡可能多的有效（非零的）數(shù)據(jù)點(diǎn)，以減小數(shù)據(jù)擬合誤差，土柱試驗(yàn)通常進(jìn)行到底部達(dá)到一定擊穿濃度才終止。

圖1(b)為溶質(zhì)運(yùn)移9 d后的濃度剖面線，此時(shí)土柱已被擊穿，可見在試樣上半部分各邊界對應(yīng)解析解的濃度基本接近，在試樣下半部分各解差異較大，其中零濃度邊界對應(yīng)的底部濃度最低，零濃度梯度邊界對應(yīng)的濃度最高，半無限邊界、柯西參數(shù)μ/λ取10和100的邊界對應(yīng)的濃度居中。如前所述，零濃度邊界假定底部能充分?jǐn)U散，與式（21）中μ/λ=∞的解相對應(yīng)；半無限邊界假定底邊界處存在有限的擴(kuò)散作用，本算例中與μ/λ=57的解答一致；零濃度梯度邊界假定擴(kuò)散為0，與μ/λ=0的解相對應(yīng)。因此，計(jì)算所得底部邊界處的濃度隨底部濃度梯度的增大而減小，也隨柯西邊界參數(shù)μ/λ的增大而減??；μ/λ越大，代表底邊界處擴(kuò)散在溶質(zhì)遷移中所占比例越大，亦即底部濃度梯度越大。

試驗(yàn)采取何種邊界條件解析解擬合，應(yīng)根據(jù)試驗(yàn)時(shí)底部邊界處的具體情況，即綜合考慮對流-擴(kuò)散比例和濃度梯度情況來確定。若底部保持為純水，則擊穿后濃度梯度很大，接近零濃度邊界的情形；若擴(kuò)散為主導(dǎo)，土樣底部接觸出流溶液，測出流溶液的累積濃度，此時(shí)邊界濃度梯度大于下為無限土層的情形，則采用μ/λ取較大值（如圖1中取100）的邊界；若土樣底部接觸出流溶液，出流溶液移除頻率較高，則對應(yīng)半無限邊界；若土樣底部不接觸出流溶液，土的滲透性較好，底部孔隙中的溶液不積聚，則采用μ/λ取較小值（如圖1中取10）的邊界；若底部邊界處擴(kuò)散可以忽略，即符合自由出流條件，則采用零濃度梯度邊界。

如某土柱試驗(yàn)為黏土試樣，土樣底部為透水石接出流管，則屬前述底部接觸出流溶液且出流溶液能及時(shí)移走的情形，對應(yīng)（或接近）半無限邊界。此處目的在于分析使用不同邊界條件解析解擬合造成的誤差，為了排除試驗(yàn)操作及數(shù)據(jù)本身離散性的誤差影響，以半無限邊界解析解數(shù)據(jù)為基準(zhǔn)。取土柱高度H=5 cm，采用惰性溶質(zhì)（取Rd=1），已知孔隙水流速v= 1 .7× 1 0-6cm/s，試驗(yàn)20 d后的溶質(zhì)濃度分布剖面如圖2中圓點(diǎn)所示，要擬合求水動(dòng)力彌散系數(shù)D。擬合需滿足所有基準(zhǔn)點(diǎn)處濃度擬合值與基準(zhǔn)值之差的平方和最小，編制簡單程序采用步進(jìn)搜索方法[15]求得最佳擬合值。使用無限邊界的解析解擬合得到D=1.45 × 1 0-6cm2/s（基準(zhǔn)解）。若用零濃度邊界擬合，得到D=1.67 × 1 0-6cm2/s，偏大15%，且擬合殘差較大；若用零濃度梯度邊界擬合，得到D=1.32×10-6cm2/s，雖然擬合殘差很小，但擬合值仍偏小 9%，誤差不容忽視。因此，根據(jù)試驗(yàn)的邊界條件合理選擇解析解進(jìn)行數(shù)據(jù)擬合，不僅是理論上的需要，也是數(shù)據(jù)處理準(zhǔn)確性的需要。

圖2 不同底部邊界條件解析解擬合結(jié)果Fig.2 Fitting curves of analytical solutions for different boundary conditions

5 結(jié) 論

（1）土柱試驗(yàn)數(shù)據(jù)可用零濃度邊界、零濃度梯度邊界、半無限邊界對應(yīng)的解析解來模擬，每一種邊界條件都與實(shí)際的試驗(yàn)條件相對應(yīng)。

（2）柯西邊界可看做各種邊界條件的統(tǒng)一形式，此邊界條件下的微分方程無法直接求解，通過引入變換式和輔助問題，求解得到適用于各種底邊界條件的有限厚度土層中考慮對流、擴(kuò)散、吸附的污染物運(yùn)移統(tǒng)一形式解析解。

（3）算例表明，土柱未被擊穿時(shí)各邊界條件解析解無差別；土柱被擊穿后，越靠近土柱底面，各種邊界條件解析解計(jì)算結(jié)果差異越大，底面處濃度值隨柯西參數(shù)μ/λ的增大而減小，μ/λ越大，代表底邊界處擴(kuò)散在溶質(zhì)遷移中所占比例越大，亦即底部濃度梯度越大，這是選擇何種邊界條件解析解進(jìn)行試驗(yàn)數(shù)據(jù)擬合的依據(jù)。

（4）擬合試驗(yàn)數(shù)據(jù)求水動(dòng)力彌散系數(shù)時(shí)，根據(jù)試驗(yàn)時(shí)土柱底部的對流-擴(kuò)散比例和濃度梯度具體情況，選擇相應(yīng)的邊界條件解析解，不僅是理論上的需要，也是顯著減小誤差的需要。

[1] SHACKELFORD C D. Cumulative mass approach for column testing[J]. Journal of Geotechnical Engineering,1995, 121(10): 696－703.

[2] PARKER J C, VAN GENUCHTEN M T. Flux-averaged and volume-averaged concentrations in continuum approaches to solute transport[J]. Water Resources Research, 1984, 20(7): 866－872.

[3] SHACKELFORD C D. Transit-time design of earthen barriers[J]. Engineering Geology, 1990, 29(1): 79－94.

[4] ACAR Y B, HAIDER L. Transport of low-concentration contaminants in saturated earthen barriers[J]. Journal of Geotechnical Engineering, 1990, 116(7): 1031－1052.

[5] SHACKELFORD C D. Critical concepts for column testing[J]. Journal of Geotechnical Engineering, 1994,120(10): 1804－1828.

[6] SHACKELFORD C D, PATRICK L R. Solute break through curves for processed kaolin at low flow rates[J].Journal of Geotechnical Engineering, 1995, 121(1): 17－32.

[7] BARRY D A, ANDERSON S J. Discussion of ‘Critical concept for column testing'[J]. Journal of Geotechnical Engineering, 1996,122 (1): 84－87.

[8] BARRY D A, ANDERSON S J. Discussion of ‘Solute break through curves for processed kaolin at low flow rates'[J]. Journal of Geotechnical Engineering, 1996,122 (9): 781－785

[9] BARRY D A, SPOSITO G. Application of the convection-dispersion model to solute transport in finite soil columns[J]. Soil Science Society of American Journal, 1988, 52(1): 3－9.

[10] RABIDEAU A, KHANDELWAL A. Boundary condition for modeling transport in vertical barriers[J]. Journal of Environmental Engineering, 1998, 124(11): 1135－1139.

[11] 余闖, 方冬芳, 杜廣印, 等. 均質(zhì)土中考慮溶質(zhì)蛻變效應(yīng)的污染物遷移模型的解析解[J]. 東南大學(xué)學(xué)報(bào)(自然科學(xué)版), 2012, 42(6): 1206－1210.YU Chuang, FANG Dong-fang, DU Guang-yin, et al.Analytical solution on pollutant migration considering solute decay in homogeneous soil[J]. Journal of Southeast University ( Natural Science Edition) , 2012,42(6): 1206－1210.

[12] 謝海建, 唐曉武, 陳云敏, 等. 原始土層影響下成層介質(zhì)污染物一維擴(kuò)散模型[J]. 浙江大學(xué)學(xué)報(bào)(工學(xué)版),2006, 40(12): 2191－2195.XIE Hai-jian, TANG Xiao-wu, CHEN Yun-min, et al.One-dimensional model for contaminant diffusion through layered media[J]. Journal of Zhejiang University (Engineering Science), 2006, 40(12): 2191－2195.

[13] 詹良通, 曾興, 陳云敏. 惰性污染物在低滲透性黏土中彌散作用離心模擬的相似研究[C]//第七屆全國巖土工程物理模擬學(xué)術(shù)研討會(huì)論文集. 杭州: [s. n.], 2013: 131－137.

[14] PARLANGE J Y, STARR J L, VAN GENUCHTEN M T,et al. Exit condition for miscible displacement experiments[J]. Soil Science, 1992, 153(3): 165－171

[15] 張文杰, 賈文強(qiáng), 張改革, 等. 黏土-膨潤土屏障中氯離子對流擴(kuò)散規(guī)律研究[J]. 巖土工程學(xué)報(bào), 2013,35(11): 2076－2081.ZHANG Wen-jie, JIA Wen-qiang, ZHANG Gai-ge, et al.Research on advection and dispersion of Cl-in claybentonite barriers[J]. Chinese Journal of Geotechnical Engineering, 2013, 35(11): 2076－2081.