阻尼色散波的人工边界条件①

黄良易 王贤明

(浙江工业大学特种装备制造与先进加工技术教育部重点实验室 杭州 310014)

0 引言

近几十年来,随着计算机能力和数值算法的发展,研究者大量地采用计算机数值模拟来研究材料和结构的力学性能和行为。一方面,碳纳米管、纳米半导体材料等新兴材料层出不穷,新的材料结构使其往往具有新的力学性能。另一方面,在非常规极端环境下服役的材料及其结构的力学行为更加复杂。这两个趋势推动着材料力学性能和行为的微观机理研究,人们大量地使用分子动力学模拟研究此类问题。实际的物理系统通常具有庞大的原子数目,采用完全的原子模拟计算量巨大。因此,在数值模拟中通常采用截断原子计算区域的办法来减小计算量与存储量。这样就引入了截断区域的计算边界,这种人为划定的计算边界,往往不能随意构造,需要针对分子动力学的运动特征设计边界条件。这种边界条件被称为人工边界条件(artificial boundary conditions)[1-4]。

不准确的人工边界条件往往会产生显著的数值反射,从而影响计算的结果,甚至会带来数值不稳定性而导致计算结果发散。波动问题的人工边界条件,也被称为吸收边界条件(arbsorbing boundary conditions)[1-3,5-7]。对于无界域连续介质波传播问题的人工边界条件,已有大量的研究[1,3,8-12]。由于晶格动力学的离散和色散特性,这些方法大多数无法直接应用,还需要进一步发展[13-15]。例如,文献[13]和[16]分别将连续介质波传播的完美匹配层方法[3]推广到分子动力学模拟中。

对于晶体的原子模拟,目前也已经发展了多种人工边界条件[6,11-12,17-32]。最初由文献[17]提出了分子动力学模拟的一种精确人工边界条件。文献[21,23-26,28,29,33,34]系统地将这种人工边界条件推广到处理二维和三维晶格,称之为时间历史积分方法(time history kernel treatment),并以之为基础发展出桥接多尺度方法(bridging scale method)。时间历史积分方法有以下局限性:一方面,该方法依赖核函数,而核函数的解析求解困难,尤其是对于多维问题不易实现[35],通常采用数值方法[18,26];另一方面,时间历史积分的数值卷积计算量较大,为了提高求解效率,通常进行时间截断,而这又带来了误差[30]。最后,基于基本解的时间历史积分在时间和空间上都是非局部的,不易推广应用到非线性晶体的原子模拟中。因此,除了全局人工边界条件之外,人们发展出了精度稍低、计算量小、便于应用的局部人工边界条件。以文献[1]为基础,通过最小化反射系数的方式,文献[11,12,36-37]构造了晶格动力学的局部人工边界条件。通过对离散晶格动力学方程进行泰勒展开、算子分裂的方式,文献[27]构造了一种主要处理长波的多层的速度界面条件。考虑到晶格动力学色散波动特征,文献[31]和[38]提出了基于匹配色散关系的匹配边界条件(matching boundary conditions,MBCs)。匹配边界条件方法采用边界附近原子速度-位移的线性组合形式,组合系数通过匹配晶格行波的色散关系确定。其结构紧凑,构造简单,吸收效果好,已运用到简单晶格[39-40]以及复式晶格[31]中。

本文旨在设计出能够有效处理阻尼色散波的局部人工边界条件。阻尼是材料的固有属性,在研究材料的动力学性能时考虑阻尼是必要的。离散晶格中波传播具有色散性,晶体中不同频率波的传播速度不同。由于阻尼的存在,使得晶体中的波传播具有边传播边衰减的特性。其人工边界条件构造的难点在于要同时处理不同波速、不同空间衰减率的阻尼色散波。针对晶体中的阻尼色散波,本文将基于匹配边界条件方法[31]通过匹配以频率而非波数为自变量建立的阻尼波色散关系,来构造人工边界条件,并利用反射系数分析和数值算例来验证所提人工边界条件的有效性。

1 人工边界条件的引入

本文以如图1 所示的阻尼一维单原子链为例,考察阻尼色散波的波传播问题。阻尼一维单原子链由集中质量、弹簧以及阻尼元件组合而成。其中,每个原子质量为m,相邻原子之间采用弹簧系数为k的弹簧相连,各个原子独立地与阻尼相连,阻尼系数为c。

图1 阻尼单原子链模型

在数值仿真中,由于无法直接求解无限长链,需要进行区域截断。图1 中,原无限长链的原子编号为j=-∞,…,-2,-1,0,1,2,…,+∞,在计算中只求解实心表示的右侧原子,左侧虚线表示的部分则被截去而不求解。截断后,0 号原子的时间演化无法按原动力学方程求解,需要构造人工边界条件。人工边界条件的构造并不是任意的,最理想的情况是能反映出原问题中被截去的左侧半无限长链的动力学运动特征。

若只考虑最近邻相互作用,在无外力作用下,第j个原子的归一化运动方程为

为了模拟其动力学特性,匹配边界条件方法[30-31]采用匹配原子链的色散关系来构造人工边界条件。本文首先研究阻尼单原子链的色散关系。

1.1 色散关系

考察满足动力学方程式(1)的如下形式的稳态解:

式中,ω为无量纲角频率,k=kR+ikI为复数值的无量纲波数。角频率一般取非负值,本文简称频率。将上式带入到动力学方程式(1)中,可得到特征方程:

将上式实部与虚部分离,可得:

以频率ω为自变量,即可得到:

这种单色波的频率ω与波数k之间需要满足的表达式,称为色散关系[41]。

下面对满足阻尼单原子链运动方程的稳态解式(2)作一个物理解释。uj(t)=Aei(ωt+kRj) e-kIj所表示的运动情况为:无限长链中所有原子都以同一个频率ω做周期振动,原子振幅随空间指数衰减;在同一时刻,相邻原子间的相位差为kR,且右侧原子的相位滞后于其左侧相邻原子。也就是说,在上述波数和频率取正值的规定下,uj(t)=Aei(ωt+kRj) e-kIj表示了右行的衰减单色波。反之,uj(t)=Aei(ωt+kRj) ekIj则表示了左行的衰减单色波。

倘若按常规的以波数为自变量建立色散关系,则可求得复数值的频率为

其中,波数是实数值。这样给出的特征波解uj(t)=Aei(ωRt+kj)e-ωIt表示的运动情况是简谐波形中的原子振幅在空间上分布大小相同,在时间上以同一衰减率衰减。

这2 种方式表示的特征波不同。在人工边界条件设计时,所求解问题中的扰动波源一般位于所关心的截断后保留的有限计算区域内部。相比于以波数为自变量的色散关系,以频率为自变量的色散关系所表示的阻尼波物理特征更吻合数值计算的实际情况。采用后者,更便于人工边界条件的有效处理。

采用频率作为自变量的色散关系式(5),能够方便地表示出阻尼波的传播和空间衰减特性。此外,对于截止频率以上的倏逝波区域,无阻尼单原子链与阻尼单原子链的波数皆为复数形式。因此,采用频率作为自变量表示色散关系更统一。

图2为不同阻尼系数下阻尼单原子链的色散关系。各子图分别对应阻尼系数c=0,0.01,0.1,1。图中,横坐标为频率,纵坐标为波数。ω∈[0,2] 为行波区域,ω=2 为截止频率,ω∈[2,4] 为倏逝波区域。实线表示波数实部,虚线表示波数虚部。与无阻尼链色散关系图2(a)不同,阻尼链色散关系在行波区域,存在非零的波数虚部;在倏逝波区域,阻尼链的波数实部随频率变化,其相邻原子之间的相位差随频率改变,而无阻尼链的波数实部为定值,其相邻原子之间的相位差保持恒定不变;此外,由于阻尼的存在,阻尼链波数虚部的绝对值更大。阻尼系数越大,波数虚部绝对值越大,倏逝波振幅的空间衰减越快。

图2 阻尼单原子链的频率-波数关系

2 人工边界条件的构造

本文利用上述色散关系,构造能同时处理不同波速、不同空间衰减率阻尼波的人工边界条件。以图1 所示的左侧人工边界为例,演示人工边界条件的构造方法。原无限长阻尼单原子链被截断后,在左侧边界原子u0处,本文采用如式(8)所示的边界条件。

该边界条件利用了内部N个原子位移和速度的线性组合。待定系数bj和cj通过匹配前述色散关系式(5)和(6)来确定。

作为不失一般性的示例,本文用以下条件:

确定边界条件式(8)中的待定系数。其中,定义了匹配边界条件的“频域残差函数”[30]

这里,色散关系由式(5)和(6)给出。由文献[30,38]引入的频域残差函数,改以频率作自变量,表示了人工边界条件对于原无限长链色散关系和阻抗特性的逼近程度。对于特定频率ω=ωb,若有Δ(ωb,k(ωb))=0,这意味着该频率相应的原子链的特征稳态简谐解uj=Aei(ωbt+jk(ωb)),既满足运动方程,同时也满足本文所提的人工边界条件式(8)。于是,此特征波解也是满足包含了人工边界条件的截断区域问题的精确解。在这种情况下,相应频率的波被边界条件完全吸收而无任何反射。

条件式(9)和(10)在一共N+1 个特殊频率处匹配了色散关系。其中,在0 频率处有1 个条件,在其余特殊频率处各有方程的实部、虚部2 个条件。于是,上述条件的数目刚好可以用来确定2N+1 个待定系数(给定c0=1)。

由此确定的人工边界条件,是文献[30]中匹配多个波数的泰勒-牛顿型匹配边界条件的推广。在文献[30,31]中,匹配边界条件通过同时匹配多个波数的行波的色散关系来构造。但是,以波数为自变量的色散关系无法反映具有复波数的阻尼行波。本文采用前述以频率为自变量、包含波数实部和虚部的完备形式的色散关系,通过“匹配多个频率”,构造能同时有效处理若干个行波的人工边界条件。为了区别于通过“匹配波数”的匹配边界条件[30],本文采用了“匹配频率”的方式构造处理阻尼波的人工边界条件。特别地,将由条件式(9)和(10)所确定的人工边界条件,记为MBCN∗(0,ω1,ω2,…,ωN)。其中,附加星号∗,以与之前匹配波数的匹配边界条件相区别;N表示涉及到的内部原子数;括号内为匹配的各个频率。

总之,对于阻尼单原子链,本文采用了以频率为自变量的统一色散关系,由此可以方便地构造人工边界条件。例如,对于阻尼系数c=0.01 的阻尼单原子链,表1 给出了MBC2∗(0,0.01,0.5)、MBC3∗(0,0.01,0.5,1)和MBC4∗(0,0.01,0.5,1,1.5)的系数。

表1 阻尼单原子链MBCN∗的系数

3 反射系数分析

当数值求解无限域阻尼单原子链的波动问题时,原始无限域模型被包含人工边界的有限区域模型所替代。模型替代后所带来的数值误差可视为人工边界处的虚假的反射波。一般情况下,波在介质界面或表面处的反射现象,常用反射系数加以分析[11,30-31]。

下面通过反射系数来分析上述人工边界条件的吸收效果。

对于左端人工边界条件,考虑任意频率的阻尼波入射下的稳态总波场:

其中,反射系数R(ω) 为频率ω的入射波在边界引起的反射波的相对波幅。将总波场式(12)代入边界条件式(8)中,可得:

式中的色散关系和频域残差函数由式(5)和(11)给出。

图3 为阻尼单原子链匹配边界条件(MBCN∗)的反射系数。图中,横坐标为频率ω,纵坐标为反射系数的模|R|。从图中可以看到,MBCN∗的反射系数在总体上都小于1。这意味着所提人工边界条件可以很好地抑制行波的反射。此外,如前所述,所提人工边界条件可完全吸收任意匹配的频率的入射波,这一特点可从对数坐标图(图3(b))中反射系数的向下尖峰上体现出来。对于涉及更多原子数目的边界条件MBC2∗(0,0.5)、MBC3∗(0,0.5,1)、MBC4∗(0,0.5,1,1.5),随着匹配频率数目的增加,能够在宽频范围上一致地减小反射系数。

图3 MBCN∗反射系数

总之,反射系数分析的结果验证了所提人工边界条件可有效抑制阻尼波的反射。

4 数值算例验证

本节通过数值算例来验证所提人工边界条件的有效性。

4.1 单波包数值算例

本节对一条无限长阻尼单原子链进行数值模拟。在数值模拟中,选取包含81 个原子的片段进行计算求解,其原子编号为-40～40,原子质量为1,弹簧的刚度系数为1,阻尼系数为0.1。左、右两侧原子u-40和u40处施加人工边界条件。每次计算时,左、右两端施加相同的边界条件MBCN∗。作为对比,本算例还求解了一条很好地近似了原始无限域的足够长片段的数值解,截取其-40 号到40 号的部分原子作为参考解。本算例采用中心差分法进行数值计算,时间步长Δt取1/128。

本算例取初始静止、如下形式的波包作为初始条件。

其中,波包主频率为1.4,对应波数实部k取值1.5508。

下面来观察单波包算例的数值结果。图4 为两侧边界施加MBC3∗(0,0.5,1)时,不同时刻数值解与参考解的对比。各子图分别对应t=0、t=25、t=45 时刻。在t=0 时刻,初始波包位于参考解中心位置,如图4(a)所示。图4(b)中,中心波包分解成左、右两部分波包向两端传播,向左、右端传播的波包逐渐变宽。由于阻尼单原子链具有色散性,长波的主峰传播速度较快,较短的短波紧随其后。与此同时,由于阻尼的存在,向左右两端传播的波包幅值不断减小。此时,所有解的波形与参考解吻合得很好。随着左、右波包到达边界后,波包逐渐被MBC3∗吸收,图4(c)中计算中心位置波包已完全分解,左、右波包的整体宽度增加,波包的幅值明显减小,整体数值解与参考解非常吻合。可见,MBC3∗能有效吸收阻尼色散波。

图4 两边界施加MBC3∗,在不同时刻下的数值解

4.2 双波包算例

为了进一步考察所提边界条件在较宽频率范围上的吸收效果,取初始静止、如下形式的中心波包作为初始条件进行计算。

它包含了2 个频率的波包叠加,第1 个波包主频率ω为1,对应的波数实部k1为1.0472;第2 个波包主频率ω为1.5,对应的波数实部k2为1.6952。在此算例中,选取包含201 个原子的片段,阻尼系数为0.01。

图5 为两边界采用MBC3∗(0,0.5,1)在t=0、60、150 时刻下的数值解。如图5(a)所示,初始波包在阻尼单原子链的中心位置,中心波包由2 个不同主频率的波包叠加而成。随着时间的推移,中心波包分成左右两部分进行传播。由于阻尼单原子链的色散性,中心波包中不同主频率波包的传播速度不同,中心波包分成了2 个传播较快的长波波包和2 个传播较慢的短波波包,如图5(b)所示。随着波包逐步到达边界后,被MBC3∗吸收。从图5(c)中可见,数值解与参考解吻合得很好。对于宽频域的波传播问题,MBCN∗仍具有良好的吸收效果。

图5 两边界施加MBC3∗,在不同时刻下的数值解

图6 为t=150 时刻下,不同边界条件下的误差。MBC2∗的数值误差明显可见,部分行波的反射已经影响到内部-50、50 号原子附近。随着匹配的阶数的提高,MBC3∗和MBC4∗的数值反射显著减小。可见,随着匹配阶数的提高,MBCN∗的吸收效果增强。

图6 不同人工边界条件下的数值误差

5 结论

本文设计了针对阻尼色散波的人工边界条件。本文采用匹配边界条件方法,通过匹配色散关系构造了匹配边界条件。这是一类可以有效处理宽频带上阻尼色散波的局部人工边界条件。反射系数分析与数值算例结果都验证了所提人工边界条件的有效性。

晶体动力学人工边界处理的主要困难之一是晶格波动的色散性,而阻尼进一步增加了处理的复杂性。晶格中的阻尼波具有边传播边衰减的特征,其空间衰减率不仅取决于阻尼系数,还与波的频率有关。基于阻尼单原子链的运动特征,本文通过以频率而非波数为自变量建立了阻尼单原子链的色散关系,可以方便地反映出阻尼波复数值的波数,更具一般性。在此基础上,发展了匹配边界条件方法,使它能够处理一般的晶格波。

总之,本文所提人工边界条件形式紧凑、构造方便、匹配频带宽、吸收效果好、数值稳定性好,具有良好的适用性。作为一种高效、实用的局部人工边界条件,可以推广到多维晶格的动力学问题计算中。