高聚物黏结炸药(polymer bonded explosive, PBX)是以高能单质炸药为主体，加入黏结剂等物质的固体炸药，在常规武器战斗部中具有广泛应用^[1]。PBX在制造、运输及发射等过程中经历的复杂应力过程，影响着武器弹药的安全和使用。因此，对PBX的力学行为进行有限元模拟研究具有重要意义^[2]。PBX力学行为的有限元模拟依赖于本构模型的正确性，一个合适本构模型中参数的准确性对于PBX的有限元模拟至关重要。材料本构的参数识别属于带有约束条件寻找最优解的优化问题，针对该问题，目前已有一系列的研究报道^[3-4]。从研究现状来看，PBX炸药本构模型具有强非线性特点，其参数无法通过对标准拉伸等实验数据直接处理获得。另外，参数空间具有非凸性，使参数搜索可能陷入局部最优值。所以，传统优化方法不能保证获得待识别参数的全局最优解。

随着智能分析方法的发展，遗传算法、蚁群算法等被不断引入参数识别中。其中，遗传算法在处理复杂非线性问题具有独特的优势，得到了广泛的应用。B.M.Chaparro等^[5]、G.H.Majzoobi等^[6]通过遗传算法，对多种本构模型参数进行了识别；P.A.Muñoz-Rojas等^[7]结合自编有限元程序和遗传算法，进行参数识别；陈炳瑞等^[8]将遗传算法和神经网络结合，进行岩体流变参数的识别；高军等^[9]通过遗传算法，对PBX本构模型参数进行了识别。然而遗传算法虽然简单、通用，但算法本身受遗传操作参数、群体规模等条件影响，客观存在着早熟收敛等问题，最终可能无法获得最优解。解决这些问题的一个有效途径是调整遗传算法的搜索方式。随着大规模并行机的快速发展，结合遗传算法固有的并行性，遗传算法的多种群并行化成为解决这些问题的有效方法^[10]。

本文中，采用一种改进的多种群并行遗传算法来取代传统的遗传算法，并结合有限元软件ABAQUS，建立材料本构模型参数识别方法。并利用该方法，对PBX炸药黏弹性损伤本构模型参数进行模拟识别，与标准遗传算法识别结果进行比较。

1.   数学模型

本构模型的参数识别是，在已知材料本构模型的基础上，结合同种材料试件的相关实验数据，合理地选择材料本构模型参数，使材料模型可以真实地模拟材料的物理特性。本构模型中的待识别参数，可以表示为^[11]：

式中：m是待识别参数的个数。

通过实验测量可以获得选定测点处的实验值，如应变、位移等。测点处实验值可表示为：

式中：u₁⁰、u₂⁰、…、u_n⁰为测点处实测值；n为测量数据的个数。

测点处值是材料参数变量X=(x₁, x₂, …, x_m)^T的函数，则有：

测点处的模拟计算值，可以在给定待识别参数X=(x₁, x₂, …, x_m)^T的初始值时通过数值计算获得，表示为：

式中：u₁、u₂、…、u_n为测点处模拟值。

参数识别就是寻找一组材料参数，使模拟计算结果和测量结果之间的差值最小。假设计算两者之间差值的目标函数为：

则参数识别的过程就是寻找一组材料参数，使下式成立：

这组材料参数X^*即为本构模型参数的最优解。通常，材料参数在一定的允许范围内约束条件为：x_{i, min}≤x_i≤x_{i, max}(i=1, 2, …, m)，x_{i, min}和x_{i, max}为材料参数的允许范围。

2.   多种群遗传算法

遗传算法(standard genetic algorithm, SGA)是根据生物界自然选择和自然遗传机制发展而来的随机搜索算法。它基于群体运算的方法，对每个个体进行的各种运算都具有一定的相互独立性，所以具有一种隐含的并行性^[12]。根据这种并行性，遗传算法的并行结构分为单种群主从并行结构、细粒度并行结构、多种群并行结构和分级混合并行结构。其中多种群并行结构计算效率高，结构实际操作简单，因此该结构被广泛应用^[13-14]。通过多种群并行结构对遗传算法进行并行化处理，形成多种群遗传算法(multiple-population genetic algorithm, MPGA)。

多种群遗传算法的算法结构如图 1所示。图中种群1~N均是基本的遗传算法，N个种群均运行基本的遗传算法，同时进行优化搜索。进化过程中，通过移民操作进行各种群之间的信息交换，将源种群中的最优个体利用移民算子定期地迁移到目标种群，并同目标种群的最劣个体比较后淘汰较差个体。移民操作将各种群相互联系，否则MPGA相当于独立进行多次SGA计算。

图  1  MPGA算法结构图

Figure  1.  Structure of MPGA

下载: 全尺寸图片 幻灯片

为了更好地模拟生物进化过程，并行进化时对各种群选择不同的遗传操作参数，实现不同的搜索目的。在SGA中，交叉和变异概率的恰当选择决定着全局和局部搜索能力的均衡。一般建议取较大的交叉概率(0.7~0.9)和较小的变异概率(0.001~0.1)。对于不同的取值组合，优化结果也具有较大的差异。MPGA通过多个赋以不同遗传操作参数的种群协同进化, 在取值范围内随机产生多个遗传参数组合，兼顾算法的全局搜索和局部搜索。

进化过程中，利用人工选择算子选出各种群的最优个体放入精华种群进行保存，并从中选取全局最优个体作为算法的最优解。精华种群同时作为判断算法终止的条件，采用最优个体最少保持代数作为终止判据^[15]。

3.   算法收敛速度改进策略

基本遗传算法中，交叉和变异概率均为一个固定常量，不同适应度的个体被选择进行交叉和变异操作的概率相同，使子代个体的产生随机性较强，父代的优良基因易被破坏，对算法的收敛和计算速度产生不利的影响。本文中，设计了一种与父代种群整体适应度相关联的自适应交叉和变异概率：

式中：P_c0和P_m0为各种群中预先给定的交叉和变异概率；P_{c, min}和P_{m, min}为最小交叉和变异概率，这里分别取0.7和0.001；F_avg为父代种群中个体的平均适应度，F_max为父代种群中个体的最大适应度。F_avg和F_max计算公式如下：

这种自适应变化概率能够根据遗传进化过程中父代种群的整体适应度不断地自动调整交叉和变异概率，可以增大遗传算法获得全局最优解的可能性，同时可以提高遗传算法的收敛速度。将自适应交叉和变异概率引入多种群遗传算法中，形成改进的多种群遗传算法(improved multiple-population genetic algorithm, IMPGA)。

4.   参数识别的实现

基于多种群遗传算法的本构模型参数识别方法由正演和反演两个部分组成。正演部分是在参数识别过程中，计算一组参数所对应的模拟值。本文中，通过有限元软件ABAQUS的二次开发接口UMAT，将本构模型嵌入到ABAQUS进行计算。反演部分由多种群遗传算法组成，本文中采用MATLAB语言进行计算程序的编写。具体实现步骤如下。

(1) 初始化。种群中个体采用实数编码，在各待识别参数范围内，随机产生种群数N组初始种群。种群数和种群中群体规模根据实际情况选择。每个染色体个体(x₁, x₂, …, x_m)为m维向量，其中x_i(i=1, 2, …, m)为待识别参数。随机产生N组交叉和变异概率，分别赋予各组种群。

(2) 正演计算。调用ABAQUS命令程序流，利用真实模拟试验件几何尺寸和实验条件的有限元模型，分别计算N组种群中每个参数个体对应的模拟计算值。

(3) 个体适应度计算。利用测点处实验值和模拟计算值，通过式(5)计算每个个体的目标函数值。本文中参数识别为求最小值问题，因此适应度亦可取式(5)，则适应度为：

式中：u_i⁰是实验值，u_i是模拟值，X是模拟计算对应的参数，i是数据个数。

(4) 移民操作。根据适应度大小判断个体优劣，将源种群中的最优个体迁移到目标种群，并同目标种群的最劣个体比较，再淘汰较差个体。同时根据适应度，通过人工选择算子选出各种群的最优个体放入精华种群MaxChrom中。

(5) 终止判断。选出精华种群MaxChrom中的最优个体，同上一代中选出的最优个体进行比较。若两者相同，最优个体保持代数gen=gen+1，否则令gen=0，同时将本次选出的最优个体代替全局最优个体进行保存。最优个体保持代数gen同设定的最少保持代数MINGEN比较，若gen>MINGEN则识别过程终止，此时的全局最优个体为参数识别的最优解；否则，继续下一步。

(6) 对各种群中的个体分别执行选择、交叉、变异操作，生成下一代种群，再返回到步骤(2)进行新一轮的优化。遗传操作中采用轮盘赌选择、实数交叉和随机变异方法。

具体计算流程如图 2所示。

图  2  参数识别的计算流程

Figure  2.  Calculation process of parameter identification

下载: 全尺寸图片 幻灯片

5.   参数识别方法的比较

通过基于多种群遗传算法MPGA和IMPGA的参数识别方法，对PBX黏弹性损伤本构模型参数进行模拟识别，验证该方法的正确性和可靠性，并同基于基本遗传算法SGA的参数识别方法的结果进行比较。

5.1   PBX黏弹性损伤本构模型

Visco-scram方程被广泛应用于PBX炸药的力学行为描述。在该方程中，引入表示内部细观缺陷的损伤变量c，表示微裂纹的特征尺寸(平均尺寸)，则PBX炸药黏弹性损伤本构模型^[9]为：

式中：σ_m为体积应力，S为偏应力，η为黏性系数，G为剪切模量，K为体积模量，β_m与β_s是和微裂纹扩展相关的两个常数。基于文献[16]，损伤变量c的演化方程具有幂律形式：

式中：χ与n是两个材料常数。本构模型中，共涉及7个材料参数，分别为杨氏模量E、泊松比ν、黏性系数η及4个常数β_m、β_s、χ与n。通过ABAQUS二次开发用户材料子程序UMAT，将PBX黏弹性损伤本构模型嵌入ABAQUS中，应用于PBX材料。

参数识别正演部分分析模型采用PBX炸药圆柱体试件模拟压缩，试件长度6 mm，直径10 mm，单元类型为C3D10M四面体实体单元。一端固定，另一端通过位移加载形式进行压缩。压缩加载最大位移0.2 mm，加载时间1 s。取待识别参数的一组合理值为真实值，如表 1。利用这组参数进行模拟压缩，得到选定测点的数据，将此数据作为实验数据进行本构模型的参数识别。通过比较参数识别结果和参数真实值的差距，分析识别方法的优劣。

表  1  参数识别结果

Table  1.  Results of parameter identification

识别参数	真实值	取值范围	SGA识别			MPGA识别			IMPGA识别
			结果	误差/%		结果	误差/%		结果	误差/%
E/MPa	496.1	400~600	507.142	2.22		504.931	1.78		505.306	1.86
ν	0.38	0.2~0.6	0.372	2.17		0.387	1.93		0.373	1.84
η	260	200~300	252.460	2.90		263.38	1.30		263.256	1.25
βm	0.6	0.3~0.9	0.619	3.17		0.613	2.09		0.587	2.16
βs	0.4	0.1~0.7	0.384	3.86		0.414	3.40		0.413	3.25
χ	155.5	100~200	148.814	4.30		149.358	3.95		149.249	4.02
n	0.628	0.3~0.9	0.643	2.38		0.638	1.73		0.639	1.75

下载: 导出CSV 

| 显示表格

参数识别中，避免多解的一般方法是，预先确定待识别参数的取值范围，将参数取值约束在变化范围内。但是，在实际操作中取值范围难于直接确定，只能根据理论分析或工程经验来预估。通常可以在符合理论前提下暂定较大的取值范围，通过参数识别确定一组参数，再进行其他工况下的验证。PBX本构模型参数取值范围，见表 1。

参数识别中，实验数据的类型和质量直接影响参数识别问题的适定性，即解的存在性、稳定性和唯一性问题。因此，实验数据需要尽量选择对本构参数敏感度较大的，另外测点位置也要根据实际问题进行选择。本文中，选取圆柱体加载端表面一点为测点，取模拟载荷值作为实验数据用于参数识别，该点处载荷-位移曲线如图 3所示。

图  3  测点处的载荷-位移曲线

Figure  3.  Load-displacement curve of measurement point

下载: 全尺寸图片 幻灯片

5.2   算例结果比较

分别采用基于SGA的参数识别方法、本文中建立的基于MPGA和IMPGA的参数识别方法，对PBX黏弹性损伤本构模型进行参数识别。

基于SGA参数识别方法中的操作参数取值分别为：种群规模20，迭代次数20次，交叉概率0.8，变异概率0.1。图 4为SGA方法运行5次得到的种群最优解适应度随进化代数的变化曲线。从图 4可以看出，5次得到的适应度均不相同，相应的识别得到的参数也各不相同。其中一次识别过程中适应度较大时已停止进化，说明算法陷入某个局部收敛区。更改遗传操作交叉和变异概率后进行多次识别计算，最终适应度均不能收敛到一个稳定的范围。所以，对于待识别参数较多且具有强非线性的本构模型，通过标准遗传算法进行参数识别时结果并不稳定，需要计算多次取相对最优解。

图  4  SGA方法适应度

Figure  4.  Fitness value with SGA

下载: 全尺寸图片 幻灯片

采用MPGA和IMPGA方法进行参数识别时, 取种群数N为10, 每个种群的规模均为10。各种群的初始交叉和变异概率分别在区间(0.7~0.9)和(0.001~0.1)中随机选取，IMPGA方法中的交叉和变异概率在进化过程中进行自适应调整，算法终止条件是最优个体最少保持代数为10。采用两种方法计算运行5次得到最优解的适应度变化曲线，如图 5~6所示。从图 5~6可以看出，两种方法运行5次的适应度结果都能趋近于稳定值, 说明稳定性很好。和SGA方法相比，多种群遗传算法在更少的进化次数下收敛到一个稳定值，收敛速度快而且寻优能力更强。IMPGA方法和MPGA方法相比，适应度从进化初期快速下降，进化代数也更少，基本在10代内已趋于稳定。IMPGA方法收敛速度更快，寻优能力更强，计算效率也更高。

图  5  MPGA方法适应度

Figure  5.  Fitness value with MPGA

下载: 全尺寸图片 幻灯片

图  6  IMPGA方法适应度

Figure  6.  Fitness value with IMPGA

下载: 全尺寸图片 幻灯片

对5次采用SGA、MPGA和IMPGA方法的参数识别结果中最优解进行对比，见表 1。从表中可以看出，采用3种方法的结果均满足工程计算基本要求，但是采用多种群遗传方法的结果明显优于采用SGA方法的，MPGA和IMPGA方法的结果较接近。将通过3种方法获得的参数代入正演分析模型中，计算参数对应的测点载荷和位移，图 7为计算结果和实验数据的载荷-位移曲线对比。从图 7可以看出，3种方法识别得到的参数模拟曲线和实验曲线整体上都比较吻合, 但MPGA和IMPGA方法的结果优于SGA方法的。采用SGA的参数识别方法在优化结果上不稳定，但是在进行多次识别计算后，也可以获得较优解。而采用MPGA和IMPGA的参数识别方法，运行5次的结果最优值近乎相等，算法稳定性更好，识别效率更高。

图  7  载荷-位移曲线对比

Figure  7.  Contrast of load-displacement curves

下载: 全尺寸图片 幻灯片

可见，多种群遗传算法MPGA中，通过多个种群同时对参数空间进行搜索，兼顾了全局和局部搜索能力，降低了算法对遗传操作参数的敏感性，对克服未成熟收敛有显著的效果，同时算法寻优结果比标准遗传算法SGA的稳定性更好，计算效率更高。多种群遗传算法中，IMPGA比MPGA算法收敛速度更快、寻优能力更强。

6.   结论

通过多种群并行结构对标准遗传算法进行并行化处理，并引入移民算子和精华种群形成多种群遗传算法。利用移民算子实现各种群之间的信息交换，采用精华种群作为最优个体储存空间和算法终止判据。针对算法收敛速度的改进，设计了与父代种群整体适应度相关联的自适应交叉和变异概率。自适应交叉和变异率随着遗传进化过程不断的自动调整，可以增大遗传算法获得全局最优解的可能性。

基于改进的多种群遗传算法和ABAQUS软件，建立了材料本构模型参数识别方法。通过对PBX黏弹性损伤本构模型参数进行模拟识别，验证了该方法的正确性和可靠性，并同基于标准遗传算法的参数识别方法的结果进行比较。算例证明，基于改进多种群遗传算法的参数识别方法识别结果的稳定性更好，对克服未成熟收敛有显著的效果。同时，该方法收敛速度更快，寻优能力更强，计算效率更高，适用于复杂非线性材料参数的识别。