本發(fā)明涉及一種壩基帷幕防滲性能衰減的數(shù)值模擬方法，屬于大壩安全監(jiān)測技術(shù)領(lǐng)域。

背景技術(shù)：

帷幕作為壩基防滲體系的重要組成部分，對大壩安全運行及防滲起著極其重要的作用。水利工程中有眾多事故是由滲流所引發(fā)的，統(tǒng)計結(jié)果表明我國1981年241座大型水庫的100宗事故中，由滲流問題而導(dǎo)致水庫失事的占到30％～40％左右。因此，壩基帷幕的防滲問題得到越來越廣泛的重視。

目前對帷幕防滲的研究主要集中在以下兩個方面：(1)帷幕不同性態(tài)參數(shù)條件下對壩基滲流的影響：如通過數(shù)值模擬手段研究壩基帷幕不同深度、厚度、位置等因素對防滲效果的影響，以及通過建立帷幕隨機缺損的數(shù)學模型模擬防滲體缺損大小、分布隨機性對滲流場的影響；(2)對帷幕防滲效果及其耐久性進行評價：如通過單一指標，如水質(zhì)等微觀效應(yīng)量，基于水化學圖示法、飽和指數(shù)模型等方法對不同時期的帷幕防滲性能進行評價，以及通過建立幕后地下水位、排水量、揚壓力等宏觀指標和壩基地下水溫度、pH值及其他水化學特征等微觀指標相結(jié)合的多指標評價體系進行評價。但申請人檢索發(fā)現(xiàn)，目前對水庫運行過程中帷幕防滲性能衰減的研究還較為少見，尤其是考慮帷幕中水泥結(jié)石等組分的溶失對帷幕細觀結(jié)構(gòu)的改變等還未涉及。

技術(shù)實現(xiàn)要素：

本發(fā)明要解決技術(shù)問題是：提供一種模擬大壩基礎(chǔ)帷幕防滲性能衰減的方法，該方法可以為大壩的安全運行和維護提供依據(jù)從而減少因滲流引起的水壩失事事故。

為了解決上述技術(shù)問題，本發(fā)明提出的技術(shù)方案是：一種壩基帷幕防滲性能衰減的數(shù)值模擬方法，包括以下步驟：

1)根據(jù)待模擬的壩基的實際物理參數(shù)設(shè)置邊界條件，并得到模擬前所述壩基的滲流速度U₀、孔隙度θ₀、目標組分濃度C₀、滲透系數(shù)K₀、彌散系數(shù)D₀以及溶解速率R₀；所述壩基包括基巖和帷幕兩部分,所述目標組分濃度包括鈣離子濃度和氫氧根離子濃度。如何設(shè)置邊界條件為現(xiàn)有技術(shù)，本領(lǐng)域技術(shù)人員在看到壩基的實際物理參數(shù)后可以很容易地得到邊界條件，不再贅述。

2)根據(jù)步驟1)設(shè)置的邊界條件，建立滲流方程、溶質(zhì)遷移方程和固相介質(zhì)時變方程；

所述滲流方程為：

$\frac{\∂}{\∂x}\left(K\frac{\∂H}{\∂x}\right)+\frac{\∂}{\∂y}\left(K\frac{\∂H}{\∂y}\right)=\mathrm{\μ}\frac{\∂H}{\∂t}$

式中，H為水頭，K為滲透系數(shù)，t為時間參數(shù)，μ為給水度，x、y分別為壩基的橫坐標、縱坐標；

所述溶質(zhì)遷移方程為：

$\frac{\∂\left(\mathrm{\θ}C\right)}{\∂t}=\frac{\∂}{\∂x}\[D_L\frac{\∂\left(\mathrm{\θ}C\right)}{\∂x}\]+\frac{\∂}{\∂y}\[D_T\frac{\∂\left(\mathrm{\θ}C\right)}{\∂y}\]-\frac{\∂\left(u_x\mathrm{\θ}C\right)}{\∂x}-\frac{\∂\left(u_y\mathrm{\θ}C\right)}{\∂y}+\mathrm{\θ}\·R$

式中，R為可溶組分的溶解速率，C為目標組分濃度；滲流速度U分解為水平方向的滲流速度u_x和垂直方向的滲流速度u_y；θ為孔隙度；彌散系數(shù)D分解為縱向彌散系數(shù)D_T和橫向彌散系數(shù)D_L；x、y分別為壩基的橫坐標、縱坐標；t為時間參數(shù)；

所述固相介質(zhì)時變方程為：

$\frac{\∂\mathrm{\θ}}{\∂t}=R\·M\·\mathrm{\θ}$

式中，M為Ca(OH)₂的摩爾體積，R為可溶組分的溶解速率，θ為孔隙度，t為時間參數(shù)。

3)將待模擬時段劃分為n個子時段，其中第i個子時段的初始時刻為終了時刻為1≤i≤n，n≥2；

根據(jù)步驟1)得到的所述壩基帷幕的滲流速度U₀、孔隙度θ₀、目標組分濃度C₀、滲透系數(shù)K₀、彌散系數(shù)D₀以及溶解速率R₀，對步驟2)建立的滲流方程、溶質(zhì)遷移方程和固相介質(zhì)時變方程進行求解，依次計算得到第j個子時段的終了時刻的滲流速度、孔隙度以及目標組分濃度，1≤j≤n；同時利用滲透系數(shù)K與孔隙度θ之間的函數(shù)、彌散系數(shù)D與滲流速度U、孔隙度θ之間的函數(shù)以及可溶組分的溶解速率R與目標組分濃度C之間的函數(shù)，得到第j個子時段的終了時刻的滲透系數(shù)、彌散系數(shù)以及溶解速率；并將第j個子時段的終了時刻的滲流速度、孔隙度、目標組分濃度、滲透系數(shù)、彌散系數(shù)以及溶解速率作為第j+1個子時段的初始時刻的相應(yīng)值；

其中，所述滲透系數(shù)K與孔隙度θ之間的函數(shù)為：

K＝K₀·(θ/θ₀)³[(1-θ₀)/(1-θ)]²

式中，K₀和θ₀分別為模擬前所述壩基帷幕的滲透系數(shù)和孔隙度；

所述彌散系數(shù)D與滲流速度U、孔隙度θ之間的函數(shù)為：

D＝α·|U|+D_m·θ

式中，α為彌散度，D_m為鈣離子在水中的擴散系數(shù)；

所述可溶組分的溶解速率R與目標組分濃度C之間的函數(shù)為：

$R=A\·{(1-\frac{C_{{\mathrm{Ca}}^{2+}}\·C_{{\mathrm{OH}}^{-}}^2}{K_{sp}})}^n$

式中，A為反應(yīng)動力學系數(shù)，n為動力學指數(shù)，K_sp為Ca(OH)₂的溶度積常數(shù)，為鈣離子濃度，為氫氧根離子濃度；A，n，K_sp均為預(yù)設(shè)系數(shù)。

模擬前所述壩基的孔隙度θ₀可由壩基的幾何參數(shù)直接得到。步驟1)中以上游庫水位作為初始水頭，在所述邊界條件的約束下對壩基進行滲流和溶質(zhì)運移的穩(wěn)定態(tài)模擬，獲得滲流穩(wěn)定態(tài)時的滲流場和化學場分布情況，從而得到模擬前所述壩基的滲流速度U₀、目標組分濃度C₀、滲透系數(shù)K₀、彌散系數(shù)D₀以及溶解速率R₀。在所述壩基的孔隙度θ₀已知的前提下，對壩基進行滲流和溶質(zhì)運移穩(wěn)定態(tài)模擬為現(xiàn)有技術(shù)，不再贅述。

本發(fā)明的原理及帶來的有益效果如下：

大壩在長期運行條件下，帷幕受到壩基水的侵蝕作用，其礦物組成的溶失使得帷幕細觀結(jié)構(gòu)發(fā)生改變，引起防滲性能衰減，可能造成滲流安全隱患。為反映帷幕防滲性能衰減機理及其造成的影響，本發(fā)明依據(jù)地下水動力學、溶質(zhì)運移以及化學動力學的相關(guān)理論，建立了多物理場耦合的數(shù)學模型，以易發(fā)生溶解且作為帷幕主要組成的Ca(OH)2作為目標組分，研究其溶解引起的帷幕防滲性能變化，以及由此引起的幕后滲流、離子濃度等的時空演變。本發(fā)明的模擬結(jié)果表明，模型能較好地反映帷幕孔隙度等細觀結(jié)構(gòu)的變化，溶蝕部位主要發(fā)生在帷幕頂部和底部，其中頂部以貫穿型侵蝕為主，而底部則主要為表層侵蝕。同時，模型還可預(yù)測由帷幕引起的壩基滲流、離子濃度等的變化。

通過使用本發(fā)明的方法模擬大壩基礎(chǔ)帷幕在某一時間段內(nèi)的防滲性能衰減值，包括帷幕的滲流速度和孔隙度等，這樣可以反映大壩基礎(chǔ)帷幕防滲性能在該時段內(nèi)的衰減過程及衰減程度，從而可以刻畫該時段內(nèi)壩基帷幕細觀結(jié)構(gòu)和防滲性能等隨時間演變及空間分布情況。

本發(fā)明提出的壩基帷幕防滲性能衰減的數(shù)值模擬方法具有良好的精度，為直觀、方便地了解壩基帷幕防滲性能的衰減發(fā)生的部位及程度提供可能，解決了因溶失引起帷幕防滲性能衰減的量化解析問題。本發(fā)明通過模擬帷幕中可溶組分溶解引起的帷幕孔隙度增大以及由此引起的防滲性能變化，可以為大壩的安全運行和維護提供依據(jù)從而減少因滲流引起的水壩失事事故。

附圖說明

下面結(jié)合附圖對本發(fā)明作進一步說明。

圖1是本發(fā)明實施例中待模擬壩基的結(jié)構(gòu)示意圖。

圖2為模擬結(jié)束后帷幕的孔隙度二維空間分布示意圖。

圖3為模擬結(jié)束后帷幕典型斷面孔隙度的一維空間分布示意圖。

圖4為模擬時段內(nèi)帷幕防滲性能衰減程度的斷面示意圖。

具體實施方式

實施例

本實施例的待模擬的壩基二維模型如圖1所示，大壩高150m，壩底寬60m，帷幕深100m，帷幕寬5m。其中，由ABCDEFIHG組成的區(qū)域Ω1為壩基基巖區(qū)域，由GHIF組成的區(qū)域Ω2為壩基帷幕。以大壩上游側(cè)底部與壩基相連接處作為坐標原點(x＝0m,y＝0m)，指向下游側(cè)為x軸正向，指向壩頂方向為y軸正向，建立坐標系。

本實施例中壩基礦物組成可視為惰性，并僅考慮帷幕組成中Ca(OH)₂的溶解作用。

本實施例的壩基帷幕防滲性能衰減的數(shù)值模擬方法，包括以下步驟：

1)根據(jù)待模擬的壩基的實際物理參數(shù)設(shè)置邊界條件，并得到模擬前所述壩基的滲流速度U₀、孔隙度θ₀、目標組分濃度C₀、滲透系數(shù)K₀、彌散系數(shù)D₀以及溶解速率R₀，邊界條件的設(shè)定為現(xiàn)有技術(shù)，不再贅述。所述壩基包括基巖和帷幕兩部分，所述目標組分濃度包括鈣離子濃度和氫氧根離子濃度。

模擬前所述壩基的孔隙度θ₀可由壩基的幾何參數(shù)直接得到。在實施時，限定區(qū)域Ω1和Ω2中的初始水頭為上游庫水位，在邊界條件的約束下對壩基進行滲流和溶質(zhì)運移的穩(wěn)定態(tài)模擬，獲得滲流穩(wěn)定態(tài)時的滲流場和化學場分布情況，從而得到模擬前所述壩基的滲流速度U₀、目標組分濃度C₀、滲透系數(shù)K₀、彌散系數(shù)D₀以及溶解速率R₀，即第1個子時段的初始時刻的滲流速度、孔隙度、目標組分濃度、滲透系數(shù)、彌散系數(shù)以及溶解速率。

在所述壩基的孔隙度θ₀已知的前提下，對壩基進行滲流穩(wěn)定態(tài)模擬為現(xiàn)有技術(shù)，不再贅述。

2)根據(jù)步驟1)設(shè)置的邊界條件，建立滲流方程、溶質(zhì)遷移方程和固相介質(zhì)時變方程；

所述滲流方程為：

$\frac{\∂}{\∂x}\left(K\frac{\∂H}{\∂x}\right)+\frac{\∂}{\∂y}\left(K\frac{\∂H}{\∂y}\right)=\mathrm{\μ}\frac{\∂H}{\∂t}$

式中，H為水頭，K為滲透系數(shù)，t為時間參數(shù)，μ為給水度，x、y分別為壩基的橫坐標、縱坐標。

滲流方程中邊界條件設(shè)定為：邊界AG和邊界DE為第一類邊界條件，即給定水頭邊界，分別等于上游庫水位和下游庫水位；邊界AB、BC、CD、EF、FG作為第二類邊界條件中的隔水邊界。

所述溶質(zhì)遷移方程為：

$\frac{\∂\left(\mathrm{\θ}C\right)}{\∂t}=\frac{\∂}{\∂x}\[D_L\frac{\∂\left(\mathrm{\θ}C\right)}{\∂x}\]+\frac{\∂}{\∂y}\[D_T\frac{\∂\left(\mathrm{\θ}C\right)}{\∂y}\]-\frac{\∂\left(u_x\mathrm{\θ}C\right)}{\∂x}-\frac{\∂\left(u_y\mathrm{\θ}C\right)}{\∂y}+\mathrm{\θ}\·R$

式中，R為可溶組分的溶解速率，C為目標組分濃度；滲流速度U分解為水平方向的滲流速度u_x和垂直方向的滲流速度u_y；θ為孔隙度；彌散系數(shù)D中涉及到的縱向彌散度α_L＝40m和橫向彌散度α_T＝5m；x、y分別為壩基的橫坐標、縱坐標；t為時間參數(shù)。

溶質(zhì)遷移方程中邊界條件設(shè)定為：邊界AG為第一類邊界條件，即給定濃度邊界，等于上游庫水中Ca²⁺濃度；邊界DE為第二類邊界，AB、BC、CD、EF、FG作為零通量邊界。

所述固相介質(zhì)時變方程為：

$\frac{\∂\mathrm{\θ}}{\∂t}=R\·M\·\mathrm{\θ}$

式中，M為Ca(OH)₂的摩爾體積，M＝33cm³mol^-1，R為可溶組分的溶解速率，θ為孔隙度，t為時間參數(shù)；

固相介質(zhì)時變方程中邊界條件設(shè)定為：邊界GH、HI、IF、FG為第二類邊界條件。

3)將待模擬時段劃分為n個子時段，其中第i個子時段的初始時刻為終了時刻為1≤i≤n，n≥2；本實施例模擬大壩在運行一百年期間帷幕防滲性能的衰減過程，令時間間隔為10天，則總的模擬時段可劃分為3650個子時段，即n＝3650。

根據(jù)步驟1)得到的所述壩基帷幕的滲流速度U₀、孔隙度θ₀、目標組分濃度C₀、滲透系數(shù)K₀、彌散系數(shù)D₀以及溶解速率R₀，對步驟2)建立的滲流方程、溶質(zhì)遷移方程和固相介質(zhì)時變方程進行求解，依次計算得到第j個子時段的終了時刻的滲流速度、孔隙度以及目標組分濃度，1≤j≤n；同時利用滲透系數(shù)K與孔隙度θ之間的函數(shù)、彌散系數(shù)D與滲流速度U、孔隙度θ之間的函數(shù)以及可溶組分的溶解速率R與目標組分濃度C之間的函數(shù)，得到第j個子時段的終了時刻的滲透系數(shù)、彌散系數(shù)以及溶解速率；并將第j個子時段的終了時刻的滲流速度、孔隙度、目標組分濃度、滲透系數(shù)、彌散系數(shù)以及溶解速率作為第j+1個子時段的初始時刻的相應(yīng)值；

其中，所述滲透系數(shù)K與孔隙度θ之間的函數(shù)為：

K＝K₀·(θ/θ₀)³[(1-θ₀)/(1-θ)]²

式中，K₀和θ₀分別為模擬前所述壩基帷幕的滲透率和孔隙度；

所述彌散系數(shù)D與滲流速度U、孔隙度θ之間的函數(shù)為：

D＝α·|U|+D_m·θ

式中，α為彌散度，D_m為鈣離子在水中的擴散系數(shù)，D_m＝1e^-9m²s^-1；

所述可溶組分的溶解速率R與目標組分濃度C之間的函數(shù)為：

$R=A\·{(1-\frac{C_{{\mathrm{Ca}}^{2+}}\·C_{{\mathrm{OH}}^{-}}^2}{K_{sp}})}^n$

式中，A為反應(yīng)動力學系數(shù)，n為動力學指數(shù)，K_sp為Ca(OH)₂的溶度積常數(shù)，為鈣離子濃度，為氫氧根離子濃度；A，n，K_sp均為預(yù)設(shè)系數(shù)。本實施例中n取4.5，A取1×10^-8mol·L^-1·s^-1，K_sp取1×10^-5.33。

模擬結(jié)束后，參照圖2可以看出，帷幕溶解主要發(fā)生在兩個以下兩個部位：1)帷幕淺層，即帷幕上部靠近壩體部位，帷幕溶失現(xiàn)象較為明顯，自上游側(cè)至下游側(cè)孔隙度皆相對較大，形成一個較為明顯的“貫穿型”滲流通道，且上游側(cè)溶失面積較大；2)帷幕底部溶失也較為明顯，但較頂部相對較弱，未形成自上游至下游的貫穿型滲流通道，僅沿帷幕表面向內(nèi)部發(fā)生淺層溶蝕。從圖2可以看出，模擬的結(jié)果與帷幕溶解的實測結(jié)果非常接近。

另外，分別選取帷幕內(nèi)x＝1m、x＝2.5m和x＝4m處作為帷幕靠近上游側(cè)、中心線和靠近下游側(cè)的典型斷面，作模擬時段結(jié)束時此三條斷面處帷幕孔隙度的分布，如圖3所示可以看出：1)帷幕頂部和底部處孔隙度相對較大，而帷幕中部孔隙度相對較?。?)靠近帷幕頂部部位處，無論是上游側(cè)、中心線還是下游側(cè)，孔隙度皆接近0.2，表明此部位的水泥結(jié)石已全部溶失，且自上游至下游側(cè)溶失面積逐漸減??；3)帷幕中部，孔隙度自上游側(cè)值下游側(cè)逐漸減小，表明帷幕中部所受溶蝕自上游至下游逐漸減?。?)帷幕底部部位，上、下游側(cè)的孔隙度明顯大于帷幕中心線位置。

本實施例還可以計算獲得帷幕防滲性能的衰減程度，選取帷幕中心線x＝2.5m處，作y＝0～-250m的典型斷面，通過該斷面的水流量如圖4所示，通過壩基典型斷面處的流量隨時間呈明顯增大的趨勢，初始時刻通過斷面流量為3.261356×10^-4m²/s，而在模擬時段末期流量為3.818986×10^-4m²/s，較初始時刻增大了約17.1％。

本發(fā)明不局限于上述實施例所述的具體技術(shù)方案，除上述實施例外，本發(fā)明還可以有其他實施方式。凡采用等同替換形成的技術(shù)方案，均為本發(fā)明要求的保護范圍。