李同春，樊舒婕，劉曉青，趙蘭浩

（1.河海大學(xué)水利水電學(xué)院，江蘇省南京市 210098；2.水資源高效利用與工程安全國家工程研究中心，江蘇省南京市 210098）

0 引言

在水利工程中存在著許多局部不連續(xù)的現(xiàn)象，例如拱壩的橫縫、巖體的裂隙、地基的軟弱面等。這些局部不連續(xù)區(qū)域在地震等動力荷載作用下的張開、閉合和摩擦滑移關(guān)系著結(jié)構(gòu)穩(wěn)定與工程安全?？紤]動載下不連續(xù)區(qū)域交界面的動態(tài)接觸問題對于結(jié)構(gòu)動力分析有著十分重要的意義。對于此類問題，已有的數(shù)值分析方法包括傳統(tǒng)的有限元法FEM[1-3]、離散元法DEM[4，5]和非連續(xù)變形分析方法DDA[6-8]，以及FEM-BEM[9]、DEM-FEM[10-13]、DEM-BEM[14]等混合算法。

本文針對動力情況下非連續(xù)變形的接觸非線性問題，基于文獻(xiàn)[15]的力學(xué)模型提出了一種新的動力迭代求解方法。推導(dǎo)了包含接觸力、塊體形心剛體加速度和接觸結(jié)點(diǎn)總位移在內(nèi)的非線性方程，將該非線性方程與塊體形心動力平衡方程進(jìn)行耦合，建立了以界面接觸力和塊體形心剛體加速度為混合變量與界面結(jié)點(diǎn)總位移進(jìn)行迭代求解的動力界面迭代方程，給出了混合解法的實(shí)施步驟并通過典型算例驗(yàn)證了動力混合解法的正確性。

1 IBE-PFE動力混合算法的力學(xué)表達(dá)

IBE-PFE動力混合算法將求解系統(tǒng)分解為若干塊體與非連續(xù)界面，每個塊體同時作為系統(tǒng)有限元的一個分區(qū)，以塊體結(jié)點(diǎn)加速度作為分區(qū)有限元計(jì)算基本變量，以塊體剛體形心加速度作為剛體運(yùn)動變量，當(dāng)接觸界面上接觸點(diǎn)對的接觸內(nèi)力增量已知，則可通過動力有限元法求出塊體的加速度增量，進(jìn)而根據(jù)Newmark法[16]求解各結(jié)點(diǎn)的應(yīng)力和位移。建立接觸界面上力與點(diǎn)對總位移方程，在塊體形心處建立接觸內(nèi)力、塊體形心加速度增量引起的力的平衡方程，將以上方程進(jìn)行耦合，塊體形心接觸內(nèi)力、加速度作為混合迭代變量和接觸界面點(diǎn)對總位移的求解可以通過迭代方程求得。

2 IBE-PFE動力混合算法

2.1 動力分區(qū)有限元方程

其中，Ki為剛度矩陣；Di為阻尼矩陣，這里采用瑞利比例阻尼假定Di=αMi+βKi，α、β為瑞利比例阻尼系數(shù)；Mi為質(zhì)量矩陣，F(xiàn)i為作用于第i個塊體上的外荷載；Ri為由不連續(xù)界面上由內(nèi)力傳輸?shù)暮奢d向量，可由下式表示：

其中，Ti為該塊體內(nèi)結(jié)點(diǎn)上的作用力與接觸點(diǎn)對上的接觸力轉(zhuǎn)換系數(shù)矩陣，F(xiàn)cL為局部坐標(biāo)系下單元接觸點(diǎn)對之間的接觸力。

2.2 接觸塊體的剛性運(yùn)動及動力平衡方程

如果第i塊的形心位移被設(shè)為剛體位移變量，那么根據(jù)剛體運(yùn)動法則，可視為塊體內(nèi)的每個點(diǎn)的剛體位移是與此相關(guān)、且遵守剛體運(yùn)動法則的。在每一個接觸塊體形心點(diǎn)上分別建立塊體力的平衡方程，對于每個塊體，離散單元的動力平衡方程為：

其中，ωk為第k個結(jié)點(diǎn)相對于形心的轉(zhuǎn)換矩陣。mi為與第i塊塊體的質(zhì)量，Ii為與第i塊塊體相對于形心的轉(zhuǎn)動慣量為質(zhì)量阻尼系數(shù)，、為第i塊塊體的剛體加速度和剛體速度。npi為第i塊塊體的結(jié)點(diǎn)總數(shù)，F(xiàn)k為作用在第k個結(jié)點(diǎn)的外荷載，Rk為作用在不連續(xù)界面上第k個接觸點(diǎn)的接觸內(nèi)力。

2.3 時域方程離散化

在時域中，采用廣義Newmark法[16]離散化有限元平衡方程（1）與動力平衡方程（3）。將離散后的（3）回帶入方程式（3）可得時間步為n+1時塊體形心的動力平衡方程如下，

2.4 局部不連續(xù)邊界的相對位移方程

局部不連續(xù)界面的相對位移方程可寫作：

假設(shè)第i塊塊體中，不連續(xù)界面上接觸點(diǎn)對的數(shù)目為nbi，對于每一個塊體i，當(dāng)時間步為n+1時各個點(diǎn)的整體位移為：

Ti為該塊體內(nèi)結(jié)點(diǎn)上的作用力與接觸點(diǎn)對上的接觸力轉(zhuǎn)換系數(shù)矩陣，等式兩邊分別乘以整理可得：

其中，令

2.5 包含剛體轉(zhuǎn)動位移和不連續(xù)邊界接觸力的混合求解方程

聯(lián)立界面結(jié)點(diǎn)總位移方程和塊體力的平衡方程可以得到的IBE-PFE混合算法的混合迭代方程。通過引入塊體力的平衡方程后把接觸力和塊體形心剛體位移作為混合變量與結(jié)點(diǎn)總位移進(jìn)行迭代，求解接觸力增量與接觸結(jié)點(diǎn)總位移增量之間的非線性關(guān)系得以實(shí)現(xiàn)。將式（1）、式（2）代入方程式（8）：

這里引入k來描述不連續(xù)界面之間的間隙狀態(tài)，當(dāng)k等于0時，Δδn+1，L為0，表示兩接觸體之間無位移差。式（9）進(jìn)一步可改寫為：

參考方程式（2），系統(tǒng)塊體形心滿足動力平衡方程，則對所有塊體應(yīng)用式（4），可得：

聯(lián)立方程式（10）、式（11）可得：

式（12）即為IBE-PFE動力混合解法的迭代方程。求解時將形心加速度增量視為混合迭代變量與界面結(jié)點(diǎn)動力總位移進(jìn)行混合迭代求解。顯然，該方法的迭代方程組具有對稱性，可優(yōu)化求解步驟，具有較高的求解效率。

3 算例驗(yàn)證

本算例以預(yù)設(shè)切口的混凝土三點(diǎn)彎曲梁受沖擊荷載開裂為例，首先參考文獻(xiàn)[17]建立有限元模型如圖1所示，混凝土梁尺寸為228.6mm×76.2mm×25.4mm，2個支撐點(diǎn)之間的距離為203.2mm，預(yù)設(shè)裂紋長度為17.0mm，有限元模型共有414個單元、456個結(jié)點(diǎn)。在梁的中心點(diǎn)處施加沖擊荷載，動力荷載施加方向垂直于梁，作用點(diǎn)位于跨中點(diǎn)，規(guī)定沖擊荷載加載在試件中點(diǎn)處，因此混凝土梁受力狀態(tài)為對稱的，有限元模型可簡化如圖2所示，接觸界面由14組接觸點(diǎn)對模擬?；炷亮旱牟牧蠀?shù)取值如表1所示，表中ρ為密度，E為彈性模量，v為泊松比，σρ為抗拉強(qiáng)度，ε為彈簧剛度。采用PFE-IBE方法計(jì)算該含缺口的混凝土簡支梁受動力荷載作用下裂縫擴(kuò)展的數(shù)值模擬解。

表1 混凝土材料參數(shù)Tab.1 Material parameters of the concrete specimen

圖1 預(yù)設(shè)切口的三點(diǎn)彎曲梁有限元模型Fig.1 Finite element model of a notched three-point bending beam

圖2 梁左側(cè)結(jié)構(gòu)及接觸點(diǎn)對示意圖Fig.2 Configuration and finite element mesh of the left part

將IBE-PFE并與FEM結(jié)果[17]及Gopalaratnam和Shah試驗(yàn)[18]結(jié)果比較，圖3中實(shí)驗(yàn)結(jié)果的測量存在著一定延時，但與數(shù)值計(jì)算結(jié)果對比延時極?。籌BE-PFE開裂時程曲線的斜率與實(shí)驗(yàn)結(jié)果較為一致，說明模擬的開裂速率與實(shí)驗(yàn)較吻合。圖4為同樣計(jì)算條件下FEM與IBE-PFE的撓度時程曲線對比，結(jié)果吻合較好。

圖3 混凝土梁開裂時程曲線Fig.3 Solutions of the crack extension history for the concrete beam

圖4 混凝土梁撓度時程曲線Fig.4 Solutions of the deflection history for the concrete beam

4 結(jié)束語

本文基于IBE-PFE理論提出了一種新的動力迭代求解方法，并應(yīng)用于求解局部非連續(xù)變形的動接觸問題；詳細(xì)闡述了IBE-PFE動力混合算法的力學(xué)表達(dá)和公式推導(dǎo)，實(shí)現(xiàn)了分區(qū)有限元與界面剛體元的結(jié)合，將求解系統(tǒng)分為若干塊體和接觸面，由接觸狀態(tài)形成界面柔度矩陣，非線性迭代僅在接觸面上進(jìn)行，提高了計(jì)算效率的同時保證了計(jì)算精度。數(shù)值算例的結(jié)果也證明了該方法的準(zhǔn)確性和有效性。

[1] C.S.Desai，M.M.Zaman. Thin layer element for interfaces and joints[J]. International Journal for Numerical and Analytical Methods in Geomechanics 1984，8（1）：19-43.

[2] A.Gens，I.Carol，E.E.Alonso. An interface element formulation for the analysis of soil-reinforcement interaction[J]. Computers and Geotechnics 1989，7：133-151.

[3] R.Buczkowski，M.Kleiber. Elasto-plastic statistical model of strongly anisotropic rough surfaces for finite element 3D-contact analysis[J]. Computer Methods in Applied Mechanics and Engineering 2006，195：5141-5161.

[4] B.C.Burman. A numerical approach to the mechanics of discontinue[D]. James Cook University of North Qeensland，Townsville，Australia，1971.

[5] P.A.Cundall，O.D.L. Strack A discrete numerical model for granular assemblies[J]. Geotechnique 1979，29（1）：47-65.

[6] G.H.Shi，R.E.Goodman. Two dimensional discontinuous deformation analysis[J]. International Journal for Numerical and Analytical Methods in Geomechanics，1985，9：541-556.

[7] L.Jing，Y.Ma，Z.Fang. Modeling of fluid flow and solid deformation for fractured rocks with discontinuous deformation analysis（DDA）method[J]. International Journal of Rock Mechanics and Mining Sciences，2001，38：343-355.

[8] J.Liu，X.J.Kong，G.Lin. Formulations of the three-dimensional discontinuous deformation analysis method[J]，Acta Mechanica Sinica，2004，20（3）：270-282.

[9] Humberto Breves Coda. Dynamic and static non-linear analysis of reinforced media ：a BEM/FEM coupling approach[J].Computers and Structures，2001，79 ：2751-2765.

[10] J.L.Xu，Z.P.Tang. Combined Discrete Finite Element Multiscale Numerical Method and Its Application[J]. Chinese Journal of Computation Physics，2003，6 ：477-482.

[11] N.Guo，J.D.Zhao. A coupled FEM/DEM approach for hierarchical multiscale modelling of granular media[J]. International Journal for Numerical Methods in Engineering，2014，99（11）：789-818.

[12] Mark Michael，F(xiàn)rank Vogel，Bernhard Peters. DEM–FEM coupling simulations of the interactions between a tire tread and granular terrain[J]. Computer Methods in Applied Mechanics and Engineering，2015，289：227-248.

[13] X.Q. Liu，T.C. Li，L.H. Zhao. An Iteration Algorithm based on Mixed FEM and DEM for Safety Factor of Slope Stability with Determined Sliding Surface[J]. Applied Mathematics and Athematics and Information Sciences，2012，6（3）：721-725.

[14] S.G Chen，J. Zhao. Modeling of Tunnel Excavation Using a Hybrid DEM/BEM Method[J]. Computer-Aided Civil and Infrastructure Engineering，2002，17（5）：381-386.

[15] T.C. Li，X.Q. Liu，L.H. Zhao. An interactive method of interface boundary elements and partitioned finite elements for local continuous discontinuous deformation problems[J].International Journal for Numerical Methods in Engineering，2014，100（7）：534-554.

[16] WANG Xiaojian，LI Tongchun，HE Jinwen，et al. Application of IBE-PFE method to analysis on stability of slope around flood discharge tunnel outlet of a hydropower station [J].Water Resources and Hydropower Engineering，2017，48（3）：140-145，157.

[17] M.G.Katona，Zienkiewicz. A unified set of single-step algorithms，Part3：the beta-m method，a generalization of the Newmark scheme[J]. International Journal for Numerical Methods in Engineering，1985，21：1345-1359.

[18] Jiaji Du，Albert S. Kobayashi，Neil M. Hawkins.FEM Dynamic Fracture Analysis OF Concrete Beams[J].Journal of Engineering Mechanics，1989，115（10）：2136-2149.