一種適合基于原子尺度有限元方程的低通濾波邊界條件

任國武， 湯鐵鋼， 李慶忠

（中國工程物理研究院流體物理研究所，四川綿陽 621999）

一種適合基于原子尺度有限元方程的低通濾波邊界條件

任國武， 湯鐵鋼， 李慶忠

（中國工程物理研究院流體物理研究所，四川綿陽 621999）

發(fā)展一種基于原子尺度的有限元動力學方程，其動態(tài)行為跨越了從原子尺度到宏觀尺度.理論計算該方程的色散關系和動力學散射特征.在此基礎上，基于濾波器的概念，設計低通濾波邊界條件用于減少網(wǎng)格增大帶來的高頻反射波，同時又保證低頻振動波傳播不受影響.通過一維數(shù)值模擬，計算了能量反射和透射系數(shù)，展示低通濾波邊界條件可以吸收高頻反射波而不影響低頻波的傳播.

有限元；濾波；反射系數(shù)

0 引言

多尺度模擬是一個復雜而非?；钴S的研究領域，材料斷裂、湍流現(xiàn)象等研究引起了廣泛的關注.過去十多年間，不同研究者發(fā)展了多種多尺度方法［1］，跨越了從原子尺度到宏觀尺度，在探索材料的斷裂、位錯等現(xiàn)象［2－12］的物理機制等方面取得了許多富有性的成果.多尺度方法的思想是在包含有非線性變形的區(qū)域采用精細的模擬方法，如第一性原理［4，10］或分子動力學，而在外圍變形較小的區(qū)域采用近似的模擬方法，如有限元方法.然而，不同模擬尺度上物理模型的差異必然要帶來耦合的問題，從而造成一些非物理的現(xiàn)象發(fā)生，比如虛假波反射［13］.

當前主要的目標是發(fā)展有限溫度多尺度模擬方法模擬動力學問題，由于耦合邊界造成的聲子反射，從而造成邊界區(qū)存在溫度梯度，以至整個耦合體系很難達到熱力學平衡.因此在耦合邊界考慮有效的邊界條件用于吸收這些虛假的高頻反射很有必要.已有一些方法［14－21］解決這個問題，如精確邊界條件、吸收邊界條件等等.但是這些方法整體上是吸收來源于原子區(qū)域的所有振動模，包括在有限元區(qū)域體現(xiàn)長程相互作用的低頻波.因而在整個模擬過程中還是限制于準靜態(tài)過程，難以建立帶有熱力學耦合的多尺度模擬方法.

本工作中，我們發(fā)展一種基于原子尺度的有限元方法，利用發(fā)展的低通濾波邊界條件來解決網(wǎng)格增加造成的虛假反射問題.數(shù)值模擬證實該方法能有效消除高頻的反射波而不影響低頻波傳播.

1 模型及方法

我們建立的動力學有限元模型，如圖1所示.考慮到初始整個模擬體系由無窮多個原子構(gòu)成，可采用經(jīng)典的Lagrangian力學來描寫.在不考慮外力的作用下，該體系的Lagrangian量等于其動能減去勢能，表示為

其中u是原子位移場，MA是原子質(zhì)量mμ的對角矩陣，V（u）是原子體系的多體相互作用勢.

現(xiàn)在把不同長度的單元分配到整個原子體系，由于外載作用下有些區(qū)域出現(xiàn)非線性變形，甚至有斷裂鍵，采用細單元來描述，而有些區(qū)域則變形很小，可包含多個原子的大單元，由此整個原子體系的自由度就被減少.圖1清晰展示了該體系局部某個區(qū)域內(nèi)有些單元長度剛好為原子間距離，有些單元則包含多個原子.對于圖1所示的一維模型，我們采用兩節(jié)點的單元.為方便理論推導，每個單元的長度hl是原子距離aμ的整數(shù)倍，即hl＝nlaμ．依據(jù)動力學約束條件，每個單元內(nèi)原子的位移場可以表示為u＝Jd，d是節(jié)點的位移場，J是一插值函數(shù).利用該條件能實現(xiàn)整個體系自由度的減少.通過插值計算得到單元內(nèi)所有原子的位移場.

圖1 一維有限元模型示意圖Fig.1 Schematic illustration of one-dimensional FE model

把原子和節(jié)點位移場關系代入式（1）后，求解在受約束后的有限元動力學方程，

其中M和K分別是質(zhì)量矩陣和剛度矩陣，

從式（3）中可看出，質(zhì)量矩陣與離散原子質(zhì)量有關，剛度矩陣與原子間相互作用勢能有關.因此，我們稱式（2）為基于原子尺度的動力學方程（atomic-based finite element method，簡寫為AFEM）.

為了推導式（2）對應于圖1有限元模型的具體形式，J選擇為線性插值函數(shù)，一維情形下在節(jié)點k處為

其中，xk為節(jié)點k的位置，它是一個局域的函數(shù)，在節(jié)點處的值為1，而在次近鄰節(jié)點以及次近鄰單元之外其值為0.因此我們將利用J的定義來推導AFEM的質(zhì)量和剛度矩陣形式，并與傳統(tǒng)有限元求得矩陣作比較.

通過理論推導首先獲得集成的AEFM質(zhì)量矩陣形式為

其中，ζl＝（nl－1）（nl＋1）／6nl，ηl＝nl－2ζl，是一個三對角的質(zhì)量矩陣，矩陣元值依賴于單元長度與原子間距比值nl.

同時計算了相應的剛度矩陣，從形式可以看出，它是依賴于原子體系的勢函數(shù)，因而為了方便理論推導，選取的勢函數(shù)只考慮了最近鄰和次近鄰的簡諧相互作用

其中f1和f2分別為最近鄰和次近鄰的力常數(shù).上面給出的剛度矩陣表達式具體寫為

式（6）代入式（7）得到集成剛度矩陣為

由此得到了AFEM的質(zhì)量和剛度矩陣，與網(wǎng)格單元的長度有關.當單元長度與原子間距離相等時，質(zhì)量矩陣只含有對角元，剛度矩陣為五對角元，這與原子體系下獲得動力學方程是相同的.當單元的長度為無窮大時，質(zhì)量矩陣對應為連續(xù)情形下的一致質(zhì)量矩陣，剛度矩陣退化為連續(xù)情形，此時AFEM運動方程過渡為連續(xù)彈性理論離散后的有限元方程.因此我們發(fā)展的一維AFEM模型建立了從原子區(qū)域過渡到連續(xù)區(qū)域，在靠近原子長度的單元，節(jié)點的相互作用展現(xiàn)了非局域的特性，節(jié)點的相互作用除了最近鄰外，還包括了次近鄰的作用；單元長度趨向無窮大時，此時節(jié)點作用表現(xiàn)為局域性，具有與傳統(tǒng)有限元方法相同的性質(zhì).

2 AFEM方程的動力學行為

2.1 色散關系

基于理論推導得到的AFEM方程，可計算動力學色散關系.單元大小都相同且長度為hl，每個單元中包含的原子數(shù)目為nl個.在引入正則的平面波解di＝d0ei（qx－ωt）代入運動方程（2）后，得到對應的動力學方程.

所以，AFEM模型的頻譜為

為了比較，我們也求得了在此條件下的傳統(tǒng)有限元方法，由集中質(zhì)量矩陣近似和連續(xù)彈性常數(shù)描述的動力學方程，

其中cμ＝（f1＋4f2）／nl．注意上式中的dμ為節(jié)點位移.對應的頻譜表示為

原子模型的運動方程為

對應的聲子譜是

在圖2，我們給出以上三種情況的頻譜.計算中我們選擇了有限元區(qū)域內(nèi)每個單元的原子數(shù)目為20.在考慮周期性邊界條件后，其相應Brillouin邊界的波矢為π／20.由圖中比較可以看出，在長波區(qū)域，也就是波矢q靠近0的區(qū)域，這三個方法得到頻譜曲線相近.但是在靠近Brillouin邊界時，AFEM模型計算的譜相對于原子模型的誤差是低于傳統(tǒng)有限元方法的，表明AFEM包括更多的與原子區(qū)域相同的振動模，并且AFEM模型在Brillouin內(nèi)所得到的頻譜是優(yōu)于傳統(tǒng)有限元方法.

2.2 反射系數(shù)

接下來分析AFEM方程的散射問題.由于單元長度的變化，在能量守恒的有限元區(qū)域，必然會有反射的發(fā)生，短波模不允許在整個區(qū)域中傳播，因此這些非物理的反射波在動力學模擬過程中會造成局域的熱，在界面產(chǎn)生溫度梯度，而這些熱會讓整個體系處于非平衡的狀態(tài).我們將用散射方法來定量分析這個問題，即計算反射系數(shù)和透射系數(shù).

為了方便計算，只考慮單元增大一倍的情形.入射波在細單元的區(qū)域，在節(jié)點k處單元增加，在界面出發(fā)生反射，一部份會繼續(xù)傳播到粗單元，而一部分會反射回細單元，因此在細單元和粗單元的位移為平面波形式

圖2 AFEM、原子模型以及集中質(zhì)量近似的有限元模型色散關系Fig.2 Dispersion relations of AFEM，atomic model and FE model with lumped mass approximation

反射系數(shù)和透射系數(shù)則表達成R＝｜r｜2，T＝｜t｜2．波矢q與頻率的關系經(jīng)由動力學矩陣得到.至于系數(shù)A只是表示波的振幅大小，而反射系數(shù)和透射系數(shù)是個比值，因此可選擇A等于1.雖然位移場是由多個頻率組合而成，但由于方程是簡諧性的，所以只考慮單頻就可以了.把式（14）代入AFEM的運動方程（2）得到

式（15）是代表節(jié)點的運動方程.變換uj（t）為

注意節(jié)點的質(zhì)量只有最近鄰有作用，而節(jié)點的力作用項除了最近鄰外還包括有次近鄰作用.代入式（14）有

式（19）需要注意的是M項，只有在｜k｜＜2才不為零.表面上看，解決上式的v項需要求解所有的節(jié)點數(shù)才行，但是由于在遠離單元突變的區(qū)域外，這些節(jié)點的v項是相等的，只是存在相位的差別.利用數(shù)值解得到相應的反射系數(shù).

圖3給出了隨波矢變化的反射系數(shù)，同時也計算了集中質(zhì)量近似下傳統(tǒng)有限元方法的結(jié)果.從圖可以看出，在采用AFEM后反射系數(shù)降低了很多.AFEM發(fā)生全反射時的最低頻率為π／2，而傳統(tǒng)有限元方法的最低頻率為π／6.因此很顯然，采用AFEM后意味著有更多的低頻模能傳播到粗單元的區(qū)域，降低了聲子散射.然而，大量的高頻振動模仍然被反射了，其攜帶的聲子能量反射回細網(wǎng)格區(qū)域，造成這個區(qū)域的溫度局域身高，從而導致整個包含有細和粗網(wǎng)格區(qū)域間出現(xiàn)能量梯度，體系處于非平衡熱力學狀態(tài).這是一種非物理的熱能，需要被消除.為此我們發(fā)展了一個基于低通濾波邊界條件的方法用于耗散這些高頻反射波.

圖3 AFEM和集中質(zhì)量近似的有限元模型反射系數(shù)Fig.3 Reflection coefficients of AFEM and FE with lumped mass approximation

3 低通濾波邊界條件

濾波器方法在信號處理中有著廣泛的應用，用來處理不同頻率信號波的傳播，在實際應用中加入一些適合的電子元件（如RLC電路），就可以實現(xiàn)哪些波通過，哪些波被消除，這種方法的理論指導就是在輸入信號和輸出信號之間構(gòu)建一個傳輸函數(shù).它的作用就是對輸入信號進行過濾，得到需要的輸出信號.這個方法一直被應用于電路中，近些年來，已經(jīng)有學者把它應用到了計算模擬之中，獲得生物分子中有用的原子振動信息.

3.1 理想傳輸函數(shù)

在動力學模擬中，每個原子的振動是隨時間變化的，它們也可以被看成是一個信號，稱為原子信號.原子的相互作用就是不同原子信號間的相互傳播.因此在原子系統(tǒng)可構(gòu)造在原子信號間傳遞的傳輸函數(shù).

考慮最簡單的一維模型，原子與原子間的相互作用只考慮最近鄰簡諧勢來描述，其中力常數(shù)為k1，運動方程為

其中ω0

2＝k1／mμ，mμ為原子的質(zhì)量．對于n個原子的體系，這樣對應的方程有n個．在以下的推導中，令k1和mμ為單位1.經(jīng)由Laplace變換，方程（20）可變換為

對于無窮長的原子鏈，且其中原子不受外力作用下，可近似θj＝θj＋1，得到

在頻域空間，H（iω）的模為1，表明所有的振動?？蔁o阻來回傳播.注意到雖然H（s）的推導是基于無窮長原子鏈，但從定義看出它描述局域的近鄰原子間振動傳輸關系.

3.2 低通濾波耗散邊界條件

以上我們得到了理想原子模型情況下的傳輸函數(shù)，希望在這個模型的基礎上設計一個與原子信號相關的低通濾波器來處理提出的虛假反射問題.所設計的濾波器依賴于AFEM得到的反射系數(shù)，目的是過濾掉在波矢大于π／nl時的全反射，而讓低頻的通過.在此情況中，通過色散關系變換反射系數(shù)與波矢的關系為反射系數(shù)與聲子振動頻率的關系.需要注意的是發(fā)生全反射時對應的頻率為ωt.所以待設計的低通濾波器的要求是過濾掉頻率大于ωt的波，而不影響低于這個頻率波的傳播.

通過一系列的設計和分析，得到滿足要求的低通濾波器，其數(shù)學形式為

事實上，在得到這個結(jié)果前，我們已經(jīng)做了其他耗散方法的低通特性測試，比如Lagevin耗散，獲得其傳輸函數(shù)，在引入多個耗散層的條件下，得到了整體傳輸函數(shù)隨頻率的關系，給出其濾波特性，結(jié)果證實了添加一定數(shù)目的耗散項后，所有來源于原子區(qū)域的波都能被吸收掉.通過得到的傳輸函數(shù)，模擬中能定量分析如何加入耗散層過濾掉所有反射的波.

在頻域空間，傳輸函數(shù)H（s）模隨頻率的變化如圖4.在圖4中，虛線是理想的濾波器設計，也就是如果能設計出大于ωt的波被完全過濾的濾波器，而不會影響低頻波的傳播，那將是最完美的，但是由于物理模型的限制，實際設計的濾波器由實線表示，顯而易見它是一個低通的濾波器，雖然它也過濾掉了一些不該濾掉的低頻波，但是對高頻的過濾基本可滿足要求，大大減少高頻的反射.對我們而言，最為關心的是傳輸函數(shù)的模，事實上在實際模擬中，它對應于能量多少，而相位作用會帶來色散效應，在這里暫時略去這方面的分析.其中的參數(shù)β雖說是一個可調(diào)參數(shù)，但是在分析中發(fā)現(xiàn)，它的大小不會影響濾波的特性，仍保持了低通濾波器的性質(zhì).在下面我們將給出詳細的說明.

由低通濾波的數(shù)學形式，基于Laplace逆變換獲得對應時空域下的運動方程，表示為

圖4 傳輸函數(shù)與振動頻率關系（虛線為理想的傳輸函數(shù)，實線為依賴實際原子模型的傳輸函數(shù).）Fig.4 Schematic illustration of transfer function and vibrated frequency（The dashed line is an ideal transfer function while the solid line is a realistic one relying on atomic model.）

從等式（27）發(fā)現(xiàn)運動方程中多余的部分是近鄰原子的速度二階差分，由于這個方法是過濾高頻的聲子，而不影響低頻聲子傳播，因此這個方法稱為聲子濾波器方法（PFM）.從運動方程中看出，這個低通濾波耗散方法只是比Lagevin耗散多了兩個次近鄰項，造成對低頻的耗散影響變小.參數(shù)β的選擇是非常重要的，在實際的應用中，可選擇多個β項來耗散高頻的反射波，一般而言，最好選擇其值逐漸增加的原則，如果值太大，會出現(xiàn)過衰減的情況，反而會增強反射.

圖5 波矢K＝0.1和K＝1.0的波包在一維模型中的傳播，（a）集中質(zhì)量近似的有限元方法，（b）AFEM，（c）AFEM加濾波耗散方法Fig.5 Wave-packet propagation of wave vectors of K＝0.1 and K＝1.0 in one-dimensional model，（a）lumped mass FE；（b）AFEM；（c）AFEM with PF method

通過理論的分析，以上得到的低通濾波耗散方法滿足了要求，雖然與理想的情況差別很大，但是從本質(zhì)上，這個傳輸函數(shù)展示了低通性質(zhì).而β項的選擇，我們只是做了簡單地分析.事實上，可以通過數(shù)值優(yōu)化β來定出所需的傳輸函數(shù)滿足理想情況，達到完全過濾掉高頻波而不影響低頻波的傳播.

4 一維數(shù)值模擬結(jié)果

我們將實施一維有限元動力學模擬，模型如圖1.在最左邊是均勻大小的粗單元，然后往右逐漸變小至細單元，長度為原子間距.最長單元的長度是原子間距的8倍.在實施模擬開始時，整個區(qū)域的質(zhì)量矩陣和剛度矩陣需要提前生成，而且在以后都不發(fā)生變化.速度－Verlet算法用于整個模型的時間積分.

在模擬過程中，我們將采用Gaussian波包作為測試信號，其形式

其中K是平均波矢.把這個波包作為初始位移加入細單元區(qū)域，x0是波包的中心位置，σ波包的寬度.

為了比較，模擬條件選擇了三種：①傳統(tǒng)有限元方法，集中質(zhì)量近似；②AFEM；③AFEM＋聲子濾波方法.采用了兩個不同的波矢K＝0.1和K＝1.0作為測試輸入.測試的結(jié)果如圖5所示，對于圖5中的每個子圖，虛線表示的是網(wǎng)格逐漸增大的過渡區(qū)域，左邊為粗網(wǎng)格區(qū)域，右邊為細網(wǎng)格區(qū)域.每個圖示的上下兩個波包圖，對應于兩個時刻用于說明波包在從初始的原子區(qū)域以及經(jīng)過過渡區(qū)域后到達的粗單元區(qū)域.首先分析波矢K＝0.1時，它位于圖形的左邊.從圖可以看出在采用AFEM后，低頻的波包完全地進入了粗網(wǎng)格區(qū)域，而條件（1）的模擬方法，仍然有部分被返回.同時也發(fā)現(xiàn)引入低通濾波耗散方法后，對這個波矢影響很小.對于圖形的右邊波矢K＝1.0的模擬，如果不采用濾波方法，顯然波完全被反射回細網(wǎng)格區(qū)域，而采用濾波方法后，這個反射波基本被消除了.因此，圖5展示不同條件下的波包傳播過程，證實AFEM加上低通濾波方法保證了K＝0.1波能傳播到組單元區(qū)域，而過濾掉了K＝1.0的波包傳播.

通過定義在細網(wǎng)格區(qū)域和粗網(wǎng)格區(qū)域間能量占有率，我們計算在某個波矢K下對應的反射系數(shù)和透射系數(shù)，結(jié)果如圖6所示.從反射系數(shù)可以看出，如果不考慮低通濾波方法，由之前分析可知全反射發(fā)生位置大約在π／8處.所以很顯然，發(fā)生全反射的波矢寬度遠大于需要透射的波矢寬度.在AFEM下加入濾波耗散方法后，大部分的虛假反射波消除了.同時從透射系數(shù)可以看到，對低頻傳播影響很小.所以對整個Brillouin區(qū)計算得到的反射系數(shù)和透射系數(shù)，進一步證實我們設計的耗散方法具有低通濾波性質(zhì).雖然也看到在反射系數(shù)中K＝0.5的位置有稍強的反射，原因主要是選擇的耗散參數(shù)β沒有很好地被優(yōu)化，這將是后續(xù)的工作.

圖6 三個不同模擬條件的反射系數(shù)和透射系數(shù)Fig.6 Reflection and transmission coefficients with three simulation conditions

5 結(jié)論

我們發(fā)展了一個基于原子尺度的有限元動力學方程.詳細的理論推導給出一維情形下AFEM的質(zhì)量矩陣和剛度矩陣具體形式，展示了其跨越從非局域的原子模型到局域的連續(xù)模型.比較其他物理模型，色散關系的計算證實AFEM擁有更多的振動頻率，散射計算結(jié)果表明AFEM也具有更高的截止全反射波矢.但是由于單元增加造成的無法避免虛假波反射，引入了一種不同以往的耗散方法.這個耗散方法基于濾波器的原理，目的是為了消除高頻的反射波，而不影響低頻的透射波，從而保證粗單元區(qū)域具有低頻的動力學行為.數(shù)值模擬不同波矢的傳播特性，以及反射和透射系數(shù)的計算證實了低頻濾波耗散方法有效地消除了大多數(shù)的反射波，而對低頻波影響很小，同時在低頻區(qū)域AFEM基本沒有影響低頻波的衰減.

原子尺度下的聲子是應力波及熱量傳播的載體，分別對應于低頻和高頻振動模.宏觀上則體現(xiàn)于連續(xù)力學本構(gòu)關系和熱傳導方程.耦合原子和連續(xù)區(qū)域的模型后，為保證原子區(qū)域的物理信息能有效地傳播到連續(xù)區(qū)域，原子區(qū)域的所有振動模需被分離為低頻和高頻.低頻的振動模需要保證傳播到連續(xù)區(qū)域體現(xiàn)力學特性而高頻則需要被收集作為能量體現(xiàn)熱傳播.建立多尺度模型后，由于大網(wǎng)格的存在導致高頻振動模被反射回原子區(qū)域，造成該區(qū)域局域溫升，體系處于非熱力學平衡.本工作的目的是消除高頻振動模而不影響低頻振動傳播，這個思想在近來發(fā)展的熱力學耦合多尺度方法可以看到［22］，本文的低通濾波耗散邊界更容易實現(xiàn).因此，我們發(fā)展的這個多尺度算法為發(fā)展熱力學耦合的多尺度方法提供了一條新的途徑.

［1］Miller R E，Tadmor E B.A unified framework and performance benchmark of fourteen multiscale atomistic／continuum coupling methods［J］.Modell Simul Mater Sci Eng，2009，17：1－51.

［2］Kohlhoff S，Gumbsch P，F(xiàn)ischmeister H F.Crack propagation in bcc crystals studied with a combined finite-element and atomistic model［J］.Phil Mag A，1991，64：851－878.

［3］Tadmor E B，Ortiz M，Phillips R.Quasicontinuum analysis of defecs in solids［J］.Phil Mag A，1996，73：1529－1563.

［4］Broughton J Q，Abraham F F，Bernstein N，Kaxiras E.Concurrent coupling of lengh scales：Methodology and application ［J］.Phys Rev B，1999，60：2391－1403.

［5］Rudd R E，Broughton J Q.Coarse-grained molecular dynamics：Nonlinear finite elements and finite temperature［J］.Phys Rev B，2005，72：144104.

［6］Wanger G J，Liu W K.Coupling of atomistic and continuum simulations using a bridging scale decomposition［J］.J Comp Phys，2003，190：249－274.

［7］Miller R，Tadmor E B，Phillips R，Ortiz M.Quasicontinuum simulation of fracture at the atomic scale［J］.Modelling Simul Mater Sci Eng，1998，6：607－638.

［8］Park H S，Liu W K.An introduction and tutorial on multiple-scale analysis in solids［J］.Comput Methods Appl Mech Eng，2004，193：1732－1772.

［9］Warner D H，Curtin W A，Qu S.Rate dependence of crack-tip processes predicts twinning trends in fcc metals［J］.Nature Mater，2007，6：876－881.

［10］Peng Q，Zhang X，Hung L，Carter E A，Lu G.Quantum simulation of materials at micro scales and beyond［J］.Phys Rev B，2008，78：054118.

［11］徐云，陳軍，陳棟泉.加載條件微裂紋動力學行為的多尺度模擬［J］.北京理工大學學報，2010，30：357－360.

［12］邵宇飛，王紹青.基于準連續(xù)介質(zhì)方法模擬納米晶體Ni中裂紋的擴展［J］.物理學報，2010，59：7258－7265.

［13］Bazant Z P.Spurious reflection of elastic waves in nonuniform finite element grids［J］.Comput Meth In Appl Mech Eng，1978，16：91－100.

［14］Zhou S J，Lomdahl P S，Thomson R，Holian B L.Dynamic crack processes via molecular dynamics［J］.Phys Rev Lett，1996，76：2318－2321.

［15］Qu S，Shastry V，Curtin W A，Miller R E.A finite-temperature dynamic coupled atomistic／discrete dislocation method［J］. Modell Simul Mater Sci Eng，2005，13：1101－1118.

［16］Doll J D，Myers L E，Adelman S A.Generalized Langevin equation approach for atom／solid-surface scattering：Inelastic studies［J］.J Chem Phys，1975，63：4908－4914.

［17］Cai W，Koning M，Bulatov V V，Yip S.Minimizing boundary reflections in coupled-domain simulations［J］.Phys Rev Lett，2000，85：3213－3216.

［18］E W，Huang Z.Matching conditions in atomistic-continuum modeling of materials［J］.Phys Rev Lett，2001，87：135501.

［19］Yang J Z，Li X.Comparative study of boundary conditions for molecular dynamic simulations of solids at low temperature［J］. Phys Rev B，2006，73：224111.

［20］To A C，Li S.Perfectly matched multiscale simulations［J］.Phys Rev B，2005，72：035414－1－8.

［21］Jones R E，Kimmer C J.Efficient non-reflecting boundary condition constructed via optimization of damped layers［J］.Phys Rev B，2010，81：094301.

［22］Mathew N，Picu R C，Bloomfild M.Concurrent coupling of atomistic and continuum models at finite temperature［J］.Comput Methods Appl Mech Engrg，2011，200：765－773.

Low-pass Filter Boundary Condition for Atomic-based Finite Element Method

REN Guowu，TANG Tiegang，LI Qingzhong
（Institute of Fluid Physics，China Academy of Engineering Physics，Mianyang，Sichuan 621999，China）

An atomic-based finite element method is developed，spanning from atomic scale to macroscopic one.With theoretical deduction dispersion relation and dynamic scattering behavior are calculated.For spurious reflection caused by non-uniform grid，lowpass filter boundary condition is designed to effectively eliminate the reflection of high-frequency phonons，while keeping low-frequency ones transparent.These schemes are demonstrated with numerically calculating reflection and transmission coefficients in onedimensional modeling.

finite element；filter；reflection coefficient

date：2013－09－12；Revised date：2013－12－19

O242

A

2013－09－12；

2013－12－19

國家自然科學基金（11102191）及流體物理研究所發(fā)展基金（SFZ20120402）資助項目

任國武（1981－），男，助理研究員，主要從事材料力學計算研究，E-mail：guowu.ren＠gmail.com

1001-246X（2014）05-0514-09