1. 晶圆表面颗粒物仿真，难点到底在哪

晶圆表面颗粒物生成与沉积，是半导体前道和后道都会遇到的顽固问题。一片12英寸晶圆，表面只要有几十颗亚微米级颗粒落在关键图形区，后续光刻的对准、薄膜的附着、键合的界面就可能出问题。产线上通常靠清洗、吹扫、洁净室等级控制去压这个指标，但真要回答“这些颗粒是从哪来的、为什么会落在这个位置、多大的颗粒最难去除”，光靠显微镜数颗粒是不够的。这正是EDEM这类离散元仿真工具能插进来的地方：把颗粒当成一颗一颗独立的个体去追踪，把晶圆当成几何体，把颗粒与晶圆、颗粒与颗粒之间的碰撞、粘附、滚动全部按力学模型算出来。

我用EDEM做这个入门案例，目的不是复现某条真实产线，而是搭一套最小可用的算例框架：晶圆摆在那儿，上方颗粒工厂源源不断生成颗粒，颗粒在重力（或简化气流）作用下往下走，落到晶圆表面后靠粘附力留下来，最后统计表面颗粒的数量、面密度和粒径分布。这套流程走通一遍，你后面想加静电场、加流场、加晶圆旋转或翻转，都是在同一个骨架上加模块。适合谁看？刚接触离散元、手里有软件但不知道第一个案例该建什么的同学；也适合做洁净控制、设备搬运、晶圆传片机构设计的工程师，想粗略估一下“这个动作会不会扬尘、颗粒落下来多少”。

先说一句可能让人不太舒服但必须提前说明的话：这个案例的物理尺度极其不友好。晶圆直径300 mm，污染颗粒粒径常在0.1 μm到5 μm之间，差了三到四个数量级；而DEM求解的时间步长又必须小于颗粒自身弹性波传播的特征时间。这两件事叠加，意味着“按真实尺度老老实实算一秒”在普通工作站上是不可行的。所以入门案例的核心不是把参数填得多准，而是先搞明白哪些量可以缩放、哪些量绝对不能乱动。下面我按自己实际搭这套算例的顺序，把设计思路、参数来源、操作细节和踩过的坑一层层铺开。

2. 案例搭建前的整体设计

动手点鼠标之前，得先把物理问题翻译成计算模型。这一步做不好，后面所有的参数调整都是在错的框架里打补丁。我通常会把一个DEM算例拆成四个问题：计算域多大、颗粒代表什么、时间算多长、输出什么指标。这四个问题回答完，几何、工厂、求解器设置基本就自动确定了。

2.1 建模范围与简化取舍

真实场景里，晶圆躺在一个片盒、传片手臂或者工艺腔体里，周围有气流、有机械振动、有静电。全都要，算不动；全都不要，结果没意义。我的取舍原则是：只保留对这个指标有直接贡献的因素。入门案例里，我只保留晶圆本体、一个简化包围盒（当作计算域边界，同时模拟腔体侧壁）、颗粒工厂、重力场。气流先不加，因为加了就得耦合，耦合就得引入流体网格和时间尺度匹配，入门阶段容易把自己绕进去。

晶圆本体我用直径100 mm的缩比圆盘代替真实的300 mm，厚度按同比例缩到0.3 mm左右。为什么缩比？因为颗粒数量、接触检测的开销跟几何尺寸和颗粒数量强相关。把直径砍到三分之一，同样面密度下颗粒数少了将近九成，跑起来轻松很多，而颗粒在表面的沉积行为、粘附判据、粒径分布规律不会因为晶圆变小而改变。只要在后处理里把“颗/mm²”这个归一化指标算出来，缩比模型和全尺寸模型是可比的。这里必须明确：这是基于常见仿真实践做的等效处理，不是真实尺寸复现，写报告时要把缩放关系写清楚，别让人误读成实物数据。

包围盒我做成一个开口长方体，底面就是晶圆所在平面稍微往下一点，四周壁面设成可以吸收颗粒或者直接让颗粒离开计算域，避免颗粒堆在角落里反复弹跳浪费时间。颗粒工厂放在晶圆正上方约80 mm处，做成一个薄片状的虚拟体积，让颗粒从近似一个面上往下释放，这样落点分布比较均匀，不会全都堆在圆心。

2.2 输入参数清单与来源

离散元算例最容易被质疑的就是“你这参数哪来的”。我习惯在动手前先列一张表，把每个参数、取值、来源和敏感度写清楚，跑完再回头看哪些参数其实不敏感可以粗放，哪些必须标定。

这张表里有两个值我要特别强调，因为它们直接决定了算例能不能跑、结果像不像。第一个是剪切模量。真实二氧化硅的剪切模量在30 GPa量级，如果老老实实填进去，Rayleigh时间步会小到0.5 ns左右，算1秒需要二十亿步，根本不现实。常规做法是把剪切模量人为降低一到两个数量级，时间步可以放大约3倍，而颗粒的宏观堆积、沉降、碰撞行为变化不大——这个技巧在离散元圈子里叫“模量缩放”，是被广泛接受的工程近似手段，但前提是颗粒的变形量不能大到穿透，且碰撞速度不能太高。第二个是JKR表面能，它是控制“颗粒落上去粘不粘”的总开关，后面我会专门讲怎么标定。

2.3 三招让算例跑得动

尺度矛盾摆在面前，必须做取舍。我的做法是三招叠加，而不是单靠某一招硬扛。

第一招是粒径放大。把代表污染颗粒的粒径从亚微米放大到20~50 μm。这样一来时间步能放大几十倍，颗粒数量也降下来了。代价是失去了真实小颗粒的布朗运动效应和范德华力主导的物理图像，所以这个算例只能用来研究“宏观沉降与粘附趋势”，不能用来预测真实亚微米颗粒的沉积率。这一点必须在结论里写明白，否则就是误导。

第二招是时间窗截断。不去算“颗粒从生成到稳态沉积几个小时”，只算“释放后0.5~1秒内颗粒的落点分布”。在这个时间窗里，颗粒基本完成了自由下落和初次碰撞，粘附与否已经能看出来。稳态覆盖率这类指标，可以用多个短时间窗的结果做统计外推，或者用“一次性释放N颗”的方式代替“持续释放”。

第三招是接触检测优化。EDEM的邻居搜索靠网格分箱，网格尺寸设得合适能省大量时间。我一般让网格边长略大于最大颗粒直径，太小会导致一个颗粒跨多个网格反复检索，太大会让每格里的候选颗粒过多。另外把“只在需要时才开启的接触模型”关掉，比如入门阶段不需要传热、不需要电荷、不需要破碎，这些模块挂着不用也吃内存。

三招叠完之后，一台普通多核工作站跑几千颗颗粒、一秒物理时间，大概几个小时能出结果，属于可接受范围。如果你连这个都嫌慢，那就把颗粒数再砍一半，先看趋势，别一上来就追求统计意义上的漂亮曲线。

3. 几何建模与颗粒工厂实操

几何和颗粒工厂是这个案例里最“手工”的部分，也是最容易出问题的地方。软件里点几下就能建个圆柱，但要让颗粒乖乖从工厂里出来、按预期落下去、落在晶圆上不弹飞，细节全在这儿。

3.1 晶圆与容器的几何处理

晶圆我用软件自带的圆柱体基元建，直径100 mm、高度0.3 mm，位置放在原点，上表面z=0。这里有个细节：晶圆的侧面和上下表面要分别命名。EDEM里几何体的每个面（section）可以单独指定接触参数和运动状态，如果整个圆柱只当成一个面，后面想单独统计“落在上表面的颗粒”就很麻烦。我在导入或创建几何之后，会在几何编辑器里把上表面单独拎出来命名成“Wafer_Top”，下表面和侧壁另起名字。这个习惯在后期做接触统计时省事很多。

包围盒我用长方体基元拼四面墙加一个底板，底板放在晶圆下方2 mm处，避免和晶圆共面导致接触判定混乱。四面墙的高度大概150 mm，顶部敞开。墙体的接触参数可以设成“颗粒一碰就消失”（用颗粒删除或者让颗粒飞出边界后自动移除），这样颗粒不会在角落里反复弹跳浪费计算资源。

注意：几何体之间千万不要有重叠或者共面。共面会导致接触法向计算出现奇异，表现为颗粒在某个平面上疯狂抖动或者直接穿过去。我最早搭这个算例的时候，底板和晶圆下表面放在同一个z坐标，结果颗粒穿模穿得一塌糊涂，排查了大半天才发现是共面问题。错开0.5 mm以上的间隙就能解决。

几何单位也要确认。软件默认可能是米，也可能是毫米，导入的CAD模型和手动建的基元如果单位不一致，会出现“晶圆大得像足球场”或者“颗粒小得看不见”的情况。我一般统一用米制SI，所有尺寸在建模时就换算好，避免中途改单位导致参数混乱。

3.2 颗粒工厂的生成逻辑

颗粒工厂是这个案例的名字来源——“颗粒物生成”。EDEM里的工厂分动态和静态两种，静态工厂在仿真开始瞬间一次性把颗粒全部生成完，动态工厂则在指定时间内按速率持续生成。研究沉降落点分布，我更倾向用动态工厂，因为一次性生成几千颗颗粒会在初始时刻造成严重的重叠和爆炸式排斥，颗粒初速度乱飞，落点完全失真。

工厂的核心设置有三项：生成体积、生成速率、粒径分布。生成体积我做成一个位于晶圆正上方80 mm处、尺寸120 mm×120 mm×2 mm的薄盒。薄一点的好处是颗粒初始位置集中在一个近似平面内，落点分布更接近“面源沉降”而不是“体源扩散”。生成速率设成每秒钟500颗，总生成2000颗，物理时间上大概4秒生成完——但因为我只算1秒，实际是分多次短算例跑，每次生成几百颗，累加统计。

粒径分布用正态分布，均值35 μm，标准差8 μm，并给上下截断（20~50 μm）。截断很重要，正态分布理论上会出现负半径，软件虽然会拦，但生成失败会在日志里刷一堆警告，看着心烦还影响排查真正的错误。初速度我给了向下的0.1 m/s，纯粹是为了让颗粒快点离开工厂区域，避免在工厂体积内互相碰撞干扰。这个速度远小于重力在1秒内产生的速度增量（约9.8 m/s），所以对最终落点影响可以忽略，属于“给个初速度省时间”的工程习惯。

提示：工厂生成颗粒时如果位置和已存在的颗粒重叠，EDEM会延迟生成或者直接报错。薄盒工厂在小粒径、高生成速率下特别容易出现这个问题。解决办法要么把生成速率降下来，要么把工厂体积在厚度方向加大一点，给颗粒留出互不干涉的空间。

3.3 接触模型与JKR参数设定

接触模型选错，整个算例就白搭。颗粒与晶圆、颗粒与颗粒之间，我用的是Hertz-Mindlin（无滑移）基础模型叠加JKR粘附模型。基础模型负责法向弹性力、切向摩擦力和滚动阻力，JKR负责在颗粒与表面分离时施加一个额外的拉脱力，模拟范德华力主导的粘附。

JKR模型的核心输入就是表面能。软件里填的这个值，物理含义是单位面积上的粘附能，实际计算拉脱力时按“拉脱力 = 1.5π·γ·R*”估算，其中R*是等效半径。我可以给你算一笔账，让你直观感受这个力有多大：取颗粒半径20 μm、表面能0.05 J/m²，拉脱力大约是1.5×3.14×0.05×2×10⁻⁵ ≈ 4.7×10⁻⁶ N。而同样一颗半径20 μm、密度2650 kg/m³的颗粒，重力大约只有 4.4×10⁻¹⁰ N。也就是说，粘附力比重力大了四个数量级。这就是为什么颗粒一旦贴上晶圆，靠重力、靠轻轻吹气根本弄不下来，也是为什么后道清洗工序要用超声、兆声、刷洗这些强制手段。

这个数量级关系直接决定了参数设定的思路：只要JKR表面能给到“粘得住”的量级，颗粒落地后基本就是粘住，不会弹跳。所以恢复系数、摩擦系数在这个算例里的敏感度其实不高，真正要调的是表面能。我在初次调试时把表面能设得偏小（0.005 J/m²），结果颗粒落地后弹几下就滚走了，表面一颗都留不住，看着像“仿真失败”，其实是粘附力不够，属于参数问题不是模型问题。

颗粒与颗粒之间我也开了JKR，因为落下来的颗粒如果聚成团，团块的沉降速度和粘附行为跟单颗粒不一样。但如果你的算例颗粒数很多，颗粒间粘附会让它们迅速抱团成一个大球，反而不利于观察单颗粒沉积。这种情况可以把颗粒间表面能设小一点，或者干脆关掉颗粒间粘附，只保留颗粒-晶圆的粘附。这是典型的“按研究目的取舍”，没有标准答案。

4. 求解器设置、运行与监控

参数填完，点开始之前，求解器这几项必须过一遍。很多新手跑出来的结果诡异，问题就出在这几个默认值上。

4.1 时间步长与计算域网格

时间步长我不用软件自动推荐的默认值，而是自己算一遍。Rayleigh时间步的估算公式是 t_R = πR·√(ρ/G) / (0.1631ν + 0.8766)，取它的20%作为实际时间步。以半径20 μm、密度2650 kg/m³、剪切模量3×10¹⁰ Pa、泊松比0.22代入，Rayleigh步长约在10⁻⁷ s量级，20%就是2×10⁻⁸ s。

这个数意味着什么？算1秒物理时间需要五千万步。即便颗粒数只有两千，也是极其漫长的。所以我前面说的“时间窗截断”不是偷懒，是必然选择。我的实际做法是把总时长设成0.2秒，几百万步，跑一晚上能出结果。想看更长时间的沉积，就分段跑，把上一段结束时的颗粒状态作为下一段的初始状态（EDEM支持从上次结果继续），或者干脆用不同的随机种子多跑几个短算例做统计平均。

计算域网格（邻居搜索网格）的边长，我设成最大颗粒直径的1.2倍左右，也就是60 μm。设太小，一个颗粒会横跨多个格子，接触检索时反复遍历；设太大，每个格子里塞的候选颗粒太多，检索量爆炸。这个值在颗粒粒径分布很宽的时候要特别小心，如果最大颗粒是最小颗粒的十倍，网格尺寸就很难兼顾，通常按最大粒径设，接受一点效率损失。

保存间隔也要设合理。设得太密，结果文件几个G，后处理打开都卡；设得太疏，颗粒落地瞬间的动态过程就丢了。我一般按总时长的百分之一设保存间隔，0.2秒的算例每2毫秒存一帧，一百帧，足够看清沉降过程。

4.2 计算过程中的观察点

点了开始不等于就可以去干别的。前几分钟必须盯着看，因为绝大多数建模错误会在最初几帧就暴露出来。

我通常关注三件事。第一，颗粒从工厂出来的时候是不是“炸开”的。如果一生成就四散飞射，说明生成位置重叠严重，颗粒间产生了巨大的排斥力，得回去调工厂体积或生成速率。第二，颗粒下落过程中有没有穿模。正常的轨迹是一条平滑的抛物线，如果某颗粒突然消失或者出现在晶圆下方，就是穿模，多半是时间步太大或者几何法向有问题。第三，颗粒落到晶圆表面后的行为。正常的现象是接触瞬间有轻微反弹，然后在一两个时间步内被粘附力“拽住”，之后基本静止或者缓慢滚动。如果颗粒一直弹、一直滚，说明粘附不够；如果颗粒接触瞬间就“定住”且位置诡异，可能是接触刚度太大导致数值刚性问题。

软件的实时监控窗口可以看颗粒总数、接触数、平均速度这些曲线。我一般会盯“与晶圆接触的颗粒数”这条曲线，它应该是阶梯式上升的，每一级台阶对应一批颗粒落地。如果曲线上下剧烈跳动，说明有颗粒在反复脱离和重新接触，粘附参数或者时间步需要复查。

5. 后处理，把“落在晶圆上”变成可量化数据

仿真跑完，画面好看没用，得拿出能写进报告的数和图。这一步我一般分三层：直接读数、空间统计、力学估算。

5.1 接触统计与区域划分

软件自带的分析模块里，几何体的接触数据可以直接导出。我关心的第一个指标是与晶圆上表面接触的颗粒数，以及这个数随时间的变化。这条曲线能告诉你沉积速率和饱和趋势。

第二个指标是空间分布。我用网格分箱工具把晶圆上表面划分成若干个小格，统计每个格子里接触颗粒的数量，换算成面密度（颗/mm²）。网格划分要跟晶圆尺寸匹配，100 mm直径的晶圆，划10 mm见方的格子，大概七八十个有效格，太少看不出分布规律，太多每格样本太少噪声大。

第三个指标是粒径分布。把落在晶圆上的颗粒按粒径分档统计，你会发现一个很有意思的现象：大颗粒落得快、动能大、容易反弹或滚走；小颗粒落得慢、粘附力相对重量更占优势、更容易留下来。最终表面上的粒径分布和释放时的粒径分布并不一致，这个“选择性沉积”恰恰是仿真的价值所在，靠实验中数颗粒很难数清楚。

5.2 用脚本算面密度和径向分布

软件导出的颗粒位置数据是CSV，字段包括ID、坐标、速度、半径。我用Python做二次统计，比在软件里点来点去灵活得多。

import pandas as pd import numpy as np # 读取导出的颗粒结果数据 df = pd.read_csv("wafer_particles.csv") r_wafer = 0.050 # 晶圆半径 50 mm，单位与导出数据保持一致 z_tol = 0.0005 # 表面判定带，0.5 mm z_top = 0.0 # 晶圆上表面高度 # 第一步：筛出高度上贴近上表面、且水平投影在圆内的颗粒 on_top = df[(df["Z"] >= z_top) & (df["Z"] <= z_top + z_tol)].copy() on_top["r"] = np.hypot(on_top["X"], on_top["Y"]) on_top = on_top[on_top["r"] <= r_wafer] # 面密度 area = np.pi * r_wafer ** 2 density = len(on_top) / area print(f"表面颗粒数: {len(on_top)}") print(f"面密度: {density:.6f} 颗/mm^2") # 第二步：按径向分环统计，看是不是边缘多、中心少 bins = np.linspace(0, r_wafer, 11) on_top["ring"] = pd.cut(on_top["r"], bins, labels=False) ring_stat = on_top.groupby("ring").agg( count=("ID", "count"), mean_d=("Radius", lambda x: 2 * x.mean()) ) print(ring_stat)

这段脚本里，z_tol的选取是个经验活。取得太小，颗粒因为接触变形微微陷入晶圆，坐标落在表面下方，就被漏掉了；取得太大，把正在下落的颗粒也算进去，沉积数虚高。我的经验是取最大颗粒半径的十分之一左右，配合接触数据交叉验证一下数量对不对得上。

径向分布那张表往往能看出问题。如果边缘环的颗粒数明显高于中心，可能是工厂的初始分布不均，也可能是包围盒侧壁反弹把颗粒“推”向了边缘。这种系统性偏差不是物理现象，是建模瑕疵，得回去改。

5.3 粘附力与去除难度的快速估算

后处理做到这里，其实还可以再往前推一步：算一下每颗颗粒受到的粘附力与自身重力的比值。这个比值直接对应“这颗颗粒有多难去掉”。比值越大，清除越困难。

做法是把拉脱力量级估算和颗粒重量估算对每一档粒径都算一遍，做成表。你会发现，粒径越小，这个比值越夸张；粒径增大到几十微米以上，重力才开始和粘附力同量级，颗粒才可能出现“抖一抖就掉”的情况。这个结论对工艺的意义很直接：小颗粒一旦落上，几乎不可能靠物理振动去除，必须靠化学或者流体力学手段削弱它的附着力。这就是为什么清洗工艺里表面活性剂、兆声波这些手段重要——它们不是“吹掉”颗粒，是降低界面能、削弱粘附力。

6. 常见问题与排查备忘

这套算例我从零搭到能出图，中间踩过的坑不少。下面按症状归类，给一份可以直接对照的排查清单。

6.1 颗粒穿模、飞散、莫名消失

穿模是所有DEM算例的头号问题，症状是颗粒穿过晶圆或者墙壁出现在计算域外。原因通常有三个：时间步太大、接触刚度太高、几何法向错误。排查顺序我建议从时间步开始，先把步长砍到原来的四分之一，如果穿模消失，就是步长问题；如果还在穿，检查几何体之间有没有重叠或共面；都没问题，再去看接触刚度的设置是不是超出了材料合理范围。

颗粒生成时“炸开”是另一类常见现象，本质是初始位置重叠。工厂体积相对于颗粒粒径太小、生成速率太快，都会导致新颗粒生在老颗粒体内，软件强行把它们推开，产生巨大的初始速度。把生成速率砍半，或者把工厂在生成方向上的厚度加大，通常能解决。

6.2 颗粒一颗都不粘，或者全粘成一坨

这两个极端都是粘附参数的问题，但方向相反。一颗不粘，要么是JKR表面能设小了，要么是恢复系数太大导致颗粒一碰就弹飞。可以先试试把表面能提高一个数量级看趋势，如果立刻粘住了，就说明量级不对，再回来细调。我最初就吃过这个亏：按直觉设了个“看起来很小”的表面能，结果什么也粘不住，后来一算拉脱力和重力的比值，发现小了两个数量级。

全粘成一坨，是颗粒间粘附太强，颗粒释放后迅速抱团，落下来是一个团块而不是分散颗粒。解决办法是把颗粒间的表面能单独调小甚至置零，只保留颗粒-晶圆的粘附。很多新手以为接触模型参数必须两两一致，其实完全可以根据研究目的区别对待。

6.3 算到一半卡死或者慢得离谱

先看是不是时间步突然被拉小了。有些情况软件会根据最大速度自动减小步长，比如有颗粒获得了异常大的速度（穿模前的征兆），步长就会被压得极小，计算直接卡住。遇到这种情况，去看日志里的最大速度，找出异常颗粒的位置，基本能定位到几何问题。

再看颗粒数是不是超了预期。工厂按速率持续生成，如果没设总量上限，颗粒会一直累积，跑着跑着就上万颗，速度自然崩。我习惯在工厂里明确设置总生成量，并在求解器里开启“超出计算域自动移除”，把飞出去的颗粒清掉。

6.4 常见问题速查表

提示：表格里每一条都是我在实际调试中真真切切遇到的，不是抄来的清单。真跑算例的时候，把这几个现象和参数对应的关系记熟，能省下大量瞎试的时间。

最后说点个人的使用体会。这个算例最大的价值不在结果本身——缩比模型算出来的面密度数字，没人会拿去做工艺判定。它的价值在于把“颗粒怎么来、怎么落、怎么粘”这条因果链条在计算机里跑通一遍，让你对每个环节的关键参数建立起手感。等以后遇到真的问题，比如某台设备搬片时扬尘、某个腔体里颗粒回吸到晶圆上，你能迅速判断该建什么模型、该重点调哪个参数、结论的可信边界在哪。入门案例就该这么用：不求准，求通。