★ FEATURED ARTICLE Abaqus 二次开发教程14VUMAT 与显式用户材料——按块推送、没有材料雅可比但本构刚度直接决定稳定时间增量版本声明块工具/软件Abaqus 20262026 FD01Abaqus/Explicit显式求解器首次出现给全称自 2026 FD01 起在 Linux 系统的 NVIDIA GPU 上支持 GPU 加速见官方 Release Notes语言/环境Fortran 用户子程序固定格式Intel FortranoneAPI 2025.0 起 ifort 移除须用 ifx脚本侧用 Python 3.10.5本文目标读懂 VUMAT 的接口契约谁在什么时候把哪些数组交给它、它必须回写什么写出一个可编译、可提交、可与内置材料对照校验的 VUMAT一句话结论VUMAT 没有收敛问题它不参与 Newton 迭代、不回写DDSDDE接口里根本没有DDSDD这个通道但它按块接收NBLOCK个材料计算点并必须回写STRESSNEW/STATENEWAbaqus/Explicit 默认从本构响应数值估计有效剪切模量DTARRAY(2:NBLOCK1)与有效体积模量DTARRAY(NBLOCK2:2*NBLOCK1)来反推稳定时间增量所以本构写硬会把稳定增量压小只有指定EFFMOD时用户写入的模量才真正参与时间增量计算。〇、本篇要解决的认知问题VUMAT 与 UMAT 的接口契约到底差在哪调用粒度、必须回写的数组、是否必须提供材料雅可比矩阵material Jacobian记为DDSDD、剪切分量按张量分量还是工程分量存储为什么 VUMAT 没有收敛问题却仍有精度与稳定性问题显式的稳定时间增量与 VUMAT 里那份本构刚度是什么关系RHS、BIGH、BIGD、DINC、KMAX、K这些名字现在还在 VUMAT 参数里吗它们与当前接口的STRAININC、RELSPININC、JINFOARRAY怎么对应状态变量solution-dependent state variables在显式里怎么传递删除标志为什么必须用状态变量而不是返回值向量化KMAX/K/KSIZE那类块循环与性能的关系以及什么样的写法会把向量化打回逐点循环VUMAT 能挂哪些单元类型哪些单元类型不传入拉伸张量与变形梯度显式 GPU 加速2026 FD01 起对 Fortran 用户子程序意味着什么一、机制解析1.1 契约对照表先看清谁负责什么UMATAbaqus/Standard隐式与 VUMATAbaqus/Explicit显式看似都是用户材料子程序但责任划分完全不同隐式把迭代的责任交给用户显式把时间步的责任交给用户。契约项UMATAbaqus/StandardVUMATAbaqus/Explicit调用粒度逐个积分点调用按块调用一次传入NBLOCK个材料计算点一个块内可含多个单元、层与截面点材料雅可比必须提供DDSDDE【铁律】第 5 条留空或置零会静默给出错误刚度接口中没有DDSDD也没有任何雅可比数组必须回写DDSDDE、STRESS、STATEVSSE/SPD/SCD/RPL按需STRESSNEW、STATENEW可选ENERINTERNNEW、ENERINELASNEWDTARRAY段按需增量语义必须把应力与状态变量更新到增量末值【铁律】第 6 条同样必须是增量末值差别在于调用者不知道末在哪——它按当前显式增量被调用剪切分量NSHR个工程剪切分量NSHR个张量剪切分量不是工程值——同一份 J2 公式搬过来会错一个因子 2分量顺序按单元库约定顺序对称张量NDIRNSHR3D 为 11,22,33,12,23,31非对称张量NDIR2*NSHR变形梯度与 UMAT 的DFGRD0(3,3)形状不同客观应力率由材料/单元选项决定Green-Naghdi格林-纳格迪客观应力率内置连续体单元多用 Jaumann有限转动伴随有限剪切时二者结果会明显不同参考基局部取向若指定随材料旋转全部应力、应变、拉伸与状态量都在随材料旋转的局部材料基corotational数据检查阶段KINC0的首次调用用于求初始雅可比data check 阶段用虚构应变调用一次STEPTIMETOTALTIME0.0用于校验本构并算初始等效材料常数与初始弹性波速只读数组可否改写STRAN、PROPS、DFGRD0/1等只读PROPS、DENSITY、CHARLENGTH、STRAININC、RELSPININC与各*OLD只读官方明确改写passed in for information的变量后果不可预测状态变量STATEV就地更新跨增量靠求解器保存STATEOLD读入、STATENEW回写逐增量传递删除标志必须占用一个状态变量退出方式工具函数XITXPLB_EXIT不得用STOP精度由安装参数决定由VABA_PARAM.INC决定默认单精度Fortran includeINCLUDE ABA_PARAM.INCINCLUDE VABA_PARAM.INC1.2 逐增量逐积分点推送数据是怎么到达 VUMAT 的Abaqus/Explicit 主循环每个显式增量 │ │ 把同一材料 同一单元类型 同一层/截面点位置的材料计算点 │ 打包成一个 block长度 NBLOCK ▼ ┌──────────────────────────────────────────────────────────┐ │ VUMAT(NBLOCK, NDIR, NSHR, NSTATEV, NFIELDV, NPROPS, │ │ JINFOARRAY, STEPTIME, TOTALTIME, DTARRAY, CMNAME) │ │ ┌────────────────────────────────────────────────────┐ │ │ │ km 1 到 NBLOCK │ │ │ │ 读STRESSOLD(km,*) STRAININC(km,*) STATEOLD │ │ │ │ 算本构刚度只作为物理量存在不作为矩阵 │ │ │ │ 写STRESSNEW(km,*) STATENEW(km,*) │ │ │ │ EFFMOD 时DTARRAY(1km) DTARRAY(NBLOCK1km)│ │ │ └────────────────────────────────────────────────────┘ │ └──────────────────────────────────────────────────────────┘ ▼ 主程序用有效模量算膨胀波速 c_d sqrt((K 4G/3)/ρ) → 稳定时间增量 Δt_stable 仅当 USER MATERIAL 指定 EFFMOD用户写入的模量才参与时间增量计算这张图解释了 VUMAT 与 UMAT 最根本的差别UMAT 里的DDSDDE是求解器用来组装方程的输入VUMAT 里没有这个通道求解器只认两件事——用户回写的应力以及可选用户回写的有效模量。所以 VUMAT 里写错本构刚度不会收敛失败它只会以另一种方式暴露要么时间增量被压小导致算不完要么应力更新在物理上不一致导致能量异常。1.3 没有收敛问题但有精度与稳定性问题显式是条件稳定conditionally stable算法单步误差不会通过迭代被修正误差只能靠步长足够小来控制。因此 VUMAT 的三个稳定性出口是稳定增量与材料刚度的关系Δt_stable由材料膨胀波速与单元特征长度CHARLENGTH决定【经验法则】本构过硬会压低稳定增量——同一套网格下若 VUMAT 的等效模量比内置材料高m倍稳定增量按1/sqrt(m)缩小。官方警告本构高度非线性时数值估计的有效模量可能不够保守不保守意味着时间增量偏大此时应在EFFMOD下直接写入有效剪切与体积模量。【经验法则】工程上更常见的做法是刻意写入偏软的模量换取速度但物理上不能低于真刚度的三分之一否则误差靠步长控制不住这是用速度换精度的交易不是免费午餐。能量一致性显式的能量平衡靠ENERINTERNNEW内能与ENERINELASNEW耗散能不写这两个值时能量平衡无法给出有意义结果。塑性耗散在显式里是硬算出来的不能指望迭代收敛后能量自然平衡。单精度累积误差显式默认单精度长程历史变量等效塑性应变、损伤在几十万个增量上累积舍入误差。【最佳实践】把历史变量存进STATENEW而不是每次重算并在单元素对比测试里把与内置材料的容差控制在 1% 以内第 2.3 节。过渡句稳定性由本构刚度这一物理量决定性能由块这一数据形状决定。下面两节就是 VUMAT 工程化的两个日常代价。1.4 向量化一个条件分支能把整段代码打回逐点循环VUMAT 是矢量接口vector interface显式按块传入数据一次收到NBLOCK个材料点的应力、应变、状态变量用户必须写DO KM 1, NBLOCK这样的块循环这就是历史命名KMAX/K的角色——块大小与块内索引。向量化收益来自整块执行同一串算术运算代价是条件分支破坏向量化一旦IF只对块内部分材料点成立例如若 Mises 应力超过屈服面则做径向返回编译器无法把这段代码向量化退化为逐点循环。【最佳实践】把分支写在整块层面全块同一状态或按弹性/屈服分桶单点分支不可避免时接受性能损失而不是追求代码更清晰。第 2.1 节代码里这处分支是有意保留的教学示例。STOP、文件 IO、调用外部子程序都会破坏向量化显式里禁止用STOP用XPLB_EXIT调试输出只走 Fortran 单元 6写入显式.log或 15–18 与大于 100 的单元【经验法则】其余单元号归 Abaqus 自用写文件必须用完整路径名且显式并行下每个处理器都会生成一个文件。用户自写的子程序与COMMON块名必须以字母 K 开头官方约定。向量化与并行不是一回事显式的并行按域分解切分单元集块是子程序内的数据视图多线程/MPI 环境下COMMON块变量不能跨线程共享除非按官方并行约定做保护。1.5KSIZE、JINFOARRAY与内存布局陷阱状态变量数组的基址寻址是显式子程序最容易出事的一层契约。历史接口里用户按BEGSUB到BEGSUB KSIZE - 1的范围定维与访问BEGSUB是本次调用数组在公共存储区里的起始位置KSIZE是可读写的元素个数。当前文档给出的写法是JINFOARRAY(*)加POINTER/LOC定位JINFOARRAY(I_INFO_ELEMNUMSTARTLOC)给出块内单元编号的起始位置本质是同一类基址 长度的寻址模型PARAMETER常量的具体取值以官方文档为准官方接口示例里I_INFO_EFFMODDEFN 5。这条契约最容易出问题的三层越界即污染求解器内存JINFOARRAY(*)是假设大小数组assumed-size不是用户私有内存越界写入会破坏 Abaqus 自身的数据现象是结果不对但没有任何报错。【铁律】相关只能改写官方标注为variables to be defined / can be updated的数组【铁律】第 10 条脚本与代码必须可复制即运行含完整声明。大数组必须从堆分配显式里 VUMAT 每个增量都可能被调用栈上定维的大数组会 stack overflow 直接退出必须用可分配数组allocatable从堆分配。精度由 include 文件决定主程序与子程序之间传变量时VABA_PARAM.INC必须在所有共享变量的子程序里出现否则精度不一致导致类型/数量不匹配的平台相关链接或运行错误【铁律】第 7 条编译前必须先加载setvars.batcompile_fortran/link_sl/link_exe三处一致。1.6 单元类型限制与某些分量必须由用户定义单元类型分量规模用户必须额外定义什么shell / plane stressNDIR3, NSHR14 分量厚度方向应变增量STRAININC(*,3)靠平面应力或不可压缩约束求出若本构用到拉伸张量U还要给STRETCHNEW(*,3)、DEFGRDNEW(*,3)shell 横向剪切刚度无法自动得到必须在材料里定义弹性横向剪切模量或在截面里定义横向剪切刚度beam拉伸张量与变形梯度不可用空间梁需给STRAININC(*,2)、STRAININC(*,3)STRAININC(*,4)是扭转剪切厚度应力STRESSNEW(*,2)、(*,3)被假定为零写入值被忽略pipe拉伸张量与变形梯度不可用环向应力由压力载荷PE/PI/HPE/HPI/PENU/PINU预定义、不得修改用户必须给轴向应力STRESSNEW(*,1)与扭转剪切应力STRESSNEW(*,4)continuum实体单元NDIR3, NSHR36 分量无额外义务【最佳实践】写 VUMAT 前先确定它要挂哪些单元类型。UMAT 里所有单元都能跑的写法在显式里不成立同一个STRESSNEW(KM,3)在实体单元里是第三个直接应力在 shell 里是厚度方向在 beam 里被忽略。1.7 显式 GPU 加速2026 FD01 起与 Fortran 用户子程序官方 Release Notes 明确2026 FD01 起 Abaqus/Explicit 在 Linux 系统的 NVIDIA GPU 上支持 GPU 加速适用于多种动力学与准静态结构分析可在单节点单卡或多卡上运行。这一条对 VUMAT 的实际影响——GPU 求解路径与 Fortran 用户子程序能否共存、以什么形式共存须以官方文档为准。工程上可确定的判断是GPU 加速的价值在内置材料与大模型的矩阵/向量运算上一旦把本构搬进 Fortran 子程序性能由第 1.4 节的向量化条件决定而不是由硬件决定。【经验法则】评估 GPU 对含 VUMAT 模型的实际收益先跑同一模型的两个版本内置材料 vs VUMAT再对比墙钟时间与稳定时间增量不要凭宣传判断。1.8 历史命名与当前接口的对应关系以下名字在早期显式文档中是 VUMAT 参数当前文档已改用右列命名逐条以官方文档为准历史名常见含义当前接口中的对应RHS应力更新方程的右端项当前接口不再传入该数组应力由用户自行更新到STRESSNEWBIGH与拉伸张量U相关的量STRETCHOLD(KM,*)/STRETCHNEW(KM,*)BIGD与速率/变形梯度相关的量DEFGRDOLD(KM,*)/DEFGRDNEW(KM,*)DINC应变增量STRAININC(KM,*)注意是张量剪切分量KMAX/K块大小与块内材料点索引NBLOCK与循环索引KMJINFOARRAY提供积分点/层/截面点信息【最佳实践】从旧代码迁移时先核对参数数量与类型官方明确参数的数量与类型必须与文档一致不匹配会造成平台相关的链接错误或运行错误——这类错误常常表现为能编译、跑起来结果不对。二、完整代码与逐行剖析2.1 完整 VUMATFortran 固定格式可直接编译保存为vumat_mises.f。本构各向同性线弹性 各向同性硬化 Mises 塑性径向返回radial return。状态变量 1 为等效塑性应变状态变量 2 为单元删除标志1 有效、0 删除。C C VUMAT : Mises isotropic hardening, radial return C Abaqus/Explicit only; fixed form; Fortran (ifx) C PROPS(1)E PROPS(2)nu PROPS(3)sigma_y0 PROPS(4)Etan C NSTATEV must be 2 (DEPVAR), STATEOLD(1)eq plastic strain, C STATEOLD(2)deletion flag (1 active, 0 deleted) C SUBROUTINE VUMAT(NBLOCK, NDIR, NSHR, NSTATEV, NFIELDV, NPROPS, 1 JINFOARRAY, STEPTIME, TOTALTIME, DTARRAY, CMNAME, COORDMP, 2 CHARLENGTH, PROPS, DENSITY, STRAININC, RELSPININC, 3 TEMPOld, STRETCHOLD, DEFGRDOLD, FIELDOLD, 4 STRESSOLD, STATEOLD, ENERINTERNOLD, ENERINELASOLD, 5 TEMPNEW, STRETCHNEW, DEFGRDNEW, FIELDNEW, 6 STRESSNEW, STATENEW, ENERINTERNNEW, ENERINELASNEW) INCLUDE vaba_param.inc C JINFOARRAY constants: values in the official interface example; C cross-check with the installed include file PARAMETER (I_INFO_EFFMODDEFN 5) INTEGER NTENS, KM, I REAL TR, SINV, J2, SM, YH, R, DPEQ, WPD, DEV, DEVP REAL E, NU, G, KULK, SIGY0, ETAN DIMENSION DEV(10) C array shapes must match the official VUMAT interface exactly DIMENSION DTARRAY(2*NBLOCK1), PROPS(NPROPS), DENSITY(NBLOCK), 1 COORDMP(NBLOCK,*), CHARLENGTH(NBLOCK), 2 STRAININC(NBLOCK,NDIRNSHR), RELSPININC(NBLOCK,NSHR), 3 TEMPOld(NBLOCK), STRETCHOLD(NBLOCK,NDIRNSHR), 4 DEFGRDOLD(NBLOCK,NDIRNSHRNSHR), FIELDOLD(NBLOCK,NFIELDV), 5 STRESSOLD(NBLOCK,NDIRNSHR), STATEOLD(NBLOCK,NSTATEV), 6 ENERINTERNOLD(NBLOCK), ENERINELASOLD(NBLOCK), TEMPNEW(NBLOCK), 7 STRETCHNEW(NBLOCK,NDIRNSHR), DEFGRDNEW(NBLOCK,NDIRNSHRNSHR), 8 FIELDNEW(NBLOCK,NFIELDV), STRESSNEW(NBLOCK,NDIRNSHR), 9 STATENEW(NBLOCK,NSTATEV), ENERINTERNNEW(NBLOCK), 1 ENERINELASNEW(NBLOCK), JINFOARRAY(*) CHARACTER*80 CMNAME E PROPS(1) NU PROPS(2) SIGY0 PROPS(3) ETAN PROPS(4) G E / (2.0 * (1.0 NU)) KULK E / (3.0 * (1.0 - 2.0 * NU)) NTENS NDIR NSHR C ---- data check call: STEPTIME TOTALTIME 0.0 with fictitious C ---- strains; publish the effective moduli, echo the input ---- IF (STEPTIME .EQ. 0.0 .AND. TOTALTIME .EQ. 0.0) THEN IF (JINFOARRAY(I_INFO_EFFMODDEFN) .EQ. 1) THEN DO KM 1, NBLOCK DTARRAY(1 KM) G DTARRAY(NBLOCK 1 KM) KULK END DO END IF DO KM 1, NBLOCK DO I 1, NTENS STRESSNEW(KM, I) STRESSOLD(KM, I) END DO DO I 1, NSTATEV STATENEW(KM, I) STATEOLD(KM, I) END DO ENERINTERNNEW(KM) ENERINTERNOLD(KM) ENERINELASNEW(KM) ENERINELASOLD(KM) END DO RETURN END IF C ---- main increment loop: this loop is the vectorized body ---- DO KM 1, NBLOCK C deleted point: solver passes zero stress and zero strain IF (NSTATEV .GE. 2 .AND. STATEOLD(KM, 2) .EQ. 0.0) THEN DO I 1, NTENS STRESSNEW(KM, I) 0.0 STRETCHNEW(KM, I) STRETCHOLD(KM, I) END DO IF (NSTATEV .GE. 2) STATENEW(KM, 2) 0.0 IF (NSTATEV .GE. 1) STATENEW(KM, 1) STATEOLD(KM, 1) ENERINTERNNEW(KM) ENERINTERNOLD(KM) ENERINELASNEW(KM) ENERINELASOLD(KM) CYCLE END IF C trial stress hydrostatic part 2*G*deviatoric strain increment TR 0.0 DO I 1, 3 TR TR STRAININC(KM, I) END DO SINV KULK * TR DO I 1, NTENS IF (I .LE. 3) THEN DEV(I) 2.0 * G * (STRAININC(KM, I) - TR / 3.0) ELSE DEV(I) 2.0 * G * STRAININC(KM, I) END IF STRESSNEW(KM, I) DEV(I) END DO DO I 1, 3 STRESSNEW(KM, I) STRESSNEW(KM, I) SINV END DO J2 0.5 * (DEV(1)**2 DEV(2)**2 DEV(3)**2) DO I 4, NTENS C shear components are tensor components: single weight, not 2 J2 J2 DEV(I)**2 END DO SM SQRT(3.0 * J2) C ---- branch below applies to single points: it breaks C ---- vectorization (see section 1.4) ---- YH SIGY0 ETAN * STATEOLD(KM, 1) IF (SM .LE. YH) THEN DO I 1, NSTATEV STATENEW(KM, I) STATEOLD(KM, I) END DO DPEQ 0.0 ELSE R YH / SM DO I 1, NTENS IF (I .LE. 3) THEN STRESSNEW(KM, I) SINV R * DEV(I) ELSE STRESSNEW(KM, I) R * DEV(I) END IF END DO DPEQ (SM - YH) / (3.0 * G) DO I 1, NSTATEV STATENEW(KM, I) STATEOLD(KM, I) END DO STATENEW(KM, 1) STATEOLD(KM, 1) DPEQ WPD SM * DPEQ C inelastic heat fraction is not defined here: temperature is C not integrated inside VUMAT, dissipated energy is reported ENERINELASNEW(KM) ENERINELASOLD(KM) WPD END IF C plastic part of the strain increment feeds the stretch tensor C (the stretch tensor is not available for beam and pipe elements) DO I 1, NTENS IF (SM .GT. YH) THEN DEVP DEV(I) / SM STRETCHNEW(KM, I) STRETCHOLD(KM, I) 1 ((1.0 - R) * DEVP) * (3.0 * DPEQ / 2.0) ELSE STRETCHNEW(KM, I) STRETCHOLD(KM, I) END IF END DO ENERINTERNNEW(KM) ENERINTERNOLD(KM) 1 0.5 * (SM * DPEQ KDOTPROD(DEV, NTENS)) IF (JINFOARRAY(I_INFO_EFFMODDEFN) .EQ. 1) THEN DTARRAY(1 KM) G DTARRAY(NBLOCK 1 KM) KULK END IF END DO RETURN END C helper: user-written routines must start with the letter K REAL FUNCTION KDOTPROD(A, N) DIMENSION A(10) INTEGER N, I KDOTPROD 0.0 DO I 1, N KDOTPROD KDOTPROD A(I) * A(I) END DO RETURN END逐段看决定成败的点INCLUDE VABA_PARAM.INC必须是参数表之后的第一条语句且所有共享变量的子程序都要带上它——显式的精度由这个文件决定缺它会在链接期出现类型/数量不匹配【铁律】第 7 条的编译环境必须先行加载。数组形状必须逐字对应官方接口DEFGRDOLD(NBLOCK, NDIRNSHRNSHR)是非对称张量的形状NDIR2*NSHR与 UMAT 的DFGRD0(3,3)不同STRAININC(NBLOCK, NDIRNSHR)里的剪切是张量分量。把 UMAT 的写法直接搬过来结果不会报错只会错。J2的权重是本文最容易写错的一行直接分量权重 0.5、剪切分量权重 1因为显式存的是张量剪切分量而 UMAT 存的是工程剪切分量。这就是第 1.1 节对照表里那一行的落地后果。IF (SM .LE. YH)是单点分支向量化代价写在注释里。想保住向量化要么把本构写成无条件更新的闭合形式要么按整块是否屈服分桶第 1.4 节。DTARRAY的下标第 1 个元素是时间增量本身DTARRAY(2:NBLOCK1)是各点剪切模量、DTARRAY(NBLOCK2:2*NBLOCK1)是各点体积模量——本代码按1 KM与NBLOCK 1 KM写入。只有EFFMOD被指定时用户写入的模量才参与时间增量计算否则显式自己从本构响应数值估计这正是本构过硬会压低稳定增量的实现路径。JINFOARRAY(I_INFO_EFFMODDEFN) 1表示必须更新模量I_INFO_*常量由 include 文件的PARAMETER给出官方接口示例里I_INFO_EFFMODDEFN 5本文按该示例硬编码实际取值以官方文档为准。删除标志只能靠状态变量表达STATENEW(KM,2)置 1 表示有效、置 0 表示删除块结构在分析中不变已删除的点不会被移出块求解器对已删除点传入零应力与零应变增量且一旦标记删除不能再激活。【最佳实践】把删除标志固定为最后一个状态变量并在输入文件用*DEPVAR预留。终止分析必须用XPLB_EXIT显式不能用STOPSTOP会让文件没有正常关闭.log/.odb可能被截断。2.2 编译环境加载与提交命令C: \Program Files (x86)\Intel\oneAPI\setvars.bat set ABA_PARALLEL_DEBUG1 C ---- abaqus_v6.env环境配置文件三处必须一致compile_fortran / link_sl / link_exe---- compile_fortran ifx -fpp -auto -assume modify link_sl ifx -lib dpc -assume modify link_exe ifx -assume modify -defaultstr:1024 fortran_type fixed提交一个显式作业abaqus jobvumat_check inputsingle_check.inp uservumat_mises.f cpus1ABA_PARALLEL_DEBUG1会把编译与链接命令打印出来排查编译通过但链接缺符号时先看这里【铁律】第 7 条。【最佳实践】把setvars.bat与abaqus_v6.env放在同一次会话里加载环境变量不会跨会话保留新开一个终端直接提交作业时编译器路径缺失是最常见的链接失败原因。2.3 与内置材料对照的 Python 侧可复制即运行单元素对照测试同一个 C3D8R 单元用内置弹性材料跑一遍再用*USER MATERIAL挂 VUMAT 跑一遍比较同一增量下的S11。importosimportsubprocess WORK_DIRos.getcwd()INP_NAMEsingle_check.inpINP*Heading ** single element check: built-in elasticity vs VUMAT *PART Name: block, ELEMENT: C3D8R 0., 0., 0. 1., 0., 0. 1., 1., 0. 0., 1., 0. 0., 0., 1. 1., 0., 1. 1., 1., 1. 0., 1., 1. *MATERIAL, NAMEREF_E *ELASTIC 21000.0, 0.3 *DENSITY 1.0e-9 *MATERIAL, NAMEUSER_M *USER MATERIAL, CONSTANTS4, EFFMOD 21000.0, 0.3, 200.0, 700.0 *DEPVAR 2 *ELEMENT, TYPEC3D8R, ELSETE_ALL 1, 1, 2, 3, 4, 5, 6, 7, 8 *SOLID SECTION, ELSETE_ALL, MATERIALREF_E *STEP, INC200000, NLGEOMON *STATIC *BOUNDARY 1, 1, 0.0 2, 2, 0.0 3, 3, 0.0 5, 1, 0.001 6, 1, 0.001 *OUTPUT, FIELD, FREQUENCY1 *ELEMENT OUTPUT S, EEQPL, SDV *END STEP withopen(INP_NAME,w,encodingutf-8)asf:f.write(INP)cmd[abaqus,jobvumat_check,inputINP_NAME,uservumat_mises.f,cpus1]print(提交命令, .join(cmd))subprocess.run(cmd,checkTrue)logvumat_check.logifos.path.exists(log):withopen(log,r,encodingutf-8,errorsreplace)asf:textf.read()hits[ln.strip()forlnintext.splitlines()iftime incrementinln.lower()]print(稳定时间增量相关行判定依据,hits[:3])print(用户材料是否被调用,USER MATERIALintextorvumatintext.lower())print(工作目录,WORK_DIR)判定要点客观vumat_check.log里必须出现稳定时间增量的数值把*USER MATERIAL换成内置*ELASTIC重跑同一增量下的S11相对误差应小于 1%弹性段进入塑性段后SDV(1)与内置材料的EEQPL在 2% 以内。【铁律】第 1 条本例整套单位是 mm–N–sE 21000即 N/mm²DENSITY 1.0e-9即 t/mm³E、sigma_y0、密度与长度必须同一套单位。三、常见报错与排查现象编译通过、链接失败报缺少符号或精度不匹配。根因没有先加载setvars.bat或compile_fortran/link_sl/link_exe三处不一致也可能是参数数量与类型与文档不符官方明确这会造成平台相关的链接或运行错误。解法先加载 oneAPI 环境再把三处统一到 ifx用set ABA_PARALLEL_DEBUG1看实际执行的命令逐条核对参数类型整型参数不能写成实参。现象跑起来结果不对但没有任何报错。根因改写了只读数组如STRAININC、PROPS、各*OLD或按BEGSUB/KSIZE基址越界访问JINFOARRAY——用户只能改写variables to be defined的数组越界会污染 Abaqus 自身内存。解法只在STRESSNEW/STATENEW以及可选的ENERINTERNNEW/ENERINELASNEW、DTARRAY段里写入需要额外存储时用可分配数组从堆分配不要用固定大数组占栈栈溢出会直接退出。现象显式作业的时间增量被自动压到很小墙钟时间是内置材料的数倍。根因VUMAT 的本构在数值上比内置材料硬例如把工程剪切分量当张量分量使用导致放大或弹性段用了明显偏大的有效模量。解法【经验法则】稳定增量按1/sqrt(模量比)变化先用单元素模型把本构的等效模量与内置材料对照1% 以内再决定是否显式指定EFFMOD。指定EFFMOD后用户写入的模量才真正参与时间增量计算。现象性能比预期差同一份本构在 UMAT 里没问题在 VUMAT 里慢。根因块内单点条件分支破坏了向量化或在 VUMAT 里做了文件 IO、STOP、外部库调用。解法把条件判断上提到整块用XPLB_EXIT代替STOP调试输出只走 Fortran 单元 6写入显式.log或 15–18 与大于 100 的单元并用完整路径名。现象shell 或 beam 单元上的结果与手算不符厚度方向应力看起来丢了。根因shell/plane stress 的厚度方向应变增量必须由用户提供STRAININC(*,3)而 beam 的STRESSNEW(*,2)、(*,3)被假定为零、写入值被忽略shell 的横向剪切刚度必须由材料或截面定义否则无法计算。解法按第 1.6 节的分量表逐单元类型改写shell 必须在材料里给弹性横向剪切模量或在截面里给横向剪切刚度。现象开启 GPU 加速后含 VUMAT 的作业没有加速或无法启动 GPU 路径。根因GPU 加速的生效范围与用户子程序的兼容性由官方文档口径决定2026 FD01 起为 Linux 上的 NVIDIA GPUFortran 用户子程序不在 GPU 计算路径内这一点须以官方文档为准。解法先跑内置材料的对照版本确认 GPU 收益基线再决定把本构换成内置材料或用 CPU 路径跑 VUMAT。四、动手练习单元素弹性段对照跑第 2.3 节的输入文件把材料从*USER MATERIAL改成内置*ELASTIC其余参数完全相同重跑一次。判定标准两次同一增量下的S11相对误差 1%vumat_check.log中出现稳定时间增量的数值把PROPS(1)从 21000.0 改成 210000.0硬 10 倍后稳定时间增量相对第一次缩小到约1/sqrt(10) ≈ 0.316倍允许 ±10% 偏差。径向返回的正确性判定在同一单元上施加单轴应变增量 0.001第一段再施加 0.002第二段SIGY0取 200.0、ETAN取 700.0。判定标准第一段结束后SDV(1)必须为 0弹性段第二段后SDV(1)与手算DPEQ (SM - YH)/(3G)一致相对误差 2%应力的偏量部分必须落在屈服面上SQRT(3.0*J2)等于SIGY0 ETAN*SDV(1)误差 2%。向量化代价测量可观察输出把第 2.1 节的单点IF (SM .LE. YH)改成整块判断版本先扫一遍块内是否全部弹性再执行无分支的更新路径同一模型同一增量数跑两次记录墙钟时间。判定标准生成两个计时记录例如 CSV 含 3 列版本、增量数、秒数行数 ≥ 2分支上提后的墙钟时间必须小于单点分支版本同一模型同一增量数且两次SDV(1)的最大相对差 0.5%证明改动没有改变结果。五、小结与下一篇预告本篇把 VUMAT 还原成一份契约它按块接收NBLOCK个材料点的数据必须回写STRESSNEW与STATENEW不提交任何材料雅可比显式对它的两个要求是物理上正确的应力更新和可用于时间增量的有效模量仅EFFMOD生效时用户写入的模量参与计算。没有收敛问题不等于没有精度问题条件稳定的算法靠步长控制误差本构过硬会按1/sqrt(模量比)压低稳定增量单精度下历史变量会累积舍入误差条件分支会破坏向量化。BEGSUB/KSIZE与JINFOARRAY这一层是内存布局契约越界写入污染的是求解器自己的数据。下一篇第 15 篇用户子程序调试方法DQN 与数值检查会把本篇的单元素对照测试做成可复用的调试工作流单精度与双精度的一致性、STRESS与DDSDD的一致性检查UMAT 侧、显式里用 Fortran 单元 6 与XPLB_EXIT的安全退出路径以及 DQN 参与显式求解的边界条件第 15 篇详述。本篇认知问题回显FAQQ1VUMAT 与 UMAT 的接口契约到底差在哪调用粒度、必须回写的数组、是否必须提供材料雅可比矩阵material Jacobian记为DDSDD、剪切分量按张量分量还是工程分量存储AVUMAT 按块调用一次传入NBLOCK个材料计算点必须回写STRESSNEW与STATENEW接口里没有DDSDDUMAT 逐个积分点调用且必须提供DDSDDEUMAT 的剪切分量是工程分量VUMAT 的剪切分量是张量分量。Q2为什么 VUMAT 没有收敛问题却仍有精度与稳定性问题显式的稳定时间增量与 VUMAT 里那份本构刚度是什么关系A显式不参与 Newton 迭代误差只能靠步长控制所以没有收敛失败稳定时间增量由膨胀波速sqrt((K4G/3)/ρ)与单元特征长度决定本构过硬会把稳定增量压到约1/sqrt(模量比)倍且只有指定EFFMOD时DTARRAY(2:NBLOCK1)与DTARRAY(NBLOCK2:2*NBLOCK1)里用户写入的模量才参与时间增量计算。Q3RHS、BIGH、BIGD、DINC、KMAX、K这些名字现在还在 VUMAT 参数里吗它们与当前接口的STRAININC、RELSPININC、JINFOARRAY怎么对应A这些是早期显式文档中的参数名当前接口用STRAININC应变增量、STRETCHOLD/STRETCHNEW、DEFGRDOLD/DEFGRDNEW、NBLOCK与循环索引KM、JINFOARRAY积分点/层/截面点与单元编号起始位置承接逐条对应关系以官方文档为准。Q4状态变量solution-dependent state variables在显式里怎么传递删除标志为什么必须用状态变量而不是返回值向量化KMAX/K/KSIZE那类块循环与性能的关系以及什么样的写法会把向量化打回逐点循环A状态变量靠STATEOLD读入、STATENEW回写逐增量传递删除标志必须占用一个状态变量1 有效、0 删除且不可重新激活向量化要求整块执行同一串算术运算只对块内部分材料点成立的IF分支会退化为逐点循环。Q5VUMAT 能挂哪些单元类型哪些单元类型不传入拉伸张量与变形梯度显式 GPU 加速2026 FD01 起对 Fortran 用户子程序意味着什么Ashell/plane stress 时NDIR3, NSHR1用户必须给厚度方向应变增量STRAININC(*,3)beam 与 pipe 拿不到拉伸张量与变形梯度2026 FD01 起 Abaqus/Explicit 在 Linux 的 NVIDIA GPU 上支持 GPU 加速其对 Fortran 用户子程序的影响以官方文档为准。
阅读完成 · 觉得有帮助？
