《晶格系统多尺度计算的一种新方法》由会员分享，可在线阅读，更多相关《晶格系统多尺度计算的一种新方法（7页珍藏版）》请在金锄头文库上搜索。
1、晶格系统多尺度计算的一种新方法肖瑶，唐少强( 北京大学工学院北京)麓要：本文提出了一种基于P M L 算法的多尺度模拟方法。在本文的计算工作中，将被研究的区域分成三个部分：M D ( M o l e c u l a rD y n a m i c s ) 区域、M C ( M a e r o s C o p i c ) 区域、以及P M L 区域。粗网格方程用于整个区域中，而M D 计算仅限于M D + P M L 区域。P M L 区域用于吸收M D 区域与M C 区域之间数值边界上的人为反射。对于粗网格计算，我们使用一种由匹配微分算子方法得到的方程作为粗网格的控制方程，理论上说，该粗网格方程
2、可以达到直至谱精度任意精度。本文使用该算法模拟几个一维晶格系统，取得了比较满意的结果。关键词：多尺度计算P M L 算法1 引言“连续性假设”是传统的固体力学及流体力学的一个基本假设，诸如位移、应变、应力等物理量，都是在连续介质力学体系中引入的。这一假设是基于所研究的力学现象的尺度远远大于分子或原子尺度而作出的近似假设。随着现代科学技术的不断发展，微观尺度上的物理、力学研究和应用不断加强，特别是M E M S 和纳米器件等方兴未艾，传统的连续介质力学对这些应用有较大的局限性。由于连续性假设忽略了物质微观结构的非连续性，因此传统固体力学在精确的描述一些现象( 如晶体生长、裂纹传播) 时遇到了根本
3、性困难。这时连续介质力学不再适用，需要引入量子力学、分子动力学等理论。基于量子力学或分子动力学理论建立的计算方法，又会由于计算复杂度等多方面原因，在解决实际问题时受到很大限制，因此需要建立多尺度模型将不同尺度算法耦合起来以回避单独使用某种算法所遇到的困难。本文中涉及的微观系统将通过分子动力学( 以下简称M D ) 方法来描述。多尺度计算中两个关键的问题在于，首先，在已有的M D 描述基础上，怎样获得粗尺度下的控制方程；其次，不同尺度算法所处理的区域之间存在数值边界，怎样让各种物理量在边界上平滑过度而不引入人为反射。本文将介绍一种基于P M L ( P e r f e c t l yM a t
4、c h e dL a y e r ) 算法的、实现起来比较容易的多尺度模拟方法( 我们称之为P M L M D 0 4 算法) ，并利用几个实例验证算法的有效性。该算法的基本特征是a 引入P M L 区域来处理宏观与微观尺度之间的数值边界问题；b 对粗网格的模拟，该算法中将使用S h a o q i a n gT a n g 等人提出的粗网格M D O 一4 方程。2 算法基本思想将我们所研究的区域Q 分成三个部分：M D 区域、P M L 区域、M C 区域。其中M D 区域是我们选取的特别关心的区域，通常很小，但区域上非线性很强，或者其中存在发生缺陷的原子。P M L 区域是位于M D 区
5、域与M C 区域之间的较为狭长的区域，用来吸收不同算法区域之间的人为反射。也就是 说，M D 与M C 区域并不直接接触，两者被P M L 区域间隔开。三者位置关系如图1 所示。这里，我们在区域Q 上建立粗网格，并在整个Q c 中进行粗网格计算。仅在M D 与P M L 区域中，我们需要以M D 方法计算每个原子的位移。由于M C 区域中不计算原予位移，因此在P M L M C 边界上会遇到原子位移的边条件问题。我们将通过对粗网格节点位移的插值得到这个的边条件。另外，我们将对P M L + M D 区域中的原子位移做粗细尺度分解，在P M L 区域将原子位移的粗尺度部分以粗网格节点位移的插值来
6、替换，而在M D 区域，则将原子位移的粗尺度部分赋值给粗网格节点。5 7考虑到原子振荡的时间尺度比粗网格的时间步长小得多，对粗、细尺度应当取不同的时间步长来进行计算。图1 计算区域划分示意图算法可以简单概括如下 轧对每个细时间步长1 通过对粗网格节点位移的插值，得到原子位移在P M L M C 的边界条件；2 在M D + P M L 区域更新原子位移、速度、加速度；b 。对每个粗时间步长1 在区域M C 上，更新节点位移、速度、加速度：2 在P M L 区域，以粗网格节点位移及速度的插值替换原子位移、速度的粗尺度部分；在M D 区域中，用原子变量向粗网格变量重新赋值。上述算法由A l b e
7、 r tC T o 与S h a o f a nL i 的P M M S ( P e r f e c t l yM a t c h e dM u l t i s c a l eS i m u l a t i o n s ) 算法【2 】改编得到。本文算法( 我们称为P M L - M D 0 4 算法) 与P M M S 算法的主要区别在于粗网格方程：前者使用M D O - 4 方程，而后者使用L E 方法来获取粗网格方程。这两种方程后文都会介绍。3P M L 算法P M L 算法最初由B 4 r e n g e r 提出，用以解决电磁波传播的问题3 1 。其后该算法被用于多尺度模拟，以吸收数
8、值边界上的反射。P M L 算法不仅可以很快的将任意频率( 除零以外) 的波动耗散掉，而且能够与其他的区域非常好的匹配，因此被用于解决很多不同的问题。P M L 算法同时也适用于各项异性及非均匀介质，并且实现起来相当简单。M D 区域中的原子运动满足牛顿第二定律：艺：聊i i 。：一3 U - ( r n ) ，( 1 )O l g a其中下标口表示第口个原予，F 、n 、r 分别表示原子所受外力、位移、所处位置。U 为原子问作用力的势能函数，一表示与第口个原子有相互作用的原予的位置。在笛卡尔坐标系中，不同坐标轴方向的力的分量相互独立，因此我们不妨仅讨论x 方向。于是由( 1 ) 得到吃：一掣
9、盟( 2 )聊口“口2 一i 、_ 二( 2 ) 毋口这里我们直接给出P M L 方程，其导出过程可参见 2 】，屹：一垦竽竺一d ( ) 2 u 口一2 J ( 皿渤。其中耗散函数矗( 讳) 的取值我们按C o l l i n o 和T s o g k a t 3 1 的半经验公式给出，也薯小( 去 芸，其中，艿是P M L 的厚度，R 是小于1 的自由参数，V 是波速，表示原子与M D P M L 边界的距离引入一个辅助变量，将( 3 ) 转化成与之等价的方程组形式，即可看出P M L 方程对波动的耗散作用，4 粗细尺度分解与重新赋值前面我们提到，每一个粗时间步长的计算结束时，需要对P M
10、 L 区域以及M D 区域中的变量重新赋值，在这个重新赋值的过程中，就涉及到原予位移的粗细尺度分解。我们把原子位移分解为两部分U = U + U ，( 5 )其中，U 表示原子运动中的低频部分；U 则表示原子运动中的高频部分，可以理解为原子在U 附近的细微扰动。设d 为粗网格节点位移，令u = N d ，其中N 为型函数。定义M = N r M _ N ，M 一为质量矩阵，为使得u ”M U 取得最小值，我们有d = M 一1 N r M 4 U ( 6 )令Q = I N M 一1 N F M ，则有u = Q u ，于是得到U = N d + Q u ( 7 )我们分别利用方程( 6 )
11、( 7 ) 对M D 区域中d 的以及P M L 区域中的U 重新赋值。5 租网格方程5 1M D O 4 方程这一节中，我们将介绍微分算子匹配法得到的粗网格M D O 一4 方程【1 1 。对于一维谐振子，适当选取单位尺度后，其M D 方程为5 9口坚加、p(盘彳一一。堡g ( 一蹦虹魄旦以瓦i i = 见u ，l 一2112 。 见2J1ll1- 2扛参 卜铲。卜其中P 为粗粒化系数，详细的推导过程参见【l 】。对于原子间作用力满足L e n n a r d J o n e s 势的一维振子，M D 方程为玩= _ 4 8 l ( + + 。一) 。3 一( r o + 一U n 一。)
12、。1 3+ 2 4 ( r o + + 1 - - U n ) 一( r o + 飞一。) 7 一粗网格的M D O 4 方程为 水+ 竽广卜半门 + 水+ 竽n + 掣) - 7 ，+ ( p2 _ 1 ) ( _ 5 2 r o - , , + 1 4 t o - s ) ( d ，一2 4 a , 一l + 6 4 4 田+ 1 + 西+ 2 )5 2L E 方法得到的租爵格方程L E 方法是应用广泛的有限元算法中的一种，即线性元。通过此方法，可以得到一组粗网格方程。L E 方法的基本想法时，通过对粗网格变量的插值，得到原子位移的近似值，代入M D 方程，得到原子加速度，再通过粗网格变量
13、于原子变量的关系，将原子加速度转化为粗网格节点的加速度。对于原子运动满足( 1 ) 的系统，L E 方法得到的粗网格方程为M a = 一V d U ( N d ) ( 1 2 )6 算法实现本文中我们采用V e r l e t 算法对原子的位移、速度、加速度等变量进行更新，q i + l = q ，+ f 彭+ 三( f ) 2 科，d 7 + l 2 = 茸+ 去f i i ，噜，：一：O U 一一B 2 q ，“一2 B 彭川2 ，咖q i + 1 = 一+ 1 2 - t - 1 2 耐“，( 1 3 )( 1 4 )( 1 5 )( 1 6 )其中，B 是由对应不同原子的耗散系数d (
14、 ) 构成的对角矩阵，在M D 区域中B = o 。对于更新粗网格节点的位移、速度、加速度，我们采用如下算法，：+ A t d V + i 1 ( & ) 2 a j ，每“= F ( d ，“) ，“= + 昙，( 彬+ 每+ 1 ) ，其中( 1 9 ) 为粗网格的M D O - 4 方程。6 1 一维谐振子( 1 7 )( 1 8 )( 1 9 )我们研究一个包含6 3 1 个原子的维谐振子模型，顺序编号，取第2 6 2 到3 9 0 个原子所占据的区域为M D 区域，第2 4 2 个到第2 6 1 个原子、第3 9 1 到第4 0 0 个原子所占据的区域位P M L 区域。我们取原子间
15、距离吃= 0 0 0 5 ，粗粒化倍数P = 1 0 ，粗网格节点间距离吃= 砘= O 0 5 ，于是确定粗网格节点数为6 4 ，粗时间尺度A t = 0 0 5 ，细时间尺度f = 0 0 0 5 。初条件由U n ( 0 ) = “o ( x n ) 给定，其中“。( x ) ：j 。5 爷( 1 + 0 1 c o s ( 8 。万x ) ) ，当l x l 。2 5 ；( 2 。)1 0 ，对其他情况在后面的图例中，M C 表示粗网格节点位移，M D 表示M D + P M L 区域中的原子位移，精确解则表示通过完全M D 方法计算得到的结果。图2P M L M D 0 4 算法与P
16、M M S 模拟得到的一维谐振子，P M L - M D 0 4 ：( a ) “O ，1 0 0 ) ，( b ) u ( x ，1 3 0 ) ，( c ) u ( x ，1 5 0 ) ；P M M S ：( d ) u ( x ，1 5 0 )6 1图3P M L M D 0 4 与P M M S 模拟模拟得到的L c n n a r d J o n e s 振子：P M L - M D 0 4 ：( a ) t = 1 0 ，( b ) 户2 0 ；P M M S ：( c ) t = 1 0 ，( d ) t = 2 0在该模型中，我们考察一列6 0 1 个原子，原子间距吃= 2 “6 ，处于中间位置的1 4 1 个原子所占据的区域为M D 区域，其两侧各3 0 个原子的区域为P M L 区域，余者为M C 区域，粗粒化倍数P = 1 0 ，所以全域中有粗网
