考虑成形损伤的抗撞零件仿真设计方法

摘要

本发明公开了一种考虑成形损伤的抗撞零件仿真设计方法，以解决现有技术存在的车身高强度抗撞零件仿真设计时忽略零件不同部位成形损伤继承特性所导致的后续仿真计算精度不足及零件结构设计不合理的问题，步骤为：建立车身高强度抗撞零件热成形损伤准则；车身高强度抗撞零件热成形模拟；车身高强度抗撞零件抗撞性能仿真评估。本方法综合考虑了成形温度、成形速率、摩擦等因素的零件仿真成形后的“成形损伤值”作为零件抗撞性能评估时考虑的因素，显著提高了后续仿真计算的精度，减少了反复修改的次数，确保了零件结构设计的合理性，大大降低了抗撞零件仿真设计的难度以及对设计人员的要求。

主权项

考虑成形损伤的抗撞零件仿真设计方法，其特征在于步骤如下：步骤一、建立车身高强度抗撞零件热成形损伤准则，具体过程为：1）材料高温单向拉伸试验利用热模拟试验机对硼钢试件进行一系列高温单向拉伸试验，试验前利用点焊机将K型热电偶丝（1）的一端焊在每个试件（2）上表面的中央，K型热电偶丝（1）另一端保持自由释放，试验中先将试件（2）夹紧于热模拟试验机的夹具内，同时将K型热电偶丝（1）自由释放的一端与热模拟试验机相连，随后，对热模拟试验机内部空间进行抽真空处理，利用电阻加热方式实现试件（2）的加热过程，而冷却过程中则通过调节压缩空气的流量控制试件（2）的冷却速度，具体试验方案如下：（1）将试件（2）以5℃/s的加热速率加热至925℃后保温3min，确保试件（2）的微观组织完全奥氏体化；（2）以50℃/s的冷却速率使试件（2）依次降至变形温度600℃、700℃、800℃，并在各变形温度下保温5s使试件（2）温度均匀稳定；（3）在设定的变形温度600℃、700℃、800℃和变形应变率0.01s^‑1、0.1s^‑1、1s^‑1、10s^‑1下对试件（2）进行拉伸，直至断裂破坏，断裂后对得到的拉伸后试件（8）进行空冷，整个拉伸过程中热模拟试验机会同时记录载荷随时间变化、温度随时间变化的曲线，整个高温单向拉伸试验包含3个变形温度和4个变形应变率组合成的12组试验条件，分别为：变形温度600℃和变形应变率0.01s^‑1、变形温度600℃和变形应变率0.1s^‑1、变形温度600℃和变形应变率1s^‑1、变形温度600℃和变形应变率10s^‑1、变形温度700℃和变形应变率0.01s^‑1、变形温度700℃和变形应变率0.1s^‑1、变形温度700℃和变形应变率1s^‑1、变形温度700℃和变形应变率10s^‑1、变形温度800℃和变形应变率0.01s^‑1、变形温度800℃和变形应变率0.1s^‑1、变形温度800℃和变形应变率1s^‑1、变形温度800℃和变形应变率10s^‑1，每组试验条件下进行一次高温拉伸试验，将试验中热模拟试验机测得的载荷随时间变化的曲线F(t)按照公式（1）换算为试件（2）的名义应力随时间变化的曲线σ_nom(t)，将CCD摄像机（3）测得的试件（2）的标距段长度随时间变化的曲线ΔL(t)按照公式（2）换算为试件（2）的标距段名义应变随时间变化的曲线ε_nom(t)，按照公式（3）将试件（2）的名义应力随时间变化的曲线σ_nom(t)换算为试件（2）的真实应力随时间变化的曲线σ_true(t)，按照公式（4）将试件（2）的标距段名义应变随时间变化的曲线ε_nom(t)换算为试件（2）的标距段真实应变随时间变化的曲线ε_true(t)，并消去两曲线中的时间变量t，以真实应变ε_true为自变量，真实应力σ_true为因变量，得到每一试验条件下的真实应力‑应变曲线σ_true(ε_true)：${\sigma }_{nom}\left(t\right)=\frac{F\left(t\right)}{A_0}---\left(1\right)$式中：F(t)为载荷随时间变化的曲线；A₀为试件标距段原始横截面积；σ_nom(t)为试件的名义应力随时间变化的曲线。${ϵ}_{nom}\left(t\right)=\frac{\Delta L\left(t\right)}{L_0}---\left(2\right)$式中：ΔL(t)为试件的标距段长度随时间变化的曲线；L₀为试件标距段原始长度；ε_nom(t)为试件的标距段名义应变随时间变化的曲线。σ_true(t)=σ_nom(t)(1+ε_nom(t))     (3)式中：σ_nom(t)为试件的名义应力随时间变化的曲线；ε_nom(t)为试件的标距段名义应变随时间变化的曲线；σ_true(t)为试件的真实应力随时间变化的曲线。ε_true(t)=ln(1+ε_nom(t))      (4)式中：ε_nom(t)为试件的标距段名义应变随时间变化的曲线；ε_true(t)为试件的标距段真实应变随时间变化的曲线。2）建立基于成形损伤的本构方程：（1）建立基于成形损伤的本构方程，以考虑热成形时材料的损伤，具体表达式如下：${\stackrel{·}{ϵ}}_e^P={\left(\frac{{\sigma }_e/(1-f_{d1})-H-k}{K}\right)}^{n_1}{(1-f_{d1})}^{-{\gamma }_1}---\left(5\right)$${\stackrel{·}{ϵ}}_{ij}^P=\frac{3S_{ij}}{2{\sigma }_e}{\stackrel{·}{ϵ}}_e^P---\left(6\right)$$\stackrel{·}{\stackrel{‾}{\rho }}=A(1-\bar{\rho })\left|{\stackrel{·}{ϵ}}_e^P\right|-C{\bar{\rho }}^{n_2}---\left(7\right)$$H=B{\bar{\rho }}^{n_0}---\left(8\right)$${\stackrel{·}{f}}_{d1}={\mathrm{D\sigma }}_e\left|{\stackrel{·}{ϵ}}_e^P\right|/{(1-f_{d1})}^{{\gamma }_2}---\left(9\right)$${\sigma }_{ij}=D_{ijkl}(1-f_{d1})({ϵ}_{kl}^T-{ϵ}_{kl}^P)---\left(10\right)$$D_{ijkl}=\frac{E}{2(1+v)}({\delta }_{il}{\delta }_{jk}+{\delta }_{ik}{\delta }_{jl})+\frac{Ev}{(1+v)(1-2v)}{\delta }_{ij}{\delta }_{kl}---\left(11\right)$式中：是成形时的等效塑性应变率；σ_e是成形时的等效应力；H是成形时由位错引起的应变强化；成形损伤变量f_d1，其变化范围为0～1，f_d1＝0时表示成形时材料没有损伤，f_d1＝1时表示成形时材料完全失效；为成形时的塑性应变率分量；S_ij为成形时的偏应力分量；ρ_i为材料初始状态下的位错密度，ρ_m为成形时材料可达到的最大位错密度，且ρ_i≤ρ≤ρ_m，即σ_ij是成形时的应力张量分量；是成形时的总应变张量分量；是成形时的塑性应变张量分量；D_ijkl是四阶刚度张量分量；E是杨氏模量；υ是泊松比；δ_ij为克罗内克因子，下标i、j、k、l变化范围为1～3，重复下标遵循爱因斯坦求和约定。参数k、K、n₁、B、C、D、E是与温度相关的材料参数，定义如下：k＝k₀exp(Q_k/RT)        (12)K＝K₀exp(Q_K/RT)          (13)n₁＝n₁₀exp(Q_n/RT)         (14)B＝B₀exp(Q_B/RT)          (15)C＝C₀exp(‑Q_C/RT)          (16)D＝D₀exp(Q_D/RT)          (17)E＝E₀exp(Q_E/RT)          (18)式中：R为通用气体常数；T为温度；Q为激活能。（2）确定本构方程中的材料参数：首先，建立求解问题的目标函数，然后对目标函数应用进化规划算法进行优化，最终确定本构方程中的材料参数，这里，需要确定的材料参数总共有20个，依次为：A、n₂、γ₁、γ₂、k₀、n₀、K₀、n₁₀、B₀、C₀、D₀、E₀、Q_k、Q_K、Q_n、Q_B、Q_C、Q_D、Q_E、R。根据拟合曲线与试验数据之间的距离建立目标函数：$f\left(x\right)=\frac{1}{n}\underset{i=1}{\overset{n}{\Sigma }}{w}_i{r}_i^2=\frac{1}{n}\underset{i=1}{\overset{n}{\Sigma }}{w}_i({\mathrm{\Delta \sigma }}^2+{\mathrm{\Delta ϵ}}^2)=\frac{1}{n}\underset{i=1}{\overset{n}{\Sigma }}{w}_i({({\sigma }_i^c-{\sigma }_i^e)}^2+{({ϵ}_i^c-{ϵ}_i^e)}^2)---\left(19\right)$式中：f(x)是关于20个未知材料参数所构成矢量x＝(A,n₂,...,Q_E,R)的实值函数；n为试验数据的总容量；w_i为第i数据点的权重值；分别为第i个试验数据对应的应力和应变值；分别为与第i个试验数据相对应的拟合曲线上的应力和应变值为使拟合曲线应变区间向试验数据收敛，将加入f(x)得到：$f\left(x\right)=\frac{1}{n}\underset{i=1}{\overset{n}{\Sigma }}{w}_i({({\sigma }_i^c-{\sigma }_i^e)}^2+{({ϵ}_i^c-{ϵ}_i^e)}^2)+W{({ϵ}_n^c-{ϵ}_n^e)}^2---\left(20\right)$式中：W为权重系数；为降低实际使用过程中确定权重值的难度及克服应力应变单位不一致的问题，最终建立的目标函数如下：$f\left(x\right)=\frac{1}{n}\underset{i=1}{\overset{n}{\Sigma }}(w{1}_i{\left(\ln \right(\frac{{\sigma }_i^c({ϵ}_i^e{ϵ}_n^c/{ϵ}_n^e)}{{\sigma }_i^e}\left)\right)}^2+w{2}_i{\left(\ln \right(\frac{{ϵ}_i^c({ϵ}_i^e{ϵ}_n^c/{ϵ}_n^e)}{{ϵ}_i^e}\left)\right)}^2)+{({ϵ}_n^c-{ϵ}_n^e)}^2---\left(21\right)$式中：$w{1}_i={\mathrm{nϵ}}_i^e/\underset{i=1}{\overset{n}{\Sigma }}{ϵ}_i^e(w{1}_1+w{1}_2+...+w{1}_n=n)---\left(22\right)$$w{1}_2={\mathrm{nϵ}}_i^e/\underset{i=1}{\overset{n}{\Sigma }}{ϵ}_i^e(w{2}_1+w{2}_2+...+w{2}_n=n)---\left(23\right)$采用改进的快速进化规划算法对建立的目标函数进行优化，确定全部材料参数，进化规划算法把目标函数作为生物种群，通过突变，选择产生新一代种群；重复这一过程，直到获得合乎要求的种群或规定的进化时限，详细的进化规划是一个迭代的过程：[1]取迭代计数k＝1，随机生成μ个种群，即随机输入μ组矢量对(x_i,η_i)，其中η_i为进化规划自适应策略参数，i＝1,2,3,...,μ；[2]对于每个个体矢量对(x_i,η_i)，计算f(x_i)；[3]对于每个父代矢量对(x_i,η_i)，生成两个子代矢量对和其中：$\left\{\begin{array}{c}{x}_i^1\left(j\right)=x_i\left(j\right)+{\eta }_i\left(j\right)N_j(0,1)\\ x_i^2\left(j\right)=x_i\left(j\right)+{\eta }_i\left(j\right){\delta }_j\\ {\eta }_i^1\left(j\right)={\eta }_i\left(j\right)\exp \left({\tau }^{1}N\right(0,1)+{\mathrm{\tau N}}_j(0,1))\end{array}\right.---\left(24\right)$计算并比较和的大小，取两者较小者所对应的矢量对，记为(x′_i,η′_i)，其中x_i(j)，x′_i(j)，η_i(j)，η′_i(j)分别为矢量x_i，x′_i，η_i，η′_i的第j个分量（j＝1,2,...,n，n为待优化材料参数的个数）；N(0,1)为服从一维标准正态分布的随机数；N_j(0,1)为服从一维标准正态分布对应于第j个分量的随机数；δ_j为服从柯西分布对应于第j个分量的随机数；参数τ¹和τ分别取和标准正态分布和柯西分布的密度函数分别为：$\begin{array}{cc}{f}_N\left(x\right)=\frac{1}{\sqrt{2n}}{e}^{-x^2/2},& -\infty <x<+\infty \end{array}---\left(25\right)$$\begin{array}{cc}{f}_{\delta }\left(x\right)=\frac{1}{\pi }\frac{1}{1+x^2},& -\infty <x<+\infty \end{array}---\left(26\right)$[4]对于所有的i＝1,2,3,...,μ，将所有的父代矢量对(x_i,η_i)和子代矢量对(x′_i,η′_i)作为一个整体，取出q个矢量对，然后，将所有父代和子代矢量对中的任意一个矢量对与取出的q个矢量对作比较，比较矢量对所对应的目标函数值，如果该矢量对小于q个矢量对中的某一个，则该矢量对得分加1，所有矢量对的最高得分为q，最低得分为0；[5]从2μ个矢量对中选出得分最高的μ个矢量对，作为下一次迭代的父代矢量对；[6]判断迭代结束条件是否满足；如不满足，则k＝k+1，并重复上述过程。（3）损伤本构方程的软件实现利用Fortran语言将已确定材料参数的成形损伤本构方程编写为Ls‑dyna用户自定义材料子程序，通过用户自定义材料子程序接口嵌入到有限元软件Ls‑dyna中。步骤二、车身高强度抗撞零件热成形模拟，具体过程为：1）利用Hypermesh有限元软件建立车身高强度抗撞零件热成形仿真模型，模型中的对象包括板料（6）、凸模（4）、凹模（7）及压边圈（5）；2）分别为模型中各对象赋予截面属性和材料特性，所有对象均采用壳单元截面属性，并在其中定义各对象的厚度t，而凸模（4）、凹模（7）及压边圈（5）采用刚体物理材料，并在其中定义各对象的密度ρ、弹性模量E及泊松比υ，而板料（6）则采用上述基于成形损伤的用户自定义物理材料，且所有对象均采用各向同性热材料，并在其中定义各对象的比热容HC及热传导系数TC；3）设置模型中各对象间的接触关系，定义各对象具有的温度场、运动特性、约束条件及模型计算所需的控制卡片；4）利用Ls‑dyna软件求解车身高强度抗撞零件热成形仿真模型，得到成形后的零件模型及其成形损伤云图和厚度分布云图。步骤三：车身高强度抗撞零件抗撞性能仿真评估：1）力学使用特性数据库（1）材料高温损伤试验利用Ls‑dyna软件对试件（2）进行虚拟试验，整个试件（2）采用壳单元截面属性，并在其中定义其厚度t，采用考虑成形损伤的用户自定义物理材料及各向同性热材料，并在其中定义其比热容HC及热传导系数TC，同时，在试件（2）中间区域定义l₀＝25mm长的均温段，在该段内施加恒定温度场，而在该段左右两侧区域则按照实际试验中试件（2）温度分布施加相应的温度场，虚拟试验时，约束试件（2）左端节点的全部6个自由度，在右端施加沿轴向随时间变化的强制位移，确保均温段以恒定的应变率进行变形，试验中使均温段取典型的试验条件：成形温度750℃，变形应变率0.1s^‑1。在该条件下，进行一系列虚拟试验，使均温段内单元的成形损伤值分别达到目标成形损伤值0,α₁,α₂,......,α_n，即得到一系列均温段具有不同成形损伤值的拉伸后试件（8），并分别测出每一试件（2）均温段长度的变化值0,利用试验机对试件（2）进行实际试验，实际试验条件与虚拟试验相同，试验中利用CCD摄像机（3）进行拍摄，并通过ARAMIS光学变形系统实时反馈试件（2）均温段的变形量，当各试件（2）均温段长度变化分别达到0,时，停止拉伸，此时，实际试验中均温段处材料达到与虚拟试验中同等损伤程度，试验中同一均温段长度变化值重复m次试验。随后，对各试件（2）均温段快速淬火冷却至室温，即得到一系列均温段具有不同成形损伤值的淬火后马氏体相的拉伸后试件（8）。（2）试件加工对材料高温损伤试验中获得的一系列淬火后的拉伸后试件（8）进行线切割加工，得到用于材料室温单向拉伸试验所用子试件（9），随后，用细砂纸对所有子试件（9）各表面进行轻微打磨，去除上面的氧化皮，并记录打磨后各子试件（9）标距段的最小截面积，作为下述材料室温单向拉伸试验中子试件（9）标距段的原始横截面积A₀。（3）材料室温单向拉伸试验利用电子万能试验机对材料高温损伤试验中获得的均温段具有不同损伤值的子试件（9）进行m种不同服役应变率的单向拉伸试验，所有子试件（9）均在室温下拉伸直至断裂，整个拉伸过程中试验机会记录载荷随时间变化的曲线F(t)，并将其按照公式（1）换算为子试件（9）的名义应力随时间变化的曲线σ_nom(t)，而采用拉伸引伸计测出子试件（9）的标距段长度随时间变化的曲线ΔL(t)，并将其按照公式（2）换算为子试件（9）的标距段名义应变随时间变化的曲线ε_nom(t)，按照公式（3）将子试件（9）的名义应力随时间变化的曲线σ_nom(t)换算为子试件（9）的真实应力随时间变化的曲线σ_true(t)，按照公式（4）将子试件（9）的标距段名义应变随时间变化的曲线ε_nom(t)换算为子试件（9）的标距段真实应变随时间变化的曲线ε_true(t)，并消去两曲线中的时间变量t，以真实应变ε_true为自变量，真实应力σ_true为因变量，得到每一试件（9）的真实应力‑应变曲线σ_true(ε_true)。最终，得具有不同成形损伤值、不同服役应变率的试件（9）的真实应力‑应变曲线σ_true(ε_true)，用于建立下述考虑成形损伤的服役本构方程。2）建立考虑成形损伤的服役本构方程建立考虑成形损伤的服役本构方程，具体表达式如下：${\stackrel{·}{ϵ}}^p={\left(\frac{\sigma /(1-f_{d2})-F{{\bar{\rho }}_s}^{n_a}-y}{Y}\right)}^{n_c}{\left(\frac{1}{1-f_{d2}}\right)}^{{\gamma }_3}---\left(27\right)$${\stackrel{·}{ϵ}}_{ij}^p=\frac{3S_{s-ij}}{2\sigma }{\stackrel{·}{ϵ}}^p---\left(28\right)$${\stackrel{·}{\stackrel{‾}{\rho }}}_s=Z(1-{\bar{\rho }}_s)\left|{\stackrel{·}{ϵ}}^p\right|-C{{\bar{\rho }}_s}^{{\gamma }_4}---\left(29\right)$${\stackrel{·}{f}}_{d2}={\beta }_1\left(\sigma {\stackrel{·}{ϵ}}^p\right){f_{d2}}^{{\gamma }_5}+\frac{{\beta }_2{\left({\stackrel{·}{ϵ}}^p\right)}^{{\gamma }_6}\cosh \left({\beta }_3{ϵ}^p\right)}{{(1-f_{d2})}^{{\gamma }_7}}---\left(30\right)$${\sigma }_{s-ij}=D_{s-ijkl}(1-f_{d2})({ϵ}_{kl}^t-{ϵ}_{kl}^p)---\left(31\right)$$D_{s-ijkl}=\frac{L}{2(1+v)}({\delta }_{il}{\delta }_{jk}+{\delta }_{ik}{\delta }_{jl})+\frac{Lv}{(1+v)(1-2v)}{\delta }_{ij}{\delta }_{kl}---\left(32\right)$式中：是服役时的等效塑性应变率；σ是服役时的等效应力；服役损伤变量f_d2，其变化范围为0～1，f_d2＝0时表示服役时材料没有服役损伤，f_d2＝1时表示服役时材料完全失效；为服役时的塑性应变率分量；S_s‑ij为服役时的偏应力分量；ρ_si为服役前材料的位错密度，ρ_sm为服役时材料可达到的最大位错密度，且ρ_si≤ρ_s≤ρ_sm，即σ_s‑ij是服役时的应力张量分量；是服役时的总应变张量分量；是服役时的塑性应变张量分量；D_s‑ijkl是四阶刚度张量分量；L是杨氏模量；υ是泊松比；δ_ij为克罗内克因子，下标i、j、k、l变化范围为1～3，重复下标遵循爱因斯坦求和约定。参数y、Y、F、G、L、β₁、β₂、β₃、γ₅、γ₆是与成形损伤相关的材料参数，定义如下：y＝y₀exp(W_y/f_d1)           (33)Y＝Y₀exp(W_Y/f_d1)          (34)F＝F₀exp(W_F/f_d1)            (35)G＝G₀exp(‑W_G/f_d1)           (36)L＝L₀exp(W_L/f_d1)           (37)${\beta }_1={\beta }_{10}\exp (W_{{\beta }_1}/f_{d1})---\left(38\right)$${\beta }_2={\beta }_{20}\exp (W_{{\beta }_2}/f_{d1})---\left(39\right)$${\beta }_3={\beta }_{30}\exp (W_{{\beta }_3}/f_{d1})---\left(40\right)$${\gamma }_5={\gamma }_{50}\exp (W_{{\gamma }_5}/f_{d1})---\left(41\right)$${\gamma }_6={\gamma }_{60}\exp (W_{{\gamma }_6}/f_{d1})---\left(42\right)$确定考虑成形损伤的服役本构方程中材料参数及服役本构方程软件实现所用的方法与建立基于成形损伤本构方程时所用的方法完全相同，而这里，需要确定的材料参数总共有25个，依次为：Z、γ₃、γ₄、γ₇、n_c、y₀、Y₀、F₀、G₀、L₀、β₁₀、β₂₀、β₃₀、γ₅₀、γ₆₀、W_y、W_Y、W_F、W_G、W_L、3）车身高强度抗撞零件虚拟试验（1）利用Hypermesh软件建立车身高强度抗撞零件虚拟试验模型，虚拟压溃试验模型中的对象包括：刚性壁障(10)、压溃薄壁抗撞零件（11)；虚拟弯曲试验模型中的对象包括：1号刚性圆滚（12）、弯曲薄壁抗撞零件（13）、2号刚性圆滚（14）、3号刚性圆滚（15）；（2）分别为模型中各对象赋予截面属性和材料特性，所有对象均采用壳单元截面属性，抗撞零件的厚度分布继承热成形后零件的厚度分布，刚性壁障(10)及所有刚性圆滚均采用刚体物理材料，并在其中定义各对象的密度ρ、弹性模量E及泊松比υ，而压溃薄壁抗撞零件（13）和弯曲薄壁抗撞零件（13）则采用考虑成形损伤的用户自定义服役本构物理材料；（3）设置模型中各对象间的接触关系，定义各对象的运动特性、约束条件及模型计算所需的控制卡片；（4）利用Ls‑dyna软件求解车身高强度抗撞零件虚拟试验模型，输出抗撞零件随时间变化的变形趋势云图、接触力曲线及吸能量曲线，并与设计目标值进行对比，以确认该抗撞零件是否符合设计要求。