NiTi形状记忆合金 （shape memory alloy, SMA）因其独特的形状记忆效应（shape memory effect, SME）和良好的超弹性，被广泛应用于航空航天、机械、电子和医疗器械等领域^[1-3]，是目前形状记忆合金中研究和应用最多的一种^[4]。近年来，NiTi合金由于其优异的吸能和减振特性，在武器装甲防护以及建筑物的减振防护领域被认为具备重要的应用潜能，因此其在冲击与高速切削等极端环境下的高应变率动态力学行为逐渐成为研究热点。

为了研究和应用NiTi形状记忆合金的超弹性与形状记忆效应，国内外学者近几十年来开展了深入的理论和试验研究，通过不同方法提出了多种形式的NiTi合金本构模型。其中基于自由能驱动力理论建立的Tanaka^[5]、Liang^[6]与Brinson^[7]系列本构模型在形式上比较相似，主要区别在于采用了不同的马氏体相变动力学模型。该系列本构模型形式简单，避开了自由能等难以测量的参量，得到了广泛的应用。Lagoudas等^[8-9]基于自由能与耗散势概念提出了可以描述超弹性与形状记忆效应的本构模型，该模型可以很好地描述应力诱发下的正、反方向的相变行为，但形式复杂计算量大，难以用于实际工程。Auricchio等^[10-12]以相变中马氏体百分比含量和马氏体百分比变化率作为两个独立的内变量发展了可描述超弹性的一维与三维本构关系，该模型宏观材料常数易测量，便于与有限元法结合并嵌入大型商用有限元软件。

上述的本构模型主要关注NiTi合金准静态条件下的形状记忆效应与超弹性力学行为，但对NiTi合金在高应力、高应变水平下的动态力学行为研究报道还比较有限。为了适应NiTi合金在冲击载荷等极端环境下的应用，本文基于不可逆热力学理论框架，确定了表征相变、塑性行为的两个内变量，提出了NiTi合金三维动态本构模型的主控方程，用于材料的动态响应数值模拟。由于材料动态响应模拟的准确性取决于其本构模型的准确性以及本构模型参数的精度。本文构建的动态本构模型共22个参数，若采用传统的单因素分析方法，对全部参数的识别计算量会非常大，并且容易忽略本构参数间的相互作用，进而影响结果的准确性。

本文采用拉丁超立方抽样方法对基于不可逆热力学理论框架确定的NiTi合金动态本构模型参数在整个参数空间中抽样，使用非参数统计方法中的Spearman秩相关分析法对本构参数随机输入样本集与其对应的目标函数输出结果集作相关性分析，建立Spearman秩相关系数等效求解参数敏感度的表达式，实现本构参数敏感度的整体性分析。并且利用改进遗传算法对NiTi合金动态本构模型的材料参数进行识别，提高参数识别的可靠性，为NiTi合金的实际工程应用奠定基础。

1 NiTi合金动态本构模型

以不可逆热力学理论框架为基础构建本构关系，首先需要确定内外变量，然后假设自由能函数（内、外变量的函数）^[13]。本文选用应变张量与温度作为外变量，因NiTi合金变形过程存在2个不可逆过程，假定2个内变量：$ \xi $表征应力诱发马氏体相变行为的内变量，$ \eta $表征塑性行为的内变量。单位质量的Helmholtz自由能函数Φ可表示为：

$ \mathit{\qquad\varPhi}=\mathit{\varPhi}\left(\varepsilon_{ij}^e,T,\xi,\eta\right) $

(1)

式中：$ \varepsilon _{ij}^{e} $为弹性应变张量；$ T $为热力学温度。

在等温条件下，由Clausius-Duhem不等式得到

$ \qquad\sigma_{ij}:\dot{\varepsilon}_{ij}-\rho\dot{\mathit{\varPhi}}\geqslant0 $

(2)

式中：$ \sigma_ij $为应力张量；$ \rho $为密度；总应变$ \varepsilon _{ij}^{} $为弹性应变与非弹性应变之和，而对于NiTi合金的非弹性应变由塑性应变$ \varepsilon _{ij}^{\rm{p}} $与相变应变$ \varepsilon _{ij}^{\rm{tr}} $组成。

将式（1）代入Clausius-Duhem不等式可得：

$\qquad \begin{split} &\left(\sigma_{i j}-\rho \dfrac{\partial \varPhi}{\partial \varepsilon_{i j}^e}\right) \dot{\varepsilon}_{i j}^e+\left(\sigma_{i j} \dot{\varepsilon}_{i j}^{\rm{t r}}-\rho \dfrac{\partial \varPhi}{\partial \xi} \dot{\xi}\right)+\\ &\qquad \left(\sigma_{i j} \dot{\varepsilon}_{i j}^p-\rho \dfrac{\partial \varPhi}{\partial \boldsymbol{\eta}} \dot{\boldsymbol{\eta}}\right) \geqslant 0 \end{split}$

(3)

假定加载过程中弹性应变、相变应变与塑性应变相互独立，要使式（3）恒成立须有：

$\qquad \sigma_{i j}=\rho \dfrac{\partial \varPhi}{\partial \varepsilon_{i j}^e} $

(4)

$ \qquad {\sigma }_{ij}\dot{\varepsilon }_{ij}^{\rm{tr}}-\rho \frac{\partial \varPhi }{\partial \mathbf{\xi }}\mathbf{\dot{\mathbf{\xi }}}\geqslant 0 $

(5)

$\qquad {\sigma }_{ij}\dot{\varepsilon }_{ij}^{\rm{p}}-\rho \frac{\partial \varPhi }{\partial \mathbf{\eta }}\mathbf{\dot{\mathbf{\eta }}}\geqslant 0 $

(6)

以上为NiTi合金的动态本构框架，式（4）～式（6）分别表征了弹性应力应变关系、相变演化规律与塑性演化规律。

根据材料的本构框架，相变行为内变量演化必须满足式（5）。本文类比经典Chaboche塑性本构模型^[14]构造与NiTi合金相变行为相关的广义力$ {A }\text{=}-\rho \frac{\partial \varPhi }{\partial \mathbf{\xi }} $与势函数$ \Theta \text{=}\Theta ({\sigma }_{ij},{A }) $，推导得到NiTi合金的相变演化规律如下，具体推导过程参见文献^[15]：

$\qquad \dot{\varepsilon }_{ij}^{\rm{tr}}={\dot{\lambda }}_{1}\frac{\partial \Theta }{\partial {\sigma }_{ij}}\text=\frac{3}{2}{\left(\frac{F_{y}^{\rm{tr}}}{{Z}_{1}}\right)}^{{{n}_{1}}}\frac{{s}_{ij}-{{{A}^{\prime}}}_{ij}}{\sigma _{{}_{\rm{eq}}}^{\rm{tr}}} $

(7)

$\qquad F_{y}^{\rm{tr}}={\sigma }_{\rm{eq}} - \sigma _{s}^{\rm{tr}}(n) = \sqrt{\frac{3}{2}({s}_{ij} - {{{A}^{\prime}}}_{ij})\colon ({s}_{ij} - {{{A}^{\prime}}}_{ij})} - \sigma _{s}^{\rm{tr}}(n) $

(8)

$\qquad \sigma _{s}^{\rm{tr}}(n)\text=(\text{1}-n)\sigma _{s}^{\rm{tr}}+n\sigma _{f}^{\rm{tr}} \text{，} n=\varepsilon _{\rm{eq}}^{\rm{tr}}/{\varepsilon }_{m} $

(9)

$\qquad \varepsilon_{\rm{eq}}^{\rm{t r}} = \sqrt{\frac{2}{3} \boldsymbol{\varepsilon}^{\rm{t r}}: \boldsymbol{\varepsilon}^{\rm{t r}}}, \quad \dot{\varepsilon}_{i j}^{\rm{t r}} = \dot{\lambda}_1 \frac{\partial \Theta}{\partial \sigma_{i j}} = \frac{3}{2}\left(\frac{F_y^{\rm{t r}}}{Z_1}\right)^{n_1} \frac{s_{i j} - A_{i j}^{\prime}}{\sigma_{\rm{eq}}^{\rm{t r}}} $

(10)

$\qquad {\dot{A}}_{ij}={k}_{1}\dot{\varepsilon }_{ij}^{\rm{tr}}-{k}_{2}\dot{\varepsilon }_{\rm{eq}}^{\rm{tr}}{A}_{ij} $

(11)

式中：$ \dot{\varepsilon }_{ij}^{\rm{tr}} $为相变应变率；$ {s}_{ij} $为偏应力张量；$ {{{A}^{\prime}}}_{ij} $为广义力$ {A}_{ij} $的偏量；$ {Z}_{1} $为材料参数；$ F_{y}^{\rm{tr}} $为相变屈服面方程；$ n $为马氏体体积分数；$ \varepsilon _{\rm{eq}}^{\rm{tr}} $为等效相变应变；$ {\varepsilon }_{m} $为单向加载条件下的最大相变应变值；$ \sigma _{s}^{\rm{tr}} $为相变的起始应力；$ \sigma _{f}^{\rm{tr}} $为相变的结束应力；$ {k}_{1} $、$ {k}_{2} $为广义力$ {A}_{ij} $演化相关的材料参数。

塑性行为内变量演化必须满足式（6），类似地构造与塑性行为相关的广义力$ {B}\text{=}-\rho \dfrac{\partial \psi }{\partial \mathbf{\eta }} $与势函数$ \varOmega \text{=}\varOmega ({\sigma }_{ij},{B}) $，推导得到NiTi合金的塑性演化规律如下：

$\qquad \dot{\varepsilon }_{ij}^{\rm{p}}={\dot{\lambda }}_{2}\frac{\partial \varOmega }{\partial {\sigma }_{ij}}\text=\frac{3}{2}{\left(\frac{F_{y}^{\rm{p}}}{{Z}_{2}}\right)}^{{{n}_{2}}}\frac{{s}_{ij}-{{{B}^{\prime}}}_{ij}}{\sigma _{{}_{\rm{eq}}}^{\rm{p}}} $

(12)

$\qquad F_{y}^{\rm{p}}={\sigma }_{\rm{eq}} - \sigma _{y}^{\rm{p}}(q)=\sqrt{\frac{3}{2}({s}_{ij} - {{{B}^{\prime}}}_{ij})\colon ({s}_{ij} - {{{B}^{\prime}}}_{ij})} - \sigma _{y}^{0} - R $

(13)

$\qquad \dot{R}=m({R}_{1}-R)\dot{\varepsilon }_{\rm{eq}}^{\rm{p}} \text{，} \dot{\varepsilon }_{\rm{eq}}^{\rm{p}}\text=\sqrt{\frac{2}{3}{\dot{\varepsilon }}^{\rm{p}}\colon {\dot{\varepsilon }}^{\rm{p}}} $

(14)

$ \qquad {\dot{B}}_{ij}={k}_{3}\dot{\varepsilon }_{ij}^{\rm{p}}-{k}_{4}\dot{\varepsilon }_{\rm{eq}}^{\rm{p}}{B}_{ij} $

(15)

式中：$ \dot{\varepsilon }_{ij}^{\rm{p}} $为塑性应变率；$ {s}_{ij} $为偏应力张量；$ {{{B}^{\prime}}}_{ij} $为广义力$ {B}_{ij} $的偏量；$ {Z}_{2} $为材料参数；$ F_{y}^{\rm{p}} $为塑性屈服面方程；$ \sigma _{y}^{0} $为初始塑性屈服应力；$ R $为屈服面演化半径；$ \dot{\varepsilon }_{\rm{eq}}^{\rm{p}} $为等效塑性应变增量；$ m $、$ {R}_{1} $为屈服面演化半径演化相关的材料参数，$ {k}_{3} $、$ {k}_{4} $为广义力$ {B}_{ij} $演化相关的材料参数。

2 动态本构参数敏感度分析 2.1 本构参数识别模型

基于不可逆热力学理论框架确定的NiTi合金动态本构模型，总共包含22个待识别参数。根据该本构模型的特点将其分为4个阶段分别进行优化识别，这4个阶段的本构模型待定参数可采用以下4个向量进行表示：

奥氏体弹性变形阶段：

$\qquad {X}_{A}\text={[{{E}_{A}}]}^ {\mathrm{T}}$

(16)

马氏体相变阶段：

$ \quad X_{\rm{tr}} = [\sigma_{s0}^{\rm{tr}},\sigma_{f0}^{\rm{tr}},\varepsilon_m,k_1,k_2,Z_1,n_1,C_1,C_2,m_1,m_2]\mathrm{^{\mathrm{\mathit{\mathrm{T}}}}} $

(17)

马氏体弹性阶段：

$ \qquad X_M=[E_M]\mathrm{^T} $

(18)

塑性屈服阶段：

$ \qquad X_p=[\sigma_{y0}^{\rm{p}},k_3,k_4,Z_2,n_2,R_1,m,C_3,m_3]^{\mathrm{T}} $

(19)

NiTi合金动态本构参数识别的目的是分别寻找这4个阶段的最优解$ {X}^{*} $，使模拟结果与试验测试结果之间的误差为最小，即寻找一组本构参数，使如下目标函数成立：

$ \qquad \left\{ \begin{array}{l} obj\colon \min W(\boldsymbol{X})=\sqrt{\sum\limits_{i\text{=1}}^{n}{({{\sigma }_{i}}({\varepsilon _{i}^{\rm{p}}})-{{\tilde{\sigma }}_{i}}({\varepsilon _{i}^{\rm{p}}}))}^{2}}\\ \boldsymbol{X}=[{{x}_{1}},{{x}_{2}},{{x}_{3}},{{x}_{4}},{{x}_{5}},{{x}_{6}},{{x}_{7}},{{x}_{8}},{{x}_{9}},{{x}_{10}},{{x}_{11}},\\{{x}_{12}}, {{x}_{13}},{{x}_{14}},{{x}_{15}},{{x}_{16}},{{x}_{17}},{{x}_{18}}]^{\text{T}}\\ =[{{E}_{A}},{{E}_{M}},{\sigma _{s0}^{\rm{tr}}},{\sigma _{f0}^{\rm{tr}}}, {\sigma _{y0}^{\rm{p}}},{{Z}_{1}},{{Z}_{2}},{{n}_{1}},{{n}_{2}},{{c}_{1}},\\ {{c}_{2}},{{c}_{3}},m,{{R}_{1}},{{k}_{1}},{{k}_{2}},{{k}_{3}},{{k}_{4}}]^{\text{T}}\\ s.t.\colon {\tilde{x}}_{i\min }\leqslant {x}_{i}\leqslant {\tilde{x}}_{i\max } \end{array}\right. $

(20)

式中：$ W({X}) $为目标函数；$ {X} $为离散设计变量；$ \tilde{\sigma}_i(\varepsilon_i^{\mathrm{p}}) $为塑性应变$ \varepsilon _{i}^{\rm{p}} $处所对应的流变应力实验测量值，$ \sigma_i(\varepsilon_i^{\mathrm{p}}) $为通过动态本构模型计算得到的塑性应变$ \varepsilon _{i}^{\rm{p}} $处所对应的流变应力模拟计算值。此时目标函数$ W({X}) $取得极小值，相应模拟计算用的材料参数$ {\boldsymbol{X}}^{*}={[{x_{1}^{*}},{x_{2}^{*}},\cdots,{x_{18}^{*}}]}^{\text{T}} $即为NiTi合金动态本构模型参数的最优解。

2.2 本构参数敏感度分析方法

在参数优化识别过程中，若对所有参数取值区间进行精密离散势必会产生巨大的计算工作量，因此需要开展参数对计算结果的敏感度分析。通过传统的单因素敏感度分析方法通常会忽略本构参数之间的相互作用对计算结果的影响。此外，由试验结果可知，NiTi合金的本构关系为非线性的，各本构参数之间存在相关性，能够对本构关系的拟合精度产生共同影响，因此需要开展参数敏感度的整体性分析，分析所有参数对目标函数的敏感程度^[16]。

对所有参数进行整体敏感度分析之前，首先必须对全参数空间抽取样本，本文通过适用于复杂多维参数空间抽样的拉丁超立方抽样（latin hypercube sampling, LHS）方法^[17-18]进行参数抽样。LHS的主要操作步骤及其过程示意如图1所示。

（1）设有$ m $个本构参数构成的参数集$ ({x}_{1}, {\mathrm{x}}_{2},\cdots,{x}_{m}) $，将每个参数的取值范围分为等概率的$ n $个互不重叠的区间；

（2）在每个参数$ {x}_{i} $的各个区间内随机抽取一个样本代表$ x_{i}^{k} $，即每个参数$ {x}_{i} $含$ n $个样本值；

（3）将每个参数$ {x}_{i} $的$ n $个样本值均随机进行排列，共生成$ m $个随机排列；

（4）每次从$ m $个排列中顺序提取各参数的一个样本值，即形成一个样本参数集$ {X}_{k} $；依次提取$ n $次，即得$ n $个样本参数集$ \left\{{X}_{1},{X}_{2},\cdots,{X}_{n}\right\} $。

通过LHS抽样方法得到NiTi合金本构关系的材料常数样本参数集，然后通过非参数统计方法中的Spearman秩相关分析法对NiTi合金本构的材料常数进行敏感度整体性分析，其基本实现过程如下：

（1）设有样本集$ \left\{{A}_{1},{A}_{2},\cdots ,{A}_{n}\right\} $，把它们按从大到小的顺序排列，构成一个有序序列$ A $；若$ {A}_{i} $在该序列$ A $中占第$ {R}_{i} $个位置，则称$ {R}_{i} $为$ {A}_{i} $在$ \left\{{A}_{1},{A}_{2},\cdots ,{A}_{n}\right\} $中的“秩”；

（2）将$ n $个样本参数集$ \left\{{X}_{1},{X}_{2},\cdots ,{X}_{n}\right\} $代入目标函数可计算得目标函数解集$ \left\{{W}_{1},{W}_{2},\cdots ,{W}_{n}\right\} $，即：每个参数$ {x}_{i} $的样本集$ \left\{{x}_{i}{}^{1},{x}_{i}{}^{2},\cdots ,{x}_{i}{}^{n}\right\} $所对应的目标函数解集为$ \left\{{W}_{1},{W}_{2},\cdots ,{W}_{n}\right\} $；

（3）设参数$ {x}_{i}{}^{k} $在$ \left\{{x}_{i}{}^{1},{x}_{i}{}^{2},\cdots ,{x}_{i}{}^{n}\right\} $中的“秩”为$ {\alpha }_{k} $，$ {W}_{k} $在$ \left\{{W}_{1},{W}_{2},\cdots ,{W}_{n}\right\} $中的“秩”为$ {\beta }_{k} $，则参数$ {x}_{i} $的随机输入样本集$ \left\{{x}_{i}{}^{1},{x}_{i}{}^{2},\cdots ,{x}_{i}{}^{n}\right\} $与其对应的输出结果集$ \left\{{W}_{1},{W}_{2},\cdots ,{W}_{n}\right\} $的Spearman秩相关系数的绝对值即为参数$ {x}_{i} $的敏感度$ {r}_{i} $，其表达式为：

$\qquad {r}_{i}=\left| \frac{n\sum\limits_{k=1}^{n}{\alpha }_{k}{\beta }_{k}-\sum\limits_{k=1}^{n}{\alpha }_{k}\sum\limits_{k=1}^{n}{\beta }_{k}}{\sqrt{n\sum\limits_{k=1}^{n}{({{\alpha }_{k}})}^{2}-{\left(\sum\limits_{k=1}^{n}{\alpha }_{k}\right)}^{2}}\sqrt{n\sum\limits_{k=1}^{n}{({{\beta }_{k}})}^{2}-{\left(\sum\limits_{k=1}^{n}{\beta }_{k}\right)}^{2}}}\right| $

(21)

若通过上式求得的$ {r}_{i} $数值越大，则说明该本构参数对目标函数输出结果的敏感度越高。根据参数的敏感度分析结果，将所有参数的取值区间进行不同疏密程度的离散划分，上述操作能够有效减少后续参数优化识别的计算量，并保证参数识别结果的可靠性^[19]。

2.3 NiTi合金动态本构参数敏感度分析

NiTi合金动态本构参数及每个参数的LHS抽样范围如表1所示，在每个参数范围内均匀抽样20次。按照上一节Spearman 秩定义方法可以得出NiTi合金动态本构参数秩值，然后根据式（21）可以计算出NiTi合金本构参数Spearman相关系数，如表2所示。从计算结果可以确定出$ {Z}_{1} $、$ {Z}_{2} $、$ {n}_{1} $、$ {n}_{2} $、$ {c}_{1} $、、$ {c}_{3} $、$ \sigma_{s0}^{\mathrm{tr}} $、$ \sigma_{f0}^{\mathrm{tr}} $、$ \sigma_{y0}^{\mathrm{p}} $、$ {R}_{1} $、$ {E}_{A} $、$ {E}_{M} $的敏感度值相对较大，即参数对输出值的影响程度较大，在识别NiTi合金动态本构参数时需重点关注这些参数的取值范围，而$ {k}_{1} $、$ {k}_{2} $、$ {k}_{3} $、$ {k}_{4} $敏感度值相对较低，即参数在整个计算过程中对输出值的影响相对较小，成为NiTi合金动态本构次要影响因素。

根据Spearman秩相关分析法理论可知，为了保证遗传算法的准确性，需分别调整各参数的取值范围和所有参数的抽样次数，使各参数敏感度值彼此接近且均小于0.2。若某一参数的敏感度比其他参数的敏感度过于偏小，则说明该参数的求解精度过高，与其他参数的求解精度不协调，会无谓地增加计算量；若某一参数的敏感度大于0.2，则该参数的求解精度不够。

3 动态本构参数的优化识别

基于本构参数敏感度分析结果，采用遗传算法可快速实现对本文构建的NiTi合金动态本构模型22个待定参数的优化识别。遗传算法是模拟自然界生物进化机制的一种算法，通过编码组成初始群体后，对个体按照它们对环境的适应度施加一定操作，从而实现优胜劣汰的进化过程。遗传算法包括以下3个基本遗传算子：选择、交叉、变异。与传统的优化算法相比，遗传算法具有更好的全局搜索能力，其基本特征是利用群体进化，即在求解的过程中，通过种群不断优化，从而找到最优解。但基本遗传算法（simple genetic algorithm, SGA）在应用中存在诸多缺陷，主要包括：无法同时满足精度、可靠性和计算时间三方面的要求，且容易产生早熟现象，局部寻优能力较差等。

本文针对基本遗传算法的诸多缺陷进行有效改进，采用改进小生境算法、可疑峰值点判断策略和局域精确搜索技术，提高了遗传算法的全局搜索能力。其关键实现要点为：

（1）对小生境算法的改进：

建立最优个体保存策略。各小生境独立进化，每代进化完毕后保存当前最优个体。保存策略为：在每代进化开始前首先复制当前最优个体作为副本，然后让所有个体平等地参与交叉、变异运算；该代进化完毕后，若该最优个体仍保留在群体中，则不做任何处理，否则就用当前最优个体的副本替换掉当代最差个体。该方法能确保最优个体不在进化中被淘汰，还能加快算法收敛速度。

在各小生境间建立数据交换机制。每进化一代都会在各小生境间进行最优个体的交换，即用第1个小生境的最优结果替换第2个小生境的最差结果，用第2个的最优替换第3个的最差，以此类推。该方法能在保证群体多样性的前提下提高优等个体的比例，加快收敛速度。

（2）建立可疑峰值点判断策略：

每个小生境独立进化完毕后，均收敛于1个峰值点处。当小生境数目较大时，求得的所有峰值点包含全局最优点和局部最优点。在最终确定全局最优点前，对求得的所有峰值点均视为“可疑峰值点”。所有可疑峰值点中的最优点A必为全局最优点，记其目标函数值为F^*，对应的变量为$ \{{x}_{1}{}^{*},{x}_{2}{}^{*},\cdots ,{x}_{m}{}^{*}\} $。设某一可疑峰值点B的目标函数为F，其对应的变量为$ \{{x}_{1},{x}_{2},\cdots ,{x}_{m}\} $，则B点为相异于A点的全局最优点的判断条件为：

$\qquad \left\{\begin{array}{l} \left| \dfrac{F-{F}^{*}}{{F}^{*}}\right| <\alpha \\ {d}_{i}=\sqrt{\displaystyle\sum\limits_{k=1}^{n}{({{x}_{i}}-{x_{i}^{*}})}^{2}} >{\beta }_{i}(i=1,2,\cdots ,m) \end{array}\right. $

(22)

式中：$ \alpha $为0～1之间的常数，指最优点附近的目标函数值域相对搜索范围；$ {\beta }_{i} $为判断两个体是否“相邻”的常数；$ {d}_{i} $为两个体的海明距离，若$ {d}_{i}\leqslant {\beta }_{i} $，则说明A、B的相隔距离较近，经精确搜索后必收敛于同一个峰值点。式（4）的第1个不等式用于确定该点为全局最优点，第2个不等式用于确定B点为相异于A点的全局最优点。

（3）对所有相异的全局最优点作局域精确搜索：

设某个全局最优点对应的变量集为$ \{{x}_{1}{}^{*}, {x}_{2}{}^{*},\cdots ,{x}_{m}{}^{*}\} $，其各个变量的取值范围为$ {x}_{il}\leqslant {x}_{i}\leqslant {x}_{iu} $，则将各个变量的上下限边界条件改为：

$\quad {x}_{i}\in [{x}_{il},{x}_{iu}]\cap [x_{i}^{*}-{\gamma }_{i},x_{i}^{*}+{\gamma }_{i}]\quad (i=1,2,\cdots ,m) $

(23)

式中：$ {\gamma }_{i} $为变量$ {x}_{i} $附近的搜索邻域。

然后以各变量新的上下限进行遗传进化操作，直到达到预先设定的最大进化代数为止，这就进一步提高了最优解的精度。改进遗传算法的程序流程图如图2所示。

利用改进遗传算法对本构参数进行优化识别时，可将基于不可逆热力学理论框架构建的NiTi合金动态本构模型分为4个阶段（奥氏体弹性变形阶段，相变阶段，马氏体弹性变形阶段和塑性流动阶段）分别建立目标函数进行本构参数识别求解。在求解时设置交叉概率为0.5，变异概率为0.09，终止代数为500，组成各变量的二进制码有9位，最终识别得到NiTi合金动态本构参数的优化结果如表3所示。

4 本构模型验证与结果分析

将表3中的本构参数优化结果代入NiTi合金相变与塑性统一的动态本构模型，考虑到复杂的本构演化方程和硬化法则，采用全隐式应力积分方法求解非弹性应变增量比较困难，本文采用一种较为简单的半隐式应力积分方法，即半隐式迭代法求解非弹性应变增量，在下一个增量步对其进行显式更新。主要步骤如下：

（1）输入应力初始值、应变增量、状态变量初始值、时间步增量以及初始弹性刚度阵；

（2）假设应变增量全部为弹性应变，按弹性方法计算应力并得到等效应力，进行屈服条件判断，若发生屈服则分别进入相应的相变或塑性迭代步骤；若不满足屈服条件则进入相应的弹性阶段（即母相弹性或马氏体弹性段）；

（3）非弹性段（即相变或塑性段）的统一迭代过程，通过计算下一增量步的非弹性应变增量和新的等效非弹性应变增量值，进行迭代收敛性判断进而更新所有的状态变量；

（4）更新非弹性应变及其增量、弹性应变及其增量。广义力和总应力等，进入下一加载步。

具体的NiTi合金本构模型数值求解流程示意图如图3所示。

根据图3数值求解流程与表3中NiTi合金动态本构参数的优化结果，可以得到NiTi合金在500/s、1 500/s、2 100/s、3 000/s应变率下的应力-应变数值计算结果，与试验数据^[15]对比如图4所示。

由对比结果可知，不同应变率条件下NiTi合金动态本构模型的数值计算结果与实验数据吻合良好。在整个冲击模拟过程中材料先后出现了4个变形阶段：奥氏体弹性变形、马氏体相变、马氏体弹性变形及塑性屈服阶段，所构建的相变与塑性统一本构模型对这4个阶段均能够较好地拟合。同时也说明利用Spearman秩相关分析法获得NiTi合金动态本构参数的敏感度，并通过改进遗传算法得到NiTi合金动态本构参数的合理性。本文中的参数敏感度整体性分析方法和改进遗传算法对其他工程材料本构参数的优化识别同样具有重要参考价值。

此外，在较高的应变率条件下，塑性屈服极限的实验结果与数值模拟结果出现一定差异，推断该现象由以下原因引起：（1）在高应变加载条件下，材料发生绝热变形，其内部因塑性功转化的热量会出现局部温升，温升会导致流变应力出现一定程度下降，即温度软化效应^[20-21]；（2）可能是动态冲击实验过程中的试件摩擦、应变片测试信号误差等因素综合作用导致的。上述推断将在后续的研究工作中开展进一步的分析与验证。

5 结　论

（1）本文以不可逆热力学理论为基础，构建了NiTi合金的相变与塑性统一动态本构模型。采用拉丁超立方抽样方法在整个本构参数空间中抽样，并用非参数统计方法中的 Spearman 秩相关分析法对本构参数随机输入样本集与其对应的目标函数输出结果集作相关性分析，建立了用 Spearman 秩相关系数等效求解参数敏感度的表达式，进而实现了参数敏感度的整体性分析，根据 Spearman相关系数大小来反映NiTi合金动态本构参数对目标函数的影响程度, 从而确定了各个参数的主次关系，为多因素系统分析提供了一种行之有效的手段。该方法可以克服单参数分析法的缺点，是一种更接近实际的方法。

（2）基于Spearman秩相关分析法得到的NiTi合金动态本构参数敏感度值，确定本构参数的识别范围。利用改进遗传算法对NiTi合金动态本构参数进行识别，减少本构模型参数识别的工作量，并且实现了快速、精确、可靠地搜索到本构参数的最优解。

（3）将本构参数的优化识别结果代入NiTi合金动态本构模型，并采用半隐式应力积分方法更新非弹性应变增量，得到材料在不同应变率下的应力−应变数值计算结果与试验数据吻合良好，验证了本文的参数整体敏感度分析与优化识别方法的准确性，所提方法对其他工程材料本构参数的优化识别同样具有重要参考价值。