接触面损伤模型应用于三维有限元分析.doc

接触面损伤模型应用于三维有限元分析胡黎明，濮家骝，王 刚(清华大学 水利水电工程系)摘要：本文简要介绍了作者提出的接触面损伤模型，推导了接触面单元的三维有限元格式。将损伤模型接触面单元应用于三维有限元计算，对滑块问题进行了数值试验，比较了不同接触面模型的计算结果，验证了接触面损伤模型分析土-结构相互作用问题的合理性和优越性。文中对一围堰工程进行了有限元分析，计算结果合理可靠，表明本文提出的损伤模型接触面单元可以方便地应用于实际工程问题的数值模拟。关键词：土与结构物接触面；应变软化；损伤模型；有限元分析作者简介：胡黎明(1974-)，男，山西文水人，清华大学讲师，主要从事岩土工程研究。不同材料接触面的力学特性及其数值模拟是岩土工程长期研究的重要课题。Goodman1对岩石节理裂隙的特性进行了研究，认为可以用4节点无厚度节理单元来模拟节理的滑移和张裂等变形特性。Clough和Duncan提出非线性弹性模型描述应力与相对位移之间的关系2。Ghaboussi提出了基于弹性理论分析岩石节理相互作用面的单元，但由于采用相对自由度作为基本未知量，实际应用不方便，不适合推广使用3。Desai认为两种材料接触面剪应力传递和剪切带的形成均发生在接触面附近的一薄层土体中，从而提出薄层单元的概念，用于模拟土与结构物接触面及岩石节理的粘接、滑动、张开和闭合等各种接触状态4。殷宗泽提出接触面单元的变形由两个部分组成，一部分是土体的基本变形，其变形特性与周围土体单元相同；另一部分是接触面上的破坏变形，包括滑动破坏和拉裂破坏，总的变形是两者的叠加5。胡黎明建立了一个能够反映接触面变形机理和力学特性的接触面损伤模型，可以模拟其剪切过程中的应变软化和剪胀现象6。本文将损伤模型应用于有限元计算，验证其分析土与结构相互作用问题的合理性和有效性。1 损伤模型接触面单元1.1 接触面非线性弹性模型 以往的接触面本构关系不考虑切向和法向的耦合作用，接触面应力和相对位移关系为：(1)Clough和Duncan提出的接触面非线性弹性模型为：当接触面受压，n0时Kss=Ks(n/Pa)ns(1-Rfs/ntans+cs)2, Knn=106MPa(2)当接触面受拉，n0时Knn=Kss=0.01MPa(3)接触面模型参数有5个：剪切劲度参数Ks、ns，接触面强度参数cs、s和破坏比Rfs。1.2 接触面损伤模型 胡黎明进行了系统的土与结构物接触面的直剪试验，研究了接触面力学性质和变形机理6，7。试验表明，粗糙接触面剪切破坏时在接触面附近产生应变局部化形成剪切带，并伴有应变软化和剪胀现象。剪切带厚度t5D50。剪切应变和剪胀均发生在剪切带内。根据试验结果建立了描述接触面应变软化和剪胀现象的损伤力学模型，认为某一状态接触面由无损伤状态和临界状态组合而成，损伤演化过程为无损伤状态向临界状态的逐渐转化过程。1.2.1 无损伤状态 无损伤状态接触面采用弹塑性本构模型来描述，其屈服函数为F=i2/Ds+ps(i-sign(i)5tg0)=0(4)式中：0为摩擦角，ps为塑性剪切应变，Ds为剪切弹性模量，上标i表示无损伤状态。采用不相适应流动法则，设塑性势函数Q=i，即法向无塑性变形，只有弹性变形。采用接触面塑性剪切应变为硬化参数，即H=ps.1.2.2 临界状态 临界状态接触面可用理想塑性模型来描述。其应力-应变关系为(5)式中：r为残余摩擦角；cv为临界状态体积应变；0v为土的最大体应变；K为体积应变参数，与粗糙度R有关，与正应力无关；Dr为土体初始相对密度，0为常数，一般取为10kPa；上标c表示临界状态。1.2.3 损伤状态变量演化律 土与结构物接触面剪切应变是造成损伤的直接原因。假设损伤演化过程只与接触面塑性剪切应变有关，损伤状态变量为D=1-exp(-apsb)(6)式中：a、b为与R相关的参数，与正应力无关。根据损伤力学的基本原理，接触面应力-应变状态变量为无损伤和临界状态的线性叠加，即(7)对接触面应力-应变状态变量之间的关系式取微分，可推导得到增量应力-应变关系：d=Depd(8)接触面应力-应变关系矩阵，即弹塑性损伤矩阵为：(9)式中：，分别为接触面上的正应力和剪应力；n、s分别为接触面上的法向应变和相对剪切应变；Dsn和Dns可以反映接触面法向和切向的耦合作用。弹塑性损伤矩阵中各元素的详细推导见文献6。该模型包括9个参数：Mohr-Coulomb强度参数0、r，接触面弹性剪切刚度Ds和弹性法向刚度Dn，损伤演化参数a、b，土的最大体应变参数0v、相对密度Dr以及土的临界状态体变参数K.模型参数可以根据接触面剪切试验成果确定。1.3 空间问题接触面单元有限元格式 三维问题中的接触面单元如图1所示。接触面单元的两个面可以是平面也可以是曲面，每个面由4个节点来定义。接触面单元为8节点等参单元，其几何插值公式为图1 空间问题接触面单元示意(10)式中：x是接触面上一点的整体坐标，xj是第j和j+4节点的节点坐标。插值形函数用自然坐标表示为：Nj=1/4(1+rjr)(1+sjs)(11)(j=1,2,3,4)式中：rj、sj为节点的自然坐标。如果将节点1到4组成的面定义为下表面，而节点5到8组成的面定义为上表面，则接触面单元两个面之间的位移可分别写成：ubot=N1u1+N2u2+N3u3+N4u4(12)utop=N1u5+N2u6+N3u7+N4u8(13)式中：uj为节点j的位移。进一步，将两个面位移转换到自然坐标系中，bot=aubot, top=autop(14)式中：a是整体坐标和自然坐标间的坐标转换矩阵。接触面单元两面间的相对位移为：v=top-bot,即：v=Bu(15)同样，8节点接触面单元也可以推广为其它更高精度的16节点、24四节点等参单元，更方便地与周围实体单元相结合或得到更精确的解答。接触面的应力-应变关系可以用下式表示。=D(16)式中：D为接触面单元本构矩阵。接触面单元厚度为t，相对位移向量可以用节点位移向量表示为v=t=Bu.可以得到接触面单元刚度矩阵为：(17)根据上述推导，本文编制三维有限元计算程序，对滑块问题作数值试验，论证接触面损伤模型的合理性和优越性，并对某围堰工程进行有限元计算，验证其应用于分析土与结构相互作用实际工程问题的可行性。2 三维混凝土滑块试验 置于混凝土基座上的滑块模型如图2所示。为了分析方便，设基座和滑块的容重均为零，基座固支，滑块自由，滑块与基座之间为中砂夹层，用空间接触面单元模拟。首先给滑块上表面施加100kPa的竖直方向初始应力，然后分级控制滑块的水平位移，直至接触面滑动破坏。接触面单元本构关系分别采用损伤模型和lough-Duncan非线性模型，模型参数分别见表1和表2.对损伤模型的计算结果与Clough-Duncan非线性模型相比较。不同模型计算得到的剪应力、相对法向位移和相对剪切位移的关系曲线分别如图3和图4所示。图2 滑块试验计算模型(单位：mm)表1 接触面损伤模型计算参数D50/mmDrDs/MPaDn/Mpaab0/ r/ K0v1.00.904.0400.03.171.6054.031.00.00510.22表2 接触面单元Clough-Duncan非线性弹性模型计算参数Cs/Mpas/ RfsKsns0.0043.30.48193931.00由滑块问题数值试验结果可以看出，接触面损伤模型可以很好地模拟土与结构物接触面和有夹层的岩体节理的力学与变形特性，并且可以描述应变软化和剪胀等岩土材料特有的现象，比以往的非线性弹性模型有较大的优越性。图3 剪应力和相对剪切位移关系曲线图4 相对法向位移和相对剪切位移关系曲线3 某围堰工程三维有限元分析某围堰工程全长1440m，高近90m，施工水深60m，堰顶高程88.5m.围堰采用砂砾石料、花岗岩风化砂及石碴混合料填筑，其防渗结构采用塑性混凝土防渗墙上接土工合成材料防渗心墙结构。其中河床深槽段为长约150m的双墙，墙厚1.0m，两道墙中心相距6.0m；位于基岩40m以上的两岸附近漫滩部分为单墙，墙厚0.8m.堰体深槽段典型剖面如图5所示。图5 某围堰深槽段典型剖面塑性混凝土防渗墙的应力变形是衡量围堰安危的重要指标。三维计算模型中考虑了单、双墙结合段材料的空间变化和防渗墙各槽段施工时间的差异，较为合理地模拟实际施工和运行工况。三维有限元计算共选取了21个断面，划分节点2460个，单元2123个。堰体有限元网格和防渗墙网格如图6所示。其中塑性混凝土防渗墙两侧与风化砂的接触面处设置接触面单元，采用损伤模型模拟。根据风化砂的三轴试验成果，确定损伤模型接触面单元计算参数如表3所示。表3 接触面单元损伤力学模型计算参数D50/mmDrDs/MPaDn/Mpaab0/ r/ K0v5.00.754.0400.03.171.6045350.0120.261围堰大规模抽水期和运行期是两种典型工况，抽水期上游水位为69m高程，下游水位为基岩面。运行期上游水位升至设计高水位85m高程。三维有限元计算分析得到各期的塑性混凝土上游防渗墙的水平位移、大主应力、小主应力和应力水平等值线分别如图7和图8所示。图6 围堰三维计算网格图7 围堰大规模抽水期塑性混凝土防渗墙计算成果由三维有限元计算成果可以看出，防渗墙体位移分布合理，墙体一般无拉应力，剪应力动用水平较低，在0.4以下。围堰运行期防渗墙体最大压应力和应力水平较小，受力状况好于围堰大规模抽水期，这是由于在运行期围堰堰体由湿容重变为浮容重，防渗墙所受拖曳力减小的缘故。围堰运行期防渗墙水平位移较大，这是由于上游水位升高、作用于防渗墙的水压力增大引起的。此结果与以往的分析成果6，8相一致。由于本问题防渗墙周围风化砂体积很大，接触面的应变软化和剪胀对整个结构的影响较小，因此损伤模型计算结果与传统非线性模型的分析成果6差别不大。上述计算可以看出，本文提出的损伤模型接触面单元能够方便地应用于实际工程问题的有限元分析，计算结果合理可靠。 4 结论图8 围堰运行期塑性混凝土防渗墙计算成果本文简要介绍了接触面损伤模型，推导了空间接触面单元的有限元格式，编制了三维岩土工程问题有限元计算程序，对土与结构相互作用问题进行了有限元计算分析。滑块问题数值模拟结果表明，损伤模型能够反映土与结构物接触面剪切过程中的应变软化和剪胀特性，较非线性弹性模型有显著的优越性。文中对一实际围堰工程进行了有限元分析，计算结果合理可靠，表明本文提出的损伤模型接触面单元可以方便地应用于实际工程问题的数值模拟。参 考 文 献：1 Goodman R E,Taylor R L,Brekke T L.A model for the mechanics of jointed rockJ。Journal of the Soil Mechanics and Foundations Division,ASCE,1968,94(3):637-659.2 Clough G W,Duncan J M.Finite element analysis of retaining wall behaviorJ。Journal of the Soil Mechanics and Foundations Division,ASCE,1971,97(12):1657-1672.3 Ghaboussi J,Wilson E L,Isenberg J.Finite element for rock joints and interfacesJ。Journal of the Soil Mechanics and Foundations Division,ASCE,1973,99(10):833-848.4 Desai C S,Rigby D B.Cyclic interface and joint shear device including pore pressure effectsJ。Journal of Geotechnical and Geo-en