張 文, 丁子榮, 沈 啟, 王海濤， 阮周生

(1.東華理工大學 理學院,江西 南昌 330013; 2.合肥鐵路工程學校,安徽 合肥 230000;3.東華理工大學 核應用技術研究所,江西 南昌 330013)

砂巖型鈾礦是我國主要的鈾礦床類型之一，研究深部砂巖型鈾礦原地浸出的溶質遷移滲流規(guī)律，對于豐富和提高地浸采鈾的基礎理論認識和技術發(fā)展具有重要的理論意義與應用價值。

由于砂巖型鈾礦儲層的滲透系數較低，導致固液及液液界面的作用力較大，溶浸液很難遷移滲流。針對低滲透多孔介質滲流理論方面，許多學者開展了不同程度的研究并取得了一定成果(Mudd，2001；魏恒等，2013；闕為民等，2002；尹升華等，2006，2008，2013；李超等，2007；曾晟等，2011；Mellado et al.，2009；趙賀永，2015；李衡等，2019；Zhang et al.,2016；焦養(yǎng)泉等，2020；何潤發(fā)等，2020)。溶質在多孔介質中的遷移一般以對流與機械擴散為主，基于Fick擴散定律的對流擴散方程是描述溶質運移行為模式的經典模型(Zheng et al.，2002)，即溶質濃度在遷移過程中滿足的動力學模型為：

(1)

式中，u(x,t)表示在空間位置x處、t時刻的溶質濃度，k為擴散系數(k>0)，v為對流系數(v>0)，f(x,t)為源項，f(x,t)>0時為源，f(x,t)<0時則為匯。然而，用對流擴散方程模型來描述溶質的運移過程時，并不能很好地模擬穿透曲線的提前穿透或拖尾現象。若用微觀粒子的擴散過程來解釋這種不足，則表現為Fick定理描述的是微觀粒子向相鄰一個單位空間的擴散，但是在各向異性多孔介質中，微觀粒子的擴散可能向相鄰多個單位空間擴散。從隨機過程的角度去看，經典擴散為微觀粒子的局域性運動，是中心極限定理的直接結果，遷移過程遵守Fick定律，統(tǒng)計規(guī)律表現為均方位移與時間呈線性依賴關系；反常擴散現象表現為微觀粒子的非局域性(時間和空間)運動，是一種復雜系統(tǒng)的擴散過程，遷移過程不再遵守Fick定律，均方位移與時間呈現出非線性的依賴關系。描述反常擴散現象的方法和模型主要有分數布朗運動、廣義的擴散方程、連續(xù)時間隨機游走模型等。近年來，在很多復雜系統(tǒng)中均觀察到這種反?，F象，如金融市場以及多孔介質中的擴散等(張繼偉，2021；Podlubny,1999)。

針對不滿足布朗運動的反常擴散行為，通過Eulerian推導提出分數階Fick定理(Schumer et al.，2001)。在此基礎上，關于流體在多孔介質中遷移的分數階對流擴散方程也被推導出來，其能夠很好地描述分形幾何、冪律現象及記憶過程等反常擴散現象(Wheatcraft et al.，2008)。由于分數階微積分適合描述反常擴散過程，在各種科學和工程領域得到了關注，但其非局部性質的缺點導致很難計算出分數微分方程的解析解，于是尋求更簡單高效的數值方法來求解分數階微分方程成為眾多學者研究分數階微分方程的核心問題之一。筆者結合Caputo-Fabrizio時間分數階導數與對流擴散方程，研究深部砂巖型鈾礦地浸過程中溶質遷移的數學模型與數值解法。

1 數學模型

假設含鈾溶液的濃度u(x,y,z,t)是一個與三維空間(x,y,z)和時間都有關的量，考慮含鈾溶液在砂巖型鈾礦地浸遷移的機理。根據質量守恒定律，小區(qū)域中濃度變化由物質流入流出引起，用連續(xù)性方程的微分形式表示為：

(2)

F=-ku

(3)

(4)

當考慮流體中的靜水壓力時，由外力、重力、浮力、毛細管力等共同作用引起的壓力差會導致流體介質整體的流動，對濃度變化也將做出通量貢獻。根據質量守恒定律可得對流擴散方程：

(5)

(6)

(7)

為使得表達形式更為簡潔，考慮空間方向為一維形式，則：

(8)

(9)

(10)

2 時間分數階導數的數值計算

2.1 Caputo-Fabrizio時間分數階導數的計算格式

根據Caputo-Fabrizio導數的定義(Caputo et al.，2016)

(11)

式中，a表示區(qū)間[a,t]的起始端點滿足相容性條件，即：

當α→0+時有：

(12)

以及當α→1-時有：

(13)

(14)

令

(15)

(16)

(17)

(18)

令

(19)

易知hj-i(t)為一個r次多項式，且滿足關系

(20)

于是有

(21)

令

(22)

(23)

則

(24)

(25)

且逼近誤差為：

(26)

(27)

[f(t0)，f(t1)，…，f(tr-1)]T

(28)

[f(t0)，f(t1)，…，f(tn)]T

(29)

于是有

(30)

若記系數矩陣

In·Ωn,r=Gn=(gn,gn-1,…,g0)(n+1)

(31)

(32)

其中

2.2 Caputo-Fabrizio時間分數階導數的數值模擬

為了直觀展示分數階導數數值計算格式(32)的精確性，有利與Li等(2016)的計算結果進行比較，下述2個算例函數均選自Li等(2016)。不同的是， Li等(2016)中討論的為奇異核函數(t-s)-α分數階導數，采用了更為復雜的4次或5次多項式以達到4階或5階精度。

算例1，選取函數f1(t)=t6，易知分數階導數的精確解為：

算例2，選取函數f2(t)=e2t-2t-2t2-

表1和表2分別給出了算例1和算例2的計算誤差和數值精度。以上2個算例表明，選取2次或3次Lagrange插值多項式進行逼近時，能達到3階或4階精度，因此在進行Caputo-Fabrizio時間分數階導數對流擴散方程初邊值問題的數值模擬時，采用2次Lagrange插值多項式進行模擬。

表1 算例1的數值結果

續(xù)表

表2 算例2數值結果

3 數學模型在地浸采鈾中的應用

將分數階導數的計算格式嵌入式(10)得出數學模型(10)的計算格式，并應用于我國西部某鈾礦采冶區(qū)，結合鈾礦采冶區(qū)的探測數據，進一步驗證數學模型(10)的有效性。

3.1 溶質遷移數學模型的離散格式

(33)

(34)

于是，式(10)中第一個方程的數值計算格式可表示為：

(35)

式中，

初始條件(m=0,1,…,M)：

(36)

邊界條件(n=0,1,…,N)：

(37)

最終，式(10)的離散數值格式為:

(38)

整理后便有(m=1,2,…,M-1,n=1,2,…,N-1)：

(39)

其矩陣形式為：

HU=F

(40)

式中，

由求解三對角方程的追趕法易解出U，進一步便能逐層求解出un(n=0,1,2,…,N)。

3.2 數學模型砂巖型鈾礦地浸中的應用

根據國際原子能機構的統(tǒng)計數據(World Nuclear Association，2021)，原位浸出模式(圖1)已在鈾礦開采中廣泛使用，如加拿大的麥克阿瑟河鈾礦、哈薩克斯坦的特爾庫特蒙庫姆鈾礦以及澳大利亞的奧林匹克壩鈾礦。適合原位浸出模式的礦床通常具有2個特點:①具有一定滲透性的疏松砂巖型鈾礦床；②礦層上、下層具有較穩(wěn)定的頂、底板隔水層。目前，地浸采鈾技術主要有酸法浸出和堿法浸出。以美國為首的西方國家基本采用以二氧化碳和氧氣為體系的堿法浸出工藝，俄羅斯、哈薩克斯坦等國家基本采用以硫酸為體系的酸法浸出工藝，我國根據砂巖型鈾礦不同的礦床地質與水文地質條件來選擇酸法或堿法浸出工藝。

我國西部某砂巖型鈾礦采冶區(qū)采用以硫酸為體系的酸法浸出工藝，在注液井和抽液井之間放置中子測井儀，可以獲取礦體中的剩余鈾含量的監(jiān)測數據。為方便對監(jiān)測數據進行了歸一化處理，下面通過式(10)模擬在一維空間下鈾元素的遷移過程。

選取源項函數：

(41)

則易知式(10)所對應的精確解為u(x,t)=exsinx,φ(x)=0,φ(t)=sint,ψ(t)=e·sint。根據監(jiān)測數據估計出環(huán)境參數: 分數階α為0.001，擴散系數k為10-8和對流系數v為10-6。

通過計算數值解與精確解的相對誤差:

可以看出，數值格式(40)計算穩(wěn)定、精度高，該數學模型能較好地模擬原位浸出問題。

4 結論

基于含鈾溶液在砂巖型鈾礦地浸遷移的機理，建立了溶質遷移的時間分數階對流擴散方程數學模型，利用數值方法推導了Caputo-Fabrizio分數階導數對流擴散方程精度較高的數值格式，并通過算例驗證了格式的高精度和高效率。將所建立的數學模型應用于我國西部某砂巖型鈾礦地浸采鈾中，結果顯示數學模型是可靠的、數值算法是高效的。由于監(jiān)測井孔代價高昂，監(jiān)測數據非常有限，導致反演環(huán)境參數(分數階α、擴散系數k和對流系數v)并不理想，下一步課題組將引進正則化等方法，針對如何準確反演環(huán)境參數、預測剩余鈾含量等方面展開研究。