相场法数值 模拟.. 共26页
相场法数值 模拟
四、相场方程 (phase-field equations)
Ginzburg-Landau 方程 f 0 k (r , t ) F Lk k k Lk t k (r , t ) k
Cahn-Hilliard 方程
2、凝固-单相场变量
ϵ 2 * 2 F f 0 ( xB , , T ) (xB ) ( ) dr V 2 2
等温凝固，假设摩尔体积不变，即组成梯度项不考虑 ϵ则=0
均质自由能密度
f 0 ( xB , , T * ) f p ( xB , , T * ) g ( )
图1（a）性能不连续
（b）性能连续
N. Moelans, B. Blanpain , P. Wollants, "An introduction to phase-field modeling of microstructure evolution", CALPHAD -Computer Coupling of Phase Diagrams and Thermochemistry, 32, 268-294, 2008
相场量(phase-fields)
两相
多相
p相，相应的变量 k 在系统中任一点 r
1表示在固相中 0表示在液相中 0 1在固液界面
k 1
p
k
1，k 0, k
三、热力学势函数 (thermodynamic energy functional)
经典热力学
F ( xB ,k ) f 0 ( xB ,k ) 1 xB (r , t ) M M xB (r , t ) Vm t xB (r , t ) xB
制冷剂中冰晶生长的数值模拟
制冷剂中冰晶生长的数值模拟
制冷剂中含有水份会导致冰堵从而影响制冷系统的的正常工作。
为探究水份在制冷剂中凝固的过程,根据凝固理论以及相场法,在微观尺度下对水的凝固过程进行研究。
采用数值模拟的方法,研究了各向异性系数和过冷度对冰晶生长的影响。
结果表明,各向异性系数大小对冰晶形貌大小有较大的影响,取值越大,冰晶的二次分枝越发达。
同时过冷度也明显影响冰晶的生长和形貌,过冷度越大,冰晶生长速度越快,二次分枝越发达。
模拟计算的冰晶生长结果与实际显微镜下观察结果较接近吻合。
三元合金等温凝固过程的相场法模拟
文章编号:167325196(2007)0420001204三元合金等温凝固过程的相场法模拟王智平1,2,冯 力1,路 阳1,朱昌盛1,肖荣振1(1.兰州理工大学材料科学与工程学院,甘肃兰州 730050;2.兰州理工大学甘肃省有色金属及复合材料工程技术研究中心,甘肃兰州 730050)摘要:采用相场与浓度场耦合的相场法,以Fe2C2P合金为例,模拟三元合金等温凝固的生长过程.结果表明,在一定的温度区间里,生长速度是随着温度的降低而减小的,随着温度的降低,二次枝晶生长越发达.模拟结果还表明,二次枝晶之间有大量的溶质富集;各向异性系数的改变会影响二次枝晶的形貌.关键词:相场法;三元合金;等温凝固;模拟中图分类号:T G248;T G244 文献标识码:ASimulation of isothermal solidif ication of ternary alloy with phase2f ield model WAN G Zhi2ping1,2,FEN G Li1,U Yang1,ZHU Chang2sheng1,XIAO Rong2zhen1(1.College of Materials Science and Engineering,Lanzhou Univ.of Tech.,Lanzhou 730050,China;2.Research Centre of Gansu Non2fer2 rous and Composite Materials,Lanzhou Univ.of Tech.,Lanzhou 730050,China)Abstract:The p hase2field model in which t he concent ration2field and p hase2field were coupled was applied to simulate micro st ruct ural evolution of isot hermal dendritic growt h of Fe2C2P ternary alloy.The calcula2 tion result indicated t hat,in a definite temperat ure range,t he growt h speed decreased wit h temperat ure re2 ducing,and,instead,t he growt h of side2branches became flourishing.The result also indicated t hat,a2 mong t he side2branches,t he solute concent rated;t he anisotropy coefficient could have effect on t he mor2 p hology of dendrite growt h.K ey w ords:p hase2field met hod;ternary alloy;isot hermal solidification;simulation 凝固过程微观组织的控制是铸造工作者致力研究的课题,过去的组织优化研究是对不同条件下制备的试样进行金相检验,以得到其组织形成的规律,这样要耗费大量的人力、物力和财力,并且有一定的盲目性.对铸件凝固过程和微观组织进行计算机模拟不仅可以减少无谓的劳动,而且做少量实验就可以达到控制微观组织的目的,并可获得主要工艺参数与铸件凝固组织的定量关系,为通过工艺控制改善微观组织提供可靠依据.相场法是一种用于描述在非平衡状态相变中复杂相界面演变的强有力工具.相场方程的解可描述金属系统中固、液界面的形态、曲率和界面的移动,从而避免跟踪复杂的固、液界面.目前该方法是微观组织模拟中的热点. 本文采用K im等人在WBM相场模型基础上 收稿日期:2006211217 基金项目:甘肃省自然科学基金(3ZS0612A252032) 作者简介:王智平(19562),男,山东荷泽人,教授,博导.改进得到的KKS模型[1～5]及耦合浓度场方程模拟了Fe2C2P三元合金在等温凝固过程中枝晶生长的演化及溶质的分布.1 相场模型[3～6]1.1 相场控制方程用稀释溶液近似来定义自由能.合金组元的化学势分别用μ1,μ2和μ3来表示,其中μs1(c1s)=μ0s1+R TV mln(γs1c1s)(1)μl1(c1l)=μ0l1+R TV mln(γl1c1l)(2)μs2(c2s)=μ0s2+R TV mln(γs2c2s)(3)μl2(c2l)=μ0l2+R TV mln(γl2c2l)(4)μs3(c1s,c2s)=μ0s3+R TV mln(1-(c1s+c2s))(5)μl3(c1l,c2l)=μ0l3+R TV mln(1-(c1l+c2l))(6)第33卷第4期2007年8月兰 州 理 工 大 学 学 报Journal of Lanzhou University of TechnologyVol.33No.4Aug.2007式中:下标1和2分别代表两种不同的溶质成分,下标3代表溶剂成分.上标l,s和0分别代表液相,固相和标准状态.γ是活性系数,在以后的控制方程推导中会被约去.f s=c1sμs1+c2sμs2+(1-(c1s+c2s))μs3(7) f l=c1lμl1+c2lμl2+(1-(c1l+c2l))μl3(8) f(φ,c1,c2)=h(φ)f s+(1-h(φ))f l+W g(φ)(9)式中:f s、f l和W分别代表固、液相的自由能和双阱势高;h(φ)为势函数,h(φ)=φ3(10-15φ+6φ2); g(φ)为剩余自由能,g(φ)=φ2(1-φ)2.相场控制方程可以表示为5φ5t=M[ε2 2φ-fφ](10)式中:M为相场参数,fφ表示自由能密度对相场的一阶导数,ε是与界面能有关的参数,可表示为ε(θ)=ε(1+υco s(kθ))(11)式中:k为各向异性的模数,通常取4;υ为各向异性强度系数;θ为界面法线方向与优先生长方向间的夹角,tanθ=φy/φx.1.2 溶质场扩散方程在合金的模拟中耦合了溶质扩散方程,其表达式为5c15t= D1(φ)f c1c1f c1(12)5c2t= D2(φ)f c2c2f c2(13)式中:D1(φ)和D2(φ)分别代表了两种溶质的扩散系数,f c1、f c2、f c1c1和f c2c2为自由能密度f分别对两个溶质浓度的一阶、二阶偏微分.在界面区域的溶质浓度c是固相和液相的质量分数的和,并且在两相平衡时,界面区域中任意点的固相和液相的化学势相等,即 c j=h(φ)c s j+(1-h(φ))c l j (j=1,2)(14) μs j(c s j(x,t)=μl j(c l j(x,t))(15)2 相场参数的确定及模拟计算2.1 相场参数的确定[6～8]相场迁移率参数M的表达式为M-1=ε3σ2W1D1iζ1(c e1l,c e1s)+1D2iζ2(c e2l,c e2s)(16)ζj=R TV m(c e jl-c e j s)2×∫10h(φ)[1-h(φ)][1-h(φ)]c e jl(1-c e jl)+h(φ)c e j s(1-c e j s)dφφ(1-φ)(17)ε0=6λσ2.2(18)W=6.6σλ(19)式中:c e jl和c e j s是溶质在平衡状态下的浓度分布.为了模拟实际凝固过程中界面处的随机起伏,计算时需要加入一种扰动,本文在相场方程中加入扰动:f5φ5t=5φ5t+16g(φ)χ (20)式中:χ是-1到1之间的随机数, 是与时间有关的相扰动强度参数.2.2 材料物性参数为了方便计算,假设在三元合金系统中的平衡分配系数和二元合金中的平衡分配系数相同[6].表1给出了计算中要使用的材料物性参数,其中C的摩尔分数为0.5%,P的摩尔分数为0.02%.表1 Fe2C2P合金的物性参数.T ab.1 Physical parameters of F e2C2P alloy物性参数C P Fe平衡分配系数0.2040.102—液相线斜率/(K・mol-1)18021836—扩散系数(液)/(m2・s-1)2×10-8 1.7×10-9—扩散系数(固)/(m2・s-1)6×10-9 5.5×10-11—界面能/(J・m-2)——0.204熔点T m/K——1810摩尔体积/(m3・mol-1)——7.7×10-62.3 初始条件和边界条件假设初始晶核半径为R,则 x2+y2≤R2时, φ=1,T=T m-ΔT(21) x2+y2>R2时, φ=0,T=T m-ΔT(22)式中:x、y分别是横坐标和纵坐标,T是有量纲温度,ΔT是过冷度.在计算区域的边界上,φ、c采用绝热边界条件,即Zero2Neumann边界条件.2.4 数值计算方法选用上述合金Fe2C2P为研究对象,采用显示有限差分法同时求解式(10、12、13).计算的时间步长受浓度场计算的限制,即ΔT<Δx2/(4Dl)式中:D l为液相中的溶质扩散系数.枝晶轴生长主轴对应着直角坐标系的x轴和y 轴,相场和溶质场的计算网格数为800个×800个,网格的尺寸为1×10-8m,初始晶核设为R=10×10-8m的球.・2・ 兰州理工大学学报 第33卷3 结果分析3.1 过冷度对枝晶形貌和溶质分布的影响随着枝晶的长大,固、液界面向液相中推移,由于溶质的再分配导致溶质从固相中析出而富集于固、液界面前沿.图1为t =0.14ms 、ν=0.08、ΔT =16K 时的枝晶形貌和溶质分布.t 代表晶粒生长的时刻.图1a 是溶质C 的溶质浓度分布,图1b 是溶质P 的溶质浓度分布,图1c 是晶粒的相场分布图.图上的颜色标尺分别给出了图中颜色代表的溶质质量分数和固、液相的比例分数.比较图1a 和图1b ,可以看出C 的溶质扩散层比P 的扩散层要厚.这是因为C 在Fe 中的扩散系数大于P 在Fe 中的扩散系数,所以C 在界面上的扩散距离大于P 在界面上的扩散距离.从图1c 中可以看出在二次枝晶上生长出三次枝晶.在二次枝晶之间,由于凝固析出的溶质富集,造成这个区域的溶质浓度高,又因为三次枝晶的生长加长了溶质扩散的液相通道,使溶质原子更难扩散到液态区域中去.由于溶质浓度高造成合金溶液局部熔点降低,因此从图1a 和图1b 中可以看出这个区域的溶质浓度最高,在图1c 中可以看出相同的区域有很多三次枝晶在二次枝晶臂上熔断.图2给出了各向异性强度系数为0.03时,在不同过冷度条件下凝固的枝晶形貌.其中图2a 是过冷熔体为1750K 时的枝晶形貌,图2b 和图2c 的过冷熔体的温度分别是1760、1765K.很明显,图2a 中的二次枝晶生长得很发达,不同方向的二次枝晶发生碰撞,长成为一体.对比图2b 和图2c 可以发现,随着过冷度的下降,二次枝晶越不发达.因为早期凝固的二次枝晶臂间距要在凝固后期粗化[9],在图2a 中可以看出部分二次枝晶出现了合并式的枝晶间距粗化.过冷熔体的温度对枝晶尖端的生长速度有很大的影响,根据模拟结果做出过冷熔体的温度和枝晶尖端生长速度的关系如图3所示.可以看出,随着温度的升高,枝晶的生长速度总体上是下降的,可是在下降的过程中有一个上升的部分.从式(16、17)中可以看到,当温度T 升高时,相场迁移速率参数M 会减小,可是温度T 的升高又对c e j s 和c e jl 有影响,会造成M 的增大,综合这两种因素,就会出现如图3中状况,枝晶的生长速率随着温度升高而下降的过程中有部分阶段是随温度升高枝晶生长速率也升高.在枝晶的生长过程中,当过冷熔体的温度较低而形(a )溶质C 的分布(b )溶质P 的分布(c )枝晶生长形貌图1 ΔT =16K 、ν=0.08、t =0.14ms 时的枝晶生长形貌和溶质分布Fig.1 Morpology of dendrite grow th and distribution of concentration 2f ield in case of ΔT =16K,ν=0.08.t =0.14ms(a )1750K (b )1760K (c )1765K图2 ν=0.03时在不同过冷度条件下的枝晶形貌Fig.2 Morpology of dendrite grow th due to different temperatures of undercooling of alloy with ν=0.03・3・第4期 王智平等:三元合金等温凝固过程的相场法模拟 成较大热过冷的时候,枝晶可以长入过冷区而具有较大的生长速度,但同时,枝晶的生长是液体中图3 过冷熔体和枝晶生长速度的关系Fig.3 R elationship of dendrite grow th speed to tempera 2ture of undercooling alloy melt原子向固液界面陆续堆砌的过程,低的温度会影响到原子的扩散速度进而影响枝晶的生长速度[9].3.2 各向异性强度系数对枝晶形貌的影响各向异性强度系数ν表示界面表面张力、界面厚度及界面动力学各向异性的程度.它对枝晶的产生有重要的作用,对晶粒生长的影响不能忽视.在模型中,它由枝晶的最大半径和最小半径的差与和的比来定义[10].图4给出了过冷度为20K ,各向异性强度系数为0.05、0.06、0.07时的枝晶形貌.从图4可以发现各向异性强度系数越大,二次枝晶生长得越发达.而且二次枝晶臂也逐渐变得细小,二次枝晶臂间距也相应地减小.这是因为随着ν的增大,相场迁移速率参数M 的变化增大,界面前沿将变的更加不稳定,引起分枝不断地分裂所致.(a )ν=0.05(b )ν=0.06(c )ν=0.07图4 不同各向异性强度系数时的枝晶形貌Fig.4 Morpology of dendrite grow th due to different anisotropy coeff icient of alloy4 结论在各向异性强度系数较大的条件下凝固,二次枝晶上会长出细小的三次枝晶.由于溶质的大量富集,造成那个区域的熔点降低,细小的三次枝晶会发生重熔、缩颈等现象,造成三次枝晶从二次枝晶臂上熔断.过冷熔体的温度对枝晶尖端生长速度的影响并不是单调的,在一定的温度区间里,生长速度是随着温度的降低而减小的.随着温度的降低,二次枝晶生长得越发达.各向异性强度系数会通过改变相场迁移速率参数来影响枝晶的生长形貌,它的变化可以影响到二次枝晶的生长形貌和二次枝晶臂间距.参考文献:[1] 刘小刚,王承志,莫春立.Al 2Cu 合金等温凝固的相场法模拟[D ].沈阳:沈阳工业学院,2003.[2] 路 阳,王 帆,朱昌盛,等.基于KKS 模型二元合金等温凝固过程的相场法模拟[J ].兰州理工大学学报,2006,32(1):22225.[3] TOSHIO S ,MACHIKO O ,KIM S G ,et al.Phase 2field modelof dendritic growt h [J ].Journal of Crystal Growt h ,2002(2372239):1252131.[4] KIM S G ,KIM W T ,SU ZU KIA T.Phase 2field model for bina 2ry alloys [J ].Physical Review E ,1999,60(6):718627197.[5] KIM S G ,KIM W T ,SU ZU KIA T.Interfacial compositions ofsolid and liquid in a phase 2field model wit h finite interface t hickness for isot hermal solidification in binary alloys [J ].Physical Review E ,1998,58(3):331623323.[6] ODE M ,L EE J S ,KIM W T ,et al .Phase 2field model for solidi 2fication of ternary alloys [J ].ISI J International ,2000,40(9):8702876.[7] 朱昌盛,王智平,荆 涛,等.二元合金非等温凝固相场法模拟[J ].稀有金属材料工程,2005,34(10):156521568.[8] 朱昌盛,肖荣振,王智平,等.定向凝固微观组织的相场法模拟研究进展[J ].铸造,2006,55(1):43246.[9] 李庆春.铸件形成理论基础[M ].哈尔滨:哈尔滨工业大学出版社.1980.[10] KOBA YASSI H ,MACHIKO O ,KIM S G ,et al .Phase 2fieldmodel for solidification of ternary alloys coupled wit h t hermo 2dynamic database [J ].Scripta Mater ,2003,48:6892694.・4・ 兰州理工大学学报 第33卷。
材料科学-相场模拟简介
编辑ppt
t=1
t=10
t=50
14
相场模拟磁畴生成及畴界演化
相场方法模拟二级相变
编辑ppt
t=55
t=65
t=100
15
相场模拟磁畴生成及畴界演化
相场方法模拟二级相变
编辑ppt
t=300
t=800
t=10000
16
相场模拟磁畴生成及畴界演化
相场方法模拟凝固过程
相场方法模拟凝固现象可得到一般的微观组织形貌演化过程； 相场方法在凝固模拟中的应用包括纯物质的凝固，合金凝固，
编辑ppt
相场方法模拟调幅分解
Fe-Mo合金的调幅分解，进一步耦合弹性应力场
忽略弹性应力 的组织演化
考虑弹性应力 的组织演化
弹性应力对CMo=0.5合金调幅分解的影响（ T=500℃ ）
(a)t=5000；(b)t=10000；(c)t=20000；(d)t=50000
24
编辑ppt
相场模拟的进展
定向凝固过程等等； 除相场动力学方程，还需要考虑传热方程，传质方程和流体力
学方程，以及各项异性问题；
编辑ppt
相场法模拟等轴晶生长 （自适应有限元法）
相场法模拟树枝晶生长 17
编辑ppt
相场方法模拟凝固过程
要从固相和液相的过渡态找出一个序参量作为 过渡态的表征·······
纯物质固液相的区别
18
来推算相图； 在许多实际问题中得到应用。
编辑ppt
相场法模拟纤维状共晶合金凝固[2]
21
[2] M.Apel et al. Journal of Crystal Growth [J] 237-239 2002:154-158.
基于相场模型及涡量流函数形式的一种多相流数值模拟方法
基于相场模型及涡量－流函数形式的一种多相流数值模拟方法1）黄军杰2），王时龙重庆大学机械传动国家重点实验室，重庆400044摘要：本文提出一种基于相场（Phase-Field）模型及涡量－流函数形式的二维多相流模拟数值方法。
基于相场的多相流数值模拟需求解两组方程：流动控制方程（具体为不可压Navier-Stokes方程）及界面演化方程（这里使用Cahn-Hilliard方程）。
与常见方法不同的是，对于流动控制方程，本文采用其涡量－流函数形式，并且给出涡量－流函数形式下包含界面张力作用的涡量演化方程及适当的边界条件。
两组方程都使用有限差分法进行空间离散，采用四阶龙格库塔（Runge-Kutta）法进行时间推进。
另外本方法采用空间交错网格，涡量和流函数定义于网格节点，而相场变量（包括相序参数和化学势）定义于网格中心。
对于相场变量的空间导数，本文尝试使用了一般二阶中心差分以及各向同性的九点差分格式（借鉴格子玻尔兹曼方法（Lattice-Boltzmann Method, LBM））。
通过三个基本的多相流算例（平整界面，静止液滴的表面张力平衡以及接触角），本方法得以初步验证；此外，以一种格子玻尔兹曼方法作为参考，本文亦对一般中心差分以及各向同性的差分格式作以比较。
本文认为涡量－流函数形式的多相流模拟方法有其一定的优势，可作为现有常见基于压力-速度形式方法的一种替代来研究某些多相流问题（特别是二维和轴对称问题）。
关键词：多相流；数值模拟；相场模型；涡量－流函数引言很多工业问题中（如石油，化工，食品，化妆品及制药等）都涉及到多相流。
对于不可混合的液－液两相流系统的研究在理论和工程应用中都有重要意义。
基于相场模型的多相流模拟方法近年来发展迅速，颇受关注。
这类方法基于流体临界点附近的理论，使用狭窄但具有有限厚度的区域来代表界面，可以更自然的处理界面的拓扑变化（如液滴融合和分离）[1,2,3,4]。
现有相场多相流模拟方法大多采用基于速度－压力形式的流体控制方程[3,4]。
烧结过程相场法模拟
pq
f0 f (Cm ) f (C) f (C,i ) f (C,i ) f (C, )
f
(ik
,
k j
)
i1
i1
k i1 j1
p
q
f (i , ) f (i , )
i1
i1
晶粒处局部最小值
1.06 107 exp( 24000) RT
Table 2. Parameters used in the simulation
Parameters Interface energy, Initial alloy composition Initial temperature Initial cooling rate Molar volume, Vm Grid size, x Thermal conductivity
感谢感各谢位各专位家专老家师老的师指的导指！导！
,C(r))
kc 2
(C(r))2
k1 2
n1 i1
(i11)2
k2 2
n2 i2
(i22 )2 ...
kn 2
nn in
(inn )2 ]dV
f (C) n Ai i1 4
n
(C Cj )4
ji
f
(C
,ijj
)
Bj 2
n ji
(C
、
、
3
等为唯象系数
控制方程
di
dt
p ji
Lij
(
F
i
F
j
)
p j 1
T型微通道内气液两相流数值模拟
T型微通道内气液两相流数值模拟王琳琳;李国君;田辉;叶阳辉【摘要】采用相场方法,在湿壁面条件下对T型微通道内不可压缩气液两相流动进行了模拟研究.数值模拟中得到微通道内特有气泡——Taylor气泡形成的过程,发现其形成共经历了4个阶段,气泡在形成过程中主要受到挤压力、表面张力、黏性力的作用,分析各阶段这3种力对气泡的作用,得到Taylor气泡形成机理.计算结果表明,挤压力和表面张力在气泡整个形成过程都起作用,表面张力在气泡脱离前达到最小值,黏性力仅在气泡形成前两阶段起作用并在阻塞阶段取得最大值.这些基本规律为有效控制微通道内气泡尺寸和微通道系统设计提供了一定的依据.%The numerical simulation was performed to investigate the incompressible gas-liquid two-phase flow in a T-junction micro-channel with wetted walls using a phase field method. The formation process including four stages of a well defined Taylor bubble was achieved. The Taylor bubble is exerted by three forces in formation process, I. e. , extrusion force, surface tension force and shear stress force of the continuous phase, and the formation mechanism of the Taylor bubble was obtained. The numerical results show that both the extrusion force and the surface tension force act in the formation process of the Taylor bubble, and the surface tension force reaches a minimum value before detachment The shear stress force acts only in the preceding two stages of the formation process, and it has a maximum value in the blocking stage【期刊名称】《西安交通大学学报》【年(卷),期】2011(045)009【总页数】5页(P65-69)【关键词】微通道;两相流;数值模拟;相场法;气泡【作者】王琳琳;李国君;田辉;叶阳辉【作者单位】西安交通大学能源与动力工程学院,710049,西安;西安交通大学能源与动力工程学院,710049,西安;西安交通大学能源与动力工程学院,710049,西安;西安交通大学能源与动力工程学院,710049,西安【正文语种】中文【中图分类】TQ021.120世纪 80年代,随着微电子机械系统(M EM S)的提出,MEMS内微通道两相流引起了学者的广泛关注.微通道的特征尺度一般小于1mm,管道内的流动几乎可忽略重力作用,因此微通道内的流动性质和传热特性可“等价于”微重力条件下的情况[1].另外,微型气泡或液滴混合[2]、纳米级粒子的合成[3]、蛋白质结晶[4]、生物鉴定[5]和 DNA 分析[6]等都属于微通道两相流动,使微通道内两相流动问题成为研究的热点.常见的两相混合微通道分为T型、Y型和十字型等,其中T型微通道使用率较高.在T型微通道的两个入口处分别注入气体和液体,当毛细数C a小于10-2时,可规律地产生大小和间距一致且不聚合的Taylor气泡[7],这种气泡流能够降低气泡的轴向混合,加速液柱的循环流动,提高换热和传质能力,减少物质损失[8].Tay lor气泡在形成过程中会受到通道形状、大小的限制,它是微通道内特有的气泡类型,与常规尺度通道内气泡运动相比,微通道内Tay lor气泡在形成过程中阻碍液相流动,使液相压强增大,液相中持续增加的压强是趋使气泡在颈部断裂的主要原因,同时Tay lor气泡还受表面张力、黏性力的影响,这些都使微通道内流动复杂性增加.Taylor气泡的长度影响换热率,许多应用领域要求精确控制微通道内气泡的长度、体积和产生率.因此,研究Tay lor气泡的形成机理,阐明影响因素,对有效控制气泡尺寸等方面具有重要意义.由于微通道尺度极小,通过实验精确观察或测量气泡的变形、脱离,以及气泡脱离后的长度、体积和通道内流场与压强分布难度大[7].随着计算机技术的迅速提高,采用数值模拟研究微通道内的流动成为主要手段,并取得了较大的进展.目前,研究两相流界面运动常用的数值方法有水平集法(LSM)、流体体积函数法(VOF)、标记粒子单元法(MAC)、相场法(PFM)等,这些方法为深入研究微通道内两相流动提供了有效和便利的途径.基于流体自由能量模型的相场方法借助Cahn-Hilliard方程,利用化学势描述两相界面,以分离两种不同的流场,在相场模型中采用标准对流技巧,使得它在使用非结构化网格或有限元技巧模拟中相比其他方法更容易实施,并能够模拟能量耗散流动[9].本文采用相场方法,借助质量守恒方程、N-S方程和Cahn-Hilliard方程,对T型微通道内密度比接近1∶1 000的气水两相流动进行数值模拟,得到Taylor气泡的形成过程,通过分析 Tay lor气泡产生过程中的受力情况和通道内压强分布,总结出微通道内流动特性和气泡产生机理,为微通道气液两相流系统的设计和气泡尺寸控制提供参考.1 物理模型及数值方程1.1 物理模型本文以 T型微通道作为物理模型,如图1所示.气体通道延伸进水平通道的部分称为气液混合区.空气和水分别从管道的上部和左侧进入,在混合区混合,从出口流出,通道直径D为111μm,空气通道长3D,水通道长6D,主管道长28D,混合区宽度为111μm.水和空气的物理性质见表1.图1 T型微通道示意图表1 水和空气的物理性质流体密度/kg·m-3黏性系数/Pa·s表面张力系数/N·m-1水 998 0.001 2空气 1.2 0.000 018 0.072 81.2 控制方程和边界条件1.2.1 控制方程由于水和空气的入口速度不高,Re小,故认为微通道内流动为不可压缩黏性层流流动.连续方程和动量方程分别为式中:U是速度向量;p是压强;密度ρ和动力黏性系数μ由以下方程得到其中下标1、2分别表示空气和水.流体体积分数V1为归一化的相场变量,在-1到1之间变化.空气、水对应的φ分别是-1,1,因此气液界面就是φ从-1变化到 1的区域.φ通过Cahn-Hilliard方程求解式中:相场辅助变量ψ =-▽·ε2▽φ+(φ2-1)φ;γ是迁移率;λ是混合能量密度;ε是两相界面厚度,通常设ε=h2c,γ=ε2,h c为网格的特征尺寸. λ满足以下方程式中:σ是表面张力系数.表面张力Fσ通过以下方程求得式中是化学势.本文工作针对二维微通道两相流进行,计算中采用三角形网格.经验证,当网格数大于11 000时,计算结果与网格数无关,故计算网格数取32 472,时间步长取0.000 2 s.1.2.2 边界条件空气和水的入口速度分别为0.04、0.042m/s,方向垂直于入口边界;出口压强为0 Pa;水和壁面接触会浸润壁面,在气液固交界处形成接触角θ,其变化范围为0～π,θ=0表示壁面与水完全浸润,θ=π表示壁面与水完全不浸润.气液界面单位法向量式中:n w、t w分别是壁面单位法向量和单位切向量.计算中取接触角为36°.1.2.3 初始条件设初始时刻的气体通道内充满空气,整个水平通道内充满水,且微通道内流体都处于静止状态.2 数值结果与机理分析2.1 气泡形成过程气相在气液混合区拐角处断裂,形成的气泡称为Taylor气泡.T型微通道通常分为对流接触通道和错流接触通道两类.对于对流接触的T型微通道,文献[10]将其内 Tay lor气泡的形成分为3个阶段.本文模拟错流接触的T型微通道,图2给出此通道内Taylor气泡的形成过程.通过分析可见,Tay lor气泡在通道内的形成可归纳成4个阶段:(1)气泡进入气液混合区阶段(t=0～0.002 s),随着进气量增加,气泡前端进入气液混合区;(2)阻塞阶段(t=0.002 2～0.004 s),继续增大的气泡几乎阻塞气液混合区,只有少量水从气泡底部流过;(3)塌陷阶段(t=0.004 2～0.008 s),气泡底部和壁面接触,气泡主体部分向主通道下游运动,气泡后端的部分流体向气体通道内运动,气泡被逐渐拉长,在气泡后端出现颈部,随着气泡向下游运动,颈部变细;(4)脱离阶段(t=0.008 2 s),气泡在颈部处断裂,此刻一部分气体收缩回气体通道,另一部分气体收缩回主通道,在主通道内形成两端呈圆形的Taylor气泡.图2 微通道内气泡形成过程2.2 气泡表面作用力的变化特性及机理分析Tay lor气泡在形成过程中主要受到3种力:液体的挤压力、黏性力、表面张力[11],其中挤压力和黏性力对气泡表面起破坏作用,而表面张力起维持作用.2.2.1 压力变化特性及机理分析水的挤压力对T型微通道中Tay lor气泡的形成起到非常重要的作用.下面分析水对气泡作用力的变化规律,水平通道中轴线上的压强分布如图3所示.图3 气泡形成过程中水平通道中轴线上的压强分布由图3可见:在气泡形成的第1阶段,水和气泡内压强逐渐增大,水中压强从 300 Pa 增大到500 Pa,气泡内压强从400 Pa增大到1.55 kPa;在阻塞阶段,水中压强加速增大到1 kPa,气泡内压强出现最大值,约为1.68 kPa,比第1阶段小幅增长,压强最大的区域变宽,说明气泡沿水平通道向下游增长;在塌陷阶段,水中压强较第2阶段继续增大到1.4 kPa左右,气泡内压强开始减小,最终下降到1.43 kPa,这一阶段内气液压强变化小,气液界面附近压差减小,气泡后端开始向下游移动,并且气泡沿水平通道继续向下游增长;在脱离阶段,气泡从颈部脱离,气泡内压强保持脱离前大小,与塌陷阶段相同,水中压强回落至接近第1阶段大小,气泡后端较其前端压强偏高,气泡内外压差约为1.08 kPa.气泡脱离前,液体通道内压强达到最大值,气泡后端内外压差达到最小值.水平通道中轴线上压强变化的原因是:随着气泡前端进入气液混合区,气泡开始阻挡水的流动,水中压强增大,在表面张力作用下,使气泡内外存在压差,出现压强跳跃区;随气泡体积继续增大,气泡对水的流动阻挡作用增强,使水中压强迅速增大;最终增大的压强能够推动气泡主体部分向主通道下游运动,小部分向气体通道内运动,在气泡后端形成颈部,气泡后端表面的曲率半径增大,这样增大了气泡的表面自由能,但消耗了气泡内的压能,使气泡内压强降低,随着气泡向下游运动,缓解了气泡对水的阻塞作用,使水中压强增长缓慢,气泡颈部在水的挤压下变细;最后,气泡在颈部断裂,并在表面张力作用下两端收缩为圆形,气泡内压强保持脱离前的大小,此时水不再受气泡的阻塞,水中压强骤降,回落至初始时刻大小.为考察气泡在产生过程中气体和液体通道内压强的变化情况,在这两通道内分别取线段A-A、B-B(两线段与T型交汇位置的距离均为D,见图1),这两条线段上的压强变化如图4所示.对比图4a、4b中气液通道内的压强变化,可见气体通道内的压强在1.44～1.67 kPa之间变动,压强波动幅度小,与文献[12]中的结论相同.观察图2,气泡始终向水中凸出,即曲率半径 r ＜0,由Young-Lap lace方程Δ P=P g-P l=-σ知,气体r通道内的压强始终高于液体通道内的压强,这与图4中的计算结果一致,克服了气体通道内压强忽高忽低与液体通道内压强的不足.液体通道内压强波动幅度大,与文献[10]结论一致,约为1.2 kPa.2.2.2 气泡表面黏性力作用特性及分析黏性剪切应力与速度梯度成正比,对气泡表面变形起重要作用.气泡在形成过程中的速度分布如图5所示,可看到在气泡形成的第1、2阶段,气泡顶部水的速度方向与壁面平行,由于流道变窄,水的速度增大;在塌陷阶段以后,气泡表面和壁面接触,气泡后端水的速度方向与气泡表面垂直.根据黏性剪切应力与速度梯度之间的关系,后两阶段切应力可忽略,因此气泡表面切应力持续时间约占气泡形成总时间的37.5%.图4 气液通道内的压强变化图5 气泡形成过程中的速度分布为比较前两阶段气泡表面所受切应力的大小,图6给出该两阶段气泡表面到下壁面的速度分布.由此图可知,速度梯度随时间的延长而增大,对应的切应力由5.23 Pa增大到17.8 Pa,最大值在阻塞阶段取得.图6 气泡顶部液体速度分布利用文献[11]中的切应力估算公式τ=μu gap/ε0(其中u gap=Q w/ε0,ε0是气泡表面到壁面的距离,Q w是水的入口流量),可估算得到前两阶段气泡附近剪切应力从1.25 Pa增大到了28.89 Pa.比较前面的计算结果,发现利用切应力估算式得到的数据量级与模拟计算得到的数据量级相同,这说明数值计算结果是合理的.2.2.3 表面张力作用特性及分析表面张力对气泡表面起维持作用,观察图2可知,气泡后端曲率半径逐渐增大,根据Young-Laplace方程知气泡后端内外压差随时间的延长而减小.利用气泡表面的曲率半径,可估算出气泡受到的表面张力.气泡形成的前两阶段,气泡曲率半径接近气体通道的半径,表面张力约为1 324 Pa;在塌陷阶段,气泡后端出现颈部,其后端曲率半径继续增大,脱离前达到最大值,此刻后端所受表面张力最小.3 结论本文采用相场方法,对T型微通道中特有的气泡——Tay lor气泡在水中的形成过程进行数值模拟,得到气泡形成中经历的4个阶段.通过分析对气泡表面起主要作用的挤压力、表面张力和黏性力,发现微通道内流体挤压力和表面张力在气泡形成的整个过程都起作用,水对气泡的挤压力随时间的延长而增大,在气泡脱离前达到最大值;气泡后端所受表面张力随时间递减,塌陷阶段达到最小值;黏性力仅在气泡形成前两阶段起作用,较前两种力小2～3个量级,在气泡形成后两阶段可忽略不计,水的流动受到扩张气泡阻碍,产生逐渐增大的挤压力是使得气泡从颈部脱离的主要原因.本文工作为设计微通道系统和控制微通道内气泡尺寸提供了新的依据.参考文献:【相关文献】[1] GALBIATIL,ANDREIN IP.Flow pattern transition for horizontal air-water flow in capillary tubes:a m icrogravity equivalent system simulation[J].International Communications in Heat and Mass Transfer,1994,21(4):461-468.[2] TICE D,SONG H,LYON D,et al.Formation of drop lets and m ixing in mu ltiphasem icrofluidics at low values of the Reyno lds and the capillarynumbers[J].Langmuir,2003,19(22):9127-9133.[3] GUNTHER A,KHANS A,THALMANN M,et al.Transport and reac tion in m icroscale segmented gasliquid flow[J].Lab Chip,2004,4(4):278-286.[4] ZH ENG B,TICE D,ISM AGILOV F.A microf luidic system for screening p rotein crystallization conditions inside nanoliter drop lets with on-chip X-ray diffraction[J].M icro Total Analysis Systems,2005,2(297):145-147.[5] SRINIVASAN V,PAM ULA K,FA IR B.An integrated digitalm icro fluidic lab-on-a-chip for clinical diagnostics on human physiological fluids[J].Lab Chip,2004,4(4):310-315.[6] BURNS A,JOHNSON N,BRAHMASANDRA N,et al.An integrated nano liter DNA analysis device[J].Science,1998,282(5388):484-487.[7] 王昆,王嘉骏,顾雪萍,等.微通道内 Tay lor流的计算流体力学数值模拟研究进展[J].化工进展,2010,29(10):1806-1810.WANG Kun,WANG Jiajun,GU Xueping,et al.Progress in modeling of Taylor flow within microchannels by computational fluid dynam ics[J].Chem ical Industry and Engineering Progress,2010,29(10):1806-1810.[8] LIU H,ZHANG Y.D roplet formation in a T-shaped microfluidic junction[J].Journal o f App lied Physics,2009,106(3):1-6.[9] JACQM IN D.Calculation of tw o-phase Navier-Stokes flow s using phase-fieldmodeling[J].Journal of Computational Physics,1999,155(1):96-127.[10]DA I L,CA IW,XIN F.Numerical study on bubble formation o f a gas-liquid flow in a T-junction m ic rochannel[J].Chem ical Engineering and Technology,2009,32(12):1984-1991.[11]GARSTECK IP,FUERSTMAN J,STONE A,et al.Formation of drop lets and bubbles in am icrofluidic T-junction:scaling and mechanism of break-up[J].Lab on a Chip,2006,6(3):437-446.[12]KASH ID N,RENKEN A,K IWI-M INSKER L.CFD modelling o f liquid-liquid mu ltiphase m ic rostructured reactor:slug flow generation[J].Chemical Engineering Research and Design,2010,88(3):362-368.。
微观组织数值模拟——相场法与元胞自动机
微观组织的数值模拟——相场法与元胞自动机法相场法和元胞自动机法是材料科学与工程研究中常用的两种数值模拟方法。
相场模型是一种建立在热力学基础上，考虑有序化势与热力学驱动力的综合作用来建立相场方程描述系统演化动力学的模型。
其核心思想是引入一个或多个连续变化的序参量，用弥散界面模型代替传统的尖锐界面来描述界面。
相场法的不足是计算量巨大，可模拟的尺度较小（最大可达几十个微米）。
元胞自动机法是一种用来描述复杂系统在离散空间-时间上演化规律的数学算法。
元胞在某一时间步的状态转变由一定的演化规则来决定，并且这种转变是随时间推移对体系各元胞同步进行的。
元胞的状态受其相邻元胞状态的影响，同时也影响着相邻元胞的状态。
局部之间相互作用，相互影响，通过一定的规则变化而整合成一总体行为。
相场法相场法的起源与发展相场法PFM(Phase Field Method)的提出是针对具有十分复杂的界面结构的问题时，用经典尖锐界面模型去跟踪界面演化，会遭遇到严重的数值困难。
并且真实材料中的相界或晶界实际上并不是严格的零厚度界面，而是具有一定厚度（纳米尺度）的边界层，这层厚度控制材料相变动力学，由此引入一个序参量场Φ来区分两相（如固相和液相），通常称之为相场。
在相场中，Φ在固/液界面的一侧从一个常值逐渐过渡至界面另一侧的某一常值，将这个扩散界面层定义为界面，因此，在相场法中的固/液界面为弥散型界面。
Φ的主要目的是跟踪两相不同的热力学状态，可以不严格地将其理解为结晶程度的度量。
相场模型的想法最初由Langer(1978, 1986)提出的，Collin和Levine (1985)也引入了类似的相场模型(Phase field model)。
Caginalp(1985-1991)分析了这些相场模型，证明它们在界面层厚度趋于零时可以还原为尖锐界面的自由边界模型，这就从数学上证明了Langer 等人相场模型的有效性。
Fix(1983),Kobayashi(1991)等采用相场模型对具体凝固过程进行数值模拟。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
固相 转变
f0()fdisAB2C3D4 f0()4(f0)max122144
f
(1,2,3)
1 2
3
A i2
k1
1 4
3
B i4
k1
1C
6
3
i2
k1
3
各向异性
Anisotropy
界面能各向异性通过序参量的梯度项引入到自
两相
1表示在固相中 0表示在液相中 0 1在固液界面
多相
p相，相应的变量 在系统中任一点
k r
p
k 1， k 0, k
k 1
三、热力学势函数 (thermodynamic energy functional)
经典热力学
FF bu lF kin tF e l F fys
图5 自由能与浓度的关系
Nele Moelans.Phase field method to simulate microstructural evolution (June 2019)
1、固态相变-对称性降低
反相位结构
立方转变为四方相
(anti-phase domain structure) (cubic to tetragonal transformation)
图1（a）性能不连续
（b）性能连续
N. Moelans, B. Blanpain , P. Wollants, "An introduction to phase-field modeling of microstructure evolution", CALPHAD -Computer Coupling of Phase Diagrams and Thermochemistry, 32, 268-294, 2019
相场法数值模拟 phase-field modeling
内容
介绍 (Introduction)
相
场
相场变量(Phase-field variables)
法
数
值
相场方程(Phase field equations)
模
拟
热力学势函数(thermodynamic energy functional)
一、介绍
,k
)
ϵ
2
(xB
)2
k
k
2
(k
)2 dr
ϵ 和 k——梯度能量系数
ϵ F in (x B t,k) V f0 (x B ,k) 2 ( x B )2k 2 k( k)2 d r
均质与非均质体系 Homogeneous versus heterogeneous systems
序参量 (order parameters)
1或 1表示有序域 0表示无序域
图3 反相位结构
非保守场(non-conserved variables)
图4 立方结构转化成四方结构有三个等同取向
非保守场(non-conserved variables)
相场量(phase-fields)
假设C组分体系
n x i i n tot
n ci
i
V
xi Vm
C
i1
xi
C
1, ci
i1
ntot 1 V Vm
1
Vcidr
Vm
Vxidrni
xi 摩尔分数
ci 摩尔浓度
图2 两种不同组成区域
非保守场(non-conserved variables)
基本类型
1.连续相场法：扩散方程 驰豫方程
2.微观相场法：实际是 Cahn-Hilliard 方程的微 观离散格点形式。Khachatuyran 引入微观场, 用于描述由原子占据晶 格位置的几率作为场变 量来描述微结构变化
于志生, 刘平, 龙永强.基于Ginzburg-Landau 理论的相 场法研究进展[J].材料热处理技术,2019,37(16):94~98
• 相场模型是一种建立在热力学基础上，考 虑有序化势与热力学驱动力的综合作用来 建立相场方程描述系统演化动力学的模型。
核心思想
引入一个或多个连续变化的序参量，用弥散 界面模型代替传统的尖锐界面来描述界面
尖锐界面与弥散界面 sharp-interface versus diffuse-interface
相场法-热力学
体积自由能 界面能
弹性应变能 电磁相互作用能
(bulk free energy) (interfacial energy) (elastic strain energy)
F (xB ,k ) V f (xB ,k , xB , k )
V
f0
(xB
α相扩散到 β相的溶质 扩散方程
尖锐界面
c D2c t
c D2c t
(c,in t c,in)tDcDc r1 r1
(c,in)t(c,in)t
相场法原理
相场法是以GinzburgLandau理论为基础, 用微 分方程来体现扩散、有序 化势和热力学驱动的综合 作用, 它是建立在GinzburgLandau 唯象理论之上的 一种近代方法。
由能表达式中，如：
ϵ F V f0 (x B ,k ) 1 2 i,3 j 1ij x r B i x r B j 1 2 k ,l p 1 i,3 j 1ij k r il k r il d r
相场模型
该方法自提出 后，迅速成为 微观组织模拟
的热点优Biblioteka 点1.通过相场与温度场、溶质 场及其它外部场的耦合，能 有效地将微观与宏观尺度结 合起来。 2.由于不需要追踪晶界位置 能方便处理晶界上溶质聚集 和第二相析出问题，并能将 晶界能和晶界迁移率的各向 异性方便地考虑进去，还能 够较大程度避免点阵的各向 异性。
1.计算量巨大，可模拟的 尺度较小（最大可达几十 个微米）。 2.相场参数不容易确定。
二、相场变量(phase-field variables)
指那些满足局域守恒条件的场变量
保守场 如人们最熟悉的浓度序参量c
非保 守场
指那些不满足局域守恒条件的场变量 如长程序参量η
保守场(conserved variables) 成分变量