摘 要：为分析简单晶体多尺度有限元计算的能量构成，利用能量最小原理得到在统一理论框架下多尺度有限元计算的统一格式，表明有限元计算可以在微观原子尺度下和在宏观连续介质尺度下进行多尺度有限元计算.基于简单晶体变形的特点说明过渡单元设计应遵守的原则，并指出理想过渡单元应该是类似于晶体结构的单元，对于较复杂的晶体，则应该利用空间群方法充分研究具有230种空间群的过渡单元的性质.引用EIDEL等的纳米压痕计算结果作为算例，表明在计算中无虚拟插值点，多级晶胞单元具有与原单胞相同的点群操作且位移场插值受晶体中原子间键长的
约束.关键词：有限元； 多尺度有限元计算； 能量最小原理； 过渡单元； 单晶体中图分类号：O241.82；O152.8 文献标志码：
A
Multiscale finite element method based on unified energy frame
HUANG Junping
1, 2, PENG Xianghe1
(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 singlecrystal inmultiscale 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 multiscale analysis under both microatomic scale or macrocontinuum scale. Based on the deformation characteristics of singlecrystal, 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 nanoindentation computation results obtained by EIDEL et al is cited to verify that no virtual interpolation point is required in the simulation and multilevel cell elements could have the same pointgroup 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; multiscale finite element computation; minimized energy principle; transitional element; singlecrystal收稿日期：2010[KG*9〗08[KG*9〗17 修回日期：2010[KG*9〗10[KG*9〗25作者简介： 黄均平(1963—)，男，四川乐山人，教授，博士研究生，研究方向为材料多尺度力学性质，(Email)jphuang_63@163.com 0 引 言
目前，多尺度模拟计算方法被广泛应用于研究不同尺度材料的力学行为和缺陷.
［1］相对于原子模拟方法，如分子动力学方法和从头算起的第一性原理方法，多尺度方法能更有效地研究多尺度系统，从而在原子尺度上找出材料破坏的本质特征.用分子动力学方法和第一性原理方法模拟大规模原子系统需对大量的原子区进行重复计算，占用时间多，且对计算机硬件要求高，操作系统须非常稳定，必须在集群计算机上并行运算，计算成本高昂，但计算的材料体积却十分有限.
［23］
为降低计算成本，能在普通PC机上进行模拟计算，出现拟连续介质多尺度分析方法
［4］、跨第一原理/原子/宏观多尺度分析方法
［56］和广义质点动力学方法
［7］等.上述方法都建立在统一的理论框架内，均需考虑区域的过渡以达到解的平滑.但是，在晶体结构的计算中除广义质点动力学方法考虑到晶格结构外，多数采用一维过渡单元作为计算算例用于如碳纳米管的模拟计算.
［813］广义质点动力学方法对体心立方晶格和面心立方晶格的单元级别有较好的描述，但还不能描述密集六角结构，因而不可能模拟出堆垛层错结果.1 对晶体的能量描述
对于简单晶格，在变形前的参考构形中，晶体中每个原子的位置坐标
［1415］都可描述为[WTHX]R[WTBX]=L1[WTHX]A[WTBX]1+L2[WTHX]A[WTBX]2+L3[WTHX]A[WTBX]3=Li[WTHX]A[WTBX]i（1）式中:([WTHX]A[WTBX]1,[WTHX]A[WTBX]2，[WTHX]A[WTBX]3)为未变形晶体晶格的逆变基矢量，其任意1组整数(L1,L2,L3)描述晶格中任意1个阵点的位置.变形后，晶体在当前构形下同一原子的坐标位置可描述为[WTHX]r[WTBX]=l1[WTHX]a[WTBX]1+l2[WTHX]a[WTBX]2+l3[WTHX]a[WTBX]3=li[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)由原子间形成的化学键所决定的晶体减聚能，记为U1，其键长和键角决定键能的大小，而键长和键角由2个原子间的相对位置决定，其键能为U1=Nk (2)能量是由未形成化学键的最邻近和次邻近原子之间形成的作用势能，记为U2.对于最邻近原子之间的作用势，可由LennardJones势 ［16］和Morse势 Jones势可改写为U′2=4ε′σ′r′ jk 12-σ′r′ jk6（8）(3)由外力作用在晶体上，在原子间产生附加力和位移引起的势能，记为U3，可表述为U3=-[WTHX]F[WTBX]-·[WTHX]u[WTBX]=-Nj=1F-j·uj[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]^为参考系原点的位移.因此，系统的总能量 UT=U1+U2+U3=ET+U3=Nk jk 12-σr jk6-Nj=1F-j·uj[ZK)][JY]（11）式中：ET=U1+U2.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，则每个晶胞的应变能为E0=E/N，晶胞中图 1 体心立方晶格图 2 单胞中的弹簧个数的应变能由各原子的化学键分配.以体心立方晶格为例，将原子间的化学键模拟为刚度为k1 （<100>方向）和k2（<111>方向）的弹簧，则晶体构成见图1.对于单个晶胞，拥有的原子个数为 8×1/8+1=2个，即每个角点原子由8个晶胞共有，每个晶胞有8个角点原子，再加上1个体心原子；弹簧个数为12×1/4+ 8×1/8=4个，即考虑次邻近原子之间的作用为3个晶格边弹簧和1个体心弹簧，见图2.设晶体沿<100>方向产生应变为[WTBX]ε 100，键长为l0，沿<111>方向产生应变为ε 111，则存储在各弹簧中的能量分别为E 100=32k1(ε 100l0)2[JY]（19）E 111=12k2ε 11132l02=38k2ε 111l02（20）得单个晶胞的总能量 E0=E 100+E 111=32k1(ε 100l0)2+38k2ε 111l02=12(3k1)+34k2υ2(ε 100l0)2（21）式中：υ=ε 111ε 100.因此，1个具有l0边长的晶体，单位体积能量E v00=E0v=12l012k1+3k2υ24(ε 100)2（22）由此可得多级晶胞的能量，而这多级晶胞可构成过渡单元.式（22）表明：可基于第一性原理方法计算多级晶胞能量，从而建立基于能量最小的刚度矩阵.4 基于晶格的有向线元变形计算 根据CauchyBorn法则 ［1819］，将晶体晶格视为无限小材料矢量，按变形梯度进行变换，有[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]=dLi[WTHX]A[WTBX]i，经过变形后在当前构形下有线元d[WTHX]r[WTBX]=dLi[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-r3r=x-c[JB)][JY]（35）式中：参数A为压力强度；θ()为分步函数；R=16a0，为压头半径；而c为压头中点位置.在模拟中 A[WTBX]= 2 000 eV/3，a0=4.032 .在计算中采用基于能量的完全非局部拟连续介质方法（QCeFNL），原子族半径分别取为Rc=1.0a0（取样原子数为19）， a0/2和22a0（取样原子数为381）.由图3 ［22］可知，所构成的各级单元均为几何相似的FCC形式，属于具有反射操作[WTHX]σ[WTBX]的[WTHX]O[WTBX]点群.图 3 FCC铝单晶<001>方向截面 ［22］模拟首先考虑初始弛豫问题，即在未作用压头情况下平衡构形主要受表面效应影响，非局部QC逐渐受到非常明显的表面效应作用.对于处于顶角的代表原子，由于三面相交，故为高能原子.图4为用晶格静力学和不同的QCeFNL 方法模拟二者在初始弛豫后z方向的位移，Rc=1.0 a0原子族在z方向的收缩计算过大，其值为uz=-3.3，比用晶格静力学方法计算的uz=1.3 大2.5倍.但是，对于 Rc=22a0，uz=-1.1 ，接近于晶格静力学参考值.(a)晶格静力学结果,uz=1.3 (b)QCeFNL连同杂交力校正结果(c)QCeFNL Rc=1.0a0结果，u z min=-3.3[JP2](d)QCeFNL Rc=22a0结果，[JP]u z min=-1.1图 4 初始弛豫位移云图 ［22］ 图 5 晶格静力学和QCeFNL方法力—压深（Fh）曲线对比 ［22］模拟中的另一个问题是压头下铝单晶材料的位错形核.在对各种原子族半径的QC模拟中，力—压深（Fh）曲线中任何半径原子族相应的加载水平都比晶格静力学小，见图5.当 Rc=22a0时，与原子决定算法符合得较好.对族半径Rc=1.0a0，取样原子数在自适应精细化过程中，从初始48 000个原子增加到86 000个原子，原子数占晶格静力学的8%，但QC模拟比晶格静力学方法快8倍. ［22］图6为由晶格静力学方法模拟的位错微结构，可看到4个沿111方向滑移的位错环，即 1.根据在111平面上的塑性滑动方向，可知形核位错为1611 2111Shockley不完全位错. ［22］图 6 晶格静力学模拟位错微结构 ［22］在QC模拟中，当Rc=1.0a0时与晶格静力学符合得较好，且在111上可观察到4个位错环，见图7.用力校正法进行的高精度模拟或增加原子族半径所得结果与晶格静力学相比，有某些偏离. ［22］(a)Rc=1.0a0(b)Rc=2a0(c)Rc=22a0(d)杂交力校正模拟图 7 QCeFNL模拟位错微结构 ［22］6 结束语 如果本文算例采用广义质点法，就可构成多级单元体，而各级单元体从几何构造上看，结构相似. ［7］KARPOV等 ［23］也提出相似的算例.从晶体学观点分析，周期性规则排列的晶体均可进行类似的质量集中，构成各种适用于不同晶体结构的过渡单元.这些单元有如下性质：(1）无虚拟插值点；(2）多级晶胞单元具有与原单胞相同的点群操作，3个四重旋转轴3C4垂直于100面［100］，［010］和［001］，4个三重旋转轴4C3<111>，即沿8个晶向的4个三重旋转轴，［111］，［1 1］，沿<110>的4根二重旋转轴4C2；(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): 20852106. ［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): 755787. ［3］ LIU B, HUANG Y, JIANG H, et al. The atomicscale finite element method［J］. Comput Methods Appl Mech Eng, 2004, 193(1720): 18491864. ［4］ MILLER R E, TADMOR E B. The quasicontinuum method: overview, applications and current directions［J］. J ComputerAided Mat Des, 2002, 9(3): 203239. ［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): 783787. ［6］ PARK H S, LIU Wingkam. An introduction and tutorial on multiplescale analysis in solids［J］. Comput Methods Appl Mech Eng, 2004, 193(1720): 17331772. ［7］ 范镜泓. 材料变形与破坏的多尺度分析［M］. 北京: 科学出版社, 2008: 198. ［8］ LU Q, BHATTACHARYA B. The role of atomistic simulations in probing the smallscale aspects of fracture:a case study on a singlewalled carbon nanotube［J］. Eng Fracture Mech, 2005, 72(13): 20372071. ［9］ ZHANG P, JIANG H, HUANG Y, et al. An atomisticbased continuum theory for carbon nanotubes: analysis of fracture nucleation［J］. J Mech & Phys Solids, 2004, 52(5): 977998. ［10］ JIANG H, FENG X Q, HUANG Y, et al. Defect nucleation in carbon nanotubes under tension and torsion: StoneWales transformation［J］. Comput Methods Appl Mech Eng, 2004, 193(30): 34193429. ［11］ QIAN Dong, WAGNER G J, LIU Wingkam. A multiscale projection method for the analysis of carbon nanotubes［J］. Comput Methods Appl Mech Eng, 2004, 193(1720): 16031632. ［12］ GRIEBEL M, HAMAEKERS J. Molecular dynamics simulations of the elastic moduli of polymercarbon nanotube composites［J］. Comput Methods Appl Mech Eng, 2004, 193(1720): 17731788. ［13］ ZHANG H W, WANG J B, GUO X. Predicting the elastic properties of singlewalled carbon nanotubes［J］. J Mech & Phys Solids, 2005, 53(9): 19291950. ［14］ 黄昆, 韩汝琦. 固体物理学［M］. 北京: 高等教育出版社, 1988: 10. ［15］ 黄克智, 薛明德, 陆明万. 张量分析［M］. 北京: 清华大学出版社, 1986: 5354. ［16］ BASKES M I. Manybody effects in FCC metals: a LennardJones embeddedatom potential［J］. Phys Rev Lett, 1999, 83(13): 25922595. ［17］ KOMANDURI R, CHANDRASEKURAN N, RAFF L M. Molecular dynamics(MD) simulation of uniaxial tension of some singlecrystal cubic metals at nanolevel［J］. Int J Mech Sci, 2001, 43(10): 22372260. ［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 atomisticbased finite deformation membrane for single layer crystal films［J］. J Mech & Phys Solids, 2002, 50(9): 19411977. ［20］ HAROLD S. P, LIU W K. An introduction and tutorial on multiplescale analysis in solids［J］. Comput Methods Appl Mech Eng, 2004, 193(1720): 17331772. ［21］ 陶瑞宝. 物理学中的群论［M］. 上海: 上海科学技术出版社, 1986: 142166. ［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): 87108. ［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): 63596379.(编辑 陈锋杰)
