统 仿 真 学 报? Vol. 18 No. 10
2006年10月 Journal of System Simulation Oct., 2006 第18卷第10期 系
多尺度仿真方法研究综述
孙西芝，陈时锦，程 凯，初文江
(哈尔滨工业大学精密工程研究所, 哈尔滨 150001)
摘 要：多尺度仿真方法结合了分子动力学方法和有限元原理，可以跨越从宏观到微观的多个尺度，具有分子动力学的精确性，而且计算量相对较小，扩展了分子动力学的应用范围。介绍了国外的多尺度仿真方法的发展及应用概况，阐述了几种较成熟的多尺度仿真方法的基本原理，指出了各自特点及不足，最后对其发展方向作了展望。
关键词：多尺度；分子动力学；有限元；计算机仿真
中图分类号：TP391.9 文献标识码：A 文章编号：1004-731X (2006) 10-2699-04 Overview of Multiscale Simulation Method Research
SUN Xi-zhi, CHEN Shi-jin, CHENG Kai, CHU Wen-jiang
(Precision Engineering Research Institute, Harbin Institute of Technology, Harbin 150001, China)
Abstract: Multiscale simulation method, combing molecular dynamics and finite element method, not only has the MD’s
precision, but also exceeds its range of application from microscale to marcroscale and owns small calculating amount relatively, expanding the range of application of
molecular dynamics. The development and application of several multiscale simulation methods were introduced, focusing on the basic principles of some comparatively advanced ones. Then their advantages and disadvantages were discussed. Finally, the future development of this method was predicted. Key words: multiscale; molecular dynamics; finite element; computer simulation
纹扩展进行了仿真，得到了较为理想的结果。随后Satashi Izumi[5]等人根据Si晶体的复杂性，对Kohlhoff的方法又作了进一步改进。并且，Tabar[6]等人在研究脆性材料的断裂问题上也采用了多尺度方法，并成功地再现了裂纹速度的振荡。与此同时，Robert[7]等人对硅和石英材料的微谐振器的振动特性进行了计算。Smirnova[8]等人应用类似方法对应力波的传播进行了研究，并与传统的分子动力学方法进行了比较，得出了误差很小的结论(5%)。还有Pillai[9]等人对Bi晶体的晶面缺陷的仿真及Tadmor
[10-11]等人对PbTiO3 晶体的磁滞现象的研究等等。
引 言近年来，研究材料的微观力学特性的分子动力学方法(MD)发展迅速，它建模简单，程序短小，可计算的原子体系大大超过第一原理等方法，在解释一些用理论分析和实验观测等方法都难以了解的微观现象上起到了不小的作用。
随着计算机计算能力的不断提高以及算法的改进，分子动力学方法可处理的原子已经数以亿计[1],但仍达不到仿真实际系统的要求，在时间和空间尺度上受到极大的限制。为解决这一难题，有人提出了多尺度仿真(Multiscale Modeling)，即把原子模型嵌入到连续介质模型中，采用分子动力学方法计算我们感兴趣的微小区域，而其他区域采用连续介质力学方法(多为FEM)计算，不仅减小了计算量，而且使计算尺度得到了极大的扩展。 2 多尺度仿真基本原理
原子区及连续介质区的计算方法都已相对成熟，主要难题集中在如何处理原子模型和连续介质的结合区上，对此许多学者各自采用了不同的方法。 1 多尺度仿真方法的发展及应用
早在上个世纪80年代初，Baskes和Mullins等人就已经运用类似的方法对晶格的裂纹扩展问题进行了仿真，但受条件制约，他们所建立的模型与真实物理条件有很大差异。后来，Kohlhoff[4]等人改进了他们的方法，对b.b.c晶体的裂 [2] [3]
2.1 FEAt方法
Mullins[3]提出，在结合区域，原子与节点一一对应，直接应用原子作用力转化成集中载荷作用于有限元节点上。这种想法简单明了，但存在许多问题，如由于原子间作用力的长程及非局部特性，很难解决区域边界处的受力平衡。而如果用节点(原子)位移来代替力的直接作用，则可避免这一问题。由Kohlhoff等人提出的FEAt方法是这一应用的典范[4]。
在他们提出的方法中，整个模型由四个部分组成(Ⅰ,Ⅱ,Ⅲ,Ⅳ)，如图1所示。 收稿日期：2005-07-25 修回日期：2005-11-28
作者简介：孙西芝(1978-), 女, 山东金乡人, 博士生, 研究方向为纳微米切削加工过程机理的建立和仿真；陈时锦(1964-), 男, 安徽肥东人, 副教授, 研究方向为超精密加工装备的设计; 程凯(1961-), 男, 黑龙江哈尔滨人, 教授, 博导, 研究方向为3D精密与超精密表面质量和功能性的控制与加工技术, e-制造技术: 模式、定量分析方法及系统实现。
·2699·
第18卷第10期 Vol. 18 No. 10 2006年10月 系 统 仿 真 学 报 Oct., 2006
过渡区，完整的硅晶格与一个有限单元相对应。为了计算晶
MD
格内原子的运动对单元产生的作用，引入了内位移的概念，即原子的位移是晶格位移和内位移的和，并引用Martin公式计算内位移的值。在确定了内位移的计算方法后，便得出了
FEM
图1 MD区与FEM区的边界关系 仿真的具体计算步骤如下：
(1) 初始化MD区的原子位置及FEM区的边界条件。 (2) 计算MD区的原子受力，从而得到区域Ⅱ内的原子位移，接着计算内位移。
(3) 根据MD区的计算结果固定区域Ⅱ的节点，计算FEM区的反作用力。被固定节点位移由对应原子位移与内位移的差确定。
(4) FEM区域计算，得到所有节点的位移。 (5) 根据区域Ⅲ的节点的位移计算内位移。
(6) 根据FEM区计算结果固定区域Ⅲ的原子，计算MD区全部原子的位移。被固定原子的位移由对应的节点位移与内位移的和确定。 (7) 检验收敛性，如不收敛则返回(2)进行计算。
MD区域从Ⅰ到Ⅲ，FEM区域从Ⅱ到Ⅳ。Ⅱ,Ⅲ为过渡区，其中的原子与节点一一对应，Ⅱ区原子为FEM提供边界条件，节点随原子运动；Ⅲ区的节点为MD区
域提供边界约束，原子随节点运动。模型采用了非局部弹性理论来描述过渡区。在线性变化范围内，晶格的弹性变形能可表示为：
E= 1 nm
∑∑φik(xn?xm)(uin?uim)(uk) (1) ?uk2nm 其中xn表示原子n 的位置，uin
表示原子n离开平衡位置的
位移，φik为能量算子。其相对应的均匀各向同性连续介质的变形能可表示为： 1E=∫∫φik(x?x')[ui(x)?ui(x')][uk(x)?uk(x')]dV'dV[0,1] (2) 2
vv2V0'可见，以上两式将离散与连续两种能量的计算联系了起来。(2)式可变换为应变形式：
E= 1
∫∫cijkl(x?x')εij(x)εkl(x')dV'dV (3) 2vv' MD FEM
式中：cijkl为连续介质的二阶弹性常量，εij(x)为位置x处的应变。将弹性能的应变形式泰勒展开：
E(ε)=E(0)+(?E/?εij)0εij+ 1216
[?2E/(?εij?εkl)]|0εijεkl+
[?E/(?εij?εkl?εmn)]|0?εij?εkl?εmn+.... (4) 3
图2 Si晶格与单元的对应关系 2.2 QC方法
QC(Quasicontunuum)方法由Tadmor和Ortiz[12]等人于1996年提出，其目的是构建一个完整的原子系统而不必显式地对每一个原子进行计算。它的原理与FEAt方法有很大的不同，体现在不再是由一特殊的过渡区来连接原子区域和连续区域，而是以原子计算为基础，通过引入有限单元概念来减小原子的计算量。 (5)
其中位移ui在仅原子的新位置xi可由xi=Xi+ui表示，
考虑其物理意义时可被看作是连续区域u(X)在Xi的值即
为了保证晶格与对应单元的应力平衡，应使式中的各系数在两种介质的计算中分别对应相等，因此做如下定义：
cij=(?E/?εij)|0cijkl=[?2E/(?εij?εkl)]|0cijklmn=[?E/(?εij?εkl?εmn)]|0 3
u(Xi)=ui，今后在公式中出现的u表示这一区域内所有原子位移的集合{u1,u2,...uN}。
为了获得系统的总能量，我们在原子仿真过程中通常总是要对每一个原子的能量进行计算，即
Etot=∑Ei(u) (6) i=1N
这样使连续介质的弹性常量由原子间的势能函数定义，达到了应力平衡的目的。 FEAt方法是开发较早的多尺度方法之一，它直接把离
散量和连续量结合了起来，体现了多尺度仿真的早期思想，并且它被应用在了一些结构简单晶体(如b.b.c, f.c.c)的仿真上，取得了不错的效果。但是由于方法本身的粗糙性，其对于复杂结构的晶体(如硅晶体)不再适用，因为原子与节点的一一对应关系很难保证，晶格内原子的位移给弹性模量的确定带来困难。因此，后来的很多学者在FEAt方法的基础上进行了新的改良，如Satashi Izumi等人的研究就是一例，他们把这种方法应用在了硅晶体的仿真上[5]，如图2所示。在
其中Ei表示原子i在位移u作用下的能量，它的准确性取决于势函数的选择。在变形相对平滑的区域对每个原子的位移都进行计算是不必要的。基于这种考虑，QC采用了近似方法，即选择一些有代表性的原子repatoms(representive atoms)计算其ui，而其他原子的位移则通过线性插补近似得到，这里恰好利用了有限元方法的插补原理，如图3所示。 ·2700·
第18卷第10期 Vol. 18 No. 10 2006年10月 孙西芝, 等：多尺度仿真方法研究综述 Oct., 2006
精确又简便。一种以力学计算为基础的QC方法也在研究之中，详细内容见参考文献[14]。
QC方法已被广泛应用在了结构破坏、微缺陷、裂纹扩展等一系列问题的处理上，是目前较完善的多尺度仿真方法，它避免了使用不够准确的理论方法将分子动力学和有限元的直接耦合，而是将有限元概念融合到了分子计算中，使
多尺度仿真方法研究综述解读
