基于非線性耦合的黏土心墻壩滲透系數(shù)的反分析

基于非線性耦合的黏土心墻壩滲透系數(shù)的反分析

0壩體材料參數(shù)的反演分析中國水庫超過90%是水庫。水庫是用松散的土壤、沙子和石頭填充的儲水建筑物。由于水庫蓄水形成的上游和下游水位差，水流將不可避免地從水庫和水庫底部滲透到下游。主要問題是洪水破壞，這對水庫的安全非常重要。而且很多土壩蓄水運行多年,在長期運行中,壩體結(jié)構(gòu)形態(tài)不斷進行調(diào)整。因此有必要對壩體主要參數(shù)進行反演分析,以評價大壩的工作性態(tài)。有的學(xué)者以土石壩滲流場中各測點水位的實測值與計算值的最優(yōu)化擬合準則,建立滲透系數(shù)的反演模型,并以青山水庫主壩滲流分析為例,根據(jù)多個滲流觀測點數(shù)據(jù),反演出大壩各個區(qū)域的滲透系數(shù);有的學(xué)者采用遺傳算法并結(jié)合Biot固結(jié)有限元數(shù)值法對一假設(shè)路基的雙層地基的一些主要土性參數(shù)進行反演;文獻基于遺傳神經(jīng)網(wǎng)絡(luò)的方法對一重力壩進行了滲透系數(shù)的反演。以往,在分析滲流場時通常忽視壩體應(yīng)力應(yīng)變對其滲透水流的作用,在分析應(yīng)力應(yīng)變時則將滲透水流作為外荷載加以處理。事實上,壩體的應(yīng)力應(yīng)變將改變土體的孔隙比,從而使壩體各部位的滲透系數(shù)發(fā)生變化,影響其滲流場的分布;與此同時,滲透水流的改變,其滲透力和浮力又影響著壩體的應(yīng)力應(yīng)變狀態(tài)。所以,在對大壩的有關(guān)材料參數(shù)進行反演計算時必須考慮應(yīng)變場和滲流場的共同作用。本文提出的基于應(yīng)變場和滲流場耦合的滲透系數(shù)的反演方法能夠更真實的反映大壩的工作狀態(tài),為科學(xué)評價大壩安全性提供有力的依據(jù)。1方程與接力學(xué)的基本方向和解算方法的研究1.1對復(fù)合方程的研究1.1.1土體結(jié)構(gòu)及滲透系數(shù)的分析在土顆粒及水不可壓縮的條件下,平面問題多孔介質(zhì)的滲流支配方程如下:??x[kx(εij)?h?x]+??z[kz(εij)?h?z]=St?h?t(1)??x[kx(εij)?h?x]+??z[kz(εij)?h?z]=St?h?t(1)式中:kx(εij),kz(εij))分別為土體水平、鉛直方向的滲透系數(shù),是應(yīng)變場的函數(shù);h為土體中任一點的水頭高度;St為貯水率。影響多孔巖土介質(zhì)的滲透性能的因素主要有兩個方面,一方面來源于土體中流體的性質(zhì),包括流體的密度和粘度等;另一方面是土體骨架的性質(zhì)。經(jīng)過大量的實驗研究表明,黏性土的滲透系數(shù)的對數(shù)lgk與孔隙比e之間存在線性關(guān)系,具體的表達式為:e=a+blgk(2)a=10b(3)b=0.01Ιp+δ(4)e=a+blgk(2)a=10b(3)b=0.01Ip+δ(4)式中:a、b為常數(shù),且都與塑性指數(shù)有關(guān);Ip為土體的塑性指數(shù);δ為與土體的類型有關(guān)的常數(shù),其平均值一般可以取0.05。1.1.2土體本構(gòu)模型的雙場耦合在滲流作用下,以張量形式表示的土體的平衡方程為:σij,j+fi(h)+Xi=0(5)σij,j+fi(h)+Xi=0(5)式中:σij,j為σij對空間坐標的偏導(dǎo);σij為有效應(yīng)力張量場;fi(h)為體積力,是滲流場水力梯度的函數(shù),與壩體外的滲透水頭有關(guān);Xi為荷載張量。將滲流場數(shù)值模型、應(yīng)力場數(shù)值模型及滲流與應(yīng)變的經(jīng)驗關(guān)系組合在一起,并以矩陣的形式表示,可得雙場耦合的數(shù)值模型為:[k]{h}+[s]{?h?t}={f}(6){dσ}=[D]{dε}=[D][B]{dw}(7)[Κ]{w}={F}+{X}(8)[k]{h}+[s]{?h?t}={f}(6){dσ}=[D]{dε}=[D][B]{dw}(7)[K]{w}={F}+{X}(8)式中:[k]為與滲透系數(shù)有關(guān)的矩陣;[S]為與貯水率St有關(guān)的貯水矩陣;{f}為已知常數(shù)項,由已知水頭結(jié)點得出;{dσ}為與土體自重有關(guān)的應(yīng)力增量矩陣;[D]為彈性矩陣(或彈塑性矩陣);[B]為應(yīng)變矩陣;[K]為剛度矩陣;{w}為結(jié)點位移;{F}為與滲透坡降有關(guān)的滲透力矩陣;{X}為荷載矩陣。上述式(6)描述應(yīng)力場影響下的滲流場;式(7)描述土體的應(yīng)力一應(yīng)變關(guān)系,可以選擇不同的土體本構(gòu)模型代入,式(8)描述滲流場影響下的力的平衡方程。當(dāng)采用Duncan模型時,式(2)和式(6)-(8)就構(gòu)成了土體在滲流作用下的非線性彈性雙場耦合模型,給定邊界條件和初始條件,便可確定滲流場和應(yīng)力場的耦合效應(yīng)。1.2間接耦合法與直接耦合法在雙場耦合分析方面,目前大致有兩類方法,一是將兩場分開計算,然后,通過兩場的交叉迭代達到耦合的目的,即間接耦合法;另一種是通過建立以應(yīng)變場和滲流場為未知值的數(shù)學(xué)模型,通過求得解析解達到完全耦合的目的,即直接耦合法。根據(jù)本文所研究對象的復(fù)雜性和特點,聯(lián)合求解上述式(6)-(8)有較大的困難,故采用間接耦合法。2基于合作土壤滲透系數(shù)的解決方案2.1土體初始孔隙率的確定由于滲透系數(shù)k是量級較小的數(shù)字(如10-8～10-4),為了給下面的反演計算提供方便,同時提高反演的精度,可由式(2)通過數(shù)學(xué)推導(dǎo)得到:k=c?k0?10αn0+εv1-n0-εv(9)式中:n0為土體的初始孔隙率;k0為土體的初始滲透系數(shù);εv為土體的體積應(yīng)變;α為與土體的塑性指數(shù)有關(guān)的常數(shù);c為由試驗確定的常數(shù)。這樣,α和v即為反演的主要參數(shù)。2.2土石壩滲流量反演模型根據(jù)上面的控制方程和定解條件,并將滲流量數(shù)據(jù)和浸潤線數(shù)據(jù)結(jié)合起來,一起作為反演的控制條件,同時對控制條件進行修正,可建立基于耦合的土石壩的滲透系數(shù)的反演模型:F(α,c)=ω11nn∑i=1|hh*-1|+ω2|QQ*-1|=min(10)式中:h是測壓管水位預(yù)測值;h*是測壓管水位觀測值;Q是滲流量預(yù)測值;Q*是滲流量觀測值;ω1和ω2是權(quán)重。2.3遺傳算法程序參數(shù)反演采用正反分析法,即把參數(shù)反分析問題轉(zhuǎn)化為一個目標函數(shù)的尋優(yōu)問題。問題的正算采用滲流場計算的有限單元法,反算用遺傳算法。土的本構(gòu)模型選用鄧肯-張模型,土壩壩體的應(yīng)力和滲流計算都采用有限元法。將反演程序GA與土壩的穩(wěn)定滲流程序正演程序SEEP和非線性有限元分析程序FEACD相結(jié)合,構(gòu)成基于流固耦合的土石壩土體滲透系數(shù)反分析程序。GA作為整個計算的主程序,SEEP和FEACD是其中的兩個子程序。程序采用Fortran語言編寫,FortranPowerstation4.0編譯,在Windows9.0以上操作環(huán)境下運行。程序的運行過程:耦合分析的初始應(yīng)力場和初始變形場選用土壩竣工時的應(yīng)力場和變形場,土壩竣工時的滲透系數(shù)k0作為初始滲透系數(shù),由k0確定的滲流場即為初始滲流場。首先由遺傳算法程序GA隨機產(chǎn)生幾組參數(shù)值,調(diào)用基于耦合的計算的子程序即滲流計算程序SEEP和應(yīng)力應(yīng)變計算程序,用這幾組參數(shù)值進行應(yīng)變場和滲流場耦合分析,通過迭代算出在各組滲透系數(shù)下浸潤線上各點的縱坐標值hi和滲流量Qi,然后返回GA程序中,用func=ω11nn∑i=1|hh*-1|+ω2|QQ*-1|計算各組參數(shù)的適應(yīng)度。其中h*是浸潤線上各點縱坐標的觀測值,即測壓管水位;計算浸潤線和觀測浸潤線上各點的橫坐標對應(yīng)相等。以適應(yīng)度為控制條件,對這幾組參數(shù)進行選擇、交換、突變等操作,得到新的幾組參數(shù),再調(diào)用SEEP和FEACD子程序。如此反復(fù)循環(huán),使待反演參數(shù)的適應(yīng)度逐漸增高,使?jié)B透系數(shù)逐漸逼近于真實滲透系數(shù)。3工程實例3.1壩體材料特性分析某水利樞紐主壩為黏土心墻壩,主壩壩頂高程為659.00m,最大壩高39.00m,壩頂長458.00m。主壩基礎(chǔ)主要地質(zhì)自上而下分別是凝灰風(fēng)化頁巖、凝灰沙礫巖與凝灰?guī)r。風(fēng)化巖層深淺不一,兩岸深約10m,河床部分深約5m。施工中,已將風(fēng)化巖層挖除,用黏土夯填,并將左岸邊坡削成斜坡。1961年10月通過竣工驗收,于1998年9月1日實行自動化觀測。在主壩共布置了3個觀測斷面,共13支測壓管,用于觀測壩體浸潤線,其中主壩0+041m斷面4支測壓管,編號為A1～A4;主壩0+066m斷面5支測壓管,編號為A5～A9。本文以典型斷面0+066為計算斷面。斷面及測壓管布置見圖1。典型時日水位測值和多年平均水位測值見表1。本文不考慮由于流固耦合引起的土體的力學(xué)參數(shù)的變化,假設(shè)土體的有關(guān)的力學(xué)參數(shù)是已知的常數(shù),主要反算和土體的滲透性能有關(guān)的參數(shù)。為了減少反算參數(shù)數(shù)量,提高計算精度,將對浸潤線影響較小的壩殼滲透系數(shù)作為已知條件,主要反演心墻的滲透系數(shù)。從式(9)的推導(dǎo)過程得反算參數(shù)的取值區(qū)間為α∈[4,6],c∈[1,10],公式(10)中權(quán)重ω1和ω2都取為0.5。壩殼的滲透系數(shù)為1×10-4cm/s,黏土心墻的初始滲透系數(shù)為1×10-6cm/s。通過大型三軸試驗,得到壩體材料的計算參數(shù)如表2。選擇初始種群規(guī)模為20,進化代數(shù)400代時,目標函數(shù)為0.0028,計算出α=5.21,a=1.56。將反算結(jié)果帶入公式(9),基于耦合進行滲流場的有限元正向計算,計算結(jié)果見圖2。從圖2可以看出,用考慮耦合的反算算結(jié)果得到的正算值比較接近實測值。3.2反分析結(jié)果的統(tǒng)計檢驗反分析完成后還需對反演結(jié)果進行檢驗與評價,具體方法是將反演結(jié)果輸入有限元正分析模型,計算出2006年7月20日的水頭計算值與該時刻測得的實際水頭值進行比較,然后檢驗二者的符合程度。采用后驗差法對反分析結(jié)果進行統(tǒng)計檢驗。假設(shè)實測水頭序列為hm(i),預(yù)測水頭序列為hc(i),二者之間的殘差為:ε(i)=hm(i)-hc(i)(i=1,2,?,n)(11)式中:n為特征點的個數(shù)。若S1為量測位移的均方差,S2為殘差的均方差,則有:S21=1nn∑i=1(hm(i)-ˉh)2(12)S22=1nn∑i=1(ε(i)-ˉε)2(13)其中:ˉh=1nn∑i=1hm(i)(14)ˉε=1nn∑i=1ε(i)(15)后驗比值c=s1s2(16)誤差概率p=p{|ε(i)-ˉε|＜0.6475S1}(17)根據(jù)上述公式計算得到后驗比值c為0.466,誤差概率p為0.894,說明反分析結(jié)果是有效的、合格的。4滲透特性的參數(shù)本文通過公式推導(dǎo)得到應(yīng)變場和滲流場耦合的關(guān)系式,可以避開對量級較小的滲透系數(shù)的求解;通過比較對土壩的滲流的正向計算采用間接耦合法,基于耦合,用遺傳算法對土壩的滲透特性的有關(guān)參數(shù)進行反演,通過實例計算得知