李同春,樊舒婕,劉曉青,趙蘭浩
(1.河海大學(xué)水利水電學(xué)院,江蘇省南京市 210098;2.水資源高效利用與工程安全國(guó)家工程研究中心,江蘇省南京市 210098)
在水利工程中存在著許多局部不連續(xù)的現(xiàn)象,例如拱壩的橫縫、巖體的裂隙、地基的軟弱面等。這些局部不連續(xù)區(qū)域在地震等動(dòng)力荷載作用下的張開(kāi)、閉合和摩擦滑移關(guān)系著結(jié)構(gòu)穩(wěn)定與工程安全??紤]動(dòng)載下不連續(xù)區(qū)域交界面的動(dòng)態(tài)接觸問(wèn)題對(duì)于結(jié)構(gòu)動(dòng)力分析有著十分重要的意義。對(duì)于此類(lèi)問(wèn)題,已有的數(shù)值分析方法包括傳統(tǒng)的有限元法FEM[1-3]、離散元法DEM[4,5]和非連續(xù)變形分析方法DDA[6-8],以及FEM-BEM[9]、DEM-FEM[10-13]、DEM-BEM[14]等混合算法。
本文針對(duì)動(dòng)力情況下非連續(xù)變形的接觸非線性問(wèn)題,基于文獻(xiàn)[15]的力學(xué)模型提出了一種新的動(dòng)力迭代求解方法。推導(dǎo)了包含接觸力、塊體形心剛體加速度和接觸結(jié)點(diǎn)總位移在內(nèi)的非線性方程,將該非線性方程與塊體形心動(dòng)力平衡方程進(jìn)行耦合,建立了以界面接觸力和塊體形心剛體加速度為混合變量與界面結(jié)點(diǎn)總位移進(jìn)行迭代求解的動(dòng)力界面迭代方程,給出了混合解法的實(shí)施步驟并通過(guò)典型算例驗(yàn)證了動(dòng)力混合解法的正確性。
IBE-PFE動(dòng)力混合算法將求解系統(tǒng)分解為若干塊體與非連續(xù)界面,每個(gè)塊體同時(shí)作為系統(tǒng)有限元的一個(gè)分區(qū),以塊體結(jié)點(diǎn)加速度作為分區(qū)有限元計(jì)算基本變量,以塊體剛體形心加速度作為剛體運(yùn)動(dòng)變量,當(dāng)接觸界面上接觸點(diǎn)對(duì)的接觸內(nèi)力增量已知,則可通過(guò)動(dòng)力有限元法求出塊體的加速度增量,進(jìn)而根據(jù)Newmark法[16]求解各結(jié)點(diǎn)的應(yīng)力和位移。建立接觸界面上力與點(diǎn)對(duì)總位移方程,在塊體形心處建立接觸內(nèi)力、塊體形心加速度增量引起的力的平衡方程,將以上方程進(jìn)行耦合,塊體形心接觸內(nèi)力、加速度作為混合迭代變量和接觸界面點(diǎn)對(duì)總位移的求解可以通過(guò)迭代方程求得。
其中,Ki為剛度矩陣;Di為阻尼矩陣,這里采用瑞利比例阻尼假定Di=αMi+βKi,α、β為瑞利比例阻尼系數(shù);Mi為質(zhì)量矩陣,F(xiàn)i為作用于第i個(gè)塊體上的外荷載;Ri為由不連續(xù)界面上由內(nèi)力傳輸?shù)暮奢d向量,可由下式表示:
其中,Ti為該塊體內(nèi)結(jié)點(diǎn)上的作用力與接觸點(diǎn)對(duì)上的接觸力轉(zhuǎn)換系數(shù)矩陣,F(xiàn)cL為局部坐標(biāo)系下單元接觸點(diǎn)對(duì)之間的接觸力。
如果第i塊的形心位移被設(shè)為剛體位移變量,那么根據(jù)剛體運(yùn)動(dòng)法則,可視為塊體內(nèi)的每個(gè)點(diǎn)的剛體位移是與此相關(guān)、且遵守剛體運(yùn)動(dòng)法則的。在每一個(gè)接觸塊體形心點(diǎn)上分別建立塊體力的平衡方程,對(duì)于每個(gè)塊體,離散單元的動(dòng)力平衡方程為:
其中,ωk為第k個(gè)結(jié)點(diǎn)相對(duì)于形心的轉(zhuǎn)換矩陣。mi為與第i塊塊體的質(zhì)量,Ii為與第i塊塊體相對(duì)于形心的轉(zhuǎn)動(dòng)慣量為質(zhì)量阻尼系數(shù),、為第i塊塊體的剛體加速度和剛體速度。npi為第i塊塊體的結(jié)點(diǎn)總數(shù),F(xiàn)k為作用在第k個(gè)結(jié)點(diǎn)的外荷載,Rk為作用在不連續(xù)界面上第k個(gè)接觸點(diǎn)的接觸內(nèi)力。
在時(shí)域中,采用廣義Newmark法[16]離散化有限元平衡方程(1)與動(dòng)力平衡方程(3)。將離散后的(3)回帶入方程式(3)可得時(shí)間步為n+1時(shí)塊體形心的動(dòng)力平衡方程如下,
局部不連續(xù)界面的相對(duì)位移方程可寫(xiě)作:
假設(shè)第i塊塊體中,不連續(xù)界面上接觸點(diǎn)對(duì)的數(shù)目為nbi,對(duì)于每一個(gè)塊體i,當(dāng)時(shí)間步為n+1時(shí)各個(gè)點(diǎn)的整體位移為:
Ti為該塊體內(nèi)結(jié)點(diǎn)上的作用力與接觸點(diǎn)對(duì)上的接觸力轉(zhuǎn)換系數(shù)矩陣,等式兩邊分別乘以整理可得:
其中,令
聯(lián)立界面結(jié)點(diǎn)總位移方程和塊體力的平衡方程可以得到的IBE-PFE混合算法的混合迭代方程。通過(guò)引入塊體力的平衡方程后把接觸力和塊體形心剛體位移作為混合變量與結(jié)點(diǎn)總位移進(jìn)行迭代,求解接觸力增量與接觸結(jié)點(diǎn)總位移增量之間的非線性關(guān)系得以實(shí)現(xiàn)。將式(1)、式(2)代入方程式(8):
這里引入k來(lái)描述不連續(xù)界面之間的間隙狀態(tài),當(dāng)k等于0時(shí),Δδn+1,L為0,表示兩接觸體之間無(wú)位移差。式(9)進(jìn)一步可改寫(xiě)為:
參考方程式(2),系統(tǒng)塊體形心滿(mǎn)足動(dòng)力平衡方程,則對(duì)所有塊體應(yīng)用式(4),可得:
聯(lián)立方程式(10)、式(11)可得:
式(12)即為IBE-PFE動(dòng)力混合解法的迭代方程。求解時(shí)將形心加速度增量視為混合迭代變量與界面結(jié)點(diǎn)動(dòng)力總位移進(jìn)行混合迭代求解。顯然,該方法的迭代方程組具有對(duì)稱(chēng)性,可優(yōu)化求解步驟,具有較高的求解效率。
本算例以預(yù)設(shè)切口的混凝土三點(diǎn)彎曲梁受沖擊荷載開(kāi)裂為例,首先參考文獻(xiàn)[17]建立有限元模型如圖1所示,混凝土梁尺寸為228.6mm×76.2mm×25.4mm,2個(gè)支撐點(diǎn)之間的距離為203.2mm,預(yù)設(shè)裂紋長(zhǎng)度為17.0mm,有限元模型共有414個(gè)單元、456個(gè)結(jié)點(diǎn)。在梁的中心點(diǎn)處施加沖擊荷載,動(dòng)力荷載施加方向垂直于梁,作用點(diǎn)位于跨中點(diǎn),規(guī)定沖擊荷載加載在試件中點(diǎn)處,因此混凝土梁受力狀態(tài)為對(duì)稱(chēng)的,有限元模型可簡(jiǎn)化如圖2所示,接觸界面由14組接觸點(diǎn)對(duì)模擬?;炷亮旱牟牧蠀?shù)取值如表1所示,表中ρ為密度,E為彈性模量,v為泊松比,σρ為抗拉強(qiáng)度,ε為彈簧剛度。采用PFE-IBE方法計(jì)算該含缺口的混凝土簡(jiǎn)支梁受動(dòng)力荷載作用下裂縫擴(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)對(duì)示意圖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é)果的測(cè)量存在著一定延時(shí),但與數(shù)值計(jì)算結(jié)果對(duì)比延時(shí)極??;IBE-PFE開(kāi)裂時(shí)程曲線的斜率與實(shí)驗(yàn)結(jié)果較為一致,說(shuō)明模擬的開(kāi)裂速率與實(shí)驗(yàn)較吻合。圖4為同樣計(jì)算條件下FEM與IBE-PFE的撓度時(shí)程曲線對(duì)比,結(jié)果吻合較好。
圖3 混凝土梁開(kāi)裂時(shí)程曲線Fig.3 Solutions of the crack extension history for the concrete beam
圖4 混凝土梁撓度時(shí)程曲線Fig.4 Solutions of the deflection history for the concrete beam
本文基于IBE-PFE理論提出了一種新的動(dòng)力迭代求解方法,并應(yīng)用于求解局部非連續(xù)變形的動(dòng)接觸問(wèn)題;詳細(xì)闡述了IBE-PFE動(dòng)力混合算法的力學(xué)表達(dá)和公式推導(dǎo),實(shí)現(xiàn)了分區(qū)有限元與界面剛體元的結(jié)合,將求解系統(tǒng)分為若干塊體和接觸面,由接觸狀態(tài)形成界面柔度矩陣,非線性迭代僅在接觸面上進(jìn)行,提高了計(jì)算效率的同時(shí)保證了計(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.