基于统一强度模型的圆形硐室应力分析.docxVIP

基于统一强度模型的圆形硐室应力分析 1 基于abaqus的数值分析 中国力学专家俞茂宏从滑动位移横截面的力学模型出发，考虑到作用在两个剪切位移体上的所有应力分量对材料的屈服或破坏的不同影响，提出了一套可适用于各种岩石和材料的统一强度理论和统一形式的数学方案。Mohr-Coulomb强度理论和双剪强度理论均为其特例,并且还包含了可以比D-P准则更合理的新的计算准则,以及可以描述非凸极限面试验结果的新的非凸强度理论。同时统一强度理论可以很好地考虑中间主应力对岩体强度的影响,并能与岩土材料的真三轴试验结果相吻合,从而在岩土工程领域得到广泛的应用。通用有限元软件ABAQUS可以分析各种固体力学、结构力学,特别是能够处理材料非线性、几何非线性、接触非线性以及这三种非线性各种组合的高度非线性问题。ABAQUS本身带有较多材料非线性弹塑性模型,包括D-P模型、Mohr-Coulomb模型、剑桥模型等,但岩土工程领域中采用的统一弹塑性本构模型并没有在ABAQUS中得以实现,这使得ABAQUS在岩土工程数值分析中难以充分发挥,加入了统一弹塑性本构模型,则可以充分发挥ABAQUS软件的计算能力。为了弥补这一不足,本文利用用户材料子程序(UMAT)接口,通过Fortran或者VC编程来开发了统一弹塑性本构模型,对单轴压缩试验的数值模拟及圆形硐室问题进行了弹塑性分析,并与ABAQUS软件自带模型及解析解进行比较,以验证编制的接口程序的正确性。 2 弹塑性刚度矩阵及点 将统一强度理论表示成应力不变量的形式为: 式中:c0为内黏聚力;?为内摩擦角;b为中间主应力系数。 如图1所示,统一强度理论包含了多个强度理论:当0≤b≤1时为统一强度理论(外凸理论);当b0或b1时为双剪非凸强度理论(非凸理论);当b=0时统一强度理论退化为Mohr-Coulomb强度理论。本文研究仅考虑0≤b≤1的外凸理论。 由屈服函数F(F′)和塑性势函数Q(Q′)决定的弹塑性刚度矩阵可以写成: 式中:A与H′为硬化函数;为等效塑性应变。 F(F′)是统一强度理论的表达式。如果采用Q=F(Q′=F′),由式(5)可以得到相关联流动的统一强度弹塑性刚度矩阵。 文献给出了相关联流动的统一弹塑性流动矢量以及本构模型中奇异点的数学处理。定义统一强度理论的流动矢量为 同理,定义屈服面F′流动矢量为 如图1所示,统一弹塑性本构模型的屈服面应具有单一的硬化参数和流动方向,但在A、B、C点存在奇异性,定义A、B、C点为奇异点。奇异点是指流动矢量{a}在该点不能唯一确定的点。这里定义两类奇异点,针对不同的奇异点,采用不同的处理方法。 (1)对于在θ=θb,如图1中B点产生的奇异性,采用矢量平均的办法,即 (2)当b=1时,在点θ=0°和θ=60°处,如图1中A和C点产生的奇异性,采用数学极限的方法,确定流动矢量。 当θ=0°时, 当θ=60°时, 当b≠1时,在点θ=0°和θ=60°处产生的奇异,采用物理的方法,确定流动矢量。 当θ=0°时, 当θ=60°时, 3 比例修正因子确定 施加荷载增量后,首先计算试探应力: 将上述试探应力带入屈服条件(1),如果它不满足,表明此时材料的行为是弹性的。硬化参数A保持不变,第r个增量步的应力就等于试探应力。 如果满足屈服条件,则表示在当前荷载步内,积分点处达到塑性条件,因此,应按照塑性规律进行计算,即该点的应力应在屈服面上移动。 分以下两种情况: 第1种,如图2所示,应力状态由A点穿越屈服面达到B点。此时 为保持应力状态在屈服面上,则C点的位置,可由下式决定: 式中:R为比例修正因子 第2种,如图2所示,应力状态由屈服面A′点达到B′点。此时 如果荷载增量比较大,且应力点位于屈服面的大曲率附近,则上述过程仍存在较大的误差。为了提高计算的精度,可以根据超出屈服面的应力大小动态设定n个等分数,依次将迭代应力拉回到屈服面,最后再采用上述的比例修正法进行修正,就可得到较为准确的结果。如图3所示,将超过屈服面应力点分成n等分,则经过n个循环后,应力点将返回偏离屈服面的E点,然后再通过比例修正使其返回屈服面的E′点。 4 弹塑性本构模型求解过程 ABAQUS的用户材料子程序(UMAT)通过与其求解器Standard的接口来实现数据交流。UMAT有自己的书写格式与一些规范,与主程序共享的变量必须在子程序开头予以定义。而主程序通过ABAQUS输入文件(.inp)中的关键字“USER MATERIAL”来判断用户是否使用了自定义材料本构模型,从而扩展了它的适用性和应用空间。 统一弹塑性本构模型主要求解过程:每一个增量加载步开始时,ABAQUS主程序在单元的积分点上调用UMAT子程序,并传入应变增量、时间步长及荷载增量,同时也传入当前已知状态的应力、应变及其他与求解过程相关的变