顆粒離散元模擬邊坡動力響應(yīng)的研究進展

常文斌，王 平，2，于一帆，王會娟

(1.中國地震局蘭州地震研究所 黃土地震工程重點實驗室， 甘肅 蘭州 730070；2.西安理工大學 巖土工程研究所， 陜西 西安 710048)

邊坡是一個復雜多變的系統(tǒng)，其穩(wěn)定性受多種條件的影響，如自然條件、邊坡結(jié)構(gòu)、人類活動等，這些影響因素大多是隨機的、可變的，其不確定性增加了邊坡穩(wěn)定性的研究難度。

邊坡在動荷載作用下的穩(wěn)定性分析較為復雜，它包含了工程地震學和巖體力學等多學科的交叉以及諸多技術(shù)問題。動荷載包括地震荷載、風荷載等自然荷載及機械荷載、爆炸荷載等人為荷載。與其他動荷載不同的是，地震荷載并非一種直接的力學作用，而是地面發(fā)生位移時坡體結(jié)構(gòu)在慣性作用下產(chǎn)生了一定的運動差，導致邊坡產(chǎn)生嚴重的變形。同時，地震荷載的頻譜特性及傳播方向?qū)ζ麦w的破壞均有很大的影響，加之需要考慮巖土體的動本構(gòu)關(guān)系及動強度，導致邊坡動力穩(wěn)定性問題變得更加復雜。然而隨著計算機科學技術(shù)的發(fā)展，數(shù)值模擬技術(shù)已成為邊坡動力分析領(lǐng)域不可或缺的手段之一。當下主要的數(shù)值方法包括有限差分法、有限元法、離散元法[1]及邊界元法等。現(xiàn)階段，有限差分法及有限元法是邊坡動力分析常用的兩種方法，由于這兩種方法將坡體視為連續(xù)介質(zhì)，因而在動力分析方面有著獨特的優(yōu)勢。但將巖土體視為連續(xù)介質(zhì)則不能真實的體現(xiàn)出介質(zhì)材料的流動變形特征，造成了有限差分法及有限元法在模擬邊坡變形時的局限性。

1971年美國學者 Cundall[2]提出了離散單元法。隨后Cundall等[3-4]首次把顆粒流法定為一種特殊的離散單元法。在此基礎(chǔ)上，美國Itasca 公司開發(fā)出了二維顆粒流軟件PFC2D(Particle Flow Code in 2Dimensions, PFC2D)[5]和三維顆粒流軟件PFC3D(Particle Flow Code in 3Dimensions, PFC3D)。PFC中，球和簇的運動遵守牛頓第二運動定律，墻體的運動通過用戶來指定。顆粒之間的接觸可以選用不同相互作用定律的接觸模型來施加，在受力體系下，通過顆粒間內(nèi)力和位移的不斷更新來模擬介質(zhì)體的變形[6]。顆粒流方法很好的解決了常規(guī)數(shù)值模擬方法在分析大變形巖土工程破壞時的局限性，能夠更好的揭示微觀及細觀介質(zhì)材料的破壞機理，因此備受廣大學者青睞。

顆粒流離散元法解決了巖土工程大變形及介質(zhì)的微觀破壞問題，可更方便的處理非連續(xù)介質(zhì)力學問題，但是對于邊坡的動力響應(yīng)研究由于啟動較晚，目前仍有諸多問題需要解決。本文總結(jié)歸納了顆粒流進行邊坡動力響應(yīng)研究時，阻尼參數(shù)標定、動荷載輸入方式和動力分析參量這三個方面的研究現(xiàn)狀及存在的問題。分析過程中基于地震動力作用的復雜及特殊性，著重對離散元顆粒流方法在邊坡地震動力響應(yīng)方面存在的關(guān)鍵問題進行了分析和探討。

1 顆粒流中阻尼參數(shù)標定

1.1 PFC程序能量耗散機制

PFC程序中，單元之間的能量是通過摩擦滑動來進行耗散的。由于單元之間的滑動摩擦并非主動發(fā)生，因此在所需的計算時步內(nèi)，僅靠摩擦滑動系統(tǒng)很難達到一個穩(wěn)定的狀態(tài)。鑒于此，PFC中提供了3種阻尼耗能機制：局部阻尼、組合阻尼、黏性阻尼。其中，局部阻尼通過對每一個顆粒單元施加一個與其速度方向相反的阻尼力。通過局部阻尼理論推導部分可知，施加于顆粒上的阻尼力于顆粒的不平衡力成正比，兩者之間的比值即為局部阻尼系數(shù)。PFC中局部阻尼力的施加則通過設(shè)置局部阻尼系數(shù)來實現(xiàn)，默認值為0.7。組合阻尼同時考慮了顆粒的速度與不平衡力，但其能量耗散率較低，一般不被考慮使用。黏性阻尼施加于PFC模型中每一個接觸上，包括法向阻尼和切向阻尼。黏性阻尼的方向與接觸處的相對速度方向相反，其定義通過指定黏性阻尼比參數(shù)來實現(xiàn)。

許多學者對PFC中阻尼參數(shù)的選取做了相關(guān)研究[7-9]，結(jié)果表明：在靜力情況下，PFC中采用局部阻尼以吸收單元之間的動能，其阻尼取值僅影響系統(tǒng)達到穩(wěn)定解的時間。而在動力問題中，應(yīng)力波在巖土體的傳播過程中，阻尼取值直接影響著巖土體的振動，從而直接影響結(jié)算結(jié)果，因此必須確保阻尼取值的科學性。當下PFC在巖土工程中的模擬研究多集中于顆粒細觀參數(shù)等對模擬結(jié)果的影響，阻尼系數(shù)取值的重要性被多數(shù)學者忽略，無形中影響了模擬結(jié)果的準確性。

1.2 顆粒阻尼參數(shù)的選取

顆粒流模擬過程中，阻尼參數(shù)的選取往往存在主觀經(jīng)驗性及隨意性。為解決這一問題，鐘文鎮(zhèn)等[7]采用聲頻采樣技術(shù)對鋼球自由落體物理模型中的阻尼系數(shù)進行了標定，得出的阻尼系數(shù)(0.5左右)與實際較為相符。同時，通過顆粒在容器中隨機堆積的模擬試驗驗證了所標定系數(shù)的可靠性。楊冰[8]在PFC2D中采取不同的局部阻尼系數(shù)及黏性阻尼系數(shù)進行了顆粒的自由落體試驗，計算不同阻尼系數(shù)所對應(yīng)的能量衰減率，擬合出了顆粒自由落體時能量衰減率與阻尼系數(shù)的關(guān)系，驗證了顆粒材料通過設(shè)置相應(yīng)的阻尼系數(shù)也可用于工程的模擬。周國慶等[9]通過鋼球自由下落和懸臂梁的剪彎受力模擬試驗，明確了顆粒流中局部阻尼系數(shù)和黏性阻尼系數(shù)的適用范圍。結(jié)果表明靜態(tài)問題可采用局部阻尼減少計算時間，動態(tài)問題需根據(jù)不同情況選取局部阻尼或黏性阻尼，并通過鋼球的自由落體試驗確定了鋼球的法向黏性阻尼比為0.16。

由此可見，PFC中阻尼參數(shù)的選取對動力問題計算結(jié)果影響較大，在進行邊坡動力分析模擬工作時更需要注意。以地震作用為例，地震作用下，顆粒集合體邊坡發(fā)生滑動后，多數(shù)顆粒為單向運動，此時阻尼類型的選取尤為關(guān)鍵。如果對顆粒設(shè)置局部阻尼，由于局部阻尼是對顆粒速度相反方向施加一個阻尼力，只要顆粒存在速度就會產(chǎn)生作用。這會造成此后計算的每一時步都會產(chǎn)生與顆粒運動相反的阻尼力，導致顆粒的運動速度急劇減小，能量耗散迅速，從而影響最終的模擬結(jié)果。因此局部阻尼不適合地震作用下的顆粒流模型，部分學者在研究過程中很容易忽視這一點，對地震邊坡模型采用局部阻尼且系數(shù)過大。相比之下，黏性阻尼比較適合地震邊坡模型，它相當于在顆粒間每一個接觸點上施加法向和切向阻尼器，阻尼力由系統(tǒng)的有限質(zhì)量及自振頻率等決定，并在計算循環(huán)時添加到接觸力上。因此在顆粒流邊坡動力分析過程中，根據(jù)不同的動力荷載及研究問題，選取合適的阻尼類型是確保模擬計算成功的前提，值得重視。

2 顆粒流中動荷載的輸入方式

在離散元顆粒流中，巖土體介質(zhì)材料通過顆粒及其之間的接觸特性來表示。應(yīng)力波在彈性狀態(tài)的連續(xù)性介質(zhì)中的傳播過程，同樣可以在顆粒流中進行模擬。

2.1 動邊界條件

在顆粒流方法中，應(yīng)力波可以通過對顆?；蛘邏w施加來模擬動荷載。無論應(yīng)力波施加于顆粒還是墻體，均要提及一個問題：應(yīng)力波傳播的動邊界條件。

對于數(shù)值模型一般采用的單面邊坡，在動力問題中，模型邊界條件對動力分析結(jié)果有較大影響。分析模型的范圍越大，分析結(jié)果就越貼近事實，但大的模型由于顆粒數(shù)目較多會導致計算量變大，這時候需要設(shè)置一定的人工邊界。全局人工邊界與局部人工邊界是人工邊界中最常用的兩類邊界。其中全局人工邊界是一種較為精確的人工邊界，其邊界參數(shù)與相鄰介質(zhì)的特性相聯(lián)系，藕聯(lián)性較強，對于顆粒離散元這種沒有控制方程的計算方法并不適用。局部人工邊界主要有黏性邊界[10]、黏彈性邊界[11]、吸收層邊界[12]與透射邊界[13]。黏性邊界、黏彈性邊界以及透射邊界在有限元及有限差分法中有大量的應(yīng)用，以往用離散元顆粒流進行邊界處理時大多參考這些數(shù)值方法的邊界處理方法[14-16]。

顆粒流中動邊界存在多種形式，其主要目的都在于吸收掉入射波的動能，防止其在邊界處發(fā)生反射，從而達到模擬無限介質(zhì)的效果。吸收層邊界通過在模型外施加一層具有耗能作用的緩沖區(qū)吸收透射波(對邊界吸收層顆粒施加一定的阻尼系數(shù))，該吸收區(qū)又叫做阻尼層[12]。該邊界最大的優(yōu)點就是邊界的幾何形狀可以任意設(shè)置，適用于任意入射角、各種頻率的入射波。但是這種對模型邊界顆粒施加阻尼系數(shù)的高阻尼吸波方法并非完美。石崇等[17-18]在對顆粒流中應(yīng)力波的傳播進行研究時，對黏滯透射邊界和高阻尼吸收邊界兩種邊界的合理性進行了驗證。結(jié)果表明應(yīng)力波傳播過程中高阻尼邊界容易導致邊界附近的應(yīng)力峰值明顯降低，有悖于自然界中應(yīng)力波的傳播規(guī)律。而黏滯透射邊界情況下雖然應(yīng)力波傳播距離略有延遲，但應(yīng)力波的振動規(guī)律較為接近。周興濤等[19-20]在二維顆粒流軟件中分別建立了黏性、黏彈性及自由邊界等人工動力模型，探討了邊界顆粒分布模式對邊界作用的影響，驗證了黏彈性邊界的正確性。當下顆粒流邊坡動力分析中，黏滯透射邊界多被學者使用，其模擬結(jié)果也與振動臺及離心機試驗結(jié)果較為相符。

當下連續(xù)性介質(zhì)方法與非連續(xù)介質(zhì)方法耦合計算方法已被廣泛使用，包括離散元與有限差分法[21-23]，離散元與有限元等[24-28]。不同數(shù)值方法的耦合取長補短，以PFC與FLAC耦合為例，李祥龍[29-30]采用顆粒離散元來模擬邊坡的潛在破壞區(qū)域，用FLAC模擬邊坡其他部位，進行了邊坡的地震動力響應(yīng)模擬。這種方法不僅解決了顆粒流模擬大型邊坡計算量過大的問題，同時避免了離散元中動邊界條件設(shè)置不成熟的問題。該方法對于離散元與其他數(shù)值方法耦合模擬邊坡動力響應(yīng)的研究提供了新的思路。

2.2 應(yīng)力波的作用方式

應(yīng)力波的作用方式包括應(yīng)力波的作用方向(不同入射方向)和作用性質(zhì),同時應(yīng)力波的頻率強度、持時及類型對數(shù)值模擬結(jié)果均存在一定的影響。通常研究中施加的爆破荷載[31-33]、機械荷載[34-35]等為動力荷載在模擬時應(yīng)力波的作用方式相對明確，影響較小。但對于類似地震波的不可控荷載，其作用方式則存在較多的爭議。

目前，在邊坡地震顆粒流值模擬分析中，地震動普遍采用一致性輸入。實際中由于遠近場的作用，地震動在經(jīng)過不同路徑傳播到邊坡底部時，邊坡底部不同位置接收到的地震動是存在一定差異的，所以當下一致性輸入地震動的方法是不完全正確的。由于目前邊坡地震動力響應(yīng)數(shù)值分析仍集中于有限元及有限差分法等方法，地震動作用方式的相關(guān)研究也都是基于這些方法展開。許多學者基于地震波的不同入射方向及不同震相的性質(zhì)做了較多的研究[36-37]。國內(nèi)以孫進忠教授為代表的團隊在這方面做了相當多的研究工作[38-40]。結(jié)果表明不同的地震波作用方式對于邊坡的動力響應(yīng)的確存在著一定的影響。這些研究成果對于在顆粒流方法中進行應(yīng)力波輸入時有著很大的啟發(fā)。

即便如此，采用地震臺站記錄的典型強震地震動數(shù)據(jù)對模型直接輸入或進行簡單調(diào)整后作為地震動輸入(顆粒流中對顆?；驂w施加速度)，雖然輸入的是實際地震動，但由于巖土體參數(shù)、傳播路徑、震源參數(shù)等的差異，所得分析結(jié)果也只能進行邊坡動力破壞機理的研究，并不能反映實際的動力響應(yīng)情況[41]。因此，顆粒流中應(yīng)力波(尤其地震波)的作用方式這是一個相當復雜問題，如何利用震源和巖土介質(zhì)特性來確定適合于數(shù)值邊坡動力分析的地震動類型，其入射方向、幅值、頻譜和持時該如何確定，仍有待于更深層次的研究。

3 顆粒流中邊坡動力分析方式

3.1 顆粒流中邊坡動力分析對象

動荷載作用下，邊坡動力響應(yīng)分析的對象主要是速度、位移、加速度、應(yīng)力等。在顆粒離散元邊坡動力分析中，現(xiàn)有的文獻大多都是對振動臺模型試驗、有限元、有限差分法等研究方法與研究結(jié)果的重復，多數(shù)學者關(guān)注點仍局限于邊坡不同部位的速度、位移及加速度變化情況[42-46]，在此不再贅述。誠然，這些參量可以在一定程度上體現(xiàn)邊坡的動力響應(yīng)，但沒有充分體現(xiàn)顆粒離散元的優(yōu)越性。例如顆粒流中，介質(zhì)的破裂能夠定義巖土體損傷、力鏈可以用來考察能量傳輸?shù)?。顆粒流邊坡模擬的過程中，可以對顆粒間的不同類型的接觸力變化及坡體內(nèi)裂紋的分布情況進行監(jiān)測[47-48]，坡體中不同類型的裂隙(拉裂隙、剪裂隙)展布及裂隙數(shù)目可以清楚的展示出來。Wu等[49]采用PFC2D對含裂縫黃土邊坡的地震作用下破壞過程進行模擬，邊坡破壞過程中裂紋分布與振動臺模擬結(jié)果類似。同時，顆粒流中每個計算步都會累積計算球、墻和接觸的能量，諸如動能、摩擦能、應(yīng)變能等能量；測量圓方法可以對小范圍內(nèi)顆粒的應(yīng)力、孔隙率、應(yīng)變率等進行監(jiān)測。丁建輝[50]和蘇燕等[51]采用二維顆粒流軟件對地震降雨型滑坡形成機理進行了研究，通過對地震及降雨作用下坡體內(nèi)的速度場、位移場、應(yīng)力場、孔隙率、顆粒接觸數(shù)及滑動比進行分析，揭示了地震降雨滑坡的形成機制。當下，將顆粒流的這些細觀分析優(yōu)勢與邊坡動力響應(yīng)結(jié)合起來，而不是僅僅局限于速度、位移等傳統(tǒng)的參數(shù)分析，有助于從不同角度分析邊坡的破壞機理。

3.2 顆粒流中邊坡加速度響應(yīng)相關(guān)問題

在顆粒流中進行邊坡動力響應(yīng)分析時，由于不能直接監(jiān)測坡體內(nèi)測點處顆粒的加速度，目前多數(shù)學者采用的做法是將記錄到的速度時程曲線進行微分，對得到的加速度時程曲線進行分析[15-16]。由于顆粒離散元的特殊性，在動荷載作用的過程中，邊坡某一監(jiān)測點的顆粒在小范圍內(nèi)處于一個不斷運動的狀態(tài)，這種運動狀態(tài)受到其周圍顆粒的碰撞和應(yīng)力波的共同影響。因此在以往的研究中，通過微分速度時程曲線得出顆粒的加速度時程曲線與原始應(yīng)力波曲線存在較大差異的現(xiàn)象多有發(fā)生。那么在顆粒流中基于這種方法所得加速度動力響應(yīng)的可靠性就值得商榷。如何有效地解決或避免當下加速度響應(yīng)分析存在偏差的這一問題，是顆粒流邊坡加速度響應(yīng)分析工作的關(guān)鍵。劉一飛[12]采用Fish語言編程實現(xiàn)對顆粒加速度的定義，即遍歷顆粒的在某一瞬間的不平衡力，通過不平衡力除以顆粒質(zhì)量得到顆粒的加速度，但其研究結(jié)果表明通過不平衡為計算的加速度分布有很大的隨機性和離散性，其峰值加速度(PGA)并不準確，參考意義并不是很大。筆者建議嘗試在顆粒流建模過程中，首先盡量設(shè)置顆粒的粒徑較小，以減小顆粒在小范圍內(nèi)的運動。其次，在保證某一監(jiān)測點處局部顆粒的總體動力響應(yīng)情況一致的前提下，通過對監(jiān)測點小范圍內(nèi)的顆粒的監(jiān)測數(shù)據(jù)求取平均值，以減少偶然性誤差。

在此基礎(chǔ)上，應(yīng)力波頻譜分析在動力響應(yīng)中的應(yīng)用較少，并不是很成熟。地震頻譜可以定量的揭示地震動的動力特性，以不同頻率分量的表現(xiàn)來表現(xiàn)地震動及其對結(jié)構(gòu)的作用。傅里葉頻譜[52]和反應(yīng)頻譜[53]是地震波頻譜特性分析時常用的頻譜，其中，傅里葉頻譜反應(yīng)了地震波在坡體內(nèi)部傳播過程中自身頻率各分量的變化情況，反應(yīng)頻譜表征了地震動對上部自由彈性結(jié)構(gòu)的影響特征。將頻譜特征響應(yīng)引入顆粒流邊坡動力分析中，結(jié)合顆粒流邊坡的宏觀破壞機制，分析結(jié)果對于邊坡的防護治理及坡體上部臨坡建筑物的抗震設(shè)計意義重大。然而，顆粒流中邊坡頻譜響應(yīng)特征的分析需要基于監(jiān)測點處加速度時程曲線提取方法的成熟，需要進一步的探索。

4 結(jié)論及展望

當前離散元方法已成為研究邊坡變形破壞的主要手段之一，開展在顆粒離散元基礎(chǔ)上的邊坡動力響應(yīng)研究，對于現(xiàn)有的邊坡動力分析研究成果是一個極大的補充。本文從阻尼參數(shù)標定、動荷載輸入方式與動力分析三個方面，對離散元顆粒流在邊坡動力響應(yīng)領(lǐng)域的研究進展進行了總結(jié)與分析。目前值得注意的關(guān)鍵問題如下：

(1)顆粒流邊坡動力破壞模擬過程中阻尼類型的選取及阻尼系數(shù)的標定應(yīng)引起重視。離散元模型中，組合阻尼由于能量耗散率較低，一般不建議使用；局部阻尼適用于靜力作用下吸收顆粒單元之間的動能減少系統(tǒng)達到穩(wěn)定的時間，而在動力作用下對能量的耗散較大，因此不適應(yīng)動力分析時使用；黏性阻尼的阻尼力在振動過程中是不斷循環(huán)添加至單元上的，因此適用于動力響應(yīng)分析時使用。因此在模擬過程中，應(yīng)根據(jù)不同情況選取對應(yīng)的阻尼類型及阻尼系數(shù)，以提高計算的準確性。

(2)如何結(jié)合巖土體材料特性及震源機制等對邊坡模型合理的輸入對應(yīng)的地震波。在模擬邊坡動力破壞的過程中，地震動參數(shù)的選取往往至關(guān)重要。不同三要素(幅值、頻譜、持時)的地震動數(shù)據(jù)對模型計算結(jié)果的影響不可忽視；其次地震動的遠近場輸入、斜入射、垂直入射等輸入方式對模擬結(jié)果也存在一定的影響。如何根據(jù)邊坡的實際場地條件及震源機制選取合理的地震動輸入方式仍需進一步的研究。

(3)離散元邊坡加速度響應(yīng)分析中顆粒加速度的獲取方式。顆粒離散元的離散特性使其在加速度響應(yīng)分析時無法直接導出顆粒的加速度值。目前兩種常用的方法均存在一定的誤差，不平衡力轉(zhuǎn)化法得到的結(jié)果隨機性及離散性較大，速度時程微分法得到的顆粒速度時程曲線受多個因素影響。因此更準確的獲取模型破壞過程中顆粒的加速度值的方法有待進一步研究。

最后，顆粒離散元模擬邊坡動力響應(yīng)的研究工作仍有很大的發(fā)展空間。首先，將離散元顆粒流可模擬巖土體大變形及細觀破壞的優(yōu)勢與動力分析相結(jié)合。例如離散元模擬過程中可檢測坡體內(nèi)裂紋數(shù)量、顆粒間接觸的破壞數(shù)量、接觸力強度的變化情況等。將顆粒流的這些細觀分析優(yōu)勢與邊坡動力響應(yīng)結(jié)合起來，而不是僅僅局限于傳統(tǒng)的參量的分析，有助于從不同角度分析邊坡的破壞機理。其次，顆粒離散元與其他數(shù)值方法耦合模擬邊坡動力響應(yīng)的相關(guān)研究已初露頭角，如顆粒離散元與有限元、有限差分法、光滑粒子流方法等的耦合可以使邊坡動力響應(yīng)的數(shù)值方法面向更多的研究，也為邊坡的穩(wěn)定性分析提供了新的途徑。