1、基于元胞自动机算法的三维玻璃微流控芯片湿法刻蚀模拟研究摘要：以二维动态元胞自动机理论为基础，提出一种针对玻璃各向同性刻蚀特性的三维动态元胞自动机（Cellular automata）算法。利用创建刻蚀链表的方式代替传统的元胞搜索，提高了算法效率，并通过改进元胞信息的存储方式，使得三维元胞阵列得以存储于二维矩阵中。最后利用OpenGL 技术对刻蚀结果进行三维成像。本文在个人计算机上实现了较高分辨率（40004000 像素）的三维刻蚀模拟与显示。通过对比玻璃刻蚀的实验结果和前人的研究评估模型，模拟结果完全能够实现对刻蚀结果的预测。本文的研究表明：建立的三维元胞自动机算法可以有效地应用于玻璃刻蚀过程
2、的模拟。关键字：玻璃芯片，湿法刻蚀，元胞自动机，OpenGL 三维成像，计算机模拟Simulation on 3-D Glass Isotropic Chemical Etching Process for microfludic chip Based onCellular Automata ModelAbstract(：For Isotropic chemical etching properties of glass, an algorithm based on improved continuous cellular automata (CA) model is produced, wh
3、ich develops the traditional CA model. It can be used for wet etching of glass in three-dimensional (3-D) space. To improve the efficiency of algorithm, structure with linked list data is used to store etched cell. With the improvement of the storage mode, a 3-D cellular matrix can be stored in 2-D
4、Array. The result of simulation is displayed with OpenGl 3-D technology. The glass etching simulation program based on common personal computer platforms has been realized with high resolution（40004000 Pixels）and visualizing results in 3-D space. The results of etching experiments are used as compar
5、ison to assess the implementations, which show that the model is able to forecast the result of etching glass. This study shows that 3-D CA model is an effective solution to the glass etching process task. Key Words: glass chip; wet etching; cellular automata; OpenGL; computer simulation1 引言1994 年 O
6、.Than 等人成功地将 CA 算法应用于硅晶体的各向异性刻蚀仿真 1。经Zhenjun Zhu 等人的改进 2，该模拟算法已较为成熟。随着微细加工技术的发展，90 年代初，玻璃芯片的出现弥补了硅芯片光电特性不足，且由于其廉价的成本，十多年来玻璃微流控芯片的研究得到重视，已经被广泛应用于电泳和色谱分离 3。由于玻璃为非晶体，相比于硅晶体各项异性的刻蚀过程，玻璃的刻蚀速率各向相等，且 O.Than 等人提出的刻蚀模型仅适用于分子结构排列整齐的晶体。1995 年Karafyllidis I 等人建立了基于元胞自动机算法的光刻胶刻蚀模型 4，2005 周再发等人对其进行改进，建立了二维动态元胞自动机
7、刻蚀模型 5，该模型不同于晶体硅的刻蚀模型，并没有将原子作为元胞单元，而是将被刻蚀工件细分为元胞，因此该模型便可以适用于非晶体材料的刻蚀模拟。但由于模型的限制，算法仅实现了对二维刻蚀过程的模拟。针对玻璃各向同性的刻蚀过程，本文在较成熟的二维动态元胞自动机算法的基础上，建立一套完整的三维动态元胞自动机各向同性刻蚀模型。由于采用了动态分配内存等优化技术，节约了可观的内存空间；结合链表存储刻蚀表面等优化方法，大大地提高了模拟算法的计算效率。这些优化处理使得原本占用较大内存空间的玻璃三维刻蚀模拟在个人计算机上的计算成为可能。2 三维刻蚀模型2.1 元胞及其状态将被刻蚀的玻璃按边长为 a 的正方体细分成
8、三维阵列，每个正方体称为 CA 的一个元胞。它是仿真运算中的基本单位。玻璃的刻蚀是自顶而下的，且各向同性。因此在本文算法中，我们将每一纵列的元胞组合成高为 h 的长方体（如图 1 所示） ，而刻蚀过程是对每一个长方体至顶而下的刻蚀。这样就将三维元胞阵列存储于二维矩阵中，用高度 h 表示该纵列剩余元胞的个数及状态。在刻蚀模拟过程中，刻蚀总是从二维矩阵单元顶部元胞开始的。该元胞的刻蚀受到位于其顶部及四周的刻蚀液影响（图 1） 。在某一时刻 t，元胞的状态 Ci,j(t)定义为此时元胞剩下体积 Vl 与一个完整元胞体积 Vc 的比值：(1),lijctt这样在刻蚀过程中某一个元胞的状态会在“1”（未
9、被刻蚀）和“0” （完全刻蚀）之间分布。因此每个长方体（阵列单元）均可用一个浮点数表示其高度 h，它是所有该纵列元胞状态值之和。即： (2),()()ijijtCt某一时刻 t，刻蚀掉的体积 Ve(t)为位于其顶部刻蚀液刻蚀掉的体积 VTop及四周刻蚀液刻蚀掉的体积 VRound之和：(3)eTopRundt(4)(1)llett图 1 纵列元胞所组成的阵列单元2.2 时间补偿值在刻蚀过程中，当元胞接触刻蚀液时便开始刻蚀模拟。这通常是由于其周围某个元胞被完全刻蚀，刻蚀液流入并与其接触造成的。但由于元胞并不总在时间步长 T 的整数倍的时间内被完全刻蚀。如果元胞在 t 时刻和 t+T 时刻之间被刻
10、蚀，就会对周围元胞的刻蚀起始时间判断造成影响，如果忽略这种影响就会导致运算结果的不准确。为解决这一问题我们引入一时间补偿值的概念，即将时间步长 T 完全刻蚀元胞后剩余的时间 Tc(t+T)补偿到下一个时间步长用于刻蚀其邻居元胞。如果在某一时刻 t 到 t+T 内长方体纵列的元胞数量没有减少，那么该元胞对应时间补偿Tc(t+T)为零。否则需要对它的时间补偿值进行计算。计算方法如下：已知时间步长 T 内刻蚀掉体积 Ve(t)，那么设完全刻蚀掉当前元胞需要时间 te，当前元胞剩余体积为 Vl(t)，那么很容易得到如下关系式：, (5)()elettT()letT因此时间补偿为：(6)()()1lce
11、eVtttT可统一表示为下式：(7)0,()()1(clehttTtVtTT2.3 刻蚀效果的计算计算来自元胞顶部刻蚀效果时，一个时间步长 内刻蚀掉的体积为刻蚀的线速率 R 与元T胞纵向刻蚀面积 S（S 的取值将在 2.5 中讨论）和时间步长 T 的乘积：(8)opVR计算来自元胞四周的刻蚀效果时，需要考虑在三维空间与其相邻的 24 个元胞对它的刻蚀作用，即水平摩尔邻居位（如图 2 元胞位（i,j）周围白色与灰色的元胞）上的纵向相邻的三个连续元胞。图 3 为元胞位（i,j）处某元胞受到水平邻位（i,j）的三个元胞的刻蚀作用。图 2 元胞（i,j）与其摩尔邻居图 3 来自水平邻位（i,j）的刻蚀
12、在摩尔邻居中，来自四角邻居(图 2 周围灰色元胞) 的刻蚀效果与冯诺依曼邻居（图 2周围白色元胞）不同，所以计算时将两部分的影响分开求解。设某一时刻 t 至 t+T 时段内四角的刻蚀体积为 Ve(t+T)diagonal，冯诺依曼邻居的刻蚀体积为 Ve(t+T)adjacent。那么总的刻蚀体积 Ve(t+T)round 为：(9)()()roundediagonleadjcenttTtTtT为了计算刻蚀掉的体积，我们首先计算刻蚀掉区域的水平截面积 Se(t+T)round（假设截面处处相等） 。(10)()()()eroundediagonleadjcentStTtTtT2.4 冯诺依曼邻居
13、和四角邻居的刻蚀效果冯诺依曼邻居和四角相邻元胞的刻蚀效果需要分别计算，并且由于有可能出现刻蚀区域的交叠现象（图 4 黑色区域） 。为了提高计算的精度，通常采用几何方法解决该问题 5。图 4 被刻蚀的元胞水平截面产生的刻蚀交叠效果通过几何方法得到完整的冯诺依曼邻居的刻蚀效果见公式（11） 。其中刻蚀线速度 Ri,j 取决于水平邻位（i,j）元胞高度情况，可由公式（13）求得。21,(1,),1(,1),()max(0(eadjcent ijijijijStTRTcttc(11)对角相邻元胞（i+1,j-1）的刻蚀效果见公式（12） ，其中参数 D 用于描述间接流入与直接流入之间的差异。同样为了提
14、高运算精度也要避免某些区域被刻蚀多次。处理方法为：如果与某对角邻居元胞相邻的两个冯诺依曼邻居（如图 3 右中元胞（i,j-1）和（ i+1,j）相对于（i+1,j-1 ） ）同时存在刻蚀液，那么参数 Di+1,j-1取 0；如果只有一个存在刻蚀液，参数 Di+1,j-1取 D/2；如果都没有刻蚀液，参数 Di+1,j-1 取D。参数 D 取值为 0.184,6。221,1,1,()()ediagonlijijijStTRTct(12)下面讨论刻蚀线速率 Ri,j的选取，在纵向界面内同样采用类似摩尔邻居的讨论，刻蚀线速率 Ri,j 可表示为公式（13） 。即当水平邻居元胞（i,j）和（ i,j）
15、正上方元胞位都未被刻蚀，该方位无刻蚀液，刻蚀速率为 0；当只有（i,j ）上方元胞位为空时，刻蚀速率为；当水平邻居元胞（i,j）和（i,j）正/2D上方元胞位都为空时，刻蚀速率为 ；当水平R邻居元胞（i,j ） ， （i,j ）正上方元胞， （i,j ）正下方元胞位都为空时，刻蚀速率为。(1/2)R, ,021(1)ijijijijij ijijijijijhDRhh (13)2.5 元胞纵向与横向刻蚀交叠的处理考虑单一元胞，其水平方向的刻蚀效果和纵向方向的刻蚀效果之间同样存在交叠，如图5（a）所示，深色区域为纵向和侧向的交叠部分。而实际上纵向的刻蚀并非仅仅体现在当前被刻蚀的元胞上，考虑图 5
16、 中（b）的情况，当纵向和横向的刻蚀同时进行时。纵向的刻蚀面积应为当前刻蚀元胞的横截面 S1 与正下方元胞裸露的表面积 S2 的和。即纵向的刻蚀面积 S始终为初始元胞的底面积 a2。图 5 (a)单一元胞纵向和横向的刻蚀交叠效果。(b)为纵列元胞刻蚀的纵向刻蚀面积 S=S1+S2 =a2 在刻蚀区域的水平截面积 Se(t+T)round 转换成刻蚀区域体积 Ve(t+T)round 的过程中，由于刻蚀后的元胞可以认为是等水平截面的立柱，因此只需要将刻蚀区域的水平截面积 Se(t+T)round乘以其刻蚀高度 he 即可。在刻蚀过程中，来自元胞顶端的刻蚀效果使元胞高度始终发生变化。为提高算法精度
17、可将 he 取高度变化过程的平均值，即 。2)()Tte2.6 元胞边长的选取在玻璃刻蚀的元胞自动机模拟算法中，元胞边长 a 的选取会影响模型计算的精度。由于在模型中采用一个浮点数表示长方体的高度值h，所以边长 a 的选取不会影响模拟计算垂直方向的精度，但却会影响水平方向的精度。边长 a 的选取同样会影响到模拟计算的最大水平范围，而水平精度和最大水平范围两个因素相互矛盾，即如果降低边长 a 的尺度，就会提高水平精度，但能够模拟的水平范围就降低了；反之，则会牺牲水平精度来提高模拟的水平范围。如果边长 a 值选为 0.5m，因此其水平精度为 0.5m。当计算中矩阵规模为。因此其水平范围为 。201
18、m2.7 时间步长的选取由于模型中计算刻蚀过程是以时间步长 T为运算的时间单位。T 正比于刻蚀算法的速度，在满足算法精度的前提下 T 的选值应尽可能大。采用时间补偿的动态元胞自动机算法可以解决相邻元胞在时间步长 T 的非整数倍的时间内被刻蚀的问题，但在一个时间步长 T 内如果某一元胞上刻蚀的体积超过其初始体积 a3，就有可能导致刻蚀效果在一个时间步长内的传递超出相邻元胞，而目前算法并不能对超出相邻范围的元胞进行时间补偿。所以应尽可能使任意时刻每个元胞上产生的刻蚀效果 Ve 都不大于元胞初始体积 a3。考虑一种极端情况，当某元胞除了正下方元胞以外的 25 个临位元胞处均存在刻蚀液，这时 Ve 有
19、可能达到极大值，其表达式为： 22(1,)(1,)3(,)(,)max(1) 02eTopRundijijijijaTRaDcttca (14)其中 D=0.18，元胞的 Tc 最大取值不应超过 T，因此公式（14）可化简为：(15)2231()max4.36,0eVRTRa假设 ，其中 k 为待求系数，则公k式（15）可转化为关于 k 的不等式：(16)32234.36()(ax,0)a ak解得： 。但上述这种极端情况在6.19k实际运算中很少出现，为了均衡运算效率考虑，进行模拟实验讨论 K 的取值。在实验中将单个时间步长内，某个元胞的刻蚀体积超过 a3 的情况认为会对模拟精度造成影响，并
20、统计该情况出现的概率，如图 6 所示。-0.51.53.55.57.59.511.513.515.517.51.5 2.5 3.5 4.5 5.5 6.5K误差百分比（%）图 6 不同 K 值对运算精度的影响当 K 取值大于 4.0 后对精度影响已经小于 0.0001%，完全可以满足工程模拟需要，所以在实验中折中选取为：(17)4aTR3 算法实现及优化根据面向对象编程思想，将模拟算法的实现部分抽象出两个类:Lattice 类、Mass 类。Lattice 类用于记录和操作一纵列元胞组成的长方体，包括：记录着长方体高度 的私有成h员（private）float m_Depth，记录时间补偿值
21、Tc 的私有成员 float m_Tc，记录长方体在二维阵列中坐标的 short int i,j 等成员；刻蚀方法 void Lattice:erode(float v)完成在每个时间步长 T 结束后对刻蚀元胞的状态 C 和长方体高度 h 的计算与更新，并依据公式(7)计算出时间补偿值 Tc。除此之外，还有一些用于操作私有成员的方法。Mass 类表示所有 Lattice 组成的被腐蚀的玻璃器件。为了实现高分辨率的刻蚀仿真，节约刻蚀过程占用的内存空间是必不可少的。传统的方法是将全部元胞的信息保存于二维矩阵内，但通常情况下，被刻蚀的区域仅占整个工件的一小部分，如果对整个工件的元胞都分配内存空间显然
22、浪费了大量的内存空间。本文的算法是创建成员 Lattice* *lattice，用于记录二维平面上所有 Lattice 实例的地址,并且仅对被刻蚀的元胞分配 Lattice 类实例空间，其余指针均赋值为 NULL。最先被加入链表中的是掩模版镂空区域的元胞地址。由于指针类型占用的内存空间（4 Byte）远小于 Lattice 类实例占用的内存空间（16 Byte） ，因此节约了可观的内存空间。在试验测试中，对于大多数掩模版的刻蚀模拟，采用新的内存分配机制可以节约 70%以上的内存开销。运算效率是算法另一个重要的属性，传统的元胞自动机刻蚀算法采用标定每个元胞的方法确定需要被刻蚀的元胞。这样每个仿真
23、周期都需要对全部元胞阵列搜索一遍，对于仅有部分元胞被刻蚀的阵列来说，这样的搜索明显浪费了大量的时间。本文采用创建一个线性链表成员 std:vector vecSurface，用于记录所有正在被刻蚀的 Lattice 的地址，每次刻蚀过程只需将链表中记录的 Lattice 刻蚀，并将新增的刻蚀区域加入链表，提高了刻蚀过程的效率。类方法 void Mass:simulation()用于实现玻璃器件的刻蚀模拟流程。算法流程见图7。根据掩模版、刻蚀时间、线速度等信息初始化模拟参数将线性链表中记录着的全部 Lattice 进行刻蚀将新产生的被刻蚀Lattice 加入链表是否完成 否输出刻蚀模拟结果是图
24、7 刻蚀模拟算法流程4 模拟与讨论在试验中，分别对模型的水平方向性和横向纵向刻蚀速率比进行了测试。前者采用将圆形镂空的掩模版进行刻蚀模拟，并观察刻蚀结果的水平截面是否为圆形，刻蚀结果见图 8，可见刻蚀模型有着较好的水平方向性；后者采用对直线形镂空掩模版进行刻蚀模拟，并分别观察横向的刻蚀距离和纵向的刻蚀深度。经过测试得到横向与纵向刻蚀速率之比约为 1.35。在实际的实验中横向与纵向刻蚀速率之比会随着刻蚀液不同有所改变，其中当 HNO3 含量为10%，HF 的含量为 13%时纵向的刻蚀速率为1.22m/min，且横向与纵向刻蚀速率之比为1.377，这与模拟数据是非常接近的。图 8 模拟结果的三维视
25、图（a）玻璃刻蚀表面的完整视图， （b）被刻蚀圆形区域的放大图， （c）右上角四分之一圆局部放大图采用一个相对复杂的掩模版图片（图 9） ，图像的分辨率为 。并设置初始化20参数:刻蚀速率为 1.5m/min，刻蚀时间为1000 秒，并动态监视刻蚀过程。模拟过程最大占用 49.6 MB 空间，模拟耗时 64 秒（采用 P4 3.0 GHz 的 CPU） 。图 11 为模拟刻蚀过程的快照图。图 9 掩模版图形图 10 模拟过程的动态显示快照图，其中（a,b,c,d）分别对应的刻蚀时间为 100 秒，300 秒， 700 秒，900秒。图 11 为最终刻蚀结果的三维全景图和局部的特写图，用户可以采
26、用任意视角对刻蚀结果进行观测。刻蚀后的玻璃器件形貌清晰可见，一些细微特征也被很好地体现出来。图 11 模拟结果的三维视图：（ a）玻璃刻蚀表面完整视图， （b）局部区域俯视图， （c ）局部区域侧视图下面是一些玻璃刻蚀的实验照片(光学显微镜照片) 与采用本文模型进行的模拟计算得出的刻蚀三维视图的对比图（图 12）和成品玻璃微流控芯片与采用本文模型进行的辅助设计图（图 13） 。图 12 实验与模拟对比图(a,c)为玻璃刻蚀的实验照片 8，（b,d）为计算机模拟结果图 13 玻璃芯片的计算机辅助设计： (a)玻璃微流控芯片照片 8, （b,c,d,e）基于模拟算法的辅助设计图5 结论玻璃的各向同
27、性模拟是预测玻璃芯片三维结构的关键。随着玻璃微流控芯片的制作工艺不断向成熟、快捷、精准方向发展，其刻蚀模拟技术显得更加重要。作为第一个可以运行于个人计算机上的三维玻璃的各向同性刻蚀模拟程序，本文以成熟的二维连续 CA 模型及其动态方法为基础，建立了三维各向同性刻蚀的模型。通过对算法进行优化，使得动态 CA 模型的模拟速度得到了较大地提升，并使内存的占用量显著地下降。这也为该算法的广泛应用奠定了基础。运用已有的玻璃刻蚀的实验结果和前人对玻璃湿法刻蚀的研究结论验证了本文提出的三维动态 CA 算法的效果，该算法有稳定性好、运算速度快的优点。模拟结果与实验结果的一致性表明了三维动态 CA 算法可以有效
28、地用于玻璃各向同性刻蚀模拟。通过模型可以准确地描述不同工艺过程、工艺参数对最终刻蚀结果的影响。这也有助于进一步实现 MEMS 加工过程的模拟及计算机辅助设计。参考文献1 THAN O. ; BTTGENBACH S. Simulation of anisotropic chemical etching of crystalline silicon using a cellular automata model. Sensors and actuators: A, Physical, 1994, 45(1):85-892 Zhu Zhen-Jun, Liu Chang. Micromachini
29、ng process simulation using a continuous cellular automata method. Journal of Microelectromechanical Systems, 2000, 9(2):252-2613 方肇伦.微流控分析芯片的制作及应用.北京：化学工业出版社,2005 年4 Karafyllidis I., Thanailakis A. Simulation of two-dimensional photoresist etching process in integrated circuit fabrication using cel
30、lular automataModeling Simulation in Materials Science and Engineering, 1995, 3(5):629-6425 周再发,黄庆安,李伟华等.光刻胶刻蚀过程模拟的元胞自动机算法研究.中国集成电路.2005,74(7):68-736 Zhou Zai-Fa, Huang Qing-An, et al. A novel 2D dynamic cellular automata model for photoresist etching process simulation. Journal of Micromechanics an
31、d Micoengineering, 2005, 15(3):652-6627 周健,闫桂珍.Pyrex 玻璃的湿法刻蚀研究.微细加工技术,2004,4(12):41-448 黎永前.电动微流体运动特性实验研究.大连：大连理工大学，大连理工大学博士后出站报告,2005 年BackgroundThe work is supported by the National Natural Science Foundation of China under grant No. 50605006: “Fracture Mechanism of Metal Film Microelectrodes unde
32、r Spatial Stress and Solving Methods”. During the course of thermal bonding of polymer microfluidic chips with electrochemical detection, metal film microelectrodes on the chip are liable to fracture easily on the boundaries between the cover plate and the substrate plate. The fracture mechanism is
33、studied by the fractographic analysis method. In order to resolve the fracture problem, a novel thermal bonding method is first put forward and studied. To the best of our know, there are no similar reports. The research of this project will accelerate the mass production and practicality of polymer microfluidic chips with electrochemical detection, and hold great theoretical and practical values for the packaging of polymer micro-electromechanical systems.
