林雪松， 王 東， 王來貴， 王永吉， 李俊岑

(1.遼寧工程技術(shù)大學 理學院，遼寧 阜新 123000; 2.阜新高等專科學校 工程系，遼寧 阜新 123000; 3.遼寧工程技術(shù)大學 力學與工程學院， 遼寧 阜新 123000)

1 研究背景

作為具有重大安全隱患的人造泥石流危險源，尾礦壩一旦發(fā)生潰決將對下游人民的生命財產(chǎn)和生存環(huán)境造成重大破壞[1-3]，因此，一直以來都是研究的重點。近年來，由于尾礦顆粒越來越細[4-5]，固結(jié)越來越困難，形成壩體潰決的概率呈不斷上升的趨勢[6]，該趨勢使尾礦壩安全問題引起許多學者的關(guān)注。水的分布狀態(tài)是影響尾礦壩安全的重要因素，尾礦壩中有相當一部分水分是存在于非飽和尾礦砂中的，在壩體的安全狀態(tài)評估中，非飽和部分內(nèi)部水分分布的影響必須予以考慮[7]，因此針對該部分水分分布的計算研究便具有重大意義。多年來國內(nèi)外學者針對非飽和滲流計算問題進行過很多有意義的研究。Neuman[8]首先提出通過有限元方法解決二維飽和-非飽和滲流問題，彭華等[9]對滲流問題的有限元分析方法加以改進，提出了修正容水度和加速迭代收斂技術(shù)的新方法，消除了非飽和滲流計算中存在的數(shù)值彌散現(xiàn)象，提高了收斂速度。王鐵行[10]考慮含水率和密度影響，給出了非飽和黃土路基水分場的數(shù)值計算模型，并探討和確定了模型參數(shù)。劉杰等[11]研究了非飽和土路基毛細作用的數(shù)值與解析方法。李毅等[12]運用變單元滲透系數(shù)法，對FLAC3D軟件進行二次開發(fā)，準確得出了連續(xù)光滑的自由水面，求解出了飽和-非飽和滲流場。劉豪杰等[13]通過自己建立的某深厚覆蓋層土石壩三維滲流有限元數(shù)值模型計算，進行了針對該土石壩滲流控制方案的合理優(yōu)化。劉夢茹等[14]運用非飽和滲流理論得到了包氣帶一維液相運動方程解析解。

縱觀以往的研究，非飽和滲流計算的研究重點均在數(shù)值計算方法方面，解析方法研究的較少，即便使用解析方法也僅限于處理一維問題。數(shù)值方法的實際工程意義是有目共睹的，但重視數(shù)值方法的同時也不應該完全忽視解析方法，同時還應重視解析法與其他方法的結(jié)合。

本文擬從非飽和滲流的Richards方程出發(fā)，得到二維非飽和尾礦砂穩(wěn)態(tài)滲流含水率的解析解，然后利用已有試驗數(shù)據(jù)，通過參數(shù)優(yōu)化的方法計算解析解中引入的待定系數(shù)，從而使解析解完整，最后對含水率計算值與實測值進行對比，以此考察解析解的實際應用效果。本文的主要目的是通過研究得出非飽和尾礦砂內(nèi)的水分分布規(guī)律以及含水率計算的思路與方法，為考慮非飽和部分的尾礦壩安全評估工作提供基礎(chǔ)資料。

2 解析模型的建立

Richards方程是解決非飽和滲流問題的基本工具，若以體積含水率θ為求解的主要參量，坐標系z軸以豎直向上為正，則三維問題中的Richards方程的形式如公式(1)所示。

(1)

式中：θ為體積含水率，mm3/mm3；D(θ)為非飽和尾礦砂中水的擴散系數(shù)，mm2/s；K(θ)為非飽和尾礦砂的導水率， mm/s。

(2)

θ|x=0=A1e-A2yθ|x=q=A3e-A4y

(3)

θ|y=0=A5e-A6xθ|y=h=A7e-A8x

(4)

式中：α、β和A1～A8均為引入的無量綱待定系數(shù)，在應用模型之前，需要先將所有參數(shù)的取值確定；q為求解區(qū)域x方向的最大坐標值，mm；h為求解區(qū)域y方向的最大坐標值，mm。

3 解析模型的求解

3.1 解析解的得出

令θ=w+v，v=Ax+B，w、v、A、B均為任意函數(shù)，假設v滿足邊界條件公式(3)，將v代入邊界條件公式(3)解得：

(5)

將θ=w+v代入方程和邊界條件公式(2)～(4)得到關(guān)于w的方程和邊界條件如公式(6)～(10)所示。

(6)

w|x=0=0

(7)

w|x=q=0

(8)

(9)

A1e-A2h

(10)

利用分離變量和傅立葉級數(shù)法可得到w，然后得到含水率空間θ分布的表達式為：

(11)

公式(11)中Cn、Dn的形式為：

(12)

(13)

很明顯最終的結(jié)果具有無窮級數(shù)的形式，但因系數(shù)均具有很強的收斂性，所以在實際應用當中只取其中前幾項即可。

3.2 模型參數(shù)的計算

3.1節(jié)中雖得出了解析解的具體形式，但其中涉及共計11個待定系數(shù)，在使用模型之前必須得到所有待定系數(shù)的值，否則模型是不完整的。因此，在結(jié)合試驗數(shù)據(jù)的基礎(chǔ)上進行參數(shù)優(yōu)化不失為一種簡單高效的方法。參數(shù)優(yōu)化的方法很多，通過多方對比，最終決定采用混沌粒子群優(yōu)化算法(Chaotic Particle Swarm Optimization, CPSO)。

3.2.1 運算的基本原理 算法的基本公式為：

vij(t+1)=ωvij(t)+c1r1j(t)(pij(t)-xij(t))+

c2r2j(t)(pgj(t)-xij(t))

(14)

xij(t+1)=xij(t)+vij(t+1)

(15)

zn+1=μzn(1-zn) (n=0,1,2,…)

(16)

式中：下標“i”為粒子的序數(shù)，本文取80個粒子；下標“j”為未知參量序數(shù)，如前述共有11個未知參數(shù)；t為迭代次數(shù)，本文在參考試算經(jīng)驗的基礎(chǔ)上將迭代總次數(shù)取為1 000；c1、c2為兩個加速常數(shù)，取值為c1=c2=1.70，c1與c2之間必須無任何關(guān)聯(lián)；r1j(t)和r2j(t)為兩個在0～1之間取值的隨機數(shù)，每次迭代中均由計算機重新給出；ω為速度的慣性權(quán)數(shù)，取為1；vij(t+1)和xij(t+1)分別為t+1次迭代的速度和位置；vij(t)和xij(t)分別為t次迭代的速度和位置；pij(t)和pgj(t)分別為t次迭代的局部最優(yōu)和整體最優(yōu)。公式(16)為Logistic系統(tǒng)迭代公式，當初值z0滿足0≤z0≤1時，通過(16)式可形成一個混沌集合，使系統(tǒng)完全處于混沌狀態(tài)，避免計算陷入局部最優(yōu)。公式(16)中μ為控制參量，取值為4。以上各參量均無量綱。

3.2.2 運算的基本過程 首先在位置、速度的取值范圍內(nèi)隨機給出二者的初值，將c1、c2、ω和最大迭代次數(shù)等參數(shù)賦值，然后進入循環(huán)。每個循環(huán)中的運算步驟與內(nèi)容如下：

(1)由公式(11)計算和對比分別得到局部最優(yōu)和整體最優(yōu)。

(2)將整體最優(yōu)從原解空間映射到[0,1]范圍內(nèi)，然后作為初值通過公式(16)迭代29次得到一組混沌向量，再將該組向量映射回原解空間，從中計算出最優(yōu)值，隨機替代當前粒子群中任意一個粒子。

(3)按公式(14)和(15)更新速度和位置，并替換速度和位置中超出取值范圍的數(shù)值，從第(1)步開始進行新的計算，直到達到最大迭代次數(shù)。

4 模型的有效性驗證

4.1 已有試驗結(jié)果及模型擬合效果

為計算模型參數(shù)和驗證模型的適用性，筆者進行了尾礦砂非飽和滲流試驗，獲得了非飽和尾礦砂穩(wěn)態(tài)含水率的分布結(jié)果與特點。試驗采用的尾礦物料粒徑很小，完全可稱為細尾礦砂[15-16]。試驗裝置正面及非飽和滲流穩(wěn)定后狀態(tài)如圖1所示，裝置主體部分為蓄水池(100 mm×100 mm×1 000 mm)和尾礦物料池(1 000 mm×100 mm×1 000 mm)。蓄水池的作用是為尾礦物料的非飽和滲流提供水頭恒定的水源，試驗中保持水頭高度50 mm不變。尾礦物料池用于盛裝干尾礦物料，物料池正面面板上設置的網(wǎng)格點是用來標記需測量含水率的位置的，在水平和豎直方向上，格點間距均為100 mm。圖1中AA′表示穩(wěn)定后的潤鋒，BB′表示實驗中尾礦砂試件內(nèi)出現(xiàn)的一個裂縫。從圖1中可看出AA′清晰的將試驗尾礦物料分為干尾礦砂和濕尾礦砂。圖1中矩形區(qū)域OCFE是一個位于濕尾礦砂內(nèi)部的矩形，可用來作為一個求解域，在矩形域上分別建立了x軸和y軸，矩形域內(nèi)共計包含了50個數(shù)據(jù)點，可用來擬合參數(shù)。矩形域長900 mm，寬400 mm，則前述公式(11)中q=900 mm，h=400 mm。測量得到的豎直向含水率分布曲線如圖2所示，水平向含水率分布曲線如圖3所示，由圖(2)～(3)可知，無論在豎直還是在水平方向，含水率均逐漸降低，降低的趨勢具有明顯的非線性，可以近似的用指數(shù)函數(shù)描述。鑒于曲線的分布特點，前述模型中的邊界條件均采用指數(shù)形式。在圖3中還給出了通過模型計算得到的各格點含水率，并將計算值與實測值進行了比較?？傮w來說，計算值與實測值是比較接近的，說明模型完全具有實際價值，另外，計算得到的分布曲線一般要比實測得到的曲線稍顯平滑。

圖1 試驗裝置與求解區(qū)域

圖2 不同水平距離含水率隨高度的變化規(guī)律

圖3 不同高度含水率實測、計算值在

4.2 模型參數(shù)的計算結(jié)果

4.1節(jié)中用模型計算各格點含水率之前，就已知模型中的參數(shù)值，所有參數(shù)值均利用3.2節(jié)所述的原理與過程，通過Matlab7.0編程計算得到，運算中誤差的變化如圖4所示，從圖4中看到，開始時收斂速度較快，然后逐漸變慢，隨著迭代次數(shù)的增加最后誤差趨于恒定值。由計算得到公式中的未知參數(shù)取值如下：

α=6.839088379,β=20000.0,

A1=44.22930231,A2=0.0960269445,

A3=979.9613199,A4=20000.0，

A5=4465.81149,A6=0.1204775985,

A7=21.48347685,A8=16768.46475,

D=16177.80965

圖4 運算誤差演化曲線

5 結(jié) 論

(1)通過解析方法得到了描述二維非飽和尾礦砂穩(wěn)態(tài)滲流含水率空間分布公式，求解過程從Richards方程的簡化開始，方程簡化主要包括兩個方面，第一是將方程形式簡化，簡化過程中包括降低維度、由瞬態(tài)變?yōu)榉€(wěn)態(tài)、將未知參數(shù)簡化為常數(shù)或固定函數(shù)形式等。第二是將復雜的求解區(qū)域簡化為矩形，在矩形區(qū)域里只要給出邊界條件是可以得到方程解析解的。

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