黏塑性本構(gòu)計(jì)算的穩(wěn)定性分析*

劉明濤，李永池，胡秀章，章 杰

(1.中國(guó)工程物理研究院流體物理研究所,四川 綿陽(yáng) 621999; 2.中國(guó)科學(xué)技術(shù)大學(xué)近代力學(xué)系,安徽 合肥 230027)

黏塑性本構(gòu)計(jì)算的穩(wěn)定性分析*

劉明濤1,2，李永池2，胡秀章2，章 杰2

(1.中國(guó)工程物理研究院流體物理研究所,四川 綿陽(yáng) 621999; 2.中國(guó)科學(xué)技術(shù)大學(xué)近代力學(xué)系,安徽 合肥 230027)

提出本構(gòu)方程計(jì)算方法的穩(wěn)定性問(wèn)題，針對(duì)黏塑性本構(gòu)計(jì)算的顯式精確算法的穩(wěn)定性進(jìn)行分析，發(fā)現(xiàn)該算法并非無(wú)條件穩(wěn)定，使用小擾動(dòng)方法給出了其計(jì)算穩(wěn)定的必要條件，穩(wěn)定性條件對(duì)數(shù)值計(jì)算中的時(shí)間步長(zhǎng)提出限制要求。通過(guò)有限元算例驗(yàn)證了分析的正確性，計(jì)算結(jié)果也表明理論推導(dǎo)得到的穩(wěn)定性公式能夠準(zhǔn)確預(yù)測(cè)滿(mǎn)足計(jì)算穩(wěn)定性條件要求的最大時(shí)間步長(zhǎng)與各參數(shù)之間關(guān)系。

本構(gòu)關(guān)系；算法穩(wěn)定性；黏塑性本構(gòu)；顯式精確算法

一直以來(lái)，塑性和黏塑性材料的本構(gòu)方程及其數(shù)值算法的研究是計(jì)算力學(xué)的核心問(wèn)題。D.C.Drucker等[1-2]提出了著名的Drucker公設(shè)，該公設(shè)是經(jīng)典塑性理論的基石。利用Drucker公設(shè)可直接推導(dǎo)出塑性流動(dòng)的正交法則和屈服面的外凸性，但其缺點(diǎn)是只適用于穩(wěn)定材料。A.C.Palmer[3]指出針對(duì)軟化的非穩(wěn)定材料也可以得到塑性流動(dòng)的正交法則和屈服面的外凸性。李永池等[4]進(jìn)一步發(fā)展了A.C.Palmer的思想，提出了廣義Drucker公設(shè)，將穩(wěn)定材料和非穩(wěn)定材料統(tǒng)一在一個(gè)框架之內(nèi)。

Drucker公設(shè)推導(dǎo)出的塑性流動(dòng)法則是塑性本構(gòu)計(jì)算的基石，其指出塑性應(yīng)變必須沿著當(dāng)前加載面的法線方向發(fā)展。迄今為止，有關(guān)沖擊動(dòng)力學(xué)問(wèn)題的計(jì)算程序中，本構(gòu)算法大部分采用傳統(tǒng)的半徑回歸法[5-7]，這是人們最初針對(duì)理想塑性材料提出的一種本構(gòu)更新算法。對(duì)于具有應(yīng)變率效應(yīng)的黏塑性材料，采用這種算法會(huì)帶來(lái)較大的誤差，尤其是當(dāng)材料由彈性狀態(tài)進(jìn)入屈服狀態(tài)時(shí)。針對(duì)此問(wèn)題，李永池等[8]提出了一種新的顯式本構(gòu)計(jì)算方法，稱(chēng)為顯式精確算法?；趶V義Drucker公設(shè)，李永池等[8]理論推導(dǎo)指出當(dāng)黏塑性材料進(jìn)入屈服狀態(tài)后，利用材料的實(shí)時(shí)狀態(tài)量即可唯一確定塑性流變過(guò)程中的塑性流動(dòng)因子。若在一階精度下，利用上時(shí)刻的狀態(tài)量即可求出現(xiàn)時(shí)刻的塑性流動(dòng)因子，繼而可計(jì)算現(xiàn)時(shí)刻的塑性變形，而后利用胡克定律可得現(xiàn)時(shí)刻的應(yīng)力增量。

通常情況下，數(shù)值計(jì)算的穩(wěn)定性是針對(duì)連續(xù)方程、動(dòng)量方程和能量方程的差分格式而言。本文中探討本構(gòu)方程計(jì)算方法的穩(wěn)定性問(wèn)題。對(duì)李永池等[8]發(fā)展的顯式精確算法進(jìn)行穩(wěn)定性分析，并通過(guò)對(duì)單個(gè)單元的有限元算例驗(yàn)證其正確性。

1 顯式精確算法及其計(jì)算穩(wěn)定性

顯式精確算法的詳細(xì)理論推導(dǎo)過(guò)程見(jiàn)文獻(xiàn)[8]，文獻(xiàn)[8]中是從一般性的熱黏塑性本構(gòu)關(guān)系出發(fā)推導(dǎo)的，具有一定的普適性。本文中對(duì)其推導(dǎo)過(guò)程進(jìn)行一定的簡(jiǎn)化，從目前常見(jiàn)的Mises屈服準(zhǔn)則下的黏塑性本構(gòu)關(guān)系出發(fā)推導(dǎo)。

1.1顯式精確算法

針對(duì)Mises類(lèi)黏塑性本構(gòu)關(guān)系，屈服函數(shù)可設(shè)為：

(1)

其中，Mises等效應(yīng)力和等效塑性應(yīng)變率計(jì)算公式分別為：

(2)

(3)

根據(jù)廣義Drucker公設(shè)可知，塑性流動(dòng)的正交法則為：

(4)

將式(1)～(2)代入式(4)得：

(5)

將式(5)代入式(3)，得：

(6)

由黏塑性本構(gòu)方程式(1)，可反解出等效塑性應(yīng)變率為：

(7)

式(6)、(7)聯(lián)立，得：

(8)

由式(8)可知，塑性流動(dòng)因子可由材料當(dāng)前的應(yīng)力狀態(tài)唯一確定。

將式(8)代入式(5)得：

(9)

根據(jù)材料的胡克定律，有：

(10)

(11)

1.2計(jì)算穩(wěn)定性分析

式(11)是一個(gè)張量表達(dá)式，不利于計(jì)算的穩(wěn)定性分析。下面對(duì)此進(jìn)行簡(jiǎn)化，選擇一種特殊情況進(jìn)行計(jì)算穩(wěn)定性理論分析。

(12)

(13)

(14)

將式(12)～(14)代入式(11)，得：

(15)

由上式可計(jì)算得：

(16)

(17)

(18)

將函數(shù)g(s(n)+Δs)在s(n)處進(jìn)行泰勒展開(kāi)，可得：

g(s(n)+Δs(n))=g(s(n))+Δs(n)g′(s(n))+O(Δs(n)2)

(19)

將式(19)代入式(18)，并略去二階小量，得：

Δs(n+1)=Δs(n)[1-3Gdtg′(s(n))]

(20)

計(jì)算穩(wěn)定性要求擾動(dòng)的放大因子小于1，即|1-3Gdtg′(s(n))|<1，解得：

0<3Gdtg′(s(n))<2

(21)

式(21)即為推導(dǎo)得出的 顯式精確算法的計(jì)算穩(wěn)定性條件，顯式精確算法的計(jì)算穩(wěn)定性對(duì)計(jì)算時(shí)間步長(zhǎng)dt提出了限制要求。還需特別指出的是，式(21)是在式(12)、(13)的特殊情況下推導(dǎo)得到的，并不能保證本構(gòu)計(jì)算絕對(duì)穩(wěn)定。

若選擇黏塑性材料的具體屈服準(zhǔn)則形式為：

(22)

將式(22)代入式(21)，可得本構(gòu)方程的計(jì)算穩(wěn)定性要求為：

(23)

傳統(tǒng)Courant穩(wěn)定性條件為：

(24)

式中：L為單元特征長(zhǎng)度，C為材料絕熱聲速，α為安全因數(shù)，通常取α=0.9。式(21)、(23)表明，本構(gòu)計(jì)算的穩(wěn)定性與材料的本構(gòu)參數(shù)和實(shí)時(shí)塑性應(yīng)變率密切相關(guān)，而傳統(tǒng)Courant穩(wěn)定性條件僅與材料的聲速和單元的尺寸相關(guān)。

2 數(shù)值計(jì)算算例

四邊形單元的運(yùn)動(dòng)過(guò)程設(shè)定為：在整個(gè)變形過(guò)程中4個(gè)節(jié)點(diǎn)在x方向均固定不動(dòng)，同時(shí)節(jié)點(diǎn)1、2在y方向也固定不動(dòng)，首先，節(jié)點(diǎn)3和4以速度-50 m/s沿著y方向勻速運(yùn)動(dòng)20 μs，四邊形單元的y方向長(zhǎng)度由20 mm壓縮至19 mm；而后，節(jié)點(diǎn)3和節(jié)點(diǎn)4再以50 m/s勻速運(yùn)動(dòng)返回至初始位置，四邊形單元的y方向長(zhǎng)度由19 mm回復(fù)至初始時(shí)的20 mm。

四邊形單元在上述變形過(guò)程中經(jīng)歷的加卸載路徑較復(fù)雜，如圖1所示。從狀態(tài)1到狀態(tài)2，經(jīng)歷了彈性加載、塑性加載；從狀態(tài)2到狀態(tài)3，經(jīng)歷了彈性卸載、反向彈性加載、反向塑性加載，共5個(gè)階段。在第2階段和第5階段，材料發(fā)生了塑性應(yīng)變，塑性應(yīng)變率均約為1 700 s-1。根據(jù)顯式精確算法的穩(wěn)定性條件式(2)，可知時(shí)間步長(zhǎng)dt需滿(mǎn)足：0

2.1不同時(shí)間步長(zhǎng)算例

2.1.1時(shí)間步長(zhǎng)dt=0.3×10-7s時(shí)

計(jì)算結(jié)果如圖2所示，黑色實(shí)線為等效應(yīng)力歷史曲線，藍(lán)色點(diǎn)劃線為等效塑性應(yīng)變歷史曲線。從圖中可以看出，當(dāng)選取的時(shí)間步長(zhǎng)滿(mǎn)足計(jì)算穩(wěn)定性條件式(23)時(shí)，利用顯式精確算法計(jì)算的結(jié)果能夠準(zhǔn)確描述四邊形單元在整個(gè)變形過(guò)程中所經(jīng)歷的復(fù)雜加卸載過(guò)程。

(1)彈性加載階段：隨著變形的增加，等效應(yīng)力勻速增加，沒(méi)有發(fā)生塑性變形。

(2)塑性加載階段：由于節(jié)點(diǎn)3和節(jié)點(diǎn)4以勻速運(yùn)動(dòng)，因此在整個(gè)加載過(guò)程中等效塑性應(yīng)變率保持不變，因此根據(jù)本構(gòu)關(guān)系式(22)，等效應(yīng)力也保持不變，而圖2中等效應(yīng)力在該階段呈現(xiàn)1個(gè)平臺(tái)段，計(jì)算結(jié)果正確。

(3)彈性卸載階段： 20 μs以后，節(jié)點(diǎn)3和節(jié)點(diǎn)4沿原路徑勻速返回，材料進(jìn)入彈性卸載階段，該過(guò)程中不產(chǎn)生塑性應(yīng)變?cè)隽浚虼说刃苄詰?yīng)變呈現(xiàn)平臺(tái)段，等效應(yīng)力由材料的屈服點(diǎn)勻速降低至零，計(jì)算結(jié)果正確。

(4)反向彈性加載階段：當(dāng)?shù)刃?yīng)力降低為零后，重新開(kāi)始勻速增大，材料進(jìn)入反向彈性加載階段，在此階段等效塑性應(yīng)變不增大，呈現(xiàn)平臺(tái)段，計(jì)算結(jié)果正確。

(5)反向塑性加載：當(dāng)材料由于反向加載再次進(jìn)入屈服后，塑性變形重新開(kāi)始累積，由于在該階段節(jié)點(diǎn)3和節(jié)點(diǎn)4勻速運(yùn)動(dòng)，因此塑性應(yīng)變率保持為恒定值，所以塑性應(yīng)變線性增大，等效應(yīng)力出現(xiàn)第2個(gè)平臺(tái)段，計(jì)算結(jié)果正確。

2.1.2時(shí)間步長(zhǎng)dt=1.0×10-7s時(shí)

該時(shí)間步長(zhǎng)不滿(mǎn)足顯式精確算法計(jì)算得出的穩(wěn)定性條件，計(jì)算結(jié)果如圖3所示?？梢钥闯?，當(dāng)材料屈服后，計(jì)算得到的等效應(yīng)力出現(xiàn)了上下抖動(dòng)現(xiàn)象，此時(shí)顯式精確算法的計(jì)算不穩(wěn)定，計(jì)算結(jié)果錯(cuò)誤。

2.1.3時(shí)間步長(zhǎng)dt=1.5×10-7s時(shí)

進(jìn)一步增大了時(shí)間步長(zhǎng)，其遠(yuǎn)不滿(mǎn)足本構(gòu)計(jì)算的穩(wěn)定性條件。計(jì)算結(jié)果如圖4所示。從圖4可以看出，當(dāng)材料屈服后，等效塑性應(yīng)變歷史曲線呈臺(tái)階式上升，四邊形單元在彈性與塑性狀態(tài)之間來(lái)回跳動(dòng)，等效應(yīng)力也變得極不穩(wěn)定，計(jì)算結(jié)果與真實(shí)值相差巨大。但是當(dāng)采用傳統(tǒng)的近似算法(半徑回歸法)并仍取時(shí)間步長(zhǎng)為dt= 1.5×10-7s時(shí)，數(shù)值計(jì)算結(jié)果穩(wěn)定收斂,具體結(jié)果如圖5所示。

綜合上述算例可以看出，顯式精確算法的確存在計(jì)算穩(wěn)定性問(wèn)題，取同樣的時(shí)間步長(zhǎng)dt=1.5×10-7s，當(dāng)采用傳統(tǒng)的近似算法(半徑回歸法)得到了正確的結(jié)果，而采用顯式精確算法時(shí)結(jié)果失穩(wěn)。當(dāng)時(shí)間步長(zhǎng)減小至滿(mǎn)足本構(gòu)計(jì)算穩(wěn)定性條件式(21)時(shí)，顯式精確算法給出的計(jì)算結(jié)果也穩(wěn)定收斂。

2.2最大時(shí)間步長(zhǎng)與材料參數(shù)和塑性應(yīng)變率關(guān)系

進(jìn)一步分析滿(mǎn)足本構(gòu)計(jì)算穩(wěn)定的最大時(shí)間步長(zhǎng)與材料參數(shù)和塑性應(yīng)變率關(guān)系，分別研究各個(gè)參數(shù)與最大時(shí)間步長(zhǎng)的關(guān)系。分為4組計(jì)算，每組只變化1個(gè)變量，其余參數(shù)的值取表1中的參數(shù)值，數(shù)值模擬得到滿(mǎn)足本構(gòu)計(jì)算穩(wěn)定的最大時(shí)間步長(zhǎng)，其隨各參數(shù)的變化如圖6～9所示。

可以看出，本構(gòu)計(jì)算穩(wěn)定性準(zhǔn)則式(21)、(23)與模擬結(jié)果符合的較好，最大時(shí)間步長(zhǎng)與各參量的依賴(lài)關(guān)系為：與應(yīng)變率敏感因子β成正比、與靜態(tài)條件下屈服強(qiáng)度Y*成正比、與剪切模量G成反比、與塑性應(yīng)變率成反比。但需要特別指出的是，本構(gòu)計(jì)算穩(wěn)定性準(zhǔn)則式(21)只是顯式精確算法計(jì)算穩(wěn)定的必要性條件。在實(shí)際工程計(jì)算過(guò)程中，為提高其可靠性，可以取一個(gè)安全系數(shù)。

3 結(jié) 論

首先對(duì)最常用的Mises類(lèi)黏塑性材料，重新推導(dǎo)顯式精確算法的計(jì)算公式和流程，然后通過(guò)理論推導(dǎo)得到顯式精確算法的穩(wěn)定性條件。通過(guò)數(shù)值算例，取不同的時(shí)間步長(zhǎng)來(lái)驗(yàn)證對(duì)顯式精確算法的穩(wěn)定性分析。

數(shù)值模擬結(jié)果表明，當(dāng)時(shí)間步長(zhǎng)過(guò)大，不滿(mǎn)足本構(gòu)計(jì)算的穩(wěn)定性條件時(shí)，計(jì)算得到的等效屈服應(yīng)力出現(xiàn)了不穩(wěn)定現(xiàn)象；而當(dāng)時(shí)間步長(zhǎng)滿(mǎn)足穩(wěn)定性條件時(shí)，計(jì)算結(jié)果準(zhǔn)確地描述了材料的復(fù)雜變形過(guò)程為：彈性加載、塑性加載、彈性卸載、反向彈性加載和反向塑性加載。進(jìn)一步的數(shù)值模擬結(jié)果表明：推導(dǎo)得到的穩(wěn)定性條件可正確預(yù)測(cè)滿(mǎn)足本構(gòu)計(jì)算穩(wěn)定的最大時(shí)間步長(zhǎng)與各參數(shù)之間的關(guān)系。

[1] Drucker D C. A more fundamental approach to plastic stress-strain relations[M]. Division of Applied Mathematics, Brown University, 1951.

[2] Drucker D C, Prager W, Greenberg H J. Extended limit design theorems for continuous media[J]. Quarterly of Applied Mathematics, 1952,9(4):381-389.

[3] Palmer A C, Maier G, Drucker D C. Convexity of yield surfaces and normality relations for unstable materials or structural elements[J]. Journal of Applied Mechanics, 1967,34(2):464-470.

[4] 李永池,唐之景,胡秀章.關(guān)于Drucker公設(shè)和塑性本構(gòu)關(guān)系的進(jìn)一步研究[J].中國(guó)科學(xué)技術(shù)大學(xué)學(xué)報(bào),1988,18(3):339-345.

Li Yongchi, Tang Zhijing, Hu Xiuzhang. Further study on the drucker postulate and plastic constitutive relations[J]. Journal of China University of Science and Technology, 1988,18(3):339-345.

[5] Hageman L J, Walsh J M. Help, a multi-material Eulerian program for compressible fluid and elastic-plastic flows in two space dimensions and time. Volume 1: AD0726459[R]. Systems Science and Software, La Jolla, California, 1971.

[6] Autodyn theory manual revision 4.3[R]. Century Dynamics Limited, Management Consulting Services, Horsham, United Kingdom, 2000.

[7] Hallquist J O. LS-DYNA theory manual[M]. Livermore: Livermore Software Technology Corporation, 2006.

[8] 李永池,譚福利,姚磊,等.含損傷材料的熱粘塑性本構(gòu)關(guān)系及其應(yīng)用[J].爆炸與沖擊,2004,24(4):289-298.

Li Yongchi, Tan Fuli, Yao Lei, et al. Thermo-viscoplastic constitutive relation of damaged materials with application[J]. Explosion and Shock Waves, 2004,24(4):289-298.

Abstract: At first, we analyzed the numerical stability of the explicit exact algorithm developed for the viscoplastic material, and then found that the explicit exact algorithm is not absolutely stable, deduced a necessary criterion that the time step should be kept below a certain value to guarantee the constitutive calculation stability. A series of numerical examples were presented to validate the reliability of the present stability analysis on the explicit exact algorithm. The results of the numerical examples show that the effective stress is unstable while the stability criterion for the constitutive calculation is not satisfied, but a complex deformation process including the elastic load, the plastic load, the elastic unload, the reverse elastic load and the reverse plastic load is accurately described while the stability criterion is satisfied. Further numerical results indicate that the stability criterion can accurately predict the relationships between the maximum time step and each parameter.

Keywords: constitutive relation; numerical stability; viscoplastic constitutive; explicit precise algorithm

(責(zé)任編輯 王易難)

Thenumericalstabilityoftheconstitutivecalculationonviscoplasticmaterials

Liu Mingtao1,2, Li Yongchi2, Hu Xiuzhang2, Zhang Jie2

(1.InstituteofFluidPhysics,ChinaAcademyofEngineeringPhysics,Mianyang621999,Sichuan,China; 2.DepartmentofModernMechanics,UniversityofScienceandTechnologyofChina,Hefei230027,Anhui,China)

中國(guó)工程物理研究院流體物理研究所發(fā)展基金項(xiàng)目(SFZ201401(04)02)

O345國(guó)標(biāo)學(xué)科代碼1301520

A

10.11883/1001-1455(2017)05-0969-07

2015-06-29；

2015-10-08

國(guó)家自然科學(xué)基金項(xiàng)目(11602250)；

劉明濤(1986— )，男，博士，副研究員，liumingtao@caep.cn。