阻尼色散波的人工邊界條件①

黃良易 王賢明

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

0 引言

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

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

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

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

1 人工邊界條件的引入

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

圖1 阻尼單原子鏈模型

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

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

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

1.1 色散關(guān)系

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

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

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

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

這種單色波的頻率ω與波數(shù)k之間需要滿足的表達(dá)式,稱為色散關(guān)系[41]。

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

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

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

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

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

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

圖2 阻尼單原子鏈的頻率-波數(shù)關(guān)系

2 人工邊界條件的構(gòu)造

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

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

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

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

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

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

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

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

表1 阻尼單原子鏈MBCN?的系數(shù)

3 反射系數(shù)分析

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

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

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

其中,反射系數(shù)R(ω) 為頻率ω的入射波在邊界引起的反射波的相對波幅。將總波場式(12)代入邊界條件式(8)中,可得:

式中的色散關(guān)系和頻域殘差函數(shù)由式(5)和(11)給出。

圖3 為阻尼單原子鏈匹配邊界條件(MBCN?)的反射系數(shù)。圖中,橫坐標(biāo)為頻率ω,縱坐標(biāo)為反射系數(shù)的模|R|。從圖中可以看到,MBCN?的反射系數(shù)在總體上都小于1。這意味著所提人工邊界條件可以很好地抑制行波的反射。此外,如前所述,所提人工邊界條件可完全吸收任意匹配的頻率的入射波,這一特點可從對數(shù)坐標(biāo)圖(圖3(b))中反射系數(shù)的向下尖峰上體現(xiàn)出來。對于涉及更多原子數(shù)目的邊界條件MBC2?(0,0.5)、MBC3?(0,0.5,1)、MBC4?(0,0.5,1,1.5),隨著匹配頻率數(shù)目的增加,能夠在寬頻范圍上一致地減小反射系數(shù)。

圖3 MBCN?反射系數(shù)

總之,反射系數(shù)分析的結(jié)果驗證了所提人工邊界條件可有效抑制阻尼波的反射。

4 數(shù)值算例驗證

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

4.1 單波包數(shù)值算例

本節(jié)對一條無限長阻尼單原子鏈進(jìn)行數(shù)值模擬。在數(shù)值模擬中,選取包含81 個原子的片段進(jìn)行計算求解,其原子編號為-40～40,原子質(zhì)量為1,彈簧的剛度系數(shù)為1,阻尼系數(shù)為0.1。左、右兩側(cè)原子u-40和u40處施加人工邊界條件。每次計算時,左、右兩端施加相同的邊界條件MBCN?。作為對比,本算例還求解了一條很好地近似了原始無限域的足夠長片段的數(shù)值解,截取其-40 號到40 號的部分原子作為參考解。本算例采用中心差分法進(jìn)行數(shù)值計算,時間步長Δt取1/128。

本算例取初始靜止、如下形式的波包作為初始條件。

其中,波包主頻率為1.4,對應(yīng)波數(shù)實部k取值1.5508。

下面來觀察單波包算例的數(shù)值結(jié)果。圖4 為兩側(cè)邊界施加MBC3?(0,0.5,1)時,不同時刻數(shù)值解與參考解的對比。各子圖分別對應(yīng)t=0、t=25、t=45 時刻。在t=0 時刻,初始波包位于參考解中心位置,如圖4(a)所示。圖4(b)中,中心波包分解成左、右兩部分波包向兩端傳播,向左、右端傳播的波包逐漸變寬。由于阻尼單原子鏈具有色散性,長波的主峰傳播速度較快,較短的短波緊隨其后。與此同時,由于阻尼的存在,向左右兩端傳播的波包幅值不斷減小。此時,所有解的波形與參考解吻合得很好。隨著左、右波包到達(dá)邊界后,波包逐漸被MBC3?吸收,圖4(c)中計算中心位置波包已完全分解,左、右波包的整體寬度增加,波包的幅值明顯減小,整體數(shù)值解與參考解非常吻合?？梢?MBC3?能有效吸收阻尼色散波。

圖4 兩邊界施加MBC3?,在不同時刻下的數(shù)值解

4.2 雙波包算例

為了進(jìn)一步考察所提邊界條件在較寬頻率范圍上的吸收效果,取初始靜止、如下形式的中心波包作為初始條件進(jìn)行計算。

它包含了2 個頻率的波包疊加,第1 個波包主頻率ω為1,對應(yīng)的波數(shù)實部k1為1.0472;第2 個波包主頻率ω為1.5,對應(yīng)的波數(shù)實部k2為1.6952。在此算例中,選取包含201 個原子的片段,阻尼系數(shù)為0.01。

圖5 為兩邊界采用MBC3?(0,0.5,1)在t=0、60、150 時刻下的數(shù)值解。如圖5(a)所示,初始波包在阻尼單原子鏈的中心位置,中心波包由2 個不同主頻率的波包疊加而成。隨著時間的推移,中心波包分成左右兩部分進(jìn)行傳播。由于阻尼單原子鏈的色散性,中心波包中不同主頻率波包的傳播速度不同,中心波包分成了2 個傳播較快的長波波包和2 個傳播較慢的短波波包,如圖5(b)所示。隨著波包逐步到達(dá)邊界后,被MBC3?吸收。從圖5(c)中可見,數(shù)值解與參考解吻合得很好。對于寬頻域的波傳播問題,MBCN?仍具有良好的吸收效果。

圖5 兩邊界施加MBC3?,在不同時刻下的數(shù)值解

圖6 為t=150 時刻下,不同邊界條件下的誤差。MBC2?的數(shù)值誤差明顯可見,部分行波的反射已經(jīng)影響到內(nèi)部-50、50 號原子附近。隨著匹配的階數(shù)的提高,MBC3?和MBC4?的數(shù)值反射顯著減小?？梢?隨著匹配階數(shù)的提高,MBCN?的吸收效果增強(qiáng)。

圖6 不同人工邊界條件下的數(shù)值誤差

5 結(jié)論

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

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

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