1種改進的流變模型在層狀覆蓋層壩基上的應用*

王 璠，毛海濤,，侍克斌，王曉菊

(1.新疆農業大學 水利與土木工程學院，新疆 烏魯木齊 830052；2.山西農業大學 城鄉建設學院，山西 太原 030000)

0 引言

河谷中的第四紀松散沉積物主要成分為顆粒粒徑較大的漂卵礫石、塊碎石以及粉細砂等，粗粒土分布相對更加廣泛，且呈強弱相間的層狀分布，深度從幾十米到幾百米不等。在這類地基上建壩，地基承載力不足、穩定性難以滿足要求、基礎沉降和不均勻沉降過大等是常見問題。除了施工措施原因外，忽略整體覆蓋層尤其是粗粒土層的流變性也是原因之一。

土石料的流變模型中主要有Maxwell模型[1]、Burgers模型[2]、Kelvin模型[3]、Bingham模型[4]、理想黏彈塑性體[5]、西原模型[6]等。這些模型通過基本元件的“串聯”和“并聯”來描述巖土介質的流變特性，建立巖土應力-變形-時間的本構模型，能反應巖土線性黏彈塑性性質，但不能表達巖土的加速蠕變特性；元件模型元件越多，模型越復雜，參數就越多，模型數值計算就越困難[7]。在以往研究巖土變形的流變模型中，H-K模型能較全面地反應巖土彈性、蠕變、松弛變形特性，如蔡新等[8]采用流變學理論，選取Maxwell模型和H-K模型模擬單一巖土料的流變特性，探討土石料流變特性對土石壩應力、位移的影響；李習平等[9]采用非線性擬合方法，比較了H-K模型和Burgers模型的擬合值與實驗數據的關系，驗證了H-K模型研究地基流變的適用性；胡剛等[10]將空區頂板-礦柱體系簡化為彈性矩形薄板H-K流變力學模型，來分析空區處理對圍巖及地表的影響規律。但H-K模型對于單一巖土體的流變研究具有優勢，而對性狀各異的多層土模擬結果往往與實際有較大出入，覆蓋層地基大多由多層巖土體組成，需要研究適合其變形的多層土流變模型。

綜上，本文對H-K模型進行改進，研究新的計算模型以適應層狀覆蓋層的特點，并以青灣壩為例，采用Comsol建立層狀覆蓋層上大壩流變及流固耦合數值模型，結合實測數據以確定模型的準確性，進而全面評價層狀覆蓋層流變對壩體、壩基以及主要結構的影響，以期為深厚覆蓋層上大壩的穩定計算提供支持。

1 層狀覆蓋層流變模型

1.1 單一土層的H-K模型

H-K模型從其元件構造上能直觀地識別彈性和黏彈性變形分量，其數學表達式能客觀地描述蠕變、應力松弛及穩定變形等。

該模型由1個彈模為E1的彈性元件(H)和Kelvin(K)體串聯而成如圖1(a)，其變形隨時間的關系如圖1(b)所示。在應力σ作用下，應力-應變關系如式(1)所示：

圖1 H-K的蠕變模型Fig.1 H-K creep model

(1)

式中：E1，e1分別為H彈簧體及Kelvin 體的彈性模量，Pa；η1為Kelvin 體的黏性系數，m2·s-1；σ是應力，Pa；而總應變ε=εH+εK，其中εH，εK分別為H彈簧體及Kelvin體的應變。

1)當?σ/?t=0，即應力σ不隨時間變化，由初始應變ε0=σ/E1可得，如式(2)所示：

(2)

由公式(2)可知，應變隨時間持續遞增(t→∞時，ε→σ/E1+σ/e1)，其應變速率由起始時的最大值逐漸趨近于零。

2)當?ε/?t=0，即應變ε不隨時間變化，由初始應力σ0=E1ε可得，如式(3)所示：

(3)

由公式(3)可知，應力隨時間增大而減小(t→∞時，σ→E1e1ε/(E1+e1))。

上述分析可知，H-K是時間上的衰減模型。模型中關鍵參數為E1，e1，η1，可通過室內外巖土體流變試驗獲得。此外，還與泊松比μ，滲透系數k等有關。對于層狀覆蓋層，各層巖土體特性差異較大，蠕變變形也各不相同，采用H-K模型不能較好反應覆蓋層整體蠕變特性。

1.2 多層土H-KS模型的建立

在層狀壩基中將每層土作為1個單元，考慮各層之間存在力的傳遞，各土層性質不同，因此，可將每層土視為1個H-K體，各層總體為串聯關系如圖2(a)所示。由H-K體串聯得層狀壩基模型如圖2(b)。

由于模型元件構造復雜不利于分析計算，結合模型中各元件所代表的含義，將圖2(b)中各層不隨時間變化的彈性模量E1、E2、…，Ei合并為模型整體彈性模量E，關系如式(4)所示：

圖2 多土層壩基串聯模型Fig.2 Multi-soil dam foundation series model

(4)

式中：E為彈性模量，Pa；Ei對應各層巖土變形初期的彈性模量，Pa。

得到簡化后的模型如圖3，其中E為與時間因素無關的整體變形模量，各層與時間因素相關的蠕變模型不變化。

圖3 簡化后層狀土層串聯模型Fig.3 Simplified series model of layered soil

簡化后形成層狀壩基流變計算新模型(以下簡稱H-KS模型)，在通過試驗獲取參數方面，該模型因試驗次數減少而降低了誤差，提高了計算結果精度。

總應變εt等于各H-K元件的應變之和，表達如式(5)所示：

(5)

式中：ε1、ε2、…、εn分別為對應各層巖土的變形量，m。

將式(2)代入式(5)得方程，如式(6)所示：

(6)

式中：σ為恒定應力，Pa；t為加載時間，s；Ei對應各層巖土變形初期的彈性模量，Pa；ηi對應各層黏彈變形的黏滯系數，m2·s-1。

聯立式(4)和式(6)可得到模型的變形方程，如式(7)所示：

(7)

式中：E為彈性模量，Pa。

壩基沉降除了受巖土體流變作用外，還受流固耦合效應的影響。因此，壩基沉降需要將流變與流固耦合模型聯合計算。

2 流固耦合模型

2.1 滲流偏微分方程

描述多孔介質飽和流體流動的達西定律質量守恒建立流體的連續性方程[11]，如式(8)所示：

(8)

式中：ρ為液體密度，kg·m-3；u為達西速度，m/d；ε?(ρεp)/?t為多孔存儲項，m/d；Qm為液體源匯項，m/d。

其中多孔存儲項表達如式(9)所示：

(9)

式中：S為多孔介質存儲系數；P為水頭，m。

達西速度表達如式(10)所示：

(10)

式中：k為滲透率，md；μ為液體動力黏度，N·s·m-2；Pf為孔隙壓力，Pa；g為重力加速度，m·s-2；D為位置水頭，m。

Qm為液體源匯項，由于孔隙空間的膨脹率(?εvol/?t)增加，使液體的體積分數增加，從而引起液體匯，如式(11)所示：

(11)

式中：αB為Biot系數；εvol為多孔介質體應變。

將式(9)～(11)代入式(8)可由質量守恒方程得式如式(12)所示：

(12)

式中：H為總水頭，m。

2.2 固體骨架變形方程

孔隙彈性將多孔介質中流體流動及顆粒變形聯系起來，得到應力、應變和孔隙壓力的本構關系，由Biot理論可得流體含量增量、體積應變和孔隙壓力的本構關系[12]，聯立如式(13)所示：

(13)

式中：σ是柯西應力張量，Pa；C為彈性矩陣；ε是應變張量；αB是Biot系數；Pf是流體孔隙壓力，Pa；I為單位矩陣；ζ為多孔介質中流體含量的變化值，m3。

固體顆粒的體積變化Pm的表達式如式(14)所示：

Pm=-Ksεvol+αBPf

(14)

式中：Ks為固體的彈性模量，Pa。

2.3 耦合的實現

采用COMSOL多場耦合來實現流固耦合，其具體流程如圖4所示。該流程能完成某一時刻的數值模擬，下一個時刻，應力張量可作為1個附加的各項同性項參與計算[13]，并視為下一時刻的1個初始應力，則該時刻總的初始應力為Pt1，其表達如式(15)所示：

圖4 壩基沉降過程數值模擬流程Fig.4 Numerical simulation procedure of dam foundation settlement process

Pt1=Pt0+σ

(15)

式中：Pt1、Pt0分別為t1及t0時刻的初始應力，Pa；σ為t0時刻計算的應力張量，Pa。

整個耦合過程在每一個時間段內反復迭代求解，直到沉降完成為止。

此外，為了解壩基各層流變對其應力應變及滲流的影響，將常用的非流變模型—分層總和法結合流固耦合模型(以下簡稱HSM模型)與本模型計算結果進行對比分析。

3 工程應用

3.1 工程概況

青灣均質土壩最大壩高50 m，壩長200 m，水庫正常蓄水高度45 m。上下游坡度均為1∶2。壩基覆蓋層約80 m，根據顆粒組成及物理力學性質等差異，可將其自下而上分為5層：1)以粗粒為主的冰水積含漂卵礫石層，厚11.5 m；2)以細粒為主的堰塞堆積粉細砂層，厚13 m；3)以粗粒為主的沖積含漂砂卵礫石層，厚17 m；4)以細粒為主的堰塞堆積粉細砂層，厚15 m；5)以中粗粒為主的沖洪堆積含漂砂卵礫石層，厚25 m。覆蓋層各土層的物理力學參數見表1。

表1 壩體及各層巖土物理力學屬性Table 1 Physical and mechanical properties of dam body and each layer of rock and soil

3.2 建模計算

在Comsol中通過自動劃分網格，得到1 086個域單元，239個邊界元，求解自由度為7 069個，具體見圖5。模型底部和兩側設置為固定約束且為不透水邊界；上游側為透水邊界，大壩上下游水位分別為45 m和0 m，淹沒在水位以下的上游河谷及壩面設置為總水頭邊界，其它邊界設置為潛在出滲邊界。

圖5 大壩數值模型網格剖分單元Fig.5 Grid division unit of dam numerical model

為了監測壩基沉降、應力和滲流等相關指標，大壩設置了完整的監測系統。為監測壩基沉降，在基建面及截水槽底部共設置5組水管式沉降儀；壩體位移監測包括水平和沉降2項，分別在壩底面、高度為20 m和35 m水平面安裝沉降和水平位移設備；在滲流監測方面，壩體內部設置測壓管并在壩后設置量水堰。此外，還在壩軸線(斷面2)、截水槽底部等重要部位設置應變及滲壓計等，大壩監測點布置如圖6。

圖6 壩體監測系統布置Fig.6 Layout of dam monitoring system

4 結果分析

4.1 大壩沉降分析

4.1.1 覆蓋層各土層沉降

采用HSM及H-KS模型分別對大壩沉降進行計算，其結果如圖7所示。

圖7 2種工況下沉降等直線Fig.7 Contour lines of settlement under two working conditions

由圖7可知，H-KS模型計算的沉降量總體上高于前者。壩基各土層的變形增量自上至下分別為11.3，2.5，17.8，4.3，22.5 mm。對比發現，以粗粒土為主的1)，2)，3)土層沉降增量較大，占總沉降量的88.4%，單位厚度增量分別為1.0，1.1，0.9 mm/m。其中5)層厚度較大且位于底部，其變形總量也相應較大，占總增量的38.5%；細粒土為主的土層2)，4)總體變形較小，沉降單位厚度增量僅為0.2，0.3 mm/m。由此可見，考慮流變后的沉降不僅與各土層的性質有關，還與其厚度及所在位置相關。

4.1.2 結果驗證

為驗證計算結果的準確性，結合圖6中壩基沉降監測數據對比結果如表2所示，并以斷面1位置為例繪制對比圖如圖8。

表2 壩基各層沉降量對比Table 2 Comparison on settlement of each layer of dam foundation mm

圖8 大壩監測斷面1處沉降對比Fig.8 Comparison on settlement of dam monitoring section 1

由表2和圖8可知，H-KS模型計算的沉降量在各土層及整體上均與實測值更接近，而HSM模型計算結果與實測值誤差較大。說明H-KS模型計算結果與實際沉降更吻合。對于該粗細相間的層狀壩基，各層土的沉降量在2種算法下也不盡相同。其中，HSM模型在1)，3)，5)粗粒土層中計算誤差較大，如SG1處誤差分別為27%，10.6%，16.1%；而細粒土層2)，4)計算結果誤差相對較小，誤差僅為為0.89%和2.6%；此外，H-KS模型在各土層的沉降計算誤差均在5%以內。

綜上，HSM模型在流變性較小的細粒土沉降計算中基本符合實際，但在粗粒土層壩基沉降計算方面與實際有較大出入，而考慮流變的H-KS模型計算結果與實際沉降基本吻合。

4.2 應力分析

壩基巖土流變對大壩整體強度會有一定影響，尤其是粗細相間的層狀壩基，各土層流變對大壩應力的影響尚不明確，需要進一步探索。

4.2.1 應力整體分布

2種模型下的大壩應力分布如圖9所示。

圖9 應力分布等值線Fig.9 Contour lines of stress distribution

由圖9可知，H-KS模型計算總體應力水平高于前者，壩體和壩基應力分布均是自上而下逐漸增大，并在壩基底部出現應力極值，分別為0.937，1.32 MPa。受壩基流變的影響，壩體整體應力水平有所增長，但增量率較小。值得注意的是，壩踵和壩趾處的應力均有明顯增大，增量達到0.027，0.03 MPa。

4.2.2 覆蓋層各土層應力變化

取壩軸線處各層應力計算結果如圖10，需要說明的是，由于壩軸線處的截水槽將1)層全部替換，替換后土料和均質土壩一致，為細粒黏性土，因此該處計算結果不反應原來的粗粒土。

圖10 壩軸線處覆蓋層各層應力變化Fig.10 Change in stress of each layer of overburden at dam axis

由圖10對比發現，以粗粒土為主的3)，5)層應力增量明顯高于以細粒土為主的1)，2)，4)層。自上而下各層的應力增量分別約為0.03，0.06，0.1，0.07，0.17 MPa。其中，性質單一、顆粒最細的1)層(截水槽處)應力差別最小，誤差在1.2%，2)，4)土層應力增長也較小，均在5.75%左右；而5)層差別最大，總體應力水平增長率約為15.2%。可見，在流變性大的粗粒土層應力方面增量也較大，而在流變性較小的細粒土層應力增量相對較小，說明土層流變性影響著應力變化。

4.2.3 應力結果驗證

結合壩體監測數據，取圖6中壩體中壩軸線位置應力監測值與2種模型計算結果進行對比如圖11所示。

圖11 壩軸線處壩體應力變化Fig.11 Change in stress of dam body at dam axis

由圖11可知，壩軸線處壩體應力自上而下變化規律不同，自壩頂向下16m應力值隨高度基本無變化，甚至出現變小趨勢，本模型應力值與實測值吻合且大于HSM模型計算結果；16～26 m處應力與高度呈非線性變化，應力由0.18 MPa增至0.23 MPa，2模型計算結果均與實測值能較好吻合；26 m至壩底應力與高度基本呈線性變化，擬合線性關系分別為：y=-91.1x+43.7、y=-74.8x+40.2，HSM隨著高度下降，偏差逐漸加大，H-KS模型計算結果與監測值基本一致。因此，本模型的正確性得到驗證，相對于傳統模型優勢顯著。

4.3 壩體變形和應力變化原因分析

綜合4.1～4.2計算結果，將2種模型得到的壩體應力、變形極值對比如表3。

表3 壩體應力變形極值對比Table 3 Comparison on extreme values of dam stress and deformation

H-KS模型壩體最大沉降較前者增大15.7%，最大應力值也增大10%。分析其原因，主要是覆蓋層各層持續的流變變形對壩體產生的影響。巖土體變形可分為2個階段：外力變形階段和流變變形階段。其中，外力變形階段大壩在施工碾壓過程中，主要受到上部碾壓作用壓迫，顆粒棱角逐漸破碎、發生位移，孔隙被填充和壓密。這一過程隨外力碾壓功的增加，壩基料孔隙逐漸壓密，但上部碾壓功影響有限，本階段巖土料變形并沒有徹底完成，大壩填筑完成后，壓縮也基本結束。流變變形階段的施工外力影響完畢，該階段變形主要是在自重及水荷載作用下的進一步變形，顆粒之間孔隙逐漸減小并密實，這是一個隨時間推移變形逐漸趨于穩定的過程。顯然，壩體的變形受到壩基覆蓋層變形的影響，傳統的HSM模型僅計算了外力變形階段，忽略了流變變形階段，故而壩體變形相對較小，這與實際情況是不相符的。H-KS模型計算結果相對偏大，但更加符合實際壩體的變形情況。

4.4 滲流分析

除了上述應力場的變化，覆蓋層中巖土體流變也影響著滲流場的改變，也是1個流固耦合過程。

4.4.1 大壩孔隙水壓力分布

分別計算HSM、H-KS模型得到大壩孔隙壓力分布如圖12所示。

圖12 孔隙水壓力等值線Fig.12 Contour lines of pore water pressure

由圖12可知，H-KS模型計算的各層土體孔壓均小于前者，其中流變性大的粗粒層1)，3)，5)孔壓減小量較大，平均減小量分別為0.039 MPa，0.043 MPa，0.06 MPa；而流變性較小的細粒料層2)，4)孔壓減小量相對較小，僅為0.013 MPa和0.015 MPa?？梢妿r土體流變對粗細土層的影響也不同，以5)，2)層為例，前者減小約50.6%，而后者僅減小2.48%。分析其原因，主要是各層土在流變效應下，土骨架被壓縮而孔隙體積減小，使原來由孔隙水承擔的部分壓力轉由土骨架承擔，孔隙水壓力減小。其中，粗粒土相對細粒土流變性較大，其孔隙水壓力減小量也更大。顯然，顆?？紫扼w積的變化加速了流固耦合過程，使壩基滲流場與應力場趨于穩定的時間變短。

4.4.2 滲流結果驗證

大壩滲流量實測結果與2模型計算結果對比如表4所示。

表4 壩后滲流總量對比Table 4 Comparison of total seepage amount behind dam

由表4可知，滲流穩定后采用HSM模型計算的滲流量比實測值增加了0.08 m3/s，增量為21.6%；采用H-KS模型計算的滲流量比實測值減少了0.03 m3/s，減小率為8.1%，與實測結果更吻合。

5 討論

5.1 深厚覆蓋層考慮流變的重要性

研究表明[14]，我國西南山區河谷覆蓋層自下而上分為1)含泥砂卵碎石層(Ⅰ)；2)漂卵石、含泥砂碎塊石、粉細砂互層(Ⅱ)；3)現代河流漂卵石層(Ⅲ)共3大層，覆蓋層中粗粒含量相對細粒占比要大得多。本文反映了由于粗粒土流變性較大[15]，大壩運行后其變形、應力、滲透性等變化主要受到粗粒土流變的影響。常規計算方法忽略了覆蓋層流變尤其是粗粒土流變的影響，結果與實際情況不相符。忽略覆蓋層流變會導致土石壩沉降等變形超過預期，壩體內部因不均勻沉降出現裂縫、結構脫落甚至破壞，影響壩體的安全穩定。深厚覆蓋層上大壩的出險加固難度相對較大，前期的設計和壩基處理若未考慮其流變性，也會為后期留下隱患。因此，在深厚覆蓋層上建壩，科學對待壩基的流變性并采取合理的工程措施是至關重要的。

5.2 流變對流固耦合的影響

計算結果表明，考慮壩基流變相比不考慮流變模型的沉降量增加了約75 mm，也就是說考慮了壩基流變性后，巖土體的實際孔隙壓縮較快且會更小，其滲透性更加弱化，增強了壩基的滲透穩定性。若不考慮流變，流固耦合過程和時間會更長，最終結果也有一定偏差。所以說，流變加速了流固耦合，縮短了流固耦合時間，滲流計算結果更加符合實際情況，為滲流控制與分析提供更合理的理論支持。

6 結論

1)對于巖性、物理力學性質差異較大的層狀覆蓋層壩基，采用H-KS流變模型能較好反應大壩應力、變形和滲流實際情況，與監測值對比，誤差在5%以內，模型在模擬層狀壩基方面具有明顯的優勢。

2)H-KS模型計算結果相比不考慮流變模型的應力、變形分別增大了10%和15.7%，部分增量會影響到大壩結構的安全穩定，應予以重視。

3)流變會導致大壩整體滲透性降低，粗細粒層孔隙水壓力分別減小了約50%和2.5%，加速了流固耦合過程，對于重新評估大壩的滲透性有重要參考意義。

4)由于流變性不同，深厚覆蓋層中粗粒土層對變形、應力和滲流方面影響最大，細粒土的貢獻相對較小；除了與巖性相關之外，還與土層所處位置、厚度等相關。