黃良易 王賢明

(浙江工業(yè)大學特種裝備制造與先進加工技術教育部重點實驗室 杭州 310014)

0 引言

近幾十年來,隨著計算機能力和數值算法的發(fā)展,研究者大量地采用計算機數值模擬來研究材料和結構的力學性能和行為。一方面,碳納米管、納米半導體材料等新興材料層出不窮,新的材料結構使其往往具有新的力學性能。另一方面,在非常規(guī)極端環(huán)境下服役的材料及其結構的力學行為更加復雜。這兩個趨勢推動著材料力學性能和行為的微觀機理研究,人們大量地使用分子動力學模擬研究此類問題。實際的物理系統(tǒng)通常具有龐大的原子數目,采用完全的原子模擬計算量巨大。因此,在數值模擬中通常采用截斷原子計算區(qū)域的辦法來減小計算量與存儲量。這樣就引入了截斷區(qū)域的計算邊界,這種人為劃定的計算邊界,往往不能隨意構造,需要針對分子動力學的運動特征設計邊界條件。這種邊界條件被稱為人工邊界條件(artificial boundary conditions)[1-4]。

不準確的人工邊界條件往往會產生顯著的數值反射,從而影響計算的結果,甚至會帶來數值不穩(wěn)定性而導致計算結果發(fā)散。波動問題的人工邊界條件,也被稱為吸收邊界條件(arbsorbing boundary conditions)[1-3,5-7]。對于無界域連續(xù)介質波傳播問題的人工邊界條件,已有大量的研究[1,3,8-12]。由于晶格動力學的離散和色散特性,這些方法大多數無法直接應用,還需要進一步發(fā)展[13-15]。例如,文獻[13]和[16]分別將連續(xù)介質波傳播的完美匹配層方法[3]推廣到分子動力學模擬中。

對于晶體的原子模擬,目前也已經發(fā)展了多種人工邊界條件[6,11-12,17-32]。最初由文獻[17]提出了分子動力學模擬的一種精確人工邊界條件。文獻[21,23-26,28,29,33,34]系統(tǒng)地將這種人工邊界條件推廣到處理二維和三維晶格,稱之為時間歷史積分方法(time history kernel treatment),并以之為基礎發(fā)展出橋接多尺度方法(bridging scale method)。時間歷史積分方法有以下局限性:一方面,該方法依賴核函數,而核函數的解析求解困難,尤其是對于多維問題不易實現[35],通常采用數值方法[18,26];另一方面,時間歷史積分的數值卷積計算量較大,為了提高求解效率,通常進行時間截斷,而這又帶來了誤差[30]。最后,基于基本解的時間歷史積分在時間和空間上都是非局部的,不易推廣應用到非線性晶體的原子模擬中。因此,除了全局人工邊界條件之外,人們發(fā)展出了精度稍低、計算量小、便于應用的局部人工邊界條件。以文獻[1]為基礎,通過最小化反射系數的方式,文獻[11,12,36-37]構造了晶格動力學的局部人工邊界條件。通過對離散晶格動力學方程進行泰勒展開、算子分裂的方式,文獻[27]構造了一種主要處理長波的多層的速度界面條件。考慮到晶格動力學色散波動特征,文獻[31]和[38]提出了基于匹配色散關系的匹配邊界條件(matching boundary conditions,MBCs)。匹配邊界條件方法采用邊界附近原子速度-位移的線性組合形式,組合系數通過匹配晶格行波的色散關系確定。其結構緊湊,構造簡單,吸收效果好,已運用到簡單晶格[39-40]以及復式晶格[31]中。

本文旨在設計出能夠有效處理阻尼色散波的局部人工邊界條件。阻尼是材料的固有屬性,在研究材料的動力學性能時考慮阻尼是必要的。離散晶格中波傳播具有色散性,晶體中不同頻率波的傳播速度不同。由于阻尼的存在,使得晶體中的波傳播具有邊傳播邊衰減的特性。其人工邊界條件構造的難點在于要同時處理不同波速、不同空間衰減率的阻尼色散波。針對晶體中的阻尼色散波,本文將基于匹配邊界條件方法[31]通過匹配以頻率而非波數為自變量建立的阻尼波色散關系,來構造人工邊界條件,并利用反射系數分析和數值算例來驗證所提人工邊界條件的有效性。

1 人工邊界條件的引入

本文以如圖1 所示的阻尼一維單原子鏈為例,考察阻尼色散波的波傳播問題。阻尼一維單原子鏈由集中質量、彈簧以及阻尼元件組合而成。其中,每個原子質量為m,相鄰原子之間采用彈簧系數為k的彈簧相連,各個原子獨立地與阻尼相連,阻尼系數為c。

圖1 阻尼單原子鏈模型

在數值仿真中,由于無法直接求解無限長鏈,需要進行區(qū)域截斷。圖1 中,原無限長鏈的原子編號為j=-∞,…,-2,-1,0,1,2,…,+∞,在計算中只求解實心表示的右側原子,左側虛線表示的部分則被截去而不求解。截斷后,0 號原子的時間演化無法按原動力學方程求解,需要構造人工邊界條件。人工邊界條件的構造并不是任意的,最理想的情況是能反映出原問題中被截去的左側半無限長鏈的動力學運動特征。

若只考慮最近鄰相互作用,在無外力作用下,第j個原子的歸一化運動方程為

為了模擬其動力學特性,匹配邊界條件方法[30-31]采用匹配原子鏈的色散關系來構造人工邊界條件。本文首先研究阻尼單原子鏈的色散關系。

1.1 色散關系

考察滿足動力學方程式(1)的如下形式的穩(wěn)態(tài)解:

式中,ω為無量綱角頻率,k=kR+ikI為復數值的無量綱波數。角頻率一般取非負值,本文簡稱頻率。將上式帶入到動力學方程式(1)中,可得到特征方程:

將上式實部與虛部分離,可得:

以頻率ω為自變量,即可得到:

這種單色波的頻率ω與波數k之間需要滿足的表達式,稱為色散關系[41]。

下面對滿足阻尼單原子鏈運動方程的穩(wěn)態(tài)解式(2)作一個物理解釋。uj(t)=Aei(ωt+kRj) e-kIj所表示的運動情況為:無限長鏈中所有原子都以同一個頻率ω做周期振動,原子振幅隨空間指數衰減;在同一時刻,相鄰原子間的相位差為kR,且右側原子的相位滯后于其左側相鄰原子。也就是說,在上述波數和頻率取正值的規(guī)定下,uj(t)=Aei(ωt+kRj) e-kIj表示了右行的衰減單色波。反之,uj(t)=Aei(ωt+kRj) ekIj則表示了左行的衰減單色波。

倘若按常規(guī)的以波數為自變量建立色散關系,則可求得復數值的頻率為

其中,波數是實數值。這樣給出的特征波解uj(t)=Aei(ωRt+kj)e-ωIt表示的運動情況是簡諧波形中的原子振幅在空間上分布大小相同,在時間上以同一衰減率衰減。

這2 種方式表示的特征波不同。在人工邊界條件設計時,所求解問題中的擾動波源一般位于所關心的截斷后保留的有限計算區(qū)域內部。相比于以波數為自變量的色散關系,以頻率為自變量的色散關系所表示的阻尼波物理特征更吻合數值計算的實際情況。采用后者,更便于人工邊界條件的有效處理。

采用頻率作為自變量的色散關系式(5),能夠方便地表示出阻尼波的傳播和空間衰減特性。此外,對于截止頻率以上的倏逝波區(qū)域,無阻尼單原子鏈與阻尼單原子鏈的波數皆為復數形式。因此,采用頻率作為自變量表示色散關系更統(tǒng)一。

圖2為不同阻尼系數下阻尼單原子鏈的色散關系。各子圖分別對應阻尼系數c=0,0.01,0.1,1。圖中,橫坐標為頻率,縱坐標為波數。ω∈[0,2] 為行波區(qū)域,ω=2 為截止頻率,ω∈[2,4] 為倏逝波區(qū)域。實線表示波數實部,虛線表示波數虛部。與無阻尼鏈色散關系圖2(a)不同,阻尼鏈色散關系在行波區(qū)域,存在非零的波數虛部;在倏逝波區(qū)域,阻尼鏈的波數實部隨頻率變化,其相鄰原子之間的相位差隨頻率改變,而無阻尼鏈的波數實部為定值,其相鄰原子之間的相位差保持恒定不變;此外,由于阻尼的存在,阻尼鏈波數虛部的絕對值更大。阻尼系數越大,波數虛部絕對值越大,倏逝波振幅的空間衰減越快。

圖2 阻尼單原子鏈的頻率-波數關系

2 人工邊界條件的構造

本文利用上述色散關系,構造能同時處理不同波速、不同空間衰減率阻尼波的人工邊界條件。以圖1 所示的左側人工邊界為例,演示人工邊界條件的構造方法。原無限長阻尼單原子鏈被截斷后,在左側邊界原子u0處,本文采用如式(8)所示的邊界條件。

該邊界條件利用了內部N個原子位移和速度的線性組合。待定系數bj和cj通過匹配前述色散關系式(5)和(6)來確定。

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

確定邊界條件式(8)中的待定系數。其中,定義了匹配邊界條件的“頻域殘差函數”[30]

這里,色散關系由式(5)和(6)給出。由文獻[30,38]引入的頻域殘差函數,改以頻率作自變量,表示了人工邊界條件對于原無限長鏈色散關系和阻抗特性的逼近程度。對于特定頻率ω=ωb,若有Δ(ωb,k(ωb))=0,這意味著該頻率相應的原子鏈的特征穩(wěn)態(tài)簡諧解uj=Aei(ωbt+jk(ωb)),既滿足運動方程,同時也滿足本文所提的人工邊界條件式(8)。于是,此特征波解也是滿足包含了人工邊界條件的截斷區(qū)域問題的精確解。在這種情況下,相應頻率的波被邊界條件完全吸收而無任何反射。

條件式(9)和(10)在一共N+1 個特殊頻率處匹配了色散關系。其中,在0 頻率處有1 個條件,在其余特殊頻率處各有方程的實部、虛部2 個條件。于是,上述條件的數目剛好可以用來確定2N+1 個待定系數(給定c0=1)。

由此確定的人工邊界條件,是文獻[30]中匹配多個波數的泰勒-牛頓型匹配邊界條件的推廣。在文獻[30,31]中,匹配邊界條件通過同時匹配多個波數的行波的色散關系來構造。但是,以波數為自變量的色散關系無法反映具有復波數的阻尼行波。本文采用前述以頻率為自變量、包含波數實部和虛部的完備形式的色散關系,通過“匹配多個頻率”,構造能同時有效處理若干個行波的人工邊界條件。為了區(qū)別于通過“匹配波數”的匹配邊界條件[30],本文采用了“匹配頻率”的方式構造處理阻尼波的人工邊界條件。特別地,將由條件式(9)和(10)所確定的人工邊界條件,記為MBCN?(0,ω1,ω2,…,ωN)。其中,附加星號?,以與之前匹配波數的匹配邊界條件相區(qū)別;N表示涉及到的內部原子數;括號內為匹配的各個頻率。

總之,對于阻尼單原子鏈,本文采用了以頻率為自變量的統(tǒng)一色散關系,由此可以方便地構造人工邊界條件。例如,對于阻尼系數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 反射系數分析

當數值求解無限域阻尼單原子鏈的波動問題時,原始無限域模型被包含人工邊界的有限區(qū)域模型所替代。模型替代后所帶來的數值誤差可視為人工邊界處的虛假的反射波。一般情況下,波在介質界面或表面處的反射現象,常用反射系數加以分析[11,30-31]。

下面通過反射系數來分析上述人工邊界條件的吸收效果。

對于左端人工邊界條件,考慮任意頻率的阻尼波入射下的穩(wěn)態(tài)總波場:

其中,反射系數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 數值算例驗證

本節(jié)通過數值算例來驗證所提人工邊界條件的有效性。

4.1 單波包數值算例

本節(jié)對一條無限長阻尼單原子鏈進行數值模擬。在數值模擬中,選取包含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 結論

本文設計了針對阻尼色散波的人工邊界條件。本文采用匹配邊界條件方法,通過匹配色散關系構造了匹配邊界條件。這是一類可以有效處理寬頻帶上阻尼色散波的局部人工邊界條件。反射系數分析與數值算例結果都驗證了所提人工邊界條件的有效性。

晶體動力學人工邊界處理的主要困難之一是晶格波動的色散性,而阻尼進一步增加了處理的復雜性。晶格中的阻尼波具有邊傳播邊衰減的特征,其空間衰減率不僅取決于阻尼系數,還與波的頻率有關。基于阻尼單原子鏈的運動特征,本文通過以頻率而非波數為自變量建立了阻尼單原子鏈的色散關系,可以方便地反映出阻尼波復數值的波數,更具一般性。在此基礎上,發(fā)展了匹配邊界條件方法,使它能夠處理一般的晶格波。

總之,本文所提人工邊界條件形式緊湊、構造方便、匹配頻帶寬、吸收效果好、數值穩(wěn)定性好,具有良好的適用性。作為一種高效、實用的局部人工邊界條件,可以推廣到多維晶格的動力學問題計算中。