非連續變形接觸問題的分區有限元-塊體界面元動力混合算法

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

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

0 引言

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

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

1 IBE-PFE動力混合算法的力學表達

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

2 IBE-PFE動力混合算法

2.1 動力分區有限元方程

其中，Ki為剛度矩陣；Di為阻尼矩陣，這里采用瑞利比例阻尼假定Di=αMi+βKi，α、β為瑞利比例阻尼系數；Mi為質量矩陣，Fi為作用于第i個塊體上的外荷載；Ri為由不連續界面上由內力傳輸的荷載向量，可由下式表示：

其中，Ti為該塊體內結點上的作用力與接觸點對上的接觸力轉換系數矩陣，FcL為局部坐標系下單元接觸點對之間的接觸力。

2.2 接觸塊體的剛性運動及動力平衡方程

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

其中，ωk為第k個結點相對于形心的轉換矩陣。mi為與第i塊塊體的質量，Ii為與第i塊塊體相對于形心的轉動慣量為質量阻尼系數，、為第i塊塊體的剛體加速度和剛體速度。npi為第i塊塊體的結點總數，Fk為作用在第k個結點的外荷載，Rk為作用在不連續界面上第k個接觸點的接觸內力。

2.3 時域方程離散化

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

2.4 局部不連續邊界的相對位移方程

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

假設第i塊塊體中，不連續界面上接觸點對的數目為nbi，對于每一個塊體i，當時間步為n+1時各個點的整體位移為：

Ti為該塊體內結點上的作用力與接觸點對上的接觸力轉換系數矩陣，等式兩邊分別乘以整理可得：

其中，令

2.5 包含剛體轉動位移和不連續邊界接觸力的混合求解方程

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

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

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

聯立方程式（10）、式（11）可得：

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

3 算例驗證

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

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

圖1 預設切口的三點彎曲梁有限元模型Fig.1 Finite element model of a notched three-point bending beam

圖2 梁左側結構及接觸點對示意圖Fig.2 Configuration and finite element mesh of the left part

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

圖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 結束語

本文基于IBE-PFE理論提出了一種新的動力迭代求解方法，并應用于求解局部非連續變形的動接觸問題；詳細闡述了IBE-PFE動力混合算法的力學表達和公式推導，實現了分區有限元與界面剛體元的結合，將求解系統分為若干塊體和接觸面，由接觸狀態形成界面柔度矩陣，非線性迭代僅在接觸面上進行，提高了計算效率的同時保證了計算精度。數值算例的結果也證明了該方法的準確性和有效性。

[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，Frank 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.