郝 林，陈子雄

(上海飞机设计研究院，上海 201210)

形状记忆合金(shape memory alloys, SMAs)是一种强非线性功能材料，因其具有超弹性效应、形状记忆效应、高阻尼能力和生物相容性等特性被广泛应用于航空航天器结构中[1-3]。其特殊的性质给仿真分析带来一定的困难，诸多学者曾尝试对形状记忆合金建立其本构模型以求在仿真过程中能真实地再现其热力响应力学特性[4]，如Lagoudas团队提出的SMAs三维唯象本构模型能较好地模拟三维形状记忆合金结构的热力学特性[5-6]。然而，由于热载荷引起相变时形状记忆合金材料的横向剪切刚度不是常数[7]，因此采用软件自带的壳体单元无法满足描述形状记忆合金薄壳结构的热力位移响应的精度要求，而大型薄壳结构若采用三维单元进行仿真分析则代价非常大。

本文基于一阶剪切变形理论(FSDT)建立了一种新的适用于形状记忆合金薄壳结构分析的4节点、20自由度壳单元，采用Lagoudas本构模型模拟形状记忆合金材料特性，运用Newton-Raphson方法迭代求解单元矩阵。最后针对给定初始相变应变条件的形状记忆合金薄壳结构进行分析，验证新单元的有效性。

1 形状记忆合金本构方程

基于材料内部和外部的状态变量，用热力学势函数来表征材料一点的热力学状态。基于Lagoudas本构模型[8]，形状记忆合金的Gibbs自由能与应力张量、温度和马氏体体积分数的关系如式(1)所示：

(1)

式中：G(σ,T,ξ,εt)为材料中一点的Gibbs自由能；σ,T,T0,εt,ξ分别为应力张量、温度、参考温度、相变应变张量和马氏体体积分数；S,α,c,ρ,s0,u0分别为柔度张量、热胀系数、比热、密度、熵和参考温度下的内能；f(ξ)为反映相变应变演化的硬化函数。f(ξ)的具体形式如式(2)所示：

(2)

式中：Δs0,Ms,Mf,As,Af,Δu0分别为熵增、马氏体相变起始温度、马氏体相变结束温度、马氏体逆相变起始温度、马氏体逆相变结束温度以及内能变化量。

根据Gibbs自由能方程可得到形状记忆合金材料的应力-应变关系[5]，如式(3)所示：

(3)

式中：ε为应变。

类比于塑性流动法则，假设形状记忆合金相变应变与马氏体体积分数有如下关系：

(4)

式中：Λ为用以描述相变应变演化方向与幅值的张量。根据Lagoudas提出的观点[5,9],形状记忆合金材料相变演化的起始与结束时刻所遵循的路径严格遵循式(5)所示的Clausius-Planck不等式：

(5)

式中:π为相变驱动力。

于是形状记忆合金发生相变时的相变函数可假定为以下形式:

(6)

式中：Φ为形状记忆合金相变函数；Y为相变热力学力的临界值。马氏体体积分数演化的约束条件采用Kuhn-Tucker条件，具体表示为：

(7)

式中：Φ(σ,T,ξ)为相变函数。

所给出的形状记忆合金本构方程式(3)依然是一组非线性方程组，求解依旧存在很大的困难，为此采用位移增量法进行求解。

为能够得到形状记忆合金的切线刚度张量，对形状记忆合金本构关系关系式进行微分，写作如下形式：

dσ=R∶dε+N ∶dT

(8)

式中：R为切线刚度张量；N为热弹性张量。

以增量形式改写式(3)并代入式(4)，得到：

dσ=S-1∶{dε-αdT-[ΔS∶σ+Δα(T-T0)+Λ]dξ}

(9)

式中：ΔS为柔度；Δα为热胀系数。

考虑到相变函数具体形式，则有：

(10)

式中：∂σΦ为相变函数对应力导数。

形状记忆合金发生相变时有Φ=0。相变函数全微分形式为：

dΦ=∂σΦ∶dσ+∂TΦdT+∂ξΦdξ=0

(11)

通过消除马氏体体积分数dξ，并使用式(10)、(11)改写式(9)即可得到切线刚度张量以及热弹性张量具体形式，如式(12)所示：

(12)

式中：∓代表马氏体正/逆相变过程。

2 壳单元模型

形状记忆合金壳单元采用4节点单元，每个节点有5个自由度，即4节点20自由度单元。对于每个自由度，使用线性形函数。为了使单元适应复杂边界，采用如图1所示的等参数变换。

图1 壳单元等参变换示意图

2.1 位移场假设

基于一阶剪切变形假设，形状记忆合金壳单元中任意一点位移模式具有如下形式：

(13)

式中：u(x,y,z)、v(x,y,z)、w(x,y,z)分别为该点在X,Y,Z方向的总位移；符号(•)0表示中性面的值；Ψx，Ψy分别为过点的法线绕X，Y轴转角。基于一阶剪切理论变形示意图如图2所示。

图2 基于一阶剪切理论变形示意图

于是可通过几何方程得到应变为：

(14)

2.2 单元建立

单元的平衡方程可通过虚功原理建立：

(15)

式中：U为广义位移；u为位移向量；ΨU为广义力；f为单元体力向量；τ为面力向量；δ是变分符号；Γ为边界；Ω为体积。

形状记忆合金壳单元在每个节点上有5个自由度，采用形函数Ni(x,y), 单元的自由度表示为如下形式：

(16)

式中：nodes为具体点的编号。

3 非线性分析过程

由公式(15)可知，形状记忆合金本构方程求解为强非线性问题。综合考虑公式(15)及(16)，应力平衡方程的弱形式可表述为：

Ψ4}TdA+F)=0

(17)

其中，

(18)

式中：σx，σxy，…为不同方向应力；elm表示单元信息；nelm为单元数；Ae为单元面积。

式(17)大括号中的表达式表示单元的内力, 其中F包含了外面施加的面力与场力。

由于形状记忆合金本构方程求解是强非线性响应问题，因此本文采用Newton-Raphson迭代法以及向后差分积分方法求解。分析时，载荷以增量形式进行施加，每个增量步施加后均需求解式(19)：

-t+Δt[KU]t+Δt{ΔU}=t+Δt{ΨU}

(19)

式中：KU为单元刚度阵；Δt为时间步增量；ΔU为时间步增量Δt下广义位移U的增量。

(20)

4 基于ABAQUS软件UEL建立

为了在ABAQUS中实现所开发的形状记忆合金壳单元，基于ABAQUS软件用户单元子程序(UEL)功能将开发单元数学模型按ABAQUS软件UEL要求编写成FORTRAN子程序。为了可视化有限元分析结果，使用了虚拟单元技术和子程序UVARM。虚拟单元和展开单元具有相同的几何特征，可以如实反映UEL分析所得结果。

5 应用分析以及结果讨论

为了检验所开发单元的实际应用能力，本文将其应用于一个形状记忆合金薄壳结构进行仿真分析。形状记忆合金薄壳模型示意图如图3所示，其几何尺寸为长140 mm、宽20 mm、厚0.55 mm，初始时室温下的挠度约为5 mm。设薄壳沿着X方向存在大小为-0.5%的初始相变应变，形状记忆合金初始处于完全马氏体相，此时具有升温鼓起的形状记忆效应。表1、表2分别给出了材料参数以及部分模型信息。

图3 形状记忆合金薄壳模型

仿真过程中假定形状记忆合金薄壳结构与外界绝热，结构温度从260 K线性上升至290 K，并保持恒定，温度加载历程如图4所示。形状记忆合金薄壳结构随着温度的升高在Y轴方向上的挠度增大，如图5所示，加载结束后，薄壳中心点变形挠度达到2.57 mm。图6所示为形状记忆合金薄壳中心点Y方向挠度随温度变化历程，由图可知，随着薄壳温度升高，尽管结构因热胀冷缩效应产生变形，但在Y方向上挠度变化并不明显。当温度高于262 K时，此时薄壳开始发生相变(图7)，形状记忆合金开始向高温形状转变，Y方向上变形挠度迅速增加。当温度达到280 K时，薄壳完全相变为奥氏体相，随后因热胀效应其挠度有所增加。

表1 材料参数

图4 加载历程

表2 模型信息

图5 温度加载结束后位移分布

图6 薄壳中心点Y方向挠度随温度变化历程

图7 薄壳马氏体体积分数随温度变化历程

图8所示为加载过程中形状记忆合金薄壳中性面形状变化过程。从图中可以看出，随着温度逐渐升高，薄壳在Y方向上产生的变形挠度也在增大，且挠度最大位置始终在薄壳中心区域。

图8 薄壳随温度变形过程

6 结束语

本文提出了一种4节点20自由度的形状记忆合金壳单元，并采用一阶剪切理论描述其位移场，使该壳单元能够有效地模拟形状记忆合金的热力行为。基于UEL功能将所开发的壳单元纳入ABAQUS软件实现仿真应用，用所开发的形状记忆合金壳单元对悬臂薄壳结构进行数值分析，结果表明，该单元能够很好地预测SMA薄壳结构的时间响应问题。