砂巖型鈾礦溶質(zhì)遷移的數(shù)學(xué)模型及數(shù)值算法

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

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

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

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

(1)

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

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

(5)

(6)

(7)

為使得表達(dá)形式更為簡(jiǎn)潔，考慮空間方向?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ù)對(duì)流擴(kuò)散方程初邊值問(wèn)題的數(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)用于我國(guó)西部某鈾礦采冶區(qū)，結(jié)合鈾礦采冶區(qū)的探測(cè)數(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)

式中，

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

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

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

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

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

(41)

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

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

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

4 結(jié)論

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