DEFORM3D二次开发：USRMSH实现热锻组织演化与晶粒预测

简介：围绕DEFORM3D二次开发实现塑性成形过程组织演化模拟的文档，面向材料成形、锻造工艺及有限元二次开发方向的学习者与工程人员，重点回应商用软件偏重宏观热力计算、微观组织预测能力不足的问题。文档结合刚粘塑性热力耦合有限元与动态再结晶、静态再结晶及晶粒长大模型，说明组织模型离散、迭代步计算和USRMSH子程序变量修改思路，并以20CrMnTi钢摇臂轴镦锻为算例，讨论温度场、应变、应变速率与晶粒尺寸的对应关系。资源包共1个doc文件，大小约158KB，内容紧凑，适合作为论文式技术材料查阅；已有654人学习。读者可从中了解再结晶分数判断、临界应变与晶粒长大衔接等处理方式，也可对照算例理解热力参数分布与内部晶粒度变化之间的联系，为热成形工艺制定和二次开发提供参考。

1. 云图好看、性能不达标：从 DEFORM3D 二次开发补上组织演化这一环

锻件下模架之后做金相，晶粒粗大，力学性能压不住——可回去翻模拟结果，温度场、等效应变场一片均匀，找不出毛病。这类翻车在热锻工艺开发里太常见了，问题出在商业软件默认只解宏观场量，晶粒尺寸要么不给，要么内置一个和你材料对不上的通用模型。这份资源盯的就是这个缺口：在 DEFORM3D 的热力耦合刚粘塑性有限元框架上做二次开发，通过用户子程序 USRMSH 把 20CrMnTi 的动态再结晶、静态再结晶和晶粒长大模型塞进每一次迭代，让再结晶分数和晶粒尺寸跟着温度、应变、应变率一起算出来。它适合已经会用 DEFORM 跑基础成形模拟、想往组织预测走一步的工艺开发，也适合做材料成形数值模拟的人拿它当二次开发的第一份完整样板——代码结构不复杂，但链路是通的。

2. 架构拆解：刚粘塑性热力耦合框架怎么容纳一套 Avrami 组织模型

2.1 为什么把组织模型挂在刚粘塑性有限元上

体积成形（镦粗、挤压、模锻）的塑性变形量远大于弹性变形量，弹性卸载那点回弹对晶粒演化的影响可以忽略，所以用刚粘塑性本构就够，省掉弹性迭代，单步收敛快、内存占用低，这是选它的第一层原因。第二层原因是热力耦合：塑性功转化的热量直接进温度场，而温度又反过来决定再结晶动力学里的 Zener-Hollomon 因子，两者必须同一个时间步内交替求解。组织模型不是独立的一层，它是挂在热力耦合迭代上的一个附加计算环节——主程序每个迭代步把温度、等效应变、等效应变率、塑性应变增量算完，回调一次用户子程序，你在里面把组织量更新完写回去，下一次迭代继续。这种耦合方式的好处是组织计算不改动求解器的刚度矩阵和残差组装，出问题也只影响附加量，不会让整个求解崩掉，代价是组织演化的时间积分精度受限于主程序的时间步长。

常见做法是在前处理阶段就把材料的热物性参数（导热系数、比热、密度随温度变化）和流变应力曲线一起建进材料库，20CrMnTi 这种结构钢商业库里不一定有完整数据，得用热力模拟试验的数据自己回归。流变应力曲线不对，应变率场就是错的，后面所有组织计算都建立在错误输入上，这一步偷懒没有后悔药。

2.2 组织模型的数学骨架：Avrami 分数方程加晶粒尺寸经验式

热变形过程中的组织演化拆成三件事：动态再结晶（变形过程中发生）、静态再结晶（道次间或卸载后发生）、晶粒长大（再结晶完成后的保温或缓冷阶段）。三者的数学形式都是 Avrami 型方程回归出来的准经验模型，形式大致如下，具体系数靠试验数据拟合：

动态再结晶分数

X_drx = 1 - exp[ -0.693 · ((ε - εc) / ε0.5)^k ]

其中 εc 是临界应变，ε0.5 是再结晶分数达到 50% 时的应变，两者都由 Zener-Hollomon 因子 Z = ε̇·exp(Q/(R·T)) 和初始晶粒尺寸 d0 决定：

εc = a_c · d0^h_c · Z^m_c ε0.5 = a_5 · d0^h_5 · Z^m_5

再结晶后的晶粒尺寸同样是经验式：

d_drx = a_d · d0^h_d · ε^n_d · ε̇^m_d · exp(Q_d/(R·T))

晶粒长大用经典的幂律时间关系：

d^n = d0^n + k_g · t · exp(-Q_g/(R·T))

静态再结晶的形式和动态再结晶几乎一样，区别在于驱动力用道次间的残余应变，时间用道次间隔时间，且只有在前一道次变形停止后才启动。这套模型在文献里被反复验证过，参数依赖材料，20CrMnTi 和别的钢种不能混用，这也是为什么必须做二次开发而不是直接用软件内置模型。

2.3 一次求解的完整数据链路

从主程序回调进 USRMSH 开始，到结果写回数据库，中间的数据流是这样的：主程序把当前迭代步的单元温度 TEMP、温度增量 DTMP、等效应变率 EFEPS、塑性应变增量 TEPS 传进来；子程序第一步判断当前单元处于哪个阶段——动态再结晶分数小于 0.95 时先看等效应变有没有超过临界应变，超了就算动态再结晶；分数超过 0.95 说明再结晶基本完成，转入晶粒长大计算；如果当前步是道次间或者变形已停止，走静态再结晶分支。算完之后，通过 VF、DFV、DFN、GRNSZ 这几个数组把再结晶分数、分数增量、静态再结晶分数和晶粒尺寸写回去，主程序下一次写数据库时就会把它们当作场变量输出。

这里容易忽略的一点是状态量的持久化。再结晶分数不是每步从零开始算的，它依赖上一步的值做累积；晶粒尺寸也是。DEFORM 主程序在网格重划分时会把宏观场量插值到新网格上，但你自己定义的普通数组不会被自动迁移，必须用 USRVE 这类用户状态变量槽位把历史值存住，否则一重划分组织场就清零，曲线会出现莫名其妙的跳变。

3. 跑通一次摇臂轴镦锻：网格、子程序编译与求解提交

3.1 工程文件与目录怎么摆

一个能复现的最小工程，文件结构大致是下面这样。不同版本的目录名和启动脚本会有差异，按自己安装路径对齐即可。

把模型参数放在外部.dat里而不是硬编码进 Fortran，是我强烈建议的做法。换钢种时只改参数文件，不用重新编译，批量试算的时候差别很大。

3.2 前处理：对称面、网格与边界条件

摇臂轴这类长轴类锻件，几何上存在平面对称，取一半建模能把单元数砍掉一半。参考做法是坯料按六面体单元离散，单元数在数千到数万量级，上下模各自独立划网格，模具当成刚性体处理，只参与接触和传热计算。文献里给出的坯料数字是五万余单元、七千余节点，这个比例按六面体八节点单元的拓扑关系反推其实对不上，更像是排版笔误，实际建模别照抄数字，按自己几何的尺寸重划，然后在关键变形区（飞边、圆角过渡处）做局部加密——晶粒尺寸的梯度在那些位置最陡，网格粗了会算出一片糊。

工序设置上，一步成形就把上模下行速度、坯料初始温度、模具预热温度、环境换热系数定死。摩擦用剪切摩擦模型，摩擦因子对等效应变分布影响直接，进而影响再结晶分数，这个参数别随手填 0.3，有条件用圆环镦粗试验标一下。求解控制里时间步长建议按模具行程的百分比给，每步行程控制在最小单元尺寸的十分之一以内，太大会让组织量的时间积分失真。

3.3 编译用户子程序并提交计算

Fortran 子程序需要编译成求解器能加载的动态库。下面是 Linux 环境下的典型流程，Windows 下用对应的编译器生成 DLL，思路一致。

# 1) 进入用户子程序目录，路径按自己安装位置改 cd $DEFORM_DIR/usr/def_usr # 2) 备份原始模板，改动前的最后一份后悔药 cp def_usr.f def_usr.f.bak # 3) 编译。makefile 里一般已经封装好编译器和优化选项， # 手工等价命令大致是 -c 出目标文件，-shared 出动态库 make -f makefile.linux # 等价于： # gfortran -c -O2 -fPIC def_usr.f # gfortran -shared -o def_usr.so def_usr.o # 4) 提交计算，把数据库和刚编出来的库一起交给求解器 $DEFORM_DIR/bin/def_solver -b shaft_billet.DB -usr ./def_usr.so > run.log 2>&1 & # 5) 盯日志，前几步就能看出子程序有没有被真正调用 tail -f run.log | grep -i -E "usr|error|warn"

参数说明几个关键点：-b后面跟数据库文件，工序、材料、边界全在里面；-usr指定用户子程序库，路径写错求解器不会报错退出，而是静默使用默认空子程序，算完一看组织场全是零——这个坑后面还会提。run.log里如果出现USRMSH called之类的打印（需要你自己在子程序里加），说明回调链路通了；一条都没有，先去查库有没有被加载，别急着改公式。

3.4 后处理：三张图看懂结果

结果判读的顺序有讲究，别一上来就看晶粒尺寸云图。先看温度场，确认坯料心部和表层温差在合理区间、模具预热温度有没有生效；再看等效应变和应变率场，确认大变形区落在你想让它落的位置，飞边桥口是不是过度变形；最后才看再结晶分数和晶粒尺寸云图，两相对照才能判断某个区域的细晶是"变形充分再结晶"带来的，还是"温度低没长大"侥幸留下的。这两种细晶的稳定性完全不同，前者可复现，后者换个批次料就没了。

导出路径上的数据做定量对比更靠谱。DEFORM 支持把场变量沿一条用户定义的路径导出成 CSV，取轴心线和表层两条路径，画晶粒尺寸随位置的曲线，比看云图颜色直观得多：

import pandas as pd import matplotlib.pyplot as plt # 从后处理导出的路径数据，列名按实际导出表头改 df = pd.read_csv("path_centerline.csv") fig, ax1 = plt.subplots(figsize=(7, 4)) # 晶粒尺寸随距轴心距离的变化 ax1.plot(df["distance_mm"], df["grain_size_um"], "o-", label="grain size") ax1.set_xlabel("distance from center (mm)") ax1.set_ylabel("grain size (um)") ax2 = ax1.twinx() # 叠一条温度曲线，看晶粒粗化是不是跟着高温区走 ax2.plot(df["distance_mm"], df["temp_C"], "--", color="gray", label="temperature") ax2.set_ylabel("temperature (C)") plt.title("grain size vs temperature along centerline") plt.tight_layout() plt.savefig("centerline_check.png", dpi=150)

这段脚本本身没什么技术含量，价值在于把"云图好看"变成"曲线可核对"：晶粒峰值出现在哪个温度区间、对应多大的应变率，一眼能读出来。工艺窗口就是从这个对应关系里反推的。

4. USRMSH 逐段拆解：一次迭代里组织演化的四步计算

4.1 参数表：哪些是输入，哪些是输出

先把子程序声明里跟组织模拟相关的参数捋清楚，剩下的接触、磨损、相变槽位这次用不上，但声明不能删，参数个数必须和主程序约定的一致，少一个就链接失败。

四个输入量里，真正驱动组织计算的是温度、等效应变率和等效应变，温度增量只用于判断升降温方向。GRNSZ 是双向的——进来时是上一步的晶粒尺寸，出去时是本步算完的新值，写错了就会每步都从同一个初值重新起算，曲线变成锯齿。

4.2 第一步：把热力场量取出来

进入单元循环之前先把材料常数准备好。参数文件只在第一次调用时读，读进 COMMON 块缓存，避免每个时间步都做一次磁盘 IO——一次镦锻模拟动辄几百上千步，每步几千个单元，IO 开销能把总时长拉高一个量级。

C ---- 材料常数缓存：首次调用读文件，之后走 COMMON ---- REAL*8 QACT, RGAS, ACRIT, NCRIT, A5, N5, KAV, DGROW COMMON /MATPAR/ QACT, RGAS, ACRIT, NCRIT, A5, N5, KAV, DGROW LOGICAL LINIT COMMON /MATINI/ LINIT IF (.NOT. LINIT) THEN OPEN(11, FILE='mat_20crmnti.dat', STATUS='OLD') READ(11,*) QACT ! 再结晶激活能 J/mol READ(11,*) ACRIT ! 临界应变系数 READ(11,*) NCRIT ! 临界应变 Z 指数 READ(11,*) A5 ! 50% 再结晶应变系数 READ(11,*) N5 ! 50% 应变 Z 指数 READ(11,*) KAV ! Avrami 指数 k READ(11,*) DGROW ! 晶粒长大速率常数 CLOSE(11) RGAS = 8.314D0 ! 气体常数 LINIT = .TRUE. END IF

这段逻辑很朴素但很关键：LINIT用 COMMON 保存，保证整个求解过程只读一次文件。读入顺序必须和.dat文件里的行序完全对应，中间插一行注释行就会错位，后续所有计算结果都是垃圾——参数文件要么严格按行读，要么改成带关键字的KEY = VALUE格式自己解析，后者更稳，代码多几十行而已。

4.3 第二步：动态再结晶的判据与分数更新

单元循环里第一件事是把该单元的场量取成局部变量，然后算 Z 因子和临界应变。文献里的流程是：先看当前再结晶分数是否小于 0.95，小于则判断等效应变是否超过临界应变，超过才启动动态再结晶计算；分数一旦超过 0.95，说明这一轮再结晶完成，改走晶粒长大分支。

DO IEL = 1, NUMEL C ---- 取局部量，避免长表达式反复索引 ---- TT = TEMP(IEL) + 273.15D0 ! 统一到 K，别用摄氏度算 exp ER = MAX(EFEPS(IEL), 1.0D-6) ! 应变率下限，防除零和 Z=0 EP = EPS(IEL) D0 = GRNSZ(IEL) XOLD = VF(1,IEL) C ---- Zener-Hollomon 因子 ---- ZN = ER * EXP(QACT / (RGAS * TT)) C ---- 临界应变与 50% 应变 ---- EPSC = ACRIT * D0**NCRIT * ZN**(1.0D0/3.0D0) EPS50 = A5 * D0**0.5D0 * ZN**(1.0D0/5.0D0) IF (XOLD .LT. 0.95D0) THEN C ---- 未完成再结晶：判据过了才更新分数 ---- IF (EP .GT. EPSC) THEN XNEW = 1.0D0 - EXP(-0.693D0 * & ((EP - EPSC) / MAX(EPS50, 1.0D-8))**KAV) C ---- 单调不减，且截断到 [0,1]，Avrami 公式会溢出 ---- XNEW = MIN(MAX(XNEW, XOLD), 1.0D0) DVF(1,IEL) = XNEW - XOLD VF(1,IEL) = XNEW END IF END IF

几个细节值得说清楚。温度先加 273.15 转开尔文，因为激活能 Q 是按开尔文标定的，用摄氏温度算exp(Q/RT)会偏出几倍甚至几十倍，结果要么再结晶分数恒为 1，要么恒为 0。MAX(EFEPS,1e-6)是防除零和防 Z 因子取零的保护，刚体接触点上应变率可能是零。XNEW = MIN(MAX(XNEW, XOLD), 1.0)这一行是血泪经验：Avrami 公式在应变刚过临界值时增量很小，数值噪声会让新值比旧值还低，分数出现"倒退"，后处理曲线上表现为抖动；截断到 [0,1] 又防止指数项溢出把分数顶到 1.2 这种非法值。

4.4 第三步：静态再结晶与晶粒长大

动态再结晶完成后转晶粒长大，另外在道次间隙（ER接近零且CURTIM跨过变形结束时刻）走静态再结晶分支。两者都是时间驱动的，用DTMAXC或相邻两步的时间差做积分步长。

C ---- 分支一：再结晶完成后，晶粒长大 ---- IF (XOLD .GE. 0.95D0) THEN DDT = MAX(DTMAXC, 1.0D-3) ! 时间步，秒 GRNSZ(IEL) = (D0**DGROW + & 1.0D-2 * DDT * & EXP(-QACT / (RGAS * TT)))**(1.0D0/DGROW) END IF C ---- 分支二：道次间静态再结晶，靠应变率骤降判定 ---- IF (ER .LT. 1.0D-4 .AND. USRVE(1,IEL) .GT. 0.5D0) THEN XSRX = 1.0D0 - EXP(-0.693D0 * & (USRVE(1,IEL) / MAX(EPS50, 1.0D-8))**KAV) DFN(1,IEL) = MAX(DFN(1,IEL), XSRX) ! 静态分数只增不减 END IF C ---- 状态量落盘，网格重划分时才不会丢 ---- USRVE(1,IEL) = EP ! 累计应变快照 USRVE(2,IEL) = GRNSZ(IEL) ! 晶粒尺寸快照 USRVE(3,IEL) = VF(1,IEL) ! 再结晶分数快照 END DO

晶粒长大这一段用的是幂律式的简化积分，DGROW是长大指数，一般取 2 到 4 之间。注意它只在再结晶分数超过 0.95 之后才执行，否则一旦升温就开始粗化，前面所有细化努力全白费。静态再结晶分支用应变率骤降作为触发条件，这在单道次镦锻里不明显，多道次锻造、有传送等待时间的工况下必须打开，否则算出来每道次晶粒都只细化不粗化，和实测对不上。

最后把EP、晶粒尺寸、再结晶分数写进USRVE的三个槽位。这三个槽位需要在辅助变量对象里注册，注册之后重划分网格会自动做场量映射，不注册就是普通数组，一切网格全丢。

4.5 第四步：写回与状态量持久化

写回这一步没有花哨操作，就是把算好的值放回主程序约定的数组位置。真正值得盯的是 NPHASE 这个维度——多相材料里第二相、第三相的槽位是分开的，VF(1,IEL)和VF(2,IEL)含义完全不同，索引写错不会报错，只会静默算出一个看起来合理的错值。单相材料也建议在代码里用命名常量IPH_DRX = 1代替裸数字 1，改起来省心。

5. 避坑记录：二次开发里五个必翻车的地方

5.1 现象：算完了，组织场全是零

原因有两种。最常见的是动态库根本没被加载，-usr路径写错、文件权限不对、编译器和求解器位数不匹配（32 位库喂给 64 位求解器），这三种情况求解器都不一定报错。另一种是子程序被调用了，但数组下标或 NPHASE 维度搞错，写回的位置没有任何单元在用。排查办法很土但有效：在子程序第一行加一句无条件打印，带上 IEL 和计算值，跑三步就停，看日志里有没有输出、值是不是合理。先验证调用链，再验证算法。

5.2 现象：网格重划分之后晶粒尺寸曲线出现断崖

原因是历史量没进 USRVE。DEFORM 在单元畸变超限时会触发局部重划分，宏观场量（温度、应变、应力）会自动插值到新网格，但你的普通数组不在映射清单里，重划分那一刻组织场归零，下一步又从初始晶粒尺寸重新起算。解决办法是把再结晶分数、晶粒尺寸、累计应变全部注册成用户状态变量，每个单元每步都写快照。代价是数据库文件变大，几百兆很正常，属于必要成本。

5.3 现象：再结晶分数恒为 1 或者恒为 0，怎么调参数都没用

温度单位和激活能量纲不一致。温度用摄氏度、激活能用 kJ/mol、气体常数用 8.314，三种单位混着算，exp(Q/RT)直接溢出到无穷大或下溢到零。统一的做法是全流程 SI：温度开尔文、激活能 J/mol、气体常数 8.314 J/(mol·K)、应变率 1/s。改完之后顺手在日志里打印一次 Z 因子，量级在 10^10 到 10^18 之间算正常，偏出去两三个数量级一定是单位错了，别去怀疑模型公式。

5.4 现象：多道次模拟里晶粒只细化不粗化

静态再结晶分支从来没被触发。触发条件是应变率骤降，但很多工况下模具回程和坯料转移阶段被简化掉了，时间轴上根本没有"空载"区间，程序看不到等待时间。解决办法是在工序定义里显式插入道次间冷却或者传热工序，把间隔时间给出来；或者在子程序里用CURTIM对比预设的道次时间表，主动判定当前处于间隙阶段。这根时间轴不对，静态再结晶和晶粒长大的计算结果全是摆设。

5.5 现象：单步计算时间从几秒涨到几分钟

两种可能。一是把参数文件的读取放在了单元循环里，几千单元乘几百步，磁盘 IO 直接把性能吃干净——LINIT缓存必须做。二是网格量级失控，文献里那种五万单元、七千节点的组合在六面体拓扑下本来就不成立，照着抄会让求解器在接触判断上疯狂回溯。按对称面取半、关键区局部加密、非关键区放粗，把总单元数控制在万级以内，单步时间能回到可接受范围。二次开发阶段建议先用一两千单元的粗模型把链路跑通，再换精细网格出正式结果。

6. 换成自己的材料：参数标定、模型替换与结果验证

拿到这份资源之后真正要干的活，是把 20CrMnTi 的模型换成你自己的钢种。改造路线按投入产出比排，第一步永远是参数外置，第二步才是模型形式。

先做参数化。把 COMMON 里硬编码的常数抽到外部.dat，按KEY = VALUE格式解析，一个钢种一个文件，换材料只改文件名。这一步半天就能做完，但省下来的是后面几十次重新编译。紧接着把 Avrami 的指数 k 和临界应变系数从固定值改成可配置——不同文献回归出来的 k 值差异不小，取值 1.5 到 2.5 之间结果能差出一倍晶粒尺寸，必须允许调。

然后是验证，验证方法比模型本身更重要。实验室热力模拟试验机上按几个典型变形条件（不同温度和应变率组合）做单道次压缩，快冷保留组织，测晶粒尺寸，和同工况的模拟结果对比。误差带按 ±15% 看，超出这个范围八成不是模型形式的问题，而是流变应力曲线或者摩擦参数先错了——宏观场量错了，组织模型再准也救不回来。这一轮对完，再去做多道次验证，把静态再结晶和晶粒长大那段补上。

工艺窗口的扫描值得单独说一句。把加热温度、变形量、变形速率、模具预热温度做成四因素正交，用.key文件脚本批量生成算例，跑完自动提取轴心线的晶粒尺寸均值，画成响应面。这一步不需要任何新的算法，纯粹是把已经跑通的链路用起来，但它产出的东西才是工艺员真正要的——不是一张漂亮的云图，而是一句"这个温度区间加这个变形量，晶粒能稳在某个范围内"。

从那以后我每次改子程序，都强制先跑一个两三千单元的小算例，只跑五步，把温度、应变率、再结晶分数、晶粒尺寸四个量在日志里打出来对一遍，确认量级和趋势都对，再上正式模型。翻车翻得最多的从来不是公式写错，而是链路没通、单位没对齐、状态量没存住这类"低级"问题，而它们在小算例上暴露得最快。希望帮到你。

本文还有配套的精品资源，点击获取