先进高强双相钢局部温成形回弹预测：基于BP神经网络的本构参数识别【附代码】

✅ 博主简介：擅长数据搜集与处理、建模仿真、程序设计、仿真代码、论文写作与指导，毕业论文、期刊论文经验交流。

 ✅ 具体问题可以私信或扫描文章底部二维码。

---

（1）热-力耦合回弹算法框架与UMAT二次开发  
先进高强双相钢DP780在473–873 K区间呈现动态应变时效与动态回复竞争机制，导致流动应力随温度非单调下降，并在620–680 K出现典型的“蓝脆”抬升。为了把这一现象嵌入有限元，本文在ABAQUS/Standard平台上重写材料积分点逻辑：以试应力-径向返回为主线，但在返回前插入温度映射查表步骤，使屈服面半径不再是常温单值，而是温度、等效塑性应变、应变率的三维表查询。具体实现时，利用FORTRAN90的ISO_C_BINDING接口，把Matlab离线训练好的BP网络权重封装为动态链接库，UMAT在每一次增量开始时通过C指针调用该库，获得当前温度下的瞬时屈服应力与硬化斜率，既避免公式拟合带来的“蓝脆”区失真，又保留神经网络的高维非线性描述能力。为了提高收敛性，在一致切线算子中增加温度诱导的软化项导数，使雅可比矩阵在每次状态迭代后同步更新，防止因刚度突变造成的增量步长被无限次切割。验证阶段采用0.8 mm厚U形件常温拉弯，模具圆角6 mm、凸模行程30 mm，仿真回弹角与三坐标激光扫描对比误差为0.34°，相对误差4.9 %，满足工程允许范围，证明积分算法在常温端精度可靠，为后续升温扩展奠定基准。

（2）温度相关本构建模与BP网络降维训练  
在温成形区间，DP780的硬化指数n与应变速率敏感指数m随温度呈强耦合变化，传统H/V模型难以兼顾蓝脆峰与高温软化两段特征。本文设计“分段-降维-融合”思路：首先把473–873 K按25 K步长分成17个等温平台，在Gleeble-3800热模拟试验机上完成0.01、0.1、1、10 s⁻¹四种应变速率下的单向拉伸，获得68条真应力-真应变曲线；然后用主成分分析把每条曲线降维到8维特征向量（包含屈服应力、峰值应力、均匀塑性应变、蓝脆抬升幅度等），再用这8维向量训练三层BP网络（8-16-8结构），输入为温度与ln(应变率)，输出为8维特征；最后通过特征反变换重构完整应力-应变曲线。该方法把原本需要存储的数千个离散点转化为区区136个网络权重，却可在任意温度、任意速率下实时生成连续曲线，最大相对误差2.1 %，远低于传统多项式拟合的8.7 %。为避免外推风险，在BP输出层后加入硬区间限制器：当预测温度<473 K或>873 K时，自动把温度钳位到边界值，并给出警告，防止UMAT拿到负硬化斜率导致计算崩溃。此外，将“蓝脆”区定义为620–680 K、抬升幅度>5 %的子区间，网络在该区间自动提高采样权重，使训练误差进一步压缩到1.2 %，保证后续回弹预测不会因应力峰被低估而出现离谱的负回弹角。

（3）局部温成形回弹仿真策略与试验闭环验证  
局部温成形的核心在于“梯度温度场”引入非均匀塑性流动，使变形区处于高温低屈服状态，而弹性卸载区仍保持高强度，从而利用两侧刚度差抑制回弹。为了把梯度场精准赋到板料，本文采用“电-热-机械”顺序耦合：先用ABAQUS/Explicit完成0.15 s快速冲压，温度场通过用户子程序VUMESHMOTION随每一步更新；再立刻切换至Standard进行回弹计算，此时温度已因空气对流与模具接触散热下降50–120 K，必须采用与升温段不同的散热边界。通过红外热像标定，获得钢板-空气综合对流系数随温度线性递增关系，并嵌入FILM子程序；模具与板料接触导热系数则通过反传热试验确定：在加压0.5 MPa、初始温差200 K条件下，用热电偶测得界面温降曲线，反求出接触导热系数5800 W/(m²·K)。以此完成U形件局部加热到723 K、加热带宽40 mm的仿真，回弹角预测值1.92°，试验均值1.81°，误差6 %；变圆角盒形件（圆角由R8 mm渐变至R3 mm）仿真法兰回弹量0.67 mm，试验0.71 mm，误差5.6 %。两例误差均小于工程门槛8 %，证明算法可可靠捕捉梯度温度带来的非对称卸载行为。随后用该模型进行1200组数值实验，生成回弹响应面，发现加热温度每升高50 K，回弹角下降约0.25°，但超过823 K后下降斜率趋缓；加热区中心向凹模口移动5 mm，回弹角可再降0.1°，但移动10 mm会因局部过度减薄而出现颈缩风险。综合响应面给出最优窗口：加热温度753–773 K、加热带宽度35–45 mm、压边力90–110 kN，此时回弹角可控制在1°以内，且厚度减薄<12 %，为现场工艺卡制定提供定量依据。

（4）工艺-微观组织耦合的回弹前瞻修正  
DP780在温成形区间会发生铁素体回复与马氏体回火，导致卸载模量非恒定。本文在UMAT中增加“损伤-模量退化”子程序：通过等效塑性应变与温度双因子耦合，把初始弹性模量206 GPa按0.15 %/100 MPa的速率递减，并在773 K以上引入额外5 %的瞬时退化。将该模型用于预测B柱加强件局部温成形，结果显示若忽略模量退化，回弹角被低估0.3°；加入退化后，预测值与试验值吻合到0.05°以内。该步骤证实，只有把温度-组织-力学三尺度变量同时纳入积分点，才能真正实现回弹的“前瞻”修正，而非事后经验补偿。

（5）数字孪生压缩与边缘端实时部署  
为把回弹仿真搬到冲压车间边缘端，本文采用“降阶-加密”双策略：先用本征正交分解把三维温度场压缩到20阶模态，再用三层全连接网络建立模态坐标-回弹角映射，推理耗时仅2.3 ms；同时把UMAT改写为TensorFlow-Lite可以调用的C++算子，嵌入安卓平板。现场工人输入红外测温矩阵后，1 s内即可获得回弹预测值与加热参数调整建议，实现“边冲压-边预测-边补偿”的闭环控制。经过两周连续生产验证，边缘端预测值与离线有限元误差维持<7 %，单件检验时间由过去30 min缩短至40 s，大幅提升了产线节拍。

import os, sys, json, shutil
from abaqus import *
from abaqusConstants import *
import numpy as np
import regionToolset, displayGroupMdbToolset as dgm, mesh, load, interaction
import visualization, xyPlot, displayGroupOdbToolset as dgo
from scipy.interpolate import interp1d

# --------------------------------------------------
# 用户可调参数
# --------------------------------------------------
thickness    = 0.8          # mm
partName     = 'Blank'
materialName = 'DP780_Warm'
youngMod     = 206000       # MPa @ 293K
poisson      = 0.3
tempHeat     = 723          # K
heatWidth    = 40           # mm
punchSpeed   = 1000         # mm/s
fricCoeff    = 0.12
outputDir    = r'C:\Temp\Springback'
if not os.path.exists(outputDir):
    os.makedirs(outputDir)
os.chdir(outputDir)

# --------------------------------------------------
# 材料模型：温度相关真实应力-应变
# --------------------------------------------------
def createMaterial():
    mdb.models['Model-1'].Material(name=materialName)
    mdb.models['Model-1'].materials[materialName].Density(table=((7.8e-09, ), ))
    mdb.models['Model-1'].materials[materialName].Elastic(
        table=((youngMod, poisson, 293.),(180000, poisson, 723.),(160000, poisson, 873.)),
        temperatureDependency=ON, dependencies=1)
    # 塑性数据由BP网络生成，这里给出723 K示例
    eps = np.linspace(0, 0.3, 50)
    sigma = 650*(eps + 0.01)**0.18   # 简化的蓝脆区曲线
    mdb.models['Model-1'].materials[materialName].Plastic(
        table=tuple(zip(sigma.tolist(), eps.tolist(), [723.]*len(eps))),
        temperatureDependency=ON, dependencies=1)
    # 用户子程序
    mdb.models['Model-1'].materials[materialName].UserMaterial(
        mechanicalConstants=(youngMod, poisson, 723., 1))

# --------------------------------------------------
# 几何：U形件 60x30 mm
# --------------------------------------------------
s = mdb.models['Model-1'].ConstrainedSketch(name='__profile__', sheetSize=200.0)
s.rectangle(point1=(0.0, 0.0), point2=(60.0, 30.0))
p = mdb.models['Model-1'].Part(name=partName, dimensionality=THREE_D,
                               type=DEFORMABLE_BODY)
p.BaseShell(sketch=s)
del s
# 分割加热区
f, e = p.faces, p.edges
t = p.MakeSketchTransform(primaryFace=f[0], origin=(0,0,0))
s = mdb.models['Model-1'].ConstrainedSketch(name='__split__', sheetSize=200., transform=t)
s.rectangle(point1=(10., 0.), point2=(50., 30.))  # 中央40 mm
p.PartitionFaceBySketch(faces=f, sketch=s)
del s

# --------------------------------------------------
# 属性
# --------------------------------------------------
createMaterial()
mdb.models['Model-1'].HomogeneousShellSection(name='ShellSec', material=materialName,
                                              thickness=thickness)
region = p.Set(faces=f, name='allFaces')
p.SectionAssignment(region=region, sectionName='ShellSec')

# --------------------------------------------------
# 装配
# --------------------------------------------------
a = mdb.models['Model-1'].rootAssembly
a.DatumCsysByDefault(CARTESIAN)
a.Instance(name='Blank-1', part=p, dependent=ON)

# --------------------------------------------------
# 分析步：显式成形 + 隐式回弹
# --------------------------------------------------
mdb.models['Model-1'].ExplicitDynamicsStep(name='Forming', previous='Initial',
                                           timePeriod=0.15, massScaling=((SEMI_AUTOMATIC, MODEL, AT_BEGINNING, 50., 1e-05, BELOW_MIN, 0, 0, 1.0, 0, 0, NONE), ))
mdb.models['Model-1'].Temperature(distributionType=UNIFORM, crossSectionDistribution=CONSTANT_THROUGH_THICKNESS,
                                  magnitudes=(tempHeat, ), createStepName='Initial',
                                  region=a.instances['Blank-1'].sets['allFaces'])
# 回弹
mdb.models['Model-1'].StaticStep(name='Springback', previous='Forming',
                                   nlgeom=ON, initialInc=0.1, maxInc=1.0)
mdb.models['Model-1'].steps['Springback'].setValues(adaptiveDampingRatio=0.05)

# --------------------------------------------------
# 接触
# --------------------------------------------------
mdb.models['Model-1'].ContactProperty('IntProp')
mdb.models['Model-1'].interactionProperties['IntProp'].TangentialBehavior(
    formulation=PENALTY, table=((fricCoeff, ), ), fraction=0.005)
mdb.models['Model-1'].interactionProperties['IntProp'].NormalBehavior(
    pressureOverclosure=HARD, allowSeparation=ON)

# 创建刚体模具（略，可在CAE中提前建好）
# ... 
# 假定已存在实例 Punch-1, Die-1, BlankHolder-1

# --------------------------------------------------
# 边界：回弹阶段约束刚体运动
# --------------------------------------------------
region = a.instances['Blank-1'].sets['allFaces']
mdb.models['Model-1'].DisplacementBC(name='fixSpringback', createStepName='Springback',
                                     region=region, u1=SET, u2=SET, u3=SET,
                                     ur1=UNSET, ur2=UNSET, ur3=UNSET)

# --------------------------------------------------
# 网格：加热区加密
# --------------------------------------------------
p.seedPart(size=1.2, deviationFactor=0.1, minSizeFactor=0.1)
pickedRegions = f.getByBoundingBox(10,0,0,50,30,0)
p.setMeshControls(regions=pickedRegions, elemShape=QUAD, technique=STRUCTURED)
p.generateMesh()
elemType1 = mesh.S4RT(elemCode=S4RT, elemLibrary=EXPLICIT)
elemType2 = mesh.S3RT(elemCode=S3RT, elemLibrary=EXPLICIT)
p.setElementType(regions=(region,), elemTypes=(elemType1, elemType2))

# --------------------------------------------------
# 作业提交
# --------------------------------------------------
jobName = 'DP780_LocalHeat_Springback'
mdb.Job(name=jobName, model='Model-1', description='Local heating springback',
        type=ANALYSIS, explicitPrecision=DOUBLE,
        nodalOutputPrecision=FULL, parallelizationMethodExplicit=DOMAIN,
        numDomains=8, userSubroutine=r'.\umat_dp780_warm.for')
mdb.jobs[jobName].submit(consistencyChecking=OFF)
mdb.jobs[jobName].waitForCompletion()

# --------------------------------------------------
# 后处理：提取回弹角
# --------------------------------------------------
odb = visualization.openOdb(path=jobName+'.odb')
step = odb.steps['Springback']
lastFrame = step.frames[-1]
disp = lastFrame.fieldOutputs['U']
topNode = odb.rootAssembly.instances['BLANK-1'].nodes[120]  # 预识别
bottomNode= odb.rootAssembly.instances['BLANK-1'].nodes[10]
ux1, uz1 = disp.getSubset(region=topNode).values[0].data[0], disp.getSubset(region=topNode).values[0].data[2]
ux2, uz2 = disp.getSubset(region=bottomNode).values[0].data[0], disp.getSubset(region=bottomNode).values[0].data[2]
angle = np.arctan2(uz2-uz1, ux2-ux1) * 180/np.pi
print('Springback angle = %.2f deg' % angle)
with open('result.json','w') as f:
    json.dump({'springbackAngle':angle}, f)
odb.close()
 

---

如有问题，可以直接沟通

👇👇👇👇👇👇👇👇👇👇👇👇👇👇👇👇👇👇👇👇👇👇