基于三维变形体离散元法的材料受拉开裂行为研究

基于三维变形体离散元法的材料受拉开裂行为研究

在过去的20年里，科学家们进行了大量的数值模拟，并将混凝土作为具有中等困难的弯曲行为。以前的基于有限差分的数值方法可分为两种类型：分散裂缝模型和分离裂缝模型。在弥散裂缝模型中,裂缝以弥散的方式表示,无穷多的平行微小裂缝分布于整个有限单元中,裂缝通常以固定的有限单元网格模拟,裂缝的扩展以材料的刚度和强度的降低来表示,用于定义材料的应力应变本构关系呈现应变软化特性。当材料在受荷条件下存在大量微观裂缝以及随机分布的裂缝时,该方法是一种非常有效的方法。根据开裂面是否随缝端应力主轴的旋转而旋转,弥散模型又分为固定裂缝模型和旋转裂缝模型两种。分离裂缝模型基于位移场的不连续,通常采用非线性界面单元进行模拟。在计算中,当沿某一面或边的法向超过了材料的强度,则引入一界面单元,由于界面单元的引入使得该方法在有限元等数值方法中的实现面临下列问题:1)需要确定一种裂缝扩展准则;2)需要有效的网格自动重新剖分方法;3)用于转换新旧两套网格状态变量的网格映射技术。这些问题限制了分离裂缝模型在断裂力学中的应用,尤其是网格重新剖分降低了程序处理问题的能力和效率;然而,当可能开裂路径已经预知,我们也可以在分析前就引入分离裂缝而使问题变得简单。上述两类模型均在实践中得到一定的应用。此外,在断裂破坏的模拟方面,有限元方法由于其小位移假定使它在模拟裂缝的不稳定扩展直至系统崩溃等大位移问题时受到了限制,而离散单元法正是从模拟不连续介质的角度发展起来的,允许系统发生大的位移、转动甚至分离,如何将连续与离散这两类数值方法结合起来研究结构的断裂破坏过程,以期用于水利工程中坝体—地基系统破坏溃决过程的模拟正是本文的目标。迄今将离散单元法与断裂力学结合分析力学问题的工作仍很少。本文基于离散单元法结合弥散式旋转裂缝模型以及分离裂缝模型进行Ⅰ型断裂问题的研究,包括混凝土单轴拉伸的数值模拟以及三点弯曲梁的数值计算与试验结果对比,一方面验证了提出的离散元—断裂力学模型的可靠性与精度,另一方面为研究混凝土、岩石由断裂到破坏的连续—非连续系统过渡过程的模拟进行了有益的尝试。1开裂方向的应力应变力学模型Rashid于1968年首次提出弥散裂缝模型,而后Bazant基于弥散裂缝模型发展的断裂带理论实现了断裂过程的数值模拟,给出了Ⅰ型裂缝某点发生开裂后局部坐标下的应力应变全量关系。许多学者应用弥散式旋转裂缝模型在模拟混凝土开裂方面做了大量工作。该模型考虑当一点主应力分量达到抗拉强度时,材料在该方向发生开裂,在开裂方向的应力应变分量满足材料的软化力学行为,发生开裂后的不同软化阶段,裂缝方向始终与当前主应力方向一致,该模型假定混凝土总应变由材料自身应变εco和裂缝张开引起的应变εcr组成:以裂缝面法向为x轴建立局部坐标系,根据坐标转换关系有:式中N为坐标转换矩阵,ecr为局部坐标系下的开裂应变,s为同一局部坐标系下的应力,假设混凝土材料自身满足线弹性本构关系:在局部坐标系下裂缝应变与应力的关系为:Dcr为考虑裂缝受拉软化后在三个主应力方向上建立的解耦的开裂应变应力关系矩阵,综合上面几式,推导得各点有一个方向发生开裂时混凝土材料的应力~应变全量关系为:式中:si,ei(i=1,2,3)分别为三个主应力和主应变,E为弹性模量,ν为Poisson比,根据当前的应变状态确定开裂参数μ,以表征不同的开裂程度。采用旋转裂缝模型,依据Bazant断裂带理论,根据网格大小调整应力应变软化曲线的形状以保证材料的断裂能唯一,材料的开裂通过刚度和强度的降低来体现,进而可以获得分析系统的力学响应。2结构采后裂缝模型分离裂缝模型是由Hillerborg提出的虚拟裂缝模型(FictitiousCrackModel)发展而来的。学者们基于有限元、边界元方法应用该模型在模拟素混凝土以及钢筋混凝土的Ⅰ型和Ⅰ/Ⅱ混合型开裂上进行了广泛的研究,近来,Vahid基于分离裂缝模型和非正交弥散裂缝模型的耦合分析了混凝土拱坝的动力行为。分离式裂缝模型将准脆性材料的带状微裂区简化为一条分离的裂缝,在裂缝面上定义传递的应力与裂缝张开位移之间的关系,该模型作了如下基本假设:1)混凝土等材料的开裂存在一材料部分损伤的断裂过程区,沿开裂面传递的应力大小由相应裂缝的张开度w决定;2)断裂过程区外的材料仍为弹性变形区,当过程区的尖端应力超过抗拉强度后,过程区向前发展,当真实裂缝的尖端位移超过wf后,真实裂缝向前扩展;3)应力~裂缝张开度曲线下的面积即为材料的断裂能Gf,即Gf=∫0wfσ(w)dw。由此可见,分离裂缝模型包含了两个材料参数:断裂能Gf和抗拉强度ft,由于考虑了裂缝前缘的应变软化效应,引入了断裂力学的概念,采用了合理的裂缝扩展准则,能有效的描述混凝土结构中存在宏观可见裂缝的情况,使得该模型已广泛应用于静力、动力以及冲击荷载条件下混凝土开裂过程的模拟。在Ⅰ型裂缝情况下,垂直于裂缝面的应力~张开度关系曲线有多种形式,根据缝面初始刚度是有限的还是无限的,将这些曲线基本分为两类,即裂缝面本构关系是初始弹性的还是初始刚性的(图1a,图1b)。对于初始弹性的分离裂缝模型,首先随着裂缝张开度的增加,裂缝面传递的应力增加,一旦应力达到抗拉强度,则缝面传递的应力随着裂缝张开度的增加而逐渐下降到0,Xu和Needleman曾采用该本构关系研究脆性材料的动力开裂行为。对于初始是刚性的分离裂缝模型,裂缝面上初始刚度假定为无限大,只定义下降段裂缝面传递的应力与张开度的关系,Camacho和Ortiz曾应用该本构关系模拟脆性材料任意路径的冲击破坏行为。一旦材料的断裂能Gf和抗拉强度ft确定,便可以根据选取的软化曲线的形状确定裂缝最大张开度wf,本文假定缝面相比基质材料刚度而言无限大初始刚度进行分析。3裂缝端面形状与结构在弥散的旋转裂缝模型中,混凝土的开裂通过定义裂缝面的应力开裂应变的关系来实现,在分离裂缝模型中,混凝土的开裂通过定义裂缝面的应力和缝面张开位移的关系来实现,不同的软化曲线选择确定了裂缝尖端的形状,其反过来影响裂缝尖端附近的应力分布以及开裂区的大小,各种软化曲线的参数均由断裂能Gf和材料抗拉强度ft两个参数确定,最为常用的为线性软化曲线和Petersson提出的双线性软化曲线(图2)。由断裂带理论的断裂能准则,对于直线形式软化曲线,弥散裂缝模型临界开裂应变εf和分离裂缝模型中裂缝面应力为零时的裂缝张开度wf为:上式中h为由网格尺寸确定断裂带宽度,对于双线性软化曲线有:折点的位置为(2εf9,ft3)或(2wf9,ft3)。图中(wm,fm),(εm,fm)为曲线折点的位置。4离散单元的动力平衡方程三维可变形离散元将系统离散成可变形块体,块体边界作为不连续界面通过接触本构模型规定界面上的力和位移关系,块体内部差分为四面体常应变单元,该方法基于牛顿第二定律采用显式时步步进的中心差分法求解系统的运动,这一特点非常适合于渐进式破坏问题的模拟。计算的求解流程如图3所示,说明如下:(1)在每一时步初始时刻,首先对每一个节点由它所受到的合力F和初始运动状态通过动力平衡方程mu=F求解出它在下一个时步初始时刻的运动状态以及接触位移增量,此时加入运动边界条件,输入指定节点的位移、速度或加速度。(2)得到每一个节点新的运动状态后,根据应力应变本构关系求得块体内部应力,根据接触力—位移关系求解出接触力,以及其它需要考虑在内的各种力的变化。将各种力分配到节点上,再加入输入的力边界条件可以得到下一个时步初各个节点的合力,转入下一个时步的计算。在考虑质量阻尼和自重的条件下,对于离散单元中每一个差分节点,其动力平衡方程为:式中m为分配在节点上的质量,α为自适应阻尼的质量阻尼系数,F为节点受到除重力外的合力,包括块体间的接触力,积分到节点上的弹性力Fe以及施加的外载,对运动方程采用中心差分格式积分有:在变形体离散元中,块体进一步差分为四面体常应变单元,假定每一时步的变形为小变形,单元应变率为:i,j=1,2,3对应坐标x、y、z,对差分四面体单元应用高斯定理有:其中V为四面体的体积,S为四面体的外表面,nj为外表面的单位外法线向量分量。由于在每一个差分四面体单元内,应变率为常数,速度场为线性,每个面的外法线向量也是常数,由上式可推导得:设四面体节点编号由1到4,第l个顶点所对的面为第l个面,则上式中上标l表示节点l的变量,(l)表示面l的变量,求得速度的偏导数后,便可以求得应变率分量,相应的应变增量为,进而可以求出单元的总应变εij。当在离散元中实现弥散的旋转裂缝模型时,是在程序计算中修改基质材料的本构模型,一旦获得了单元的总应变εij,可以根据单元开裂的起裂条件,对于开裂单元根据公式(5)获得单元的应力。对于未开裂单元,采用常用的线弹性本构关系求解单元应力σij,最后由单元应力通过积分便可以获得等效到每一个节点上的弹性力Fe。可见弥散裂缝模型在离散元中是通过修改单元的材料本构模型来实现的,该方法对处理随机分布的裂缝的问题比较有效。在离散元中,如图4所示,块体与块体之间通过相互连接的弹簧和阻尼器传递相互之间的作用力。为了模拟结构的渐进性破坏过程,本文提出了分离裂缝模型与传统离散元法的耦合,通过块体间接触弹簧本构关系的修改来实现块体界面间的分离,进而模拟界面的裂缝扩展,而不是采用传统离散元中接触弹簧受拉脆断的方式。以线性软化关系为例,在受拉软化阶段,弹簧上传递的法向接触力nF与裂缝面上相对法向张开度之间的关系为:式中:ft为离散缝面的抗拉强度,Ac为接触面积,w为裂缝张开度,wf为由公式(6)确定的最大裂缝张开度。基于离散单元法采用分离裂缝模型实现开裂过程的模拟,其最大优点在于不仅可以获得裂缝稳定扩展阶段的破坏过程全曲线,而且可以获得材料部分或完全破坏后,裂缝扩展过程失稳时的整个溃决过程,实现全过程的仿真分析。5计算与分析5.1采用本构模型的方法计算选取尺寸为0.5m×0.5m×0.5m的立方块体,在块体底部施加竖向滚轴约束,块体顶部做匀速直线运动,速度为1.0×10-4m/s,两种开裂模型计算采用的网格剖分如图5所示。混凝土材料参数为:弹性模量E=20GPa,Poisson比v=0.2,抗拉强度ft=2.0MPa,断裂能Gf=100N/m。当采用分离裂缝模型时,假定混凝土材料为线弹性,块体中部存在一预先设定的分离裂缝。应用两种开裂模型分别用线性软化和双线性软化曲线来模拟块体断裂过程。根据断裂能的定义,荷载位移曲线所围的面积即块体断裂所需要的能量应为:根据软化曲线各个参数的定义,在混凝土单轴拉伸情况下,计算得的荷载位移曲线、裂缝面法向应力与开裂应变以及裂缝面法向应力与裂缝张开度曲线如图6所示,根据计算结果,可以得出如下结论:(1)从图中可以看出,根据断裂带理论确定选取的软化曲线的软化参数,不同开裂模型,不同软化曲线的计算结果与理论解相同,说明程序的实现是正确的。(2)采用线性与双线性软化曲线模型,其荷载位移曲线的形状不同,但其所围成的面积相同,即混凝土发生完全开裂所需要的能量相同,这与软化参数依据确定的断裂能获得相吻合。不同的开裂模型,只要采用的软化曲线相同,计算得的荷载位移曲线相同,说明在保证断裂能守恒的条件下,两种开裂模型在模拟混凝土Ⅰ型开裂问题上是等价的。5.2模拟结果与分析选取Petersson三点弯梁试验,梁的尺寸为2000mm×200mm×50mm,切口深度为半梁高,力学参数为:混凝土弹性模量E=30GPa,Poisson比v=0.2,抗拉强度ft=3.33MPa,断裂能Gf=124N/m,采用旋转裂缝模型计算时,梁切口宽度为20mm,计算中采用位移加载控制,切口上方对应点做匀速直线运动,速度为2.0×10-4m/s,两种模型网格划分情况如图7所示。采用前述的旋转裂缝模型和分离裂缝模型,分别采用线性软化、双线性软化曲线对三点弯梁试验进行模拟,依据断裂能准则确定相应软化曲线的各个参数,计算得施加的总荷载F与切口顶部节点的位移D关系曲线如图8所示,基于计算结果,可以看出:(1)对比试验结果的上下包络线与数值模拟结果,可知两种模型数值计算都与试验结果符合较好。采用线性软化曲线计算的峰值荷载稍偏大,采用双线性软化模型比线性软化模型计算的峰值荷载前段偏小,而后段则偏大,这是由于在断裂能守恒的条件下,双线性软化曲线在拐点以前的下降段较线性软化曲线陡,在拐点后下降段相比较缓所致。(2)对比弥散裂缝模型和分离裂缝模型计算结果(图8c),发现分离裂缝模型获得的峰值荷载略偏大,但随着荷载的施加,后期承载力下降的快,比弥散裂缝模型低,这是由于在弥散裂缝模型中,当开裂方向与网格线不平行时,在弥散表征裂缝的单元内随着裂缝扩展会产生虚假剪切应力,呈现出剪切闭锁现象,高估了结构的残余承载能力,而分离裂缝模型在完全对称的几何与荷载条件下,裂缝面的传递的剪切应力为零,残余承载能力较小。(3)为了保证计算结果的客观性,在保证断裂能守恒的条件下,弥散裂缝模型需要根据网格大小来调整断裂带宽度,以消除计算结果对网格尺寸的依赖性;分离裂缝模型不存在这一问题,只需在分离裂缝面上多布置一些非线性接触弹簧以较好地模拟裂缝的扩展。(4)图9给出了采用双线性软化分离裂缝模型,加载点位移D为不同值,系统变形放大100倍的形态图,可以看到,当继续加载时,预设的分离裂缝将完全脱开,系统失稳,基于离散单元法的分离式裂缝模型的实现使得跟踪系统从线弹性阶段到完全破坏阶段的模拟成为可能。6破坏阶段的荷载位移全曲线(1)本文基于变形体离散元方法分别采用弥散式旋转裂缝和分离式裂缝两种模型分析了混凝土、岩石等准脆性材料的受拉开裂行为。结果表明:模型能够合理预测裂缝的起裂和扩展,获得了分析系统在外载作用下从线弹性阶段到破坏阶段的