基于变形体离散元方法的裂缝模型及变形分析

基于变形体离散元方法的裂缝模型及变形分析

混凝土结构的负荷影响下的损害和破坏，岩质边坡的变形和稳定，以及不连续的地下孔室围岩的稳定性是当前防灾和民防工程中必须解决的问题。这些问题存在一个显著的特点,系统在外载作用下从连续或断续状态逐渐转变为非连续状态而形成独立块体,系统的力学行为最终由独立块体间的相互作用决定。采用数值方法描述系统的力学行为,实现连续介质的损伤断裂并向非连续介质的转化以致失稳是当今计算固体力学迫待研究的课题。Rashid于1968年提出弥散裂缝模型模拟材料开裂行为,Bazant发展的断裂带理论实现了断裂过程的数值模拟,Rots采用固定裂缝模型和Jirasek采用旋转裂缝模型研究了混凝土的Ⅰ型开裂,但以有限元为基础的连续介质力学方法限制于小位移问题,不能处理大变形失稳问题,基于非连续介质的离散元法在模拟连续介质损伤断裂力学方面研究较少。本文将准脆性材料的开裂模型嵌入到离散单元法中,提出一种将损伤断裂模型与变形离散元结合的方法以研究连续-非连续系统的过渡过程,文中首先讨论材料开裂模型,不连续界面的引入以及裂缝表面的接触模拟,离散块体的运动求解等问题,文末以三点弯梁断裂问题作数值解与试验比较。1三维非线性衰减理论模型依据裂缝方向的不同处理方法,基于断裂带理论的本构关系,将断裂力学模型分为固定裂缝和旋转裂缝两种。1.1土体裂裂变形模型固定裂缝模型模拟材料的一应力点的开裂方向由该点在加载过程中达到抗拉强度时主应力的方向所决定,在后续计算中不再发生改变。考虑到裂缝面上由于骨料互琐而产生的残余剪切刚度,通常引入一随着裂缝的张开而逐渐减小或常剪切抗力系数,采用应变分解方法,认为混凝土总应变ε由混凝土材料自身应变εco和裂缝张开引起的应变εcr组成,即ε=εco+εcr.(1)以裂缝面法向为x轴建立局部坐标系,根据坐标转换关系有:εcr=Νecr‚(2)s=ΝΤσ.(3)式中:N为坐标转换矩阵,开裂应变ecr包含垂直裂缝面方向的应变分量及与裂缝面平行的剪应变分量,s为同一局部坐标系下的应力。假设混凝土材料自身满足线弹性本构关系,在局部坐标系下裂缝应变与应力的关系为s=Dcrecr‚(4)式中Dcr为综合考虑裂缝拉伸和剪切软化后构成的3×3矩阵。综合式(1)—(4),可得开裂混凝土材料的应力-应变全量关系为σ=[Dco-DcoΝ(Dcr+ΝΤDcoΝ)-1ΝΤDco]ε‚(5)式中Dco为混凝土弹性刚度矩阵。固定裂缝模型裂缝方向由达到抗拉强度时的主应力方向决定,在后续的计算中不发生改变,随着裂缝的扩展,一点的材料主轴可能发生偏转,因此采用固定裂缝模型计算一点的力学行为通常与材料发生开裂后的实际力学特征不符。其后,学者们发展了固定裂缝模型,定义了一个开裂后连续旋转的应力主轴,在旋转的主轴空间中建立应力应变关系,提出了旋转裂缝模型,并应用该模型在模拟混凝土开裂方面做了大量工作。1.2局部应力应变动力系统的求解旋转裂缝模型考虑当一点主应力分量达到抗拉强度时,材料发生开裂,在开裂方向的应力应变分量满足材料的软化力学行为,发生开裂后的不同软化阶段,裂缝方向始终与当前主应力方向垂直,只需定义主轴空间中的应力应变关系,即:εcr=Τecr‚(6)ecr=gs.(7)式中:ecr为局部坐标系下3个主方向的开裂应变,s为局部坐标系下的应力,g为裂缝特性矩阵,T为坐标转换矩阵。根据旋转裂缝模型的假设,在求解过程中需要不断更新T矩阵,使得局部坐标系随主应变方向变化而旋转。以各点有一个主方向发生开裂为例:ecr1=g1s1‚(8)g1=(1-μ)μE‚(0＜μ≤1).(9)式中μ为材料割线模量相对于线弹性模量E的折减系数。综合式(1)、(6)—(9),假定混凝土材料自身满足线弹性本构关系,在由主应力方向建立的局部坐标系中应力应变全量关系为:{s1s2s3}=E1-ν-2μν2[μ(1-ν)μνμνμν1-μν21+νν(1+μν)1+νμνν(1+μν)1+ν1-μν21+ν]{e1e2e3}‚(10)式中:si、ei(i=1,2,3)分别为3个主应力和主应变,E为弹性模量,ν为Poisson比。根据当前的应变状态确定折减系数μ,以表征不同的开裂程度。1.3断裂带及材料参数对力学性能的影响在复杂应力路径作用下,准脆性材料的开裂行为趋向于以Ⅰ型裂纹为主的方向。因此,混凝土的张拉软化曲线成为上述裂缝模型的基础,依据断裂带理论,遵循断裂能Gf唯一准则,根据材料抗拉强度ft和基于网格尺寸确定的断裂带宽度h,调整材料应变软化曲线参数以准确描述材料的断裂能释放,通常采用线性和双线性软化曲线描述准脆性材料的力学响应。两种软化曲线的临界开裂应变为:εf1=2Gf/fth‚(11)εf2=3.6Gf/fth.(12)对于双线性软化曲线,软化段中折点位置由式(13)和式(14)确定:εm=2εf2/9‚(13)fm=ft/3.(14)2局部阻尼分析方法文提出的离散单元法将为不连续面所切割的系统差分为四面体常应变单元以模拟块体的变形,不连续界面通过接缝本构模型规定界面上的力和位移关系,作用力被集中在差分节点上。该模型基于Newton第二定律,采用动态松弛法求解节点运动方程,使用显式中心差分方法求解节点的位移和单元应变,各节点的动力平衡方程为m⋅⋅ui=∫sσijnjds+Fexti‚(15)式中:m为分配在节点上的质量,Fexti为节点所受的外力,s为包围积分质点域的外边界,nj为边界的单位外法向分量。由于采用动力松弛解法,对于稳态问题,通常采用粘性阻尼以消散系统的能量。对于有荷载突变或系统逐渐破坏的情况,采用局部阻尼,阻尼力与节点的不平衡力成正比,即Fdi=α|∑F(t)i|sgn(˙u(t-Δt/2)i)‚(16)式中:Fdi为节点阻尼力,∑F(t)i为除阻尼力外节点所受的合力,α为小于1的常数。对运动方程采用差分格式得:˙u(t+Δt/2)i=˙u(t-Δt/2)i+(∑F(t)i-Fdi)Δtm‚(17)u(t+Δt)i=uti+˙u(t+Δt/2)iΔt.(18)对于已有的不连续面以及材料开裂后形成的新的不连续面,通过法向和切向的弹簧传递接触面之间的相互作用力,接触模型法向采用拉断准则,切向遵循Mohr-Coulomb准则:ΔFn=-ΚnΔunAc‚(19a)ΔFs=-ΚsΔusAc‚(19b)Fn=0‚Fs=0‚当Fn＞ftAc‚(20a)Fs=Fs(f|Fn|+CAc)/|Fs|‚当|Fs|＞f|Fn|+CAc‚(20b)式中:Kn、Ks是接触法向和切向弹簧的刚度,Δun、Δus为相对位移增量的法向和切向分量,ΔFn、ΔFs为法向和切向弹簧力变化量,Fn、Fs为切向和法向弹簧力的大小,Ac为子接触面积,f为接触面的摩擦因数,C为粘聚力。可变形离散元能够很好模拟不连续介质的静动力力学行为,由于其采用增量解法,可以较容易将弥散裂缝模型引入到材料本构模型中以模拟材料开裂行为。3材料的拉拔试验对于准脆性材料,定义开裂系数Ff即沿裂缝法向的非弹性开裂应变εcrn和最大开裂应变εf的比来表征开裂程度,即Ff=εcrn/εf.(21)当Ff为0时,材料处于线弹性阶段,当0<Ff<1时,材料处于软化阶段,当Ff=1时,沿裂缝面的抗拉强度降为零,此时通过常应变四面体单元形心插入一离散裂缝,在该离散裂缝面上形成新的接触,接触面的抗拉强度和粘聚力均为0,f采用本体材料的摩擦因数,接触的法向和切向仍遵循Mohr-Coulomb准则。对于固定裂缝模型,离散裂缝方向垂直于裂纹起裂时的最大主应力方向,对于旋转裂缝模型,离散裂缝方向垂直于Ff=1时的最大主应力方向,而忽略开裂过程中裂缝方向的旋转。为此,根据开裂系数Ff的大小可以确定弥散裂缝张开的程度以及离散裂缝插入的顺序,最终形成非连续或断续系统。4网格划分的确定选取Peterson三点弯梁试验,梁的尺寸为2000mm×200mm×50mm,切口深度为半梁高,力学计算参数为:混凝土弹性模量E=30GPa,Poisson比ν=0.2,抗拉强度ft=3.33MPa,断裂能Gf=124N/m。计算采用两套网格,其主要差别在于切口宽度和切口处的网格划分,细网格的单元尺寸为粗网格(切口宽度20mm)的单元尺寸的一半,从正面看粗网格划分情况如图1所示。本文采用三维固定和旋转裂缝两种模型,分别应用线性、双线性软化曲线计算施加的总荷载F与切口顶部节点的位移D关系曲线如图2—图6所示。5结构承载能力分析依据图2—图5,按照Bazant断裂带理论,每个模型中不同粗细的网格剖分计算得到的荷载位移关系曲线基本相同,说明在非线性断裂力学模型中,当计算遵循断裂能守恒准则,根据数值模拟的断裂带宽度的大小来调整应变软化参数,可以克服网格剖分尺寸对计算结果的影响,同时也说明了细网格并不能克服由于网格线与开裂方向不一致而带来的剪切锁死现象,即目前的局部化模型存在一定程度的网格方向依赖性。对比试验结果的上下包络线与数值模拟结果,可知采用线性软化曲线较双线性软化曲线计算的峰值荷载偏大,荷载位移曲线下降较陡,双线性软化模型后期承载力偏大,这是由于对于同一断裂带宽度,在保证断裂能守恒的情况下,双线性软化曲线在拐点以前的下降段较线性软化曲线陡,在拐点后下降段相比较缓所致。从图6中可以看出,在裂缝发生的初始阶段,固定裂缝模型和旋转裂缝模型的计算结果基本相同,而后旋转裂缝模型计算的残余承载能力偏小,更接近于试验值。分析原因知随着裂缝扩展,裂缝应变在总应变分布中占主导地位,由于开裂方向与网格线的不平行造成弥散表征裂缝的单元内产生了虚假的剪切应力,裂缝单元主应力方向发生一定程度的偏转。采用固定裂缝模型使得该虚假剪应力的作用随着裂缝的进一步扩展而增强,系统刚度较实际情况偏大;采用旋转裂缝模型由于裂缝法向始终与主应力方向共轴,一定程度上避免了主应力方向偏离初始开裂方向而带来的剪切锁死现象,计算的结构承载能力更符合实际情况。图7给出了粗网格划分情况下不同断裂模型的裂缝开裂方向,旋转裂缝模型与固定裂缝模型相比开裂方向与梁所承受的对称荷载更为相符。6支护结构的变形1)本文采用非线性断裂力学模型,考虑混凝土、岩石等准脆性材料的受拉开裂,基于变形体离散元方法分析了材料受拉破坏行为。模型能够定量预测初始裂纹时的应力状态及开裂过程中的变形,获得了合理的结构荷载位移响应曲线,为连续-非