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

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

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

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

(1)

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

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

1 數(shù)學(xué)模型

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

(2)

F=-ku

(3)

(4)

當(dāng)考慮流體中的靜水壓力時(shí)，由外力、重力、浮力、毛細(xì)管力等共同作用引起的壓力差會(huì)導(dǎo)致流體介質(zhì)整體的流動(dòng)，對濃度變化也將做出通量貢獻(xiàn)。根據(jù)質(zhì)量守恒定律可得對流擴(kuò)散方程：

(5)

(6)

(7)

為使得表達(dá)形式更為簡潔，考慮空間方向?yàn)橐痪S形式，則：

(8)

(9)

(10)

2 時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)的數(shù)值計(jì)算

2.1 Caputo-Fabrizio時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)的計(jì)算格式

根據(jù)Caputo-Fabrizio導(dǎo)數(shù)的定義(Caputo et al.，2016)

(11)

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

當(dāng)α→0+時(shí)有：

(12)

以及當(dāng)α→1-時(shí)有：

(13)

(14)

令

(15)

(16)

(17)

(18)

令

(19)

易知hj-i(t)為一個(gè)r次多項(xiàng)式，且滿足關(guān)系

(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)

若記系數(shù)矩陣

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

(31)

(32)

其中

2.2 Caputo-Fabrizio時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)的數(shù)值模擬

為了直觀展示分?jǐn)?shù)階導(dǎo)數(shù)數(shù)值計(jì)算格式(32)的精確性，有利與Li等(2016)的計(jì)算結(jié)果進(jìn)行比較，下述2個(gè)算例函數(shù)均選自Li等(2016)。不同的是， Li等(2016)中討論的為奇異核函數(shù)(t-s)-α分?jǐn)?shù)階導(dǎo)數(shù)，采用了更為復(fù)雜的4次或5次多項(xiàng)式以達(dá)到4階或5階精度。

算例1，選取函數(shù)f1(t)=t6，易知分?jǐn)?shù)階導(dǎo)數(shù)的精確解為：

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

表1和表2分別給出了算例1和算例2的計(jì)算誤差和數(shù)值精度。以上2個(gè)算例表明，選取2次或3次Lagrange插值多項(xiàng)式進(jìn)行逼近時(shí)，能達(dá)到3階或4階精度，因此在進(jìn)行Caputo-Fabrizio時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)對流擴(kuò)散方程初邊值問題的數(shù)值模擬時(shí)，采用2次Lagrange插值多項(xiàng)式進(jìn)行模擬。

表1 算例1的數(shù)值結(jié)果

續(xù)表

表2 算例2數(shù)值結(jié)果

3 數(shù)學(xué)模型在地浸采鈾中的應(yīng)用

將分?jǐn)?shù)階導(dǎo)數(shù)的計(jì)算格式嵌入式(10)得出數(shù)學(xué)模型(10)的計(jì)算格式，并應(yīng)用于我國西部某鈾礦采冶區(qū)，結(jié)合鈾礦采冶區(qū)的探測數(shù)據(jù)，進(jìn)一步驗(yàn)證數(shù)學(xué)模型(10)的有效性。

3.1 溶質(zhì)遷移數(shù)學(xué)模型的離散格式

(33)

(34)

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

(35)

式中，

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

(36)

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

(37)

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

(38)

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

(39)

其矩陣形式為：

HU=F

(40)

式中，

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

3.2 數(shù)學(xué)模型砂巖型鈾礦地浸中的應(yīng)用

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

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

選取源項(xiàng)函數(shù)：

(41)

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

通過計(jì)算數(shù)值解與精確解的相對誤差:

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

4 結(jié)論

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