赵澎涛，裘进浩，于慧臣

（1.北京航空材料研究院航空材料检测与评价北京市重点实验室，先进高温结构材料重点实验室，北京100095；2.南京航空航天大学机械结构力学及控制国家重点实验室，南京210016）

0 引 言

形状记忆合金（SMA）同时具备感知和驱动功能，可以感知外界的温度变化，并能够将热能转化为机械能，其以优异的形状记忆效应、伪弹性［1-2］、生物相容性和阻尼特性，广泛应用于航空航天［3-4］、生物医药［5］、土木工程［6］等领域。但由于SMA 复杂的热力学耦合关系，其本构关系不能用通常的弹塑性模型描述，因此SMA 本构行为的描述是目前研究的热点和难点。为了实现SMA更广泛的工程应用，作者以SMA在三维自适应发动机进气道变体结构中的应用为研究背景，利用SMA所特有的形状记忆效应，使用温度控制的方式实现了鼓包结构的变形。通过建立SMA三维本构关系模型，采用回归映射应力更新算法进行有限元实现，并对SMA三维鼓包结构在热变形过程中的热力学特性响应进行有限元分析，得到了一定温度条件下该结构的变形参数，为SMA在发动机进气道中的应用提供支持。

1 SMA的本构关系模型

采用Lagoudas关于SMA的三维唯象理论模型［7-8］，即SMA在变形过程中产生的应变主要由弹性应变、热膨胀应变和相变引起的应变等三部分构成，假设SMA的Gibbs自由能（G）为

式中：σ，εt，ξ，T 以及T0分别为真实应力张量、相变应变张量、马氏体体积分数、温度以及参考温度；S，α，ρ，c，s0和u0分别为等效柔度张量、等效热膨胀系数、密度、等效比热、参考温度下的等效比熵和参考温度下的等效内能；f（ξ）为相变硬化函数，表征了马氏体相与母相或者马氏体变体之间发生相互转化而产生的应变能。

以上各等效材料特性参数可由材料中所含的马氏体体积分数以及纯晶体相的特性参数计算得到，其具体的表达式如式（2）。式中物理量右上角处的A与M分别表示为奥氏体相和马氏体相：

结合热力学相关理论，SMA的本构关系及约束条件可以从Gibbs自由能推导得到，其表达式如式（3）和式（4）：

假设在马氏体相变过程中忽略马氏体变体的重定向过程，材料任何微观状态的改变都仅仅是马氏体含量变化的结果。基于此假设，SMA相变过程中的相变应变εt和马氏体体积分数ξ的演化方程（类似于塑性力学理论中的流动法则）可表示为

式中：Λ为相变转换张量，表征了相变应变的方向。

将相变演化方程式（5）代入式（4）的本构关系约束条件中，可得到：

式中：π为与马氏体体积分数ξ共轭的热驱动力。

将Gibbs自由能的显示表达式（1）代入上式，可得到：

基于以上关于马氏体正逆相变的假设以及边界条件的表述，引入相变函数φ（相当于塑性力学中的屈服函数），其表达式为

式中：Y为相变临界屈服值。

基于不同相变过程的假设，SMA唯象理论模型可以通过不同的相变硬化函数来表征。Tanaka指数形式理论模型［9］、Liang和 Rogers余弦模型［10］以及Lagoudas多项式模型［7］的硬化函数分别见式（9），（10），（11）：

2 UMAT材料子程序的编写及验证

在Lagoudas的SMA三维本构关系模型的基础上，利用ABAQUS的二次开发功能编写FORTRAN程序［11］实现SMA的有限元模拟。

针对前面提到的SMA本构关系模型，在UMAT中定义了13个变量（1为马氏体体积分数、2～7为相变应变分量，8～13为个相变方向）和20个材料常数（1～2为弹性模量，3～4为热膨胀系数，5～8为相变温度，9～10为应力影响参数，11～16为初始相变应变，17为最大相变应变，18为泊松比，19为容差，20为初始马氏体体积分数）。

为了验证编写的UMAT子程序的准确性，在ABAQUS软件中建立受单轴加载的单元模型，如图1所示，一端固支约束，一端加载轴向力载荷。载荷加载历史如图2所示。模型中的SMA材料参数采用文献［12］中的等原子钛镍合金参数，即：奥氏体相弹性模量EA为70GPa，马氏体相弹性模量EM为30GPa，泊松比νA和νM均为0.3，奥氏体相膨胀系数αA为22×10-6K-1，马氏体相膨胀系数αM为22×10-6K-1，马氏体相变开始温度Ms为291K（18℃），马氏体相变结束温度Mf为280K（7℃），奥氏体相变开始温度As为295K（22℃），奥氏体相变结束温度Af为306K（33℃），最大相变应变H为0.05。针对超弹性的验证，模型初始相为奥氏体，初始温度为325K（52℃），单轴最大载荷为60kN。针对形状记忆效应的验证，模型初始相为马氏体，初始温度为292K（19℃），单轴最大载荷为30kN，加载结束后升温至350K（77℃），高于奥氏体相变结束温度。

图1 单轴加载模型示意Fig.1 Abridged general view of uniaxial loading model

图2 载荷加载历史示意Fig.2 Abridged general view of applied load history

从图3，4可以看出，以上基于不同硬化函数的三种本构关系模型都能够较好地描述SMA的超弹性和形状记忆效应等特性，相变过程也遵循相应的马氏体奥氏体相变模型表征函数，证明了该UMAT子程序的正确性，从而保证了在UMAT子程序基础上进行的鼓包热变形模拟结果的准确。

图3 基于不同硬化函数的SMA超弹性模型的模拟结果Fig.3 Simulated results of SMA pesudoelasticity model based different hardening functions

图4 基于不同硬化函数的SMA形状记忆效应的模拟结果Fig.4 Simulated results of SMA shape memory effect based different hardening fuctions

3 有限元模拟结果与讨论

3.1 有限元模型的建立

SMA鼓包热变形模型是三维薄壁球体结构的一部分，其尺寸参数包括：过中心点A的弧线在x-z平面上的投影长（图5中的BC线段，10cm）、鼓包厚度（1mm）、曲率半径（42cm）。该鼓包模型在经历双程记忆效应训练之后，沿过中心点的径长方向具有高温伸长、低温收缩的性质。模型网格划分采用了减缩积分的线性六面体单元（C3D8R），在厚度方向上划分5个单元，以免产生沙漏效应。材料参数与前面验证模型中采用的参数相同。

图5 三维鼓包热变形模型的网格划分Fig.5 Meshing of the three demensional bump thermal deformation model

鼓包热变形模型边界条件为周边固支约束，初始温度为291K（18℃），初始状态为马氏体相。沿过中心点的弧长方向具有初始的相变应变，根据不同的训练效果分为-0.5%，-1%，-1.5%，-2%几种情况。考虑鼓包沿厚度方向的初始相变应变，初始马氏体体积分数为1。结构的网格划分如图5所示。SMA鼓包表面加载温度载荷，从初始温度291K（18℃）加载至360K（87℃）并保持温度，增量加载步数为100步，具体过程如图6所示。

图6 温度加载历程示意Fig.6 Abridged general view of temperature loading history

3.2 模拟结果与讨论

3.2.1 定初始相变的SMA鼓包相变特性

从图7中可以看出，鼓包模型在初始状态时为纯马氏体相，体积分数为1；随着温度的逐渐升高，在鼓包中部马氏体相的体积分数开始下降，并逐步向鼓包的四周边界扩展；加载结束后鼓包的大部分区域基本完成了马氏体逆相变，转化为奥氏体相，只有在接近固支约束的边界部分仍残留有马氏体相。

从图8中可以看出，图5中的A、B两点在297K（24℃）左右开始向奥氏体相转变，A点在313K（40℃）结束了马氏体逆相变过程，马氏体体积分数为0，全部转化为奥氏体相；而在B点，由于边界条件的约束，限制了鼓包的马氏体逆相变过程，在加载结束后仍残留有部分马氏体相。

3.2.2 定初始相变应变的SMA鼓包热变形力学特性

由图9可见，A点的竖向挠度在初始阶段增长较慢，在23～40℃（296～313K）的温度区间增长较快，并且大部分的挠度变化发生在该区间内，而且该区间与马氏体逆相变的温度区间相近。可知，鼓包竖向挠度的产生主要与SMA的相变有关。

3.2.3 不同初始相变的SMA鼓包热变形结果

鼓包在加载过程中最大竖向挠度发生在加载结束温度360K。由图10可见，最大竖向挠度与初始相变应变绝对值之间基本呈线性正比关系。当初始相变应变绝对值为2%时，鼓包最大的竖向挠度达到了5mm，相当于鼓包模型BC线段长度的5%。

图7 初始相变应变为-1%时，SMA鼓包模型在不同温度下的马氏体体积分数云图Fig.7 Martensite fraction of SMA bump model at initial phase-transformation strain of-1%and different temperatures

图8 SMA鼓包模型中不同位置处在不同温度下的马氏体体积分数Fig.8 Volume fractions of martensite vs temperture for different locations in SMA bump model

图9 SMA鼓包中心A点处的竖向挠度随温度的变化Fig.9 Change of vertical deflection of the point A（in the middle of the SMA bump）with temperature：（a）vertical deflection of bump and（b）vertical deflection of point A

图10 SMA鼓包在加载过程中最大竖向挠度与初始相变应变绝对值之间的关系Fig.10 The relation ofmaximum vertical deflection and initial phase-transformation strain during SMA bump

从图11中可以看出，在一定的初始相变应变下，随着温度的升高，竖向高度曲线的形状和最大高度也发生改变；竖向高度在313K（40℃）之后的变化较小，该温度与鼓包中A点相变结束点的温度基本吻合。

图11 SMA鼓包过A点弧线各点在不同初始相变下的竖向高度曲线Fig.11 Vertical height curves of the points which pass by the point A at different initial phase-transformation strains

3 结 论

（1）通过ABAQUS软件材料子程序UMAT的接口，编写FORTRAN程序，将SMA的本构关系模型导入到ABAQUS材料库中，从而使有限元软件具备了对SMA模型进行数值模拟的能力；通过建立SMA单轴受力模型，证明了材料子程序UAMT能够同时描述SMA的超弹性和形状记忆效应特性。

（2）对于三维鼓包热变形模型，马氏体逆相变首先发生在鼓包的中部并向四周扩展，除边界部分未完成逆相变外，其余大部分完全转化为奥氏体相；鼓包的挠度最高达到5mm，相当于其弧线投影长度的5%。

致谢：

感谢南京航空航天大学机械结构力学及控制国家重点实验室的鼎力支持。

［1］MIRZAEIFAR R，SHAKERI M，SADIGHI M.Nonlinear finite element formulation for analyzing shape memory alloy cylindrical panels［J］.Smart Materials and Structures，2009，18（3）：1-14.

［2］LAGOUDAS D C.Shape memory alloys：modeling and engineering applications［M］.New York：Springer，2008：1-43.

［3］GAO X，TURNER T L，BURTON D，et al.Finite element analysis of adaptive-stiffening and shape-control SMA hybrid composites［J］.Journal of Engineering Materials and Technology，2006，128（3）：285-293.

［4］REY N，TILLMAN G，MILLER R，et al.Shape memory alloy actuation for a variable area fan nozzle［C］//Proceedings of SPIE，Smart Structures and Materials：Industrial and Commercial Application of Smart Structures Technologies.［S.l.］：［s.n］，2001：371-382.

［5］徐祖耀，江伯鸿.形状记忆材料［M］.上海：上海交通大学出版社，2000：3-6.

［6］WILLIAMS K，CHIU G，BERNHARD R.Adaptive-passive absorbers using shape-memory alloys［J］.Journal of Sound and Vibration，2002，249（5）：835-848.

［7］BO Z，LAGOUDAS D C.Thermomechanical modeling of polycrystalline SMAs under cyclic loading-Part I：theoretical derivations［J］.International Journal of Engineering Science，1999，37：1089-1140.

［8］QIDWAI M A，LAGOUDAS D C.On thermomechanics and transformation surfaces of polycrystalline NiTi shape memory alloy material［J］.International Journal of Plasticity，2000，16：1309-1343.

［9］TANAKA K，NISHIMURA F，HAYASHI T，et al.Phenomenological analysis on subloops and cyclic behavior in shape memory alloys under mechanical and/or thermal loads［J］.Mechanics of Materials，1995，19：281-92.

［10］LIANG C，ROGERS C A.The multi-dimensional constitutive relations of shape memory alloys［J］.Journal of Engineering Mathematics，1992，26：429-443.

［11］庄茁，张帆，岑松，等.ABAQUS非线性有限元分析与实例［M］.北京：科学出版社，2005.

［12］LAGOUDAS D C，BO Z，QIDWAI M A.A unified thermodynamic constitutive model for SMA and finite element analysis of active metal matrix composite［J］.Mechanics of Composite Materials ＆ Structures，1996，3：153-79.