子1外力作用之和fd
[0073] 步骤2. 4. 2 :根据式（14)及式（15)分别计算研磨颗粒S的线加速度
和角加 速虔
[0076] 其中Ms和Is分别为研磨颗粒S的质量和转动惯量，rs为研磨颗粒中心的位移矢 量。并可以根据式（16)求得研磨颗粒的速度
[0078] 步骤2. 4. 3 :通过文献[11]所述方法，采用式（17)对研磨颗粒与其他固体材料 (研磨颗粒、抛光垫、硅片）的碰撞进行模拟，获得颗粒a和b之间的相互作用力Fab :
[0080] 其中Cab为排斥力的尺度因子，通常与重力加速度同一量级；?为刚性系数，一般取 I X 10 6 ; ξ为排斥力的判定距离，一般为3 Λ X ;rab为两个研磨粒中心连线矢量，R为物体半 径。当两个固体距离小于ξ时，由Cab和?决定的强排斥力会将两个物体排开；
[0081] 当研磨粒靠近抛光垫和硅片固体边界的时候，通过边界外的镜像研磨粒w提供相 应的排斥力，关系式如下
[0083] 其中，raw是从研磨粒a中心指向边界外镜像研磨粒中心的矢量；
[0084] 本发明步骤2. 5中，通过下述子步骤，利用SPH方法计算粒子相互作用力及运动轨 迹，对特征尺寸级CMP工艺进行动态仿真：
[0085] 步骤2. 5. 1 :根据文献[12]所述方法对所有粒子进行相邻关系查找，生成邻居粒 子列表；
[0086] 步骤2. 5. 2 :对任意时间步n，采用式（19)、（20)计算〃 + $·时间步粒子的速度
和位移
[0089] 其中At为时间步长；
[0090] 步骤2. 5. 3 :利用半时间步长速度
和位移
结合式（7)求出n+1步的密度 Pn+1，并更新第n+1步的位移一+1
[0093] 步骤2. 5. 4 :根据式（21)和（22)获得的n+1步密度和位移，求解式⑶至式（10) 获得研磨液粒子的加速度，利用式（13)至式（18)获得研磨粒粒子的加速度，并更新所有粒 子第n+1步的速度
[0095] 步骤2. 5. 5 :判断当前时刻是否等于总仿真时间T，如果不相等则进行到下一个离 散时间步，重复步骤2. 5. 1至2. 5. 4 ;否则进行步骤3 ;
[0096] 本发明步骤3中，通过下述子步骤，依据仿真数值仿真结果，通过统计方法得到 MRR，进而实现对实际工艺过程的预测：
[0097] 步骤3. 1 :将每一时刻粒子i的物理量（如位移Γι、密度P i、速度V1、压力P1、剪切 力T1等）作为仿真结果输出，并获得相应时刻的研磨液粒子压力分布 Pwga (x，y)、剪切力 分布 T研磨液 (χ, y)；
[0098] 步骤3. 2 :根据研磨液粒子的剪切力分布τ W)ga (χ, y)可以获得有研磨粒时娃片 表面任意一点k的剪切力
[0100] 硅片表面所受最大剪切力为
[0102] 当\大于无研磨颗粒时硅片表面最大剪切力Tth时，即为有效剪切Tactire。统 计τ &tlTC出现的时间tartlTC占总仿真时间ttotal的比例，获得剪切频率
根据文献
[13]中MRR与研磨粒的剪切频率FT(shearfrequenct)呈现正比关系，即：
[0103] MRR Ft (26)
[0104] 即可预测MRR的比例值。
[0105] 本发明方法采用微小粒子表征化学机械抛光中的工艺材料，用NS流体力学方程 表述粒子间相互作用，通过光滑化粒子流体动力学方法进行计算获得抛光仿真结果，经测 试结果显示，不依赖经验模型，只需根据工艺材料参数即可完成仿真，具有普适性强的优 点。
[0106] 本发明的优点在于：
[0107] 1.本发明不依赖经验模型，只需根据工艺参数即可完成仿真建模，具有普适性强 的优点；
[0108] 2.本发明方法的仿真过程模拟抛光垫、抛光液、抛光颗粒与硅片之间实际物理过 程，具有精度高的优点；
[0109] 3.利用本发明方法可以对各类工艺条件进行有效仿真，为优化化学机械抛光工艺 流程和工艺参数控制提供有效参考。
【附图说明】
[0110] 图1为一种旋转式化学机械抛光机台的示意图。
[0111] 图2为碟陷（dishing)与侵蚀（erosion)示意图。
[0112] 图3为本发明方法的流程示意图。
[0113] 图4为本发明将空间离散为粒子的示意图。
[0114] 图5为本发明依据离散粒子所处材料区域分配给对应材料示意图。
[0115] 图6为本发明实施算例1的模型示意图。
[0116] 图7为本发明实施算例1中管道上表面压力分布曲线与FLUENT比较。
[0117] 图8为本发明实施算例1中管道上表面剪切力分布曲线与FLUENT比较。
[0118] 图9为本发明实施算例2的模型示意图。
[0119] 图10为本发明实施算例2中无研磨颗粒情形、研磨颗粒未通过狭缝时、研磨颗粒 通过狭缝时管道上表面剪切力分布曲线。
[0120] 图11为本发明实施算例2中不同研磨颗粒浓度变化与有效剪切频率关系曲线。
[0121] 图12为本发明实施算例2中不同研磨颗粒浓度与有效剪切频率关系曲线拟合情 况。
【具体实施方式】
[0122] 通过具体实例的实施过程，进一步描述本发明的方法。
[0123] 实施算例1
[0124] 本发明模型计算得到的无抛光颗粒研磨液流体力学特性（流场与力场）结果与商 用流体力学计算工具FLUENT之间的比较。
[0125] 在该实例中，选取二维粗糙边墙管道流进行仿真。采用如图6所示的管道，该管道 由顶部平直的边墙和底部正弦曲线形状的边墙构成，管道的长度为L，高度H的表达式为
[0127] 管道的上下边墙存在速度为Vwall的相对运动，并拖拽管道中的流体发生运动。
[0128] 本例管道长L = 3X 10 5m，正弦边界H。= I X 10 5m，正弦边界振幅Ha = 5X 10 6m, 流体密度P = l〇〇〇kg/m3,流体动力学粘度η = 1X10 3Pa · s，边墙速度Vwall = lm/s。
[0129] 通过本发明模型获得的上边墙压力和剪切力与现有技术方法之间的比较（图7及 图8)，结果显示，本发明模型的计算结果与现有技术方法FLUENT误差很小，其中最大误差 7. 7%，平均误差1. 3%，实验数据表明本发明计算方法的正确性。值得指出的是，本算例中 并未考虑液体中的悬浮固体运动，当考虑固体（抛光颗粒）在液体（抛光液）中的运动情 况时，FLUENT等主流商用工具由于采用网格方法，无法进行求解，而本方法可以统一进行处 理，本申请将在实施算例2中说明该优点。
[0130] 实施算例2
[0131] 在本实施算例中，采用高斯曲线管道流模型仿真硅片与抛光垫界面情形，并增加 了圆形的漂浮颗粒用于模拟研磨颗粒（如图9所示）。由于固体悬浮颗粒的存在，基于网格 方法的常见商用软件如Fluent等，无法对本算例进行仿真，因此本发明在处理此类问题上 具有优势。管道高度H表达式为：
[0133] 其中X(] = 4〇表示下表面尖峰x坐标，尖峰的宽度由〇决定。
[0134] 本算例的基本参数为：上下边墙最小间距Hgap = 3Χ 10 6m、上下边墙最大间距Η_ =2Χ 10 5m、管道长度8 σ = 4Χ 10 5m、研磨粒直径IX 10 6m、流体密度P = 1000kg/m3、流 体动力学粘度Π = 1X10 3Pa · s、研磨粒密度2000kg/m3、边墙速度Vwall = lm/s。本算例 对研磨粒浓度（用质量百分浓度wt%表示）从1.07%到7. 49%进行了一系列仿真。管道 中研磨粒数目与研磨粒浓度的关系式为：
t
[0136] 其中Vfluid与Vatoasire分别表不研磨液与研磨粒的体积。
[0137] 根据仿真结果，当研磨颗粒在研磨液厚度较宽区域运动时，对流体在硅片表面的 压力和剪切力作用影响很小。然而，当这些研磨粒进入硅片与抛光垫之间的狭缝时，则会引 起流体压力和剪切力的较大变化（如图10所示）。在本领域，已发现研磨粒浓度对研磨效 果会产生影响，并发现MRR与（Wt% )1/3成正比关系，并且随着Wt%的进一步增加，MRR会达 到饱和状态[14]。此外大量针对CMP移除机理的实验研究表明MRR与研磨粒的剪切频率 (shearfrequenct)F T亦呈现正比关系[13]，即
[0138] MRR Ft (30)
[0139] 图11及图12展示了研磨粒浓度（wt% )与剪切频率（Ft)之间的关系。在这两 幅图中，采用了归一化浓度（normalized wt% )作为横坐标，其中归一化浓度通过不同的浓 度值除以最大浓度值获得。从中可以发现：
[0140] a)FT随着研磨粒浓度的增长而显著地增长；
[0141] b)当研磨粒浓度达到某具体值（约5% )之后，Ft达到饱