龔興勇，王銘明，楊雙鋼，胡圣明，茍 超

(1.昆明理工大學(xué) 電力工程學(xué)院，昆明 650500；2.湖南省水利水電勘測設(shè)計(jì)規(guī)劃研究總院有限公司，長沙 410007)

尾礦庫是礦山企業(yè)儲存尾礦渣及澄清尾礦水的重要生產(chǎn)設(shè)施[1]，我國90%以上的尾礦壩均采用方式簡單、管理方便的上游式堆筑。由于特殊的筑壩方式和筑壩材料，使得尾礦壩在地震荷載作用下非常敏感[2]，在地震作用下容易壩體失穩(wěn)，對下游居民的生命和財(cái)產(chǎn)造成威脅。因此研究尾礦壩在地震作用下的動力特性尤為重要。

地震的隨機(jī)性與突發(fā)性和壩體空間結(jié)構(gòu)的非均勻性及變異性，致使尾礦壩抗震問題十分復(fù)雜，國內(nèi)外學(xué)者主要就地震液化、動位移、加速度等方面進(jìn)行了研究。AMBRASEYS[3]采用剪切楔理論進(jìn)行尾礦壩動力分析，分析中只考慮了土體材料的剪切畸變，計(jì)算結(jié)果不夠準(zhǔn)確，只能看作近似。汪聞韶[4]分析了不同機(jī)理下的砂土液化，分析得出土體處于極限平衡狀態(tài)下將先于土體液化的發(fā)生。因此，本研究以某上游式尾礦壩為研究對象，借助ABAQUS有限元軟件運(yùn)用等效線性原理進(jìn)行尾礦壩動力分析，通過Fortran語言編寫孔壓應(yīng)力模型，對尾礦壩在地震作用下壩體的加速度、動位移、動應(yīng)力及孔隙水壓力進(jìn)行計(jì)算，通過計(jì)算結(jié)果采用有效應(yīng)力判別法對尾礦壩的液化區(qū)進(jìn)行判定。

1 靜動力分析基本理論

1.1 鄧肯—張E-B

靜力分析選用國內(nèi)外許多工程實(shí)例證明有效的鄧肯—張模型[5]。其切線模量Et為：

(1)

考慮土體材料強(qiáng)度的非線性，內(nèi)摩擦角可視為隨圍壓變化的函數(shù)：

(2)

切向體積模量：

(3)

鄧肯—張E-B模型公式中共有9個(gè)模型參數(shù)：Kb、K、n、m均為實(shí)驗(yàn)常數(shù)；c為黏聚力；Δφ為對φ進(jìn)行修正的系數(shù)；Rf為破壞比；Pa為大氣壓。

1.2 等效線性模型

等效線性模型在周期荷載下的滯回性是通過黏彈性Kelvin模型來反映的[6]，其應(yīng)力應(yīng)變關(guān)系為：

τ=Gγ+ηGγ

(4)

式中：τ為剪應(yīng)力；G為剪切模量；ηG為剪切黏滯系數(shù)；γ為剪應(yīng)變。

剪切黏滯系數(shù)ηG為：

(5)

式中：λ為阻尼比，ω為圓頻率。

動力反應(yīng)分析采用沈珠江提出的如下形式：

(6)

(7)

(8)

式中：γdmax為地震過程中的最大剪應(yīng)變；Pa為大氣壓。

1.3 孔壓應(yīng)力模型

孔壓應(yīng)力模型采用徐志英和沈珠江建議的公式[7]：

(9)

將式(9)對t求導(dǎo)，并寫成增量形式得：

(10)

式中：α為土體的試驗(yàn)參數(shù)，取3.0；σ′0為有效應(yīng)力；τ0水平剪應(yīng)力；Δug為Δt時(shí)間內(nèi)孔隙水壓力增量；ΔN為Δt時(shí)間內(nèi)的振動周數(shù)；N為累計(jì)振動周數(shù)；m為決定孔隙水壓力隨著α而遞減的常數(shù)，對于尾礦砂已求得m=1.1～1.2，本文取1.1；θ為試驗(yàn)常數(shù)，取0.7；Nt為α=0時(shí)的液化周數(shù)，其表達(dá)式為：

(11)

τav=0.65τmax

(12)

式中：τav為平均剪應(yīng)力；a、b為尾礦料抗液化性能有關(guān)的常數(shù)。計(jì)算時(shí)段等效周數(shù)ΔN和累計(jì)等效周數(shù)Neq采用MARTIN等[8]研究方法，震級和振動持續(xù)時(shí)間根據(jù)表1計(jì)算得出。

表1 震級、等效振動次數(shù)和振動持續(xù)時(shí)間的關(guān)系Table 1 The relationship between magnitude，time duration，and equivalent cyclic number

1.4 液化判別法

尾礦壩的液化判別方法主要有剪應(yīng)力對比法、循環(huán)應(yīng)力法、有效應(yīng)力法等。本研究選用許多工程實(shí)例證明有效的有效應(yīng)力法來作為液化判別的標(biāo)準(zhǔn)。計(jì)算公式如下：

(13)

對于三維問題，震前平均主應(yīng)力的計(jì)算采用下式：

(14)

式中：σx、σy、σz為震前應(yīng)力分量；uw為孔隙水壓力。計(jì)算得出ks>1，判定該單元完全發(fā)生液化；反之，該單元不液化。

1.5 計(jì)算步驟

1)采用鄧肯—張E-B模型靜力計(jì)算，根據(jù)靜力結(jié)果算出各單元的最大剪切模量Gmax，取初始阻尼比λ0=0.05。

2)由初始剪切模量Gmax和初始阻尼比λ0形成各單元的剛度矩陣[K]、質(zhì)量矩陣[M]、阻尼矩陣[C]，建立動力平衡方程。

3)計(jì)算整個(gè)壩體的圓頻率ω。

將輸入地震時(shí)程劃分為若干個(gè)時(shí)段Δt。

5)將計(jì)算出的G/Gmax和λ/λmax重新代入進(jìn)行數(shù)值分析，一般迭代3～4次即可收斂。

6)計(jì)算該時(shí)段Δt內(nèi)的液化周數(shù)Nl，根據(jù)表1計(jì)算等效周數(shù)ΔN和累計(jì)等效周數(shù)Neq。

7)按式(10)計(jì)算該時(shí)段的孔隙水壓力增量。

8)如此一個(gè)時(shí)段接著一個(gè)時(shí)段計(jì)算，直至地震結(jié)束。

2 工程概況及有限元模型建立

2.1 工程概況

尾礦壩采用上游式筑壩方式，初期壩透水性良好，初期壩頂寬4 m，壩長84 m，上下游坡度均為1∶2，壩高20 m。堆積子壩外邊坡坡比1∶4.8，干灘坡比1∶100，考慮正常工況條件，干灘長度取300 m，堆積子壩高為85 m，整個(gè)尾礦壩最終壩高為105 m，堆積壩材料分區(qū)依次為尾中砂、尾細(xì)砂、尾粉砂，具體的材料分區(qū)如圖1所示，三種材料的滲透系數(shù)分別為1.25×10-5、1.85×10-6、1.25×10-6m/s。

2.2 有限元模型建立

2.2.1 邊界條件

在計(jì)算中根據(jù)建設(shè)的實(shí)體模型，固定壩體底部三個(gè)方向的位移，約束其他各面法向位移。在滲流分析時(shí)，設(shè)初期壩為透水邊界，按正常水位考慮干灘長度300 m，在上游邊界指定隨高程線性變化水頭邊界。

2.2.2 網(wǎng)格劃分

壩體共劃分2 720個(gè)單元，包括4 103個(gè)節(jié)點(diǎn)。尾礦壩單元?jiǎng)澐智闆r詳見圖1。

圖1 尾礦壩計(jì)算模型Fig.1 Calculation model of tailings dam

2.2.3 荷載

計(jì)算中所用的荷載主要有自重荷載、水壓力和地震荷載，動力時(shí)程分析中采用設(shè)計(jì)院提供的規(guī)范譜地震波，該地震波場地烈度為7級，并對其地震加速度進(jìn)行修改，將地震波加速度峰值調(diào)整為0.148 g，經(jīng)研究認(rèn)為順河向地震對壩體造成危害最大[9]，因此在下面中值輸入順河向加速度。地震加速度時(shí)程曲線如圖2所示。地震歷時(shí)30 s，整個(gè)地震過程1 s一個(gè)時(shí)段，共有30個(gè)時(shí)段，每個(gè)時(shí)段采用Wilson-θ法逐步積分，積分步長取為0.01 s，每段地震數(shù)據(jù)有100個(gè)，共有3 000個(gè)數(shù)據(jù)。計(jì)算參數(shù)見表2和表3。

圖2 地震波加速度時(shí)程曲線Fig.2 Duration curve of seismic acceleration

表2 鄧肯-張E-B模型參數(shù)Table2 Duncan-Chang E-B model parameters

3 主要計(jì)算結(jié)果

3.1 靜力計(jì)算結(jié)果

尾礦壩靜力計(jì)算結(jié)果見圖3～6。通過圖3、圖4應(yīng)力等值線圖可看出，尾礦壩的最大主應(yīng)力和最小主應(yīng)力均為負(fù)值，最大值出現(xiàn)在靠近壩基處。壩體無受拉區(qū)域，土體單元始終為受壓狀態(tài)。通過圖5可得到，最大應(yīng)力水平為0.84，其最大值出現(xiàn)在尾細(xì)砂與尾粉砂的交界處，最大值小于1，表明未出現(xiàn)剪切破壞。通過圖6可得到，孔隙水壓力為0的面為浸潤面，等值線圖可看出浸潤線的逸出位置在初期壩壩腳，沒有從堆積壩壩坡逸出，由于該初期壩屬于透水壩，從壩腳逸出是符合實(shí)際的。綜上，該尾礦壩壩體靜力穩(wěn)定。

圖3 最大主應(yīng)力(單位：kPa)Fig.3 The maximum principal stress(unit：kPa)

圖4 最小主應(yīng)力(單位：kPa)Fig.4 The minimum principal stress(unit：kPa)

圖5 壩體應(yīng)力水平Fig.5 Stress level of dam body

圖6 孔隙水壓力分布(單位：kPa)Fig.6 Pore water pressure distribution(unit：kPa)

3.2 動力響應(yīng)分析

3.2.1 模態(tài)分析

本文開發(fā)的等效線性模型參數(shù)中，ω取結(jié)構(gòu)的圓頻率，因此首先對尾礦壩進(jìn)行模態(tài)分析，提取壩體的頻率，圖7為尾礦壩第一階振型(基本振型)圖。

圖7 尾礦壩第一階振型(基本振型)圖Fig.7 First order vibration pattern(basic vibration pattern)of tailings dam

尾礦壩壩體的自振頻率計(jì)算結(jié)果如表4所示。

表4 壩體自振頻率Table 4 Self-vibration frequency of dam body

研究表明：壩體動力響應(yīng)的強(qiáng)弱并不完全取決于壩基地震的輸入，還與壩體自身結(jié)構(gòu)性質(zhì)有關(guān)，壩體結(jié)構(gòu)的自振頻率和輸入地震波頻率越接近，說明壩體的動力反應(yīng)越強(qiáng)烈。從表4中可以看出壩體的自振頻率與輸入地震波的頻率相差較大，即壩體的動力響應(yīng)不是很強(qiáng)烈。

3.2.2 動位移反應(yīng)

壩體動位移反應(yīng)如圖8～9所示。堆積壩壩體水平向最大動位移0.47 m，垂直向最大動位移為0.037 m。堆積壩最大動位移均發(fā)生在壩頂附近，且從壩頂向下動位移逐漸減小。從動位移反應(yīng)分布來看，水平向動位移相對垂直向動位移較大。

圖8 水平向位移(單位：m)Fig.8 Horizontal displacement(unit：m)

圖9 垂直向位移(單位：m)Fig.9 Vertical displacement(unit：m)

3.2.3 加速度反應(yīng)分析

垂直壩軸線中心剖面順河向加速度等值線見圖10，從圖中可以看出水平向加速度分布不具有一定的規(guī)律，加速度最大值出現(xiàn)在干灘位置，最大加速度5 m/s2，相對于輸入地震波峰值0.148 g，放大系數(shù)為3.38。

圖10 中心剖面水平方向加速度等值線圖Fig.10 Horizontal acceleration contour of central section

為了反映尾礦壩加速度隨時(shí)間的變化趨勢，通過計(jì)算得出中心截面初期壩壩頂節(jié)點(diǎn)、h/2壩高節(jié)點(diǎn)、3h/4壩坡節(jié)點(diǎn)、堆積壩壩頂節(jié)點(diǎn)四個(gè)典型位置的順河向加速度的時(shí)程曲線，如圖11所示。

圖11 加速度時(shí)程曲線Fig.11 Acceleration time history curve

由圖11可以看出尾礦壩地震反應(yīng)水平向最大加速度出現(xiàn)在5～10 s，即輸入地震荷載歷時(shí)的強(qiáng)震段，符合壩體受力的一般規(guī)律。初期壩壩頂節(jié)點(diǎn)最大值2.1 m/s2，放大系數(shù)為1.4；h/2壩高節(jié)點(diǎn)最大值2.7 m/s2，放大系數(shù)為1.8；3h/4壩高節(jié)點(diǎn)最大值3.4 m/s2，放大系數(shù)為2.2；堆積壩壩頂節(jié)點(diǎn)最大值2.8 m/s2，放大系數(shù)為1.89。其中最大值出現(xiàn)在3h/4壩高節(jié)點(diǎn)處，之后再往上加速度峰值有所減小。

3.2.4 動應(yīng)力分析

尾礦壩體動應(yīng)力反應(yīng)如圖12、13所示。從壩體中心剖面上動應(yīng)力分布可知，堆積壩壩體最大動主應(yīng)力為572 kPa，最小動主應(yīng)力為260 kPa，動主應(yīng)力分布均勻且全為壓應(yīng)力，最大動主應(yīng)力和最小動主應(yīng)力從上到下逐漸增大，極大值出現(xiàn)在壩基附近。

圖12 最小主應(yīng)力(單位：kPa)Fig.12 The minimum principal stress(unit：kPa)

圖13 最大主應(yīng)力(單位：kPa)Fig.13 The maximum principal stress(unit：kPa)

3.2.5 地震液化分析

需要說明，在浸潤線以上的區(qū)域?yàn)榉秋柡蛥^(qū)，可直接判定為不液化，浸潤線以下的區(qū)域通過孔隙水壓力和震前平均應(yīng)力的關(guān)系來判斷。為了更加直觀地觀察尾礦壩是否發(fā)生液化及判斷液化區(qū)域。本文基于FORTRAN語言編寫孔壓應(yīng)力模型，計(jì)算地震過程中壩體單元的振動孔隙水壓力。通過MATLAB語言編寫后處理程序，生成不同地震時(shí)程下的壩體液化區(qū)域見圖14所示。

圖14 不同時(shí)刻孔壓比等值線云圖 Fig.14 Contour cloud map of pore water pressure ratio at different times

為了更為直觀地判斷液化區(qū)域，圖中只給出了ks≥1和ks<1的兩個(gè)區(qū)域，圖中可以看出黃色區(qū)域?yàn)橐夯瘏^(qū)，藍(lán)色區(qū)域?yàn)榉且夯瘏^(qū)。液化結(jié)果顯示：尾礦壩在地震發(fā)生3 s時(shí)局部發(fā)生液化，液化主要發(fā)生在上游積水淺層區(qū)域，其主要原因?yàn)樯嫌畏e水淺層區(qū)域砂土處于飽和狀態(tài)，上覆有效應(yīng)力較其他區(qū)域小，在地震荷載作用下孔隙水壓力在較短時(shí)間達(dá)到上覆有效應(yīng)力，導(dǎo)致土體液化。0～18 s液化區(qū)隨著地震時(shí)間的延長而增大，18～30 s液化區(qū)范圍有小幅度的增加。綜上所述，尾礦壩連續(xù)液化單元所形成的液化區(qū)主要集中在浸潤線以下的沉積灘淺層區(qū)域，在壩頂和下游未出現(xiàn)明顯的液化區(qū)域，且只在沉積灘區(qū)域出現(xiàn)小范圍的液化，未形成滑移通道，但是小范圍的液化區(qū)域也會對尾礦壩的穩(wěn)定性造成不利影響。根據(jù)尾礦壩抗震災(zāi)害調(diào)查[10]，尾礦壩沉積灘處尾礦料松散且長期處于飽和狀態(tài)，在動荷載作用下容易發(fā)生液化，這表明本次計(jì)算結(jié)果合理。

4 結(jié)論

對尾礦壩進(jìn)行靜動力分析及壩體液化分析，得出以下結(jié)論：

1)靜力分析計(jì)算得出最大主應(yīng)力和最小主應(yīng)力均為壓應(yīng)力，無拉應(yīng)力出現(xiàn)；應(yīng)力水平最大值為0.84，小于1表明未出現(xiàn)剪切破壞；浸潤線逸出位置在初期壩壩腳，未從堆積壩壩坡逸出。故尾礦壩靜力穩(wěn)定。

2)在荷載作用下，壩體單元均處于受壓狀態(tài)；位移主要發(fā)生在壩頂附近和積水區(qū)域，順河向最大位移0.47 m，垂直向最大位移0.037 m。壩體整體位移較小。四個(gè)典型位置在地震作用下加速度呈現(xiàn)相同的變化趨勢，最大值均發(fā)生在強(qiáng)震階段，垂直方向的振動相對水平方向反應(yīng)較小，在該地震作用下反應(yīng)不大。

3)尾礦壩在地震發(fā)生3 s時(shí)局部發(fā)生液化，液化主要發(fā)生在上游積水淺層區(qū)域，其主要原因在于上游積水淺層區(qū)域砂土處于飽和狀態(tài)，上覆有效應(yīng)力較其他區(qū)域小，在地震荷載作用下孔隙水壓力在較短時(shí)間達(dá)到上覆有效應(yīng)力，導(dǎo)致土體液化。0～18 s液化區(qū)隨著地震時(shí)間的延長而增大，18～30 s液化區(qū)范圍有小幅度的增加。該尾礦壩在地震時(shí)只出現(xiàn)小范圍的液化，未形成滑移通道。