基于ABAQUS 的耐火材料细观损伤本构模型的开发与实现

基于ABAQUS的耐火材料细观损伤本构模型的开发与实现基于ABAQUS的耐火材料细观损伤本构模型的开发与实现陈瑶，王志刚，刘昌明基金项目：博士点基金（20214219110003）作者简介：陈瑶，（1988-），男，硕士研究生。主要从事结构设计与仿真和软件的二次开发。通信联系人：王志刚，（1973-），男，博士，教授，博士生导师。主要从事细观力学在耐火材料设计中的应用、金属与非金属复合结构热机械行为研究、机械设备故障诊断、结构强度分析与计算机仿真等方面的研究。E-mail:wzgwy@（武汉科技大学机械自动化学院，武汉430081）5摘要：针对耐火材料复杂的非线性力学行为及其损伤机理，在细观力学基础上提出了基于界面相模型的耐火材料细观损伤力学模型，利用ABAQUS软件的二次开发平台，开发了该本构关系的子程序.利用该本构模型，对镁碳质耐火材料的受压试验进行了数值模拟.结果表明，运用该本构模型对耐火材料的非线性损伤力学行为进行模拟与试验结果能够较好的吻合。关键词：耐火材料；本构模型；细观力学；二次开发；ABAQUS10中图分类号：TB35SecondaryDevelopmentofMicromechanicsDamageConstitutiveModelofRefractorybasedonABAQUSCHENYao,WANGZhigang,LIUChangming15(MachineryandAutomationSchool,WuhanUniversityofScienceandTechnology,WuHan430081)Abstract:Accordingtothecomplexnonlinearmechanicsbehaviorandthedamagemechanismsofrefractorymaterials,damageconstitutivemodelofrefractorymaterialswasputforwardbasedontheinterphasemodelbythemethodofmicromechanics.Utilizingthemodelproposed,a20subroutinewasdevelopedthroughthecommercialsoftwareABAQUS.Usingthesubroutine,thedamageofthemagnesiumcarbonrefractorywassimulatedundercompression.Theresultsofthesimulationareingoodagreementwiththeexperimentaldata.Keywords:refractories;constitutivemodel;micromechanics;secondarydevelopment;ABAQUS250引言耐火材料作为高温容器件的隔热保温材料，是保障高温设备安全运行的关键.热应力破坏是耐火材料主要损毁方式之一，数值模拟是研究耐火材料热应力的主要技术手段.本构模型能否真实反映其力学行为，将直接决定有限元法等数值模拟技术计算的准确性.然而耐火材料的微观非均质性使其物理性能，特别是热机械性能受温度和载荷状况的影响较大，在宏30观上表现为依赖于温度的非线性力学行为.由于耐火材料力学行为的非线性源可能来自任何可能的微观变形和损伤机理，因此，建立在连续介质假设基础上的材料本构行为的统一描述就变得非常困难.到目前为止，仍未找到一种适合于数值模拟的耐火材料力学行为的表征方法。耐火材料复杂的非线性力学行为，实际上是其微观结构与外载荷相互作用的结果.任何35材料在使用中都不可避免地会产生微损伤，而这些微损伤的存在必将对材料的力学性能产生重要的影响.因此,从材料的显微结构出发,利用细观力学方法,研究其微观结构与外载荷的相互作用及其对宏观性能的影响，对探明耐火材料非线性力学行为的产生机制具有重要的意义.目前,国内外学者利用细观力学理论研究复合材料的损伤已经取得了一些进展.他们提40出了一系列基于细观力学方法的耐火材料本构关系，为进一步开发适用于耐火材料的本构关系模型奠定了基础.本文针对细观力学方法提出的基于界面相模型的耐火材料损伤本构关系模型，采用ABAQUS有限元软件提供的二次开发平台，实现了该损伤本构关系的数值模拟.利用数值算法，对耐火材料的受压试验进行了数值模拟，对比分析解析结果和实验结果，验证了该本构45关系的有效性.1耐火材料本构模型1.1耐火材料本构模型概述材料本构方程是用来描述其应力-纵向应变关系的数学表达式，一般来说本构方程应该尽可能地简单，且需要的力学参数也要尽可能的少，但必须包括造成非线性行为的原因，如50损伤等.目前，在已有的耐火材料本构模型里，热弹性模型[1]最为简单，它是一种近似模型，没有考虑任何非线性力学行为.等人[2]提出的弹塑性模型、J.M.Robi等人[3]提出的损伤弹塑性模型和Hemandez等人[4]提出的考虑脱水的弹塑性模型，是目前使用最广泛的一种模型.但模型涉及到许多需要实验测定的材料参数和只考虑了受压时的材料特征.弥散裂纹模型[5]将拉伸状态下的裂纹扩展行为与压缩的塑性行为统一起来，较为全面地描述耐火材料的55非线性行为.MitsuoSugawara等人[6]提出的神经网络模型，只是本构方程的一种近似算法.实际上，对耐火材料合适的非线性行为分析将取决于材料的细观结构及其与载荷之间的相互作用.因此，从微观尺度出发，利用细观力学方法，研究耐火材料宏观非线性行为的细观力学机理，不失为是解决耐火材料复杂非线性行为的另一种有效途径[7].目前,国内外学者利用细观力学理论研究复合材料的损伤已经取得了一些进展.H.K.Lee和S.H.Pyo[8]在60Eshelby细观力学方法的基础上,利用Weibull概率函数来表征颗粒和基质之间的弱界面变化概率.Schjudt-Thomsen和Pyrz[9]利用Mori-Tanaka平均场理论以及改进的Eshelby张量来研究具有弱界面的短纤维复合材料.Burr等[10-13]从界面损伤的角度出发，用细观力学方法解释材料的宏观强度；但文章中只给出了材料的几种破坏假设，没有实际模拟材料的非线性损伤过程.刘昌明等人[14]在细观力学角度上，运用多尺度广义自洽模型模拟了镁碳质耐火材料的65损伤过程.1.2基于界面相模型的耐火材料细观损伤本构模型基于以上讨论，本课题组从细观力学的角度出发，提出了基于界面相模型的耐火细观损伤本构模型：在细观力学方法的基础上引入界面相模型的方法,根据界面相模型建立颗粒、基质、界面的弹性常数与耐火材料的弹性常数的关系,结合损伤力学理论,研究耐火材料的70损伤演化规律,建立宏观损伤变量与界面弹性常数的关系.在实现该本构关系的问题上，将耐火材料的行为过程看作分段的非线性弹性模型，基于材料行为过程中杨氏模量的变化，用杨氏模量的变化来表征耐火材料的损伤，结合界面相模型方法计算出的杨氏模量和界面参数的变化曲线，从而建立起在宏观和微观两个尺度上的本构关系.界面相模型即二次组装CSA模型，其等效力学性能由广义自洽模型[7]推广得到.基于界75面相模型的耐火材料损伤演化表征可以分为三个阶段：损伤起始阶段、基质损伤阶段和界面损伤阶段.损伤起始阶段是线弹性模型，杨氏模量E不变；后两个阶段里杨氏模量E变化，为分段的非线性弹性模型.为了表征损伤，以杨氏模量E的变化定义界面损伤阶段的损伤d,为01EdE=?（1）式中E0为初始时的杨氏模量值.则材料的损伤应力-纵向应变关系为80()01Edσξ=?（2）对镁碳质压缩试验测点测得的应力σ与纵向应变ξ的变化关系进行拟合，再通过解析法得到损伤d和纵向应变ξ的变化关系、杨氏模量E与纵向应变ξ的变化关系，然后结合运用界面相模型方法得到了杨氏模量E与界面参数q的变化关系，最终得到损伤d与界面参数q的变化关系.该本构模型技术路线图如图1所示.85图1理论模型技术路线图Fig1.Thetechnologyroadmapofthetheoreticalmodel每个阶段的变量间变化关系如下：90损伤起始阶段：d0,1,0.00004~0qξ===?(3)基质损伤阶段：1036231541238267.02101.33101.23100.0011.645101.04101.263101.02610400.00078~0.00004dqqqqξξξξξ????=××?××+××??=?××?××+××???××???=???(4)界面损伤阶段：1037241038261.554105.411108.288100.0866.123104.872101.3511069.480.00031~0.00078dqqqqξξξξ????=××?××+××?+??=?××?××?××????=????(5)95公式（1）、（2）、（3）中，d表示损伤，q表示界面参数，ξ表示纵向应变.2ABAQUS用户子程序开发2.1ABAQUS用户子程序功能ABAQUS10.1一共有42个用户子程序接口,15个应用程序接口,用户可以定义包括边界条件、荷载条件、接触条件、材料特性以及利用用户子程序和其它应用软件进行数值交换等.100这些用户子程序接口使得用户在解决一些问题时有很大的灵活性,同时这些接口也大大扩充了ABAQUS的功能.常用的接口程序有：用于定义材料力学行为的UMAT、用于定义场变量的USDFLD、用于定义不同单元模式的UEL和自定义硬化变量子程序UHARD等，涉及材料、单元、算法等多个方面.本文根据提出的细观力学损伤模型，介绍利用用户子程序UMAT和USDFLD来实现该本构模型，在相应的ABAQUS输入*.inp文件中，输入“*USER105MATERIAL”关键字表示采用用户自定义材料属性，计算时会自动调用用户材料子程序.UMAT子程序可以定义ABAQUS材料库中没有包含的各类材料本构模型，以扩充材料库；UMAT可以用于力学行为分析的任何分析过程，几乎可以把用户材料属性赋予ABAQUS中的任何单元；另外还可以与USDFLD一起使用，通过“USDFLD”来重新定义单元每一个积分点上传递到UMAT中的场变量的数值.1102.2本构模型的子程序实现UMAT子程序的核心内容是定义材料本构模型的雅克比矩阵，即应力增量对应变增量的变化率?σ/?ξ.为了验证本构模型的思路有效性，将它看作是分段的非线性各向同性弹性模型，所以UMAT定义的?σ/?ξ=D，弹性刚度矩阵D的定义如下：()()()()()()()()()()()()()()()()()()()()()()()()100011211211210001121121121000112112112021002100000210000021EEEEEEEEEDEEEνννννννννννννννννννννννννννννν?????+?+?+????????+?+?+????????+?+?+?=????+??????+??????+??115其中E、μ分别代表弹性模量和泊松比，整个非线性力学行为过程中，弹性模量E是不断变化的，E=E0(1-d);泊松比是不变的.每个阶段对应的损伤d、界面参数q和纵向应变ξ间的变化关系如上式（1）、（2）、（3）.对应于一个应变增量，可以求得相应的界面参数q，进一步可以得到损伤d，根据损伤d和杨氏模量的关系得到杨氏模量E，从而可以求得数值求解过程中的弹性刚度矩阵D,最后更新应力、状态变量和其他数据.120基于上面所述的材料本构模型、应力更新算法和场变量之间的传递关系，将纵向应变ξ和界面参数q定义为场变量1和场变量2，纵向应变ξ、界面参数q、损伤变量d和杨氏模量E为状态变量，在每一个增量步进行