基于偶應(yīng)力理論的高階彈塑性本構(gòu)模型的無(wú)網(wǎng)格法

王瑞昌， 孫玉周

(中原工學(xué)院 建筑工程學(xué)院， 河南 鄭州 450007)

近年來(lái)，各種特殊材料以及微/納米量級(jí)新型構(gòu)件被越來(lái)越多地應(yīng)用到工程實(shí)際中，偶應(yīng)力理論[1]和應(yīng)變梯度理論[2]等高階連續(xù)理論被廣泛用在微/納米結(jié)構(gòu)或者需要考慮微尺度影響的宏觀結(jié)構(gòu)的研究中[3]。BORST等率先提出了基于Von-Mises屈服準(zhǔn)則的Cosserat彈塑性理論[4]；FLECK等分別從統(tǒng)計(jì)存儲(chǔ)位錯(cuò)和幾何必需位錯(cuò)兩個(gè)角度出發(fā),改進(jìn)了Cosserat彈塑性模型[5]；李錫夔等基于Drucker-Prager屈服準(zhǔn)則給出了Cosserat彈塑性理論的塑性一致性算法[6]；冀賓等基于偶應(yīng)力彈塑性理論，對(duì)軟化行為的剪切帶問(wèn)題進(jìn)行了有限元數(shù)值模擬[7]；張建成等利用非線性Cosserat擴(kuò)展模型對(duì)層狀巖體的洞室開(kāi)挖過(guò)程進(jìn)行了研究[8]。整體上說(shuō)，由于數(shù)值離散時(shí)要求位移函數(shù)滿足C1連續(xù)性，高階連續(xù)理論、特別是高階彈塑性理論的數(shù)值模擬仍然是工程數(shù)值仿真領(lǐng)域的研究熱點(diǎn)。

無(wú)網(wǎng)格法是一種新型數(shù)值計(jì)算方法，與傳統(tǒng)有限元法相比，具有形函數(shù)滿足高階連續(xù)特征、不受網(wǎng)格限制等特點(diǎn)，特別適合大變形、位移不連續(xù)、裂紋擴(kuò)展和高速碰撞等問(wèn)題的數(shù)值模擬[9]。本文嘗試用無(wú)網(wǎng)格法建立基于偶應(yīng)力理論的高階彈塑性本構(gòu)模型的數(shù)值計(jì)算框架 ，用Fortran語(yǔ)言編寫計(jì)算程序，對(duì)懸臂梁進(jìn)行數(shù)值模擬，以研究無(wú)網(wǎng)格法在高階彈塑性問(wèn)題數(shù)值模擬中的應(yīng)用。

1 基于偶應(yīng)力理論的高階彈塑性本構(gòu)模型

1.1 本構(gòu)模型

偶應(yīng)力理論的幾何方程為

(1)

式中：εij為常規(guī)應(yīng)變，ui、uj為位移，下標(biāo)“，”表示對(duì)其后坐標(biāo)求偏導(dǎo)，χij為微曲率，ωj為微轉(zhuǎn)角。

增量形式的本構(gòu)方程為

(2)

式(2)也可寫為

(3)

不記體力和體力偶時(shí)，平衡方程為

σij,j=0，mij,j+ejklσkl

(4)

式中：ejkl為排列算子。

靜力邊界條件為

(5)

位移邊界條件為

(6)

1.2 增量彈塑性陣

材料屈服后可把應(yīng)變?cè)隽糠譃閺椥栽隽亢退苄栽隽績(jī)蓚€(gè)部分，即

(7)

對(duì)于塑性變形而言，一應(yīng)力增量可對(duì)應(yīng)不同的塑性變形增量。根據(jù)相應(yīng)的屈服準(zhǔn)則和流動(dòng)法則，可得

[D]ep=[D]e-[D]p=

(8)

對(duì)于平面應(yīng)力問(wèn)題，應(yīng)力增量和應(yīng)變?cè)隽糠謩e為

(9)

等效應(yīng)力可通過(guò)應(yīng)力偏量寫為

(10)

于是

(11)

偶應(yīng)力理論的彈性矩陣為[10]

(12)

式中：E為彈性模量，v為泊松比。

將式(10)-式(12)代入式(8)并經(jīng)過(guò)整理，最終得到平面應(yīng)力問(wèn)題的彈塑性陣顯式表達(dá)式為

(13)

2 無(wú)網(wǎng)格法計(jì)算框架

控制方程的弱形式可寫為

(14)

在無(wú)網(wǎng)格法中，域內(nèi)任意一點(diǎn)x的位移增量可由移動(dòng)最小二乘近似表示為[10-11]

(15)

域內(nèi)任意一點(diǎn)x處的應(yīng)變?cè)隽靠山茷?/p>

(16)

式中：B為幾何矩陣，表達(dá)式如下

(17)

則域內(nèi)任意一點(diǎn)x處的應(yīng)力增量為

(18)

由式(14)可得離散控制方程為

(19)

由于位移、轉(zhuǎn)角和曲率均由節(jié)點(diǎn)位移來(lái)插值，位移和轉(zhuǎn)角邊界條件可采用罰數(shù)法施加。

3 數(shù)值算例及分析

圖1所示為一端受集中力作用的懸臂梁，其幾何尺寸取長(zhǎng)度L=8.0 m,高度h=1.0 m，集中力F=1 N，彈性模量E=1.0×105Pa，泊松比v=0.25，屈服極限σs=25 Pa。采用線性強(qiáng)化彈塑性模型，H′=0.25E，切線模量E′=0.2E，材料服從Mises屈服條件。

圖1 受集中力作用的懸臂梁模型

采用無(wú)網(wǎng)格法計(jì)算時(shí)，均布節(jié)點(diǎn)數(shù)為81×11，相鄰節(jié)點(diǎn)間距在x和y方向上均為10 cm，背景網(wǎng)格數(shù)為80×10，與節(jié)點(diǎn)布置一致，背景網(wǎng)格內(nèi)高斯點(diǎn)個(gè)數(shù)為3×3個(gè)。采用ANSYS建模時(shí)，節(jié)點(diǎn)及網(wǎng)格布置方案與無(wú)網(wǎng)格法相同，加載步數(shù)均為100。

表1給出了經(jīng)典力學(xué)彈性解、ANSYS彈塑性解和本文方法(尺度因子l=0)計(jì)算得到的中性層各點(diǎn)撓度值?？梢钥闯觯瑹o(wú)網(wǎng)格法計(jì)算結(jié)果與有限元軟件ANSYS計(jì)算結(jié)果相差很小,且都大于經(jīng)典力學(xué)彈性解。這說(shuō)明本文建立的數(shù)值計(jì)算框架可以獲得較好的效果。

圖2所示為Von-Mises應(yīng)力云圖。當(dāng)構(gòu)件中某點(diǎn)的Mises應(yīng)力超過(guò)屈服極限(σs=25 Pa)后即進(jìn)入塑性階段，從圖中可明顯觀察到塑性區(qū)。

表1 一端受集中力作用的懸臂梁豎向位移 m

圖2 無(wú)網(wǎng)格法數(shù)值模擬懸臂梁Von-mises應(yīng)力云圖

圖3所示為懸臂梁右端中點(diǎn)的撓度與荷載的關(guān)系曲線圖。從圖中可以看出，在增量加載過(guò)程中，當(dāng)荷載大于其彈性極限值(P>0.6 N)時(shí)，懸臂梁由彈性變形階段進(jìn)入塑性變形階段。

圖3 梁右端中點(diǎn)的撓度與荷載的關(guān)系曲線圖

為了分析尺度因子對(duì)懸臂梁彈塑性彎曲變形的影響，令尺度因子l在0～2.0b(b為厚度)之間變化，分別計(jì)算懸臂梁的中性層各點(diǎn)撓度值，并將其與有限元軟件ANSYS計(jì)算結(jié)果進(jìn)行比較，如圖4所示。可以看出，當(dāng)梁厚度遠(yuǎn)遠(yuǎn)大于材料的特征長(zhǎng)度尺寸(如l<0.01b)時(shí)，梁的撓度曲線接近于ANSYS計(jì)算得到的彈塑性解，尺度效應(yīng)可以忽略。隨著材料特征長(zhǎng)度尺寸的不斷增大，由于存在偶應(yīng)力，其彎曲剛度顯著增加，對(duì)梁撓度的影響也隨之增大。

圖4 不同l/b條件下懸臂梁中性層撓度值

4 結(jié)語(yǔ)

本文以基于偶應(yīng)力理論的高階彈塑性本構(gòu)模型為對(duì)象，建立了無(wú)網(wǎng)格法數(shù)值計(jì)算框架，用Fortran語(yǔ)言編寫計(jì)算程序，并對(duì)懸臂梁的彈塑性彎曲變形進(jìn)行了數(shù)值模擬，研究了無(wú)網(wǎng)格法在高階彈塑性數(shù)值模擬中的應(yīng)用。結(jié)果表明，本文方法不僅可以求解高階彈性問(wèn)題，對(duì)于高階彈塑性問(wèn)題也是非常有效的。當(dāng)二維懸臂梁厚度接近材料特征長(zhǎng)度尺寸時(shí)，由于存在偶應(yīng)力，其彎曲剛度較之前有明顯的增加，尺度效應(yīng)較為明顯；當(dāng)其厚度遠(yuǎn)遠(yuǎn)大于材料的特征長(zhǎng)度尺寸時(shí)，材料尺度效應(yīng)非常弱，可以忽略不記。