近場動(dòng)力學(xué)算子方法1)

李志遠(yuǎn) 黃 丹 Timon Rabczuk

* (河海大學(xué)工程力學(xué)系,南京 211100)

? (魏瑪-包豪斯大學(xué)結(jié)構(gòu)力學(xué)研究所,德國魏瑪 99423)

引言

20 多年來,近場動(dòng)力學(xué)(peridynamics,PD)[1-3]在計(jì)算力學(xué)與相關(guān)工程領(lǐng)域受到了廣泛關(guān)注.它采用空間積分方程代替經(jīng)典連續(xù)介質(zhì)力學(xué)中的空間微分方程,從而避免了基于連續(xù)性假設(shè)的傳統(tǒng)局部模型在面臨不連續(xù)問題時(shí)的奇異性,在處理諸多工程領(lǐng)域?qū)嶋H問題時(shí)表現(xiàn)出較明顯優(yōu)勢[4-5],如沖擊破壞[6-7]、水力劈裂[8-9]和熱力耦合[10-12]等.

PD 模型通常可分為3 類: 鍵型近場動(dòng)力學(xué)(bond-based peridynamics,BB-PD)、常規(guī)態(tài)型近場動(dòng)力學(xué)(ordinary state-based peridynamics,OSB-PD)和非常規(guī)態(tài)型近場動(dòng)力學(xué)(non-ordinary state-based peridynamics,NOSB-PD).BB-PD[1]中物質(zhì)點(diǎn)間成對的相互作用類似于獨(dú)立的彈簧,在應(yīng)用中存在泊松比限制等缺陷.態(tài)型PD[2]模型有效地克服了這一限制.當(dāng)限定泊松比時(shí),OSB-PD 可簡化為BB-PD.NOSB-PD 提出關(guān)聯(lián)模型的概念,可重構(gòu)經(jīng)典連續(xù)介質(zhì)力學(xué)模型,但它受零能模式引起的數(shù)值振蕩影響[13].為消除振蕩,學(xué)者們又提出一些處理措施,如附加額外的力狀態(tài)[13-14]和鍵關(guān)聯(lián)的近場范圍[15]等.對偶域近場動(dòng)力學(xué)[16](dual-horizon Peridynamics,DHPD)的提出則突破了傳統(tǒng)PD 模擬中的固定近場范圍限制.

基于PD 的非局部作用思想,2016 年Madenci等[17]提出近場動(dòng)力學(xué)微分算子(peridynamic differential operator,PDDO).PDDO 求解偏微分方程(partial differential equations,PDE)是基于強(qiáng)形式,并通過拉格朗日乘子法施加邊界條件.2020 年Ren 等[18-19]提出非局部算子方法(nonlocal operator method,NOM).NOM 求解PDE 是基于弱形式,并運(yùn)用罰能量泛函來抑制零能模式.運(yùn)用PDDO 和NOM 均可以將PDE 從局部微分形式重構(gòu)為非局部積分形式.這兩類非局部算子近年來均在各種物理學(xué)方程求解中得到應(yīng)用.Li 等[20]提出用于求解PDE 的鍵關(guān)聯(lián)的PDDO 弱形式,并通過引入鍵關(guān)聯(lián)的近場范圍[15]來消除零能模式.周保良等[21-23]基于PDDO 研究了瞬態(tài)熱傳導(dǎo)[21],正交各向異性熱傳導(dǎo)[22]以及非線性熱傳導(dǎo)[23].Li 等[24-27]應(yīng)用PDDO 分析了復(fù)合材料與結(jié)構(gòu)的熱彈性[24-25]、動(dòng)力特性 [26]以及大變形[27].Ren 等[28-30]運(yùn)用NOM 求解高階梯度固體[28]、斷裂相場[29]與電磁波導(dǎo)[30]等問題.

在PDDO 和NOM 這兩類非局部算子的已有研究基礎(chǔ)上,本文進(jìn)一步提出一種更為一般的近場動(dòng)力學(xué)算子方法(peridynamic operator method,PDOM).PDOM 不僅可直接將局部微分轉(zhuǎn)化為相應(yīng)的非局部積分形式,而且同樣適用于微分的乘積.因此,PDDO 和NOM 皆可視為PDOM 的一種特例.以彈性力學(xué)問題為例,下文基于變分原理和拉格朗日方程,推導(dǎo)了適用于靜/動(dòng)態(tài)彈性力學(xué)問題的PDOM 模型,并證明了,當(dāng)分別限定相互作用域?yàn)榕c位置無關(guān)或相關(guān)的圓形域時(shí),該P(yáng)DOM 彈性模型即可簡化為通常的PD 或DH-PD 模型.并通過求解3 個(gè)典型問題: 桿的拉伸與波動(dòng)、亥姆霍茲方程和含孔板拉伸問題,說明本方法的準(zhǔn)確性、收斂性與穩(wěn)定性.

1 近場動(dòng)力學(xué)算子方法

考慮定義在M維空間Ω ?RM中的標(biāo)量函數(shù)

u(x)∈R,由N階泰勒展開可得

式中,η=u(x′)-u(x),x′=x+ξ ∈Hx ,Hx為點(diǎn)x的相互作用域,如圖1 所示.定義

圖1 PD 點(diǎn)的相互作用域Fig.1 Interaction domain of PD points

考慮定義在M維空間Ω ?RM中的Q階向量函數(shù)u(x)∈RQ,式(1)可擴(kuò)展為

式中,η=u(x′)-u(x).定義

當(dāng)Q=1時(shí),式(3)可簡化為式(1).

PD 函數(shù)構(gòu)造為

式中,wm為權(quán)函數(shù),為待定系數(shù).PD 函數(shù)具有正交性

式中,δqipi為克羅內(nèi)克符號.由式(3)與式(6)可得偏微分乘積的PDOM 構(gòu)型

將式(5)代入式(6)可得

式(8)可改寫為矩陣形式

由式(10)可得待定系數(shù),將其回代式(5),即可確定PD 函數(shù).當(dāng)Q=1時(shí),PDOM 即可退化為PDDO[17]或NOM[18-19].

接下來,具體給出兩種情況的PDOM 構(gòu)型.

情況 ?:M=2,N=1,Q=1

由式(2)可得

將式(11)代入式(4)可得

將式(11)代入式(5)可得

將式(11)和式(12)代入式(7)可得

將式(11)代入式(8)和式(9)可得

情況 П:M=2,N=2,Q=2

由式(2)可得

將式(11)和式(16)代入式(4)得

將式(11)和式(16)代入式(5)可得

將式(11)、式(16)和式(17)代入式(7)可得

將式(11)和式(16)代入式(8)和式(9)可得

2 PDOM 彈性模型

應(yīng)用PDOM 可以輕易地直接建立很多物理問題的非局部模型.限于篇幅,本節(jié)以二維線彈性固體為例,建立PDOM 彈性模型.

2.1 應(yīng)變能密度

應(yīng)變張量可表示為

式中,ui為位移.由式(21)和式(14),可得體應(yīng)變

由式(21)和式(19)可得

由式(22)和式(24),可得應(yīng)變能密度

式中,λ和 μ為拉梅常數(shù).

2.2 運(yùn)動(dòng)方程

系統(tǒng)的拉格朗日量為

式中,T為系統(tǒng)動(dòng)能,U為系統(tǒng)勢能,可表示為

式中,ρ為密度.拉格朗日方程可表示為

式中,Q(k)為廣義力,此處只含體力b(k).

式(22)和式(25a)的離散形式可表示為

式中,η(j)(k)=u(j)-u(k).由式(30)可得

由式(26)、式(28)和式(31)可得

將式(27)、式(28)和式(32)代入式(29)并轉(zhuǎn)換成積分形式可得

2.3 求解方案

式(30)可改寫為矩陣形式

其中

式中,N(k)為相互作用域H(k)中的離散點(diǎn)個(gè)數(shù).將式(36)代入式(26)可得

針對靜力學(xué)問題,可根據(jù)最小位能原理,建立能量泛函

式中,K為整體剛度矩陣,P為載荷列陣

其中,Ntotal為求解域的總離散點(diǎn)數(shù).

針對動(dòng)力學(xué)問題,可采用中心差分法,構(gòu)建時(shí)間積分

式中,Δt為時(shí)間步長.由式(34)和式(43)可得

3 與PD 模型的聯(lián)系

當(dāng)相互作用域Hx是由與位置相關(guān)的半徑 δx所定義的圓形域時(shí),式 (35 a)與DH-PD[16]中的內(nèi)力密度一致.本文模型即可退化為通常的DH-PD 模型.

當(dāng)相互作用域Hx是由與位置無關(guān)的半徑 δ所定義的圓形域時(shí),可得=Hx,式(35a)可表達(dá)為

式(45)與PD[1-2]中的內(nèi)力密度一致,本文模型即可退化為通常的近場動(dòng)力學(xué)模型.下文運(yùn)用PDOM彈性模型生成常用的OSB-PD 模型.

將w(1,0)=w(0,1)=ω(|ξ|)代入式(15),并考慮圓形作用域的Hx對稱性,可得

式(47)與OSB-PD 中的加權(quán)體積一致.將式(46)和式(13)代入式(23)可得

將式(48)代入式(22)可得

式(49)與OSB-PD 中的膨脹標(biāo)量狀態(tài)[2]一致.

將式(50)和式(18)代入式(25b)可得

將式(51)和式(25a)代入式(26)可得

式(52)與OSB-PD 中的彈性能密度一致.將式(48)和式(51)代入式(35b)可得

式(53)與OSB-PD 中的力矢量狀態(tài)[2]一致.

4 算例分析

本節(jié)以一維桿的拉伸與波動(dòng)、二維亥姆霍茲方程和含孔板拉伸為例,來說明PDOM 的求解能力.當(dāng)均勻離散求解域 Ω 時(shí),離散間距為 Δx,采用方形相互作用域

當(dāng)非均勻離散時(shí),選定N(k)個(gè)距離x(k)最近的點(diǎn)組成相互作用域H(k).定義.權(quán)函數(shù)設(shè)定為

式中,nw為權(quán)重指數(shù).

4.1 桿的拉伸與波動(dòng)

一維桿的運(yùn)動(dòng)方程可表示為

式中,u為軸向位移,A為截面面積,E為彈性模量.對于靜態(tài)拉伸問題,不考慮式(56)左端的慣性力.桿的左端固定,右端受到軸向力P作用.邊界條件可表示為

拉伸問題的精確解為

對于動(dòng)態(tài)波動(dòng)問題,設(shè)置式(58)為初始位移,初速度為零.邊界條件和初始條件可表示為

波動(dòng)問題的解析解為

該問題的能量泛函為

根據(jù)Q=1,2情況的PDOM,可構(gòu)造

通過拉格朗日方程,式(56)可改寫為

對于桿的拉伸問題,相關(guān)參數(shù)設(shè)定為EA=1,P=1,nw=3,N=2(Q=2時(shí)),N=1(Q=1時(shí)).

圖2 給出Q=1,2,m=2,3,Δx=0.02時(shí)桿拉伸的位移誤差對比.可以看出Q=1時(shí)的結(jié)果(此時(shí)得到的就是沒有引入其他額外的數(shù)值振蕩消除方法時(shí)的PDDO 或NOM 結(jié)果)有明顯的由零能模式引起的數(shù)值振蕩,而Q=2時(shí)的PDOM 結(jié)果非常穩(wěn)定.這表明本文方法可以從根本上有效避免零能模式.

圖3 給出Q=2,m=2,3,Δx=0.005,0.01,0.02和0.05 時(shí)桿拉伸的全局誤差收斂.收斂率分別為r=1.0301,1.1568,說明本方法具有良好的收斂性.其中,全局誤差由L2 范數(shù)衡量

圖3 桿拉伸的全局誤差收斂Fig.3 Convergence of global error for bar tension

對于桿的波動(dòng)問題,相關(guān)參數(shù)設(shè)定為 ρ=1,Δt=1.0×10-4,Δx=0.01,m=3.圖4 給出桿波動(dòng)的位移時(shí)空分布.圖5 給出桿波動(dòng)的位移分布對比.從圖中可以看出,PDOM 結(jié)果與精確解吻合良好.結(jié)果表明本方法可以高精度求解一維穩(wěn)態(tài)與瞬態(tài)問題.

圖4 桿波動(dòng)的位移時(shí)空分布Fig.4 Displacement in space and time for bar wave

圖5 桿波動(dòng)的位移分布對比Fig.5 Comparisons of displacement for bar wave

4.2 亥姆霍茲方程

二維亥姆霍茲方程可表達(dá)為

式中,k為波數(shù).精確解為

式中,J0是第一類零階貝塞爾函數(shù)

根據(jù)精確解,求解域四邊施加狄利克雷邊界條件.

該問題的能量泛函為

根據(jù)Q=2情況的PDOM,可構(gòu)造

將相關(guān)參數(shù)設(shè)定為m=3,nw=6,N=2,Δx=0.01.圖6 給出不同波數(shù)情況下亥姆霍茲方程的PDOM位移結(jié)果分布.圖7 給出不同波數(shù)情況下亥姆霍茲方程在 0 ≤x1≤1,x2=0上的位移分布對比.可以看出,PDOM 結(jié)果與精確解完全重合.表明了本方法求解二維問題的能力.

圖6 亥姆霍茲方程的PDOM 位移分布Fig.6 PDOM displacement for Helmholtz equation

圖7 亥姆霍茲方程的位移分布對比Fig.7 Comparisons of displacement for Helmholtz equation

4.3 含孔板拉伸

考慮均質(zhì)無限大的含孔板拉伸,如圖8(a)所示,P為拉伸均布力,a為圓孔半徑.該問題精確解為

圖8 均質(zhì)含孔板拉伸Fig.8 Homogeneous plate with a hole in tension

式中,r和 φ分別為極坐標(biāo)系極徑和極角,κ為體積模量,σφ為環(huán)向應(yīng)力.

設(shè)定方板邊長為1,圓孔半徑a=0.1,非均勻離散11152 個(gè)點(diǎn),如圖8(b)所示.彈性模量取1000,泊松比取0.3,均布力取P=1.根據(jù)精確解施加狄利克雷邊界條件.

相關(guān)參數(shù)設(shè)定為nw=6,N=2,N(k)=25.圖9給出含孔板拉伸的PDOM 位移結(jié)果分布.圖10給出含孔板拉伸在r=0.3上的位移與應(yīng)力分布對比.可以看出,PDOM 結(jié)果與精確解完全一致,表明本方法在非均勻離散情況下依然可以保證高精度.

圖9 含孔板拉伸的PDOM 位移分布Fig.9 PDOM displacement for plate with a hole in tension

圖10 含孔板拉伸的位移與應(yīng)力分布對比Fig.10 Comparisons of displacement and stress for plate with a hole in tension

5 結(jié)論

本文提出一種基于非局部思想求解物理學(xué)問題的方法,稱之為“近場動(dòng)力學(xué)算子方法”.給出了理論推導(dǎo)過程,與已有方法和模型進(jìn)行了對比,并以彈性力學(xué)問題為例,運(yùn)用變分原理和拉格朗日方程,建立了適用于靜/動(dòng)態(tài)彈性力學(xué)問題的PDOM 彈性模型,通過幾個(gè)典型算例進(jìn)行了驗(yàn)證.本文得出主要結(jié)論如下.

(1) PDOM 可將任意階局部微分及其乘積轉(zhuǎn)化為相應(yīng)的非局部積分形式.當(dāng)Q=1時(shí),PDOM 可退化為目前兩種受到關(guān)注的求解微分方程的非局部算子: PDDO 或NOM.

(2) 當(dāng)相互作用域取為與位置無關(guān)的半徑 δ或相關(guān)的半徑 δx所定義的圓形域時(shí),PDOM 彈性模型即可簡化為通常的PD 或DH-PD 模型.

(3) 通過求解桿的拉伸與波動(dòng)、亥姆霍茲方程和含孔板拉伸3 個(gè)經(jīng)典問題,表明本方法具有良好的計(jì)算精度、收斂性和數(shù)值穩(wěn)定性,能有效避免零能模式,且適用于非均勻離散.

PDOM 方法可為基于非局部思想求解微分方程組、分析不連續(xù)問題提供一種選擇.本文受篇幅限制,僅給出理論部分并以二維彈性力學(xué)問題為例開展了相關(guān)建模過程演示和分析.特別歡迎和期待相關(guān)領(lǐng)域廣大同仁采用PDOM 方法求解其他各種問題,特別是不連續(xù)問題.