三維接觸面單元彈塑性損傷數(shù)值模擬分析方法

鄧建，肖明，謝冰冰(1.武漢大學 水資源與水電工程科學國家重點實驗室，湖北 武漢，430072；2.武漢大學 水工巖石力學教育部重點實驗室，湖北 武漢，430072)

三維接觸面單元彈塑性損傷數(shù)值模擬分析方法

鄧建1,2，肖明1,2，謝冰冰1,2
(1.武漢大學 水資源與水電工程科學國家重點實驗室，湖北 武漢，430072；
2.武漢大學 水工巖石力學教育部重點實驗室，湖北 武漢，430072)

針對巖土工程中出現(xiàn)的非連續(xù)結(jié)構(gòu)面、接觸面等問題，建立一種新的模擬結(jié)構(gòu)面或接觸問題的三維非線性接觸單元模型；針對接觸面非線性變形以及損傷演化特性，提出一種可以考慮接觸面損傷演化過程的彈塑性損傷數(shù)值模擬方法；針對接觸面張開、閉合以及剪切滑移失效等變形特性，模擬接觸面的黏結(jié)、滑移和開裂等非線性現(xiàn)象。針對某水電站地下廠房巖錨吊車梁圍巖與混凝土接觸問題進行數(shù)值分析。研究結(jié)果表明：采用接觸面單元模擬混凝土結(jié)構(gòu)與圍巖的相互作用，吊車梁受力情況良好，應力、變形分布正常，基本規(guī)律與實際相符，在工程實踐中是合理可行的。

巖土工程；非線性；接觸單元；彈塑性；損傷演化；數(shù)值模擬

在巖土結(jié)構(gòu)工程中，存在著大量接觸問題，如巖體的結(jié)構(gòu)面、混凝土結(jié)構(gòu)的分縫以及結(jié)構(gòu)與基礎、不同材料的交界面等[1]。這些接觸面表現(xiàn)出不同于一般巖體的力學特性，即在荷載作用下，可能發(fā)生張開、滑移、閉合等現(xiàn)象[2]，對結(jié)構(gòu)的安全穩(wěn)定有較大影響，因此越來越受到人們的重視。長期以來，許多中外學者對接觸問題進行了多方面研究，先后基于不同的理論提出了各種接觸面計算模型。GOODMAN等[3]基于彈簧剛度提出了用于模擬節(jié)理巖體的無厚度接觸面單元，得到了廣泛應用，其主要問題是彈簧剛度系數(shù)選取的任意性及其合理取值的困難性。DESAI等[4]提出了有厚度的非線性薄層單元，在一定程度上克服了GOODMAN接觸面單元的缺點，其界面參數(shù)包括彈性模量E、剪切模量G、泊松比γ這3個參數(shù)，但在參數(shù)取值方面缺乏理論依據(jù)。BOULON等[5]提出了接觸面彈塑性模型，但是由于較復雜而未能得到推廣應用。CLOUGH等[6]在實驗基礎上提出了非線性彈性模型。武亞軍等[7]將非線性彈性理論與彈塑性理論相結(jié)合，提出了一種非線性彈性-理想塑性本構(gòu)模型。張冬霽等[8]在實驗基礎上提出了剪切錯動帶單元。盧廷浩等[9]基于剪切錯動帶提出了接觸面法向剛度與切向剛度相耦合的非線性本構(gòu)模型。高俊合等[10]通過大型單剪實驗提出了有厚度剪切滑移薄層單元。胡黎明等[11]進行了系統(tǒng)的土與結(jié)構(gòu)物接觸面的直剪試驗,根據(jù)試驗結(jié)果建立了描述接觸面應變軟化和剪脹現(xiàn)象的損傷力學模型。安關峰等[12]提出了接觸面的三維彈-粘-塑性模型。本文作者在以上研究的基礎上，提出了一種具有滑移失效及拉伸、剪切、黏結(jié)特性的三維接觸面彈塑性損傷有限元計算方法，并通過工程實例驗證該方法的合理性。

1　接觸面單元基本模型

圖1所示為接觸面單元模型。考慮圖1中三維8節(jié)點等厚度接觸單元，沿法向zˊ方向分為上、下2平面，接觸單元厚度為t，t/ b=1/100～1/10，b為接觸單元的長度。形函數(shù)滿足

在有限元計算中，整體坐標與接觸面的局部坐標往往是不一致的，在三維坐標系中，由Oxyz轉(zhuǎn)換到Oˊxˊyˊzˊ坐標系時，其坐標轉(zhuǎn)換矩陣[R]可表示為

圖1　接觸面單元模型Fig.1　Contact element model

假設薄層單元沿厚度方向即法向zˊ方向應變均勻分布，以接觸面單元對應點的相對位移表征單元的應力應變[13]，其相對位移為

則接觸面上表面任意一點的相對位移為：

則單元應變?yōu)椋?/p>

若接觸面單元的應力矩陣為[D]，則根據(jù)虛位移原理，有

考慮到厚度t盡可能小，接觸面單元整體剛度矩陣可表示為

2　接觸面單元彈塑性本構(gòu)

接觸面單元在復雜外荷載作用下其應力應變關系常表現(xiàn)出一定的非線性，在發(fā)生剪切變形和受高壓應力作用時，可能產(chǎn)生沿粗糙體表面的滑移以及粗糙體剪斷、壓碎和磨損、拉裂等現(xiàn)象，接觸面單元的變形中不僅有可逆的彈性變形，而且有不可逆的塑性變形，本文采用接觸面單元彈塑性本構(gòu)模型進行數(shù)值模擬[14-15]。

根據(jù)彈塑性理論，當空間某一點的應力狀態(tài)進入屈服后，其應力應變關系表現(xiàn)為非線性，此時應變增量可分解為彈性應變增量和塑性應變增量兩部分，即

其中：彈性部分滿足虎克定律，塑性部分滿足塑性流動法則，即

式中：F為屈服函數(shù)。

對于彈性部分，忽略法向和切向變形的耦合影響，只考慮接觸面單元法向和切向的單獨承載響應，采用彈性應力矩陣：

對于塑性部分，采用下述混合準則：當接觸面產(chǎn)生剪切破壞時，采用Mohr-Coulomb準則；當接觸面產(chǎn)生受拉破壞時，因為其抗拉強度低，承受法向拉應力時，容易開裂，所以，受拉破壞的強度條件采用最大拉應力準則，即

根據(jù)變形相容條件，有

式中：F為屈服函數(shù)，H為硬化參數(shù)。可得 λd為

對于彈塑性應力增量則可以表示為

其中，塑性應力矩陣可表示為

3　接觸面單元損傷開裂模擬

3.1接觸面損傷塑性模型

在外荷載作用下，接觸面除了表現(xiàn)出非線性變形特性外，還會逐漸產(chǎn)生微裂紋和微孔隙等損傷破壞情況，損傷的發(fā)生和積累導致結(jié)構(gòu)材料剛度劣化和強度降低。為了模擬這種隨著塑性變形的產(chǎn)生、積累而形成的損傷演化過程，本文在彈塑性本構(gòu)模型的基礎上，通過引入損傷變量[16]，建立了接觸面單元彈塑性損傷本構(gòu)模型。

在承載過程中，由于出現(xiàn)新的損傷，并出現(xiàn)新的塑性變形，應變增量可分解為3部分：

d{ε}d和d{ε}p分別是損傷應變和塑性應變增量，其滿足廣義正交法則：

引入損傷變量ω，則損傷應變增量可以定義如下：

假設接觸面的損傷變形和塑性變形是同時出現(xiàn)的，當應力狀態(tài)滿足屈服準則F(σ, ω)=0時,才有可能出現(xiàn)新的損傷變形和塑性變形。

在不考慮黏性效應的情況下，可以假定損傷變量增量dω與塑性應變增量之間滿足以下關系：

式中：Tn是一個三維矢量，它是損傷變量ω和單元應力σ的函數(shù)。通過式(21)~(23)可得

采用與不考慮損傷劣化相同的本構(gòu)矩陣形式，即式(18)所示，其中H可以表示為

10月24日，“經(jīng)協(xié)”籌備會召開座談會，提出經(jīng)濟民主的九點主張，并于11月5日在《商務日報》上發(fā)表。這年冬，“經(jīng)協(xié)”籌備會舉行歡迎中共代表團的小型集會。周恩來、鄧穎超、陸定一等參加，周恩來講了話。這次會議對“經(jīng)協(xié)”籌備會的負責人鼓舞很大。

式中：矩陣[D]和硬化參數(shù)H均含有損傷變量ω。再由式(9)和式(10)計算接觸面單元剛度矩陣。

3.2接觸面破壞狀態(tài)模擬

由于接觸面單元在平行和垂直于接觸面方向分別表現(xiàn)出不同的物理力學特性，所以，在彈塑性非線性分析時，應對平行于接觸面和垂直于接觸面2個方向的破壞特性進行判斷。

3.2.1拉裂破壞

對垂直于接觸面的破壞情況，先計算接觸面單元局部坐標下的應力：

[S]為應力轉(zhuǎn)化矩陣。然后根據(jù)局部坐標下的應力分量校核垂直于接觸面的正應力狀態(tài)，當

說明垂直接觸面的正應力超過極限抗拉強度Rt，垂直于接觸面方向?qū)㈤_裂。此時必須將超出接觸面抗拉強度的應力轉(zhuǎn)化成節(jié)點荷載，并將其轉(zhuǎn)移。對于發(fā)生拉裂破壞的單元，抗拉強度Rt設置為零，單元不再進行損傷處理，直接標記為分離狀態(tài)。

對平行于接觸面的破壞情況，應校核是否沿接觸面滑移。當接觸面單元的應力狀態(tài)滿足

式中：f為接觸面的摩擦系數(shù)；c為黏結(jié)力；k為接觸面滑動安全系數(shù)。此時平行于接觸面方向?qū)a(chǎn)生剪切滑動，同時應將超出抗剪強度的剪力：

轉(zhuǎn)化成節(jié)點荷載，并將其轉(zhuǎn)移。對于發(fā)生剪切破壞的單元，則把單元標記為滑移狀態(tài)，抗拉強度Rt修正為殘余強度。

3.2.3拉裂破壞后閉合

再接觸狀態(tài)是指已經(jīng)拉裂破壞的接觸單元再次受壓時會產(chǎn)生閉合，恢復為接觸狀態(tài)，在力學性質(zhì)上表現(xiàn)為法向剛度和切向剛度的恢復。再接觸狀態(tài)判斷可以采用如下方法：

如果單元法向位移滿足上式，那么單元處于再接觸狀態(tài)，抗拉強度Rt同樣設置為0 MPa。

4　工程實例

4.1工程概況黃登水電站位于云南省怒江州蘭坪縣境內(nèi)，為瀾滄江上游河段規(guī)劃中的第5個梯級，引水發(fā)電系統(tǒng)位于樞紐左岸，主要建筑物有主副廠房、主變室、尾水調(diào)壓井等3大洞室及電站進水口、壓力管道、機組尾水檢修閘門室、尾水隧洞及其他相關輔助洞室。整個廠區(qū)洞室縱橫交錯規(guī)模宏大，屬大型地下洞室群。其中引水隧洞(包括鋼管段)共4條，平行布置，中心間距為31 m，單條隧洞長244.54 m；主副廠房開挖長× 寬×高為228.0 m×30.3 m×76.5 m；主變室開挖長× 寬×高為153 m×18 m×30 m；尾水調(diào)壓井采用圓筒阻抗式，內(nèi)徑為26 m，高為73 m，阻抗孔直徑為8.5 m；調(diào)壓井后接兩條尾水隧洞，1號隧洞長度為422.092 m，隧洞長度為323.339 m，洞徑為15 m。

4.2計算模型和計算條件

建立吊車梁局部細化模型進行計算分析。綜合分析本工程的地形、地質(zhì)特性，取2號機組段進行巖錨吊車梁的穩(wěn)定分析，模型包括2號機組段的主廠房、引水管、母線洞、主變洞、尾水管等。模型全部采用八節(jié)點六面體單元進行離散，共剖分了70 496個單元和74 978個節(jié)點。其中澆筑混凝土單元11 200個，接觸單元1 120個。三維有限元模型開挖單元如圖2所示，巖錨吊車梁混凝土單元及接觸單元模型如圖3所示。

圖2　三維有限元模型Fig.2　Three-dimensional finite element model

圖3　混凝土及接觸單元模型Fig.3　Concrete and contact element model

初始地應力場通過實測地應力反演獲得，屬中等略低地應力。地下廠房區(qū)域巖體條件較好，以Ⅱ和Ⅲ類圍巖為主，材料力學參數(shù)取值如表1所示。

表1　材料力學參數(shù)取值Table 1　Mechanical parameters of materials

計算工況：1)采用等效連續(xù)模型分析洞室群開挖對吊車梁混凝土受力影響；2)采用接觸面單元模型分析洞室群開挖對吊車梁混凝土受力影響。

4.3計算結(jié)果分析

為了論證接觸面單元彈塑性損傷數(shù)值模擬分析方法的合理性，本文分別采用等效連續(xù)模型和接觸面單元模型進行數(shù)值模擬，選取吊車梁機組段典型斷面，對比分析地下廠房開挖完畢巖錨吊車梁混凝土應力和變形分布規(guī)律。

4.3.1巖錨吊車梁應力分布規(guī)律

分別采用等效連續(xù)模型和接觸面單元模型進行數(shù)值模擬，廠房開挖完畢后巖錨吊車梁第1、第3主應力(拉為正，壓為負)如圖4~7所示。

采用等效連續(xù)模型進行分析計算(見圖4和圖5)，開挖完畢吊車梁混凝土第1主應力分布在-0.07~-2.52 MPa之間，第3主應力分布在-0.31~1.81 MPa之間。從圖4和圖5可以看出：第1主應力不大，但第3主應力局部超過混凝土的抗拉強度，發(fā)生局部開裂。

圖4　等效連續(xù)模型第1主應力Fig.4 The 1st principle stress with equivalent continuum model

采用接觸面單元模型進行分析計算(見圖6和圖7)，開挖完畢吊車梁混凝土第1主應力分布在-0.05~ -1.31 MPa之間，第3主應力分布在-0.10~1.13 MPa之間。從圖6和圖7可以看出：吊車梁混凝土應力值較小，第3主應力均小于混凝土的抗拉強度。

圖5　等效連續(xù)模型第3主應力Fig.5　The 3rd principle stress with equivalent continuum model

圖6　接觸面單元模型第1主應力Fig.6　The 1st principle stress with contact element model

圖7　接觸面單元模型第三主應力Fig.7　The 3rd principle stress with contact element model

從圖6和圖7可以看出：采用接觸面單元模擬吊車梁混凝土與圍巖的接觸狀態(tài)，吊車梁受力情況良好，應力分布均勻，能夠有效控制后期開挖對吊車梁的擾動影響，且與工程實際更為吻合。

4.3.2巖錨吊車梁變形分布規(guī)律

分別采用等效連續(xù)模型和接觸面單元模型進行數(shù)值模擬，廠房開挖完畢巖錨吊車梁位移分布如圖8和圖9所示。

圖8　等效連續(xù)模型位移分布Fig.8　Displacements distribution with equivalent continuum model

圖9　接觸面單元模型位移分布Fig.9　Displacements distribution with contact element model

采用等效連續(xù)模型進行分析計算(見圖8)，廠房開挖完畢時，吊車梁變形在11.5~12.1 mm之間，最大變形量出現(xiàn)在吊車梁自由面，接觸面部位位移從上到下逐漸增大。吊車梁上下游位移和分布規(guī)律均相似，基本以廠房中心線呈對稱分布。總體來看，吊車梁下部區(qū)域巖體和母線洞巖體開挖對吊車梁影響較大。

采用接觸面單元模型進行分析計算(見圖9)，廠房開挖完畢時，吊車梁變形在7.6~8.8 mm之間，吊車梁變形均勻，變形較小，上下游位移和分布規(guī)律基本一致。

從圖8和圖9可以看出：采用接觸面單元模擬吊車梁混凝土與圍巖的接觸狀態(tài)，吊車梁受力良好，變形規(guī)律正常，以隨圍巖整體變形為主，變形量減小約3 mm，更加符合工程實際，有一定的模擬效果。

5　結(jié)論

1)基于Desai薄層單元的基本理論，建立了一種新的模擬接觸問題的三維非線性接觸單元模型，并基于非線性損傷演化基本理論，提出了接觸面單元彈塑性損傷開裂模型。

2)采用接觸面單元模擬混凝土結(jié)構(gòu)與圍巖的相互作用，吊車梁受力情況良好，應力、變形分布正常，基本規(guī)律與實際相符，在工程實踐中是合理可行的，并為巖土工程非連續(xù)性接觸問題的模擬提供了一種有效的分析方法。

[1]張楚漢.論巖石、混凝土離散-接觸-斷裂分析[J].巖石力學與工程學報,2008,27(2):217-220. ZHANG Chuhan.Discrete-contact-fracture analysis of rock and concrete[J].ChineseJournalofRockMechanicsand Engineering,2008,27(2):217-220.

[2]金峰,邵偉,張立翔,等.模擬軟弱夾層動力特性的薄層單元及其工程應用[J].工程力學,2002,19(2):36-37. JIN Fen,SHAO Wei,ZHANG Lixiang,et al.A thin-layer element for simulation of static and dynamic characteristics of soft interlayer and its application[J].Engineering Mechanics, 2002,19(2):36-37.

[3]GOODMAN R E,TAYLOR R L,BREKKE T L.A model for the mechanics of jointed rock[J].Journal of Soil Mechanics and Foundation Engineering Division,1968,94(3):637-660.

[4]DESAI C S,ZAMAN M M,LIGHTNER J G,et al.Thin-layer element for interfaces and joints[J].International Journal for Numerical and Analytical Methods in Geomechanics,1984,8(1): 19-43.

[5]BOULON M,NOVA R.Modeling of soil structure interface behavior:a comparison between elastoplastic and rate-type laws[J].Computers and Geotechnics,1990,17(9):21-46.

[6]CLOUGH G W,DUNCAN J M.Finite element analysis of retaining wall behavior[J].Journal of Soil Mechanics and Foundation Engineering Division,1971,97(12):1657-1673.

[7]武亞軍,欒茂田,楊敏.土與結(jié)構(gòu)間一種新的接觸單元模型[J].同濟大學學報(自然科學版),2005,33(4):432-435. WU Yajun,LUAN Maotian,YANG Ming.New contact element of interface between soil and structure[J].Journal of Tongji University(Natural Science),2005,33(4):432-435.

[8]張冬霽,盧廷浩.一種土與結(jié)構(gòu)接觸面模型的建立及其應用[J].巖土工程學報,1998,20(6):62-66. ZHANG Dongji,LU Tinghao.Establishment and application of an interface model between soil and structure[J].Chinese Journal of Geotechnical Engineering,1998,20(6):62-66.

[9]盧廷浩,鮑伏波.接觸面薄層單元耦合本構(gòu)模型[J].水利學報,2000(2):71-75. LU Tinghao,BAO Fubo.A coupled constitutive model for interface thin layer element[J].Journal of Hydraulic Engineering, 2000(2):71-75.

[10]高俊合,于海學,趙維炳,等.土與混凝土接觸面特性的大型單剪試驗研究及數(shù)值模擬[J].土木工程學報,2000,33(4): 42-46. GAO Junhe,YU Haixue,ZHAO Weibing,et al.Characteristics of interface between soil and concrete by using large size single shearapparatusandnumericalanalysis[J].ChinaCivil Engineering Journal,2000,33(4):42-46.

[11]胡黎明,濮家騮.土與結(jié)構(gòu)物接觸面損傷本構(gòu)模型[J].巖土力學,2002,23(1):6-11. HULiming,PUJialiu.Damagemodelofsoil-structure interface[J].Rock and Soil Mechanics,2002,23(1):6-11.

[12]安關峰,高大釗.接觸面彈粘塑性本構(gòu)關系研究[J].土木工程學報,2001,34(1):88-91. AN Guanfeng,GAO Dazhao.Research on elastic-visco-plastic constitution of interfaces[J].China Civil Engineering Journal, 2001,34(1):88-91.

[13]SAMADHIYA N K,VILADKAR M N,AL-OBAYDI M A. Three-dimensional joint/interface element for rough undulating major discontinuities in rock masses[J].International Journal of Geomechanics,2008,8(6):328-331.

[14]WANG J G,ICHIKAWA Y,LEUNG C F.A constitutive model for rock interfaces and joints[J].International Journal of Rock Mechanics and Mining Sciences,2003,40(1):41-53.

[15]肖明.地下洞室圍巖穩(wěn)定與支護數(shù)值方法研究[D].武漢:武漢大學水利水電學院,2002:56-61. XIAO Ming.Study on numerical analysis method of stability and supporting for underground caverns[D].Wuhan:Wuhan University.SchoolofWaterResourcesandHydropower Engineering,2002:56-61.

[16]殷有泉.考慮損傷的節(jié)理本構(gòu)模型[J].工程地質(zhì)學報,1994, 2(4):2-4. YIN Youquan.Joint constitutive model considering damage[J]. Journal of Engineering Geology,1994,2(4):2-4.

(編輯羅金花)

Three-dimensional elastic-plastic damage numerical simulation and analysis method for contact elements

DENG Jian1,2,XIAO Ming1,2,XIE Bingbing1,2
(1.State Key Laboratory of Water Resources and Hydropower Engineering Science,Wuhan University, Wuhan 430072,China;
2.Key Laboratory of Rock Mechanics in Hydraulic Structural Engineering of Ministry of Education, Wuhan University,Wuhan 430072,China)

In light of discontinuous joints and interfaces in geotechnical engineering,a new three-dimensional nonlinear contact element model for the simulation of joints and contact problems was presented.Based on the nonlinear deformation and damage evolution characteristics of contact elements,an elastic-plastic damage numerical simulation method which takes the damage evolution process into consideration was proposed.In view of the deformation characteristics of contact elements such as opening,closing and shearing slip damage,some nonlinear phenomena like bonding,slip and cracking were simulated.The contact problem between surrounding rocks and concrete crane beam of underground powerhouse was taken as an example.The results show that concrete crane beam is of better stress condition and normal deformation distribution with the contact element model,and the basic laws conform to the actual situation, which verifies the good rationality and accuracy of the present model in engineering practice.

geotechnical engineering;nonlinear;contact elements;elastic-plastic;damage evolution;numerical simulation

鄧建，博士研究生，從事地下結(jié)構(gòu)穩(wěn)定數(shù)值分析研究；E-mail:xianyudj@126.com

TU45

A

1672-7207(2016)07-2383-07

10.11817/j.issn.1672-7207.2016.07.028

2015-07-23；

2015-09-23

國家自然科學基金重大資助項目(91215301)；國家自然科學基金資助項目(51279136,51209164)(Project(91215301) supported by the Major Program of National Natural Science Foundation of China;Projects(51279136,51209164)supported by National Natural Science Foundation of China)