CN112131762B - 模拟马氏体相变的网格自适应有限元方法 - Google Patents

Abstract

本发明公开了一种模拟马氏体相变的网格自适应有限元方法，应用于形状记忆合金的工业生产模拟设计和科学研究领域。该方法通过建立马氏体相变相场有限元模型，设定合适的网格单元后验误差估计，标记需要加密/粗化的网格单元并执行加密/粗化步骤，达到有效追踪马氏体变体界面并在界面处加密网格和在非界面处粗化网格的效果，在保持计算精度的条件下，减小相场模拟计算所需求解的自由度总数目，有效提高计算效率和推动模拟研究进程。本发明方法能有效动态追踪动态变化的马氏体变体界面，所需的计算资源和计算时间更少，能有效提高计算效率和推动模拟研究进程。

Description

模拟马氏体相变的网格自适应有限元方法

技术领域

本发明涉及一种马氏体相变的相场模型的有限元方法数值计算模拟，应用于形状记忆合金的数值模拟技术领域。

背景技术

在形状记忆合金材料的工业生产和科学研究中，形状记忆合金的马氏体相变行为需要进行详细的研究。相场方法是在介观尺度下模拟材料微结构演化的一种数值模拟方法。在过去的几十年中，相场方法已经成为介观尺度下模拟材料的微观结构演化的有力工具，并成功地被应用于模拟马氏体相变过程中的材料微观结构演化现象。马氏体相变相场模型的控制方程是一组耦合程度较高的偏微分方程组，在实际模拟过程中需要通过数值计算求解偏微分方程组。常用于求解相场模型偏微分方程组的数值计算方法有傅里叶谱方法、有限元方法、有限差分方法。其中，有限元方法由于其能求解复杂的几何区域和边界条件，成为了数值求解相场模型的一种主要方法。有限元方法是基于求解域的网格划分进行计算的，因此为了确保模拟计算在动态变化的马氏体变体界面处的计算精度不得不进行全局网格加密，使得有限元模型的总体自由度大大增加，耗费大量计算资源和时间。

因此需要对有限元方法进行优化以减小马氏体相变相场模拟计算所需的计算资源和时间。网格自适应有限元通过加密/粗化局部区域的网格减小所需求解的总体自由度数目，达到在满足精度要求的情况下减小计算资源耗费的效果。目前的网格自适应有限元方法多用于求解静态问题，在求解动态演化问题时比较少见。因此研究和使用基于网格自适应有限元的马氏体相变相场数值模拟方法，对于工业生产和科学研究都是很有必要的，成为亟待解决的技术问题。

发明内容

为了解决现有技术问题，本发明的目的在于克服已有技术存在的不足，提供一种模拟马氏体相变的网格自适应有限元方法，达到有效追踪马氏体变体界面并在界面处加密网格和在非界面处粗化网格的效果，在保持计算精度的条件下，减小相场模拟计算所需求解的自由度总数目，有效提高计算效率和推动模拟研究进程。本发明方法能有效动态追踪动态变化的马氏体变体界面，所需的计算资源和计算时间更少，效率更高，成本更低。

为达到上述发明创造目的，本发明采用如下技术方案：

一种模拟马氏体相变的网格自适应有限元方法，包括以下步骤：

(1)建立马氏体相变的相场模型，同时建立相应的有限元模型；

(2)进行有限元前处理；

(3)编写相应的网格自适应有限元程序代码，将相应有限元程序代码与开源有限元程序deal.II结合调用，选取误差估计准则计算每个网格单元的误差估计，通过误差估计的大小判断网格是否需要加密/粗化，在模拟计算中实现动态划分网格，追踪马氏体变体的界面并在界面处加密网格，在非界面区域粗化网格，降低总体网格数量，从而减小模型的自由度总数达到减小模拟计算所需的计算资源的效果；

(4)编写后处理脚本，通过后处理提取所需的信息和结果。本发明模拟马氏体相变的网格自适应有限元方法，可以在动态变化的马氏体变体界面处加密网格，在非界面处粗化网格，在确保达到同样初始网格密度下的精度要求的情况下减小模拟计算所需的计算资源，为快速模拟计算研究形状记忆合金的马氏体相变行为提供更快的计算方法。

优选地，在所述步骤(3)中，选取Kelly-error-estimator作为马氏体相变相场模型的网格后验误差估计，网格单元K的Kelly后验误差估计表达式为：

其中h_K为网格单元K的网格大小，u为计算Kelly误差时的精确解，

为拉普拉斯算子；而在数值模拟计算中无法得到精确解和精确的二阶导数从而无法求得Kelly后验误差，需要对Kelly后验误差进行数值简化：

其中u_h为通过数值模拟计算求得的数值解，

为数值解u_h的梯度在网格单元边界上的点x_i上的左右极限的差值：

其中

为梯度算子，ε为任意小的正数。

优选地，在所述步骤(3)中，根据网格后验误差的大小按照自定义的规则标记应加密/粗化的网格，并执行加密/粗化网格操作，所述标记规则采用如下任意一种：

第一种标记规则：选取后验误差占总的网格后验误差的一定比例的网格最小集合进行加密标记，选取后验误差占总的网格后验误差的一定比例的网格最小集合进行粗化标记；

第二种标记规则：将网格按照网格后验误差的大小进行排序，选取一定比例的排序靠前的网格进行加密标记，选取一定比例的排序靠后的网格进行粗化标记；同时还可指定网格最多加密次数和最多粗化次数。

优选在于步骤(3)中选取位移变量的数值解计算网格后验误差估计。因为马氏体相变是无扩散的结构型相变，在马氏体变体界面附近的位移场将会被扭曲，界面附近的网格单元的Kelly后验误差将会非常大，选取位移变量的数值解计算网格后验误差估计能通过网格单元的Kelly后验误差的大小有效地追踪马氏体变体界面。

优选地，在所述步骤(1)中，相场模型中使用时间相关的金茨堡朗道方程(TimeDependent Ginzburg Landau equation，TDGL方程)来描述马氏体相变的演化过程，同时还需要满足力学平衡方程：

σ_ij，j+f＝0，

其中η_p(x，t)为相场模型的场变量也成为序参数，代表着第p个马氏体变体；L为动力学系数；F为总自由能；δF/δη_p为变分导数；f为体积力；σ_ij为应力张量，通过位移场和本构关系求得：

σ＝C(η_p)：ε^el(u)＝C(η_p)：(ε^tot(u)-ε^MT(η_p))，

其中C(η_p)为与序参数相关的弹性张量；ε^el(u)为弹性应变；ε^tot(u)为从位移场u求得的总应变；ε^MT(η_p)为与序参数相关的相变转化应变，由化学能F_ch、弹性能F_el、界面能F_gr组成：

F＝F_ch+F_el+F_gr，

F_ch＝∫_Vf_ch(η_p，T)dV，

其中f_el为弹性能密度，f_ch为化学自由能密度，β_ij为梯度能系数；

根据数学推导，上述相场模型的伽辽金弱解形式为：

其中

为对应于序参数场η_p的标量形函数；φ为对应位移场的u的向量形函数；t为应力边界Γ_t上的外加应力；TDGL方程使用隐式欧拉格式离散，非线性方程组使用牛顿迭代法进行求解。

本发明与现有技术相比较，具有如下显而易见的突出实质性特点和显著优点：

1.本发明方法在动态变化的马氏体变体界面处加密网格，在非界面处粗化网格，在确保达到同样初始网格密度下的精度要求的情况下减小模拟计算所需的计算资源，为快速模拟计算研究形状记忆合金的马氏体相变行为提供更快的计算方法；

2.本发明方法有效解决了常规有限元方法模拟计算中在动态变化的马氏体变体界面处的计算精度不得不全局加密网格，使得有限元模型的总体自由度大大增加，耗费大量计算资源和时间的计算效率问题，可以有效动态追踪动态变化的马氏体变体界面，所需的计算资源和计算时间更少，能有效提高计算效率和推动模拟研究进程。

附图说明

图1是计算方法的流程图。

图2是以TiNb合金材料为例，计算模拟温度诱导马氏体相变时的网格自适应示意图。

图3是以TiNb合金材料为例，计算模拟应力诱导马氏体相变时的应力应变曲线对比图。

图4是以TiNb合金材料为例，计算模拟应力诱导马氏体相变时的网格自适应示意图。

图5是以TiNb合金材料为例，计算模拟应力诱导马氏体相变时的系统总自由度变化图。

具体实施方式

以下结合具体的实施例子对上述方案做进一步说明，本发明的优选实施例详述如下：

实施例一：

在本实施例中，一种模拟马氏体相变的网格自适应有限元方法，包括以下步骤：

(1)建立马氏体相变的相场模型，同时建立相应的有限元模型；

(2)进行有限元前处理；

(3)编写相应的网格自适应有限元程序代码，将相应有限元程序代码与开源有限元程序deal.II结合调用，选取误差估计准则计算每个网格单元的误差估计，通过误差估计的大小判断网格是否需要加密/粗化，在模拟计算中实现动态划分网格，追踪马氏体变体的界面并在界面处加密网格，在非界面区域粗化网格，降低总体网格数量，从而减小模型的自由度总数达到减小模拟计算所需的计算资源的效果；

(4)编写后处理脚本，通过后处理提取所需的信息和结果。

本实施例模拟马氏体相变的网格自适应有限元方法能在动态变化的马氏体变体界面处加密网格，在非界面处粗化网格，在确保达到同样初始网格密度下的精度要求的情况下减小模拟计算所需的计算资源，为快速模拟计算研究形状记忆合金的马氏体相变行为提供更快的计算方法。

实施例二：

本实施例与实施例一基本相同，特别之处在于：

在本实施例中，以Ti2448合金为例，计算模拟了温度在100K时的温度诱导马氏体相变。在所述步骤(1)中，相场模型中使用时间相关的金茨堡朗道方程(Time DependentGinzburg Landau equation，TDGL方程)来描述马氏体相变的演化过程，同时还需要满足力学平衡方程：

σ_ij，j+f＝0，

其中η_p(x，t)为相场模型的场变量也成为序参数，代表着第p个马氏体变体；L为动力学系数；F为总自由能；δF/δη_p为变分导数；f为体积力；σ_ij为应力张量，通过位移场和本构关系求得：

σ＝C(η_p)：ε^el(u)＝C(η_p)：(ε^tot(u)-ε^MT(η_p))，

其中C(η_p)为与序参数相关的弹性张量；ε^el(u)为弹性应变；ε^tot(u)为从位移场u求得的总应变；ε^MT(η_p)为与序参数相关的相变转化应变，由化学能F_ch、弹性能F_el、界面能F_gr组成：

F＝F_ch+F_el+F_gr，

F_ch＝f_V，f_ch(η_p，T)dV，

其中f_el为弹性能密度，f_ch为化学自由能密度，β_ij为梯度能系数；

根据数学推导，上述相场模型的伽辽金弱解形式为：

其中

为对应于序参数场η_p的标量形函数；φ为对应位移场的u的向量形函数；t为应力边界Γ_t上的外加应力；TDGL方程使用隐式欧拉格式离散，非线性方程组使用牛顿迭代法进行求解。

在所述步骤(3)中，选取Kelly-error-estimator作为马氏体相变相场模型的网格后验误差估计，网格单元K的Kelly后验误差估计表达式为：

其中h_K为网格单元K的网格大小，u为计算Kelly误差时的精确解，

为拉普拉斯算子；对Kelly后验误差进行数值简化：

其中u_h为通过数值模拟计算求得的数值解，

为数值解u_h的梯度在网格单元边界上的点x_i上的左右极限的差值：

其中

为梯度算子，ε为任意小的正数。

按照本实施例方法，选取位移变量的数值解计算网格后验误差估计。因为马氏体相变是无扩散的结构型相变，在马氏体变体界面附近的位移场将会被扭曲，界面附近的网格单元的Kelly后验误差将会非常大，选取位移变量的数值解计算网格后验误差估计能通过网格单元的Kelly后验误差的大小有效地追踪马氏体变体界面。计算流程图如图1所示，在使用建立的网格自适应有限元方法模拟的马氏体相变过程中，能明显看到在马氏体变体的界面处网格得到加密，在非界面处得到粗化，如图2所示。在相同硬件条件下的相同的模拟过程不同方法所花费计算时间如表1所示，可以看出所建立的模拟马氏体相变方法可以有效地减小耗费的计算资源和时间。

实施例三：

本实施例与前述实施例基本相同，特别之处在于：

在本实施例中，以Ti2448合金为例，计算模拟了温度在200K时的应力诱导马氏体相变。参见图1，一种模拟马氏体相变的网格自适应有限元方法，包括如下步骤：

(1)建立马氏体相变的相场模型及有限元模型：

相场模型中使用时间相关的金茨堡朗道方程(Time dependent Ginzburg Landauequation，TDGL)来描述马氏体相变的演化过程，同时还需要满足力学平衡方程：

σ_ij，j+f＝0，

其中η_p(x，t)为相场模型的场变量也成为序参数，代表着第p个马氏体变体；L为动力学系数；F为总自由能；δF/δη_p为变分导数；f为体积力；σ_ij为应力张量，可以通过位移场和本构关系求得：

σ＝C(η_p)：ε^el(u)＝C(η_p)：(ε^tot(u)-ε^MT(η_p))，

其中C(η_p)为与序参数相关的弹性张量；ε^el(u)为弹性应变；ε^tot(u)为从位移场u求得的总应变；ε^MT(η_p)为与序参数相关的相变转化应变，不同的模型下的相变转化应变的计算方式不同。同时不同模型下的总自由能构建方法也不同，由化学能F_ch、弹性能F_el、界面能F_gr组成：

F＝F_ch+F_el+F_gr，

F_ch＝∫_Vf_ch(η_p，T)dV，

其中f_el为弹性能密度，f_ch为化学自由能密度，β_ij为梯度能系数；以上马氏体相变相场模型的本构模型、化学自由能密度函数和材料参数都是与模拟材料相关的，通过实验或者第一性原理计算获得；

根据数学推导，上述相场模型的伽辽金弱解形式为：

其中

为对应于序参数场η_p的标量形函数；φ为对应位移场的u的向量形函数；t为应力边界Γ_t上施加的外加应力；TDGL方程使用隐式欧拉格式离散，非线性方程组使用牛顿迭代法进行求解；

(2)进行有限元前处理：

设置初始网格、边界条件、初始条件，处理多种马氏体相变的诱导方式，包括温度诱导、应力诱导；

(3)编写相应的有限元程序代码，将相应有限元程序代码与开源有限元程序deal.II结合调用，编写的有限元程序代码包括：

(3.1)网格设定，边界条件，初始条件；

(3.2)材料模型及总体刚度矩阵、总体残差的计算；

(3.3)使用牛顿迭代法求解非线性的方程组；

(3.4)网格后验误差指标的设定和计算：

采用Kelly-Error-Estimator作为网格后验误差指标(Kelly，D.W.，De S.R.Gago，J.P.，Zienkiewicz，O.C.and Babuska，I.(1983)，A posteriori error analysis andadaptive processes in the finite element method：Part I-error analysis.Int.J.Numer.Meth.Engng.，19：1593-1619.doi：10.1002/nme.1620191103)，在实际计算过程中，对于网格单元K的Kelly-Error-Estimator(e_K)通过下面的公式进行计算：

其中h_K为网格大小，u为选择的变量，u_h为对应变量u的数值逼近解，

为拉普拉斯算子，

为使用有限差分方法计算的第i个点的左右梯度的差值：

由于马氏体相变是结构型相变，所以在计算网格的Kelly后验误差时选取位移为计算的变量；

(3.5)根据网格后验误差的大小按照自定义的规则标记应加密/粗化的网格，并执行加密/粗化网格操作，所述标记规则采用如下任意一种：

第一种标记规则：选取后验误差占总的网格后验误差的一定比例的网格最小集合进行加密标记，选取后验误差占总的网格后验误差的一定比例的网格最小集合进行粗化标记；

第二种标记规则：将网格按照网格后验误差的大小进行排序，选取一定比例的排序靠前的网格进行加密标记，选取一定比例的排序靠后的网格进行粗化标记；同时还；指定网格最多加密次数和最多粗化次数；

(3.6)将自由度重新分配到新的网格中；

(4)编写后处理脚本，通过后处理提取所需的信息和结果。

本实施例以Ti2448合金材料为例，在环境温度为200K时使用所建立的方法模拟计算应力诱导马氏体相变，计算流程图如图1所示，并与相同初始网格的常规有限元方法的计算结果进行对比,如图3所示本实施例中所建立的网格自适应有限元方法计算的应力应变曲线与常规有限元方法计算得到的应力应变曲线基本一致，验证所建立的方法的正确性。同时如图4所示，在使用建立的网格自适应有限元方法模拟的马氏体相变过程中，能明显看到在马氏体变体的界面处网格得到加密，在非界面处得到粗化。系统总自由度个数随无量纲演化时间变化的曲线如图5所示，可以看出随着演化进行系统的总自由度个数随演化时间进行不断变化，总体比初始总自由度个数大大减小。在相同硬件条件下的相同的模拟过程不同方法所花费计算时间如表1所示，证明所建立的方法能有效动态追踪马氏体变体界面和通过网格自适应减小耗费的计算资源和时间。

表1相同硬件条件下的相同的模拟过程耗费的计算时间对比

	温度诱导马氏体相变	应力诱导马氏体相变
			建立的网格自适应有限元方法	8810s	123000s
常规有限元方法	21000s	344000s

本实施例模拟马氏体相变的网格自适应有限元方法，应用于形状记忆合金的工业生产模拟设计和科学研究领域。该方法通过建立马氏体相变相场有限元模型，设定合适的网格单元后验误差估计，标记需要加密/粗化的网格单元并执行加密/粗化步骤，达到有效追踪马氏体变体界面并在界面处加密网格和在非界面处粗化网格的效果，可以在保持计算精度的条件下，减小相场模拟计算所需求解的自由度总数目，有效提高计算效率和推动模拟研究进程。本发明方法有效地解决了常规有限元方法进行马氏体相变相场模拟时，为了确保模拟计算在动态变化的马氏体变体界面处的计算精度不得不全局加密网格，使得有限元模型的总体自由度大大增加，耗费大量计算资源和时间的计算效率问题，可以有效动态追踪动态变化的马氏体变体界面，所需的计算资源和计算时间更少，能有效提高计算效率和推动模拟研究进程。

上面对本发明实施例结合附图进行了说明，但本发明不限于上述实施例，还可以根据本发明的发明创造的目的做出多种变化，凡依据本发明技术方案的精神实质和原理下做的改变、修饰、替代、组合或简化，均应为等效的置换方式，只要符合本发明的发明目的，只要不背离本发明的技术原理和发明构思，都属于本发明的保护范围。

Claims (4)

1.一种模拟马氏体相变的网格自适应有限元方法，其特征在于：包括以下步骤：

(1)建立马氏体相变的相场模型，同时建立相应的有限元模型；

(2)进行有限元前处理；

(3)编写相应的网格自适应有限元程序代码，将相应有限元程序代码与开源有限元程序deal.II结合调用，选取误差估计准则计算每个网格单元的误差估计，通过误差估计的大小判断网格是否需要加密/粗化，在模拟计算中实现动态划分网格，追踪马氏体变体的界面并在界面处加密网格，在非界面区域粗化网格，降低总体网格数量，从而减小模型的自由度总数，达到减小模拟计算所需的计算资源的效果；

(4)编写后处理脚本，通过后处理提取所需的信息和结果。

2.根据权利要求1所述模拟马氏体相变的网格自适应有限元方法，其特征在于，在所述步骤(3)中，选取Kelly-error-estimator作为马氏体相变相场模型的网格后验误差估计，网格单元K的Kelly后验误差估计表达式为：

其中h_K为网格单元K的网格大小，u为计算Kelly误差时的精确解，

为拉普拉斯算子；对Kelly后验误差进行数值简化：

其中u_h为通过数值模拟计算求得的数值解，

为数值解u_h的梯度在网格单元边界上的点x_i上的左右极限的差值：

其中

为梯度算子，ε为任意小的正数；由于马氏体相变是无扩散的结构型相变，选取位移变量的数值解计算网格后验误差。

3.根据权利要求1所述模拟马氏体相变的网格自适应有限元方法，其特征在于，在所述步骤(3)中，根据网格后验误差的大小按照自定义的规则标记应加密/粗化的网格，并执行加密/粗化网格操作，所述标记规则采用如下任意一种：

第一种标记规则：选取后验误差占总的网格后验误差的一定比例的网格最小集合进行加密标记，选取后验误差占总的网格后验误差的一定比例的网格最小集合进行粗化标记；

第二种标记规则：将网格按照网格后验误差的大小进行排序，选取一定比例的排序靠前的网格进行加密标记，选取一定比例的排序靠后的网格进行粗化标记；同时还指定网格最多加密次数和最多粗化次数。

4.根据权利要求1所述模拟马氏体相变的网格自适应有限元方法，其特征在于，在所述步骤(1)中，相场模型中使用时间相关的金茨堡朗道方程来描述马氏体相变的演化过程，同时还需要满足力学平衡方程：

σ_ij，j+f＝0

其中η_p(x，t)为相场模型的场变量也称为序参数，代表着第p个马氏体变体；L为动力学系数；F为总自由能；δF/δη_p为变分导数；f为体积力；σ_ij为应力张量，通过位移场和本构关系求得：

σ＝C(η_p)：ε^el(u)＝C(η_p)：(ε^tot(u)-ε^MT(η_p))，

其中C(η_p)为与序参数相关的弹性张量；ε^el(u)为弹性应变；ε^tot(u)为从位移场u求得的总应变；ε^MT(η_p)为与序参数相关的相变转化应变，由化学能F_ch、弹性能F_el、界面能F_gr组成：

F＝F_ch+F_el+F_gr，

f_ch＝∫_Vf_ch(η_p，T)dV，

其中f_el为弹性能密度，f_ch为化学自由能密度，β_ij为梯度能系数；

根据数学推导，上述相场模型的伽辽金弱解形式为：

其中

为对应于序参数场η_p的标量形函数；φ为对应位移场的u的向量形函数；t为应力边界Γ_t上的外加应力；TDGL方程使用隐式欧拉格式离散，非线性方程组使用牛顿迭代法进行求解。

Images