基于能量統一格式的多尺度有限元法

摘 要：為分析簡單晶體多尺度有限元計算的能量構成，利用能量最小原理得到在統一理論框架下多尺度有限元計算的統一格式，表明有限元計算可以在微觀原子尺度下和在宏觀連續介質尺度下進行多尺度有限元計算.基于簡單晶體變形的特點說明過渡單元設計應遵守的原則，并指出理想過渡單元應該是類似于晶體結構的單元，對于較復雜的晶體，則應該利用空間群方法充分研究具有230種空間群的過渡單元的性質.引用EIDEL等的納米壓痕計算結果作為算例，表明在計算中無虛擬插值點，多級晶胞單元具有與原單胞相同的點群操作且位移場插值受晶體中原子間鍵長的

約束.關鍵詞：有限元； 多尺度有限元計算； 能量最小原理； 過渡單元； 單晶體中圖分類號：O241.82；O152.8 文獻標志碼：

A

Multiscale finite element method based on unified energy frame

HUANG Junping

1, 2, PENG Xianghe1

(1. Dept. of Eng. Mechanics, Chongqing Univ., Chongqing 400044, China;2. School of Mechanical Eng., Chongqing Vocational Institute of Eng., Chongqing 400037, China)

Abstract： To analyze the energy of singlecrystal inmultiscale finite element computation, the minimized energy principle is used to obtain the unified formatunder unified theory frame. It is shown that finite element computation can be used for multiscale analysis under both microatomic scale or macrocontinuum scale. Based on the deformation characteristics of singlecrystal, a design principle for transitional elements is introduced and the optimal transitional elements should be analogous to crystal structure. For complicated crystal, space group method can be used for a further study on the character of transitional elements which contain 230 kinds of space groups. The paradigm of nanoindentation computation results obtained by EIDEL et al is cited to verify that no virtual interpolation point is required in the simulation and multilevel cell elements could have the same pointgroup symmetry operations as those of the primitive cell. Furthermore the displacement interpolation field must be constrained by the length of bonds between atoms in the crystal.Key words： finite element; multiscale finite element computation; minimized energy principle; transitional element; singlecrystal收稿日期：2010[KG*9〗08[KG*9〗17 修回日期：2010[KG*9〗10[KG*9〗25作者簡介： 黃均平(1963—)，男，四川樂山人，教授，博士研究生，研究方向為材料多尺度力學性質，(Email)jphuang_63@163.com 0 引 言

目前，多尺度模擬計算方法被廣泛應用于研究不同尺度材料的力學行為和缺陷.

［1］相對于原子模擬方法，如分子動力學方法和從頭算起的第一性原理方法，多尺度方法能更有效地研究多尺度系統，從而在原子尺度上找出材料破壞的本質特征.用分子動力學方法和第一性原理方法模擬大規模原子系統需對大量的原子區進行重復計算，占用時間多，且對計算機硬件要求高，操作系統須非常穩定，必須在集群計算機上并行運算，計算成本高昂，但計算的材料體積卻十分有限.

［23］

為降低計算成本，能在普通PC機上進行模擬計算，出現擬連續介質多尺度分析方法

［4］、跨第一原理/原子/宏觀多尺度分析方法

［56］和廣義質點動力學方法

［7］等.上述方法都建立在統一的理論框架內，均需考慮區域的過渡以達到解的平滑.但是，在晶體結構的計算中除廣義質點動力學方法考慮到晶格結構外，多數采用一維過渡單元作為計算算例用于如碳納米管的模擬計算.

［813］廣義質點動力學方法對體心立方晶格和面心立方晶格的單元級別有較好的描述，但還不能描述密集六角結構，因而不可能模擬出堆垛層錯結果.1 對晶體的能量描述

對于簡單晶格，在變形前的參考構形中，晶體中每個原子的位置坐標

［1415］都可描述為[WTHX]R[WTBX]=L1[WTHX]A[WTBX]1+L2[WTHX]A[WTBX]2+L3[WTHX]A[WTBX]3=Li[WTHX]A[WTBX]i（1）式中:([WTHX]A[WTBX]1,[WTHX]A[WTBX]2，[WTHX]A[WTBX]3)為未變形晶體晶格的逆變基矢量，其任意1組整數(L1,L2,L3)描述晶格中任意1個陣點的位置.變形后，晶體在當前構形下同一原子的坐標位置可描述為[WTHX]r[WTBX]=l1[WTHX]a[WTBX]1+l2[WTHX]a[WTBX]2+l3[WTHX]a[WTBX]3=li[WTHX]a[WTBX]i（2）式中:([WTHX]a[WTBX]1,[WTHX]a[WTBX]2,[WTHX]a[WTBX]3)為變形后晶體晶格的逆變基矢量.對于晶體中第j個原子，在參考構形和當前構形下其坐標位置分別為[WTHX]R[WTBX]j=L

由N個原子系構成的超晶結構是個多體結構體系，體系的總能量由3部分構成：

(1)由原子間形成的化學鍵所決定的晶體減聚能，記為U1，其鍵長和鍵角決定鍵能的大小，而鍵長和鍵角由2個原子間的相對位置決定，其鍵能為U1=Nk

(2)能量是由未形成化學鍵的最鄰近和次鄰近原子之間形成的作用勢能，記為U2.對于最鄰近原子之間的作用勢，可由LennardJones勢

［16］和Morse勢

Jones勢可改寫為U′2=4ε′σ′r′

jk

12-σ′r′

jk6（8）(3)由外力作用在晶體上，在原子間產生附加力和位移引起的勢能，記為U3，可表述為U3=-[WTHX]F[WTBX]-·[WTHX]u[WTBX]=-Nj=1F-j·uj[JY]（9）式中：F-j為作用在第j個原子上的外力；[WTHX]u[WTBX]^=

(u^1,

u^2,…,u^N)T為N個原子產生的位移，[WTHX]u[WTBX]^=[WTHX]r[WTBX]-[WTHX]R[WTBX]+[WTHX]r[WTBX]^[JY]（10）式中：[WTHX]r[WTBX]^為參考系原點的位移.因此，系統的總能量 UT=U1+U2+U3=ET+U3=Nk

jk

12-σr

jk6-Nj=1F-j·uj[ZK)][JY]（11）式中：ET=U1+U2.2 基于能量最小的統一格式

17）式中：[WTHX]K[WTBX]為剛度矩陣;[WTHX]P[WTBX]為非平衡力矢量；[WTHX]F[WTBX]-=

(F[WTBX]-1,F[WTBX]-2,…，F[WTBX]-N)T；[WTHX]u[WTBX]=[WTHX]r[WTBX]-[WTHX]r[WTBX]0.

對于非線性關系，通過迭代直到[WTHX]P[WTBX]=[WTHX]O[WTBX]達到平衡.因此，式（17）與有限元中的控制方程等同.由式（16）和（17）得剛度矩陣[WTHX]K[WTBX]=[WTHX]K[WTBX]Q+[WTHX]K[WTBX]

MD[JY]（18）式中：[WTHX]K[WTBX]Q為由第一性原理方法確立的剛度矩陣；[WTHX]K[WTBX]

MD為由分子動力學確立的剛度矩陣.因此，由式（17）和（18）可知，對于由原子區和分子動力學區建立的控制方程，實際上可在統一的理論框架內（基于能量最小原理）由相同的有限元格式求解.如果將式（17）右邊的非平衡力[WTHX]P[WTBX]理解為在連續介質中的非平衡力，則式（17）為涵蓋原子區、分子動力學區和連續介質區的統一多尺度有限元計算方程.對于某些材料，當相互作用勢遠小于化學鍵能時，可直接將原子區與有限元區進行“握手”計算.3 非局部作用勢的非連續介質單元

在宏觀連續介質條件下，設非局部區域Ω含有N個晶胞對應相應的局部區域，局部區域對應的應變能為E，則每個晶胞的應變能為E0=E/N，晶胞中圖 1 體心立方晶格圖 2 單胞中的彈簧個數的應變能由各原子的化學鍵分配.以體心立方晶格為例，將原子間的化學鍵模擬為剛度為k1

（<100>方向）和k2（<111>方向）的彈簧，則晶體構成見圖1.對于單個晶胞，擁有的原子個數為

8×1/8+1=2個，即每個角點原子由8個晶胞共有，每個晶胞有8個角點原子，再加上1個體心原子；彈簧個數為12×1/4+

8×1/8=4個，即考慮次鄰近原子之間的作用為3個晶格邊彈簧和1個體心彈簧，見圖2.設晶體沿<100>方向產生應變為[WTBX]ε

100，鍵長為l0，沿<111>方向產生應變為ε

111，則存儲在各彈簧中的能量分別為E

100=32k1(ε

100l0)2[JY]（19）E

111=12k2ε

11132l02=38k2ε

111l02（20）得單個晶胞的總能量 E0=E

100+E

111=32k1(ε

100l0)2+38k2ε

111l02=12(3k1)+34k2υ2(ε

100l0)2（21）式中：υ=ε

111ε

100.因此，1個具有l0邊長的晶體，單位體積能量E

v00=E0v=12l012k1+3k2υ24(ε

100)2（22）由此可得多級晶胞的能量，而這多級晶胞可構成過渡單元.式（22）表明：可基于第一性原理方法計算多級晶胞能量，從而建立基于能量最小的剛度矩陣.4 基于晶格的有向線元變形計算

根據CauchyBorn法則

［1819］，將晶體晶格視為無限小材料矢量，按變形梯度進行變換，有[WTHX]a[WTBX]=[WTHX]F[WTBX]·[WTHX]A[WTBX][JY]（23）式中：[WTHX]a[WTBX]=([WTHX]a[WTBX]1,[WTHX]a[WTBX]2,[WTHX]a[WTBX]3);[WTHX]A[WTBX]=([WTHX]A[WTBX]1,[WTHX]A[WTBX]2,[WTHX]A[WTBX]3)；[WTHX]F[WTBX]為2階變形梯度張量.由于從晶體結構觀點看，晶體的晶向矢量是有限長度的，而非無限小量，故式（19）為近似式.如果變形場是均勻的位移場，視晶體晶格為無限小材料矢量，則式（19）可認為精確.

將晶體變形前、后定義為參考構形和當前構形2個構形，并與1對點([WTHX]R[WTBX],[WTHX]r[WTBX])相對應，[WTHX]F[WTBX]為這2點張量.在參考構形下（即晶體未發生變形前），若取有向線元為d[WTHX]R[WTBX]=dLi[WTHX]A[WTBX]i，經過變形后在當前構形下有線元d[WTHX]r[WTBX]=dLi[WTHX]a[WTBX]i，將式（1）和（2）寫為協變基形式，

由于R分別為在參考構形下的剛度矩陣和非平衡力;[WTHX]F[WTBX]為變形梯度張量.

在式（29）表述的方程中，如果定義[WTHX]u[WTBX]^=[WTHX]F[WTBX]

-1[WTHX]u[WTBX][JY]（34）則除表達式與式（17）一致外，還由于剛度矩陣在微觀原子區和宏觀連續介質區表述的不變性，使得在模擬計算中不會出現交界處激波的反射

［20］，即從原子區到分子動力學區再到有限元區的平滑過渡.

將位移場與變形梯度聯系起來，這樣在不同的區域進行有限元計算，就類似于在連續介質力學條件下采用不同的單元進行計算.如果將原子區稱為非局部區域，將連續介質力學區稱為局部區域，那么在局部區和非局部區域交界處必須構建過渡單元.過渡單元不能在力學軟件的500多種單元庫中選取，而各種軟件基于二次開發的用戶子程序使單元的設計成為可能，如Abaqus/Standard中的用戶單元子程序UEL等.在簡單晶體的模擬計算中，充分利用晶格的結構特性設計過渡單元是實現多尺度有限元計算的最佳途徑，因為周期排列的簡單晶體晶胞是最理想的單元之一.5 基于點群的單元設計算例

基于晶格變形的過渡單元設計與常規有限元單元設計不同之處在于過渡單元須反映原簡單晶體的特性，對不同晶格具有更強的針對性，即對于確定晶格的過渡單元具有結構設計唯一性.用廣義質點動力學方法構建的多級晶胞具有至少存在1個點保持不動的對稱操作，滿足點群的定義，并且這樣的多級晶胞由于在不動點具有實際存在的質點，故在單元插值設計中不會出現虛擬質點，可在數學上簡化其方法.在晶體結構中，高階軸（[WTHX]C[WTBX]n,n≥3）多于1個的多軸點群晶體具有高度的對稱性，如[WTHX]T[WTBX],[WTHX]O[WTBX]和[WTHX]I[WTBX]點群，涵蓋正四面體、立方體、正八面體、正五角十二面體和正三角二十面體等.包含[WTHX]C[WTBX]n,[WTHX]D[WTBX]n,[WTHX]T[WTBX],[WTHX]O[WTBX]和[WTHX]I[WTBX]的第1類點群，由于不存在反演，其變換矩陣行列式不為負；而對于第2類點群，由于存在反演，其變換矩陣為負，單元插值可由廣義函數的奇偶性決定.

［21］

作為算例， EIDEL等

［22］以FCC鋁單晶材料進行納米壓痕模擬計算，3D塊體由64×64×64布拉菲點陣超晶組成.塊體側面原子在垂直于該面方向上固定，底面原子在z方向固定，但可在底面平面內移動.在圖3中，坐標系的軸相應于<001>方向.半徑為R的球型壓頭被模擬為外部作用勢[JB({]Vx=A·θR-r·R-r3r=x-c[JB)][JY]（35）式中：參數A為壓力強度；θ()為分步函數；R=16a0，為壓頭半徑；而c為壓頭中點位置.在模擬中

A[WTBX]=

2 000 eV/3，a0=4.032 .在計算中采用基于能量的完全非局部擬連續介質方法（QCeFNL），原子族半徑分別取為Rc=1.0a0（取樣原子數為19），

a0/2和22a0（取樣原子數為381）.由圖3

［22］可知，所構成的各級單元均為幾何相似的FCC形式，屬于具有反射操作[WTHX]σ[WTBX]的[WTHX]O[WTBX]點群.圖 3 FCC鋁單晶<001>方向截面

［22］模擬首先考慮初始弛豫問題，即在未作用壓頭情況下平衡構形主要受表面效應影響，非局部QC逐漸受到非常明顯的表面效應作用.對于處于頂角的代表原子，由于三面相交，故為高能原子.圖4為用晶格靜力學和不同的QCeFNL 方法模擬二者在初始弛豫后z方向的位移，Rc=1.0 a0原子族在z方向的收縮計算過大，其值為uz=-3.3，比用晶格靜力學方法計算的uz=1.3 大2.5倍.但是，對于

Rc=22a0，uz=-1.1 ，接近于晶格靜力學參考值.(a)晶格靜力學結果,uz=1.3 (b)QCeFNL連同雜交力校正結果(c)QCeFNL Rc=1.0a0結果，u

z min=-3.3[JP2](d)QCeFNL Rc=22a0結果，[JP]u

z min=-1.1圖 4 初始弛豫位移云圖

［22］

圖 5 晶格靜力學和QCeFNL方法力—壓深（Fh）曲線對比

［22］模擬中的另一個問題是壓頭下鋁單晶材料的位錯形核.在對各種原子族半徑的QC模擬中，力—壓深（Fh）曲線中任何半徑原子族相應的加載水平都比晶格靜力學小，見圖5.當 Rc=22a0時，與原子決定算法符合得較好.對族半徑Rc=1.0a0，取樣原子數在自適應精細化過程中，從初始48 000個原子增加到86 000個原子，原子數占晶格靜力學的8%，但QC模擬比晶格靜力學方法快8倍.

［22］圖6為由晶格靜力學方法模擬的位錯微結構，可看到4個沿111方向滑移的位錯環，即

1.根據在111平面上的塑性滑動方向，可知形核位錯為1611

2111Shockley不完全位錯.

［22］圖 6 晶格靜力學模擬位錯微結構

［22］在QC模擬中，當Rc=1.0a0時與晶格靜力學符合得較好，且在111上可觀察到4個位錯環，見圖7.用力校正法進行的高精度模擬或增加原子族半徑所得結果與晶格靜力學相比，有某些偏離.

［22］(a)Rc=1.0a0(b)Rc=2a0(c)Rc=22a0(d)雜交力校正模擬圖 7 QCeFNL模擬位錯微結構

［22］6 結束語

如果本文算例采用廣義質點法，就可構成多級單元體，而各級單元體從幾何構造上看，結構相似.

［7］KARPOV等

［23］也提出相似的算例.從晶體學觀點分析，周期性規則排列的晶體均可進行類似的質量集中，構成各種適用于不同晶體結構的過渡單元.這些單元有如下性質：(1）無虛擬插值點；(2）多級晶胞單元具有與原單胞相同的點群操作，3個四重旋轉軸3C4垂直于100面［100］，［010］和［001］，4個三重旋轉軸4C3<111>，即沿8個晶向的4個三重旋轉軸，［111］，［1

1］，沿<110>的4根二重旋轉軸4C2；(3）位移場插值受晶體中原子間鍵長的約束.

由上述分析可知，在統一的理論框架下應用能量最小原理進行多尺度分析計算，各區域的計算必須在過渡區實現無縫連接才能反映正確結果，跨尺度的計算才能光滑；而在非局部區和過渡區，需利用第一性原理方法對位移場插值進行約束.在各種多尺度有限元計算中，眾多學者雖然用各種方法構建相應的過渡單元，但基于群論方法研究具有230種空間群性質的各種晶體的過渡單元結構和性質，以及其插值形函數的數學方法是仍需要進一步努力的方向.參考文獻：

［1］ SHILKROT L E, CURTINA W A, MILLER R E. A coupled atomistic/continuum model of defects in solids［J］. J Mech & Phys Solids, 2002, 50(10): 20852106.

［2］ SHILKROT L E, MILLER R E, CURTIN W A. Multiscale plasticity modeling: coupled atomistics and discrete dislocation mechanics［J］.

J Mech & Phys Solids, 2004, 52(4): 755787.

［3］ LIU B, HUANG Y, JIANG H, et al. The atomicscale finite element method［J］. Comput Methods Appl Mech Eng, 2004, 193(1720): 18491864.

［4］ MILLER R E, TADMOR E B. The quasicontinuum method: overview, applications and current directions［J］. J ComputerAided Mat Des, 2002, 9(3): 203239.

［5］ ABRAHAM F, BROUGTON J J, BERNSTEIN N, et al. Spanning the continuum to quantum length scales in a dynamic simulation of brittle fracture［J］. Europhysics Lett, 1998, 44(6): 783787.

［6］ PARK H S, LIU Wingkam. An introduction and tutorial on multiplescale analysis in solids［J］. Comput Methods Appl Mech Eng, 2004,

193(1720): 17331772.

［7］ 范鏡泓. 材料變形與破壞的多尺度分析［M］. 北京: 科學出版社, 2008: 198.

［8］ LU Q, BHATTACHARYA B. The role of atomistic simulations in probing the smallscale aspects of fracture:a case study on a singlewalled carbon nanotube［J］. Eng Fracture Mech, 2005, 72(13): 20372071.

［9］ ZHANG P, JIANG H, HUANG Y, et al. An atomisticbased continuum theory for carbon nanotubes: analysis of fracture nucleation［J］. J Mech & Phys Solids, 2004, 52(5): 977998.

［10］ JIANG H, FENG X Q, HUANG Y, et al. Defect nucleation in carbon nanotubes under tension and torsion: StoneWales transformation［J］. Comput Methods Appl Mech Eng, 2004, 193(30): 34193429.

［11］ QIAN Dong, WAGNER G J, LIU Wingkam. A multiscale projection method for the analysis of carbon nanotubes［J］. Comput Methods Appl Mech Eng, 2004, 193(1720): 16031632.

［12］ GRIEBEL M, HAMAEKERS J. Molecular dynamics simulations of the elastic moduli of polymercarbon nanotube composites［J］. Comput Methods Appl Mech Eng, 2004, 193(1720): 17731788.

［13］ ZHANG H W, WANG J B, GUO X. Predicting the elastic properties of singlewalled carbon nanotubes［J］. J Mech & Phys Solids, 2005,

53(9): 19291950.

［14］ 黃昆, 韓汝琦. 固體物理學［M］. 北京: 高等教育出版社, 1988: 10.

［15］ 黃克智, 薛明德, 陸明萬. 張量分析［M］. 北京: 清華大學出版社, 1986: 5354. 

［16］ BASKES M I. Manybody effects in FCC metals: a LennardJones embeddedatom potential［J］. Phys Rev Lett, 1999, 83(13): 25922595.

［17］ KOMANDURI R, CHANDRASEKURAN N, RAFF L M. Molecular dynamics(MD) simulation of uniaxial tension of some singlecrystal cubic metals at nanolevel［J］. Int J Mech Sci, 2001, 43(10): 22372260.

［18］ TADMOR E B, ORTIZ M, PHILLIPS R. Quasicontinuum analysis of defects in solids［J］. Philos Mag A, 1996, 73(6): 1529.

［19］ ARROYO M, BELYTSCHKO T. An atomisticbased finite deformation membrane for single layer crystal films［J］. J Mech & Phys Solids, 2002, 50(9): 19411977.

［20］ HAROLD S. P, LIU W K. An introduction and tutorial on multiplescale analysis in solids［J］. Comput Methods Appl Mech Eng, 2004,

193(1720): 17331772.

［21］ 陶瑞寶. 物理學中的群論［M］. 上海: 上?？茖W技術出版社, 1986: 142166.

［22］ EIDEL B, STUKOWSKI A. A variational formulation of the quasicontinuum method based on energy sampling in clusters［J］. J Mech & Phys Solids, 2009, 57(1): 87108.

［23］ KARPOV E G, YU H, PARK H S, et al. Multiscale boundary conditions in crystalline solids: theory and application to nanoindentation［J］. Int J Solids & Structures, 2006, 43(21): 63596379.(編輯 陳鋒杰)