二維單裂隙-孔隙雙重介質(zhì)的核素遷移數(shù)學(xué)模型及參數(shù)反演

張 文， 阮周生， 邱淑芳， 何 杰

(1.東華理工大學(xué)放射性地質(zhì)與勘探技術(shù)國防重點學(xué)科實驗室，江西 撫州 344000;2.東華理工大學(xué)理學(xué)院 江西 撫州 344000)

90Sr是一種高毒性人工放射性核素，物理與生物半衰期均較長，資料顯示核素90Sr在地下水中平均遷移1 km所需時間為104 a。在過去全球性核試驗中，90Sr作為主要核污染物，各國對它在環(huán)境中的污染水平及動向進(jìn)行了觀測與研究(仵彥卿，2007;馬應(yīng)明等，2007;溫瑞媛等，2000)。近年來，隨著核能的廣泛應(yīng)用，各種核設(shè)施排放的廢水中，90Sr仍是一個重要核素。在長期科研生產(chǎn)中積累了大量的含鍶核廢物，特別是在核設(shè)施退役過程中，產(chǎn)生了大量含鍶極低放廢物(VLLW)。核素在孔隙-裂隙中隨地下水遷移問題從模型角度可分為單一介質(zhì)遷移模型與雙重介質(zhì)遷移模型;從計算順序角度又可分為兩類，一類是正問題的研究，即采用何種方法解出滿足模型要求的解。另一類則是反問題的研究，簡單地認(rèn)為是“由果溯因”思想，即根據(jù)測量數(shù)據(jù)來反推模型的特征(如水文地質(zhì)參數(shù)、應(yīng)力場等)。

關(guān)于正問題的研究，國內(nèi)外的許多學(xué)者都取得了有意義的結(jié)論和在一定條件下的解析解(王榕樹等，1994;劉金英等，1996;Ahsanuzzaman et al.，2003;Burnett，2001;Hatano et al.，1998;Gastberger et al.，2000;Borje，1986;Zhang et al.，1995;Pollner et al.，1999)。有學(xué)者較早地對核素在裂隙介質(zhì)中遷移模型進(jìn)行了研究，提出了用有限元與算子分裂迎風(fēng)均衡格式相結(jié)合的方法求解放射性核素遷移方程(謝運棉等，1992)。近年來，楊勇等(2004)、郭敏麗等(2004)、王青海等(2004)先后進(jìn)行了針對90Sr，137Cs，3H和Br－在不同條件下的遷移數(shù)值模擬研究;李尋等(2010)考慮了指數(shù)衰變注入源的基質(zhì)域-裂隙域中溶質(zhì)運移隨地下水遷移的數(shù)學(xué)模型解析解。

由于反問題所求解的滲透系數(shù)、給水度、彈性釋水系數(shù)、擴(kuò)散系數(shù)等工程水文地質(zhì)參數(shù)是工程設(shè)計的基本參數(shù)，一般通過工程地質(zhì)鉆探和抽水試驗等勘探手段或通過經(jīng)驗公式獲得，由于勘探手段代價高且受地域影響大，而經(jīng)驗公式回避了土壤結(jié)構(gòu)的復(fù)雜性，缺乏堅實的物理基礎(chǔ)，因而參數(shù)反演問題備受專家學(xué)者的關(guān)注，張東輝等(2003)、郁伯銘(2003)、劉建國(2001)、劉曉麗等(2004)將分形理論運用于多孔介質(zhì)中，討論多孔介質(zhì)的滲透率、熱導(dǎo)率、熱彌散系數(shù)、水力參數(shù)等輸運性質(zhì)參數(shù);張文等(2010)探討了核素在單裂隙孔隙-裂隙雙重介質(zhì)中的一維遷移模型的解析解以及反演排污點源的反問題。

本文探討了核素90Sr在孔隙-裂隙雙重介質(zhì)中的二維平面遷移模型正、反問題。利用Fourier變換和Saul’yev不對稱差分格式相結(jié)合求得耦合模型的數(shù)值解?；谡龁栴}的計算結(jié)果，采用遺傳算法研究了遷移模型的參數(shù)反演問題，即根據(jù)裂隙中某測量位置的核素濃度反求裂隙遷移模型的裂隙介質(zhì)水動力彌散系數(shù)D1和裂隙平均滲流速度u。數(shù)值結(jié)果表明，參數(shù)反演效果較好。

1 核素遷移的耦合模型

在張文等(2010)基礎(chǔ)上，將遷移模型進(jìn)行推廣，由一維推廣至二維平面遷移模型，同時探討遷移模型的參數(shù)反演問題?？紤]核素在一個二維平面飽和孔隙域中含有一條狹長裂隙的孔隙-裂隙雙重介質(zhì)中的遷移耦合模型，如圖1所示，假定核素在裂隙域中有對流、彌散、衰減和吸附作用;裂隙寬度遠(yuǎn)小于巖塊孔隙厚度，裂隙域為穩(wěn)態(tài)的均勻一維滲流，孔隙域中的水流靜止;含水層的吸附為物理吸附作用，符合等溫線性吸附;核素衰減符合一階衰減動力學(xué)方程;裂隙源頭排污點源以指數(shù)遞減形式注入核素。那么，核素在平面單裂隙多孔介質(zhì)核素遷移過程耦合模型的基本微分方程(仵彥卿，2007;羅嗣海等，2005;李金軒等，2001;中國人民解放軍總裝備部軍事訓(xùn)練教材編輯工作委員會，2002)為裂隙域遷移平衡方程:

為確保結(jié)果準(zhǔn)確性，對模型進(jìn)行帕克檢驗確定其是否具有異方差性。驗證結(jié)果顯示，F(xiàn)統(tǒng)計量的值為1.988844，Prob（F-statistic）=0.108443大于顯著性水平α=0.05，說明模型不存在異方差性，證明模型不存在異方差性，模型設(shè)定完全可靠，回歸分析結(jié)果可以作為分析結(jié)論以及政策建議的依據(jù)。

圖1 單裂隙核素遷移示意圖Fig.1 The schematic diagram of radionuclide migration in a single fracture

孔隙域遷移平衡方程:

其中C1，C2分別表示裂隙域中核素的濃度和孔隙域中核素的濃度;D1表示裂隙介質(zhì)水動力彌散系數(shù);D2，D3分別為孔隙介質(zhì)沿x方向和y方向的水動力彌散系數(shù)(x，y方向如圖1所示);u為平均裂隙的滲流速度;λ為核素衰變常數(shù);b為裂隙半寬;Γ為孔隙介質(zhì)與裂隙域的溶質(zhì)交換量。根據(jù)費克定律知，R1，R2分別為裂隙域的阻滯因子和孔隙域的容量因子，滿足關(guān)系式，Q為污染物總質(zhì)量，φ為巖石的有效孔隙度，ρb為巖石密度，Kf為吸附分布系數(shù)，Km為孔隙域中溶質(zhì)平均分配系數(shù)，δ(x)為狄拉克函數(shù)。該模型對均質(zhì)各向同性介質(zhì)適用。

2 正問題的求解

本文所考慮的正問題是，已知裂隙源頭排污點的核素濃度變化狀態(tài)，求出核素濃度在裂隙和孔隙中的變化規(guī)律。為了能夠求解耦合方程組，對孔隙域遷移平衡方程兩邊關(guān)于x和y進(jìn)行二維Fourier變換，求出耦合項，再對裂隙域遷移平衡方程使用Saul’yev不對稱差分格式的有限差分法求解得出數(shù)值解。

引進(jìn)Fourier變換的相關(guān)符號，記

根據(jù)Fourier變換的性質(zhì)容易得到

由已知條件

得到

聯(lián)立式(6)和式(7)，得出

進(jìn)一步使用Fourier逆變換式(4)以及位移性質(zhì)，便得到核素在孔隙域中變化規(guī)律的解析解為

下面求解核素在裂隙域中變化規(guī)律的數(shù)值解。由式(9)易知

則耦合項為

將耦合項代入裂隙域遷移平衡方程中得到定解方程

引進(jìn)Saul’yev不對稱差分格式，該差分格式是無條件穩(wěn)定的(張文生，2006)，對于函數(shù)u(x，t)來說，其差分格式為

將式(11)離散后便容易得出核素在裂隙域中的數(shù)值解。

3 反問題的求解

考慮耦合模型中2個系數(shù)的反演問題，即在裂隙介質(zhì)水動力彌散系數(shù)D1和裂隙平均滲流速度u未知情況下，通過附加某測量位置x0、不同時刻測得裂隙域中的核素濃度C1(x0，t)，T1≤t≤T2，來反求裂隙介質(zhì)水動力彌散系數(shù)D1和裂隙平均滲流速度u。

記 f(t;D1，u)=C1(x0，t;D1，u)為已知 D1和 u經(jīng)正問題求解出的濃度，設(shè)測量位置x0、不同時刻測得裂隙域中的核素濃度為h(t)，引進(jìn)泛函

則反問題可歸結(jié)為求參數(shù)D1和 u，使得泛函J(f)達(dá)到極小，即使得

成立的D1和u就是系數(shù)反演反問題的最優(yōu)解。

為探討解決此類問題的有效方法，引進(jìn)遺傳算法。眾所周知，遺傳算法根據(jù)達(dá)爾文進(jìn)化論中的“生存競爭”和“優(yōu)勝劣汰”原則，不受求導(dǎo)和函數(shù)連續(xù)性的限定，從某一初始群體出發(fā)，借助復(fù)制、交叉、變異等操作，不斷迭代計算，經(jīng)過若干代的演化后，群體中的最優(yōu)值逐步逼近最優(yōu)解，直至最后達(dá)到全局最優(yōu)(王彥飛，2007;王澤文等，2008;顧正華，2006;閔濤等，2004)。

遺傳算法描述如下:

STEP 1 產(chǎn)生初始群體。在給定的參數(shù)范圍內(nèi)均勻隨機(jī)生成實數(shù)編碼建立初始群體，種群數(shù)量為N.個體中每個十進(jìn)制數(shù)計算公式為

其中r為[0，1]服從均勻分布的隨機(jī)數(shù)，X為個體目標(biāo)變量，Xmin，Xmax分別為個體目標(biāo)變量可取的最小、最大值。

STEP 2 交叉。將選出的采用綜合交叉法將N個個體兩兩雜交，生成N個新的個體，交叉算法如下:

其中，X和Y為父代個體，X′和Y′為子代個體，r為[0，1]服從均勻分布的隨機(jī)數(shù)。

STEP 3 計算適應(yīng)值。取式(14)作為遺傳算法的適應(yīng)函數(shù)。

STEP 4 選擇。對STEP 3中得到的由父代和子代共2N個個體的適應(yīng)值進(jìn)行比較，從中挑選N個更小的個體作為新一代個體。

STEP 5 變異。對STEP 4中選出的N個個體按照式(16)進(jìn)行變異。

STEP 6、判斷新一代群體是否滿足結(jié)束條件，若滿足則停機(jī)，否則轉(zhuǎn)至STEP 2繼續(xù)計算。

4 數(shù)值算例

圖2 平面孔隙中所含核素遷移規(guī)律Fig.2 Migration of radionuclides contained in the plane pore

給定計算參數(shù)(李錄等，1998)。核素90Sr衰變系數(shù) λ =1.54 × 10－4d－1，裂隙平均滲流速度 u=0.1 m/d，裂隙域的阻滯因子和孔隙域的容量因子分別為R1=1.0，R2=0.15，裂隙介質(zhì)水動力彌散系數(shù)D1=0.351 m2/d，孔隙介質(zhì)水動力彌散系數(shù)D2=1.38 × 10－5m2/d，D3=6.9 × 10－5m2/d，裂隙半寬b=2.5 × 10－5m.

實際問題中，測量數(shù)據(jù)往往帶有測量隨機(jī)誤差。因此，通過 h(x)=C1(x，T)+ ε·C1(x，T)·R獲得測量數(shù)據(jù)，其中C1(x，T)為正問題的解，ε為噪聲水平，R為介于區(qū)間[－1，1]且服從均勻分布的隨機(jī)數(shù)。遺傳算法的參數(shù)設(shè)置如下:群體大小(30)，交叉概率(0.8)，變異概率(0.2)，參數(shù)范圍(0.01，0.8)，終止代數(shù)(500)。

正問題的計算結(jié)果見圖2～5，圖2表示平面孔隙中所含核素當(dāng)時刻T=20041 d時在平面上的分布圖，圖3表示平面孔隙中所含核素時刻T=20041 d，x=100 m時的變化規(guī)律，圖4表示單裂隙中所含核素各個位置隨時間變化的遷移規(guī)律，圖5表示單裂隙中所含核素在測量位置x=506 m時隨時間變化的遷移規(guī)律。

表1為在不同測量誤差下的參數(shù)反演結(jié)果。圖6為模擬裂隙x0=60 m，t∈[1 000，10 000]不同噪音水平ε時的測量數(shù)據(jù)。從表1可以看出，本文算法對于單裂隙核素遷移的耦合模型參數(shù)反演是有效的。

圖3 平面孔隙中所含核素遷移規(guī)律(x=100)Fig.3 Migration of radionuclides contained in the plane pore(x=100)

5 結(jié)語

圖4 單裂隙中所含核素遷移規(guī)律Fig.4 Migration of nuclides contained in a single fracture

圖6 不同噪音水平下的觀測數(shù)據(jù)(x0=60 m)Fig.6 The observed data under different noise levels(x0=60 m)

針對二維單裂隙-孔隙雙重介質(zhì)系統(tǒng)中的核素遷移耦合問題，利用Fourier變換方法與Saul’yev不對稱差分格式得到正問題的解;另一方面，首先將參數(shù)反演反問題歸結(jié)為一等價的泛函極小化問題，然后利用遺傳算法經(jīng)過交叉、選擇和變異等步驟得出反問題的解。從算例的計算結(jié)果可以看出，正問題的解析解能較好地刻畫核素的遷移規(guī)律;遺傳算法也能有效地實現(xiàn)單裂隙-孔隙雙重介質(zhì)耦合模型地質(zhì)參數(shù)的反演。對于有效地反演多裂隙以及三維裂隙-孔隙雙重介質(zhì)系統(tǒng)的遷移問題，需要展開更深入的研究。

圖5 單裂隙中所含核素遷移規(guī)律(x0=506 m)Fig.5 Migration of nuclides contained in a single fracture(x0=506 m)

表1 x0=60 m時參數(shù)反演結(jié)果Tab.1 The inversion results of the two parameters under different noise levels(x0=60 m)

顧正華.2006.廣義遺傳算法及其在水流參數(shù)反演中的應(yīng)用[J].系統(tǒng)工程理論與實踐，26(2):133-137.

郭敏麗，王金生，李書紳.2004.飽和黃土中3H和Br－遷移的延遲特征[J].環(huán)境科學(xué)學(xué)報，24(5):878-882.

李金軒，李尋.2001.基于雙重介質(zhì)理論的單裂隙核素遷移模型[J].勘察科學(xué)技術(shù)，2:7-11.

李錄，李春江，張吳平.1998.單裂隙花崗巖中核素遷移的連續(xù)數(shù)學(xué)模型研究[J].山西大學(xué)學(xué)報:自然科學(xué)版，21(1):1-9.

李尋，張土喬，李金軒.2010.單裂隙介質(zhì)中核素遷移模型的解析解及其應(yīng)用[J].哈爾濱工業(yè)大學(xué)學(xué)報，42(12):1990-1994.

梁冰，劉磊.薛強(qiáng)，等.2007.核素滲漏對地下水污染的數(shù)值仿真研究[J].系統(tǒng)仿真學(xué)報，19(2):261-263.

劉建國，聶永豐.2001.非飽和土壤水力參數(shù)預(yù)測的分形模型[J].水科學(xué)進(jìn)展，12(1):99-106.

劉金英，王新民，彭澤洲.1996.雙重介質(zhì)中核素遷移模型的解析解及其應(yīng)用[J].世界地質(zhì)，15(1):68-73.

劉曉麗，梁冰，薛強(qiáng)，等.2003.多孔介質(zhì)滲透率的分形描述[J].水科學(xué)進(jìn)展，14(6):769-773.

羅嗣海，錢七虎，李金軒，等.2005.高放廢物深地質(zhì)處置中的多場耦合與核素遷移[J].巖土力學(xué)，26(Z1):264-270.

馬應(yīng)明，金玉仁.2007.钚在地下水中的存在形態(tài)分析[J].東華理工學(xué)院學(xué)報，30(3):236-240.

閔濤，周孝德，張世梅，等.2004.對流－擴(kuò)散方程源項識別反問題的遺產(chǎn)算法[J].水動力學(xué)研究與進(jìn)展(A輯)，19(4):520-524.

王青海，王蘭生，李曉紅.2004.基巖裂隙水中90Sr遷移的數(shù)值模擬[J].吉林大學(xué)學(xué)報，34(3):405-409.

王榕樹，馮為.1994.放射性核素在地質(zhì)介質(zhì)中的遷移研究[J].核化學(xué)與放射化學(xué)，16(2):117-121.

王彥飛.2007.反演問題的計算方法及其應(yīng)用[M].北京:高等教育出版社.

王澤文，邱淑芳.2008.一類流域點污染源識別的穩(wěn)定性與數(shù)值模擬[J].水動力學(xué)研究與進(jìn)展(A 輯)，23(4):364-371.

仵彥卿.2007.多孔介質(zhì)污染物遷移動力學(xué)[M].上海:上海交通大學(xué)出版社.

謝運棉，楊天行，王運國等.1992.放射性核素在裂隙巖石中的遷移模型研究[J].巖石學(xué)與工程學(xué)報，11(1):1-24.

楊巍，楊亞新，曹龍生，等.2011.某鈾尾礦庫中放射性核素對環(huán)境的影響[J].東華理工大學(xué)學(xué)報:自然科學(xué)版，34(2):231-235.

楊勇，苑國琪，張東.2004.90Sr，137Cs在某種包氣帶土壤中的遷移研究[J].四川環(huán)境，23(3):85-89.

郁伯銘.2003.多孔介質(zhì)輸運性質(zhì)的分形分析研究進(jìn)展[J].力學(xué)進(jìn)展，33(3):333-346.

張東輝，施明恒.2003.多孔介質(zhì)擴(kuò)散、導(dǎo)熱、滲流分形模型的研究[D].南京:東南大學(xué).

張文，王澤文，樂勵華.2010.雙重介質(zhì)中的一類核素遷移數(shù)學(xué)模型及其反演[J].巖土力學(xué)，31(2):553-558.

張文生.2006.科學(xué)計算中的偏微分方程有限差分法[M].北京:高等教育出版社.

中國人民解放軍總裝備部軍事訓(xùn)練教材編輯工作委員會.2002.地下核爆炸現(xiàn)象學(xué)概論(下冊)[M].北京:國防工業(yè)出版社.

Ahsanuzzaman A，Kolar R，Zaman M.2003.Limiting Source Dimensions of Three-Dimensional Analytical Point Source Model for Solute Transport[C].Hydrology Days Proceedings:1-15.

Pollner L，Zagyvai P.1999.Subsurface transport modeling of longlived nuclides[J].Waste Management，19:45-54.

Zhang L X，Zhang M S.Tu G R.1995.A field study of tritium migration in groundwater[J].The Science of the Total Environment.173/174:47-51.

Gastberger M，Steinhausler F，Gerzabek M H，et al.2000.Soil-toplant transfer of fallout caesium and strontium in Austrian lowland and Alpinepastures[J].J.Environ Radioactivity，49:217-233.

Borje T.1986.Migration of the fission products strontium.tech-netium.lodine and cesium in clay[J].Radiochim Acta，39:97-104.

Burnett W C.2001.Elzerman A W.Nuclide migration and theenvironmental radiochemistry of florida phosphogyp sum[J].J Environ Radioactivity，54:27-51.

Hatano Y，Hatano N，Amano H.et al.1998.Aerosol migration near chernobyl:longterm data and modeling[J].A tmospheric Environment，32(14-15):2587-2594.