《接触面损伤模型应用于三维有限元分析》由会员分享，可在线阅读，更多相关《接触面损伤模型应用于三维有限元分析（8页珍藏版）》请在金锄头文库上搜索。
1、接触面损伤模型应用于三维有限元分析胡黎明，濮家骝，王 刚(清华大学 水利水电工程系)摘要：本文简要介绍了作者提出的接触面损伤模型，推导了接触面单元的三维有限元格式。将损伤模型接触面单元应用于三维有限元计算，对滑块问题进行了数值试验，比较了不同接触面模型的计算结果，验证了接触面损伤模型分析土-结构相互作用问题的合理性和优越性。文中对一围堰工程进行了有限元分析，计算结果合理可靠，表明本文提出的损伤模型接触面单元可以方便地应用于实际工程问题的数值模拟。关键词：土与结构物接触面；应变软化；损伤模型；有限元分析作者简介：胡黎明(1974-)，男，山西文水人，清华大学讲师，主要从事岩土工程研究。不同材料接
2、触面的力学特性及其数值模拟是岩土工程长期研究的重要课题。Goodman1对岩石节理裂隙的特性进行了研究，认为可以用 4 节点无厚度节理单元来模拟节理的滑移和张裂等变形特性。Clough 和 Duncan 提出非线性弹性模型描述应力与相对位移之间的关系 2。Ghaboussi 提出了基于弹性理论分析岩石节理相互作用面的单元，但由于采用相对自由度作为基本未知量，实际应用不方便，不适合推广使用 3。Desai 认为两种材料接触面剪应力传递和剪切带的形成均发生在接触面附近的一薄层土体中，从而提出薄层单元的概念，用于模拟土与结构物接触面及岩石节理的粘接、滑动、张开和闭合等各种接触状态4。殷宗泽提出接触面
3、单元的变形由两个部分组成，一部分是土体的基本变形，其变形特性与周围土体单元相同；另一部分是接触面上的破坏变形，包括滑动破坏和拉裂破坏，总的变形是两者的叠加 5。胡黎明建立了一个能够反映接触面变形机理和力学特性的接触面损伤模型，可以模拟其剪切过程中的应变软化和剪胀现象 6。本文将损伤模型应用于有限元计算，验证其分析土与结构相互作用问题的合理性和有效性。1 损伤模型接触面单元1.1 接触面非线性弹性模型 以往的接触面本构关系不考虑切向和法向的耦合作用，接触面应力和相对位移关系为：(1)Clough 和 Duncan 提出的接触面非线性弹性模型为：当接触面受压，n0 时Kss=Ks ( n/Pa)n
4、s(1-Rfs/ ntan s+cs)2, Knn=106MPa (2)当接触面受拉，n0 时Knn=Kss=0.01MPa (3)接触面模型参数有 5 个：剪切劲度参数 Ks、n s，接触面强度参数 cs、 s和破坏比 Rfs。1.2 接触面损伤模型 胡黎明进行了系统的土与结构物接触面的直剪试验，研究了接触面力学性质和变形机理 6，7 。试验表明，粗糙接触面剪切破坏时在接触面附近产生应变局部化形成剪切带，并伴有应变软化和剪胀现象。剪切带厚度t5D 50。剪切应变和剪胀均发生在剪切带内。根据试验结果建立了描述接触面应变软化和剪胀现象的损伤力学模型，认为某一状态接触面由无损伤状态和临界状态组合而
5、成，损伤演化过程为无损伤状态向临界状态的逐渐转化过程。1.2.1 无损伤状态 无损伤状态接触面采用弹塑性本构模型来描述，其屈服函数为F= i2/Ds+ ps( i-sign( i)5tg 0)=0 (4)式中： 0为摩擦角， ps为塑性剪切应变，D s为剪切弹性模量，上标 i 表示无损伤状态。采用不相适应流动法则，设塑性势函数 Q= i，即法向无塑性变形，只有弹性变形。采用接触面塑性剪切应变为硬化参数，即 H= ps.1.2.2 临界状态 临界状态接触面可用理想塑性模型来描述。其应力-应变关系为(5)式中： r为残余摩擦角； cv为临界状态体积应变； 0v为土的最大体应变；K 为体积应变参数，
6、与粗糙度 R 有关，与正应力 无关；D r为土体初始相对密度， 0为常数，一般取为 10kPa；上标 c 表示临界状态。1.2.3 损伤状态变量演化律 土与结构物接触面剪切应变是造成损伤的直接原因。假设损伤演化过程只与接触面塑性剪切应变有关，损伤状态变量为D=1-exp(-a ps b) (6)式中：a、b 为与 R 相关的参数，与正应力无关。根据损伤力学的基本原理，接触面应力-应变状态变量为无损伤和临界状态的线性叠加，即(7)对接触面应力-应变状态变量之间的关系式取微分，可推导得到增量应力-应变关系：d=D epd (8)接触面应力-应变关系矩阵，即弹塑性损伤矩阵为：(9)式中：， 分别为接
7、触面上的正应力和剪应力； n、 s分别为接触面上的法向应变和相对剪切应变；D sn和 Dns可以反映接触面法向和切向的耦合作用。弹塑性损伤矩阵中各元素的详细推导见文献6。该模型包括 9 个参数：Mohr-Coulomb 强度参数 0、 r，接触面弹性剪切刚度 Ds和弹性法向刚度Dn，损伤演化参数 a、b，土的最大体应变参数 0v、相对密度 Dr 以及土的临界状态体变参数 K.模型参数可以根据接触面剪切试验成果确定。1.3 空间问题接触面单元有限元格式 三维问题中的接触面单元如图 1 所示。接触面单元的两个面可以是平面也可以是曲面，每个面由 4 个节点来定义。图 1 空间问题接触面单元示意接触面
8、单元为 8 节点等参单元，其几何插值公式为(10)式中：x 是接触面上一点的整体坐标，x j是第 j 和 j+4 节点的节点坐标。插值形函数用自然坐标表示为：Nj=1/4(1+rjr)(1+sjs)(j=1,2,3,4) (11)式中：r j、s j为节点的自然坐标。如果将节点 1 到 4 组成的面定义为下表面，而节点 5 到 8 组成的面定义为上表面，则接触面单元两个面之间的位移可分别写成：ubot=N1u1+N2u2+N3u3+N4u4 (12)utop=N1u5+N2u6+N3u7+N4u8 (13)式中：u j为节点 j 的位移。进一步，将两个面位移转换到自然坐标系中，bot=aubo
9、t, top=autop (14)式中：a 是整体坐标和自然坐标间的坐标转换矩阵。接触面单元两面间的相对位移为：v= top- bot,即：v=Bu (15)同样，8 节点接触面单元也可以推广为其它更高精度的 16 节点、24 四节点等参单元，更方便地与周围实体单元相结合或得到更精确的解答。接触面的应力-应变关系可以用下式表示。=D (16)式中：D 为接触面单元本构矩阵。接触面单元厚度为 t，相对位移向量可以用节点位移向量表示为v=t=Bu.可以得到接触面单元刚度矩阵为：(17)根据上述推导，本文编制三维有限元计算程序，对滑块问题作数值试验，论证接触面损伤模型的合理性和优越性，并对某围堰工程
10、进行有限元计算，验证其应用于分析土与结构相互作用实际工程问题的可行性。2 三维混凝土滑块试验 置于混凝土基座上的滑块模型如图 2 所示。为了分析方便，设基座和滑块的容重均为零，基座固支，滑块自由，滑块与基座之间为中砂夹层，用空间接触面单元模拟。首先给滑块上表面施加 100kPa 的竖直方向初始应力，然后分级控制滑块的水平位移，直至接触面滑动破坏。接触面单元本构关系分别采用损伤模型和lough-Duncan 非线性模型，模型参数分别见表1 和表 2.对损伤模型的计算结果与 Clough-Duncan 非线性模型相比较。不同模型计算得到的剪应力、相对法向位移和相对剪切位移的关系曲线分别如图 3 和
11、图 4 所示。图 2 滑块试验计算模型(单位：mm)表 1 接触面损伤模型计算参数D50/mm Dr Ds/MPa Dn/Mpa a b 0/ r/ K 0v1.0 0.90 4.0 400.0 3.17 1.60 54.0 31.0 0.0051 0.22表 2 接触面单元 Clough-Duncan 非线性弹性模型计算参数Cs/Mpa s/ Rfs Ks ns0.00 43.3 0.48 19393 1.00由滑块问题数值试验结果可以看出，接触面损伤模型可以很好地模拟土与结构物接触面和有夹层的岩体节理的力学与变形特性，并且可以描述应变软化和剪胀等岩土材料特有的现象，比以往的非线性弹性模型有
12、较大的优越性。图 3 剪应力和相对剪切位移关系曲线 图 4 相对法向位移和相对剪切位移关系曲线3 某围堰工程三维有限元分析某围堰工程全长 1440m，高近 90m，施工水深 60m，堰顶高程 88.5m.围堰采用砂砾石料、花岗岩风化砂及石碴混合料填筑，其防渗结构采用塑性混凝土防渗墙上接土工合成材料防渗心墙结构。其中河床深槽段为长约 150m 的双墙，墙厚 1.0m，两道墙中心相距 6.0m；位于基岩 40m 以上的两岸附近漫滩部分为单墙，墙厚 0.8m.堰体深槽段典型剖面如图 5 所示。图 5 某围堰深槽段典型剖面塑性混凝土防渗墙的应力变形是衡量围堰安危的重要指标。三维计算模型中考虑了单、双墙
13、结合段材料的空间变化和防渗墙各槽段施工时间的差异，较为合理地模拟实际施工和运行工况。三维有限元计算共选取了 21 个断面，划分节点 2460 个，单元 2123 个。堰体有限元网格和防渗墙网格如图 6 所示。其中塑性混凝土防渗墙两侧与风化砂的接触面处设置接触面单元，采用损伤模型模拟。根据风化砂的三轴试验成果，确定损伤模型接触面单元计算参数如表 3 所示。表 3 接触面单元损伤力学模型计算参数D50/mm Dr Ds/MPa Dn/Mpa a b 0/ r/ K 0v5.0 0.75 4.0 400.0 3.17 1.60 45 35 0.012 0.261围堰大规模抽水期和运行期是两种典型工况
14、，抽水期上游水位为 69m 高程，下游水位为基岩面。运行期上游水位升至设计高水位 85m 高程。三维有限元计算分析得到各期的塑性混凝土上游防渗墙的水平位移、大主应力、小主应力和应力水平等值线分别如图 7 和图 8 所示。图 6 围堰三维计算网格 图 7 围堰大规模抽水期塑性混凝土防渗墙计算成 果由三维有限元计算成果可以看出，防渗墙体位移分布合理，墙体一般无拉应力，剪应力动用水平较低，在0.4 以下。围堰运行期防渗墙体最大压应力和应力水平较小，受力状况好于围堰大规模抽水期，这是由于在运行期围堰堰体由湿容重变为浮容重，防渗墙所受拖曳力减小的缘故。围堰运行期防渗墙水平位移较大，这是由于上游水位升高、
15、图 8 围堰运行期塑性混凝土防渗墙计算成果作用于防渗墙的水压力增大引起的。此结果与以往的分析成果 6，8 相一致。由于本问题防渗墙周围风化砂体积很大，接触面的应变软化和剪胀对整个结构的影响较小，因此损伤模型计算结果与传统非线性模型的分析成果 6差别不大。上述计算可以看出，本文提出的损伤模型接触面单元能够方便地应用于实际工程问题的有限元分析，计算结果合理可靠。4 结论本文简要介绍了接触面损伤模型，推导了空间接触面单元的有限元格式，编制了三维岩土工程问题有限元计算程序，对土与结构相互作用问题进行了有限元计算分析。滑块问题数值模拟结果表明，损伤模型能够反映土与结构物接触面剪切过程中的应变软化和剪胀特
16、性，较非线性弹性模型有显著的优越性。文中对一实际围堰工程进行了有限元分析，计算结果合理可靠，表明本文提出的损伤模型接触面单元可以方便地应用于实际工程问题的数值模拟。参 考 文 献：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
