基于參數優化的非飽和尾礦砂滲流場解析解研究

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

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

1 研究背景

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

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

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

2 解析模型的建立

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

(1)

式中：θ為體積含水率，mm3/mm3；D(θ)為非飽和尾礦砂中水的擴散系數，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均為引入的無量綱待定系數，在應用模型之前，需要先將所有參數的取值確定；q為求解區域x方向的最大坐標值，mm；h為求解區域y方向的最大坐標值，mm。

3 解析模型的求解

3.1 解析解的得出

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

(5)

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

(6)

w|x=0=0

(7)

w|x=q=0

(8)

(9)

A1e-A2h

(10)

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

(11)

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

(12)

(13)

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

3.2 模型參數的計算

3.1節中雖得出了解析解的具體形式，但其中涉及共計11個待定系數，在使用模型之前必須得到所有待定系數的值，否則模型是不完整的。因此，在結合試驗數據的基礎上進行參數優化不失為一種簡單高效的方法。參數優化的方法很多，通過多方對比，最終決定采用混沌粒子群優化算法(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”為粒子的序數，本文取80個粒子；下標“j”為未知參量序數，如前述共有11個未知參數；t為迭代次數，本文在參考試算經驗的基礎上將迭代總次數取為1 000；c1、c2為兩個加速常數，取值為c1=c2=1.70，c1與c2之間必須無任何關聯；r1j(t)和r2j(t)為兩個在0～1之間取值的隨機數，每次迭代中均由計算機重新給出；ω為速度的慣性權數，取為1；vij(t+1)和xij(t+1)分別為t+1次迭代的速度和位置；vij(t)和xij(t)分別為t次迭代的速度和位置；pij(t)和pgj(t)分別為t次迭代的局部最優和整體最優。公式(16)為Logistic系統迭代公式，當初值z0滿足0≤z0≤1時，通過(16)式可形成一個混沌集合，使系統完全處于混沌狀態，避免計算陷入局部最優。公式(16)中μ為控制參量，取值為4。以上各參量均無量綱。

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

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

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

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

4 模型的有效性驗證

4.1 已有試驗結果及模型擬合效果

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

圖1 試驗裝置與求解區域

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

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

4.2 模型參數的計算結果

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

α=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 結 論

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

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