孫璐妍，余 音，鄭曉玲

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

1 引 言

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

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

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

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

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

2.1 基本理論

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

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

(1)

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

設(shè)結(jié)構(gòu)離散后，節(jié)點總數(shù)為N，有Nx個節(jié)點qi(1≤i≤Nx)，滿足|qi-x|<2δ，將式(1)的積分展開，得到節(jié)點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ù)值方法

按照達朗貝爾原理，建立節(jié)點x的運動方程為

(t=nΔt) (4)

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

(6)

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

(7)

t+Δt時刻節(jié)點的位移為

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

本文選取的初始位移場引起的系統(tǒng)不平衡力較小，在進行系統(tǒng)平衡判定時，若選用絕對不平衡力收斂準(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)一計算設(shè)備架構(gòu)技術(shù)CUDA(Compute Unified Device Architecture)和共享存儲并行編程語言擴展OpenMP(Open Multi-Processing)，提供了CPU+GPU協(xié)同和CPU多核兩種并行計算方式。

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

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

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

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

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

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

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

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

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

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

圖3 載荷位移曲線

圖4 中心點撓度隨載荷變化的曲線

圖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，計算得到σc r=25.47 MPa，乘以加載邊截面面積得到臨界壓力10.190 kN，數(shù)值仿真結(jié)果與該值相差約5.7%，證明仿真得到的初始屈曲結(jié)果是可靠的。

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

5 總 結(jié)

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

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

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

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

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