基于態(tài)型近場動(dòng)力學(xué)的準(zhǔn)靜態(tài)數(shù)值仿真方法

孫璐妍，余 音，鄭曉玲

(1.上海飛機(jī)設(shè)計(jì)研究院，上海 201210；2.上海交通大學(xué) 航空航天學(xué)院，上海 200240)

1 引 言

2000年，Silling[1]提出了一套能夠求解多尺度不連續(xù)力學(xué)問題的近場動(dòng)力學(xué)理論，其優(yōu)勢得益于借助相同的數(shù)學(xué)模型描述物體從原子尺度到宏觀尺度的力學(xué)行為。它是在空間域中將對象離散為一系列物質(zhì)點(diǎn)，物質(zhì)點(diǎn)之間通過非局部作用產(chǎn)生聯(lián)系，物質(zhì)點(diǎn)的運(yùn)動(dòng)由空間積分方程表達(dá)，可以使用統(tǒng)一的數(shù)學(xué)模型求解多尺度力學(xué)行為，突破了空間微分方程及連續(xù)性假設(shè)的瓶頸，不連續(xù)現(xiàn)象隨求解自然發(fā)生[2]。

Silling等[3]提出的態(tài)型PD理論利用物質(zhì)點(diǎn)間的力狀態(tài)和變形狀態(tài)來構(gòu)造本構(gòu)關(guān)系，物質(zhì)點(diǎn)之間相互作用力和相對變形分解成和傳統(tǒng)力學(xué)概念類似的各向同性量和偏量，周圍環(huán)境對物質(zhì)點(diǎn)的影響通過各向同性量表達(dá)，PD模型參數(shù)和傳統(tǒng)材料力學(xué)參數(shù)通過能量推導(dǎo)建立關(guān)系，彌補(bǔ)了早期鍵型PD理論無法正確反應(yīng)某些材料屬性的不足。Silling[4]又推導(dǎo)了態(tài)型PD線性化理論，得到了模量狀態(tài)。在Silling的指引下，態(tài)型PD理論推廣應(yīng)用于瞬態(tài)動(dòng)力學(xué)問題的求解[5]、彈塑性[6]和準(zhǔn)脆性斷裂特征[7]本構(gòu)模型的描述。態(tài)型PD理論相應(yīng)的計(jì)算體系也進(jìn)一步完善，如特定條件下計(jì)算塑性變形的隱式算法[6]、態(tài)型PD材料失效判別準(zhǔn)則[8]和影響非局部作用的權(quán)函數(shù)設(shè)定[9]等。

真實(shí)的薄壁結(jié)構(gòu)往往存在釘孔連接等疲勞細(xì)節(jié)，而失穩(wěn)是薄壁結(jié)構(gòu)主要的失效模式之一，PD理論為后續(xù)研究失穩(wěn)對疲勞細(xì)節(jié)裂紋萌生與擴(kuò)展的影響開辟了一條新路。動(dòng)力松弛法是一種借助帶臨界阻尼的動(dòng)力學(xué)方程求解靜力學(xué)問題的方法，尤其適用于材料和幾何非線性問題的求解。Kilic等[10]以鍵型PD理論為基礎(chǔ)，運(yùn)用動(dòng)力松弛法求解了軸壓屈曲和熱屈曲問題。

本文將在態(tài)型PD線性化理論基礎(chǔ)上，推導(dǎo)節(jié)點(diǎn)剛度矩陣和結(jié)構(gòu)剛度矩陣的表達(dá)關(guān)系，改進(jìn) Kilic 等[10]基于鍵型PD理論的動(dòng)力松弛法，構(gòu)建態(tài)型PD理論中準(zhǔn)靜態(tài)仿真的數(shù)值計(jì)算方法，并將該方法用于金屬矩形薄板的軸向壓穩(wěn)定性預(yù)測。使用合理的離散方法建立薄板的態(tài)型PD模型以提高計(jì)算速率，通過矩形薄板拉壓的基準(zhǔn)數(shù)值實(shí)驗(yàn)與有限元仿真結(jié)果對比，驗(yàn)證薄板態(tài)型PD模型的正確性。最后，本文通過經(jīng)驗(yàn)公式證明了金屬平板軸壓穩(wěn)定性仿真結(jié)果的可靠性。

2 態(tài)型PD理論簡介

2.1 基本理論

假設(shè)物體占據(jù)某一空間域R，線性化態(tài)型PD理論[4]將域內(nèi)物質(zhì)點(diǎn)x的運(yùn)動(dòng)方程描述如下，

P(x)u(x,t)+b(x,t)

(1)

2.2 節(jié)點(diǎn)剛度矩陣與結(jié)構(gòu)剛度矩陣

設(shè)結(jié)構(gòu)離散后，節(jié)點(diǎn)總數(shù)為N，有Nx個(gè)節(jié)點(diǎn)qi(1≤i≤Nx)，滿足|qi-x|<2δ，將式(1)的積分展開，得到節(jié)點(diǎn)x的增量平衡方程如下，

(2)

從零開始編號，xi(0≤i

(0≤i

式中Ci j=C(xi,xj)，uj=u(xj)，bi=b(xi)，結(jié)構(gòu)剛度矩陣具有對稱性[4]。

3 準(zhǔn)靜態(tài)問題的數(shù)值方法

按照達(dá)朗貝爾原理，建立節(jié)點(diǎn)x的運(yùn)動(dòng)方程為

(t=nΔt) (4)

虛擬阻尼一般不高于系統(tǒng)臨界阻尼系數(shù)，cn(x)為對角矩陣，對角元素取值如下，

(6)

將式(4)進(jìn)行有限差分寫成分量形式即為

(7)

t+Δt時(shí)刻節(jié)點(diǎn)的位移為

初始位移場的選擇直接影響到算法的收斂速度，故節(jié)點(diǎn)的初始位移根據(jù)節(jié)點(diǎn)空間位置插值得到。

本文選取的初始位移場引起的系統(tǒng)不平衡力較小，在進(jìn)行系統(tǒng)平衡判定時(shí)，若選用絕對不平衡力收斂準(zhǔn)則，則迭代次數(shù)過少影響求解精度，若選用相對不平衡力收斂準(zhǔn)則，則迭代次數(shù)過多影響求解速度。本文融合了相對及絕對兩種準(zhǔn)則的優(yōu)勢，得到了一種修正的收斂準(zhǔn)則為

(8)

4 數(shù)值仿真

本文分析借助了結(jié)構(gòu)分析軟件SBPD(State -Based Peridynamics)，軟件求解器基于VC++語言，通過統(tǒng)一計(jì)算設(shè)備架構(gòu)技術(shù)CUDA(Compute Unified Device Architecture)和共享存儲(chǔ)并行編程語言擴(kuò)展OpenMP(Open Multi-Processing)，提供了CPU+GPU協(xié)同和CPU多核兩種并行計(jì)算方式。

4.1 基準(zhǔn)數(shù)值試驗(yàn)

以矩形薄板拉壓作為基準(zhǔn)試驗(yàn)，將SBPD與ABAQUS的仿真結(jié)果進(jìn)行對比。矩形薄板長為 300 mm，寬為200 mm，厚為2 mm，彈性模量為 69 GPa，泊松比為0.33，載荷形式為沿長度方向的位移。平板拉伸的邊界條件為雙邊簡支，雙邊自由；壓縮的邊界條件為四邊簡支。拉和壓施加的位移大小均為1.2 mm。

在PD理論的數(shù)值仿真中，物體通常劃分成大小一致的正方體晶格，節(jié)點(diǎn)設(shè)在晶格的正中，對于薄板問題，若厚度方向只劃分一層晶格，會(huì)使節(jié)點(diǎn)剛度矩陣中與厚度方向有關(guān)的量為0，導(dǎo)致方程無法求解，不僅如此，薄板結(jié)構(gòu)發(fā)生屈曲時(shí),面內(nèi)的應(yīng)力應(yīng)變在厚度方向上是變化的，因此在厚度方向上至少劃分兩層晶格。而薄板結(jié)構(gòu)的長寬方向的尺寸遠(yuǎn)大于厚度方向的尺寸，正方體晶格的離散方法會(huì)導(dǎo)致模型節(jié)點(diǎn)的數(shù)量過于龐大，影響求解效率。本文將結(jié)構(gòu)離散成方板晶格，在板厚方向劃分了兩層晶格，晶格長寬為10 mm，厚為1 mm，晶格數(shù)量為1200，方板晶格的對角線長為14.2 mm，即為節(jié)點(diǎn)的近場范圍，節(jié)點(diǎn)2倍近場范圍內(nèi)一般有49個(gè)節(jié)點(diǎn)與之關(guān)聯(lián)。與之相比，常規(guī)離散方法至少需要?jiǎng)澐?20000個(gè)正方體晶格，節(jié)點(diǎn)2倍近場范圍內(nèi)一般有73個(gè)節(jié)點(diǎn)與之關(guān)聯(lián)，節(jié)點(diǎn)剛度矩陣計(jì)算需要經(jīng)過兩次關(guān)聯(lián)節(jié)點(diǎn)的遍歷，其計(jì)算時(shí)間與關(guān)聯(lián)節(jié)點(diǎn)數(shù)量的平方正相關(guān)，整個(gè)模型的求解時(shí)間與節(jié)點(diǎn)數(shù)量和節(jié)點(diǎn)剛度矩陣計(jì)算時(shí)間正相關(guān)。綜上所述，本文方法將求解效率提升了約200倍。

在ABAQUS中做了類似的離散處理，單元類型為三維實(shí)體減縮積分單元C3D8R。分析得到的平板拉伸沿y方向的位移場云圖如圖1所示。

ABAQUS得到的拉壓外力值分別為112.2 kN和-124.6 kN，SBPD得到的拉壓外力值為119.9 kN和-131.9 kN，SBPD的結(jié)果比ABAQUS大6%～7%。這是由于兩者離散方法的差異造成的，PD模型的實(shí)際加載端距為290 mm，相同的位移載荷下，需作用的外力也更大，因此兩者結(jié)果差異的趨勢也是合理的，證明了金屬矩形薄板態(tài)型PD模型的正確性。

圖1 平板拉伸沿y方向的位移場云圖(單位：mm)

圖2 平板壓縮沿x方向的位移場云圖(單位：mm)

4.2 金屬平板軸壓穩(wěn)定性分析

對于對稱結(jié)構(gòu)而言，數(shù)值仿真穩(wěn)定性問題需要先引入初始幾何缺陷作為擾動(dòng)，本文引入的初始缺陷形態(tài)為ABAQUS屈曲分析的第1階模態(tài)位移場，設(shè)置的縮比因子為0.2(幅值大小為1/10的平板厚度)。

仿真得到的軸向載荷位移變化曲線如圖3所示，可以看出，平板在點(diǎn)a發(fā)生了初始屈曲，載荷為9.6 kN，軸壓剛度突變，與初始值比下降了40%，此時(shí)平板仍然具有承載能力；初始屈曲后軸壓剛度緩慢降低，在點(diǎn)b發(fā)生了二次分支型屈曲，載荷為35.2 kN，此時(shí)的平衡狀態(tài)并不穩(wěn)定；很快在點(diǎn)c發(fā)生了三次跳躍型屈曲，載荷為39.6 kN，在點(diǎn)d達(dá)到新的平衡狀態(tài)，從點(diǎn)c到點(diǎn)d的過程中，平板的軸壓剛度突變，出現(xiàn)不連續(xù)跳躍，點(diǎn)c和點(diǎn)d的軸壓剛度分別約為初始值的40%和30%；此后平板繼續(xù)承載，面內(nèi)剛度保持緩慢下降。平板中心點(diǎn)撓度隨載荷變化的曲線如圖4所示，屈曲模態(tài)變化如圖5所示，可以看出，屈曲前中心點(diǎn)撓度變化很小，初始屈曲后迅速變大，在中心區(qū)域出現(xiàn)一個(gè)半波波形，半波慢慢鼓起，在波峰撓度無法繼續(xù)變大時(shí)，發(fā)生二次屈曲，一個(gè)半波擴(kuò)展為兩個(gè)方向一致的半波，然而這個(gè)狀態(tài)并不穩(wěn)定，迅速發(fā)生三次屈曲，中心點(diǎn)撓度發(fā)生了突變，在兩個(gè)波之間出現(xiàn)第三個(gè)方向相反的半波。

圖3 載荷位移曲線

圖4 中心點(diǎn)撓度隨載荷變化的曲線

圖5 屈曲模態(tài)

常用的矩形平板壓縮穩(wěn)定性臨界應(yīng)力σc r公式[11]如下,

(9)

式中 楊氏模量E取69×103MPa，泊松比μe取 0.33，平板厚度δ取2 mm，加載邊寬度b取為 200 mm，壓縮臨界應(yīng)力系數(shù)Kc與平板長寬比和邊界約束有關(guān)，對于本文問題，Kc取4，計(jì)算得到σc r=25.47 MPa，乘以加載邊截面面積得到臨界壓力10.190 kN，數(shù)值仿真結(jié)果與該值相差約5.7%，證明仿真得到的初始屈曲結(jié)果是可靠的。

從仿真結(jié)果可知，矩形平板初始屈曲后仍具備一定的承載能力，且三倍初始屈曲載荷內(nèi)，平板的屈曲模態(tài)形式不會(huì)發(fā)生變化，這與矩形薄板軸壓穩(wěn)定性的試驗(yàn)研究成果一致[12,13]。試驗(yàn)研究發(fā)現(xiàn)金屬矩形薄板是否發(fā)生多次屈曲與矩形薄板的長寬比和邊界條件有關(guān)，而最終破壞與金屬材料進(jìn)入塑性狀態(tài)有關(guān)[13]。

5 總 結(jié)

本文提出了一種基于PD狀態(tài)理論準(zhǔn)靜態(tài)問題仿真的數(shù)值方法。

(1) 針對PD薄板模型提出了合理的離散方法以提高計(jì)算效率。

(2) 推導(dǎo)了節(jié)點(diǎn)剛度矩陣和結(jié)構(gòu)剛度矩陣的表達(dá)關(guān)系，改進(jìn)了Kilic[9]基于PD鍵理論的動(dòng)力松弛法。

(3) 結(jié)合絕對準(zhǔn)則和相對準(zhǔn)則，提出了一種新的修正迭代收斂準(zhǔn)則，并將該方法成功運(yùn)用于薄板拉壓及軸壓穩(wěn)定性的數(shù)值仿真中。

目前，態(tài)型PD理論的研究工作還處于起步階段，相信待理論框架和計(jì)算體系成熟之后，將在不連續(xù)力學(xué)領(lǐng)域發(fā)揮更大的作用。