拋石防波堤地震反應的顆粒流數(shù)值模擬

李存柱， 盛 儉， 張鴻斌， 曹蘇茜

(1. 江蘇煤炭地質(zhì)物測隊， 江蘇 南京 210046； 2. 南京工業(yè)大學交通運輸工程學院， 江蘇 南京 210009； 3. 南京水務集團有限公司， 江蘇 南京 21000)

0 引言

近年來，國內(nèi)外人工島工程屢見不鮮，著名的有迪拜棕櫚群島、日本大阪灣人工島和伊勢灣人工島、韓國首爾人工島及中國永暑礁西南島、美濟島等。歷史上關于拋石防波堤破壞的實例屢有報道， 1978年12月和1979年2月葡萄牙悉尼斯港拋石防波堤受風暴引起的波浪影響，其塊石護堤發(fā)生嚴重破壞[1]， 1964年美國阿拉斯加發(fā)生地震，致使蘇爾德和科迪亞克的防波堤嚴重破壞[2]，拋石防波堤破壞的后果不堪設想，因此保持拋石防波堤穩(wěn)定的重要性不言而喻。影響拋石防波堤穩(wěn)定性的主要因素有波浪和地震作用。目前國內(nèi)外對拋石防波堤穩(wěn)定性的研究主要基于波浪荷載的影響，對于地震作用下拋石防波堤穩(wěn)定性的研究相對較少。世界上很多拋石防波堤建造在地震高發(fā)地帶，并且已有很多拋石防波堤在破壞性地震作用下遭受破壞的工程實例[2-3]。拋石防波堤作為人工島的主要防護工程，可以防止波浪沖刷，抵擋海浪潮汐，有效地保護人工島基礎設施，減少經(jīng)濟損失和人員傷亡。

目前已有少量關于拋石防波堤在地震動作用下穩(wěn)定性的研究[3-5]，而部分研究中對于單一級配拋石防波堤地震動力響應的研究相對更少，且地震動輸入方式和顆粒粒徑對拋石防波堤動力響應的影響尚不明確。此外，拋石防波堤在地震動作用下的變形破壞機制尚不清楚。因此，研究單一級配堆石料的力學性質(zhì)和變形特征以及拋石防波堤在地震動作用下的動力響應，可以作為工程設計的依據(jù)，大大提高設計的準確性，減少因拋石工程失穩(wěn)引起的人員傷亡及財產(chǎn)損失，且有助于彌補科研空缺，為今后的相關研究提供一定的參考依據(jù)。

圖 1 魯?shù)榈卣鸺铀俣葧r程Fig.1 The acceleration time-history of Ludian earthquake

圖 2 魯?shù)榈卣鹚俣葧r程Fig.2 The velocity time-history of Ludian earthquake

顆粒流離散元被廣泛應用于巖石材料的相關研究中[6-8]。堆石料為不連續(xù)且各向異性的材料，顆粒流離散元剛好可以彌補大型振動臺試驗和有限元存在的不足，以便探明拋石防波堤在地震動作用下的動力響應以及變形破壞機制?；赑FC2D軟件平臺，對建造在剛性地基上的拋石防波堤進行數(shù)值模擬，魯?shù)榈貐^(qū)場地為基巖類似剛性地基，所以選取魯?shù)榈卣鹚涗浀募铀俣葧r程作為輸入地震動，分別分析了拋石防波堤的動力響應受輸入加速度峰值、不同輸入方式和顆粒粒徑大小的影響，以及其變形破壞機制。

1 數(shù)值模擬工況與建模

1.1 數(shù)值模擬工況

目前拋石防波堤在地震動作用下動力響應的研究主要集中于不同輸入峰值加速度對其動力響應的影響[3， 5， 9]，而關于地震動輸入方式、單級配堆石料粒徑大小對拋石防波堤在地震動作用下動力響應的影響研究較少。本文依據(jù)3種不同粒徑的顆粒共設計了3組模型，結合輸入地震動峰值加速度大小及加載方式確定了36個數(shù)值模擬工況。第一組、第二組、第三組數(shù)值模型中顆粒平均粒徑分別為1.2m、1.6m、2.0m，每一組地震動輸入方式為單一水平、單一豎向、水平和豎向組合輸入，單一水平輸入的峰值加速度分別為0.05g、0.1g、0.2g、0.4g，單一豎向輸入的峰值加速度分別為0.035g、0.07g、0.14g、0.28g，水平和豎向組合輸入的峰值加速度分別為0.05g和0.035g、0.1g和0.07g、0.2g和0.14g、0.4g和0.28g。通過上述不同工況拋石防波堤動力響應的顆粒流數(shù)值模擬，分析不同顆粒粒徑、不同輸入方式、不同輸入峰值加速度對拋石防波堤在地震動作用下動力響應的影響，以及拋石防波堤的變形破壞機制。

1.2 輸入地震動

數(shù)值模擬中選取魯?shù)榈卣鹬旋堫^山鎮(zhèn)強震動臺站記錄的南北向和豎向加速度時程作為水平向和豎向輸入地震動，其加速度峰值分別為0.73g和0.52g，圖 1為原始的加速度時程曲線。為了考慮不同輸入地震動水平的影響，通過調(diào)幅，將水平向輸入地震動峰值加速度分別調(diào)為0.05g、0.1g、0.2g、0.4g，豎向調(diào)幅與水平向等比縮小，調(diào)幅后的峰值加速度分別為0.035g、0.07g、0.14g、0.28g。

由于PFC的地震動輸入只能為速度時程，故需將原始加速度時程經(jīng)基線校正并一次積分得到速度時程，如圖 2所示。

1.3 模型構建

圖 3 模型幾何示意圖Fig.3 The geometric diagram of the model

圖 4 初始模型Fig.4 The initial model

參照文獻[6]中的幾何尺寸，構建了簡化且?guī)缀蜗嗨频膾伿啦ǖ棠Ｐ?。依?jù)拋石防波堤模型利用CAD中的多線段命令構建了相應的幾何圖形，如圖 3所示，并通過PFC2D的相關命令將幾何模型導入PFC2D中。

表 1 模型細觀接觸參數(shù)Tab.1 The meso-contact parameters of the model

以第一組數(shù)值模擬工況為例，顆粒粒徑取1.0～1.4m，平均粒徑1.2m，按照均勻分布的方式在圖 3所示的紅色矩形框中生成形狀不規(guī)則的顆粒簇，并賦予線性接觸模型參數(shù)，且設置法向和切向剛度以及摩擦系數(shù)，顆粒按照設定的孔隙率布滿紅色矩形區(qū)域，在此基礎上通過時步迭代進行應力釋放，并使之達到不平衡力極小的狀態(tài)。對模型顆粒之間及顆粒與墻體之間賦予平行黏結模型參數(shù)，在重力作用下使模型重新達到平衡狀態(tài)。最后進行削坡，并按照如圖 3所示的①→②、③→④、⑤→⑥、⑦的順序進行，每次削坡之后進行一次平衡。圖 4給出了最終的分析模型示意圖，模型頂寬20m，高30m，底寬140m，兩側坡面坡度比 1︰2。區(qū)域①中的顆粒在刪除前充當預壓的作用，使之達到相對密實的狀態(tài)。每組試驗模型中分別用不同粒徑的顆粒模擬實際拋石防波堤的堆石料，顆粒形狀仍然采用顆粒流雙軸試驗數(shù)值模擬中的不規(guī)則顆粒。在模型頂部中心位置處設置監(jiān)測點A用來監(jiān)測堤頂?shù)呢Q向位移和速度。

1.4 參數(shù)選取

顆粒流數(shù)值模型中顆粒之間的接觸參數(shù)進行標定，即通過數(shù)值模擬三軸試驗和室內(nèi)試驗的結果對比分析，具體顆粒間接觸參數(shù)詳見表 1，表中pb_kn、pb_ks分別為顆粒法向剛度和切向剛度，pb_ten、pb_coh分別為法向黏結強度和切向黏結強度，pb_r為黏結半徑，n為試樣孔隙率，ρ為顆粒密度，f為摩擦系數(shù)。

2 數(shù)值結果分析

在圖 1、2所示經(jīng)調(diào)幅的輸入地震動作用下，利用PFC2D分析得到了拋石防波堤地震反應，包括位移和加速度反應，并據(jù)此對拋石防波堤在輸入地震動作用下的變形和加速度響應特征進行分析。

2.3.2 穩(wěn)定性試驗 取“2.2.2”項下供試品溶液（批號：20170508）適量，分別于室溫下放置0、2、4、8、12、24 h時按“2.1”項下色譜條件進樣測定，以橙皮苷峰的保留時間和峰面積為參照，記錄各共有峰的相對保留時間和相對峰面積。結果，11個共有峰相對保留時間的RSD均小于0.50%，相對峰面積的RSD均小于0.79%（n＝6），表明供試品溶液于室溫下放置24 h內(nèi)基本穩(wěn)定。

圖 5 不同工況下坡頂沉降變化趨勢Fig.5 Variation of slope top settlement under different working conditions

2.1 拋石防波堤坡頂沉降

(1)不同輸入條件下坡頂沉降分析

圖 5給出了坡頂中心A處的豎向位移隨時間的變化規(guī)律。以圖 5a為例，發(fā)現(xiàn)坡頂沉降受地震動輸入方式影響較大，坡頂沉降量從大到小的順序為：水平和豎向組合輸入最大、單一水平輸入次之、單一豎向輸入最小。相關研究通常只對水平輸入地震動下的拋石防波堤進行研究[3， 4，5， 9]，由圖 5可以看出豎向輸入對拋石防波堤變形的影響不可忽略。在顆粒流數(shù)值模擬中設定終止條件為模型達到穩(wěn)定狀態(tài)，由圖 5可以看出，單一豎向加載時模型達到平衡狀態(tài)所用時間最短，單一水平加載次之，水平和豎向組合加載最長。由此可知，在輸入地震動作用下，拋石防波堤的穩(wěn)定性受水平和豎向荷載組合作用時的影響最大，水平加載次之，豎向加載影響最小。

由圖 5可以看出，3種加載方式最初的沉降速率比較一致，隨后沉降速率逐漸減小。由輸入的加速度時程曲線得知，隨著時間的推移加速度幅值越來越小，大約20s后加速度很小。從定性角度分析，沉降速率逐漸減小的原因主要有：模型自身的內(nèi)部結構隨著地震動的輸入在變化，加速度幅值在不斷地減小。原始輸入的加速度持時80s，但是所有模型記錄的位移時程持時均未達80s，這表明加速度時程中加速度幅值極小時對拋石防波堤的穩(wěn)定性幾乎沒有影響。

(2)峰值加速度對坡頂沉降的影響

依據(jù)不同工況下的拋石防波堤坡頂?shù)淖罱K沉降量，見表 2，分析了輸入峰值加速度對坡頂沉降的影響，結果表明：所有工況在輸入方式相同、顆粒粒徑相同條件下，輸入加速度峰值越大最終沉降量越大，所以拋石防波堤的穩(wěn)定性受輸入峰值加速度影響較大。

圖 6給出了拋石防波堤在地震動輸入方式相同時，坡頂最終沉降隨輸入加速度峰值的變化規(guī)律。由圖 6可以看出輸入條件相同時，峰值加速度越大，拋石防波堤最終的沉降量越大。如圖6(a)和(c)所示，單一水平輸入和水平與豎向組合輸入時，輸入峰值加速度相同，顆粒粒徑越大最終沉降量越小。

表 2 不同工況下坡頂最終沉降量統(tǒng)計結果Tab.2 The statistical results of final settlement of slope top under different working conditions

圖 6 拋石防波堤坡頂沉降隨地震動峰值加速度的變化趨勢Fig.6 Variation of slope top settlement of rubble-mound breakwater with peak acceleration of ground motion

較大粒徑顆粒組成的堆石料顆粒之間的咬合作用大于小粒徑顆粒組成的堆石料，表現(xiàn)為大粒徑堆石料整體抗剪強度大于小粒徑堆石料，由此可以看出相同輸入地震動下，顆粒流數(shù)值模擬得到的剪應力引起的大粒徑顆粒模型變形小于小粒徑顆粒模型的變化規(guī)律是合理的。然而，豎向輸入加速度峰值小于0.1g時，顆粒粒徑越小坡頂?shù)淖罱K沉降量越大； 當豎向峰值加速度0.2g時，平均粒徑為120cm的最終沉降量最大，平均粒徑200cm的次之，平均粒徑 160cm的最??； 豎向峰值加速度0.4g時，平均粒徑的最終沉降量排序為200cm>120cm>160cm。由此可知，單一豎向加載最終沉降量變化規(guī)律不同于單一水平加載、水平與豎向組合加載，說明了最終沉降量不僅與輸入峰值加速度有關，還與加載方式有關。單一豎向輸入時顆粒應力以正應力為主，單一水平輸入時顆粒應力以剪應力為主，水平和豎向同時輸入時內(nèi)部應力比較復雜，所以豎向作用引起的拋石防波堤沉降規(guī)律與其它兩種作用方式不同。

圖 7 拋石防波堤典型的破壞模式[6]Fig.7 Typical failure modes for rubble-mound breakwaters[6]

圖 8 工況12拋石防波堤初始模型與輸入地震動作用下變形模型Fig.8 Theinitial model of rubble-mound breakwater and the deformation model under input ground motion under working condition 12

2.2 拋石防波堤破壞模式

拋石防波堤在地震動作用下有3種典型的破壞模式[5]，如圖 7所示。圖 7(a)為拋石防波堤在地震作用下上部結構發(fā)生不均勻沉降； 圖7(b)為地震作用下，拋石防波堤地基沉降或地基液化引起的坡頂沉降及橫向變形或上部結構不均勻沉降； 圖7(c)為地震引起地基液化，從而導致結構失穩(wěn)或上部沉降及傾覆。顆粒流數(shù)值模擬結果表明，相對幾何尺寸較大的模型，在地震動作用下模型的變形較小，所以在分析拋石防波堤破壞模式時選取所有工況中變形最大的一組，即平均顆粒粒徑為120cm，地震動為水平和豎向組合輸入下，峰值加速度分別為0.4g和0.28g的工況。該工況下拋石防波堤初始模型與輸入地震動作用下變形模型如圖 8所示，由圖可以看出初始模型邊界顆粒緊貼模型幾何構造線邊緣，在加載達到穩(wěn)定狀態(tài)以后，上部兩側坡體及坡頂均與原始模型構造線有所偏離。從模型的破壞變形模式可以看出與圖 7(a)所示的變形規(guī)律一致，由此驗證了利用顆粒流離散元進行拋石防波堤動力響應數(shù)值分析的可行性。

圖 9 工況12拋石防波堤位移矢量圖(單位：m)Fig.9 Displacement vector diagram of rubble-mound breakwater under working conditions 12(unit： m)

表 3 坡頂峰值加速度放大系數(shù)統(tǒng)計結果Tab.3 Statistical result of peak acceleration amplification factor at slope top

圖 10 監(jiān)測點A峰值加速度放大系數(shù)隨輸入地震動峰值加速度變化趨勢Fig.10 Variation of acceleration amplification coefficient with peak acceleration of input strong motion at monitoring point A

為便于更清晰地了解拋石防波堤在地震動作用下的變形機制，將以上工況對應的位移矢量圖進行深入分析。由圖 9可以看出，區(qū)域①主要以豎向沉降為主，兩側坡體②和③分別斜向下滑動，坡面表層滑動，其主要原因是坡體內(nèi)部受側向約束力較大，因此區(qū)域①變形以豎向變形為主，而兩側坡體面向凌空面，坡面無側向約束，故而受力主要為豎向正應力和水平向的側壓力，由此引起的區(qū)域②和③既有豎向位移也有水平位移，即總位移方向斜向下。此外，盡管坡體幾何對稱，但內(nèi)部堆石顆粒結構不對稱，其材料存在各向異性，從而導致左右坡體變形不完全對稱，且左側坡面的淺層滑動略大于右側坡面。同時，拋石防波堤顆粒流數(shù)值模擬顯示兩側坡腳處顆粒有明顯向外運動的趨勢，因此在拋石防波堤工程設計中宜考慮坡腳的這種變形，采取相關措施加以控制。

2.3 加速度響應特征

依據(jù)計算峰值加速度可統(tǒng)計得到水平輸入地震動下各工況坡頂峰值加速度放大系數(shù)，統(tǒng)計結果列于表 3，峰值加速度放大系數(shù)隨輸入地震動峰值加速度的變化關系如圖 10所示。由圖 10和表 3可以看出，同一個模型，顆粒粒徑相同時，峰值加速度放大系數(shù)隨輸入地震動峰值加速度增大而逐漸減小，且峰值加速度小于等于0.10g時放大系數(shù)衰減較快，隨著輸入峰值加速度增加放大系數(shù)衰減逐漸變??； 相同輸入峰值加速度下，顆粒粒徑越大加速度放大系數(shù)越大。由此可得出，輸入峰值加速度越大，拋石防波堤發(fā)生的變形越大，耗散的地震動能量越多； 相同輸入峰值加速度下，大粒徑模型發(fā)生的變形較小，小粒徑模型發(fā)生的變形大，同樣由于小粒徑的顆粒消耗的能量更多，所以粒徑越小峰值加速度放大系數(shù)越小。

3 結論

基于顆粒流離散元方法分析了拋石防波堤在輸入地震動作用下的動力響應，在此基礎上研究了不同粒徑顆粒構成的拋石防波堤在地震動作用下的變形特征。輸入地震動引用魯?shù)榈卣鹬旋堫^山強震動臺站記錄到的南北向和豎向的加速度時程，離散元數(shù)值模擬中輸入方式分別為單一水平向輸入、單一豎向輸入及水平和豎向組合輸入。經(jīng)過對離散元數(shù)值模擬結果進行分析，得出的主要結論如下：

(1)當顆粒粒徑相同時，輸入峰值加速度相同，發(fā)現(xiàn)拋石防波堤在水平和豎向組合輸入時坡頂沉降最大，單一水平向輸入次之，單一豎向輸入最??； 當顆粒粒徑相同、輸入方式也相同時，拋石防波堤坡頂沉降隨輸入峰值加速度的增大而增大； 單一豎向輸入防波堤坡體達到穩(wěn)定狀態(tài)用時最短，水平和豎向組合輸入時用時最長，單一水平向輸入時用時介于二者之間。

(2)當顆粒粒徑相同時，拋石防波堤坡頂峰值加速度放大系數(shù)隨輸入峰值加速度的增加逐漸減小，且輸入峰值加速度越大衰減越慢； 相同輸入峰值加速度時，拋石防波堤坡頂?shù)姆逯导铀俣确糯笙禂?shù)受粒徑影響較大，顆粒粒徑越大峰值加速放大系數(shù)越大，反之，顆粒粒徑越小峰值加速度放大系數(shù)越小。

(3)本研究所得到的拋石防波堤變形特征與其在地震動作用下的經(jīng)典變形特征相符，表明通過離散元數(shù)值模擬拋石防波堤在地震動作用下的動力響應是可行的。拋石防波堤在地震動作用下的變形特征表明，剛性地基上的拋石防波堤受地震動作用主要表現(xiàn)為坡頂沉降，兩側坡體斜向下滑動，且表面存在淺層滑動變形，由于內(nèi)部結構具有不對稱性和各向異性，所以拋石防波堤變形不完全對稱。兩側坡腳有向外運動的趨勢或已發(fā)生向外變形，所以有必要在拋石防波堤工程設計和施工中采取相關措施，以防止坡腳在地震動作用下發(fā)生破壞。