不同溶蝕速率對灰?guī)r單裂隙滲流-溶解模型的影響

2019-11-15 02:11:52

人民長江 2019年10期

(1.長江水利委員會水文局長江下游水文水資源勘測局，江蘇南京 210011； 2. 南京水利科學研究院水文水資源與水利工程科學國家重點實驗室，江蘇南京 210029)

我國巖溶分布面積廣，引起的資源環(huán)境及工程問題也較為突出，如巖溶區(qū)隧洞涌突水、水庫滲漏、邊坡穩(wěn)定等問題，而此類問題多與巖溶演化發(fā)育有關[1-4]。地下水在運移過程中與可溶性碳酸鹽類發(fā)生反應，引起礦物溶解，在以裂隙為主的巖溶系統(tǒng)中使得裂隙開度發(fā)生擴張，并導致滲流場、化學場等發(fā)生不可逆的改變，此過程可通過滲流-溶解耦合模型加以反映刻畫[5]。模型主要由水流模塊、溶質(zhì)運移模塊、裂隙開度時變模塊等組成，并通過各模塊變量間的函數(shù)關系實現(xiàn)模型耦合[6-7]。其中，礦物溶蝕速率作為模型的重要組成，在以上3個模塊中均有所體現(xiàn)。已有研究表明，在不考慮CO2對反應速率影響的封閉環(huán)境中，溶蝕速率主要受表面反應和擴散遷移雙重因素的聯(lián)合控制，該兩項進程逐次進行，溶蝕速率值為兩者間相對較小者[8]。但在實際模擬應用中，由于該速率方程表達較為復雜，常僅考慮受單因素控制，即只受表面反應或擴散遷移其中某一因素控制，從而使得不同研究中溶蝕速率方程有所差異[9-11]。可以看出，以往研究多側(cè)重于研究裂隙系統(tǒng)的發(fā)育演變，以及由此帶來的滲流動態(tài)、溶液組分等的改變，而未涉及在模擬過程中，當溶蝕速率表達不同時，對模擬結果存在何等影響。

針對目前研究中存在的不足，本文旨在量化解析不同溶蝕速率方程對模型結果的影響。首先根據(jù)溶蝕方程機理的不同，分為單因素控制(表面反應或擴散遷移)以及兩者聯(lián)合控制。然后假設裂隙形態(tài)皆為平行板裂隙，分別模擬獲得不同溶蝕速率下裂隙開度、滲流性態(tài)以及溶液組分等的空間分布和時間演變規(guī)律，并進行比較獲得其異同。研究成果可為今后實際應用中根據(jù)需要選擇更為合適的模型提供指導。

1 溶蝕機理與速率模型

灰?guī)r單裂隙滲流-溶解耦合模型主要由水流模塊、溶質(zhì)運移模塊和裂隙開度時變模塊3部分組成，具體可見文獻[7，10]，在此不再贅述。由于灰?guī)r溶蝕速率是滲流-溶解耦合模型的核心內(nèi)容之一，在模型3個模塊中皆有涉及。根據(jù)其受控機理的不同，一般包括以下3種：① 僅受表面反應單一因素控制；② 僅受擴散遷移單一因素控制；③ 受表面反應和擴散遷移雙重因素聯(lián)合控制。以上各機理及速率模型如下所述。

(1) 表面反應控制，即僅受可溶性礦物的溶解動力學控制，其速率大小與溶液中對應礦物的相對飽和度有關。對于碳酸鹽類，在接近平衡態(tài)和遠離平衡態(tài)時，其速率表達略有不同。當遠離平衡態(tài)時反應相對較快，而接近飽和平衡態(tài)時，其速率陡降，該過程可表示為

(1)

式中，k1、k2為動力學反應常數(shù)；n1、n2為指數(shù)，通常n1為等于1或為接近于1的數(shù)，而n2在3～6之間；Cs為拐點濃度值，一般為0.7Ceq～0.9Ceq之間。

(2) 擴散遷移控制，即對應離子組分從固相表面進入溶液中的遷移速率，一般認為該速率不僅與溶液中離子濃度有關，還與該位置處裂隙開度b有關，表達如下：

(2)

式中，Dm為離子在水溶液中的擴散系數(shù)；ShD為Sherwood數(shù)，在層流情形下為8.24[12]；Ceq為離子飽和濃度。

需要指出的是，國內(nèi)在模擬滲流-溶解耦合問題時主要基于“薄膜假設”，認為緊貼裂隙壁有一層極薄的“薄膜”，內(nèi)部溶液濃度飽和，薄膜內(nèi)組分可進入擴散區(qū)，導致薄膜內(nèi)溶液濃度降低；此時，固壁表面發(fā)生溶解作用，離子進入“薄膜”補充其內(nèi)物質(zhì)損失，使之重新達到飽和；同時，由于固壁表面物質(zhì)的溶失，造成裂隙開度的增大[13]。該假設本質(zhì)上也為擴散遷移控制，但在該速率模型中不考慮裂隙開度變化對溶蝕速率的影響。

(3) 表面反應和擴散遷移兩因素聯(lián)合控制，即認為在灰?guī)r溶解過程中，表面反應與擴散過程逐次依序進行，因此總的溶解速率取決于各時刻以上兩過程中的最慢者，即以上式(1)與式(2)計算結果間相對較小者作為該時刻下灰?guī)r的溶解速率。

2 不同溶蝕速率對模型影響分析

2.1 模型簡介

巖溶地區(qū)裂隙相互交錯形成裂隙網(wǎng)絡，其真實形態(tài)極為復雜，因此常以對單裂隙滲流特性的研究作為認識巖體裂隙整體滲流特征的基礎。將裂隙壁簡化為光滑平板，同時隙寬相對較小時，可采用立方定律研究模型隙寬與單寬流量之間的關系，在此基礎上加以修正可推廣至粗糙裂隙及裂隙網(wǎng)絡。因此，平行板裂隙下滲流特性的研究作為真實裂隙網(wǎng)絡的基礎具有較為重要的意義[14-15]。

考慮研究區(qū)為一長為L、寬為W的單裂隙，其初始裂隙隙寬為b0，裂隙兩端水頭分別為Hin和Hout，流入溶液的濃度為Cin，裂隙內(nèi)初始水頭和初始濃度分別為H0和C0，模型示意圖如圖1所示，具體取值見表1。

圖1 模型示意Fig.1 Sketch of model

L/mW/mαLαTM/(g·mol-1)ρ/(kg·m-3)Cin/(mol·m-3)C0/(mol·m-3)b0/mmH0/mHin/mHout/m10.50.050.0210026600.0500.200.10

2.2 溶蝕速率參數(shù)及比較

表面反應和擴散遷移控制中涉及到的主要參數(shù)如表2所示。為研究不同溶蝕速率方程對模型結果的影響，分別建立以下模型，其溶蝕速率對應控制因素為：模型1為僅表面反應單一因素控制，模型2為僅擴散遷移單一因素控制，模型3為表面反應及擴散遷移兩者聯(lián)合控制。

表2 反應速率方程中相關參數(shù)Tab.2 Parameters in reaction rate equation

注：表中參數(shù)引自文獻[8]。

在有限元軟件COMSOL Multiphysics軟件中建立尺寸為L×W(長×寬)的二維研究區(qū)域，并設定相應的模型方程和定解條件。采用瞬態(tài)分離式求解器( time dependent segregated) 對以上3個模型分別進行求解。在其它條件相同情形下，模擬時段(2 000 d)內(nèi)選取500,1 000,1 500 d和2 000 d，在上述時刻時各模型裂隙開度變化如圖2所示。由圖2可以看出：

(1) 3個模型在上游進水口位置溶蝕皆較為嚴重，而向下游溶蝕逐漸減緩，從溶蝕鋒面來看，在前期(500 d和1 000 d)皆較為平直，而在中后期(1 500 d和2 000 d)出現(xiàn)溶蝕差異不均衡現(xiàn)象。

(2) 3個模型的不同主要體現(xiàn)在溶蝕“程度”和“深度”兩方面。從溶蝕“程度”來看，在研究區(qū)上游側(cè)，尤其是靠近進水口部位，模型1結果較模型2和模型3明顯偏大，而在研究區(qū)的中后部，模型1和模型2未發(fā)生明顯溶蝕，其開度b較模型3明顯偏小；從溶蝕“深度”看，模型3>模型2>模型1，其中模型3于2 000 d時在研究區(qū)靠近軸部的部位形成了一條較寬的集中滲流通道，且基本貫穿整條裂隙，模型2則形成了局部的滲流通道，達到了裂隙中部，但未貫穿，而模型1的溶蝕僅發(fā)生在靠近進水口的較短距離內(nèi)(<0.2 m)。

圖2 不同時刻下裂隙開度b的分布(由上至下依次為500,1 000 ,1 500 d和2 000 d)Fig.2 Distribution of fracture aperture at different time(from top to bottom are in 500,1 000,1 500 d and 2 000 d)

為更直觀地比較3個模型各時刻下開度在沿裂隙延伸方向上的變化，以y=0.25 m作為典型斷面，以上時刻各模型開度b在該典型斷面上的分布見圖3。

由圖3可以看出：相同時刻下，在裂隙前端靠近進水口部位(<0.1 m)，模型2和模型3的開度分布較為近似，且都明顯小于相同位置處模型1的。這表明模擬時段較長時，該部位溶解主要受擴散遷移控制，而采用表面反應控制的動力學表達，會導致該部位灰?guī)r溶蝕的過量估計。在研究區(qū)中后段(>0.4 m)，尤其是模擬時刻的后期(如2 000 d)，模型1和模型2的結果明顯小于模型3的，表明單一因素控制的模型在研究區(qū)范圍相對較大時，其位于研究區(qū)中后部的結果誤差相對較大。

可通過計算獲得模擬時段內(nèi)研究區(qū)平均隙寬隨時間的變化。平均隙寬一方面可反映裂隙形態(tài)的變化，另一方面也可反映灰?guī)r裂隙整個溶失量的變化。同時，為更直觀地比較各模型間的差異，以模型3的結果作為基準，計算獲得模型1、模型2與模型3間的比例。各模型平均隙寬及對應比例隨時間變化見圖4。

由圖4可以看出，3個模型計算獲得的實際平均隙寬在整個模擬時段內(nèi)皆呈增大的趨勢，但彼此間的大小關系在不同時刻有所不同。在該算例中，各模型間的差異性變化可分為以下4個階段：① 300 d之前，各模型的平均隙寬表現(xiàn)為模型2>模型1>模型3；② 300～510 d之間，表現(xiàn)為模型1>模型2>模型3；③ 在510～1 370 d之間，則為模型1>模型3>模型2；④ 在1 370 d之后，模型3>模型1>模型2，且這段期間模型3與其他兩個模型間的差異越來越大，在模擬結束時模型1和模型2僅為模型3值的0.739和0.593。以上4個階段可簡單概括為：前期單一因素控制下模型平均隙寬變化大于兩因素聯(lián)合控制的；而后模型3(兩因素聯(lián)合控制)的隙寬變化依次超過模型2(擴散遷移控制)和模型1(表面反應控制)，且與后兩者間的差異越來越顯著。

由以上可知，不同溶蝕速率方程對裂隙開度的變化有所影響，因此通過整個裂隙的流量也會有所差異。同樣以模型3的結果作為基準，分別計算獲得模型1、模型2與模型3間的比例，則出口位置處各模型通過裂隙的流量Q及對應比例隨時間的變化如圖5所示。

由圖5可以看出：① 3種模型中通過裂隙的流量，總體上隨時間的推移皆呈逐漸增長的趨勢，而相同時刻下，三者比較，總體上呈“模型3>模型1>模型2”。該結果表明，由單一因素控制溶解速率的模型計算獲得的流量結果普遍低于實際通過裂隙的流量，其中尤以模型2的差別相對較大。在290d時，模型2與模型3的差別即超過了10%，而模型1在1080d后才達到該數(shù)值。②模型3中流量增長出現(xiàn)“臨界點”，在該時刻前增長速率較為平緩，而該時刻后，其速率陡增，并最終達到初始時的96.57倍，此與模型3中裂隙形成貫穿型滲流通道有關；而模型1和模型2在整個模擬時段內(nèi)流量皆呈緩慢增長的趨勢，無明顯的流量激增時刻，而終了時刻時的流量也僅為初始時的3.22倍和1.83倍。以上結果表明，在對通過裂隙的滲流量進行預測時，采用單因素控制的模型會低估通過裂隙的滲流量，尤以擴散遷移控制的模型表現(xiàn)明顯。

圖3 不同時刻下y=0.25 m處裂隙開度變化Fig.3 Fracture aperture along y=0.25 m

圖6 不同時刻下研究區(qū)內(nèi)Ca2+濃度的分布(由上至下依次為500，1 000，1 500 d和2 000 d)Fig.6 Distribution of Ca2+ concentration in 500,1 000,1 500 d and 2 000 d

圖4 不同模型裂隙平均隙寬及相應比例Fig.4 Average aperture width and corresponding proportion of different models

圖5 不同模型通過裂隙流量及比例過程線Fig.5 Flow through aperture of different models and corresponding proportion curve

除裂隙開度、滲流場外，模型間化學場分布亦有所不同，其中Ca2+分布如圖6所示。由圖6可以看出各模型計算獲得的Ca2+分布具有以下異同點。

(1)對于每個模型而言，Ca2+濃度沿裂隙延伸方向逐漸增大，在進水口位置處等于流入溶液的Ca2+濃度，而在一段距離后接近飽和；且與裂隙開度分布相似，鋒面在前期較為平直，后在中后期出現(xiàn)溶蝕差異。

圖7 不同時刻y=0.25 m處Ca2+濃度分布Fig.7 Concentration of Ca2+ along y=0.25 m

(2)模型3由于在模擬時段后期在近裂隙軸部部位形成較為集中的滲流通道，且貫穿整個裂隙，因此在該段時期裂隙溶液中Ca2+濃度相對較小，而相同時段內(nèi)模型1和模型2該部位的Ca2+濃度仍相對較大。

(3)在模擬時段前期(<1 500 d)，盡管3個模型在裂隙中后部的溶液Ca2+濃度皆接近飽和，但其值大小仍有區(qū)別，其中模型2的濃度值已基本等于灰?guī)r溶解的飽和濃度值，即0.1768 mol/m3；而模型1和模型3的最大值仍小于0.16 mol/m3，這與表面反應控制在溶液達到拐點濃度后，其反應速率大大減慢有關。

為更清晰地比較各模型中Ca2+濃度沿裂隙延伸方向上的變化，以y=0.25 m作為典型斷面，則各模型在該斷面上Ca2+濃度分布如圖7所示。由圖7可以看出：

(1)無論是單一因素控制的模型1和模型2，還是雙重因素控制的模型3，Ca2+濃度在沿裂隙延伸方向上皆呈逐漸增大的趨勢，且在裂隙前部增長較快，而在裂隙中后部增長緩慢。

(2)模型1與模型3在裂隙中后部Ca2+的濃度值呈緩慢增長的趨勢，但仍未達到飽和狀態(tài)，而模型2則達到飽和平衡態(tài)。

3 結論

在建立灰?guī)r單裂隙滲流-溶解耦合模型的基礎上，研究了不同溶蝕速率對模型結果的影響，并以一平行板裂隙中滲流溶解為例，得出如下結論。

(1) 對于裂隙開度而言，從溶蝕深度來看，模擬前期3種模型基本一致，而中后期雙重因素下模型明顯更深，從而形成貫穿裂隙的滲流通道；就溶蝕程度來看，3種模型下進水口位置溶蝕更為嚴重，其中尤以表面反應單因素控制的更為明顯，其在上游位置的溶蝕明顯大于另外兩種模型。

(2) 就總的溶失量而言，前期單一因素控制下模型平均隙寬變化大于兩因素聯(lián)合控制的，但相差不大；而后兩因素聯(lián)合控制依次超過擴散遷移控制和表面反應控制，且與后兩者間的差異越來越顯著。

(3) 從通過裂隙的滲流量來看，采用單因素控制的模型會低估通過裂隙的滲流量，尤以擴散遷移控制的模型表現(xiàn)明顯。同時，溶蝕使得裂隙中可形成貫穿型集中滲流通道，從而引起流量陡增，單因素控制的模型預測的該臨界點遠遠晚于實際。

(4) 從化學場分布來看，模擬前期3種模型差別并不顯著，但隨著兩因素聯(lián)合控制的模型形成貫穿型滲流通道，其裂隙內(nèi)Ca2+濃度分布明顯小于另外另種模型。此外，對于擴散遷移控制的模型，其在裂隙下游側(cè)溶液中Ca2+濃度達到飽和，其值明顯大于另兩種模型。

(5) 當采用單因素控制對溶解速率進行簡化時，會對模型結果造成一定影響，但可根據(jù)應用中關注的內(nèi)容(如裂隙開度、溶蝕量、滲流量、Ca2+濃度等)以及模擬時段的長度，從表面反應和擴散遷移中選擇合適的模型，亦可滿足一定的精度要求。