CN102411648A - 高强度钢板温热成形数值模拟方法 - Google Patents

Abstract

本发明公开了一高强度钢板温热成形数值模拟方法，包括：根据用户输入的材料参数或模型调用指令，建立待模拟模具的分析模型；接收用户输入的工艺参数信息，并根据用户输入的工艺参数信息计算得到待模拟模具的热、力及相变耦合模型，动力显示模型，接触及摩擦模型，和/或热成形温度场模型后显示。本发明提供的高强度钢板温热成形数值模拟方法实现了对温热成形过程中，待模拟材料的数值模拟功能，可以对温热成形模具设计及工艺参数的选取进行自动分析，进而提供优化方案，相较于现有的人工选取方式，可靠性及精确度更高，节省了制程成本。

Description

高强度钢板温热成形数值模拟方法

技术领域

本发明属于热处理技术领域，尤其涉及一高强度钢板温热成形数值模拟方法。

背景技术

当前，安全、环保、节能成为汽车制造业发展的主题，采用高强度钢板冲压件制造车身既可以减轻车体重量，又能提高安全性，是同时实现车体轻量化和提高碰撞安全性的最好途径。高强度钢板，强度越高，越难成形。尤其是当强度达到1500MPa时，常规的冷冲压成形工艺几乎无法成形。因此，如何实现高强度钢板的高精度冲压成形就成为一项迫切需要解决的技术难题。

温热成形技术(即：热冲压技术)，是一种专门用于成形高强度冲压件的先进制造技术，可用来成形强度高达1600MPa的冲压件，目前成为世界上众多汽车生产厂商关注的热点。现有的温热成形技术依靠人的经验来指导温热成形模具设计及工艺参数的选取，可靠性差且浪费人力。

发明内容

本发明的目的在于提供一高强度钢板温热成形数值模拟方法，以解决现有的温热成形技术依靠人的经验来指导温热成形模具设计及工艺参数的选取，可靠性差且浪费人力的问题。

本发明实施例是这样实现的，一种高强度钢板温热成形数值模拟方法，所述方法包括以下步骤：

根据用户输入的材料参数或模型调用指令，建立待模拟模具的分析模型；

接收用户输入的工艺参数信息，并根据用户输入的所述工艺参数信息计算得到待模拟模具的热、力及相变耦合模型，动力显示模型，接触及摩擦模型，和/或热成形温度场模型后显示。

所述根据用户输入的所述工艺参数信息计算得到待模拟模具的热、力及相变耦合模型的步骤可以表示为：

$\left\{\begin{array}{c}\mathrm{\ξ}=1-\exp [-\mathrm{\θ}(M_S-T)]\\ {\stackrel{\·}{\mathrm{\ϵ}}}^{\mathrm{tp}}=\mathrm{KS}(1-\mathrm{\ξ})\stackrel{\·}{\mathrm{\ξ}}\\ {\stackrel{\·}{\mathrm{\ϵ}}}^{\mathrm{tr}}=\mathrm{\β}\stackrel{\·}{\mathrm{\ξ}}\end{array}\right.$

其中，ξ为马氏体转变量；

为相变塑性应变；S为偏应力；

为相变体积应变增量；β为相变膨胀系数；K为相变塑性系数；T为温度；Ms为马氏体相变开始温度；θ为马氏体转变速率参数、k为相变塑性系数；

所述马氏体相变开始温度Ms、马氏体转变速率参数θ、相变塑性系数k满足：

$\left\{\begin{array}{cc}{M}_s=M_{s0}& \stackrel{\&OverBar;}{\mathrm{\&sigma;}}\&le;{\mathrm{\&sigma;}}_{a0}\\ M_s=M_{s0}+\mathrm{\&chi;}\&times;(\stackrel{\&OverBar;}{\mathrm{\&sigma;}}-{\mathrm{\&sigma;}}_{a0})& {\mathrm{\&sigma;}}_{a0}<\stackrel{\&OverBar;}{\mathrm{\&sigma;}}<{\mathrm{\&sigma;}}_{p0}\\ M_s=M_{s0}+\mathrm{\&chi;}\&times;({\mathrm{\&sigma;}}_{p0}-{\mathrm{\&sigma;}}_{a0})& \stackrel{\&OverBar;}{\mathrm{\&sigma;}}\&GreaterEqual;{\mathrm{\&sigma;}}_{p0}\end{array}\right.$

$\left\{\begin{array}{cc}\mathrm{\&theta;}={\mathrm{\&theta;}}_0& \stackrel{\&OverBar;}{\mathrm{\&sigma;}}\&le;{\mathrm{\&sigma;}}_a\\ \mathrm{\&theta;}={\mathrm{\&theta;}}_0+\mathrm{\&lambda;}\&times;(\stackrel{\&OverBar;}{\mathrm{\&sigma;}}-{\mathrm{\&sigma;}}_a)& {\mathrm{\&sigma;}}_a<\stackrel{\&OverBar;}{\mathrm{\&sigma;}}<{\mathrm{\&sigma;}}_p\\ \mathrm{\&theta;}={\mathrm{\&theta;}}_0+\mathrm{\&lambda;}\&times;({\mathrm{\&sigma;}}_p-{\mathrm{\&sigma;}}_a)& \stackrel{\&OverBar;}{\mathrm{\&sigma;}}\&GreaterEqual;{\mathrm{\&sigma;}}_p\end{array}\right.$

$\left\{\begin{array}{cc}K=K_0& \stackrel{\&OverBar;}{\mathrm{\&sigma;}}\&le;{\mathrm{\&sigma;}}_a^{\&prime;}\\ K=K_0+\mathrm{\&gamma;}\&times;(\stackrel{\&OverBar;}{\mathrm{\&sigma;}}-{\mathrm{\&sigma;}}_a^{\&prime;})& {\mathrm{\&sigma;}}_a^{\&prime;}<\stackrel{\&OverBar;}{\mathrm{\&sigma;}}<{\mathrm{\&sigma;}}_p^{\&prime;}\\ K=K_0+\mathrm{\&gamma;}\&times;({\mathrm{\&sigma;}}_p^{\&prime;}-{\mathrm{\&sigma;}}_a^{\&prime;})& \stackrel{\&OverBar;}{\mathrm{\&sigma;}}\&GreaterEqual;{\mathrm{\&sigma;}}_p^{\&prime;}\end{array}\right.$

其中，

为应力，M_s0为马氏体点基点温度，χ为相变温度与应力之间的耦合系数，σ_a0、σ_p0分别为马氏体相变点与应力耦合关系的区间应力两基点；θ₀为马氏体转变速率基点，λ为相变速率与应力之间的耦合系数，σ_a、σ_p分别为马氏体相变速率与应力耦合关系的区间应力两基点；K₀为相变塑性系数基点，γ为相变塑性系数与应力之间的耦合系数。

所述动力显示模型为应力的导数张量模型，所述根据用户输入的所述工艺参数信息计算得到动力显示模型的步骤可以表示为：

${\stackrel{\&dtri;}{\mathrm{\&sigma;}}}_{\mathrm{ij}}={\stackrel{\&OverBar;}{D}}_{\mathrm{ijkl}}^{\mathrm{ep}}({\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_{\mathrm{kl}}\left(t\right)-{\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_{\mathrm{kl}}^{\mathrm{th}}\left(t\right)-{\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_{\mathrm{kl}}^{\mathrm{tr}}\left(t\right)-{\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_{\mathrm{kl}}^{\mathrm{tp}}\left(t\right))+{\stackrel{\&CenterDot;}{\stackrel{\&OverBar;}{D}}}_{\mathrm{ijkl}}^{\mathrm{ep}}({\mathrm{\&epsiv;}}_{\mathrm{kl}}\left(t\right)-{\mathrm{\&epsiv;}}_{\mathrm{kl}}^{\mathrm{th}}\left(t\right)-{\mathrm{\&epsiv;}}_{\mathrm{kl}}^{\mathrm{tr}}\left(t\right)-{\mathrm{\&epsiv;}}_{\mathrm{kl}}^{\mathrm{tp}}\left(t\right))$

其中，

是弹塑性本构矩阵，

是以相邻构型为参考构型的Green应变率张量，是以相邻构型为参考构型的热应变率张量，

是以相邻构型为参考构型的相变体积应变率张量，

是以相邻构型为参考构型的相变塑性应变率张量。

所述动力显示模型为段流动应力模型，所述根据用户输入的所述工艺参数信息计算得到动力显示模型的步骤可以表示为：

$\mathrm{\&sigma;}(\mathrm{\&epsiv;},\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}},T)=K^{\&prime;}{(b+\mathrm{\&epsiv;})}^{n\left(T\right)}{\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}^{m\left(T\right)}\exp ({\mathrm{\&beta;}}^{\&prime;}/T)$

其中，K′为强度系数，b为应变修正量，β′为温度系数，n(T)与m(T)满足：

n(T)＝n₀exp(-c_n(T_i-T₀))

m(T)＝m₀exp(-c_m(T_i-T₀))

其中，n₀、m₀、c_n、c_m分别为硬化材料常数，其值通过上述拉伸实验确定；T₀为温度初始值，T_i为温度当前值。

当所述接触及摩擦模型为动力平衡运动模型时，所述根据用户输入的所述工艺参数信息计算得到接触及摩擦模型的步骤可以表示为：

$\frac{M}{\mathrm{\&Delta;}{t}_n^2}[u_{n+1}-(1+\mathrm{\&beta;}1)u_n+\mathrm{\&beta;}1u_{n-1}]+\frac{C}{\mathrm{\&Delta;}{t}_n}(u_{n+1}-u_n)+f_{n+1}^{\mathrm{int}}=f_{n+1}^{\mathrm{ext}}+{\mathrm{\&lambda;}}_{n+1}$

其中，Δt_n＝t_n+1-t_n，β1＝Δt_n/Δt_n-1，M为质量矩阵，C为阻尼矩阵，且满足C＝αM，α为阻尼系数；λ_n+1为t_n+1时刻摩擦接触反力，

为t_n+1时刻节点内力，u_n+1为t_n+1时刻节点位移向量，

为t_n+1时刻节点外力，u_n为t_n时刻节点位移向量，u_n-1为t_n-1时刻节点位移向量。

当所述接触及摩擦模型为法向接触反力和切向摩擦力模型时，所述根据用户输入的所述工艺参数信息计算得到接触及摩擦模型的步骤可以表示为：

其中，

为接触力控制系数，且满足

其中，

为高温下待模拟材料的屈服强度，

为室温下待模拟材料的屈服强度，m为待模拟材料上一点的质量，d_N为所述点的法向穿透距离，d_T为所述点沿模具表面的切向滑移距离。

所述根据用户输入的所述工艺参数信息计算得到热成形温度场模型的步骤可以表示为：

$(\frac{11}{6\mathrm{\&Delta;t}}{C}^e+K^e)T_{t+\mathrm{\&Delta;t}}=F^e+\frac{1}{\mathrm{\&Delta;t}}{C}^e(3T_t-1.5T_{t-\mathrm{\&Delta;t}}+\frac{1}{3}{T}_{t-2\mathrm{\&Delta;t}})$

其中，K^e为热传导系数，C^e为阻尼矩阵，F^e为温度载荷矩阵，T_t+Δt为t+Δt时刻的温度值，T_t为t时刻的温度值，T_t-Δt为t-Δt时刻的温度值，T_t-2Δt为t-2Δt时刻的温度值。

在所述根据用户输入的所述工艺参数信息计算得到待模拟模具的动力显示模型的步骤之后，还可以包括对所述动力显示模型进行离散化并显示的步骤，所述对动力显示模型进行离散化的步骤表示为：

${\mathrm{\&sigma;}}_{\mathrm{ij}}^{n+1}={\widehat{\mathrm{\&sigma;}}}_{\mathrm{ij}}^{n+1}+{\left({\stackrel{\&OverBar;}{D}}_{\mathrm{ijkl}}^{\mathrm{ep}}\right)}^n({\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_{\mathrm{kl}}^{n+\frac{1}{2}}\left(t\right)-{\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_{\mathrm{th}}^{{\mathrm{th}}^{n+\frac{1}{2}}}\left(t\right)-{\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_{\mathrm{kl}}^{{\mathrm{tr}}^{n+\frac{1}{2}}}\left(t\right)-{\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_{\mathrm{kl}}^{{\mathrm{tp}}^{n+\frac{1}{2}}}\left(t\right))\mathrm{\&Delta;}{t}_n+{\left({\stackrel{\&CenterDot;}{\stackrel{\&OverBar;}{D}}}_{\mathrm{ijkl}}^{\mathrm{ep}}\right)}^n({\mathrm{\&epsiv;}}_{\mathrm{kl}}^{n+\frac{1}{2}}\left(t\right)-{\mathrm{\&epsiv;}}_{\mathrm{kl}}^{{\mathrm{th}}^{n+\frac{1}{2}}}\left(t\right)-{\mathrm{\&epsiv;}}_{\mathrm{kl}}^{{\mathrm{tr}}^{n+\frac{1}{2}}}\left(t\right)-{\mathrm{\&epsiv;}}_{\mathrm{kl}}^{{\mathrm{tp}}^{n+\frac{1}{2}}}\left(t\right))\mathrm{\&Delta;}{t}_n$

其中，Δt_n为时间增量，是弹塑性本构矩阵，

是以相邻构型为参考构型的Green应变率张量，

是以相邻构型为参考构型的热应变率张量，

是以相邻构型为参考构型的相变体积应变率张量，

是以相邻构型为参考构型的相变塑性应变率张量。

所述根据用户输入的所述工艺参数信息计算得到待模拟模具的动力显示模型的步骤之后，还可以包括计算下一时刻节点总位移矢量并显示的步骤，所述计算下一时刻节点总位移矢量的步骤表示为：

$u_{n+1}=u_n+{\stackrel{\&CenterDot;}{u}}_{n+\frac{1}{2}}\mathrm{\&Delta;}{t}_n$

其中，u_n+1为t_n+1时刻节点的位移向量，u_n为t_n时刻节点的位移向量，

为t_n+1时刻节点的速度向量，Δt_n为时间增量。

所述时间增量Δt_n可以满足：

$\mathrm{\&Delta;}{t}_n\&le;r{L}_n/c,L_n=\min (L_n^e,e-\mathrm{1,2},...\mathrm{ne})$

其中，常数r值介于0.8到0.9之间，优选为0.5，c为待模拟材料的内部声速，可由

求得，而E(T)为随温度变化的弹性模量，

为第n状态单元e的名义长度，ne为单元总数。

本发明提供的高强度钢板温热成形数值模拟方法实现了对温热成形过程中，待模拟材料的数值模拟功能，可以对温热成形模具设计及工艺参数的选取进行自动分析，进而提供优化方案，相较于现有的人工选取方式，可靠性及精确度更高，节省了制程成本。

附图说明

图1是本发明提供的高强度钢板温热成形数值模拟方法的流程图；

图2是本发明中，拉伸实验记录的板材横向宽度横向膨胀量与温度之间的膨胀量-温度曲线图；

图3是本发明中，拉伸实验中应力与马氏体相变开始温度之间耦合关系实验结果图；

图4是本发明中，拉伸实验中马氏体转变速率参数与应力耦合关系实验结果图；

图5是本发明中，拉伸实验中相变塑性系数与应力耦合关系实验结果图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

如图1所示，本发明提供的高强度钢板温热成形数值模拟方法包括以下步骤：

在步骤S101中，根据用户输入的材料参数或模型调用指令，建立待模拟模具的分析模型。

在步骤S102中，接收用户输入的工艺参数信息，并根据用户输入的工艺参数信息计算得到待模拟模具的热、力及相变耦合模型，动力显示模型，接触及摩擦模型，和/或热成形温度场模型后显示。

对于热、力及相变耦合模型，本发明以目前应用较广泛的高强度钢板热成形材料22MnB5的硼合金钢为例，对其进行拉伸实验而得出热、力及相变耦合模型。该22MnB5的硼合金钢的材料成分质量百分含量如下表一所示：

表一

实验过程为：试样以10℃/s的加热速度升温到950℃，在950℃保温5分钟，之后以大于30℃/s的冷却速率降温至室温，在此降温过程中对试样进行位移加载，位移加载的温度范围为750℃到500℃并一直保持到室温，记录实验过程中力、位移、宽度方向膨胀量及温度随时间变化的曲线，实验结束后对试样中部温区线切割后进行金相分析，确定试样马氏体转化是否充分。

如图2所示，以实验中的某一次位移加载记录的板材横向宽度横向膨胀量与温度之间的膨胀量-温度曲线为例，说明实验过程中物理量的关系及马氏体相变开始温度Ms的确定方法。其中AB段的横向收缩量是由于降温过程中的热膨胀冷缩应变ε^th与拉应力引起的应变ε^e与ε^p热力耦合作用引起，BC段的横向收缩量是由于相变引起的ε^tr、ε^tp与ε^th、ε^e、ε^p等热、力、相变耦合作用引起的，CD段的横向收缩量为相变结束后ε^th、ε^e、ε^p热力耦合作用引起。本发明中，热、力、相变耦合参数θ、k由BC段得相变过程决定。

如图3示出了应力与马氏体相变开始温度Ms之间耦合关系实验结果。当应力

小于某一值时，相变开始温度Ms基本不变，当应力

在一定范围内变化时，相变开始温度Ms随应力

的增加而增加，当应力

大于某一值时，相变开始温度Ms趋于定值。

图4示出了马氏体转变速率参数θ与应力耦合关系实验结果，图5示出了相变塑性系数k与应力耦合关系实验结果。当应力小于某一值时，马氏体转变速率参数θ、相变塑性系数k基本不变，当应力

在一定范围内变化时，马氏体转变速率参数θ、相变塑性系数k随应力

的增加而增加，当应力

大于某一值时，马氏体转变速率参数θ、相变塑性系数k趋于定值。

对上述实验结果进行拟合，建立马氏体相变开始温度Ms、马氏体转变速率参数θ、相变塑性系数k与应力

之间的关系式为：

$\left\{\begin{array}{cc}{M}_s=M_{s0}& \stackrel{\&OverBar;}{\mathrm{\&sigma;}}\&le;{\mathrm{\&sigma;}}_{a0}\\ M_s=M_{s0}+\mathrm{\&chi;}\&times;(\stackrel{\&OverBar;}{\mathrm{\&sigma;}}-{\mathrm{\&sigma;}}_{a0})& {\mathrm{\&sigma;}}_{a0}<\stackrel{\&OverBar;}{\mathrm{\&sigma;}}<{\mathrm{\&sigma;}}_{p0}\\ M_s=M_{s0}+\mathrm{\&chi;}\&times;({\mathrm{\&sigma;}}_{p0}-{\mathrm{\&sigma;}}_{a0})& \stackrel{\&OverBar;}{\mathrm{\&sigma;}}\&GreaterEqual;{\mathrm{\&sigma;}}_{p0}\end{array}\right.---\left(1\right)$

$\left\{\begin{array}{cc}\mathrm{\&theta;}={\mathrm{\&theta;}}_0& \stackrel{\&OverBar;}{\mathrm{\&sigma;}}\&le;{\mathrm{\&sigma;}}_a\\ \mathrm{\&theta;}={\mathrm{\&theta;}}_0+\mathrm{\&lambda;}\&times;(\stackrel{\&OverBar;}{\mathrm{\&sigma;}}-{\mathrm{\&sigma;}}_a)& {\mathrm{\&sigma;}}_a<\stackrel{\&OverBar;}{\mathrm{\&sigma;}}<{\mathrm{\&sigma;}}_p\\ \mathrm{\&theta;}={\mathrm{\&theta;}}_0+\mathrm{\&lambda;}\&times;({\mathrm{\&sigma;}}_p-{\mathrm{\&sigma;}}_a)& \stackrel{\&OverBar;}{\mathrm{\&sigma;}}\&GreaterEqual;{\mathrm{\&sigma;}}_p\end{array}\right.---\left(2\right)$

$\left\{\begin{array}{cc}K=K_0& \stackrel{\&OverBar;}{\mathrm{\&sigma;}}\&le;{\mathrm{\&sigma;}}_a^{\&prime;}\\ K=K_0+\mathrm{\&gamma;}\&times;(\stackrel{\&OverBar;}{\mathrm{\&sigma;}}-{\mathrm{\&sigma;}}_a^{\&prime;})& {\mathrm{\&sigma;}}_a^{\&prime;}<\stackrel{\&OverBar;}{\mathrm{\&sigma;}}<{\mathrm{\&sigma;}}_p^{\&prime;}\\ K=K_0+\mathrm{\&gamma;}\&times;({\mathrm{\&sigma;}}_p^{\&prime;}-{\mathrm{\&sigma;}}_a^{\&prime;})& \stackrel{\&OverBar;}{\mathrm{\&sigma;}}\&GreaterEqual;{\mathrm{\&sigma;}}_p^{\&prime;}\end{array}\right.---\left(3\right)$

式(1)中，M_s0为马氏体点基点温度，χ为相变温度与应力之间的耦合系数，σ_a0、σ_p0分别为马氏体相变点与应力耦合关系的区间应力两基点，其取值与待模拟材料及冷却方式有关；式(2)中，θ₀为马氏体转变速率基点，λ为相变速率与应力之间的耦合系数，σ_a、σ_p分别为马氏体相变速率与应力耦合关系的区间应力两基点，其取值与待模拟材料及冷却方式有关；式(3)中，K₀为相变塑性系数基点，γ为相变塑性系数与应力之间的耦合系数，σ′_a、σ′_p分别为相变塑性系数与应力耦合关系的区间应力两基点，其取值与待模拟材料及冷却方式有关。

根据式(1)、式(2)及式(3)计算热、力及相变耦合本构模型，表示为：

$\left\{\begin{array}{c}\mathrm{\&xi;}=1-\exp [-\mathrm{\&theta;}(M_S-T)]\\ {\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}^{\mathrm{tp}}=\mathrm{KS}(1-\mathrm{\&xi;})\stackrel{\&CenterDot;}{\mathrm{\&xi;}}\\ {\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}^{\mathrm{tr}}=\mathrm{\&beta;}\stackrel{\&CenterDot;}{\mathrm{\&xi;}}\end{array}\right.$

其中，ξ为马氏体转变量；

为相变塑性应变；S为偏应力；为相变体积应变增量；β为相变膨胀系数；K为相变塑性系数；T为温度。

对于动力显示模型，本发明中，动力显示模型包括应力的导数张量模型或阶段流动应力模型。当动力显示模型为应力的导数张量模型时，根据用户输入的工艺参数信息计算得到待模拟模具的动力显示模型的步骤可以表示为：

${\stackrel{\&dtri;}{\mathrm{\&sigma;}}}_{\mathrm{ij}}={\stackrel{\&OverBar;}{D}}_{\mathrm{ijkl}}^{\mathrm{ep}}({\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_{\mathrm{kl}}\left(t\right)-{\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_{\mathrm{kl}}^{\mathrm{th}}\left(t\right)-{\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_{\mathrm{kl}}^{\mathrm{tr}}\left(t\right)-{\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_{\mathrm{kl}}^{\mathrm{tp}}\left(t\right))+{\stackrel{\&CenterDot;}{\stackrel{\&OverBar;}{D}}}_{\mathrm{ijkl}}^{\mathrm{ep}}({\mathrm{\&epsiv;}}_{\mathrm{kl}}\left(t\right)-{\mathrm{\&epsiv;}}_{\mathrm{kl}}^{\mathrm{th}}\left(t\right)-{\mathrm{\&epsiv;}}_{\mathrm{kl}}^{\mathrm{tr}}\left(t\right)-{\mathrm{\&epsiv;}}_{\mathrm{kl}}^{\mathrm{tp}}\left(t\right))$

其中，

是弹塑性本构矩阵，是以相邻构型为参考构型的Green应变率张量，是以相邻构型为参考构型的热应变率张量，

是以相邻构型为参考构型的相变体积应变率张量，

是以相邻构型为参考构型的相变塑性应变率张量。

当动力显示模型为段流动应力模型时，根据用户输入的工艺参数信息计算得到待模拟模具的动力显示模型的步骤可以表示为：

$\mathrm{\&sigma;}(\mathrm{\&epsiv;},\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}},T)=K^{\&prime;}{(b+\mathrm{\&epsiv;})}^{n\left(T\right)}{\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}^{m\left(T\right)}\exp ({\mathrm{\&beta;}}^{\&prime;}/T)$

其中，K′为强度系数，b为应变修正量，β′为温度系数，n(T)与m(T)满足：

n(T)＝n₀exp(-c_n(T_i-T₀))

m(T)＝m₀exp(-c_m(T_i-T₀))

其中，n₀、m₀、c_n、c_m分别为硬化材料常数，其值通过上述拉伸实验确定；T₀为温度初始值，T_i为温度当前值。

进一步地，本发明在根据用户输入的工艺参数信息计算得到待模拟模具的动力显示模型的步骤之后，还可以包括对动力显示模型进行离散化并显示的步骤，其中，对动力显示模型进行离散化的步骤可以表示为：

${\mathrm{\&sigma;}}_{\mathrm{ij}}^{n+1}={\widehat{\mathrm{\&sigma;}}}_{\mathrm{ij}}^{n+1}+{\left({\stackrel{\&OverBar;}{D}}_{\mathrm{ijkl}}^{\mathrm{ep}}\right)}^n({\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_{\mathrm{kl}}^{n+\frac{1}{2}}\left(t\right)-{\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_{\mathrm{th}}^{{\mathrm{th}}^{n+\frac{1}{2}}}\left(t\right)-{\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_{\mathrm{kl}}^{{\mathrm{tr}}^{n+\frac{1}{2}}}\left(t\right)-{\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_{\mathrm{kl}}^{{\mathrm{tp}}^{n+\frac{1}{2}}}\left(t\right))\mathrm{\&Delta;}{t}_n+{\left({\stackrel{\&CenterDot;}{\stackrel{\&OverBar;}{D}}}_{\mathrm{ijkl}}^{\mathrm{ep}}\right)}^n({\mathrm{\&epsiv;}}_{\mathrm{kl}}^{n+\frac{1}{2}}\left(t\right)-{\mathrm{\&epsiv;}}_{\mathrm{kl}}^{{\mathrm{th}}^{n+\frac{1}{2}}}\left(t\right)-{\mathrm{\&epsiv;}}_{\mathrm{kl}}^{{\mathrm{tr}}^{n+\frac{1}{2}}}\left(t\right)-{\mathrm{\&epsiv;}}_{\mathrm{kl}}^{{\mathrm{tp}}^{n+\frac{1}{2}}}\left(t\right))\mathrm{\&Delta;}{t}_n$

其中，Δt_n为时间增量。

进一步地，本发明在根据用户输入的工艺参数信息计算得到待模拟模具的动力显示模型的步骤之后，还可以包括计算下一时刻节点总位移矢量并显示的步骤，其中，计算下一时刻节点总位移矢量的步骤可以表示为：

$u_{n+1}=u_n+{\stackrel{\&CenterDot;}{u}}_{n+\frac{1}{2}}\mathrm{\&Delta;}{t}_n$

其中，u_n+1为t_n+1时刻节点的位移向量，u_n为t_n时刻节点的位移向量，

为t_n+1时刻节点的速度向量，Δt_n为时间增量。

更进一步地，为了保证运算的稳定性，时间增量Δt_n优选为：

$\mathrm{\&Delta;}{t}_n\&le;r{L}_n/c,L_n=\min (L_n^e,e-\mathrm{1,2},...\mathrm{ne})$

其中，常数r值介于0.8到0.9之间，优选为0.5，c为待模拟材料的内部声速，可由

求得，而E(T)为随温度变化的弹性模量，

为第n状态单元e的名义长度，ne为单元总数。

对于接触及摩擦模型，本发明中，接触及摩擦模型包括动力平衡运动模型或法向接触反力和切向摩擦力模型。当接触及摩擦模型为动力平衡运动模型时，根据用户输入的工艺参数信息计算得到接触及摩擦模型的步骤可以表示为：

$\frac{M}{\mathrm{\&Delta;}{t}_n^2}[u_{n+1}-(1+\mathrm{\&beta;}1)u_n+\mathrm{\&beta;}1u_{n-1}]+\frac{C}{\mathrm{\&Delta;}{t}_n}(u_{n+1}-u_n)+f_{n+1}^{\mathrm{int}}=f_{n+1}^{\mathrm{ext}}+{\mathrm{\&lambda;}}_{n+1}$

其中，Δt_n＝t_n+1-t_n，β1＝Δt_n/Δt_n-1，M为质量矩阵，C为阻尼矩阵，且满足C＝αM，α为阻尼系数；λ_n+1为t_n+1时刻摩擦接触反力，

为t_n+1时刻节点内力，u_n+1为t_n+1时刻节点位移向量，

为t_n+1时刻节点外力，u_n为t_n时刻节点位移向量，u_n-1为t_n-1时刻节点位移向量。

当接触及摩擦模型为法向接触反力和切向摩擦力模型时，根据用户输入的工艺参数信息计算得到接触及摩擦模型的步骤可以表示为：

其中，

为接触力控制系数，且满足

其中，

为高温下待模拟材料的屈服强度，

为室温下待模拟材料的屈服强度，m为待模拟材料上一点的质量，d_N为该点的法向穿透距离，d_T为该点沿模具表面的切向滑移距离。

对于热成形温度场模型，本发明中，根据用户输入的工艺参数信息计算得到热成形温度场模型的步骤可以表示为：

$(\frac{11}{6\mathrm{\&Delta;t}}{C}^e+K^e)T_{t+\mathrm{\&Delta;t}}=F^e+\frac{1}{\mathrm{\&Delta;t}}{C}^e(3T_t-1.5T_{t-\mathrm{\&Delta;t}}+\frac{1}{3}{T}_{t-2\mathrm{\&Delta;t}})$

其中，K^e为热传导系数，C^e为阻尼矩阵，F^e为温度载荷矩阵，T_t+Δt为t+Δt时刻的温度值，T_t为t时刻的温度值，T_t-Δt为t-Δt时刻的温度值，T_t-2Δt为t-2Δt时刻的温度值。

本发明提供的高强度钢板温热成形数值模拟方法实现了对温热成形过程中，待模拟材料的数值模拟功能，可以对温热成形模具设计及工艺参数的选取进行自动分析，进而提供优化方案，相较于现有的人工选取方式，可靠性及精确度更高，节省了制程成本。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，根据本发明的技术方案及其发明构思加以等同替换或改变，都应涵盖在本发明的保护范围之内。

Claims (10)

1.一种高强度钢板温热成形数值模拟方法，其特征在于，所述方法包括以下步骤：

根据用户输入的材料参数或模型调用指令，建立待模拟模具的分析模型；

接收用户输入的工艺参数信息，并根据用户输入的所述工艺参数信息计算得到待模拟模具的热、力及相变耦合模型，动力显示模型，接触及摩擦模型，和/或热成形温度场模型后显示。

2.如权利要求1所述的高强度钢板温热成形数值模拟方法，其特征在于，所述根据用户输入的所述工艺参数信息计算得到待模拟模具的热、力及相变耦合模型的步骤表示为：

$\left\{\begin{array}{c}\mathrm{\&xi;}=1-\exp [-\mathrm{\&theta;}(M_S-T)]\\ {\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}^{\mathrm{tp}}=\mathrm{KS}(1-\mathrm{\&xi;})\stackrel{\&CenterDot;}{\mathrm{\&xi;}}\\ {\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}^{\mathrm{tr}}=\mathrm{\&beta;}\stackrel{\&CenterDot;}{\mathrm{\&xi;}}\end{array}\right.$

其中，ξ为马氏体转变量；

为相变塑性应变；S为偏应力；

为相变体积应变增量；β为相变膨胀系数；K为相变塑性系数；T为温度；Ms为马氏体相变开始温度；θ为马氏体转变速率参数、k为相变塑性系数；

所述马氏体相变开始温度Ms、马氏体转变速率参数θ、相变塑性系数k满足：

$\left\{\begin{array}{cc}{M}_s=M_{s0}& \stackrel{\&OverBar;}{\mathrm{\&sigma;}}\&le;{\mathrm{\&sigma;}}_{a0}\\ M_s=M_{s0}+\mathrm{\&chi;}\&times;(\stackrel{\&OverBar;}{\mathrm{\&sigma;}}-{\mathrm{\&sigma;}}_{a0})& {\mathrm{\&sigma;}}_{a0}<\stackrel{\&OverBar;}{\mathrm{\&sigma;}}<{\mathrm{\&sigma;}}_{p0}\\ M_s=M_{s0}+\mathrm{\&chi;}\&times;({\mathrm{\&sigma;}}_{p0}-{\mathrm{\&sigma;}}_{a0})& \stackrel{\&OverBar;}{\mathrm{\&sigma;}}\&GreaterEqual;{\mathrm{\&sigma;}}_{p0}\end{array}\right.$

$\left\{\begin{array}{cc}\mathrm{\&theta;}={\mathrm{\&theta;}}_0& \stackrel{\&OverBar;}{\mathrm{\&sigma;}}\&le;{\mathrm{\&sigma;}}_a\\ \mathrm{\&theta;}={\mathrm{\&theta;}}_0+\mathrm{\&lambda;}\&times;(\stackrel{\&OverBar;}{\mathrm{\&sigma;}}-{\mathrm{\&sigma;}}_a)& {\mathrm{\&sigma;}}_a<\stackrel{\&OverBar;}{\mathrm{\&sigma;}}<{\mathrm{\&sigma;}}_p\\ \mathrm{\&theta;}={\mathrm{\&theta;}}_0+\mathrm{\&lambda;}\&times;({\mathrm{\&sigma;}}_p-{\mathrm{\&sigma;}}_a)& \stackrel{\&OverBar;}{\mathrm{\&sigma;}}\&GreaterEqual;{\mathrm{\&sigma;}}_p\end{array}\right.$

$\left\{\begin{array}{cc}K=K_0& \stackrel{\&OverBar;}{\mathrm{\&sigma;}}\&le;{\mathrm{\&sigma;}}_a^{\&prime;}\\ K=K_0+\mathrm{\&gamma;}\&times;(\stackrel{\&OverBar;}{\mathrm{\&sigma;}}-{\mathrm{\&sigma;}}_a^{\&prime;})& {\mathrm{\&sigma;}}_a^{\&prime;}<\stackrel{\&OverBar;}{\mathrm{\&sigma;}}<{\mathrm{\&sigma;}}_p^{\&prime;}\\ K=K_0+\mathrm{\&gamma;}\&times;({\mathrm{\&sigma;}}_p^{\&prime;}-{\mathrm{\&sigma;}}_a^{\&prime;})& \stackrel{\&OverBar;}{\mathrm{\&sigma;}}\&GreaterEqual;{\mathrm{\&sigma;}}_p^{\&prime;}\end{array}\right.$

其中，

为应力，M_s0为马氏体点基点温度，χ为相变温度与应力之间的耦合系数，σ_a0、σ_p0分别为马氏体相变点与应力耦合关系的区间应力两基点；θ₀为马氏体转变速率基点，λ为相变速率与应力之间的耦合系数，σ_a、σ_p分别为马氏体相变速率与应力耦合关系的区间应力两基点；K₀为相变塑性系数基点，γ为相变塑性系数与应力之间的耦合系数。

3.如权利要求1所述的高强度钢板温热成形数值模拟方法，其特征在于，所述动力显示模型为应力的导数张量模型，所述根据用户输入的所述工艺参数信息计算得到动力显示模型的步骤表示为：

${\stackrel{\&dtri;}{\mathrm{\&sigma;}}}_{\mathrm{ij}}={\stackrel{\&OverBar;}{D}}_{\mathrm{ijkl}}^{\mathrm{ep}}({\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_{\mathrm{kl}}\left(t\right)-{\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_{\mathrm{kl}}^{\mathrm{th}}\left(t\right)-{\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_{\mathrm{kl}}^{\mathrm{tr}}\left(t\right)-{\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_{\mathrm{kl}}^{\mathrm{tp}}\left(t\right))+{\stackrel{\&CenterDot;}{\stackrel{\&OverBar;}{D}}}_{\mathrm{ijkl}}^{\mathrm{ep}}({\mathrm{\&epsiv;}}_{\mathrm{kl}}\left(t\right)-{\mathrm{\&epsiv;}}_{\mathrm{kl}}^{\mathrm{th}}\left(t\right)-{\mathrm{\&epsiv;}}_{\mathrm{kl}}^{\mathrm{tr}}\left(t\right)-{\mathrm{\&epsiv;}}_{\mathrm{kl}}^{\mathrm{tp}}\left(t\right))$

其中，

是弹塑性本构矩阵，

是以相邻构型为参考构型的Green应变率张量，

是以相邻构型为参考构型的热应变率张量，是以相邻构型为参考构型的相变体积应变率张量，

是以相邻构型为参考构型的相变塑性应变率张量。

4.如权利要求1所述的高强度钢板温热成形数值模拟方法，其特征在于，所述动力显示模型为段流动应力模型，所述根据用户输入的所述工艺参数信息计算得到动力显示模型的步骤表示为：

$\mathrm{\&sigma;}(\mathrm{\&epsiv;},\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}},T)=K^{\&prime;}{(b+\mathrm{\&epsiv;})}^{n\left(T\right)}{\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}^{m\left(T\right)}\exp ({\mathrm{\&beta;}}^{\&prime;}/T)$

其中，K′为强度系数，b为应变修正量，β′为温度系数，n(T)与m(T)满足：

n(T)＝n₀exp(-c_n(T_i-T₀))

m(T)＝m₀exp(-c_m(T_i-T₀))

其中，n₀、m₀、c_n、c_m分别为硬化材料常数，其值通过上述拉伸实验确定；T₀为温度初始值，T_i为温度当前值。

5.如权利要求1所述的高强度钢板温热成形数值模拟方法，其特征在于，所述接触及摩擦模型为动力平衡运动模型，所述根据用户输入的所述工艺参数信息计算得到接触及摩擦模型的步骤表示为：

$\frac{M}{\mathrm{\&Delta;}{t}_n^2}[u_{n+1}-(1+\mathrm{\&beta;}1)u_n+\mathrm{\&beta;}1u_{n-1}]+\frac{C}{\mathrm{\&Delta;}{t}_n}(u_{n+1}-u_n)+f_{n+1}^{\mathrm{int}}=f_{n+1}^{\mathrm{ext}}+{\mathrm{\&lambda;}}_{n+1}$

其中，Δt_n＝t_n+1-t_n，β1＝Δt_n/Δt_n-1，M为质量矩阵，C为阻尼矩阵，且满足C＝αM，α为阻尼系数；λ_n+1为t_n+1时刻摩擦接触反力，

为t_n+1时刻节点内力，u_n+1为t_n+1时刻节点位移向量，

为t_n+1时刻节点外力，u_n为t_n时刻节点位移向量，u_n-1为t_n-1时刻节点位移向量。

6.如权利要求1所述的高强度钢板温热成形数值模拟方法，其特征在于，所述接触及摩擦模型为法向接触反力和切向摩擦力模型，所述根据用户输入的所述工艺参数信息计算得到接触及摩擦模型的步骤表示为：

其中，

为接触力控制系数，且满足

其中，

为高温下待模拟材料的屈服强度，

为室温下待模拟材料的屈服强度，m为待模拟材料上一点的质量，d_N为所述点的法向穿透距离，d_T为所述点沿模具表面的切向滑移距离。

7.如权利要求1所述的高强度钢板温热成形数值模拟方法，其特征在于，所述根据用户输入的所述工艺参数信息计算得到热成形温度场模型的步骤表示为：

$(\frac{11}{6\mathrm{\&Delta;t}}{C}^e+K^e)T_{t+\mathrm{\&Delta;t}}=F^e+\frac{1}{\mathrm{\&Delta;t}}{C}^e(3T_t-1.5T_{t-\mathrm{\&Delta;t}}+\frac{1}{3}{T}_{t-2\mathrm{\&Delta;t}})$

其中，K^e为热传导系数，C^e为阻尼矩阵，F^e为温度载荷矩阵，T_t+Δt为t+Δt时刻的温度值，T_t为t时刻的温度值，T_t-Δt为t-Δt时刻的温度值，T_t-2Δt为t-2Δt时刻的温度值。

8.如权利要求1至7任一项所述的高强度钢板温热成形数值模拟方法，其特征在于，在所述根据用户输入的所述工艺参数信息计算得到待模拟模具的动力显示模型的步骤之后，还包括对所述动力显示模型进行离散化并显示的步骤，所述对动力显示模型进行离散化的步骤表示为：

${\mathrm{\&sigma;}}_{\mathrm{ij}}^{n+1}={\widehat{\mathrm{\&sigma;}}}_{\mathrm{ij}}^{n+1}+{\left({\stackrel{\&OverBar;}{D}}_{\mathrm{ijkl}}^{\mathrm{ep}}\right)}^n({\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_{\mathrm{kl}}^{n+\frac{1}{2}}\left(t\right)-{\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_{\mathrm{th}}^{{\mathrm{th}}^{n+\frac{1}{2}}}\left(t\right)-{\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_{\mathrm{kl}}^{{\mathrm{tr}}^{n+\frac{1}{2}}}\left(t\right)-{\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_{\mathrm{kl}}^{{\mathrm{tp}}^{n+\frac{1}{2}}}\left(t\right))\mathrm{\&Delta;}{t}_n+{\left({\stackrel{\&CenterDot;}{\stackrel{\&OverBar;}{D}}}_{\mathrm{ijkl}}^{\mathrm{ep}}\right)}^n({\mathrm{\&epsiv;}}_{\mathrm{kl}}^{n+\frac{1}{2}}\left(t\right)-{\mathrm{\&epsiv;}}_{\mathrm{kl}}^{{\mathrm{th}}^{n+\frac{1}{2}}}\left(t\right)-{\mathrm{\&epsiv;}}_{\mathrm{kl}}^{{\mathrm{tr}}^{n+\frac{1}{2}}}\left(t\right)-{\mathrm{\&epsiv;}}_{\mathrm{kl}}^{{\mathrm{tp}}^{n+\frac{1}{2}}}\left(t\right))\mathrm{\&Delta;}{t}_n$

其中，Δt_n为时间增量，

是弹塑性本构矩阵，

是以相邻构型为参考构型的Green应变率张量，

是以相邻构型为参考构型的热应变率张量，是以相邻构型为参考构型的相变体积应变率张量，

是以相邻构型为参考构型的相变塑性应变率张量。

9.如权利要求8所述的高强度钢板温热成形数值模拟方法，其特征在于，所述根据用户输入的所述工艺参数信息计算得到待模拟模具的动力显示模型的步骤之后，还包括计算下一时刻节点总位移矢量并显示的步骤，所述计算下一时刻节点总位移矢量的步骤表示为：

$u_{n+1}=u_n+{\stackrel{\&CenterDot;}{u}}_{n+\frac{1}{2}}\mathrm{\&Delta;}{t}_n$

其中，u_n+1为t_n+1时刻节点的位移向量，u_n为t_n时刻节点的位移向量，

为t_n+1时刻节点的速度向量，Δt_n为时间增量。

10.如权利要求9所述的高强度钢板温热成形数值模拟方法，其特征在于，所述时间增量Δt_n满足：

$\mathrm{\&Delta;}{t}_n\&le;r{L}_n/c,L_n=\min (L_n^e,e-\mathrm{1,2},...\mathrm{ne})$

其中，常数r值介于0.8到0.9之间，优选为0.5，c为待模拟材料的内部声速，可由

求得，而E(T)为随温度变化的弹性模量，为第n状态单元e的名义长度，ne为单元总数。

Images