1、統 仿 真 學 報© Vol. 18 No. 102006年10月 Journal of System Simulation Oct., 2006第18卷第10期 系多尺度仿真方法研究綜述孫西芝，陳時錦，程 凱，初文江(哈爾濱工業大學精密工程研究所, 哈爾濱 150001)摘 要：多尺度仿真方法結合了分子動力學方法和有限元原理，可以跨越從宏觀到微觀的多個尺度，具有分子動力學的精確性，而且計算量相對較小，擴展了分子動力學的應用范圍。介紹了國外的多尺度仿真方法的發展及應用概況，闡述了幾種較成熟的多尺度仿真方法的基本原理，指出了各自特點及不足，最后對其發展方向作了展望。關鍵詞：多尺度；分子

2、動力學；有限元；計算機仿真中圖分類號：TP391.9 文獻標識碼：A 文章編號：1004-731X (2006) 10-2699-04Overview of Multiscale Simulation Method ResearchSUN 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

3、, combing molecular dynamics and finite element method, not only has the MDsprecision, 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 s

4、everal 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;

5、 computer simulation紋擴展進行了仿真，得到了較為理想的結果。隨后Satashi Izumi5等人根據Si晶體的復雜性，對Kohlhoff的方法又作了進一步改進。并且，Tabar6等人在研究脆性材料的斷裂問題上也采用了多尺度方法，并成功地再現了裂紋速度的振蕩。與此同時，Robert7等人對硅和石英材料的微諧振器的振動特性進行了計算。Smirnova8等人應用類似方法對應力波的傳播進行了研究，并與傳統的分子動力學方法進行了比較，得出了誤差很小的結論(5%)。還有Pillai9等人對Bi晶體的晶面缺陷的仿真及Tadmor10-11等人對PbTiO3 晶體的磁滯現象的研究等等。引

6、言近年來，研究材料的微觀力學特性的分子動力學方法(MD)發展迅速，它建模簡單，程序短小，可計算的原子體系大大超過第一原理等方法，在解釋一些用理論分析和實驗觀測等方法都難以了解的微觀現象上起到了不小的作用。隨著計算機計算能力的不斷提高以及算法的改進，分子動力學方法可處理的原子已經數以億計1,但仍達不到仿真實際系統的要求，在時間和空間尺度上受到極大的限制。為解決這一難題，有人提出了多尺度仿真(Multiscale Modeling)，即把原子模型嵌入到連續介質模型中，采用分子動力學方法計算我們感興趣的微小區域，而其他區域采用連續介質力學方法(多為FEM)計算，不僅減小了計算量，而且使計算尺度得到了

7、極大的擴展。2 多尺度仿真基本原理原子區及連續介質區的計算方法都已相對成熟，主要難題集中在如何處理原子模型和連續介質的結合區上，對此許多學者各自采用了不同的方法。1 多尺度仿真方法的發展及應用早在上個世紀80年代初，Baskes和Mullins等人就已經運用類似的方法對晶格的裂紋擴展問題進行了仿真，但受條件制約，他們所建立的模型與真實物理條件有很大差異。后來，Kohlhoff4等人改進了他們的方法，對b.b.c晶體的裂232.1 FEAt方法Mullins3提出，在結合區域，原子與節點一一對應，直接應用原子作用力轉化成集中載荷作用于有限元節點上。這種想法簡單明了，但存在許多問題，如由于原子間作

8、用力的長程及非局部特性，很難解決區域邊界處的受力平衡。而如果用節點(原子)位移來代替力的直接作用，則可避免這一問題。由Kohlhoff等人提出的FEAt方法是這一應用的典范4。在他們提出的方法中，整個模型由四個部分組成(,)，如圖1所示。收稿日期：2005-07-25 修回日期：2005-11-28作者簡介：孫西芝(1978-), 女, 山東金鄉人, 博士生, 研究方向為納微米切削加工過程機理的建立和仿真；陳時錦(1964-), 男, 安徽肥東人, 副教授, 研究方向為超精密加工裝備的設計; 程凱(1961-), 男, 黑龍江哈爾濱人, 教授, 博導, 研究方向為3D精密與超精密表面質量和功能

9、性的控制與加工技術, e-制造技術: 模式、定量分析方法及系統實現。·2699·第18卷第10期 Vol. 18 No. 10 2006年10月 系統 仿 真 學 報 Oct., 2006過渡區，完整的硅晶格與一個有限單元相對應。為了計算晶MD格內原子的運動對單元產生的作用，引入了內位移的概念，即原子的位移是晶格位移和內位移的和，并引用Martin公式計算內位移的值。在確定了內位移的計算方法后，便得出了FEM圖1 MD區與FEM區的邊界關系仿真的具體計算步驟如下：(1) 初始化MD區的原子位置及FEM區的邊界條件。 (2) 計算MD區的原子受力，從而得到區域內的原子位移，接

10、著計算內位移。(3) 根據MD區的計算結果固定區域的節點，計算FEM區的反作用力。被固定節點位移由對應原子位移與內位移的差確定。(4) FEM區域計算，得到所有節點的位移。 (5) 根據區域的節點的位移計算內位移。(6) 根據FEM區計算結果固定區域的原子，計算MD區全部原子的位移。被固定原子的位移由對應的節點位移與內位移的和確定。(7) 檢驗收斂性，如不收斂則返回(2)進行計算。MD區域從到，FEM區域從到。,為過渡區，其中的原子與節點一一對應，區原子為FEM提供邊界條件，節點隨原子運動；區的節點為MD區域提供邊界約束，原子隨節點運動。模型采用了非局部彈性理論來描述過渡區。在線性變化范圍內，

11、晶格的彈性變形能可表示為：E=1nmik(xnxm)(uinuim)(uk) (1) uk2nm其中xn表示原子n的位置，uin表示原子n離開平衡位置的位移，ik為能量算子。其相對應的均勻各向同性連續介質的變形能可表示為：1E=ik(xx')ui(x)ui(x')uk(x)uk(x')dV'dV0,1 (2) 2vv2V0'可見，以上兩式將離散與連續兩種能量的計算聯系了起來。(2)式可變換為應變形式：E=1cijkl(xx')ij(x)kl(x')dV'dV (3) 2vv'MDFEM式中：cijkl為連續介質的二階彈性常

12、量，ij(x)為位置x處的應變。將彈性能的應變形式泰勒展開：E()=E(0)+(E/ij)0ij+12162E/(ijkl)|0ijkl+E/(ijklmn)|0ijklmn+. (4)3圖2 Si晶格與單元的對應關系2.2 QC方法QC(Quasicontunuum)方法由Tadmor和Ortiz12等人于1996年提出，其目的是構建一個完整的原子系統而不必顯式地對每一個原子進行計算。它的原理與FEAt方法有很大的不同，體現在不再是由一特殊的過渡區來連接原子區域和連續區域，而是以原子計算為基礎，通過引入有限單元概念來減小原子的計算量。(5)其中位移ui在僅原子的新位置xi可由xi=Xi+ui

13、表示，考慮其物理意義時可被看作是連續區域u(X)在Xi的值即為了保證晶格與對應單元的應力平衡，應使式中的各系數在兩種介質的計算中分別對應相等，因此做如下定義：cij=(E/ij)|0cijkl=2E/(ijkl)|0cijklmn=E/(ijklmn)|03u(Xi)=ui，今后在公式中出現的u表示這一區域內所有原子位移的集合u1,u2,.uN。為了獲得系統的總能量，我們在原子仿真過程中通常總是要對每一個原子的能量進行計算，即Etot=Ei(u) (6)i=1N這樣使連續介質的彈性常量由原子間的勢能函數定義，達到了應力平衡的目的。FEAt方法是開發較早的多尺度方法之一，它直接把離散量和連續量結

14、合了起來，體現了多尺度仿真的早期思想，并且它被應用在了一些結構簡單晶體(如b.b.c, f.c.c)的仿真上，取得了不錯的效果。但是由于方法本身的粗糙性，其對于復雜結構的晶體(如硅晶體)不再適用，因為原子與節點的一一對應關系很難保證，晶格內原子的位移給彈性模量的確定帶來困難。因此，后來的很多學者在FEAt方法的基礎上進行了新的改良，如Satashi Izumi等人的研究就是一例，他們把這種方法應用在了硅晶體的仿真上5，如圖2所示。在其中Ei表示原子i在位移u作用下的能量，它的準確性取決于勢函數的選擇。在變形相對平滑的區域對每個原子的位移都進行計算是不必要的。基于這種考慮，QC采用了近似方法，即

15、選擇一些有代表性的原子repatoms(representive atoms)計算其ui，而其他原子的位移則通過線性插補近似得到，這里恰好利用了有限元方法的插補原理，如圖3所示。·2700·第18卷第10期 Vol. 18 No. 10 2006年10月 孫西芝, 等：多尺度仿真方法研究綜述 Oct., 2006精確又簡便。一種以力學計算為基礎的QC方法也在研究之中，詳細內容見參考文獻14。QC方法已被廣泛應用在了結構破壞、微缺陷、裂紋擴展等一系列問題的處理上，是目前較完善的多尺度仿真方法，它避免了使用不夠準確的理論方法將分子動力學和有限元的直接耦合，而是將有限元概念融合到

16、了分子計算中，使圖3 repatoms的選擇方式仿真結果與真實物理條件更加接近，比FEAt方法更加先進，但是由于忽略了原子的振動，使其應用僅局限于零度下的靜態仿真，進一步研究成果見參考文獻15。在圖中，實心原子即repatoms,它們被作為節點將整個區域用三角形單元劃分，其他原子如A則受所在單元的限制。2.3 MAAD方法MAAD(Macro Atomistic Ab initio Dynamics)是由Abraham16-18等人提出的，用來仿真硅的裂紋擴展過程。該而且為方法不僅結合了分子動力學(MD)和有限元方法(FE)。了得到更加真實的效果，在裂紋產生處還采用緊束縛模型repatoms的

17、密度則根據具體問題來選擇，在我們最關心的地區所有的原子都是repatoms，而在較遠的地區，密度則有所下降。現在系統的總能量可表示為Etot=Ei(uh) (7)i=1N(Tight-Binding,TB)來仿真原子鍵的斷裂過程。為了方便原子區與有限元區的結合，有限單元被劃分到了原子級大小，如圖4所示。FEMD在此式中u由線性插補獲得，即uh=Su (8)=1NrephS是與repatom有關的插補函數，Nrep表示repatoms的數量，Nrep遠小于N。由于有限單元形狀函數的緊湊性，使上式的計算更加簡便。如在圖3中A的位移可由他所在單元的節點(B,C,D)的位移確定，即uh(XA)=SB(

18、XA)uB+SC(XA)uC+SD(XA)uD (9)TDMD上一步中自由度的數量雖然得到了減少，但在計算總能量時，我們仍然需要訪問每個原子。為了進一步降低計算量，引入了Cauchy-Born法則度不變，即(F)=E0(F)013圖4 FE/MD/TB的結合方式整個系統總能量的哈密頓函數表述為，即假設在任意單元中，變形梯Htot=HFE+HMD+HFE/MD+HTB+HMD/TB (13) 每一項代表各自區域及過渡區的能量函數，其中(10)HFE=11&)2d (14) (r)C(r)d+(r)(u220為單元的體積，E0為單元的變形能，這樣系統的總能量就可以表示為NeEtote(Fe

19、) (11)e=11Natom2&imiuHMD=+V (15)2i=1&為節點速度。 為應變量，C為剛度，為材料密度，u&i為第i個原子的速度，V代表MD區的勢能總和。在計算u原子勢能時，使用了Stillinger-Weber(SW) 勢函數19，它是一種專門為硅晶體設計的勢函數，提高了計算的準確性。Noccn=1Ne為單元總數。這時計算所有原子的能量就被替換為計算所有單元的變形能，每個單元只需計算一個原子，計算量大大降低了，式(11)被稱為局部QC公式。當變形并不趨于平滑時，例如物體缺陷處或裂紋處，局部QC法顯得不夠精確，所以作者又將局部QC法改進為非局部QC，其中

20、總能量的計算公式為EtotNrepHTB=n+Vrep(rij) (16)i<j該式可理解為原子間相互吸引項n及排斥項Vrep的共同作用，相關討論見參考文獻20。過渡區的能量（HFE/MD，HMD/TB）是由該區兩邊的能量函數的線性組合來確定，組合系數取決于各自函數對過渡區的能量貢獻。各能量函數采用相同的時間步計算，各區域之間的能量傳遞被忽略。nE(u) (12)=1式中僅出現了repatoms，所以n表示repatom所代替原子的數量。E通過對周圍原子的計算獲得。在接近裂紋或缺陷處，每個原子都被定義為repatom，使得結果更加精確。在實際應用中使局部和非局部方法相結合，使得計算過程既

21、MAAD方法為了得到更準確的仿真效果，使三種方法(FE/MD/TB)同時結合在一起。與FEAt的思想類似，這三種·2701·第18卷第10期 Vol. 18 No. 10 2006年10月 系 統 仿 真 學 報 Oct., 20063M Mullins, M Dokainish. Simulation of the (001) Plane Crack in Alpha-iron Employing A New Boundary Scheme J. Philos. Mag. A. (S0141-8610), 1982, 46: 771-787. 4S Khlhoff, P

22、Gumbsch, H F Fischmeister. Crack Propagation in b.b.c Crystals Studied with a Combined Finite-element and Atomistic Model J. Philos. Mag. A. (S0141-8610), 1991, 64A: 851-878.5 S Izumi, T Kawakami, S Sakai. Study of A Combined FEM-MDMethod for Silicon J. JSME Int. J. Ser. A. (S1344-7912), 2001, 44(1)

23、: 152-159.6 H Rafii-Tabar, L Hua, M. Cross. A Multi-scale Atomistic-continuumModeling of Crack Propagation in a Two-dimensional Macroscopic Plate J. J. Phys. Condens. Matter. (S0953-8984), 1998, 10: 2375-2387.7 R Robert. Coupling of Length Scales in MEMS Modeling the AtomicLimit of Finite Element J.

24、 Proceedings of SPIE (S0277-786X), 2000, 11: 16-25.8 J A Smirnova, L V Zhigilei, et al. Combined Molecular Dynamics andFinite Element Method Technique Applied to Laser Induced Pressure Wave Propagation J. Comput. Phys. Commun. (S0010-4655), 1999, 118(1): 11-16.9 A R Pillai, R E Miller. Crack Behavio

25、r at Bi-crystal Interfaces: AMixed Atomistic and Continuum Approach C/Mater. Res. Soc. Symp. Proc. (S0272-9172), 2001, 653: Z2.9.1-Z2.9.7.10 E B Tadmor, U V Waghmare, G.S. Smith and E. Kaxiras. PolarizationSwitching in PbTiO3: An ab initio Finite Element Simulation J. Acta Materialia (S1359-6454), 2

26、002, 50: 2989-3002. 11 E B Tadmor, R Phillips, M Ortiz. Mixed Atomistic and ContinuumModels of Deformation in Solids J. Langmuir (S0743-7463), 1996, 12: 4529-4534.12 E B Tadmor, M Ortiz, R Phillips. Quasicontinuum Analysis of 13 14 15 16Defects in Solids J. Philos. Mag. A. (S0141-8610), 1996, 73(6):

27、 1529. J L Ericksen. Phase Transformations and Material Instabilities in Solids M. New York, Academic Press, 1984.J Knap, M Ortiz. An Analysis of the Quasicontinuum Method J. J. Mech. Phys. Solid (S0022-5096), 2001, 49(9): 1899-1923.E B Tadmor, R E Miller. Overview EB/OL. , 2006.J Q Broughton, F F A

28、braham, N Bernstein, et al. Concurrent Coupling of Length Scales: Methodology and Application J. Phys. Rev. B. (S0163-1829), 1999, 60(4): 2391-2403.F F Abraham, J Q Broughton, N Bernstein. Spanning the Continuum to Quantum Length Scales in a Dynamic Simulation of Brittle Fracture J. Europhys. Lett.

29、(S0295-5075), 1998, 44(6): 783.F F Abraham, J Q Broughton, N Bernstein. Spanning the Length Scale in Dynamic Simulation J. Comput. Phys. (S0894-1866), 1998, 12: 538-546.F H Stillinger, T A Weber. Computer Simulation of Local Order in Condensed Phases of Silicon J. Phys. Rev. B. (S0163-1829), 1985, 3

30、1(8): 5262-5271.W M C Foulkes, R Haydock. Tight-binding Models and Density Functional Theory J. Phys. Rev. B. (S0163-1829), 1989, 39(17): 12520-12536.R E Rudd, J Q Broughton. Coarse-grained Molecular Dynamics and Atomic Limit of Finite Elements J. Phys. Rev. B. (S0163-1829), 1998, 58(10): 5893-5896.

31、R E Rudd. Coarse-grained Molecular Dynamics: Dissipation Due to Internal Modes C/Mater. Res. Soc. Symp. Proc. (S0272-9172), 2002, 695: 499-504.Continuum Model of Defects in Solids J. J. Mech. Phys. Solids (S0022-5096), 2002, 50(10): 2085-2106.方法間相互獨立，但它們間的結合方式比FEAt更加先進。雖然有限單元小到了原子級大小，有利于仿真效果，但同時也使計算量增大，因此方法的作者使用巨型機來進行計算。該方法主要用在晶體裂紋擴展的研究中，相關討論見參考文獻17,18。2.4 CGMD方法CGMD( Coarse Grained Molecular Dynamics)方法由Rudd和Broughton21等人提出。其基本思想是在需要詳細描而在遠離關鍵部分的區域，述的關鍵區域采用MD方法計算，引入有限元方法粗化計算，即挑選一定密度的原子作為節點，節點的能量由式(16)計算，而其他原子作簡化計算，降低