土木与环境工程学报  2020, Vol. 42 Issue (6): 112-118   doi: 10.11835/j.issn.2096-6717.2020.080   PDF    
形状记忆合金数值模型的不确定性分析
李金朋 1a, 陈城 2, 侯和涛 1     
1a. 山东大学 土建与水利学院, 济南 250061;
1b. 山东大学 山东省绿色建筑智能建造工程技术研究中心, 济南 250061;
1c. 山东大学 山东省工业技术研究院工业化建筑智能建造协同创新中心, 济南 250061;
2. 旧金山州立大学 工程学院, 旧金山 94132
收稿日期：2020-03-20
基金项目：国家自然科学基金（51878390）
作者简介：李金朋(1994-), 男, 主要从事结构抗震不确定性研究, E-mail:201714576@mail.sdu.edu.cn;
侯和涛(通信作者), 男, 教授, 博士, 博士生导师, E-mail:houhetao@sdu.edu.cn.
摘要：形状记忆合金（shape memory alloy，简称为SMA）具有"超弹性"，即在受到应力而发生较大变形并卸载后，可以恢复原始形状，并在这个过程中耗散能量，在建筑抗震和桥梁振动控制中具有广阔的应用前景。SMA的模型参数通常由优化方法来确定，然后用于装有SMA装置的结构地震时程响应分析中。利用Metropolis-Hasting算法（简称为MH算法）中的改进算法DRAM方法（延迟拒绝及自适应采样），基于经过"预拉伸"和热处理的SMA棒材循环拉伸试验结果，对SMA改进的Graesser & Cozzarelli模型参数进行采样，从SMA的本构模型参数和耗能能力两个方面分析了SMA材料的不确定性。建立了各参数的后验分布，并得到了参数两两之间的相关性，结果可用于概率模型的建立及基础模型数学形式的研究。研究表明，在累积概率密度为15%时，材料的能量耗散能力相对误差高达20%；累积概率密度为85%时，相对误差为10%。
关键词：形状记忆合金    不确定性分析    马尔可夫链蒙特卡罗方法    概率建模    
Probabilistic analysis of shape memory alloy modeling
Li Jinpeng 1a, Chen Cheng 2, Hou Hetao 1     
1a. School of Civil Engineering, Shandong University, Jinan 250061, P. R. China;
1b. Shandong Research Center for Green Building and Intelligent Construction, Shandong University, Jinan 250061, P. R. China;
1c. Collaborative Innovation Center for Industrialized Building and Intelligent Construction, Shandong Industrial Technology Research Institute, Shandong University, Jinan 250061, P. R. China;
2. School of Engineering, San Francisco State University, San Francisco, CA 94132, U.S.A
Received: 2020-03-20
Foundation item: National Natural Science Foundation of China (No. 51878390)
Author brief: Li Jinpeng (1994-), main research interest:uncertainty analysis, E-mail: 201714576@mail.sdu.edu.cn;
Hou Hetao(corresponding author), professor, PhD, doctorial supervisor, E-mail:houhetao@sdu.edu.cn.
Abstract: Shape memory alloy(SMA) has "super elasticity", that is, it can recover original shape after deformation and unloading due to stress, and dissipate energy in this process. It has broad application prospect in seismic control of buildings and bridge vibration. The model parameters of SMA are often determined through optimization and treated as deterministic for dynamic analysis of structures with SMA based devices. In this study, the modified Metropolis-Hasting algorithm-DRAM algorithm, which is a combination of delay rejection and adaptive sampling, is utilized to characterize the uncertainties in modified Graesser & Cozzarelli SMA model parameters. A series of SMA bars with the same geometric size and heat treatment were tested under cyclic loads. The Markov Chain Monte Carlo (MCMC) method is applied to analyze the uncertainties of SMA in terms of model parameters and energy dissipation capacity. The analysis provide insight into the underlying mathematical form of a model, suggest simplifications or modifications and begin to indicate the relative significance of individual parameters, based on a limited set of experimental data. Besides, research shows thatthe energy dissipation of the SMA bar could have up to a relative error of 20% and 10% corresponding to the CDF of 15% and 85%.
Keywords: shape-memory alloy    uncertainty analysis    Markov Chain Monte Carlo    probabilistic modeling    

形状记忆合金(SMA)具备形状记忆，这使其在经历较大幅度的变形后，可通过加热或者卸载恢复原本形状[1]。奥氏体相下的SMA受到应力而发生变形，并在卸载后恢复原始形状的行为称为“超弹性”(或“伪弹性”)。SMA的化学成分以及生产中的冶金处理过程对上述性质具有显著影响，而材料固有的不确定性导致这两个因素不易被精确控制。SMA的单轴应力-应变响应通常为典型的旗帜型滞回曲线，并且具有良好的自定心能力、能量耗散能力和循环可重复性[1-2]。这使其在结构抗震装置，如振动控制装置[3]、多跨桥梁的限位装置[4]等领域得到广泛应用，其数值模型也得到了深入的研究。现有研究中，当进行基于SMA装置的地震响应模拟时，SMA模型参数通常通过优化算法来确定[5-6]。然而，数学模型的简化、不可避免的实验误差和诸多其他因素均可能导致材料或者结构的数值模型产生不确定性，进而导致模拟结果具有局限性并可能失真[7]。

笔者提出对模型参数进行概率建模的方法，基于SMA棒材的实验数据，采用改进的Graesser和Cozzarelli模型与MCMC算法的组合来分析模型本身固有的不确定性，将模型参数视为随机变量，采用Metropolis-Hastings算法来生成样本参数集，揭示了参数的概率特性和参数之间的潜在相关性，并从模型参数的角度研究了SMA模型中固有的不确定性及其对材料能量耗散能力预测值的影响。

1 SMA的数值模型与材性试验
1.1 改进的Graesser & Cozzarelli模型

基于Ozdemir的一维滞回模型[8]，Graesser和Cozzarelli提出了描述SMA的宏观应力-应变关系的Graesser & Cozzarelli模型[9]，其数学表达式为

$ {\dot \sigma = E\left[ {\dot \varepsilon - |\dot \varepsilon |{{\left( {\frac{{\sigma - \beta }}{Y}} \right)}^n}} \right]} $ (1)
$ {\dot \beta = \alpha E\left\{ {{\varepsilon _{{\rm{in }}}} + {f_{\rm{T}}}|\varepsilon {|^c}{\mathop{\rm erf}\nolimits} (a\varepsilon )[u( - \varepsilon \dot \varepsilon )]} \right\}} $ (2)

式中：σ为一维应力；β为一维背应力；ε为一维的应变；E为材料的弹性模量；Y为材料的屈服应力；参数n控制SMA由弹性阶段向塑性阶段过渡时滞回曲线的“尖锐”程度，取值范围设计为任意正奇数；参数α控制塑性阶段滞回曲线斜率，由α=Ey/(E-Ey)得到，其中Ey为屈服后的弹性模量；fT是控制滞回曲线类型和大小的参数，当fT=0时，模型代表纯马氏体状态下的SMA，当fT>0时，模型可模拟SMA的“超弹性”行为；a是控制卸载过程中材料弹性恢复量的参数；c是控制卸载过程中应力平台段斜率的参数；(·)表示对时间的一阶导数；u(x)和erf(x)分别为单位阶跃函数和误差函数，数学表达式为

$ \begin{array}{*{20}{l}} {u(x) = \left\{ {\begin{array}{*{20}{c}} { + 1}&{x \ge 0}\\ 0&{x < 0} \end{array}} \right.}\\ {{\mathop{\rm erf}\nolimits} (x) = \frac{2}{{\sqrt {\rm{ \mathsf{ π} }} }}\int_0^{\rm{ \mathsf{ π} }} {{{\rm{e}}^{ - {t^2}}}} {\rm{d}}t} \end{array} $ (3)

但是，该模型无法模拟SMA由奥氏体相向马氏体相转化完成后出现的“硬化”现象，即材料的弹性模量突然增大，称为SMA的马氏体硬化特性。为了描述这一特性，Qian等[10]提出了改进的Graesser & Cozzarelli模型，其数学表达式为

$ {\dot \sigma = E\left[ {\dot \varepsilon - |\dot \varepsilon |\left| {\frac{{\sigma - \beta }}{Y}} \right|n - 1\left( {\frac{{\sigma - \beta }}{Y}} \right)} \right]} $ (4)
$ \begin{array}{l} \;\;\;\dot \beta = \alpha E{\varepsilon _{{\rm{in }}}} + {f_T}|\varepsilon {|^c}{\mathop{\rm erf}\nolimits} (a\varepsilon )[u( - \varepsilon \dot \varepsilon )] + \\ {f_{\rm{M}}}{\left[ {\varepsilon - {\varepsilon _{{\rm{Mf}}}}{\mathop{\rm sgn}} (\varepsilon )} \right]^m}[u(\varepsilon \dot \varepsilon )]\left[ {u\left( {|\varepsilon | - {\varepsilon _{{\rm{Mf}}}}} \right)} \right] \end{array} $ (5)

改进后，相较于式(1)，式(4)中参数n的取值范围变为任意正实数，更方便该参数的不确定性分析。另外，新引入的fM和m为控制马氏体硬化阶段曲线的参数；εMf为SMA转化为全马氏体时的应变；sgn(x)为符号函数，数学表达式为

$ {\mathop{\rm sgn}} (x) = \left\{ {\begin{array}{*{20}{l}} { + 1}&{x > 0}\\ 0&{x = 0}\\ { - 1}&{x < 0} \end{array}} \right. $ (6)

相较于式(2)，式(5)添加了描述马氏体硬化的表达式。在SMA进入马氏体硬化阶段前，这部分值为0，此时式(5)与式(2)完全相同。

1.2 SMA棒材的循环拉伸试验

采用粒子群优化(PSO)方法[11]得到使模型具备良好拟合效果的参数值作为初始参数组，通过试验得到一组SMA棒材的循环拉伸测试的实测数据，采用PSO方法得到模型的确定性参数值。试验试件的原材料为直径12 mm的镍-钛SMA棒材，其化学成分如表 1所示。棒材加工成“狗骨”形状的试样，削弱部分直径6 mm，便于循环拉伸测试，如图 1(a)所示。热处理的温度为400 ℃，持续时间为15 min。退火后，观察到试样的颜色从银色变为金色，如图 1(b)、(c)所示。试件的加载方式如图 1(d)所示。在热处理之前，先对试件进行峰值应变为7%的准静态拉伸并卸载处理，使材料内部晶体结构重新排列，有助于其性能的发挥[12]；处理结束后，试件产生了3.8%的残余应变。

表 1 SMA试件原材料的化学成分 Table 1 Chemical composition of the testing SMA bar

图 1 SMA循环拉伸试验 Fig. 1 Cyclic tests of SMA bar

将热处理后的试件进行循环拉伸试验，其加载制度如图 2(a)所示。图 2(b)为循环拉伸试验得到的滞回曲线，其数据用于后续的不确定性分析。在进行不确定性分析之前，基于峰值应变为0.08的两条滞回曲线，采用PSO方法得到使模型具备良好拟合效果的参数值作为初始参数组，有助于提高马尔可夫链蒙特卡洛(MCMC)模拟时马尔可夫链的收敛速度。由PSO方法得到的初始参数具体数值如表 2所示。从图 2(b)可以看到，PSO优化方法得到的参数虽然使得模型具有良好的拟合效果，但模拟结果与试验数据仍旧存在偏差，因此，有必要研究模型中存在的固有不确定性。

图 2 加载制度和试验及优化模型模拟结果示意图 Fig. 2 Loading protocol and comparison between experiment results and PSO methods

表 2 基于循环拉伸试验数据的PSO优化参数值 Table 2 Model parameter values from cyclic tests of SMA bars using PSO method

2 MCMC不确定性分析
2.1 MCMC方法

马尔可夫链蒙特卡洛(MCMC)是一种通过建立一条按照提议分布π(x)平滑分布的马尔科夫链来获取所需样本点的方法。MCMC方法通过沿着马尔可夫链计算$1/n\sum\limits_{i = 1}^n f \left( {{X_i}} \right)$的值，当该值趋于稳定后，将该值作为给定函数f相对于分布π(x)的期望Eπf的估计(此时认为马尔可夫链收敛)，进而得到平滑服从所给提议分布的马尔可夫链的样本点。其中，$1/n\sum\limits_{i = 1}^n f \left( {{X_i}} \right)$称为MCMC算子。通常，MCMC采样是渐进无偏且服从正态分布的[13]。

2.1.1 DRAM方法

采用MCMC方法中两个重要方法的结合，即延迟拒绝法(DR法，Delaying Rejection)和自适应采样(AM法，Adaptive Metropolis Samplers)，简称为DRAM方法[13-14]。DR法通过适当调整标准Metropolis-Hastings(MH)算法中马尔可夫链每一步中的提议分布来提高MCMC算子的效率。AM法则是基于马尔可夫链的历史来调整提议分布。当提议分布的方差非常小时，标准MH算法倾向于以小步长“遍历”目标分布，而无法有效地探索状态空间，且样本点分布偏移样本空间中心，产生偏差。Haario等[14]的研究证明了AM方法能够解决“探索范围未覆盖整个样本空间”的问题，DR方法能够解决“样本集中位置偏移样本空间中心”的问题，而DR和AM的组合，即DRAM方法可以同时解决这两个问题。

2.1.2 似然函数

似然函数定义为给定的一组参数值下，模型模拟结果与实验数据一致的概率，也可以将其视为模型预测和实验测量之间的误差概率。ss函数是似然函数的一部分，用于描述归一化后的误差[15]。似然函数及ss函数公式为

$ {L_j} = \exp \left[ { - 0.5 \times \left( {\frac{{s{s_j} - s{s_{j - 1}}}}{{{\sigma ^2}}}} \right)} \right] $ (7)
$ s{s_j} = \frac{{\sum\limits_{i = 1}^n {{{\left( {{f_j}\left( {{t_i}} \right) - {F^*}\left( {{t_i}} \right)} \right)}^2}} }}{{\sum\limits_{i = 1}^n {{{\left( {{F^*}\left( {{t_i}} \right)} \right)}^2}} }} $ (8)

式中：Lj为第j步的似然函数值；ssj为第j步的归一化误差值；n为实验数据的数据点个数；σ2为设定的容许误差，用以限制ss的大小；F*(ti)为在时间点ti的实验测量结果；fj(ti)为第j步时在时间点ti的模型输出结果。

DRAM模拟方法的具体流程如图 3所示，采用MATLAB编程并模拟。模拟次数(Nsimu)设置为2×106次，为了研究合理范围内的不确定性，容许误差σ2设为0.022。

图 3 DRAM模拟方法的流程示意图 注: Θ为参数向量(E,Y,a,c,fT,αfM, εMf,m,n)T, Θj为第j组的参数向量;γ为机器生成的[0,1]间的随机数,在每次循环中更新。 Fig. 3 Flow chart of simulation process

2.1.3 参数的先验概率

在模拟开始之前，首先要建立参数的先验概率，包括参数的范围及其在该范围内的分布。研究表明，先验分布并不是MCMC模拟得到参数所收敛的后验分布的决定性因素，而是影响收敛速度的关键因素[16]。因此，模拟中仅根据数学或物理要求对参数的范围进行合理规定，不指定其分布方式。参数的先验概率设置如表 3所示。

表 3 参数的先验概率设置 Table 3 Prior of model parameters

2.2 MCMC分析结果
2.2.1 模型参数

基于峰值应变为8%的加载循环试验数据得到的参数不确定性分析结果如表 4、图 4所示。表 4列出了10个参数的概率特性，包括均值、方差和偏度。图 4为所有参数的频率分布直方图，图 4中大多数参数的频率分布呈现出单峰结构，表明对参数采样的马尔可夫链收敛服从该参数的后验概率分布。图 4显示，除了参数E、参数Y和参数εMf外，其余参数后验分布的均值均大于PSO优化方法得到的初始参数值。这意味着PSO的优化结果是否最优有待商榷，原因可能在于其算法中粒子游走出现某一方向的偏差而错失最优路径。此外，绝大部分参数的分布峰值位置与初始值基本吻合，表明优化方法得到的参数值通常具有普遍性；然而，参数c的分布中，峰值出现的位置距离初始值有较大偏差，表明PSO优化方案应用在SMA模型中可能具有片面性以及偏差性，即优化方法提供的参数值可能并不是该参数最普遍采用的值。这表明SMA具有本身固有的不确定性，利用MCMC方法研究SMA不确定性是必要的。

表 4 参数分布的均值、标准差和偏度 Table 4 Moments of model parameters

图 4 参数频率分布直方图 注:红色实线为参数的初始值，蓝色虚线为参数分布均值。 Fig. 4 Histogram of model parameters

图 5显示了参数α、εMf、m和fM之间的成对相关性。两个参数的相关性越强，则对应椭圆的离心率越大；椭圆长轴方向为左下至右上时表示参数之间存在正相关关系，反之呈负相关关系。由图 5可见，参数m和α、a和c、m和fM两两之间存在较强正相关性；fT和α之间存在较强负相关性；fT和a、c和fM等之间不存在线性相关性或相关性较弱。参数之间的相关性可以为数学模型的研究提供参考。

图 5 参数相关性示意图 Fig. 5 Correlation between the simulated parameters

2.2.2 能量耗散

为了更好地说明不确定性研究的必要性，研究通过模型参数的概率分布建立材料耗能能力的概率特征，其概率密度示意图如图 6(a)所示。与实验数据得到的结果相比，在累积概率密度为15%时，材料的能量耗散能力相对误差高达20%；累积概率密度为85%时，相对误差为10%。图 6(b)显示了对应能量耗散的累积概率密度值分别为0.15、0.5、0.64和0.85时模型的拟合效果。结果显示，此时模型具有非常好的拟合效果，但没有PSO优化结果好。这表明更新模型参数不会消除或补偿模型模拟的偏差，但会将偏差控制在可接受的范围内。

图 6 基于能量耗散能力预测的不确定性分析结果示意图 Fig. 6 The results based on the energy dissipation analysis

能量耗散通常通过等效粘性阻尼(EVD)ξeq来衡量，ξeq是一个与尺寸无关的指数，表示为

$ {\xi _{{\rm{eq}}}} = \frac{{{W_{\rm{D}}}}}{{2{\rm{ \mathsf{ π} }}{W_{\rm{E}}}}} $ (7)

式中：WD为每个加载周期SMA的能量耗散值；WE为相应线性系统的应变能。EVD在很大程度上取决于滞回环的形状，较大的滞回环会得到较大的EVD值，滞回曲线越饱满，EVD值越大。由图 7可知，随着峰值应变的增加，EVD略微增加，当峰值应变达到8%时，EVD略有下降。这表明当加载峰值应变达到近6%时，材料达到其最佳能量耗散能力。这一结论对基于SMA的地震控制装置的设计具有借鉴意义。

图 7 基于不同应变峰值的EVD不确定性分析结果示意图 Fig. 7 Uncertainty analysis of EVD for different peak-strain

3 结论

1) 基于形状记忆合金棒材循环拉伸试验数据的DRAM算法采样得到的马尔可夫链体现出模型参数的概率特性。样本的分布特征(均值、方差等)体现出优化方法可能存在偏差，部分参数之间存在线性相关性，在进行数值模型研究时应予以重视。

2) 数值模型的不确定性也体现在模型的耗能预测上，在累积概率密度为15%时，材料的能量耗散能力相对误差高达20%；累积概率密度为85%时，相对误差为10%。加载应变峰值对材料的耗能性能有明显影响，等效粘滞阻尼分布显示，加载峰值应变为6%时，材料耗能性能较其他对比组更好。

参考文献
[1]
DESROCHES R, MCCORMICK J, DELEMONT M. Cyclic properties of superelastic shape memory alloy wires and bars[J]. Journal of Structural Engineering-ASCE, 2004, 130(1): 38-46. DOI:10.1061/(ASCE)0733-9445(2004)130:1(38)
[2]
DOLCE M, CARDONE D. Mechanical behaviour of shape memory alloys for seismic applications 2. Austenite NiTi wires subjected to tension[J]. International Journal of Mechanical Sciences, 2001, 43(11): 2657-2677. DOI:10.1016/S0020-7403(01)00050-9
[3]
XU L H, XIE X S, LI Z X. A self-centering brace with superior energy dissipation capability:development and experimental study[J]. Smart Materials and Structures, 2018, 27(9): 095017. DOI:10.1088/1361-665X/aad5b0
[4]
VARELA S, SAIIDI M S. A bridge column with superelastic NiTi SMA and replaceable rubber hinge for earthquake damage mitigation[J]. Smart Materials and Structures, 2016, 25(7): 075012.
[5]
黄宙, 李宏男, 付兴. 自复位放大位移型SMA阻尼器优化设计方法研究[J]. 工程力学, 2019, 36(6): 202-210.
HUANG Z, LI H N, FU X. Optimum design of a re-centering deformation-amplified SMA damper[J]. Engineering Mechanics, 2019, 36(6): 202-210. (in Chinese)
[6]
QIU C X, ZHU S Y. Performance-based seismic design of self-centering steel frames with SMA-based braces[J]. Engineering Structures, 2017, 130: 67-82. DOI:10.1016/j.engstruct.2016.09.051
[7]
CAICEDO J M, JIANG Z S, BAXTER S C. Including uncertainty in modeling the dynamic response of a large-scale 200 kN magneto-rheological dampe[J]. ASCE-ASME Journal of Risk and Uncertainty in Engineering Systems, Part A:Civil Engineering, 2017, 3(2): G4016002.
[8]
OZDEMIR H. Nonlinear transient dynamic analysis of yieldling structures[D]. Berkeley, California: the University of California, 1976.
[9]
GRAESSER E J, COZZARELLI F A. Shape-memory alloys as new materials for aseismic isolation[J]. Journal of Engineering Mechanics, 1991, 117(11): 2590-2608. DOI:10.1061/(ASCE)0733-9399(1991)117:11(2590)
[10]
QIAN H, LI H N, SONG G. A constitutive model of shape memory alloys with consideration of martensitic hardening effect[C]//11th Biennial ASCE Aerospace Division International Conference on Engineering, Science, Construction, and Operations in Challenging Environments. 2008.
[11]
毛湘云, 徐冰峰, 孟繁艺. PSO-SVM与BP神经网络组合预测供水系统余氯的方法[J]. 土木与环境工程学报(中英文), 2019(4): 159-164.
MAO X Y, XU B F, MENG F Y. Prediction of residual chlorine in water supply system by PSO-SVM and BP neural network combined model[J]. Journal of Civil and Environmental Engineering, 2019(4): 159-164. (in Chinese)
[12]
GENCTURK B, ARAKI Y, KUSAMA T, et al. Loading rate and temperature dependency of superelastic Cu-Al-Mn alloys[J]. Construction and Building Materials, 2014, 53: 555-560. DOI:10.1016/j.conbuildmat.2013.12.002
[13]
CUI T G, WARD N D, EVESON S, et al. Pragmatic approach to calibrating distributed parameter groundwater models from pumping test data using adaptive delayed acceptance MCMC[J]. Journal of Hydrologic Engineering, 2016, 21(2): 06015011. DOI:10.1061/(ASCE)HE.1943-5584.0001267
[14]
HAARIO H, SAKSMAN E, TAMMINEN J. An adaptive metropolis algorithm[J]. Bernoulli, 2001, 7(2): 223. DOI:10.2307/3318737
[15]
PRAJAPAT K, RAY-CHAUDHURI S. Prediction error variances in Bayesian model updating employing data sensitivity[J]. Journal of Engineering Mechanics, 2016, 142(12): 04016096. DOI:10.1061/(ASCE)EM.1943-7889.0001158
[16]
CHING J, CHEN Y C. Transitional Markov Chain Monte Carlo method for Bayesian model updating, model class selection, and model averaging[J]. Journal of Engineering Mechanics, 2007, 133(7): 816-832. DOI:10.1061/(ASCE)0733-9399(2007)133:7(816)
X