基于相场理论的固态相变模拟的有限元计算

1、 工作总结-基于相场理论的固态相变模拟的有限元计算目录1、相变机理12、相场理论32.1相场理论概述32.2相场的控制微分方程32.3体积化学自由能42.4界面能72.5弹性应变能72.6多物理场耦合83、有限元计算93.1有限元的离散化93.1.1 相场控制微分方程的化简93.1.2 相场与位移场的离散化103.2相场与位移场的动力求解113.2.1 相场的动力求解113.2.1 位移场的动力求解123.3相变模拟计算流程134、已完成工作总结144.1相场模型的研究总结144.1.1 弹性应变能的影响144.1.2 体积化学能的影响164.1.3 界面能的影响164.1.4 动力学系数17

2、4.1.5 随机噪声174.1.6 稳定结果的非均匀性174.2有限元算法的研究总结184.2.1 时域积分184.2.2 单元划分184.2.3 单元阶次194.2.4 高斯点选取194.2.5 计算效率191、相变机理相变是材料从高能量相转变为低能量相，总能量减小的过程。在材料中，各相能量的高低是由其内部状态（如应力、温度等）所决定的。例如对于NiTi形状记忆合金，在低温时，马氏体为低能量相，当温度从高温降为低温时，奥氏体转变为马氏体；在高温时，奥氏体为低能量相，当温度从低温升为高温时，马氏体转变为奥氏体。相的变体是指晶体结构相同，取向不同的同一种相1。例如对于NiTi形状记忆合金，奥氏体

3、的晶体结构为高对称的立方体（ high symmetry cubic），马氏体的晶体结构为四方体，在不同的朝向下，共有三种变体，如图1-1。图1-1 奥氏体与马氏体的晶体结构示意图在一定的内部状态下，体系生成各变体后，总能量的大小是不同的。体系会自发地选择生成使得体系能量降低最多的变体。如图1-2所示，各变体首先在使得体系能量降低最多的区域形成，后生长形成特定的位向关系。图1-2 变体产生及位向关系形成示意图2、相场理论2.1相场理论概述相场理论是一种用来模拟相变微观组织的演化过程的理论，它的主要特点是以Ginzburg-Laudau相变理论为基础，基于扩散界面模

4、型，引入连续变量来描述新旧两相，用微分方程描述系统自由能的变化。 相场理论的特点具体如下：（1）离散变量的连续化：在相变中，各相或者相的各变体是离散的物理量，各相或者相的各变体之间不存在连续变化的过程。相场模型将这些离散变量连续化，用序参量表示，便于在数值上获取导数等信息。例如，对于新相含有n个变体的相变，采用序参量1=2=n=0表示体系此处为母相所占有；序参量p=1或-1且k=0（k=1,2.p-1,p+1,n）表示此处对应为新相的第p个变体。由于连续化的处理，各相或者各变体间的界面是弥散界面（diffuse interface），而不是尖锐界面（Sharp interface），如图2-1

5、。（2）微分方程：建立相场内的平衡微分方程，描述各相或各变体的演化。（3）界面能：基于扩散界面模型建立界面能，描述相场内部序参量的变化。图2-1 尖锐界面与弥散界面示意图2.2相场的控制微分方程在相场理论中，相变微观组织演化通过求解如下的Ginzburg-Landau方程： （2-1）式中，L为相场动力学系数，控制着相场的演化速率；G为系统总的自由能，在马氏体相变中，一般来说，包括体积化学能、界面能和弹性应变能；为高斯随机噪声，满足以下条件： （2-2） （2-3）其中，为相场热起伏对空间和时间的平均；kB为Boltzmann常数，T为温度，为Kronecker Delta函数。2.3体积化学

6、自由能体积化学能是指材料本身具有的一种化学势能，各相之间的体积化学能密度与温度相关，同一相中的各变体具有相同的体积化学能密度。两相间的体积化学能密度之差促使相变从能量高的相向能量低的相转变。在相变中，体积化学自由能具有以下两点性质（1）在母相处，体积化学自由能密度为0，在新相处，体积化学自由能密度为-fchem；（2）当某点为母相或新相时，该点的体积化学能驱动力为0，处于稳定态或亚稳定态。这两点性质要求体积化学自由能密度函数满足： （2-4）在不同的文献中，体积化学能密度函数chem主要有以下两种形式：（1）2-4-6次多项式2 （2-5）A1，A2，A3为体系膨胀系数，一般取正数；n为变体总

7、数；序参量1=2=n=0表示体系此处为母相所占有；序参量p=±1且k=0（k=1,2.p-1,p+1,n）表示此处对应为新相的第p个变体。当k=0（k=1,2.p-1,p+1,n）时，体积化学能密度chem与序参量p的关系曲线见图2-2与图2-3。由图所示，该插值函数有三个稳定点：p=0,±1。在该插值函数下，相变从母相到新相转化时，需先从外界吸收一部分能量，达到跃迁所需的能量后（如图2-3中p=±0.1），才能向外界放出能量，完成相变转化。图2-2 2-4-6多项式体积化学能密度插值函数与序参量的关系曲线图2-3 2-4-6多项式体积化学能密度插值函数与序参量的

8、关系曲线（靠近0部分）（2）2-3-4次多项式3 （2-6） 与2-4-6次多项式类似，A1，A2，A3为体系膨胀系数，一般取正数；n为变体总数；序参量1=2=n=0表示体系此处为母相所占有；序参量p=1且k=0（k=1,2.p-1,p+1,n）表示此处对应为新相的第p个变体。当k=0（k=1,2.p-1,p+1,n）时，体积化学能密度chem与序参量p的关系曲线见图1-3与图1-4。由图所示，该插值函数有两个稳定点：p=0,1。与2-4-6次多项式类似，在该插值函数下，相变从母相到新相转化时，需先从外界吸收一部分能量，达到跃迁所需的能量后（如图1-4中p=0.011），才能向外界放出能量，完

9、成相变转化。 图2-4 2-3-4多项式体积化学能密度插值函数与序参量的关系曲线图2-5 2-3-4多项式体积化学能密度插值函数与序参量的关系曲线（靠近0部分）2.4界面能在相场模型中，界面能表示的是由于相场不均匀分布所产生的能量。在相变中，界面能的释放将将促使相场分布的均匀化。界面能与序参量的梯度相关，其表达式为： （2-7）其中，为界面梯度能系数；一般为了简化问题，将材料考虑为各向同性，且各变体的界面梯度能系数相等，界面能表达式可化简为： （2-8）2.5弹性应变能在相变中，由于晶体结构的改变，从母相到新相，材料将发生相变变形。在变形协调下，相变变形使得材料的变形重分布，改变材料的弹性应变

10、能。对于线弹性材料，弹性应变能密度的表达式为： （2-9）在相变中，总应变包含弹性应变与相变应变： （2-10）其中，相变应变在相场模型中采用序参量构造插值函数()来表示： （2-11）式中，为第p个变体由于晶体结构改变产生的相变应变，是材料常数。弹性应变能密度的表达式可展开为： （2-12）因此，弹性应变能驱动力为： （2-13）2.6多物理场耦合相变是从高能量相转变为低能量相的过程，由于对于材料，相的能量高低是由内部状态所决定的，相变中微观组织演化将受到其它物理场的影响。同时，相变中微观组织的演化也会改变其它物理场的分布。因此，相变本身是一个多物理场耦合的问题。在马氏体相变中，相互耦合的三

11、个物理场分别是相场、位移场与温度场，三者的关系如图2-6所示。图2-6 马氏体相变三场耦合示意图在实际模拟中，一般忽略相变过程中的发热、吸热，假定在相变过程中，温度不发生改变，原三场耦合的问题简化为两场耦合，如图2-7所示。图2-7 简化后的马氏体相变两场耦合示意图对于相场与位移场的耦合，本研究计算采用顺序耦合方法求解。具体求解策略是按照顺序进行位移场和相场的分析，采用上一时刻的相场分布，先求解位移场，得到上一时刻位移场分布，将位移场所得结果应用到相场求解中，再求解这一时刻的相场分布。3、有限元计算3.1有限元的离散化3.1.1 相场控制微分方程的化简对体积化学能、界面能和弹性应变能三种能量之

12、和的总能量泛函求变分，可得： （3-1）将式（3-1）代入式（2-1）的Ginzburg-Landau方程，可得： （3-2）式（3-2）与非稳态的温度场控制方程类似，在有限元离散化中，只需保证C0连续，只需采用常规的有限元插值。3.1.2 相场与位移场的离散化将相场中的变量序参量用有限元离散化，如图3-1，采用节点变量表示序参量，将式（3-2）离散化，可得： （3-3）（1）C由单元矩阵c组装后所得 （3-4）（2）K由单元矩阵k组装所得 （3-5）（3）Qk包含相变驱动力和随机项k （3-6）图3-1 有限元离散化示意图同理，位移场的动力方程或者静力方程可离散化为： （3-7）或：（1）M

13、为单元质量矩阵m组装后所得 （3-8）（2）K为单元刚度矩阵k组装所得 （3-9）（3）R包含外荷载fext和相变引起的等效荷载fequi3.2相场与位移场的动力求解3.2.1 相场的动力求解相场中ti+1时刻的序参量与ti时刻的序参量的关系通过下式表示： （3-10）将上述递推公式代入离散后的相场方程，整理可得 （3-11）其中等效刚度阵和等效右端项为： （3-12） （3-13）该求解方法数值稳定的最大时间步为： （3-14）为方程中的最大特征值。当=0，该方法为时域显式积分法，式（3-8）和式（3-9）可化为 （3-15） （3-16）3.2.1 位移场的动力求解（1）Newmark-法

14、：在Newmark-法中ti+1时刻的位移、速度、加速度与ti时刻的位移、速度、加速度的关系通过下式表示： （3-17） （3-18） 将式（3-17）和式（3-18）递推公式位移场的动力方程，整理后得逐步积分公式为： （3-19）其中，等效刚度阵和等效荷载向量分别为： （3-20） （3-22） Newmark-法的稳定性条件为： （3-23）当=1/2，=1/4时，算法为无条件稳定，同时具有二阶精度。（2）中心差分法：在中心差分法中ti时刻的速度、加速度与ti+1时刻和ti-1时刻的位移的关系通过下式表示： （3-24） （3-25）将递推公式代入力平衡方程，得： （3-26）其中： （3

15、-27） （3-28）中心差分法的稳定性条件为： （3-29）3.3相变模拟计算流程对于位移场和相场的耦合问题，采用的求解策略是在根据t时刻序参量分布值，先在t时刻先求解位移场，再根据得到t时刻的应力分布值，求解t+t时刻相场，得到t+t时刻的序参量分布值，具体流程如图2-8所示： 图3-1 相变模拟计算流程图4、已完成工作总结4.1相场模型的研究总结马氏体相变的相场模型包含三种能量：体积化学能、界面能、弹性应变能，三种能量相互作用，改变其中任何一种，都将影响到最终结果。4.1.1 弹性应变能的影响根据弹性应变能的表达式（2-12）与弹性应变能驱动力的表达式（2-14），弹性应变能的组成分为三

16、部分：（1）应力：应力的大小影响着相变的演化速率和演化方向。当应力的方向与相变应变一致时，应力促进相变演化；反之，抑制相变演化。从应力产生的原因来分，应力可分为外荷载作用下引起的应力和相变应变引起的应力。由于相变应变引起的应力一定与相变应变相反，因此，相变应变引起的应力总是抑制相变进一步演化。（2）弹性模量：在线弹性下，弹性模量与相变应变引起的应力大小成线性关系，因此，弹性模量越大，由相变引起的应力大小越大，相变的演化越慢，反之，相变演化越快。（算例验证）（3）泊松比：泊松比的大小与相变应变引起的应力大小相关，其影响与弹性模量相似，只是对不同的应变分量的影响程度不同。（算例验证）（4）相变应变

17、插值函数：相变应变的插值函数影响弹性应变能在序参量值介于0-1之间的大小，进而影响弹性应变能驱动力的大小，影响相变演化进程与最终结果。在选取相变应变插值函数时，需注意以下准则：I.必须满足：为保证当序参量为0时，相变应变为0；当序参量为1时，相变应变为，弹性应变能插值函数()必须满足： （4-1） II.非必须满足：为了使序参量位于0或者1时，材料处于稳定或者亚稳定态，其能量驱动力应为0，由于体积化学能驱动力在序参量位于0或者1时为0，弹性应变能驱动力也应为0。因此，相变应变插值函数()须满足： （4-2）按照以上选取准则，已进行过对比的相变插值函数如下： 表达式形式一次导数表达式线性二次三次

18、2-3次2-4次2-3-4次其中，在材料内部无缺陷且无外力作用下，对比线性、二次、三次插值函数，线性的结果好于二次和三次，原因是随着插值函数次数增大，在序参量靠近1的部分弹性应变能驱动力增大，最终稳定时，序参量离1的平均距离变远。与此同时 ，体积化学能驱动力在1附近的增大幅度大于弹性应变能驱动力，因此，应变能驱动力在总驱动力的比例降低，最终形成的位向关系变差。（算例验证）2-3次、2-4次、2-3-4次是在材料内部无缺陷且无外力作用时，结果无明显差别。对于材料受到较大外力作用，或因为自身缺陷引起应力集中，造成相变区域为大应力区域的情况，还未进行算例分析。4.1.2 体积化学能的影响在温度不变的

19、情况下，体积化学能只与自身的插值函数相关，与应力大小无关。体积化学能是否促进相变演化取决于在选取的温度下，母相的体积化学能是否较大。4.1.3 界面能的影响由于相场模型本身是将离散的相通过连续的序参量值来表示，界面能采用梯度能的描述方式，与实际情况并不相符。实际各相或者各变体间是否存在界面能，界面能的大小与什么参数相关，还需进一步阅读文献。单就相场模型中的界面能，可以从以下角度分析：（1）数学表达式：界面能表示由于相场不均匀分布所产生的能量。在相变中，界面能的释放将将促使相场分布的均匀化。（2）模型中的作用：通过梯度描述的界面能，在相场模型中，建立相变区域中相邻各点的联系，使得新相或者变体能够

20、通过“扩散”的方式生长。（3）界面能系数对模拟结果影响：界面能系数的大小影响相变模拟结果中各变体的平均尺寸和弥散界面的宽度，界面能系数越大，各变体的平均尺寸越大，弥散界面的宽度越大。（算例验证）4.1.4 动力学系数在式（2-1）的Ginzburg-Landau方程中，动力学系数控制着相变演化的速率。由于相变演化本身还与位移场相互耦合，相变演化的速率影响惯性力的大小，决定位移场分析是采用动力分析还是静力分析。4.1.5 随机噪声相场模型是一个描述相变演化过程的模型，无法描述相变的形核机理，通常要求其化学能驱动力与弹性应变能驱动力为0，在新相形核前，体系处于亚稳定状态。引入随机的噪声，主要目的是

21、引起序参量值的微小改变，从而使得相变能够演化下去。因此，在随机噪声很小的情况下，与空间坐标独立的随机噪声抽样方式和持续时间不影响相变结果的位向关系和变体的平均尺寸。（算例验证）4.1.6 稳定结果的非均匀性根据相场的控制方程，当结果稳定时，原方程可化为： （4-11） （4-12）假设变体中间区域，序参量为0或者1的均匀分布，界面能为0，方程进一步简化为： （4-13）在变体边界，序参量不均匀分布，存在界面能，弹性应变能驱动力与体积化学能驱动力之和不为0： （4-14）由于式（4-13）与式（4-14）互斥，满足式（4-13）的序参量值必然不满足式（4-14）。变体中间区域的序参量值必然不等于

22、变体边界上的序参量值，因此，中间区域的序参量值与边界上的序参量值存在一定的梯度，变体中间区域的序参量值不满足式（4-13），与假设矛盾。4.2有限元算法的研究总结4.2.1 时域积分（1）显式积分与隐式积分：在完全积分的情况下，由于求解相场驱动力与相变引起的等效荷载所需的计算时间较长，同线性回代时间差不多，因此，当每步求解过程中，刚度阵不发生改变时，无论是在相场计算还是位移场计算（动力分析）中采用显式积分不能明显地提高效率。由于显式积分的等效刚度阵为对角阵，相比隐式积分未经优化的变带宽等效刚度阵，占用内存较少。应用显式积分，可以大幅降低内存使用。（2）积分步长：步长是否合理的判断：I是否满足相

23、场数值计算的稳定性要求；II序参量是否大比例超限。在合理步长下，会造成相变结果局部区域结果的不同，但所有结果整体上都符合位向关系。同时在合理步长下，步长较大，虽然计算结果收敛时对应的相变实际时间较长，但在模拟中所需要的计算步数较少，需要的数值计算时间较短。4.2.2 单元划分在位移场中，提高网格精度，可以更准确地描述区域内应力分布；在相场中，提高网格精度，可以更准确地描述区域内序参量的梯度分布，更好地描述界面能的大小，计算结果更精确合理。同时，提高网格精度，在图形显示上，可以更好地描述变体间的边界，避免出现锯齿型区域。4.2.3 单元阶次在二维分析中，对于位移场，相比Q4单元，Q8单元可以更准确地描述单元内应力分布。对于相场，Q4单元采用的双线性插值可以保证单元内部任一点的值不超出节点值的范围，