DEFORM3D二次开发：USRMSH实现热锻再结晶与晶粒尺寸预测

简介：围绕 DEFORM3D 二次开发实现塑性成形过程组织演化模拟的技术文档，面向材料成形、锻造工艺与有限元二次开发方向的研究生、工程师及科研人员，重点解决商业软件仅能完成宏观热力耦合计算、缺乏微观组织预测能力的问题。文档以 20CrMnTi 钢摇臂轴镦锻热成形为实例，梳理动态再结晶、静态再结晶与晶粒生长的数学模型耦合流程，并给出 USRMSH 子程序接口中各参量的含义与修改方式，配合温度场、应变与应变速率、再结晶分数及晶粒尺寸的模拟结果加以说明。包内共 1 个 doc 文件，约 158KB，为完整的论文式技术文稿，含摘要、模型建立、实例分析与参考文献，可直接作为二次开发入门与工艺参数分析的参考。目前已有 655 人学习下载。

1. 从一张"只有温度没有晶粒"的云图说起

热锻件出模后硬度散差大、后续渗碳淬火变形超差，追到根上常常是晶粒没控住。可打开 DEFORM3D 跑完一套镦锻，屏幕上只有温度、等效应变、应变速率、流变应力这几张云图，晶粒尺寸那一栏是空的——软件默认只解宏观场量，组织演化那一路得自己接上去。手上这份资料讲的正是这条接法：基于 DEFORM3D 二次开发，用 USRMSH 用户子程序把 Avrami 形式的动态再结晶、静态再结晶和晶粒长大模型嵌进热力耦合刚粘塑性有限元循环，再用 20CrMnTi 钢摇臂轴镦锻做算例验证。它适合已经能独立跑通一套等温锻造模拟、想把组织预测塞进工艺评审的成形工艺工程师，也适合做塑性成形课题、需要拿出晶粒演化结果的研究生。下面按"为什么改、改哪里、参数怎么设、哪里会翻车"四层往下拆。

2. 商用成形软件的组织模拟边界：为什么非动子程序不可

2.1 宏观场量求解与微观组织演化的耦合关系

热变形过程中的组织变化不是一条连续曲线，而是分阶段切换的：变形进行时发生动态再结晶（DRX），变形停止后、温度还高时发生静态再结晶（SRX）和亚动态再结晶，随后进入晶粒长大。这三个阶段用的模型形式都差不多是 Avrami 型，但触发条件和时间基准完全不同，一个绑定应变，一个绑定保温时间。通用成形软件不区分这些阶段，因为它根本没有组织状态量。

关键在于，宏观求解器输出的温度 T、等效应变 ε、等效应变率 ε̇、以及时间 t，恰好就是组织模型的全部输入。缺的不是数据，缺的是一个"积分器"——在每个求解步末尾，把当前单元的 T、ε、ε̇ 取出来，代入经验模型算一次增量，再把新的再结晶分数和晶粒尺寸写回求解器的状态数组，让它跟着几何、温度一起被更新、被后处理读取。这就是二次开发要干的事，工作量集中在模型离散和变量管理，不在算法本身。

需要提前说清一个边界：这份资料做的是单向耦合，组织演化只读取宏观场量，不反过来影响流变应力。如果想让软化后的流变应力回灌到本构里做双向耦合，得在子程序里同时改写等效应力相关量，收敛性会立刻变成一个棘手问题，通常是另一个量级的工程。先做单向，把晶粒尺寸算对，是更稳的路径。

2.2 刚粘塑性热力耦合的求解循环与插入点

体积成形用的刚粘塑性有限元，基本假设是材料不可压缩、忽略弹性变形，用 Markov 变分原理建立泛函，体积不可压缩约束用罚函数或拉格朗日乘子处理。每个时间步的求解顺序大致是：在当前温度场下用 Newton-Raphson 迭代求速度场，收敛后更新节点坐标和等效应变、应变率，再解一次热传导方程更新温度场，然后进入下一个时间步。

USRMSH 的调用位置就在这个循环的收敛步之后。理解这一点很重要，因为它决定了三件事。第一，组织变量不是独立积分出来的，它跟随求解步走，步长就是组织模型的时间分辨率，步长给大了，再结晶分数会一步跳过去。第二，温度、应变率都是当前步的瞬时值，不是单元历史平均值，写增量公式时要自己维护累积量。第三，网格重划分会打断状态变量的连续性——新单元从哪里继承旧单元的晶粒尺寸，这一步没人替你处理，必须自己写映射，这是后面避坑章要重点讲的一条。

刚粘塑性还有一个隐含约束：材料不可压缩意味着体积应变恒为零，等效应变严格由塑性应变增量累积得到。所以子程序里累加等效应变时，一定要用塑性应变增量那一项，而不是从位移场反算的应变，否则在网格畸变区域会累积出虚假应变，直接放大再结晶分数。

2.3 组织模型的状态变量设计：VF、DFV、DFN、GRNSZ

子程序能读写的组织相关量有四组，理解它们的语义差别比记住名字更重要。

VF(NPHASE, *) 是再结晶分数，NPHASE 是相数，单相钢取 1 即可。DFV(NPHASE, *) 是再结晶分数增量。DFN(NPHASE,) 是静态再结晶分数。GRNSZ() 是晶粒尺寸。前三个的维度都按单元索引铺开，最后一个通常也按单元或节点铺开。

这里有个我踩过的判断：不要直接覆盖 VF 的绝对值，优先写 DFV 增量。原因是部分版本的求解器内部对 VF 有自己的更新策略，你在子程序里写死的绝对值有可能在后续步被求解器的逻辑覆盖掉，表现就是分数来回跳。写增量、让求解器去累加，鲁棒性明显更好。反过来，GRNSZ 一般没有内部逻辑，直接写绝对值没问题。

除了这四个专用量，每个单元还有一组用户状态变量 USRVE(NUSRVE, *)，这才是真正该用来存中间量的地方：累积等效应变、初始晶粒尺寸 d0、上一次发生静态再结晶的时刻、当前所处的阶段标志。把这些存进用户状态变量而不是每次重算，一方面省时间，另一方面重启动时能接得上。变量个数在前处理里分配，分配少了后面加不进去，我一般会一次留够 10 个以上。

2.4 三条技术路线的取舍

对绝大多数做工艺的团队，第二条是性价比最高的。自带模块的问题不在于不准，而在于它的系数是通用拟合，换一个炉批号的 20CrMnTi 就可能偏出去；而且它的阶段切换判据固定，你没法按自己的热压缩数据改临界应变条件。

3. 把再结晶模型落进 USRMSH：参数表、状态变量与代码骨架

3.1 USRMSH 形参速查与读写边界

子程序的形参有五十多个，真正用到的不到十个，但每一个的维度和排布方式必须核对清楚。下面这张表是按资料里的声明整理的，实际动手前建议先拿一个单单元模型把所有数组的 size 打印一遍验证，别凭记忆写循环上界。

顺序上的一个硬性提醒：这份资料给出的声明里同时出现了 EFSTS 和 EFEPS，正文说明 EFEPS(NUMEL) 是等效应变率、TEPS(NUMEL) 是塑性应变增量。不同版本之间形参顺序存在差异，抄代码之前必须在安装目录里找到原始的 def_usr.f，逐字对一遍声明，用关键字而不是位置去认参数。我见过有人直接把网上的代码段贴进去，编译能过，结果温度读进了应变率，晶粒尺寸算出来全场一个值，排查了半天。

3.2 动态再结晶判据与 Avrami 增量计算

模型形式是准经验的 Avrami 型，系数必须来自你自己的热力模拟试验（热压缩或扭转），资料里没有给 20CrMnTi 的具体系数，我下面用 A50、N50、M50、QACT 这类占位符，你替换成自己回归的值。别直接套文献里的数字，不同炉批、不同变形区间的差别足以让结果不可信。

离散化的核心是把连续方程写成增量形式，在每一步用当前步的 T、ε̇ 和累积应变算一次。

!------------------------------------------------------------------ ! 在每个收敛步末尾调用：遍历单元更新 DRX 分数与晶粒尺寸 ! 自由格式书写，若你的 def_usr.f 是固定格式，注意第 6 列续行 !------------------------------------------------------------------ DO IEL = 1, NUMEL TMP = TEMP(IEL) ! 温度，K ER = EFEPS(IEL) ! 等效应变率，1/s DEPS = TEPS(IEL) ! 本步塑性应变增量 EPS = USRVE(1,IEL) + DEPS ! 累积等效应变，存用户状态变量 USRVE(1,IEL) = EPS D0 = USRVE(2,IEL) ! 初始晶粒尺寸，前处理里预置 ! 变形量太小或温度太低，直接跳过，避免 exp 溢出 IF (EPS .LT. 1.0E-4 .OR. TMP .LT. 773.0) GOTO 100 ! 临界应变判据：eps_c = A_C * eps_p，eps_p 为峰值应变 EPSC = A_C * EPS_P(EPS, ER, TMP) ! 峰值应变同样由回归式给出 IF (EPS .LT. EPSC) GOTO 100 ! 未到临界，不发生 DRX ! 50% 再结晶时间：eps50 = A50 * eps^n * er^m * exp(Q/(R*T)) EPS50 = A50 * (EPS**N50) * (ER**M50) & * EXP(QACT / (RGAS * TMP)) ! Avrami：X = 1 - exp(-0.693 * (eps/eps50)^k) XNEW = 1.0 - EXP(-0.693 * ((EPS / EPS50)**AK)) IF (XNEW .GT. 0.995) XNEW = 0.995 ! 留一点余量给长大分支 IF (XNEW .LT. 0.0) XNEW = 0.0 ! 增量写法，交给求解器累加，避免覆盖内部状态 DFV(1,IEL) = XNEW - USRVE(3,IEL) USRVE(3,IEL) = XNEW ! DRX 晶粒尺寸：d = Ad * d0^h * eps^n * er^m * exp(-Qd/(R*T)) DDRX = AD * (D0**HD) * (EPS**ND) * (ER**MD) & * EXP(-QDACT / (RGAS * TMP)) ! 混合法则：d = X*d_drx + (1-X)*d0，未再结晶部分保留原始尺寸 GRNSZ(IEL) = XNEW * DDRX + (1.0 - XNEW) * D0 100 CONTINUE END DO

代码逻辑上分四层。第一层做量纲和数值保护，温度低于再结晶门槛或变形量极小时直接跳过，否则 EXP 里的指数项很容易溢出成 NaN，一崩就是全场崩。第二层是临界应变判据，这一步决定了 DRX 什么时候启动，临界应变一般取峰值应变的 0.6 到 0.8 倍，具体系数从热压缩曲线的峰值点反推。第三层算 50% 再结晶时间，这是 Avrami 的时间尺度，它对温度的敏感度极高，QACT 取 300 kJ/mol 量级时，温度差 20 K 就能让 eps50 差一倍，所以温度场本身的精度直接决定晶粒尺寸的可信度。第四层做晶粒尺寸混合，用再结晶分数加权，未再结晶的单元保留原始晶粒尺寸，这样云图上才能看出"边角没变形、晶粒没细化"的真实分布。

参数上要盯住三个：QACT 和 QDACT 的单位必须和 RGAS 匹配，统一用 J/mol 配 8.314；ER 在模具刚接触、单元还没真正变形时可能极端小，要做下限截断，我一般卡在 1e-4 量级；AK 是 Avrami 指数，通常 1 到 2 之间，取值越大曲线越陡，对步长越敏感。

3.3 静态再结晶与晶粒长大分支

变形结束后进入完全不同的物理过程：没有外加载荷，驱动力来自畸变能，时间基准从"应变"切换成"保温时间"。这段逻辑最好单独成块，用阶段标志隔开，别和 DRX 混在一个 IF 里。

!------------------------------------------------------------------ ! 卸载后分支：静态再结晶 + 晶粒长大 ! 判断依据：模具速度归零 / 接触状态解除，用阶段标志 USRVE(4,*) 锁存 !------------------------------------------------------------------ IF (ISTAGE .EQ. 2) THEN ! 已进入变形后阶段 DT = CURTIM - USRVE(5,IEL) ! 距进入该阶段的时间，s IF (DT .LE. 0.0) GOTO 200 ! 静态再结晶 Avrami：Xs = 1 - exp(-0.693*(t/t50)^ks) T50S = AS50 * (EPS**NS) * EXP(QSACT / (RGAS * TMP)) XS = 1.0 - EXP(-0.693 * ((DT / T50S)**AKS)) IF (XS .GT. 0.99) XS = 0.99 DFN(1,IEL) = XS - USRVE(6,IEL) USRVE(6,IEL) = XS ! 晶粒长大：d^n = d0^n + Kg * t * exp(-Qg/(R*T)) DN = GRNSZ(IEL)**NG DN = DN + KG * DT * EXP(-QGACT / (RGAS * TMP)) GRNSZ(IEL) = DN**(1.0 / NG) 200 CONTINUE END IF

静态再结晶的时间常数 T50S 里保留了应变项，因为变形量越大、储存畸变能越高，静态再结晶启动越快，这个耦合关系不能丢。晶粒长大用幂律形式，NG 通常取 2 到 4，取 2 时长大速率前期快后期慢，取 4 更平缓。这里最容易被忽略的是计时起点：USRVE(5, IEL) 必须在进入阶段 2 的那一步写入当前 CURTIM，而不是每步都刷新，否则 DT 永远是零，静态再结晶永远不启动。

阶段判据的写法有两种。一种是看模具速度是否归零，另一种是看单元是否还与模具接触。前者简单但对多道次锻造不友好，后者更贴近物理但接触状态的读取在各版本里不一致。我一般用"模具速度 + 上模位移不再增加"的组合判据，稳一些。

还有一个必须处理的细节：晶粒长大阶段温度还很高，如果步长给得大，长大项会一步跨过头，晶粒尺寸出现不合理的跳变。解决办法是在子程序里回写 DTMAXC，把允许步长压到合理范围，比如保温阶段最大步长不超过 1 s。这是子程序里少数几个能反向影响求解器的通道，用好了能省掉很多事后处理。

3.4 子程序编译链接与求解器挂接

写完 Fortran 只是第一步，能不能被真正调用是另一回事。DEFORM3D 的用户子程序集中放在 def_usr.f 里，编译生成自定义求解器，运行时必须显式选择这个自定义版本，否则跑的仍是发行版求解器，你的代码一行都不会执行。

# 进入用户子程序目录，先备份再改，别动发行版原文件 cd $DEFORM_DIR/user_sub cp def_usr.f def_usr.f.bak # 确认 USRMSH 在文件里的位置，只改这一段，其余保持空实现 grep -n "SUBROUTINE USRMSH" def_usr.f # 调用安装目录下的编译入口（不同版本脚本名不同，找带 makefile 的那个） ./compile_user.sh > build.log 2>&1 tail -n 30 build.log # 重点看有没有 undefined reference 与形参不匹配告警 # 运行前在前处理里指定使用自定义求解器，否则编译了也白编译

编译日志里有两类信息必须逐个看。一类是形参不匹配的警告，Fortran 对隐式接口比较宽容，参数个数对不上有时也能过，但运行时会读错内存。另一类是隐式类型警告，变量没声明就以 I-N 规则变成整型，浮点运算直接截断，算出来的晶粒尺寸永远是整数，这种错很难从结果上看出来。

挂接之后还有一步验证：在子程序开头加一行写文件或者打印，跑一个只有几步的算例，确认它确实被调用了、调用次数与步数对得上。这一步花不了五分钟，但能省掉后面几个小时的困惑。

3.5 步长、增量稳定性与收敛的相互牵制

组织模型是强非线性的指数型方程，对步长极其敏感。步长给到 0.5 s 时，一个 3 s 的变形过程只有 6 步，动态再结晶分数可能从 0 直接跳到 0.9，晶粒尺寸的云图看起来像被切了一刀。经验做法是把变形阶段的步长控制在模具每步压下量不超过最小单元尺寸的 1/10，镦锻这种大变形过程，往往要落到 0.05 s 甚至更小。

代价是计算时间。52111 个六面体单元的模型，步长减半，时间差不多翻倍。所以子程序内部要做优化：常数项提前算好，只对真正发生变形的单元做 EXP 运算，未变形区直接跳过。判断方式可以简单用应变增量阈值，也可以用单元是否在模具接触区内。我一般两个都用，先按接触筛一遍，再看应变。

4. 摇臂轴镦锻算例：52111 单元模型里的四个关键设置

4.1 对称面简化与网格规模分配

算例对象是扇形齿轮摇臂轴，一步镦锻成形，材料 20CrMnTi，模具材料 5CrNiMo。工件为平面对称，取一半建模，坯料划分 52111 个六面体单元、7085 个节点，上模 1918 个单元，下模 5823 个单元。

这里有个必须点出来的问题：52111 个六面体单元配 7085 个节点，比例明显不对称。按六面体网格的常识，单元数和节点数应该在同一量级，一比七的比例在正常划分下不可能出现。我判断资料里的数字是分开统计的（可能是坯料单元数与某一部分节点数）。抄参数的时候，网格规模按自己前处理里的实际统计来，别硬套这两个数字，否则你在前处理里怎么调都对不上。

模具按刚体处理，不做网格重划分，但模具的网格密度要足够描述型腔圆角，否则接触判断会在圆角处抖动。上模 1918、下模 5823 这个差异，通常是因为下模型腔复杂、需要更密的表面网格，这个思路是对的：坯料网格跟着变形走，模具网格跟着几何走。

对称面约束按面施加，法向位移固定。镦锻这类工艺对称面选错是常见失误，摇臂轴是非轴对称件，只有一个平面对称，别按轴对称去简化，否则整个应变场都是错的。

4.2 材料与模具的物性录入

20CrMnTi 的流变曲线要用自己的热压缩数据建库，温度点覆盖 850 到 1200 ℃，应变率覆盖 0.01 到 10 /s。刚粘塑性模型里流变应力是应变、应变率、温度的三维曲面，插值方式选对数应变率插值比线性稳，因为应变率跨两个数量级。

热物性参数同样重要：导热系数、比热、密度随温度变化。这些参数对温度场的影响是直接的，而温度场又通过 EXP(Q/RT) 直接放大到晶粒尺寸上。我见过有人流变曲线做得很细，热物性用常温常数，结果模拟出来的芯部温度和实测差 60 K，晶粒尺寸自然对不上。

模具材料 5CrNiMo 按刚性体处理，但热传导要参与，模具的导热系数和界面换热系数不能省。界面换热系数在热锻里通常取 5 到 20 kW/(m²·K) 量级，取值不同会影响坯料表层温度，进而影响表层晶粒。这个参数没有理论值可算，只能靠热电偶实测比对反调。

4.3 步长、模具速度与摩擦边界

模具速度按实际压机曲线给，不要用恒定速度。液压机在接触瞬间有速度爬升，这个阶段的应变率变化很剧烈，直接影响动态再结晶的启动时刻。如果只有恒定速度数据，至少在接触后的头几个步用小步长过渡。

摩擦模型用剪切摩擦，摩擦因子取 0.3 到 0.6。镦锻这种高压力工况，库仑摩擦会在高正应力区失效，剪切摩擦更稳。摩擦因子对等效应变分布的影响很大，边角处的大变形区范围会跟着变，而大变形区正是动态再结晶最充分的区域，所以摩擦因子实际上是间接影响晶粒分布的一个参数。

步长设置上，变形阶段按位移控制，每步压下量约为最小单元尺寸的 1/10；进入保压和冷却阶段后改用时间步长，同时通过子程序回写 DTMAXC 限制最大步长，防止晶粒长大项过冲。

4.4 结果判读：四张云图怎么对读

模拟输出四类场量：温度场、等效应变、等效应变率、再结晶分数与晶粒尺寸。单看任何一张都容易误判，必须交叉对读。

温度场先看分布形态：镦锻时坯料表层接触冷模，温度低于芯部，形成"外冷内热"。如果模拟出来芯部反而比表层冷，先查界面换热系数方向是不是反了。

等效应变和应变率要一起看：应变大不代表再结晶充分，还得看应变率。同样的应变量，高应变率下 eps50 更小，再结晶反而更快，这是 Avrami 模型里应变率指数项的直接结果。摇臂轴这类件，法兰根部应变最大、应变率也高，通常是最先完成再结晶的区域。

再结晶分数云图看的是"哪些区域完成了"。资料里提到在这种等效应变分布下坯料内各区域都发生了不同程度的再结晶，重点在"不同程度"——如果云图上是清一色的 0.99，说明模型参数偏软或者步长太大，分数被一步顶到上限了，这时的晶粒尺寸不可信。

晶粒尺寸云图最后看，并且要和应变云图对照：细化区应该和难变形区、大变形区的边界吻合。如果晶粒尺寸的分布和应变分布完全无关，八成是状态变量索引错了，或者网格重划分把晶粒尺寸重置成了初始值。这种"看起来很正常但物理上讲不通"的结果，是最需要警惕的。

5. 避坑与排查：五条把结果算飞的记录

5.1 编译通过、结果全场恒定：子程序根本没被调用

现象：算例能跑完，温度、应变云图都正常，唯独晶粒尺寸全场一个值，等于初始晶粒尺寸，再结晶分数恒为零。

原因：最常见的是运行时没有选自定义求解器，跑的还是发行版。其次是 def_usr.f 里的 USRMSH 段被注释掉了，或者改的是另一个同名文件。还有一种隐蔽情况：子程序被调用了，但循环上界写成了 NUMNP（节点数）而不是 NUMEL（单元数），数组越界读到的是相邻单元的零值。

解决：在子程序入口写一行调试输出，记录调用次数和第一个单元的输入值，跑 5 步看日志。确认调用正常后再查循环上界。同时在前处理的求解器设置里明确指定用户子程序版本，别依赖默认。

5.2 再结晶分数超过 1 或变成负数

现象：后处理里 VF 出现 1.03 或者 -0.02 这种值，晶粒尺寸随之出现不合理的小值或负值。

原因：增量式更新没有做上下限截断，Avrami 曲线在 eps 略小于 eps50 时增长很平缓，一旦 eps 跨过 eps50 就陡增，大步长下很容易一步跨过 1。另一个来源是时间步被求解器自动放大，而你只在小步长下验证过。

解决：所有分数类变量写完都做一次 0 到 0.995 的截断，晶粒尺寸做下限保护。同时在子程序里回写 DTMAXC，把最大步长卡住。别指望事后在云图上修，负的晶粒尺寸会让长大项的幂运算直接产生 NaN，一旦出现会污染后续所有步。

5.3 应变率量级离谱，晶粒尺寸整体偏小一个数量级

现象：等效应变率云图出现 1e7、1e8 这样的量级，晶粒尺寸全场偏小，细化得不像话。

原因：单位制不统一。几何用 mm、时间用 s、力用 N 的 mm-t-s 制，和 m-kg-s 制混用，应变率差 1000 倍。更隐蔽的是热激活项：QACT 用 kJ/mol 而 RGAS 用 8.314 J/(mol·K)，指数项差 1000 倍，EXP 直接溢出或者归零。

解决：开工前先把单位制写进项目备注，全套统一。QACT 和 RGAS 必须同单位，要么都用 J 配 8.314，要么都用 kJ 配 0.008314。写完代码先拿一个已知解析解的等温均匀压缩算例验一遍，量级对了再上真实模型。

5.4 一改网格，晶粒尺寸云图就跳变

现象：同一套工艺参数，网格加密 20% 后，晶粒尺寸分布形态明显改变，局部出现块状跳变。

原因：网格重划分后，用户状态变量和晶粒尺寸没有做插值映射，新单元被初始化为前处理里设定的默认值。重划分越频繁的区域，晶粒尺寸被"重置"的次数越多，表现出来就是网格越密、结果越乱。

解决：在重划分触发的那一步做变量映射，把旧单元的状态量按体积加权插到新单元上。常见做法是利用求解器提供的重划分回调，或者自己在子程序里判断单元编号变化后重新赋值。这一步没有后悔药，只能在建模阶段就把重划分准则调稳，尽量减少重划分次数。

5.5 静态再结晶在变形过程中就启动了

现象：变形还没结束，DFN 已经有值，变形阶段的晶粒尺寸曲线出现二次下降。

原因：阶段判据写成了纯温度判据（温度高于某值就开始算静态再结晶），而变形过程中温度本来就高，判据恒成立。另一个来源是计时基准写错，USRVE 里的时间戳每步都刷新，DT 恒为零，理论上不该启动，但如果用了未初始化的变量，行为就是随机的。

解决：阶段判据必须包含"变形是否停止"这个信息，用模具速度归零或者接触状态解除作为触发条件，并用一个阶段标志变量锁存，一旦进入阶段 2 就不再回退。计时基准只在进入阶段 2 的那一步写一次，之后只读不写。这类逻辑错误不会报错，只能靠单单元算例把每个分支都走一遍来验。

6. 拿什么证明晶粒尺寸不是玄学：三条验证路径与一个习惯

晶粒尺寸这种算出来没有解析解、测出来又有分布的结果，最怕的是"看着像那么回事"就收工。我一般走三条验证路径，从软到硬逐级加码。

第一条最便宜，也最能定位问题。把子程序里的那套公式单独抽出来，用热压缩试验的 T、ε、ε̇ 路径离线算一遍，和实测晶粒尺寸对比。如果这条都过不了，后面的成形模拟全是白算。

第二条是真正的落地验证。剖切位置要选在变形历史差异大的地方，法兰根部变形大、腹板中心接近不变形、飞边应变率极高但温降快，三个位置的晶粒尺寸应该拉开明显差距。如果模拟出来三个位置几乎一样，说明模型对场量不敏感，参数多半是拍脑袋的。

第三条用来校验分布形态。硬度与晶粒尺寸是 Hall-Petch 关系，虽然硬度还受相组成影响，但单调性一般成立。沿轴线打一排硬度点，和模拟出来的晶粒尺寸沿同一条线对比，形态对得上就基本可信。

从那以后我每次上新钢种，都强制先走一遍第一级：离线复算热压缩路径，偏差超过 15% 就不往下走，绝不用"成形模拟出来再看"来赌。这套流程和子程序骨架我已经整理成可复用的模板，配合资料里的形参表和算例参数，改材料系数就能跑起来，需要的按平时取资料的路径自取。希望帮到你。

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