简介面向Abaqus Explicit分析用户的复合材料三维连续介质损伤力学CDMVUMAT子程序实现资源适合从事复合材料失效仿真、需要编写或移植自定义材料本构的研发工程师。文件包共13个文件包含1个Fortran子程序源文件、6个覆盖拉伸、压缩、剪切等基本工况的inp算例、2个Python脚本分别用于原位性能换算与子程序自动验证以及使用说明、状态变量字典、License等辅助文档压缩包整体仅33KB。模型在inp中完整定义纤维/横向弹性模量、泊松比、剪切模量、强度、断裂韧性及破坏面角等属性便于替换参数后直接开展渐进损伤分析运行前需将Abaqus与Fortran编译器及Visual Studio正确链接。已有1600余人学习是快速上手VUMAT二次开发和复合材料损伤模拟的实用参考。1. VUMAT与三维连续介质损伤复合材料仿真从弹性到失效的那一步Abaqus Explicit 自带了不少复合材料损伤模型比如 Hashin 准则配合单元删除但实际做工程仿真的人都知道内置模型解决不了三个问题一是三维应力状态下层间与层内耦合失效的描述二是刚度渐进退化过程中对网格依赖的控制三是把试验测到的断裂能、损伤演化规律以自定义方程写进求解器。这时就得动用自己的材料子程序也就是 VUMAT。用 Fortran 写 VUMAT 并不是把公式搬进代码那么简单它同时要求你理解 Abaqus 的增量式求解框架、显式算法对时间增量的限制以及复合材料每一层在三维应力空间里损伤如何起始、如何演化、如何影响刚度和强度。这篇文章既不打算第一句就丢出一大段损伤张量公式也不打算通篇停留在概念层面而是按我在工程里实际开发 VUMAT 的路径走一遍从连续介质损伤力学的理论框架到可运行的 Fortran 骨架再到参数调试和排错最后给一个能拿去做单胞验证和批量计算的思路。适合已经在用 Abaqus、想从搭搭内置模型转向二次开发的工程师也适合刚接触 VUMAT 但迫切需要一个完整落点的人。2. 三维连续介质损伤力学模型的理论基础从应力更新到损伤演化2.1 连续介质损伤力学的基本框架损伤变量、有效应力、应变等效连续介质损伤力学CDM的核心思路是把材料内部的微裂纹、微孔洞用连续的损伤变量来表征而不是显式地建模每一条裂纹。对于复合材料我们通常不定义一个各向同性损伤系数而是按材料主方向定义多个损伤分量纤维拉伸、纤维压缩、基体拉伸、基体压缩甚至分层。三维模型下损伤状态由一组标量d_i表示刚度矩阵在每一次增量步里都乘以(1 - d_i)对应的退化因子。有效应力的概念是 CDM 的基石。经典形式是sigma_eff sigma / (1 - d)在复合材料三维正交各向异性材料中应力应变关系写成sigma C(d) : epsilon这里的C(d)不再是常数而是损伤变量d的函数。实现 VUMAT 时最常见的做法是引入一组内部变量SDV在每个增量步开始时读入d_i更新刚度矩阵再计算应力。Abaqus Explicit 的 VUMAT 接口要求你给出从应变增量到应力增量的显式关系所以整个模型会被嵌入一个更新流程中。提示显式求解器没有全局迭代VUMAT 必须自己处理材料非线性的数值稳定性。损伤变量突跳容易引起应力波振荡。2.2 复合材料三维本构关系平面应力只是特例三维才是常态很多内置失效准则是在壳单元或平面应力假设下给出的但真实层合板的自由边、开孔、冲击区域都是三维应力状态。三维正交各向异性本构需要 9 个独立弹性常数E1、E2、E3、nu12、nu13、nu23、G12、G13、G23。写成刚度矩阵C是 6x6 对称矩阵VUMAT 中需要根据这些工程常数构建。损伤后的刚度矩阵不能简单把某个模量减去一个百分比因为剪切项和泊松耦合项必须满足热力学一致性。常见做法是对应变能密度函数进行分解纤维方向分量、基体方向分量、剪切分量分别乘以(1 - d_i)的退化系数。一个经过验证的展开形式是C(1,1) (1.0 - df) * ( E1 * (1.0 - nu23*nu32) / delta ) C(2,2) (1.0 - dm) * ( E2 * (1.0 - nu13*nu31) / delta ) C(1,2) (1.0 - df) * (1.0 - dm) * ( E1 * (nu21 nu31*nu23) / delta )其中delta 1.0 - nu12*nu21 - nu23*nu32 - nu31*nu13 - 2.0*nu21*nu32*nu13。这个写法保留了对称性同时让每个损伤变量只作用于对应工程常数所在的刚度项。注意不要把df和dm直接乘到E1、E2上后就拿来算剪切模量。交叉项的退化会出现非物理负刚度尤其在压缩和剪切联合作用时。对于三维模型在 Abaqus 中采用 C3D8R 或 C3D20R 实体单元每个积分点上会有 6 个应变分量。VUMAT 收到的strainInc是正增量你需要把它加到总应变上再更新应力。不同单元类型对应不同的截面点组织方式这就是后文要说的多截面点问题。2.3 损伤起始准则与演化方程从 Hashin 到断裂能正则化损伤起始用三维 Hashin 准则是一步自然选择因为它区分了纤维和基体的拉伸与压缩。三维形式下纤维拉伸条件写为(df_start (sigma_11 / XT)^2 (tau_12 / S12)^2 (tau_13 / S13)^2 1)基体拉伸条件则要考虑 22、33 和剪切项。这只是起始判断真正决定模型行为的是损伤演化。单纯用应力跌落会造成网格依赖所以主流 VUMAT 实现都在演化方程里引入断裂能密度Gc和特征长度L。在 Abaqus Explicit 中特征长度由单元几何计算你可以在 VUMAT 里通过characteristicLength参数部分版本通过charLength获得。演化方程可以写成指数软化形式d_i 1.0 - exp(-(1.0 / m) * ( (strain_eq - strain_0) / strain_f )**m )这里strain_eq是等效损伤应变strain_0是损伤起始阈值strain_f由断裂能控制。在显式分析中如果时间增量过大损伤变量会在一个增量步内跳到接近 1导致应力跌落到零并引发单元畸变。所以你的 VUMAT 里最好加一个限制每个增量步的损伤增量不超过 0.05或者用上一增量步的应力状态做线性插值。3. VUMAT 接口与 Fortran 代码结构把 3D 损伤模型写进 Abaqus Explicit3.1 VUMAT 的调用约定理解 nblock、nstatev 与截面点组织VUMAT 和 UMAT 最大的不同在于它收到的是一个块block里多个积分点的数据而不是单个点。Abaqus 对 VUMAT 的签名是subroutine vumat( c 只读变量 1 nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal, 2 stepTime, totalTime, dt, cmname, coordMp, charLength, 3 props, density, strainInc, relSpinInc, 4 tempOld, stretchOld, defgradOld, fieldOld, 5 stressOld, stateOld, enerInternOld, enerInelasOld, 6 tempNew, stretchNew, defgradNew, fieldNew, 7 stressNew, stateNew, enerInternNew, enerInelasNew )其中nblock是此调用中需要处理的积分点总数它通常远大于单元数。nstatev是你自己定义的内部变量个数。charLength是每个点的特征长度。strainInc是当前增量步的应变增量数组维度为nblock * ndir在三维问题中ndir3nshr3所以每个积分点的应变分量顺序是11,22,33,12,23,13。对于普通实体单元Abaqus 已经帮你把应力应变旋转到了材料方向。但复合材料层合板如果用 solid section 并指定了 orientation那么strainInc是在材料坐标系下的分量。你要确保 Fortran 里数组顺序与 Abaqus 文档一致。3.2 最小可行的 Fortran 骨架从应变增量到应力更新下面这个骨架实现了 3D 正交各向异性弹性响应加上简单的损伤退化不涉及断裂能正则化但结构完整。我在工程开发中会先跑通这个弹性版本再逐步加入损伤起始和演化。subroutine vumat( 1 nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal, 2 stepTime, totalTime, dt, cmname, coordMp, charLength, 3 props, density, strainInc, relSpinInc, 4 tempOld, stretchOld, defgradOld, fieldOld, 5 stressOld, stateOld, enerInternOld, enerInelasOld, 6 tempNew, stretchNew, defgradNew, fieldNew, 7 stressNew, stateNew, enerInternNew, enerInelasNew ) c implicit none integer nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal real*8 stepTime, totalTime, dt character*80 cmname real*8 coordMp(nblock,*), charLength(nblock) real*8 props(nprops) real*8 density(nblock) real*8 strainInc(nblock,ndirnshr) real*8 relSpinInc(nblock,nshr) real*8 tempOld(nblock), stretchOld(nblock,*) real*8 defgradOld(nblock,*), fieldOld(nblock,*) real*8 stressOld(nblock,ndirnshr), stateOld(nblock,nstatev) real*8 enerInternOld(nblock), enerInelasOld(nblock) real*8 tempNew(nblock), stretchNew(nblock,*) real*8 defgradNew(nblock,*), fieldNew(nblock,*) real*8 stressNew(nblock,ndirnshr), stateNew(nblock,nstatev) real*8 enerInternNew(nblock), enerInelasNew(nblock) c integer i, j, k real*8 e1, e2, e3, nu12, nu13, nu23, g12, g13, g23 real*8 nu21, nu31, nu32, delta real*8 C(6,6), d(6), strain(6), stress(6) real*8 df, dm c c 假设 props 按顺序传递 9 个弹性常数和 2 个初始损伤值 c props(1)E1, props(2)E2, props(3)E3 c props(4)nu12, props(5)nu13, props(6)nu23 c props(7)G12, props(8)G13, props(9)G23 e1 props(1) e2 props(2) e3 props(3) nu12 props(4) nu13 props(5) nu23 props(6) g12 props(7) g13 props(8) g23 props(9) c nu21 nu12 * e2 / e1 nu31 nu13 * e3 / e1 nu32 nu23 * e3 / e2 delta 1.0 - nu12*nu21 - nu23*nu32 - nu13*nu31 1 - 2.0*nu21*nu32*nu13 c do i 1, nblock c 读取损伤变量这里用 stateOld(1) 和 (2) df stateOld(i,1) dm stateOld(i,2) c 构建退化后的刚度矩阵 C(1,1) (1.0 - df) * e1 * (1.0 - nu23*nu32) / delta C(2,2) (1.0 - dm) * e2 * (1.0 - nu13*nu31) / delta C(3,3) (1.0 - dm) * e3 * (1.0 - nu12*nu21) / delta C(4,4) (1.0 - df) * g12 C(5,5) (1.0 - dm) * g23 C(6,6) (1.0 - df) * g13 C(1,2) (1.0 - df) * (1.0 - dm) 1 * e1 * (nu21 nu31*nu23) / delta C(2,1) C(1,2) C(1,3) (1.0 - df) * (1.0 - dm) 1 * e1 * (nu31 nu21*nu32) / delta C(3,1) C(1,3) C(2,3) (1.0 - dm) * e2 * (nu32 nu12*nu31) / delta C(3,2) C(2,3) c do j 1, ndir nshr strain(j) stateOld(i,j3) strainInc(i,j) end do c call matvec33(C, strain, stress, ndir, nshr) c 更新应力这里先把剪切分量单独处理 do j 1, ndir stressNew(i,j) stress(j) end do stressNew(i,4) stress(4) stressNew(i,5) stress(5) stressNew(i,6) stress(6) c c 更新内部变量应变存到 state 4~9损伤变量保持原值 do j 1, ndir nshr stateNew(i,j3) strain(j) end do stateNew(i,1) df stateNew(i,2) dm end do c return end代码逻辑说明先根据工程常数计算柔度相关项构建损伤退化后的刚度矩阵。注意裁剪项C(4,4)只乘了(1-df)表示纤维剪切损伤C(5,5)乘(1-dm)对应基体剪切。实际使用时你需要根据失效模式调整退化因子组合避免剪切和拉伸产生相同退化。matvec33是矩阵向量乘自己写一个不用内建matmul是因为 Fortran 指针对数组维度的处理在 Abaqus 编译器下容易出错。参数stateOld(i,j3)用来保存总应变方便下一步做等效应变。3.3 状态变量与多截面点数据的存储策略Abaqus 可视化结果中每个积分点只能看到你通过配套VUMAT状态变量输出的值。建议一个固定排列1~2损伤变量3损伤起始函数最大值4~9应变10能量密度。这样在 Abaqus 后处理里可以直接画损伤云图。当模型里有大量实体单元尤其是一层铺层只画一个单元厚度时每个单元有多个积分点nblock会包含来自不同单元和不同截面点的混合数据。不要在 VUMAT 里通过coordMp判断单元编号来区分材料而应该依赖cmname。如果你需要对同一模型赋予不同材料方向或不同刚度就写多个VUMAT或者用props里的一个整数作为材料标识。nstatev在.inp文件的*Depvar里定义多截面点不会影响状态变量数量。4. 在 Abaqus Explicit 中调试与验证参数、稳定性与常见坑4.1 复合材料 3D 模型的参数表从单层刚度到失效强度VUMAT 的材料参数通常以props数组传入。我一般会准备一份参数表放到代码注释或一个独立文本文件里参数编号含义典型值说明1E1135 GPa纤维方向模量2E210 GPa基体横向模量3E310 GPa厚度方向通常等于 E24nu120.3主泊松比5nu130.3厚度方向泊松比6nu230.4横观各向同性面内泊松比7G125 GPa面内剪切模量8G135 GPa剪切模量9G233.5 GPa横向剪切模量10XT2200 MPa纤维拉伸强度11XC1200 MPa纤维压缩强度12YT60 MPa基体拉伸强度13YC200 MPa基体压缩强度14S1280 MPa面内剪切强度15Gf80 N/mm纤维断裂能16Gm1 N/mm基体断裂能注意props数组在 VUMAT 里是按浮点数读入的如果参数表里混有整数必须写成1.0。在.inp文件中*UserMaterial, constants16后按顺序写数值。实际调试时先只用前 9 个参数跑弹性响应确认应力应变线性关系正确再逐步加入强度参数。4.2 时间增量与稳定性为什么损伤后计算总中断显式分析的时间增量由最小单元特征长度和材料波速决定。VUMAT 引入刚度退化后材料波速变化临界时间增量也随之改变。Abaqus 默认的固定时间增量会导致单元失真甚至负体积。常见的解决方案是使用*Fixed Time Incretion的替代项*Variable Mass Scaling或者直接在跨 increat 里给 VUMAT 加最大损伤增量限制。我的经验是把损伤增量的限制与单元特征长度挂钩比如要求d_i在单个增量步中的变化不超过0.2 * charLength / (C11 * dt)的某个倍数这样能有效抑制数值振荡。需要特别强调的是Abaqus Explicit 的中断不总是报错。如果你看到The time increment is less than the minimum specified这类提示通常意味着单元刚度剧烈变化而不是单纯的时间步长设置问题。这时回到模型检查材料方向是否与单元坐标系一致尤其是曲面上的实体单元。4.3 常见运行错误libpng error、中断不了怎么办很多 VUMAT 开发者在提交计算后遇到libpng error这不是子程序的问题而是 Abaqus/CAE 在生成缩略图时与高分辨率显示器或显卡驱动冲突。解决方法是关闭后处理的缩略图预览或者用命令行abaqus jobxxx cpus4 interactive提交避免打开 CAE 界面。这与子程序本身的编译无关但会误导你半天时间。真正和 VUMAT 相关的坑是“计算中断不了怎么办”。显式分析中如果时间增量被无限缩小作业会卡在某个增量步不动。原因多半是损伤变量的状态在相邻两步间跳变导致刚度矩阵不可逆。你可以在 VUMAT 中加入状态变量检查当某个损伤变量超过 0.99 时直接把应力降至一个极小值而不是继续用退化后的刚度计算。另一个排查方法是输出每个增量步的dt和最小特征长度在.msg文件里观察是否随时间持续缩小。5. 进阶单胞验证与 GPU 加速下的批量参数反演5.1 用单元素测试验证 VUMAT 的每个损伤模式不要一上来就跑整个层合板。对每一个损伤模式写一个只含一个单元的最小模型单轴拉伸、单轴压缩、纯剪切。在.inp里施加一个光滑的位移载荷防止高频噪声。后处理时画应力-应变曲线确认损伤起始点与理论强度一致并检查软化段的斜率是否由断裂能主导。为了提高验证效率写一个 Python 脚本批量生成inp并提交多个任务。我通常会把 VUMAT 的props设计成可以通过*UserMaterial的constants读取这样不需要每次重编译 Fortran。例如把不同的XT和Gf组合放进去跑一组参数扫描。这个做法在工程上用于拟合试验数据也用于评估模型对参数的敏感性。5.2 多截面点状态映射从壳单元到实体单元的跨层对比复合材料层合板如果采用“每层一个实体单元”的建模方式截面点数量巨大。Abaqus 后处理中默认只会显示单元中心的 SOV 值如果你在 VUMAT 中输出了每个积分点的损伤变量可以通过*Element Output指定SDV在积分点上输出。进入 CAE 后用Result - Options - Integration Point就能看到不同层截面的损伤差异。当模型由壳单元和实体单元混用时状态变量的映射会导致SDV含义不一致。建议将 VUMAT 的状态变量输出固定为与单元坐标系无关的值比如应变量、损伤量而不是应力分量。因为实体单元的 material orientation 与壳单元的 section orientation 在坐标系旋转上存在细微差异直接对比应力会出现镜像效应。5.3 用 GPU 加速批量计算参数反演的应用思路Abaqus Explicit 支持 GPU 计算但 VUMAT 的 Fortran 代码运行在 CPU 端GPU 加速收益主要来自显式求解器中的接触和体单元应力计算。如果你的模型内 VUMAT 占用了大量 CPU 时间建议先用 Intel Fortran 编译器开启-O3和-xHost优化再考虑并行。实际对比中我发现单模型 VUMAT 的耗时瓶颈往往在状态变量的读写而不是浮点运算量所以不要在状态变量里存储不必要的大数组。批量参数反演时可以用一个简单的 Python 包装脚本调用abaqus job利用多核 CPU 并行跑多个小模型而不必把所有模型拼成一个大作业。这样每个作业的 VUMAT 都是独立加载不会因为某个作业中断影响其他作业。最后用 pandas 汇总所有.odb中的应力应变数据画成散点图直接观察断裂能对峰后软化速率的影响。本文还有配套的精品资源点击获取