1、基于相场理论的固态相变模拟的有限元计算2012.01-2013.06 工作总结2013.06.27 基于相场理论的固态相变模拟的有限元计算1、相变机理 22、相场理论 32.1相场理论概述 32.2相场的控制微分方程 42.3体积化学自由能 42.4界面能 72.5弹性应变能 82.6多物理场耦合 93、有限元计算 103.1有限元的离散化 103.1.1相场控制微分方程的化简 103.1.2相场与位移场的离散化 103.2相场与位移场的动力求解 113.2.1相场的动力求解 113.2.1位移场的动力求解 123.3相变模拟计算流程 144、已完成工作总结 144.1相场模型的研究总结 14
2、4.1.1弹性应变能的影响 154.1.2体积化学能的影响 164.1.3界面能的影响 164.1.4动力学系数 174.1.5随机噪声 174.1.6稳定结果的非均匀性 174.2有限元算法的研究总结 184.2.1 时域积分 184.2.2单元划分 184.2.3单元阶次 194.2.4高斯点选取 194.2.5计算效率 191相变机理相变是材料从高能量相转变为低能量相， 总能量减小的过程。在材料中，各相能量的高低是由其内部状态(如应力、温度等)所决定的。例如对于NiTi形状 记忆合金，在低温时，马氏体为低能量相，当温度从高温降为低温时，奥氏体转 变为马氏体；在高温时，奥氏体为低能量相，当
3、温度从低温升为高温时，马氏体 转变为奥氏体。相的变体是指晶体结构相同，取向不同的同一种相。例如对于NiTi形状记 忆合金，奥氏体的晶体结构为高对称的立方体( high symmetry cubic)，马氏体 的晶体结构为四方体，在不同的朝向下，共有三种变体，如图 1-1。马氏体的三种变体图1-1奥氏体与马氏体的晶体结构示意图在一定的内部状态下，体系生成各变体后，总能量的大小是不同的。体系会 自发地选择生成使得体系能量降低最多的变体。如图 1-2所示，各变体首先在使 得体系能量降低最多的区域形成，后生长形成特定的位向关系。2、相场理论2.1相场理论概述相场理论是一种用来模拟相变微观组织的演化过程
4、的理论， 它的主要特点是以Ginzburg-Laudau相变理论为基础，基于扩散界面模型，引入连续变量来描述 新旧两相，用微分方程描述系统自由能的变化。相场理论的特点具体如下：(1)离散变量的连续化：在相变中，各相或者相的各变体是离散的物理量，各相或者相的各变体之间不存在连续变化的过程。相场模型将这些离散变量连续 化，用序参量n表示，便于在数值上获取导数等信息。例如，对于新相含有 n个变体的相变，采用序参量n i = n 2= = n n=0表示体系此处为母相所占有；序参 量np=1或-1且n=0 (k=1,2.p-1,p+1, 沁表示此处对应为新相的第p个变体。由 于连续化的处理，各相或者各
5、变体间的界面是弥散界面 (diffuse in terface)，而不是尖锐界面(Sharp interfac，如图2-1。(2)微分方程：建立相场内的平衡微分方程，描述各相或各变体的演化。(3)界面能：基于扩散界面模型建立界面能，描述相场内部序参量的变化。母相新相母相新相图2-1尖锐界面与弥散界面示意图尖锐界面弥散界面2.2相场的控制微分方程在相场理论中，相变微观组织演化通过求解如下的 Ginzburg-Landau方程:式中，L为相场动力学系数，控制着相场的演化速率； G为系统总的自由能，在马氏体相变中，一般来说，包括体积化学能、界面能和弹性应变能； ；(r,t)为高斯随机噪声，满足以下条
6、件：:k r,t) = 0(2-2) 耳(r,t)Up(r,t)沖 2kBTLkp%6(r -r )芬(t -t)(2-3)其中，：:： r,t) 为相场热起伏对空间和时间的平均；kB为Boltzmann常数，T为温度，S为Kronecker Delta函数。2.3体积化学自由能体积化学能是指材料本身具有的一种化学势能，各相之间的体积化学能密度 与温度相关，同一相中的各变体具有相同的体积化学能密度。 两相间的体积化学 能密度之差促使相变从能量高的相向能量低的相转变。在相变中，体积化学自由能具有以下两点性质(1)在母相处，体积化学自由能密度为0,在新相处，体积化学自由能密度为-AChem； (2
7、)当某点为母相或新相时，该点的体积化学能驱动力为 0,处于稳定态或亚稳定态。这两点性质要求体积化学自由能密度函数满足：WchemCl =吩2 =叫=。)=。=0(kP)=1 八P (2-4)Fchem | Fchem nM=o= 虫，q=o(k#)= 0p p在不同的文献中，体积化学能密度函数 屁hem主要有以下两种形式:(1) 2-4-6次多项式-4)n /(2-5)chem hem轨 2，2 2)4( 4 4 2 4Ao 2 2 2 3(12 n) -1 - -1)6A1, A2, Ao为体系膨胀系数，一般取正数；n为变体总数；序参量n=n= = n=0 表示体系此处为母相所占有；序参量
8、n=且n=0 (k=1,2.p-1,p+1,n表示此 处对应为新相的第p个变体。当n=0 ( k=1,2.p-1,p+1,n时，体积化学能密度 艇hem与序参量n的关系 曲线见图2-2与图2-3。由图所示，该插值函数有三个稳定点： n=0, 。1在该插 值函数下，相变从母相到新相转化时，需先从外界吸收一部分能量，达到跃迁所 需的能量后(如图2-3中n= 0.1，才能向外界放出能量，完成相变转化。0.2chem( Achem)n1-1.5*-芦1111.5T-0.2/-0.4-0.6 - /-0.8 7-1 -7-1.2图2-2 2-4-6多项式体积化学能密度插值函数与序参量的关系曲线(2) 2
9、-3-4次多项式3=fA1( 2 2 2)-役(3 ；(2-6)2 3企(12 2 2)2 (oi)4与2-4-6次多项式类似，Ai, A2, A3为体系膨胀系数，一般取正数；n为变 体总数；序参量n=n= = n=o表示体系此处为母相所占有；序参量 叩=1且n=o (k=1,2.p-1,p+1,n表示此处对应为新相的第p个变体。当n=0 ( k=1,2.p-1,p+1,n时，体积化学能密度 艇hem与序参量n的关系 曲线见图1-3与图1-4。由图所示，该插值函数有两个稳定点： n=0,1。与2-4-6次多项式类似，在该插值函数下，相变从母相到新相转化时，需先从外界吸收一 部分能量，达到跃迁所
10、需的能量后(如图 1-4中n=0.011)，才能向外界放出能 量，完成相变转化。图2-4 2-3-4多项式体积化学能密度插值函数与序参量的关系曲线图2-5 2-3-4多项式体积化学能密度插值函数与序参量的关系曲线（靠近 0部分）2.4界面能在相场模型中，界面能表示的是由于相场不均匀分布所产生的能量。 在相变 中，界面能的释放将将促使相场分布的均匀化。界面能与序参量的梯度相关，其 表达式为：1 njPi pdV (2-7)D 2 p 4其中，川p为界面梯度能系数；一般为了简化问题，将材料考虑为各向同性，且 各变体的界面梯度能系数相等，界面能表达式可化简为：1 nGchem 八 p2dV(2-8)
11、D 2 P 42.5弹性应变能在相变中，由于晶体结构的改变，从母相到新相，材料将发生相变变形。在 变形协调下，相变变形使得材料的变形重分布，改变材料的弹性应变能。对于线弹性材料，弹性应变能密度的表达式为：二(:E ：f/2 (2-9)在相变中，总应变包含弹性应变与相变应变：；.=e 亠 j j j (2-10)其中，相变应变在相场模型中采用序参量 n构造插值函数0 (来)表示：n時=迟%g(p)(2-11)式中，；j0(p)为第p个变体由于晶体结构改变产生的相变应变，是材料常数。弹性应变能密度的表达式可展开为：e=二 E ： ,/2 工 ： E : 2 - : E : : E : /2 (2-
12、12)因此，弹性应变能驱动力为:2.6多物理场耦合相变是从高能量相转变为低能量相的过程， 由于对于材料，相的能量高低是由内部状态所决定的，相变中微观组织演化将受到其它物理场的影响。 同时，相变中微观组织的演化也会改变其它物理场的分布。因此，相变本身是一个多物理 场耦合的问题。在马氏体相变中，相互耦合的三个物理场分别是相场、位移场与 温度场，三者的关系如图2-6所示。相场在实际模拟中，一般忽略相变过程中的发热、吸热，假定在相变过程中，温度不发生改变，原三场耦合的问题简化为两场耦合，如图 2-7所示。图2-7简化后的马氏体相变两场耦合示意图对于相场与位移场的耦合，本研究计算采用顺序耦合方法求解。具
13、体求解策 略是按照顺序进行位移场和相场的分析， 采用上一时刻的相场分布，先求解位移 场，得到上一时刻位移场分布，将位移场所得结果应用到相场求解中，再求解这 一时刻的相场分布。3、有限元计算3.1有限元的离散化 3.1.1相场控制微分方程的化简对体积化学能、界面能和弹性应变能三种能量之和的总能量泛函求变分， 可 得:将式.:t屮的门+代佝亠E d (叭Vp,i)p=1,2,3.n(3-1)p d(3-1)代入式(2-1 )的 Ginzburg-Landau方程，可得:L二 :卩=_l（ stchem厂； p =1,2,3.n（3-2） 別p式(3-2)与非稳态的温度场控制方程类似，在有限元离散化
14、中，只需保证 C0连续，只需采用常规的有限元插值。3.1.2相场与位移场的离散化将相场中的变量序参量n用有限元离散化，如图3-1,采用节点变量表示序 参量n将式(3-2)离散化，可得：C旦止K讣二Qk 12 门小 k = 1,2,3.n (3-3) ct(1)C n由单元矩阵cn组装后所得c 二N TN d(3-4)(2)K n由单元矩阵k n 组装所得k p B TL -B d“ (3-5)(3)Qk包含相变驱动力-和随机项E k;q = - ;N T(L工 k)d(3-6)4 3图3-1有限元离散化示意图同理，位移场的动力方程或者静力方程可离散化为：M KU二R - 2 （3-7）.t或：
15、KU二R 1 . 2 n /（1） M为单元质量矩阵m组装后所得m二.NJNd（3-8）（2） K为单元刚度矩阵k组装所得k = iBTEBd（3-9）（3） R包含外荷载fext和相变引起的等效荷载fequi3.2相场与位移场的动力求解3.2.1相场的动力求解相场中ti+1时刻的序参量与ti时刻的序参量的关系通过下式表示:咕严（1-映贰+陋化+（3-10）将上述递推公式代入离散后的相场方程，整理可得k pi .1 =Q?p (3-11)其中等效刚度阵和等效右端项为：J?二C *Kc (3-12)Q?p =(1 J：tQpi JtQpi i (C(1J：tKc) pi (3-13)该求解方法数
16、值稳定的最大时间步为：2 _=t 匚 (3-14)入为方程(KJ- 4C)T =0中的最大特征值。(1-2max当3=0,该方法为时域显式积分法，式(3-8)和式(3-9)可化为冏二C(3-15)Qk =MQki +(C-加Kc)Si(3-16)3.2.1位移场的动力求解(1) Newmark-3法：在Newmark- 3法中ti+1时刻的位移、速度、加速度与ti时刻的位移、速度、 加速度的关系通过下式表示：54=5 (1 -)tu上匕 1 (3-17)* q =* 二tu(1 - Jit2* “2* 1 (3-18)将式(3-17)和式(3-18)递推公式位移场的动力方程，整理后得逐步积分
17、公式为：kui + =P?Uy Y Yui 厂阪(ui + -ui)+(1-p)ui t(1-丑)ui (3-19)- 1 1-1 -ui + =i2(ui 十ui)-(巫 T)ui其中，等效刚度阵和等效荷载向量分别为：1 y k=K t2M-C (3-20)Newmark- B法的稳定性条件为：1 1t 兰十=Tn (3-23) 兀 J2 JY2P当y1/2 , B=1/4时，算法为无条件稳定，同时具有二阶精度(2)中心差分法：在中心差分法中ti时刻的速度、加速度与ti+1时刻和ti-1时刻的位移的关系通过下式表示：将递推公式代入力平衡方程，得:肌爪二乳(3-26)其中：2 1 1? 1 =
18、Ri -(K-M)ui -(M C)uy(3-28)心t 23中心差分法的稳定性条件为：进乞( 3-29)313.3相变模拟计算流程对于位移场和相场的耦合问题，采用的求解策略是在根据t时刻序参量分布值，先在t时刻先求解位移场，再根据得到t时刻的应力分布值，求解t+ 0时刻 相场，得到t+ At时刻的序参量分布值，具体流程如图 2-8所示：t=t+ t图3-1相变模拟计算流程图4、已完成工作总结4.1相场模型的研究总结马氏体相变的相场模型包含三种能量：体积化学能、界面能、弹性应变能, 三种能量相互作用，改变其中任何一种，都将影响到最终结果。4.1.1弹性应变能的影响根据弹性应变能的表达式（2-1
19、2）与弹性应变能驱动力的表达式（2-14）, 弹性应变能的组成分为三部分：（1） 应力：应力的大小影响着相变的演化速率和演化方向。当应力的方向 与相变应变一致时，应力促进相变演化；反之，抑制相变演化。从应力产生的原因来分，应力可分为外荷载作用下引起的应力和相变应变 引起的应力。由于相变应变引起的应力一定与相变应变相反， 因此，相变应变引起的应力总是抑制相变进一步演化。（2） 弹性模量：在线弹性下，弹性模量与相变应变引起的应力大小成线性 关系，因此，弹性模量越大，由相变引起的应力大小越大，相变的演化越慢，反 之，相变演化越快。（算例验证）（3） 泊松比：泊松比的大小与相变应变引起的应力大小相关，
20、其影响与弹性模量相似，只是对不同的应变分量的影响程度不同。 （算例验证）（4） 相变应变插值函数：相变应变的插值函数影响弹性应变能在序参量值 介于0-1之间的大小，进而影响弹性应变能驱动力的大小，影响相变演化进程与 最终结果。在选取相变应变插值函数时，需注意以下准则:II.非必须满足：为了使序参量位于 0或者1时，材料处于稳定或者亚稳定态， 其能量驱动力应为0,由于体积化学能驱动力在序参量位于 0或者1时为0,弹 性应变能驱动力也应为0。因此，相变应变插值函数 9 （须满足：表达式形式一次导数表达式线性O =np p%=1二次日=H2p p儿三次6 = n3p p日p=3叫：p按照以上选取准则
21、，已进行过对比的相变插值函数如下:2-3次% =3听 2鸣日p =6睥 6听2-4次%=2p日 p=4Hp4 吧2-3-4 次2 3 4Tp = 0.07Hp +4.14p 3.07耳卩 2 3日 p =-0.14J +12.42*!p-12.28其中，在材料内部无缺陷且无外力作用下，对比线性、二次、三次插值函数,线性的结果好于二次和三次，原因是随着插值函数次数增大，在序参量靠近 1的部分弹性应变能驱动力增大，最终稳定时，序参量离1的平均距离变远。与此 同时，体积化学能驱动力在1附近的增大幅度大于弹性应变能驱动力， 因此，应 变能驱动力在总驱动力的比例降低，最终形成的位向关系变差。 （算例验证
22、）2-3次、2-4次、2-3-4次是在材料内部无缺陷且无外力作用时，结果无明显 差别。对于材料受到较大外力作用，或因为自身缺陷引起应力集中，造成相变区域 为大应力区域的情况，还未进行算例分析。4.1.2体积化学能的影响在温度不变的情况下，体积化学能只与自身的插值函数相关， 与应力大小无 关。体积化学能是否促进相变演化取决于在选取的温度下，母相的体积化学能是 否较大。4.1.3界面能的影响由于相场模型本身是将离散的相通过连续的序参量值来表示， 界面能采用梯度能的描述方式，与实际情况并不相符。实际各相或者各变体间是否存在界面能， 界面能的大小与什么参数相关，还需进一步阅读文献。单就相场模型中的界面
23、能，可以从以下角度分析：（1） 数学表达式：界面能表示由于相场不均匀分布所产生的能量。在相变中，界面能的释放将将促使相场分布的均匀化。（2） 模型中的作用：通过梯度描述的界面能，在相场模型中，建立相变区 域中相邻各点的联系，使得新相或者变体能够通过“扩散”的方式生长。（3） 界面能系数对模拟结果影响：界面能系数的大小影响相变模拟结果中各变体的平均尺寸和弥散界面的宽度，界面能系数越大，各变体的平均尺寸越大, 弥散界面的宽度越大。（算例验证）4.1.4动力学系数在式（2-1）的Ginzburg-Landau方程中，动力学系数控制着相变演化的速率。 由于相变演化本身还与位移场相互耦合， 相变演化的速
24、率影响惯性力的大小，决定位移场分析是采用动力分析还是静力分析。4.1.5随机噪声相场模型是一个描述相变演化过程的模型， 无法描述相变的形核机理，通常 要求其化学能驱动力与弹性应变能驱动力为 0，在新相形核前，体系处于亚稳定 状态。引入随机的噪声，主要目的是引起序参量值的微小改变， 从而使得相变能 够演化下去。因此，在随机噪声很小的情况下，与空间坐标独立的随机噪声抽样方式和持 续时间不影响相变结果的位向关系和变体的平均尺寸。 （算例验证）4.1.6稳定结果的非均匀性根据相场的控制方程，当结果稳定时，原方程可化为:(Gchem G GQ p(X,t)=0(4-11)花(r) -瞇假设变体中间区域，
25、序参量为0或者1的均匀分布，界面能为0,方程进- 步简化为：丄h f p（r） f P（r） （4-13）在变体边界，序参量不均匀分布，存在界面能，弹性应变能驱动力与体积化学能驱动力之和不为0:-0(4-14)由于式（4-13）与式（4-14）互斥，满足式（4-13）的序参量值必然不满足式（4-14）o变体中间区域的序参量值必然不等于变体边界上的序参量值， 因此,中间区域的序参量值与边界上的序参量值存在一定的梯度， 变体中间区域的序参量值不满足式（4-13）,与假设矛盾。4.2有限元算法的研究总结4.2.1时域积分（1） 显式积分与隐式积分：在完全积分的情况下，由于求解相场驱动力与相变引起的等
26、效荷载所需的计算时间较长，同线性回代时间差不多，因此，当每 步求解过程中，刚度阵不发生改变时，无论是在相场计算还是位移场计算（动力 分析）中采用显式积分不能明显地提高效率。由于显式积分的等效刚度阵为对角阵，相比隐式积分未经优化的变带宽等效 刚度阵，占用内存较少。应用显式积分，可以大幅降低内存使用。（2） 积分步长：步长是否合理的判断：I 是否满足相场数值计算的稳定性要求；II.序参量是否大比例超限。在合理步长下，会造成相变结果局部区域结果的不同， 但所有结果整体上都 符合位向关系。同时在合理步长下，步长较大，虽然计算结果收敛时对应的相变 实际时间较长，但在模拟中所需要的计算步数较少，需要的数值
27、计算时间较短。4.2.2单元划分在位移场中，提高网格精度，可以更准确地描述区域内应力分布；在相场中, 提高网格精度， 可以更准确地描述区域内序参量的梯度分布， 更好地描述界面能 的大小，计算结果更精确合理。同时，提高网格精度，在图形显示上，可以更好 地描述变体间的边界，避免出现锯齿型区域。4.2.3 单元阶次在二维分析中，对于位移场，相比 Q4单元，Q8单元可以更准确地描述单 元内应力分布。对于相场， Q4 单元采用的双线性插值可以保证单元内部任一点 的值不超出节点值的范围，而 Q8单元采用的二次插值则不能。当节点序参量值 落在0-1之间时，Q8单元插值给出的内部值可能会超出许可值范围。4.2
28、.4高斯点选取由于目前判断结果是否合理最重要的依据是位向关系， 而位向关系在相变演 化过程中无法体现， 应用不完全积分虽然计算不精确， 并不影响最终结果的位向 关系。4.2.5计算效率在计算中，主要耗时部分的分为四部分：相场驱动力计算（ 33%左右），相场线性回代（20%左右），相变引起的等效荷载 （22%左右），位移场线性回代（25% 左右）。其中，相场驱动力积分计算与相变引起的等效荷载积分计算可以采用 1 个高斯点进行分析，所需的时间减为原先的 1/10；相场的线性回代可采用对角 阵的显式积分，所需时间减为原先的 1/10；对于位移场，在动力分析中，可以采用对角阵的显式积分，将所需时间减为
29、原先的 1/10，但对于静力分析，目前 暂无较好的方法提高效率。参考文献1潘金生，仝健民，田民波 . 材料科学基础 M. 清华大学出版社有限公司 , 1998.2Wang Y, Khachaturyan A G. Three-dimensional field model and computer modeling of martensitic transformationsJ. Acta materialia. 1997, 45(2): 759-773.3Jin Y M, Khachaturyan A G. Phase field microelasticity theory of dislocation dynamics in a polycrystal: model and three-dimensional simulationsJ. Philosophical magazine letters.
