近場動力學系統中波傳播特性的探究1)

朱競高 任曉丹

(同濟大學土木工程學院,上海 200092)

引言

在載荷的作用下,固體的變形并不是瞬間完成的[1],而是波傳播、演化的結果.從這種意義上來說,波的傳播對固體的行為有著極其重要的影響,是固體力學中的基礎問題[2].特別是對于斷裂問題,波的傳遞更是影響著裂紋萌生、成核、分岔乃至破碎的全過程[3-5],是正確捕捉裂紋擴展的關鍵.

隨著計算機技術的發展和固體力學研究的深入,數值計算方法紛紛涌現[6-8],為研究固體力學中的科學問題提供了更為便捷的途徑.其中,近場動力學作為一類基于非局部思想的新固體力學方法[9],從鄰域內粒子間的相互作用出發,考慮了微納米尺度下的長程力,在模擬裂紋方面更具合理性.該方法采用積分形式的控制方程,利用微觀鍵的失效來描述斷裂,能夠從根本上避免求解不連續問題的奇異性和復雜性,無需特殊處理即可實現對裂紋發展的捕捉[10-12],被廣泛用于斷裂問題的求解.經過發展,目前形成了鍵基[13-14]和態基[15-17]近場動力學兩大分支.Zhou 等[18]和Zhu 等[19]還引入了剪切效應,以更準確地描述材料的力學性能.隨著對材料非線性研究的深入,Gu 等[20]和Liu 等[21]建立了彈塑性模型,王涵等則發展了熱黏塑性模型[22-23],在預測各類材料的斷裂失效問題中取得重要的進展[24-27].

然而,近場動力學作為一種非局部系統,不可避免地會引入波色散問題.波色散會改變波的傳播特性,不利于對固體行為的正確捕捉,進而影響對裂紋擴展的準確預測.Silling[13]在提出鍵基模型時,首先發現了近場動力學中的非線性色散行為.2016 年,隨著Baz?nt 等[28]從非局部性的角度對色散關系開展了進一步研究后,近場動力學中的色散問題逐漸引起了學者們的關注.Butt 等[29]對態基近場動力學中的色散關系進行了研究,提出調整鄰域半徑尺寸和影響函數來減弱波色散的方法.Wildman[30]分析了鍵基近場動力學中的色散關系,并實現了通過改變微彈模參數來匹配所需的色散方程.Gu 等[31]和Li 等[32]也開展了相關研究,探討不同衰減函數對近場動力學中色散關系的影響.除了通過調整參數來減小波色散外,還有學者[33-34]從借鑒有限差分法中色散研究的角度出發,提出有限差分增強的近場動力學方法.但是,這些研究主要局限于低頻成分,忽視了高頻成分的色散行為,同時也缺乏對波間斷性的考慮.

近年來,恐怖襲擊和局部戰爭頻發,極端載荷作用下工程結構的安全性問題引起廣泛關注,成為國防安全領域的研究重點.極端載荷,如爆炸、沖擊作用,往往伴隨著激波的產生[35-37],激波的傳播又影響著材料破碎和裂紋發展,是進行工程結構安全設計的關鍵.而激波因其高頻、強間斷的復雜性,給相關分析理論和計算方法帶來了挑戰.近場動學方法采用積分形式的控制方程,自然地適用于激波引起的材料破碎和裂紋發展模擬,但該方法對激波的準確捕捉卻不可避免地會受到高頻、強間斷成分色散特性的制約.

因此,對近場動力學的波色散特性展開全面系統的研究,納入對高頻成分和波間斷性的考慮,構成了文章的出發點.理論研究表明,相比于低頻成分,高頻成分還表現出振蕩關系和零能模式,色散問題更為突出.高頻域的色散關系會隨波的傳播方向發生變化,引起空間傳播的方向性.數值計算則發現強間斷性波的色散現象更為顯著,呈現出Gibbs 不穩定性.說明色散不僅會對波的傳播特性產生不利影響,還會引起數值不穩定性.為了解決上述問題,通過黏性引入耗散來抑制色散行為,并考慮對高頻成分的選擇性抑制,構造了相應的黏性力態.本文研究工作將為在近場動力學的框架下準確模擬波的傳播過程,正確捕捉固體行為提供重要參考.

1 近場動力學方法

1.1 基本思想

近場動力學通過引入特征長度 δ,定義了鄰域范圍

如圖1 所示,在鄰域范圍內,粒子間由鍵相連,鍵隨著粒子運動而發生變形.對某一空間域 ?,假設在某一時刻t,域內任一物質點x只與其鄰域范圍內的其他物質點x′之間存在相互作用力,則根據牛頓第二定律可得近場動力學的控制方程為[38]

圖1 近場動力學模型示意圖Fig.1 The schematic of peridynamics

其中,ρ為物質密度;u為物質點的位移;b為體力密度;f為力密度函數,包含了材料的本構關系.

根據本構關系不同,可將近場動力學分為鍵基近場動力學和態基近場動力學.鍵基近場動力學由于模型簡單、計算方便,在斷裂問題中的應用更為廣泛.因此,本文的研究將基于鍵基模型進行展開.這里,對鍵基近場動力學作簡要介紹.

1.2 鍵基近場動力學

相比于態基近場動力學,鍵基近場動力學的模型較為簡化.該模型認為鍵與鍵之間相互獨立,力密度函數只與單一鍵有關[13].若定義粒子間的相對位置和相對位移為

則未變形的鍵和變形后的鍵可以分別用 ξ和 η +ξ表示.由于采用了鍵間獨立性的假定,粒子間的相互作用可視為中心彈簧,將微觀勢函數表示為

其中,s為鍵的相對伸長量,C (∥ξ∥)為微彈性模量.對勢函數求導可得力密度函數的表達式

其中n為變形后鍵的方向

微彈性模量有多種函數形式,當取為常數c時,鍵基近場動力學就變成了微觀彈脆性模型,將被用于后面的理論分析和數值計算.該模型根據應變能等效,可得不同維度下的微彈性模量

其中,E為彈性模量,A為一維桿的橫截面積,B為二維板的厚度,二維情況只選取了平面應力問題作為代表.此時,力密度函數可寫為

2 波的傳播特性分析

這里,將從微觀彈脆性模型入手,兼顧低頻和高頻成分,對近場動力學系統中波的傳播特性展開全面的分析.采用小變形假設,微觀彈脆性模型的控制方程可以寫為[28]

2.1 一維譜分析

對于一維問題,式(9)可以寫為

其中,A為一維模型的橫截面積.為了研究波的傳播特性,有必要在頻域內展開分析.這里關注線彈性問題,因此,可對位移項進行傅里葉展開

其中,i 為虛數單位,滿足 i2=?1 ;a為幅值;v為相速度;κ為 2 π長度內的波數.將式(11)代入式(10),可約去時間項得

進一步,考慮鄰域的對稱性,可約去奇數項得

若定義與波數有關的項為

則式(13)可視為相速度v的方程

求解可得相速度為

顯然,式(16)中的虛數項為0,故近場動力學系統中沒有數值耗散.當波色散產生不利影響時,如數值振蕩、非物理變形等,系統本身無法對其進行抑制.隨著不利因素的逐漸累積,數值計算的精度將會受到嚴重的影響,無法保證對固體行為的正確捕捉.這一方面說明了進行近場動力學系統中波色散研究的必要性,另一方面為抑制色散現象提供了方向.

進一步,考慮數值計算中的空間離散,對近場動力學系統的色散特性進行分析.根據鄰域的對稱性,經過數值離散,式(14)變為

其中,m為一維鄰域內的粒子數,h為粒子間距,c為微彈性模量.采用式(7)中的微彈性模量,式(17)可寫為

將式(18)代入式(16),可得一維近場動力學系統的色散關系

其中,v0為彈性理論波速

其中,E為彈性模量,ρ為密度.由于 κ δ和 κ ξi均與 κh有關,式(19)可視為Nh=κh/(2π)的函數,表示粒子間距內的波數.若取鄰域半徑為 δ=3.015h,式(19)中的色散關系將如圖2 所示.

由圖2 可知,在一維近場動力學系統中,相速度與理論波速的比值偏離局部理論的結果,表現出色散行為.局部理論各個頻率成分的波的傳播速度相同,在傳播時會同時達到.而近場動力學中各個頻率成分的波的傳播速度不同,在傳播時會分先后達到.對于傳播速度小于理論波速的頻率成分,在傳播中會落后于其他波,表現為拖尾波.因此,定義臨界點為與局部理論色散關系的交點,即的點.可以看出,大部分頻率成分位于臨界點以下,在波的傳播中表現為拖尾波.

圖2 一維近場動力學系統的色散關系Fig.2 The dispersion relation of 1D peridynamics

為區分色散行為和研究方便,將Nh≥1定義為高頻域,0 ≤Nh<1定義為低頻域.即高頻域與低頻域的劃分與波長有關,高頻域的波長范圍為 λ ≤h,低頻域的波長范圍為 λ >h,其中h為粒子間距.顯然,低頻成分(虛線段)和高頻成分(實線段)分別呈現出不同的色散特性.在高頻域,波速不僅存在振蕩的趨勢,還在某些頻率處等于0.根據文獻[28-29],這些點的相速度為0,對應于零能模式點.另有學者將零能模式定義為群速度為0 的頻散帶,在這種定義下,首先有必要給出一維模型的頻率表達式.基于式(19),有

顯然,由于余弦函數具有周期性,頻率將呈現在Nh取整區間內的周期性變化.為了更加直觀地說明此周期性,將頻率與Nh的關系繪制成圖3.

圖3 一維近場動力學系統的頻率關系Fig.3 The frequency relation of 1D peridynamics

如圖3 中藍色線所示,頻率隨波數呈現周期性變化,并在波數取整處為0.與圖2 不同的是,這種周期性在低頻域(虛線)和高頻域(實線)呈現相同的變化規律.由于頻率呈現從0 到0 的周期性變化,對應的斜率應當呈現從正到負或從負到正的周期性變化.而頻率的斜率正比于群速度,說明存在群速度為0 的點.為了進一步驗證上述結論,基于頻率的表達式(21),求導得到群速度與波數的關系,繪制圖像進行表示.如圖3 中紅色線所示,群速度會在整波數區間內從正變為負,并在從正變為負的過程中出現零點.在這種情況下,零能模式點對應于圖3 中紅色線與y=0 水平線的交點.

綜上所述,無論是采用何種形式的零能模式定義,一維近場動力學模型中均存在零能模式點.振蕩趨勢和零能模式均是非物理響應,會加重色散行為帶來的不利影響.因此,相比于低頻成分,高頻成分的色散問題更為嚴重,應當給予關注.

2.2 二維譜分析

下面,考慮波的空間傳播,對二維系統展開分析.此時,除方向x1外,還需引入方向x2,控制方程(9)變為

其中,B為二維模型的厚度,ξ為相對位置

同時,位移項的傅里葉展開變為

其中,a為與位移u共線的幅度向量;b為波傳播方向的單位矢量;其他變量的定義同式(11).將式(24)代入式(22),可同時約去矢量項和時間項,有

同樣地,根據鄰域的對稱性,可約去奇數項得

仍可將與波數有關的項用一個變量進行表示

此時,關于波速的方程可寫為

解得相速度為

與一維的譜分析結果相同,相速度的虛數項為0,沒有數值耗散,無法抑制色散帶來的不利影響;實數項與波數有關,表現出色散行為.

進一步,考慮數值計算,對Nκ進行空間離散,有

其中,n為二維鄰域內的粒子數,h為粒子間距,ξj·b為相對位置在波傳播方向上的投影.仍采用式(7)中的微彈性模量,式(30)可寫為

將式(31)代入式(29),可以得到二維近場動力學系統的色散關系

其中,v0為式(20)定義的彈性理論波速.顯然,在二維問題中,波的色散行為不僅與頻率有關,還與波的傳播方向有關.

如圖4 所示,對于對稱鄰域,對關于x2軸對稱的傳播方向b′,一定可以同樣對稱地找到關于x2軸對稱的相對位置向量,使式(30)中的投影相等

圖4 二維鍵基模型色散關系的對稱性Fig.4 The symmetry of dispersion relation for 2D peridynamics

由于鄰域內每個相對位置向量 ξj,都存在關于x2軸對稱的相對位置向量滿足投影相等關系,故最終求和得到的式(30)也相等.同樣地,對于x1軸也滿足此對稱性

因此,可將波的傳播方向縮小至 θ=0 ～90?范圍內進行研究.顯然,在 0?～90?范圍內,關于45?方向存在對稱性

因此,二維鍵基模型的色散行為將表現出沿45?的對稱性.

基于此對稱性,對波傳播方向與x1方向的夾角在 θ=0?～45?之間的色散行為進行研究.這里,仍取Nh=κh/(2π)為橫坐標,并取波的傳播方向分別為θ=0?,30?,45?,式(32)中的色散關系將如圖5 所示.由圖5 可知,當波沿x1方向傳播時,二維近場動力學系統中色散關系的變化趨勢與一維系統較為一致,即在低頻成分和高頻成分表現出不同的色散行為,并在高頻域呈現出振蕩趨勢和零能模式.但隨著傳播方向改變,高頻成分色散關系的振蕩程度發生變化,零能模式點的位置也在移動.

圖5 二維近場動力學系統的色散關系Fig.5 The dispersion relation of 2D peridynamics

為了更加直觀地表征不同方向的頻散行為,提取高頻成分色散關系變化率的最大值作為振蕩程度的衡量指標,將 0?～360?范圍內高頻成分色散關系的振蕩程度繪制雷達圖進行對比分析.如圖6 所示,對于對稱鄰域,高頻成分色散關系的變化趨勢呈現出關于 4 5?方向的對稱性,驗證了前述分析的正確性.在 0?～45?方向范圍內,隨著波的傳播方向改變,振蕩指標的大小發生變化,與圖5 的分析結果相一致.

圖6 不同方向頻散行為的雷達圖Fig.6 Radar chart of dispersion behavior in different directions

綜合上述分析,相比于低頻成分,高頻成分的色散行為更為嚴重,表現出振蕩趨勢和零能模式.高頻域的色散關系會隨波的傳播方向變化,呈現出沿45°的對稱性.而近場動力學系統本身沒有數值耗散,無法抑制色散帶來的不利影響,會嚴重影響對固體行為的正確捕捉.

3 近場動力學方法的改進

事實上,作為一種非局部理論,近場動力學在高頻段呈現頻散行為是非局部效應的合理體現.對于非均質材料,頻散行為是真實存在的,可以通過近場動力學來模擬非均質材料中的非線性色散關系,且近場作用半徑與非均質材料的特征尺寸具有強相關性[29].但是,對于變異性不大的材料,通常可以認為是均質材料,近場動力學的頻散行為就會帶來不利影響.根據第2 章的研究,近場動力學的頻散行為在高頻域表現出振蕩趨勢.也就是說,理論上應該同時達到的不同頻率成分的波,在近場動力學的計算中傳播速度會出現差異,這種差異會在數值模擬中引起非物理變形和損傷.

對于極端載荷作用下的問題,一方面這種不利影響會加重,嚴重影響數值模擬的精度.另一方面,極端載荷作用往往伴隨著裂紋的產生,而裂紋的擴展與波的傳播緊密相關,異常的波傳播過程必定會對裂紋擴展的預測產生不利影響.而極端載荷作用下工程結構的安全性問題是國防安全領域的研究重點,近場動力學方法的優勢正是對于裂紋擴展捕捉的能力.因此,有必要抑制近場動力學中的頻散行為,提高其對波傳播的模擬精度,進而保證其對極端載荷作用下計算的可靠性,為國防安全鄰域研究提供理論指導和技術支撐.本章將從引入數值耗散的角度出發,對近場動力學體系的色散行為進行抑制.并考慮高頻域色散更為嚴重的特點,對高頻成分進行選擇性抑制.

3.1 黏性效應的研究

在固體力學中,通常認為黏性與數值耗散有關,為近場動力學中耗散機制的引入提供了方向.首先,有必要對黏性效應進行研究,證明黏性引入的有效性.這里,仍從微觀彈脆性模型入手.考慮到黏性與應變率的相關性,可將黏性項取為力密度函數的率形式

其中,γ為無量綱系數.在小變形下,有

以一維問題為例,式(38)可化簡為

其中,A為一維模型的橫截面積.沿用第2 章的譜分析方法,采用式(11)中的傅里葉展開,式(39)可寫為

其中,v′為黏性作用下的相速度.仍采用式(14)中對波數相關項的定義,式(40)可表示為相速度的方程

求解可得相速度

此時,虛數項不再為0

說明黏性項能夠引入數值耗散,這與文獻[40-43]的研究結論相一致,證明了黏性引入的有效性.

進一步,考慮高頻域色散更為嚴重的特點,對黏性項在不同頻率成分的表現進行研究.由于數值耗散的引入是為了抑制色散行為帶來的不利影響,這里關注虛數項(43)與原始相速度(16)的比值.借助式(18),比值可以表示為

其中,v0為式(20)定義的彈性理論波速,γ為無量綱系數,δ為鄰域半徑,h為粒子間距.仍取Nh=κh/(2π)為橫坐標,歸一化后式(44)中的比值關系將如圖7所示.由圖7 可知,數值耗散對色散行為的抑制作用在低頻域(虛線段)和高頻域(實線段)呈現相同的周期變化,即黏性效應對低頻和高頻成分的表現相同.實際上,式(37)本質上是一種線性黏性,無法實現對不同頻率成分的選擇性抑制.

圖7 數值耗散對不同頻率成分的抑制作用Fig.7 The dissipation effect on different frequencies

3.2 黏性效應的引入

因此,可將黏性效應引入近場動力學,抑制波色散帶來的不利影響.并有必要在黏性項的構造中,考慮對高頻成分的選擇性抑制.基于近場動力學的非局部性和鍵基模型的鍵間獨立性假設,將在鍵層次以非局部作用的方式引入黏性效應.相應地,將控制方程定義為

其中,f為力密度函數,是傳統力態;fvs為黏性力態.同時,考慮到在鍵基模型中,力密度函數的方向為變形后鍵的方向,故滿足線動量守恒

和角動量守恒

為了保持上述守恒關系,取黏性力態的方向與力密度函數的方向相同.此時,黏性力態可以分解為標量項和方向項

其中,fvs為黏性力態的標量形式,n為變形后鍵的方向.顯然,這種引入方式不僅保留了傳統近場動力學模型的理論框架,而且只需要給出黏性力態標量形式的定義.

下面,將針對體積變形,這一常見的固體變形模式,構造黏性力態的標量表達式.參照連續介質力學中體積黏性的定義[44-45]

其中,β為無量綱體積黏性系數;cd為膨脹波速;Le為單元的特征尺寸;為體積應變率

其中,u為位移場.定義體積黏性力態的標量形式為

其中,b1為無量綱體積黏性系數;v0為式(20)定義的彈性理論波速;VH為鄰域的體積;ζ為應變率.在鍵基模型中,鍵的相對伸長量s為一種應變測量,可用于應變率的計算

如圖8 所示,∥ ξ∥ 為鍵的初始長度(藍線),η ˙·n為鍵的速度在變形后鍵的方向上的投影(紅線).

圖8 應變率的運動學表示Fig.8 The kinematic expression of strain rate

基于上述定義,可得體積黏性力態的標量表達式

當b1αp=γc時,式(54)與式(39)相同,證明了3.1 中黏性效應分析的正確性.

由3.1 中黏性效應分析可知,式(53)本質上是一種線性黏性,無法實現對不同頻率成分的選擇性抑制.考慮到高頻成分波色散的嚴重性,有必要對其進行非線性改進.這里,將引入二次項,以放大對高頻成分的抑制作用

其中,b2為無量綱二次黏性系數,其余變量的定義與式(51)一致.為充分發揮一次項和二次項的黏性效應,將式(53)和式(55)進行多項式組合

其中,b1和b2分別為一次項和二次項的黏性系數,得到的黏性力態稱為多項式黏性.

式(48)給出了引入黏性的理論框架,式(56)則從體積變形出發,考慮了對高頻成分的選擇性抑制,構造了多項式黏性力態,完成了對近場動力學方法的改進.

4 數值算例

本章基于波傳播的數值算例,開展數值模擬研究.旨在對理論分析進行驗證的同時,考察波間斷性的影響.在近場動力學中,通常采用EMU 方法[46]進行數值計算,該方法對控制方程進行空間離散

其中,N為鄰域內的粒子數目,?V為離散粒子的體積.由于波的傳播為動力問題,采用顯式中心差分方法進行求解

其中,u,和分別為加速度,速度和位移;?t為滿足穩定條件的時間步長.

4.1 一維波傳播問題

對一維算例,首先考慮帶理論解的波傳播問題,進行對比研究,以驗證理論分析的正確性.再在此基礎上,考慮極端載荷下的激波傳播問題,探究波間斷性的影響.

4.1.1 帶理論解的波傳播問題

為了驗證理論分析的正確性,并說明對不同間斷性波的適用性,這里取具有解析表達式的光滑正弦波和間斷方波為研究對象,對一維波的傳播問題進行了研究.將1 m 長的彈性桿一端固定,另一端輸入目標波,提取t=77.3 μs 時沿桿長的位移分布,將鍵基模型、線性黏性和多項式黏性的計算結果與理論解進行對比分析.圖9 和圖10 分別給出的是正弦波和方波的計算結果.顯然,兩種波的鍵基模型計算結果均在波后出現了數值振蕩.結合理論分析結果,一維近場動力學系統中存在波色散,且大部分頻率成分的波速小于理論波速,表現為拖尾波的形式,第二章中譜分析的正確性得到了驗證.

圖9 帶理論解的一維正弦波傳播結果對比Fig.9 Comparative study of 1D sine wave propagation problem with theoretical result

圖10 帶理論解的一維方波傳播結果對比Fig.10 Comparative study of 1D square wave propagation problem with theoretical result

提取波后振蕩進行進一步研究,發現在兩種波形的計算結果中,黏性的引入均可以有效地抑制數值振蕩.且多項式黏性能夠在抑制波后振蕩的同時,更好地保持原始波峰和波形,得到了更為合理的波傳播結果.這驗證了第3 章所提方法的正確性.需要注意的是,相比于正弦波,方波的波后振蕩更為顯著,說明近場動力學中的色散行為與波的間斷性有關,這將在后續算例中進行研究和討論.

4.1.2 激波傳播問題

考慮1 m 長彈性桿,在兩端受爆炸載荷作用的波傳播問題.如圖11 所示,兩端炸藥距桿端均為0.5 m,且同時發生爆炸.為探究波間斷性的影響,設置小當量炸藥(W=1 kg)和大當量炸藥(W=80 kg)兩類工況進行計算.桿的材料參數為 ρ=E=193 GPa.將一維桿離散成1000 個粒子進行計算,鄰域半徑取 δ=3.015h,采用壓力時程模擬爆炸載荷,并忽略負壓區的影響.

圖11 一維波傳播算例示意圖 (單位: mm)Fig.11 The schematic of numerical example for 1D wave propagation (unit: mm)

(1)小當量工況(W=1 kg)

圖12 和圖13 分別為5 個典型時刻(t=50,100,150,200,250 μs)沿桿長方向的位移和應力分布.根據桿長和材料特性,波的傳播周期為

計算結果表明,爆炸載荷產生壓力波,壓力波由桿的兩端向中間傳播.t=100 μs 時在中點(A點)相遇后,繼續向前進行傳播.t=200 μs 時達到桿端后,反射產生拉伸波,向相反方向進行傳播.模擬得到的波的傳播周期與式(59)相符,驗證了數值計算的正確性.相比于圖12,圖13 中的應力波并不光滑,在波后出現了數值振蕩.結合第2 章的譜分析結果,一維近場動力學系統中存在波色散,且大部分頻率成分的波速小于理論波速,表現為拖尾波的形式,譜分析的正確性得到了驗證.

圖12 沿一維桿長方向的位移分布Fig.12 Displacement distribution along the 1D bar

圖13 沿一維桿長方向的應力分布Fig.13 Stress distribution along the 1D bar

進一步,提取A點的應力波時程,對黏性效應進行研究.如圖14 所示,兩端的壓力波在t=100 μs 時相遇疊加,反射產生的拉伸波在t=300 μs 時再次相遇疊加,波的傳播周期與理論解相符.對波峰處放大進行觀察,鍵基系統的計算結果在波后出現了數值振蕩,這與圖13 的計算結果相一致.施加黏性后,數值振蕩得到了有效地抑制,證明了黏性引入的有效性.

對比兩種黏性力態,多項式黏性能夠在抑制波后振蕩的同時,更好地保持原始波峰和波形,得到更為合理的波傳播結果.這是因為數值耗散的引入會消耗波的能量,從而影響波峰和波形.而多項式黏性中含有二次項,可以放大對高頻成分的抑制作用,在同等的抑制效果下減小了能量耗散,降低了對原始波峰和波形的影響.為了更直觀地從能量角度進行說明,對圖14 進行頻域分析,得到應力波的能量譜密度曲線.如圖15 所示,相比于體積黏性,多項式黏性能夠降低能量耗散,且二者的差異主要集中在105～106的頻率區間.

圖14 A 點的應力波時程Fig.14 Time history of stress wave at point A

圖15 A 點應力波的能量譜密度曲線Fig.15 Energy spectrum density of stress wave at point A

(2)大當量工況(W=80 kg)

下面,納入對波間斷性的考慮,對大當量工況進行計算.此時,t=50,100,150,200,250 μs 時刻沿桿長方向的位移和應力分布分別如圖16 和圖17 所示.顯然,相比于小當量工況,此時的位移和應力的波形變得更為尖銳,即波的間斷性變強,對應于激波的產生.與小當量的計算結果相比,位移也出現了波后振蕩,應力的數值振蕩也更加明顯.說明強間斷性波的色散行為更為顯著,呈現出Gibbs 不穩定性.因此,色散在對波傳播特性產生不利影響的同時,還將引起數值不穩定性.

圖16 沿一維桿長方向的位移分布Fig.16 Displacement distribution along the 1D bar

圖17 沿一維桿長方向的應力分布Fig.17 Stress distribution along the 1D bar

仍提取點A的應力波時程進行進一步分析,應力波時程和能量譜密度曲線分別如圖18 和圖19 所示.與圖17 的計算結果相一致,圖18 中波的間斷性變強,波后的振蕩加劇.與小當量工況相同,黏性能夠有效地抑制色散造成的數值振蕩,且多項式黏性能夠更好地保持原始波峰和波形,說明所提方法對激波同樣適用.圖19 則從能量角度說明了多項式黏性的優勢,即減小了能量耗散,從而降低了對波峰和波形的影響.

圖18 A 點的應力波時程Fig.18 Time history of stress wave at point A

圖19 A 點應力波的能量譜密度曲線Fig.19 Energy spectrum density of stress wave at point A

綜合一維計算結果,近場動力學系統中的波色散在一維波傳播問題中表現為波后振蕩,強間斷性波的色散行為更為顯著,呈現出Gibbs 不穩定性.黏性能夠有效地抑制波后振蕩,且多項式黏性能夠更好地保持原始波峰和波形.

4.2 二維波傳播問題

對二維算例,考慮內徑0.15 m,外徑0.3 m 的圓環在爆炸載荷作用下的波傳播問題.如圖20 所示,炸藥在圓環的中心位置處爆炸,距圓環內壁0.15 m.圓環的材料參數為,E=30 GPa.這里,取炸藥質量W=70 kg,對大當量爆炸下的激波傳播進行模擬.在數值計算中,粒子間距取h=10 mm,鄰域半徑取 δ=3.015h,爆炸載荷仍采用壓力時程進行模擬,并忽略負壓區的影響.

圖20 二維波傳播算例示意圖 (單位: mm)Fig.20 The schematic of numerical example for 2D wave propagation (unit: mm)

中心爆炸下圓環中波的傳播本質上是球形腔加壓問題,存在理論解[47].這里,關注壓力云圖,并采用對比研究的方式,以考察色散對波的空間傳播行為的影響.由圖21 (a1)～ 圖21(c1) 可知,爆炸理論上會產生壓力波,壓力波沿周向均勻分布,并逐漸沿徑向向外傳播.而在鍵基系統的計算結果中,如圖21 (a2)～圖21(c2)所示,壓力波的分布并不均勻,呈現出沿45°的對稱性.結合第2 章的譜分析結果,二維系統中高頻成分振蕩的色散關系具有沿45°方向的對稱性,數值計算結果驗證了譜分析的正確性.

圖21 二維圓環的壓力云圖Fig.21 The pressure cloud map of 2D ring

此外,相比于一維問題,二維情況下的拖尾波現象更為顯著,拖尾波的峰值(3 GPa)甚至超過了波前峰值(1.5 GPa).拖尾波的存在不僅會在波后引起非物理變形,還會消耗波前峰值,不利于對波傳播過程的正確捕捉.因此,在二維情況下,波色散會通過影響波的空間分布和傳播過程,對波的空間傳播行為產生不利影響.

進一步,觀察多項式黏性的計算結果,以探究黏性在二維問題中的表現.如圖21 (a3)～ 圖21(c3)所示,此時的拖尾波得到了有效的抑制,降低了對波前峰值的消耗,波的傳播過程得到了較好地捕捉.同時,壓力波的分布更為均勻,與理論解更為一致.因此,在黏性的作用下,近場動力學系統中波的空間傳播行為更為合理,證明了黏性引入的有效性.此時,壓力波的分布寬度相比理論解有所增加,這應當與特征長度和數值耗散有關.近場動力學引入了特征長度,即鄰域半徑 δ,在特征長度范圍內考慮粒子間的相互作用.特征長度的存在會在一定程度上擴大影響范圍,使波的分布寬度變大.另一方面,由一維數值模擬結果可知,數值耗散會消耗波的能量,在一定程度上增大波寬.但總體來說,黏性的引入,使得波前峰值和波的空間分布均得到了較好地捕捉,證明了所提方法同樣適用于多維問題.

5 結論

本文采用譜分析方法,對近場動力學系統中不同維度、不同頻率成分的色散特性進行了研究,并提出了相應的改進方法.同時進行數值計算,驗證理論分析結果,并探究波間斷性的影響,主要結論如下.

(1) 相比于低頻成分,高頻成分的色散關系還呈現出振蕩趨勢和零能模式,色散問題更為突出.

(2) 高頻成分的色散行為會隨波的傳播方向發生變化,呈現出沿45°的對稱性.

(3) 色散在影響波傳播特性的同時,還會引起Gibbs不穩定性,表現為強間斷性波的色散現象更為顯著.

(4) 黏性能夠引入數值耗散,抑制色散帶來的不利影響.二次項的引入可以放大對高頻成分的抑制作用,在同等的抑制效果下,減小對原始波峰和波形的影響.

(5) 色散在一維波傳播中引起波后振蕩,在二維問題中將引起更為嚴重的拖尾波和波空間分布的不均勻性.

本文研究可為減小近場動力學系統中波色散的不利影響,獲得正確的波傳播過程和固體行為提供重要參考.同時,尚存在一些值得進一步研究的問題,如在沖擊問題中的應用、對零能模式的抑制等.事實上,若想納入針對零能模式的抑制作用,需要先識別出近場動力學中的零能變形模式,再針對其進行抑制.零能模式是一種不引起外力的虛假變形模式,可以從零能模式的特點出發,對近場動力學中的零能變形模式進行研究,確定其變形特征.再在本文提出的數值黏性的基礎上,根據研究出的近場動力學中零能模式的變形特征,對變形相關項進行修改,以實現對零能模式的抑制效果.