技术领域
[0001] 本
发明属于金属凝固过程模拟领域，具体涉及一种BGA锡铅焊球凝固过程模拟的方法。
背景技术
[0002] 数值模拟在理论和实验研究间架起了一座
桥梁，极大地促进了科技进步。元胞自动机模拟作为一种通用的模拟手段日益受到广泛关注。同时，随着集成
电路特征尺寸的缩小，球栅阵列封装(BGA)的重要性日益凸显。而锡铅精密焊球的综合性能对BGA微互连的可靠性至关重要。本发明基于元胞自动机模拟对用于BGA的微米级锡铅精密焊球凝固过程进行模拟。
[0003] 无容器凝固(containerless solidification)过程可以避免由于容器壁导致的异质形核，因而可以实现材料的过深
过冷快速凝固使其具有组织更均匀、结构更细化以及强度更高等传统
铸造过程难以实现的优异性能。利用无容器快速凝固方法制备高性能的
合金材料已经引起了材料研究者的极大兴趣，并被认为是开发新材料的重要手段之一。均匀液滴喷射(UDS)法制备无铅精密焊球是典型的无容器快速凝固过程。该方法制备的锡铅焊球在尺寸均一性、表明光洁性以及内部微结构等方面具有优良特性。
[0004] 获得高性能的材料是材料学科研究的主要目的。在材料制备过程中，几乎所有的金属材料都经历了一次或多次的凝固过程。在凝固时形成的显微组织会直接影响金属后续的热加工和最终性能。正确理解凝固过程中显微组织的演化是控制凝固组织和铸件性能的关键。枝晶是一种最常见的金属凝固显微组织。对于具有一定凝固区间的合金材料，在凝固过程中会形成一个液/固相共存的糊状区。在大多数中低速的凝固过程中，枝晶组织在糊状区发生粗化是不可避免的。枝晶粗化会对最终的凝固组织、微观偏析和产品性能产生重要影响。例如，大量实验数据表明，二次枝晶臂间距(SDAS)与材料的
力学性能直接相关。
[0005] 此外，凝固也是自然界中常见并且涉及范围很广的一类
相变过程，作为凝固过程中最常被观察到的一种微观组织，枝晶则代表着自然界中一种典型的非线性、自组织模式。因此，对于凝固及其组织的研究不仅具有很强的工程应用背景，而且具有重要的理论研究价值。然而，凝固过程涉及到一系列复杂的物理化学变化。从宏观上看，凝固伴随着
传热、传质、流动；从微观上看，凝固则体现为形核、生长的动力学过程，而宏微观之间又有着复杂的相互影响。这使得理论研究和数学解析遇到了很大的困难。而且，由于发生凝固的液态金属往往处于较高的
温度下，金属本身又不透明，并处于流动状态，这使得实验研究也缺乏直接和有效的方法。
[0006] 随着计算机
硬件和数值计算技术的迅速发展，当今数值模拟已成为与实验技术和理论研究并行发展的第三种科学研究方法。它对材料科学和技术的发展起到了日益重要并且不可或缺的作用。显微组织的数值模拟可以将相变
热力学、动力学、界面能、晶体生长的
各向异性、温度场、浓度场、流场、
应力-应变场等诸多因素有效地综合起来，能够
可视化地再现合金在相变过程中显微组织和溶质偏析的演化。通过模拟研究，可以帮助人们深入理解凝固过程中显微组织形成的内在机制，为材料成分设计、改进现有工艺以及发展新工艺提供科学依据。
发明内容
[0007] 本发明的目的是为解决目前在使用元胞自动机方法模拟金属凝固过程时多采用Kurz-Givoanola-Trivedi(KGT)等平衡条件解析模型来计算界面的生长速度和过冷度之间的关系，并未考虑在非等温、非等溶质界面情况下需要使用非平衡界面系列模型来模拟过冷合金熔体快速凝固的过程，导致模拟过程不符合BGA锡铅焊球实际凝固过程。
[0008] 本发明为解决上述技术问题采取的技术方案是：
[0009] 一种基于元胞自动机的BGA锡铅焊球凝固过程模拟方法，所述方法的实现过程为：
[0010] 步骤1：确定模拟过程所需的物理性质参数：热力学参数、溶质扩散参数和形核参数，
[0011] 热力学参数包括热过冷度、合金熔点、元胞温度、
潜热、合金材料的
比热容、热扩散系数，溶质扩散参数包括液相线斜率、溶质平衡分配系数、液相溶质扩散系数，形核参数包括平均界面动力学系数、吉布斯-汤姆森系数、动力学各向异性强度；
[0012] 步骤2：将固定的二维计算区域均匀划分为正方形网格，每个正方形网格代表一个单位元胞，并且建立以i列j行的二维
坐标系，i，j为每个元胞的坐标参数；
[0013] 步骤3：初始化单位元胞，为每个单位元胞的状态参数、固相率、液相溶质浓度和固相溶质浓度赋予初始值，并通过预设时间尺度决定循环步数；其中单元胞的状态state包括三个，即液相胞L、界面胞M和固相胞S，分别对应三个状态参数值0、1和2；
[0014] 步骤4：建立枝晶生长模型，即计算单个元胞的固相率增长△fs及其剩余液相溶质浓度Cl；
[0015] 步骤5：对液相胞进行捕获，捕获方为基于冯诺依曼捕获模型，对每一个固相元胞周围上下左右四个单位液相元胞进行判定，若符合捕获规则，则将液相元胞设置为界面胞；
[0016] 步骤6：根据步骤4建立的枝晶生长模型，计算被捕获的临近界面胞M的固相增长率△fs、剩余的液相溶质浓度Cl和单位元胞的温度T；
[0017] 步骤7：被捕获的临近元胞计算完毕后，判断被捕获的临近界面胞M在此刻的单位时间步长内固相率加上固相增长率fs+△fs是否大于或等于1，若大于1，则将fs置为1，将此界面元胞的状态state置为固相胞，并执行步骤8；若否，执行步骤6；每个固相胞的四邻胞捕获过程同时、单独进行，互不影响；
[0018] 步骤8：计算固相溶质扩散，不考虑固相的溶质扩散；计算单位元胞凝固过程产生的温度扩散；
[0019] 步骤9：判断设置的i列j行坐标系内的所有元胞是否全部遍历结束，如果这一时间步长内全部元胞遍历完成，更新计算下一个时间步长，在下一个时间步长内重复步骤4至步骤8；
[0020] 步骤10：在每次计算得到的单位时间步长进行累加，若达到初始设置的总时间，计算模拟过程结束。
[0021] 所述步骤4中的枝晶生长模型具体为：
[0022] 单个元胞的固相增长率计算公式如下所示，
[0023]
[0024] 式中，S(Vn,△t)为元胞内因界面移动增加的固相面积，(△x)2为元胞面积，Vn为固液界面处的法向速度，下
角标处n表示界面处的法向方向，上角标处0或者n代表0或n时刻；
[0025] 0时刻初始生长速度Vn0先给一个很小的初始数值，得到初始计算时的△fs，得到界面胞凝固过程产生的热△T为：
[0026]
[0027] 式中，ΔHf为凝固潜热释放的热量，cp为合金材料的比
热容；
[0028] 界面胞凝固过程产生热量使自身升温后，再由热扩散计算模型将热传递给周围上下左右四个元胞，热扩散方程为：
[0029]
[0030] 式中，DT、DL分别为热扩散系数和液相溶质扩散系数；
[0031]
[0032] 根据上式则可以得到CL*的数值，由CL*可得出下一时刻的固液界面速度[0033]
[0034] 随着程序的运行，依据每次运行步数的推进，依次可得到程序运行n时刻之后的固液界面速度Vnn；
[0035] 式中，Tm为合金熔点，Tl*为元胞温度，ΓKη(θ)为
曲率过冷， 为动力学过冷，μ0为界面动力学系数， 为成分过冷；
[0036] 式中，Γ为吉布斯-汤普森系数，C0为元胞初始液相成分， 为当前时刻元胞固液界面处成分，mL为液相线斜率， 固液界面速度为Vn时的液相线斜率，表示为：
[0037]
[0038] 式中，
[0039] 其中，
[0040] 式中，ke为溶质平衡分配系数，k(Vn)为固液界面速度为Vn时的溶质瞬时分配系数，VDI为界面溶质扩散速度；
[0041] η(θ)为界面能各向异性函数，表示为：
[0042] η(θ)＝1-εcos[4(θ-θ0)]；
[0043] 式中，ε为各向异性强度，θ0为择优生长方向，θ则由界面法向量 得到；
[0044]
[0045] 式中，fs为界面胞元的固相份数；界面处曲率可由
整理得
[0046]
[0047] 式中，fx、fy表示fs的2个一阶偏导数，fxx、fxy、fyy表示fs的3个二阶偏导数；
[0048] 界面元胞中固相份数△fs的增加需要排出的多余溶质为
[0049] ΔC＝CL(1-ke)△fs；
[0050] 这些溶质将平均分配给当前界面元胞的液相临胞；当前界面元胞的固相份数fs和成分C则分别更新为：
[0051] fs＝fs0+△fs；
[0052]
[0053]
[0054] 式中， 分别为本次循环和上次循环所计算的剩余液相溶质浓度，分别为本次循环和上次循环所计算的剩余固相溶质浓度，fs为元胞的固相率，k(Vn)为固液界面速度为Vn时的溶质瞬时分配系数。
[0055] 所述步骤8中液相的溶质扩散的计算模型和热扩散的计算模型如下：
[0056] 合金凝固过程中，枝晶的生长过程由溶质扩散和热扩散控制，在不考虑自然
对流和强制对流的情况下，液相胞的溶质扩散由下式表示：
[0057]
[0058] 界面胞凝固过程产生的热△T为：
[0059]
[0060] 式中，ΔHf为凝固潜热释放的热量，cp为合金材料的比热容；
[0061] 界面胞凝固过程产生热量使自身升温后，再由热扩散计算模型将热传递给周围上下左右四个元胞，热扩散方程为：
[0062]
[0063] 式中，DT、DL分别为热扩散系数和溶质扩散系数。
[0064] 所述步骤9中时间步长的计算方程为：
[0065]
[0066] 式中， 为固液界面处的最大速度。
[0067] 模拟过程所需的物理性质参数如表1所示：
[0068] 表1物理性质参数表
[0069]
[0070] 采用上述技术方案所产生的有益效果在于：本发明提供的基于元胞自动机的BGA无铅焊球凝固过程模拟方法，发明点在于采用元胞自动机方法模拟快速非平衡条件下的BGA锡铅焊球的凝固过程。前人在进行金属凝固过程的元胞自动机模拟多采用Kurz-Givoanola-Trivedi(KGT)等平衡条件解析模型来计算界面生长速度和过冷度之间的关系。然而，沿着固液界面，界面曲率和法向速度实际上不是常数，这将导致沿着固液界面界面曲率过冷和溶质偏析的变化。这意味着界面将是非等温、非等溶质的。利用本发明的方法可研究等轴晶的生长过程及其周围的溶质偏析情况、研究不同过冷度对等轴晶生长的影响、计算柱状晶的生长情况的同时，本发明还考虑了热扩散对BGA锡铅焊球凝固过程的影响。
附图说明
[0071] 图1为本发明
实施例提供的基于元胞自动机的BGA锡铅焊球凝固过程模拟方法
流程图；
[0072] 图2为本发明实施例提供的模型网格划分示意图；
[0073] 图3为本发明实施例的元胞捕获模型示意图；
[0074] 图4为本发明实施例提供的等轴晶生长图；
[0075] 图5为本发明实施例提供的柱状晶生长图；
[0076] 图6为前人在
显微镜下Al-Si合金铸态下的枝晶α组织；
[0077] 图7为本发明实施例提供的施加随机扰动后枝晶分枝的竞争生长及成分分布。
具体实施方式
[0078] 下面结合附图对本发明的技术方案作进一步的说明，但并不局限于此，凡是对本发明技术方案进行
修改或者等同替换，而不脱离本发明技术方案的精神和范围，均应涵盖在本发明的保护范围中。
[0079] 一种基于元胞自动机的BGA锡铅焊球凝固过程模拟方法流程如图1所示，本实施例的方法如下所述。
[0080] 步骤1：从实验及枝晶生长相关文献获取计算所需的热力学参数、溶质扩散参数和形核参数，其中热力学参数包括热过冷度、合金熔点、元胞温度、潜热、合金材料的比热容、热扩散系数，溶质扩散参数包括液相线斜率、溶质平衡分配系数、液相溶质扩散系数，形核参数包括平均界面动力学系数、吉布斯-汤姆森系数、动力学各向异性强度。具体数据如表1所示。
[0081] 表1物理性质参数表
[0082]
[0083]
[0084] 步骤2：将固定的二维计算区域均匀划分为正方形网格，每个正方形网格代表一个单位元胞，并且建立以i列j行的二维坐标系，i，j为每个元胞的坐标参数；
[0085] 本实施例选用的合金为锡铅合金，原始液相溶质浓度C0＝2.86，每个网格的宽度Δx＝2μm，计算区域的大小为直径750μm的圆形区域，划分为375×375个网格，如图2所示，为保证扩散方程的
稳定性，需要对时间步长△t进行控制，即公式：
[0086]
[0087] 本实施例的方法在Visual studio 2012平台上进行自主编程、调试运算。
[0088] 步骤3：初始化单位元胞，为每个单位元胞的状态参数、固相率、液相溶质浓度和固相溶质浓度赋予初始值，并通过预设时间尺度决定循环步数；其中单元胞的状态state包括三个，即液相胞L、固相胞S和界面胞M，分别对应三个状态参数值0、1和2；初始时刻设置中心(375/2,375/2)元胞为固相胞，即为1，其余都为液相胞，即为0；对应中心(375/2,375/2)元胞的固相率为1，其余为0；中心元胞(375/2,375/2)的液相溶质浓度为0，其余都为2.86；中心元胞(375/2,375/2)的固相溶质浓度为0.4862，其余为0。本实施例中预设运行总时间为0.3秒。
[0089] 步骤4：建立枝晶生长模型，即计算单个元胞的固相率增长Δfs及其剩余液相溶质浓度Cl；
[0090] 单个元胞的固相增长率计算公式如下所示，
[0091]
[0092] 式中，S(Vn,△t)为元胞内因界面移动增加的固相面积，(△x)2为元胞面积，Vn为固液界面处的法向速度，下角标处n表示界面处的法向方向，上角标处0或者n代表0或n时刻，△t为时间步长。
[0093] 0时刻初始生长速度Vn0先给一个很小的初始数值，得到初始计算时的△fs，得到界面胞凝固过程产生的热△T为：
[0094]
[0095] 式中，ΔHf为凝固潜热释放的热量，cp为合金材料的比热容。
[0096] 界面胞凝固过程产生热量使自身升温后，再由热扩散计算模型将热传递给周围上下左右四个元胞，热扩散方程为：
[0097]
[0098] 式中，DT、DL分别为热扩散系数和液相溶质扩散系数。
[0099]
[0100] 根据上式则可以得到CL*的数值，由CL*可得出下一时刻的固液界面速度[0101]
[0102] 随着程序的运行，依据每次运行步数的推进，依次可得到程序运行n时刻之后的固n液界面速度Vn。
[0103] 式中，Tm为合金熔点， 为元胞温度，ΓKη(θ)为曲率过冷， 为动力学过冷，μ0为界面动力学系数，mLC0-mVnCL*为成分过冷；
[0104] 式中，Γ为吉布斯-汤普森系数，C0为元胞初始液相成分， 为当前时刻元胞固液界面处成分，mL为液相线斜率， 固液界面速度为Vn时的液相线斜率，表示为：
[0105]
[0106] 式中，
[0107] 其中，
[0108] 式中，ke为溶质平衡分配系数，k(Vn)为固液界面速度为Vn时的溶质瞬时分配系数，VDI为界面溶质扩散速度；
[0109] η(θ)为界面能各向异性函数，表示为：
[0110] η(θ)＝1-εcos[4(θ-θ0)]；
[0111] 式中，ε为各向异性强度，θ0为择优生长方向，θ则由界面法向量 得到；
[0112]
[0113] 式中，fs为界面胞元的固相份数。界面处曲率可由 整理得
[0114]
[0115] 式中，fx、fy表示fs的2个一阶偏导数，fxx、fxy、fyy表示fs的3个二阶偏导数。
[0116] 界面元胞中固相份数△fs的增加需要排出的多余溶质为
[0117]
[0118] 这些溶质将平均分配给当前界面元胞的液相临胞。当前界面元胞的固相份数fs和成分C则分别更新为
[0119] fs＝fs0+△fs；
[0120]
[0121]
[0122] 式中， 分别为本次循环和上次循环所计算的剩余液相溶质浓度，分别为本次循环和上次循环所计算的剩余固相溶质浓度，fs为元胞的固相率，k(Vn)为固液界面速度为Vn时的溶质瞬时分配系数。
[0123] 步骤5：对液相胞进行捕获，捕获方为基于冯诺依曼捕获模型，对每一个固相元胞周围上下左右四个单位液相元胞进行判定，若符合捕获规则，则将液相元胞设置为界面胞；判定的捕获规则为：
[0124]
[0125] 式中，Tm为合金熔点，Tl*为界面胞温度， 为动力学过冷，C0为元胞初始液相成分，CL为当前时刻元胞液相成分，mL为液相线斜率， 固液界面速度为Vn时的液相线斜率。
[0126] 步骤6：根据步骤4建立的枝晶生长模型，计算被捕获的临近界面胞M的固相增长率△fs、剩余的液相溶质浓度Cl和单位元胞的温度T；
[0127] 步骤7：被捕获的临近元胞计算完毕后，判断被捕获的临近界面胞M在此刻的单位时间步长内固相率加上固相增长率fs+△fs是否大于或等于1，若大于1，则将fs置为1，将此界面元胞的状态state置为固相胞，并执行步骤8；若否，执行步骤6；每个固相胞的四邻胞捕获过程同时、单独进行，互不影响；
[0128] 步骤8：计算固相溶质扩散，不考虑固相的溶质扩散；计算单位元胞凝固过程产生的温度扩散。
[0129] 合金凝固过程中，枝晶的生长过程由溶质扩散和热扩散控制，在不考虑
自然对流和强制对流的情况下，液相胞的溶质扩散由下式表示：
[0130]
[0131] 界面胞凝固过程产生的热△T为：
[0132]
[0133] 式中，ΔHf为凝固潜热释放的热量，cp为合金材料的比热容。
[0134] 界面胞凝固过程产生热量使自身升温后，再由热扩散计算模型将热传递给周围上下左右四个元胞，热扩散方程为：
[0135]
[0136] 式中，DT、DL分别为热扩散系数和溶质扩散系数。
[0137] 步骤9：判断设置的i列j行坐标系内的所有元胞是否全部遍历结束，如果这一时间步长内全部元胞遍历完成，更新计算下一个时间步长，在下一个时间步长内重复步骤4至步骤8；
[0138]
[0139] 式中， 为固液界面处的最大速度。
[0140] 步骤10：在每次计算得到的单位时间步长进行累加，若达到初始设置的总时间，计算模拟过程结束。结束后，将数据导入到Excel表格内进行处理，之后用画图
软件做出微观组织的生长模型。
[0141] 本实施例提供的方法可研究等轴晶的生长过程及其周围的溶质偏析情况；可研究不同过冷度对等轴晶生长的影响；可计算柱状晶的生长情况。
[0142] 如图4所示，是过冷度为10K，择优生长方向角为0°时的等轴晶生长图。可以看出等轴晶稳定的沿着择优方向生长，并且枝晶形貌呈现四重对称。如图5所示，是过冷度为10K，择优生长方向角为0°时的柱状晶晶生长图。如图7所示，是过冷度为10K，施加振幅为0.01的随机扰动后模拟得到的枝晶形貌及成分偏析。从图中可以看出因为施加了界面扰动，枝晶已不再仅仅以光滑的主干生长，而是出现了分枝结构。这些分枝之间相互竞争生长，并且由于分枝的端部具有较好的扩散条件，因而往往端部则比较粗大，而根部则呈颈缩状。在分枝之间则形成了溶质富集区，这些富集区的液相具有较高的成分，并且由于枝晶骨架的形成而进一步阻碍了溶质的扩散，因而将成为最后凝固的部分。
[0143] 最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明
权利要求所限定的范围。
