镁合金综合性能优良，在工业领域具有广阔的应用前景，温成形技术在镁合金成形中的应用也越来越受到关注，与等温、常温成形工艺相比，差温成形更有利于提高成形件质量。

数值仿真方法能有效预测成形件的缺陷^[1]，而准确的本构参数是保证数值仿真结果可靠的前提。直接通过准静态实验获得的材料参数，由于没有考虑复杂加载情况和温度的影响而会引入较大误差^[2]。采用分离式Hopkinson压杆装置能够大幅提高获取的材料本构参数的精度，但昂贵、复杂的实验设备使得这类方法并不能广泛地应用^[3]。Clausen等^[4]将数值模拟与实验结合，在优化迭代过程中不断调用有限元模型从而反求出材料参数，但是耗时很长，效率低下。乔良等^[5]采用RBF近似模型代替有限元模型，结合粒子群算法对本构参数实现了反求。李恩颖等^[6]利用支持向量机结合试验直接进行本构参数反求，其结果表明支持向量机在小样本非线性问题中具有良好的拟合精度。Wang等^[7-8]对利用近似模型反求参数及相关不确定性因素等问题进行了研究，提高了反求方法的稳健性。Goel等^[9]通过实验论述了多个近似模型组合使用解决实际问题的优势。Han等^[10]将量子计算与传统遗传算法相结合并引入量子旋转门，通过对组合优化问题求解，取得了比传统遗传算法更好的效果。本文选取NUMISHEET2011标准考题中的十字杯形件作为研究对象，采用基于有限元思想的数值方法，通过建立以本构参数为设计变量，目标温度为响应的混合近似模型，利用改进的量子遗传优化算法寻求仿真结果与试验数据之间误差函数的最小解^[11-12]，从而更准确地获取适用于差温成形过程的AZ31B镁合金的本构参数。

1 基于SVR-ELM混合近似模型的反求模型 1.1 改进的量子遗传算法

量子遗传算法是量子通信与遗传算法相结合的一种进化算法。它采用量子比特来表达一个基因，因此比经典遗传算法具有更好的特征多样性与收敛性。采用多量子比特编码p个参数的基因如下：

$ q_{^\mathit{j}}^\mathit{t} = \left( \begin{array}{l} \left| \begin{array}{l} \alpha _{11}^t\\ \beta _{11}^t \end{array} \right.\left| \begin{array}{l} \alpha _{12}^t\\ \beta _{12}^t \end{array} \right.\left| \begin{array}{l} \cdots \\ \cdots \end{array} \right.\left| \begin{array}{l} \alpha _{1\mathit{k}}^\mathit{t}\\ \beta _{1\mathit{k}}^\mathit{t} \end{array} \right.\left| \begin{array}{l} \alpha _{21}^t\\ \beta _{21}^t \end{array} \right.\left| \begin{array}{l} \alpha _{22}^t\\ \beta _{22}^t \end{array} \right.\left| \begin{array}{l} \cdots \\ \cdots \end{array} \right.\\ \left| \begin{array}{l} \alpha _{{\rm{2}}\mathit{k}}^\mathit{t}\\ \beta _{{\rm{2}}\mathit{k}}^\mathit{t} \end{array} \right.\left| \begin{array}{l} \cdots \\ \cdots \end{array} \right.\left| \begin{array}{l} \alpha _{p1}^t\\ \beta _{p1}^t \end{array} \right.\left| \begin{array}{l} \alpha _{p{\rm{2}}}^t\\ \beta _{p{\rm{2}}}^t \end{array} \right.\left| \begin{array}{l} \cdots \\ \cdots \end{array} \right.\left| \begin{array}{l} \alpha _{pk}^t\\ \beta _{pk}^t \end{array} \right. \end{array} \right) $

(1)

式中：q_j^t指第t代中第j个个体的染色体编码，α和β表示量子比特编码，p为染色体基因个数，k为编码基因的量子比特数。

量子遗传算法的迭代更新通过量子旋转门实现：

$ \left[\begin{array}{l} \alpha _i^*\\ \beta _i^* \end{array} \right] = \mathit{\boldsymbol{G}}\left( {{\theta _i}} \right)\left[\begin{array}{l} {\alpha _i}\\ {\beta _i} \end{array} \right] $

(2)

$ \mathit{\boldsymbol{G}}\left( {{\theta _i}} \right) = \left[\begin{array}{l} {\rm{cos}}{\theta _i}\;\;-{\rm{sin}}{\theta _i}\\ {\rm{sin}}{\theta _i}\;\;\;{\rm{cos}}{\theta _i} \end{array} \right] $

(3)

式中:(α_i, β_i)^T, (α_i^*, β_i^*)^T分别表示编码染色体第i个量子比特旋转门更新前后的概率幅; θ_i为量子旋转门旋转角。

利用动态量子门代替固定旋转角量子门，对量子遗传算法进行改进。根据进化进程实时调整量子旋转门旋转角θ_i的大小，迭代过程中将通过测量染色体所得的二进制解与当前最优解进行逐位比较，根据比较结果以及进化程度，调整个体中相应位量子比特，使个体更快速地向最优解方向进化^[13]。

设通过测量染色体得到的某二进制解为：

$ d = {\lambda _1}{\lambda _2} \ldots {\lambda _s} $

(4)

式中：λ₁λ₂…λ_s表示长度为s的二进制串，λ_i(i=1, 2, …，s)表示染色体上第i位编码。

为方便与最优解比较，通过反正弦函数将λ_i向角度映射，则更新后染色体第i位对应的旋转角可表示为：

$ {\theta _i} = \left( {\frac{{{\rm{arc}}\;{\rm{sin}}{\lambda _i} - \theta _i^*}}{g} + \delta } \right)\theta _i^* $

(5)

式中：θ_i^*是当前第i位编码对应的旋转角; g为当前进化代数; δ为平衡参数，避免随进化代数增加而使更新后旋转角过小导致收敛速度减慢。

1.2 SVR-ELM混合模型建立

不同近似模型针对不同工程问题的预测能力表现各异，精度高的近似模型不一定预测能力强，精度低的近似模型也可能得到更合适的最优解，选择某一种近似模型可能会存在一定的不足，而由多个近似模型组合而成的混合模型对高维度非线性问题的预测精度更高更稳定。

$ {y_{avg}} = \sum\limits_{\mathit{i} = 1}^\mathit{N} {{\phi _\mathit{i}}\left( x \right){y_i}\left( x \right)} $

(6)

$ \sum\limits_{\mathit{i} = 1}^\mathit{N} {{\phi _\mathit{i}}\left( x \right)} = 1 $

(7)

式中：y_avg为混合模型的输出响应，ϕ_i(x)为第i个模型的权系数，y_i(x)为第i模型的输出响应，N为近似模型的个数。

构造混合模型的关键是权系数的计算，根据各模型的预测精度来分配权系数。考虑到不同模型的精度有较大差别，采用启发式计算方法来计算各自的权系数。

$ {E_{avg}} = \frac{1}{N}\sum\limits_{\mathit{i} = 1}^\mathit{N} {{E_i}} $

(8)

$ \phi _i^* = {\left( {{E_i} + {\gamma _1} \cdot {E_{{\rm{avg}}}}} \right)^{{\gamma _2}}},{\gamma _1} < 1,{\gamma _2} < {\rm{ }}0 $

(9)

$ {\phi _\mathit{i}} = \phi _i^*/\sum\limits_{\mathit{j} = 1}^\mathit{N} {\phi _j^*} $

(10)

式中:E_i为第i模型的均方误差; E_avg为系统平均均方误差; γ₁, γ₂为权重控制系数, 为保证混合模型的稳定性，此处取γ₁=0.05, γ₂=-1。

基于支持向量回归机(support vector regression，SVR)与极限学习机(extreme learning machine，ELM)两种回归近似模型，提出一种SVR-ELM混合近似模型用于反求镁合金本构参数。

1) 支持向量回归机。

支持向量机是以结构风险最小化原则为基础的一种基于统计学理论的新型机器学习方法，回归型支持向量机的基本思想是寻找一个最优分类超平面, 使得所有样本到该超平面的距离最小，即超平面回归函数与样本之间的误差最小。

利用对应约束条件以及Largrange函数，通过求解其对偶形式，得到回归方程：

$ \begin{array}{l} f\left( \mathit{\boldsymbol{x}} \right) = {\mathit{\Psi }^*}\varphi \left( \mathit{\boldsymbol{x}} \right) + {\mathit{\boldsymbol{b}}^*} = \\ \sum\limits_{\mathit{i} = 1}^\mathit{l} {\left( {{\tau _i} - \tau _i^*} \right)} K({\mathit{\boldsymbol{x}}_i}, \mathit{\boldsymbol{x}}) + {\mathit{\boldsymbol{b}}^*} \end{array} $

(11)

式中：τ=[τ₁, τ₂, …, τ_l], τ^*=[τ₁^*, τ₂^*, …, τ_l^*]为最优解；K(x_i, x)=φ(x_i)φ(x)为核函数；x_i为输入响应; l为训练样本数量；b^*为包含支持向量的待定参数。

2) 极限学习机回归模型。

极限学习机是由单隐含层前馈神经网络(single-hidden layer feed forward neural network，SLFN)发展而来^[14]，相比传统神经网络, 其学习速度更快，收敛效果更好。

设有Q个不同样本(x_i, y_i)，i=1, 2, …, Q, 其中x_i=[x_i1, x_i2, …, x_in]^T∈Rⁿ, y_i=[y_i1, y_i2, …, y_im]^T∈R^m; Rⁿ, R^m分别为与x_i，y_i维数对应的实数空间。

对具有n个输入变量，m个输出变量及η个隐含层节点的单隐含层神经网络可表示为:

$ {\mathit{\boldsymbol{y}}_\mathit{j}} = \sum\limits_{i = 1}^\eta {{\beta _\mathit{i}}\mathit{g}({\omega _\mathit{i}}\cdot{\mathit{\boldsymbol{x}}_\mathit{j}} + {\mathit{\boldsymbol{b}}_\mathit{i}}), \mathit{j} = 1, 2, \cdots, Q} $

(12)

式中：β_i, ω_i，b_i分别为输出权值、输入权值和隐含层阈值；g(ω_i·x_j+b_i)为激活函数。

为训练SLFN，需要确定ω_i, b_i, β，使得:

$ \xi = \mathop {{\rm{min}}}\limits_{{\omega _\mathit{i}}, {\mathit{\boldsymbol{b}}_i}, \beta } \left\| {\mathit{\boldsymbol{H}}\left( {{\omega _1}, {\omega _2}, \cdots, {\omega _\eta }, {\mathit{b}_1}, {\mathit{b}_2}, \cdots, {\mathit{b}_\eta }} \right)\beta - \mathit{\boldsymbol{Y}}} \right\| $

(13)

式中:H为g(ω_i·x_j+b_i)构成的隐含层输出矩阵；Y=[y₁, y₂, …, y_Q]_m×Q为神经网络的期望输出。

ELM训练时随机产生输入层与隐含层间的连接权值ω及隐含层神经元的阈值b，且整个训练过程中无需调整，只需要设置隐含层神经元的个数，便可以获得唯一的最优解。由ω与b可唯一确定出隐含层输出矩阵H，进而更快速获得构建神经网络所需的ω, b, β三个参数。通过对式(13) 的转化，神经网络可简述为线性系统。

$ \mathit{\boldsymbol{H}}\beta = \mathit{\boldsymbol{Y}} $

(14)

求解该方程可得隐含层与输出层连接权值的最小二乘解β=H⁺Y，其中H⁺是H的Moore-Penrose广义逆。

1.3 SVR-ELM混合近似模型的自适应更新

近似模型训练时，样本空间的大小与多样性直接制约着近似模型的拟合精度。采用量子遗传算法对混合近似模型求解域寻优并通过自适应方法将产生的局部最优解选择性地加入各模型初始样本库，重新拟合相应的近似模型并再次组合新的混合模型，通过不断地更新样本库，使得混合近似模型的求解域更快速地逼近目标最优解。具体操作如下：

1) 通过初始样本库分别训练SVR, ELM两模型并计算各自的均方误差E_svr, E_elm作为模型精度的判断依据。

2) 建立SVR-ELM混合近似模型，并利用量子遗传算法寻优，将迭代寻优得到的局部最优样本与初始样本库暂时组合成新样本库，用于重新训练SVR与ELM，并再次计算各自的均方误差E′_svr, E′_elm以及ΔE_svr=E′_svr-E_svr, ΔE_elm=E′_elm-E_elm。

3) 若ΔE_svr＜0, ΔE_elm＜0，将寻优得到的样本加入初始样本库用于后续的SVR, ELM训练；若ΔE_svr≥0, ΔE_elm＜0，将该样本只加入ELM的初始样本库；若ΔE_svr＜0, ΔE_elm≥0，将该样本只加入SVR的初始样本库；否则舍弃该样本并返回至步骤2)。

1.4 基于SVR-ELM混合近似模型的算法流程

1) 使用拉丁超立方试验设计方法进行抽样并组成初始样本库，对样本点进行仿真计算出响应值。利用初始样本库分别训练SVR和ELM模型。

2) 建立SVR-ELM混合近似模型并利用启发式算法计算SVR和ELM的权系数。

3) 采用量子遗传算法对初始种群测量寻优，记录当前代最优个体。

4) 通过自适应方法将获得的个体选择性地加入对应的样本库。

5) 利用更新的样本库重新拟合两模型并建立混合模型，利用量子遗传算法迭代寻优，判断是否达到收敛条件，满足条件就结束，否则返回至步骤4)。

该算法流程如图 1。

1.5 SVR-ELM混合近似模型的验证

为了验证SVR-ELM混合模型的效果，选取非线性函数对其测试，并与单一SVR模型以及ELM模型对比。此处采用的测试函数为：

$ \mathit{f}\left( {\mathit{x}, \mathit{y}} \right) = \frac{{{\rm{sin}}\sqrt {{\mathit{x}^2} + {\mathit{y}^2}} }}{{{\mathit{x}^2} + {\mathit{y}^2}}} + {{\rm{e}}^{\frac{{{\rm{cos}}2\pi \mathit{x} + {\rm{cos}}2\pi \mathit{y}}}{2}}}, - 2 < x, y < 2 $

(15)

函数图像如图 2所示。

常用的近似模型精度评价指标有均方根误差(root mean square error, RMSE)和决定系数R², 定义如下：

$ {\rm{RMSE = }}\frac{1}{k}\sqrt {\sum\limits_{i = 1}^k {{{\left( {{\mathit{f}_i} - {\widehat {{\mathit{f}_\mathit{i}}}}} \right)}^2}} } $

(16)

$ {\mathit{R}^2} = \sum\limits_{i = 1}^k {{{\left( {{\widehat {{\mathit{f}_\mathit{i}}}} - {{\mathit{\bar f}}_\mathit{i}}} \right)}^2}/\sum\limits_{i = 1}^k {{{\left( {{\mathit{f}_i} - {{\mathit{\bar f}}_\mathit{i}}} \right)}^2}} } $

(17)

式中：k为样本数，f_i为真实响应值，$ \widehat {{\mathit{f}_\mathit{i}}}$为样本点对应的近似值，f_i为真实响应值均值。RMSE值越小，R²值越接近1，表明模型拟合精度越高。

通过拉丁超立方抽样方法对测试函数随机抽样16组，分别采用SVR，ELM，SVR-ELM三种模型对样本点进行拟合。利用均方根误差RMSE与决定系数R²评价模型拟合精度，结果如表 1所示，图 3为预测结果对比。

测试结果表明，相比SVR，ELM两种模型，SVR-ELM混合近似模型均方根误差更小，决定系数R²更接近1，具有更好的拟合精度。

2 SVR-ELM混合近似模型在镁合金本构参数反求中的应用 2.1 有限元模型的建立与验证

选取NUMISHEET2011中的Benchmark 2十字杯形件为研究对象，对镁合金AZ31B十字杯形件差温成形过程进行仿真。冲压成形的工件如图 4。

考虑到计算时间和十字杯形件的两向对称性，选取1/4杯形件，在ABAQUS中建立对称边界条件有限元模型，材料为镁合金AZ31B，板料尺寸如图 5所示，板厚为0.5 mm, 材料参数如表 2所示。

通过模具的作用使板料温度场差异化是差温成形的主要特征，将凹模和压边圈加热到250 ℃，通过热传导方式提高板料法兰部位的温度，从而降低材料流动阻力，板料、凸模和顶料板加热到100 ℃，使得处于拉深部位的板料拥有相对较低的塑性变形能力，从而有效控制该部位材料的减薄率。利用有限元软件ABAQUS进行差温成形仿真，设拉深深度为18 mm, 摩擦系数为0.05。成形仿真过程采用动力显式分析方法，温度—位移加载，板料网格划分如图 6所示，DE，DF代表截面位置，网格类型采用C3D8RT(八节点温度位移耦合减缩积分单元)，沿厚度方向设置4层单元以控制沙漏问题，图 7为有限元仿真模型示意图。

为了验证该有限元模型的有效性，利用表 3中镁合金AZ31B的Johnson-Cook本构模型关键参数^[15]进行数值模拟并将仿真结果中的板料DE截面温度分布情况与试验数据进行对比，仿真结果如图 8所示，对比结果如图 9所示。

由图 8可知，法兰处最高温度为229.233 ℃，与试验数据相比误差为7.44%。

由图 9对比结果可知，仿真数据与试验结果基本吻合，证明建立的有限元模型有效。

2.2 Johnson-Cook本构模型

Johnson-Cook本构模型(J-C模型)由应变强化项、应变率强化项和温度软化项三部分组成，它适用于描述金属材料在不同应变率、不同温度下的材料变形行为。表达式如下：

$ \sigma = \left( {\mathit{A} + \mathit{B}{\varepsilon ^\mathit{n}}} \right)\left( {1 + \mathit{C}\cdot{\rm{ln}}\frac{{\dot \varepsilon }}{{{{\dot \varepsilon }_0}}}} \right)\left( {1 - {\mathit{T}^*}^\mathit{m}} \right) $

(18)

$ {{\mathit{T}}^{*}}=\frac{\mathit{T}-{{\mathit{T}}_{\rm{room}}}}{{{T}_{\rm{melt}}}-{{\mathit{T}}_{\rm{room}}}} $

(19)

式中：A，B，C，m，n为待求参数；$ {{\mathit{\dot \varepsilon }}_{\rm{0}}}$为参考应变率；T, T_room, T_melt分别为实验温度、室温和材料熔点。

J-C模型适用于准静态变形到大应变率的高速变形，针对不同的金属材料，一般通过准静态实验和霍普金森杆实验来获取J-C模型本构参数，但由图 9对比结果可知用此法获得的本构模型不能很好地描述AZ31B镁合金板差温成形过程。

2.3 基于SVR-ELM混合模型的本构参数反求

鉴于NUMISHEET2011提供的十字杯形件试验数据，在本模型中选取沿板料DE截面方向温度最高为响应目标。以本构模型中A，B，C，m，n为设计变量，仿真值与试验值误差最小化为优化目标。使用拉丁超立方试验设计方法抽样20组，组成初始样本库，结果如表 4所示，对样本点进行仿真计算得到响应值。选取前15组数据作为训练集，后5组作为检验集。为更好地对比反求结果，分别采用SVR，ELM和SVR-ELM三种近似模型结合量子遗传算法对AZ31B差温成形的Johnson-Cook本构参数进行反求，算法最大迭代次数为50。

2.4 反求结果与分析

基于SVR，ELM，SVR-ELM三种不同的近似模型，通过改进的量子遗传算法寻优得到Johnson-Cook本构参数结果，如表 5所示。

将基于不同近似模型反求得到的3组数据应用到十字杯形件的差温成形仿真中，成形件温度场与应力场云图分别如图 10与图 11所示。3组仿真计算得到的温度最高值与试验值对比如表 6所示。根据温度场云图显示，成形件最高温度均分布在沿DE截面方向的法兰边缘。

结果表明在较少的算法迭代次数下，SVR-ELM混合近似模型即可获得最优解，收敛速度比2个单一近似模型快。与试验值相比，利用SVR-ELM混合近似模型反求的本构模型仿真所得的最高温度的相对误差仅为0.39%，精度最高。在基于3个不同近似模型反求的本构参数所得到的应力云图中，凹模底部圆角与侧壁半径较小的圆角附近区域应力较大，这主要是因为差温成形过程中，与凸模接触的板料部分温度较低，材料流变应力较大，随着拉深过程进行，处在圆角半径较小位置的材料还要承受较大的弯曲应力。

为进一步验证反求的本构参数的准确性，选取仿真结果中的DE，DF截面(如图 12所示)的厚度分布与试验数据进行对比，对比结果分别如图 13、图 14所示。

从图 13和图 14结果可知板料在模具侧壁位置处由于拉深产生了较大的减薄，板料法兰区域DE截面位置由于材料流动不均造成材料堆积，增厚现象比较严重。

由上述结果可以看出，成形件温度场分布非常不均匀，进一步导致应力场分布不均，从而影响成形件厚度减薄率与成形质量。运用反求得到的本构方程，通过一定方法优化相关工艺，提高成形件质量是接下来的首要任务。

3 结论

本文通过自适应方法更新样本库，采用SVR-ELM混合近似模型与量子遗传寻优算法建立了反求模型来求解镁合金AZ31B差温成形的Johnson-Cook本构参数。

1) 利用改进的量子遗传算法和自适应样本更新方法，提高了近似模型的拟合精度，加快了寻优算法向最优求解域收敛的速度。与采用单一近似模型相比，混合近似模型具有更好的泛化能力与更高的精度。

2) 基于NUMISHEET2011十字杯形件的差温成形有限元模型，利用SVR-ELM混合模型反求本构参数的结果误差较小，避免了常规准静态实验获得的材料参数不能很好描述冲压成形等高应变率问题，提高了镁合金差温成形数值模拟的精度，对预测成形质量具有重要的意义。