林雪松, 王永吉, 王 東, 王來貴, 王漢勍, 宋澤友, 王中東, 李俊岑, 馮 歡

(1.遼寧工程技術(shù)大學(xué) 理學(xué)院，遼寧 阜新 123000; 2.阜新高等?？茖W(xué)校，遼寧 阜新 123000;3.遼寧工程技術(shù)大學(xué) 力學(xué)與工程學(xué)院，遼寧 阜新 123000)

0 引言

尾礦壩一旦發(fā)生潰壩將對下游人民的生命安全和生存環(huán)境造成嚴(yán)重破壞[1]，因此一直以來都是安全生產(chǎn)領(lǐng)域的重點研究內(nèi)容。近年來，由于尾礦顆粒越來越細(xì)[2]，固結(jié)越來越困難，壩體潰決概率出現(xiàn)了不斷上升的趨勢，使得針對細(xì)尾礦砂堆筑的尾礦壩的研究開始吸引學(xué)者們的注意[3-4]。水的分布狀態(tài)是影響尾礦壩安全的重要因素，細(xì)尾礦砂堆筑的尾礦壩更是如此。尾礦壩中的水分有很大一部分是存在于非飽和尾礦砂中的，在尾礦壩的安全狀態(tài)評估中，必須考慮非飽和部分的影響[5]，而清楚地描述非飽和尾礦砂內(nèi)部水分的分布狀況是考慮該部分影響的重要前提，必須進(jìn)行深入探索。多年來國內(nèi)外學(xué)者針對非飽和滲流問題進(jìn)行過很多有意義的研究。張志軍等[6]通過尾礦砂室內(nèi)毛細(xì)水上升測試發(fā)現(xiàn)：尾礦砂中毛細(xì)水的上升規(guī)律呈初期迅速而后期緩慢的變化趨勢，是其受多種作用力的共同作用所致;魏樣等[7]通過室內(nèi)土柱試驗發(fā)現(xiàn)：不同土質(zhì)的污染土壤內(nèi)，毛細(xì)水上升高度與時間均符合良好的冪函數(shù)關(guān)系;郝瑞等[8]利用自行設(shè)計的層狀土室內(nèi)模型試驗監(jiān)測了毛細(xì)水上升過程，通過試驗結(jié)果對LW模型進(jìn)行了修正。前述工作可歸類為實驗方面的研究，在數(shù)值計算方面，Neuman[9]首先提出通過有限元方法解決二維飽和-非飽和滲流問題;湯卓等[10]基于非飽和土滲流理論建立了尾礦壩三維有限元模型，總結(jié)了上游法尾礦壩流固耦合計算規(guī)律；彭華等[11]對滲流問題的有限元分析方法加以改進(jìn)，消除了非飽和滲流計算中存在的數(shù)值彌散現(xiàn)象，提高了收斂速度 ;王鐵行[12]給出了考慮含水率和密度影響的非飽和黃土路基水分場數(shù)值計算模型;劉杰等[13]研究了非飽和土路基毛細(xì)作用的數(shù)值與解析方法;李毅等[14]運用變單元滲透系數(shù)法，通過FLAC3D的二次開發(fā)得到了飽和-非飽和滲流場;王睿等[15]通過引入時間分?jǐn)?shù)階導(dǎo)數(shù)，結(jié)合2種水力傳導(dǎo)系數(shù)對新的Richards方程進(jìn)行了數(shù)值求解。

縱上發(fā)現(xiàn)，前人在非飽和滲流問題的實驗和數(shù)值計算方面進(jìn)行了大量研究，但較少見到解析方法方面的研究，即便使用解析方法也僅限于處理一維問題。實驗和數(shù)值計算當(dāng)然是探討非飽和滲流問題的有效方法，應(yīng)該投入精力進(jìn)行研究與發(fā)展，但重視實驗和數(shù)值方法的同時也應(yīng)該重視解析方法，加之尾礦壩問題通常涉及2個維度，因此對二維問題的解析研究意義也同樣明顯。

本文擬從非飽和滲流的Richards方程出發(fā)，在對求解區(qū)域和邊界進(jìn)行適當(dāng)加工的基礎(chǔ)上，得到尾礦砂二維非飽和滲流的穩(wěn)態(tài)解析解，然后通過得到的解析解擬合細(xì)尾礦砂實驗數(shù)據(jù)，觀察解的擬合效果。主要目的是通過研究給出符合實際的非飽和細(xì)尾礦砂內(nèi)含水率分布計算的思路與方法，為考慮非飽和部分的尾礦壩安全評估工作提供基礎(chǔ)資料。

1 創(chuàng)建解析模型

Richards方程是細(xì)尾礦砂非飽和滲流所涉及的基本方程，若以體積含水率θ為待求解函數(shù)，坐標(biāo)系的z軸以豎直向上為正，三維問題中的Richards方程的形式如式(1)所示[16]。方程中的D(θ)稱為擴(kuò)散系數(shù)，D(θ)隨θ的變化而變化，因此表達(dá)為函數(shù)形式；K(θ)表示非飽和尾礦砂中水分的滲透系數(shù)，也表達(dá)成了θ的函數(shù)。在尾礦壩工程中通常研究的都是二維問題，此情況下，如將坐標(biāo)軸y以豎直向上為正方向，且僅僅研究穩(wěn)態(tài)問題不討論瞬態(tài)問題，式(1)可變化為式(2)。雖然式(2)的形式已經(jīng)比較簡單，但若想求得解析解尚需在符合實際的基礎(chǔ)上繼續(xù)簡化。接下來將D(θ)取為常數(shù)D，則式(2)又可化為式(3)。

(1)

(2)

(3)

需要指出的是，將D(θ)取為常數(shù)D并不是認(rèn)為實際中D(θ)一定為常數(shù)，而是在運算中用某種平均值來代替D(θ)，前人曾實踐過類似的簡化思路[13]。式(3)必須和一定的求解區(qū)域和邊界條件聯(lián)系在一起才能求解，若區(qū)域的形狀太過復(fù)雜是無法得到解析解的，因此應(yīng)當(dāng)將求解區(qū)域選為簡單形狀，在接下來的研究中，選擇矩形求解區(qū)域，因為矩形區(qū)域內(nèi)的方程求解是完全可以做到的，且實際中很多工程問題的求解也可歸結(jié)到矩形區(qū)域內(nèi)。接下來考慮式(3)的邊界條件。通過對已有實驗結(jié)果(論文后面會給出)的分析可看出，含水率在豎直和水平方向上都近似遵循指數(shù)變化的規(guī)律，有鑒于此，在矩形的求解區(qū)域內(nèi)，建立式(4)～(6)所示的方程和邊界條件。

(4)

θ|x=0=ae-byθ|x=q=ce-dy

(5)

θ|y=0=fe-gxθ|y=h=je-kx

(6)

2 求解模型的具體過程

方程和邊界條件建立起來之后，就涉及到求解整個模型。首先，方程和邊界條件都是非齊次的，可對邊界條件先進(jìn)行齊次化處理。令含水率θ=w+v，w和v為2個待定函數(shù)。假設(shè)v滿足邊界條件(5)，且令v=Ax+B，A和B同樣為2個待定函數(shù)。將v帶入邊界條件(5)解得：

(7)

從v的表達(dá)式中可看出：

(8)

將θ=w+v及v的表達(dá)式帶入方程和邊界條件(4)～(6)得到關(guān)于w的方程和邊界條件如式(9)～(11)所示：

(9)

w|x=0=0，w|x=q=0

(10)

(11)

利用分離變量和傅立葉級數(shù)法可得到w的表達(dá)式為：

(12)

公式中Cn，Dn的形式為：

(13)

(14)

則最終得到含水率空間分布的表達(dá)式為：

(15)

很明顯最終的結(jié)果具有無窮級數(shù)形式，但由于系數(shù)都具有很強(qiáng)的收斂性，所以在實際中只取其中前幾項就足夠了。

3 基于實驗結(jié)果的模型效果檢驗

3.1 實驗設(shè)計與結(jié)果

1)實驗裝置與物料

主要實驗裝置如圖1所示，包括儲水腔室(1 m×0.1 m×1 m)和儲料腔室(0.1 m×0.1 m×1 m)，尾礦物料全部放置在儲料腔室內(nèi)，圖1中儲料腔室上的格點表示計劃測含水率的位置。儲水腔室用于提供滲流過程中所需的水，2個腔室之間的隔板上設(shè)置有鉆孔，用來保證滲流過程中水供應(yīng)的順暢。實驗對象為某鐵礦尾礦砂，顆粒組成見表1，表1中數(shù)據(jù)顯示：50%以上顆粒粒徑小于0.019 mm，大于0.074 mm的顆粒含量低于10%，大于0.037 mm的顆粒含量低于30%，實驗用尾礦砂完全可歸類于細(xì)尾礦砂[17-18]?；蛘甙凑展こ谭诸惙绞綒w類為尾黏土。

2)實驗主要過程與參數(shù)設(shè)定

在實驗開始之前，先將尾礦砂填筑到儲料腔室內(nèi)，然后在儲水腔室中注入清水，當(dāng)腔室內(nèi)水位達(dá)到0.05 m時停止注水，并通過虹吸方式使水位保持恒定。接下來細(xì)尾礦砂開始在恒定水源水位的基礎(chǔ)上進(jìn)行飽-非飽和滲流。實驗中潤鋒的位置與形態(tài)不斷發(fā)生變化，當(dāng)潤鋒完全穩(wěn)定時認(rèn)為實驗結(jié)束。細(xì)尾礦砂滲流速度較慢，從實驗開始到潤鋒穩(wěn)定共計用了2個月時間。剛開始水平向滲流速度大于豎直向滲流速度很多，但隨著滲流的發(fā)展二者差距越來越小。穩(wěn)定后的狀態(tài)如圖1所示，圖1中AA′表示穩(wěn)定后的潤鋒，BB′表示實驗中尾礦砂試件內(nèi)出現(xiàn)的1個裂縫。滲流達(dá)到穩(wěn)態(tài)后對儲料腔室上標(biāo)記的各格點進(jìn)行挖掘，利用烘干法測含水率。該細(xì)尾礦砂飽和狀態(tài)體積含水率為49.80%。

表1 細(xì)粒尾礦砂粒徑組成Table 1 Composition of particle size of fine tailing sand

圖1 實驗結(jié)果Fig.1 Experimental results

3)實驗結(jié)果

測量得到的豎直向含水率分布曲線如圖2所示，水平向含水率分布曲線如圖3中的實測值所示。由圖2～3可知：無論豎直還是水平方向，含水率均逐漸降低，降低趨勢具明顯的非線性，可近似用指數(shù)函數(shù)描述。鑒于曲線的分布特點，前述模型中的邊界條件均采用指數(shù)形式。

圖3 水平向含水率分布Fig.3 Horizontal distribution of water content

3.2 模型參數(shù)的確定

前述解析解(15)中包含共計10個未知參數(shù)，未知參數(shù)完全可利用實驗數(shù)據(jù)通過數(shù)學(xué)方法得出?？蛇x的數(shù)學(xué)方法有很多，通過仔細(xì)調(diào)研，最終決定采用混沌粒子群優(yōu)化算法(Chaos Particle Swarm Optimization，CPSO)來計算參數(shù)。由于理論公式只能應(yīng)用在矩形區(qū)域內(nèi)，因此應(yīng)用公式的時候首先需要構(gòu)造矩形域。

為檢驗計算得到的參數(shù)對所有數(shù)據(jù)的預(yù)測結(jié)果如何，需將由理論模型計算得到的各點含水率與實際測量結(jié)果進(jìn)行比較，接下來展示出不同高度處含水率沿水平方向分布的實測值和理論值，以此將所有測量點處的實際測量值與理論計算值進(jìn)行比較。比較結(jié)果如圖3所示。如果仔細(xì)觀察圖3可看出，在高度0.1，0.3和0.4 m處，實測值與計算值符合的比較好，在高度0和0.2 m處，數(shù)據(jù)擬合結(jié)果相對要稍差一些，但相差的也不是很大，總體來說，模型計算值和實測值是比較接近的，證明了用推導(dǎo)出的公式來描述實驗數(shù)據(jù)是完全可行的。因此解析解可被認(rèn)為是含水率空間分布的理論規(guī)律，計算得到的分布曲線一般是要比測量得到的曲線稍顯平滑。

3.3 擬合結(jié)果

針對3.1節(jié)實驗構(gòu)造的矩形域如圖1中的COEF所示，在矩形域上分別建立了x軸和y軸，矩形域內(nèi)共計包含了50個數(shù)據(jù)點，可用來擬合參數(shù)。運算中在參考試算經(jīng)驗的基礎(chǔ)上，將最大迭代次數(shù)設(shè)置為1 000，迭代過程中誤差的變化如圖4所示，從圖4中看到，開始的時候收斂速度比較快，然后收斂速度逐漸變慢，隨著迭代次數(shù)的增加最后誤差趨于恒定值。由計算得到公式中的未知參數(shù)取值如下：

α=6.839 088 379,β=20 000.0,a=44.229 302 31,b=0.096 026 944 5,c=979.961 319 9,d=20 000.0，f=4 465.811 49,g=0.120 477 598 5,j=21.483 476 85,k=16 768.464 75,D=16 177.809 65。

圖4 誤差變化Fig.4 Error development

4 結(jié)論

1)通過解析方法得到了描述非飽和尾礦砂含水率空間分布的公式，求解過程表明，細(xì)尾礦砂非飽和滲流穩(wěn)態(tài)問題的解析求解可以從Richards方程的簡化開始，簡化主要沿著2個方向，第1是將方程形式簡化，簡化過程中包括降低維度，變瞬態(tài)為穩(wěn)態(tài)，將未知參數(shù)簡化為常數(shù)或固定函數(shù)形式等；第2是將復(fù)雜求解區(qū)域簡化為矩形，在矩形區(qū)域里只要給出合適的邊界條件是可以得到方程解析解的。

2)模型求解所需的邊界條件可在參考實驗數(shù)據(jù)的基礎(chǔ)上給出，以此來保證求解結(jié)果的有效性。簡化過程中引入的參數(shù)可通過CPSO算法編制程序利用測量數(shù)據(jù)得出。算法計算中誤差具有明顯的收斂性。得出的參數(shù)能夠很好地擬合全部測量數(shù)據(jù)。由計算得到的曲線較實測曲線具有明顯的光滑性。