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

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

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

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

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

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

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

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

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

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

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

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

1.1顯式精確算法

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

(1)

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

(2)

(3)

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

(4)

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

(5)

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

(6)

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

(7)

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

(8)

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

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

(9)

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

(10)

(11)

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

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

(12)

(13)

(14)

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

(15)

由上式可計算得：

(16)

(17)

(18)

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

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)

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

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

(21)

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

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

(22)

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

(23)

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

(24)

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

2 數(shù)值計算算例

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

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

2.1不同時間步長算例

2.1.1時間步長dt=0.3×10-7s時

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

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

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

(3)彈性卸載階段： 20 μs以后，節(jié)點3和節(jié)點4沿原路徑勻速返回，材料進(jìn)入彈性卸載階段，該過程中不產(chǎn)生塑性應(yīng)變增量，因此等效塑性應(yīng)變呈現(xiàn)平臺段，等效應(yīng)力由材料的屈服點勻速降低至零，計算結(jié)果正確。

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

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

2.1.2時間步長dt=1.0×10-7s時

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

2.1.3時間步長dt=1.5×10-7s時

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

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

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

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

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

3 結(jié) 論

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

數(shù)值模擬結(jié)果表明，當(dāng)時間步長過大，不滿足本構(gòu)計算的穩(wěn)定性條件時，計算得到的等效屈服應(yīng)力出現(xiàn)了不穩(wěn)定現(xiàn)象；而當(dāng)時間步長滿足穩(wěn)定性條件時，計算結(jié)果準(zhǔn)確地描述了材料的復(fù)雜變形過程為：彈性加載、塑性加載、彈性卸載、反向彈性加載和反向塑性加載。進(jìn)一步的數(shù)值模擬結(jié)果表明：推導(dǎo)得到的穩(wěn)定性條件可正確預(yù)測滿足本構(gòu)計算穩(wěn)定的最大時間步長與各參數(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].中國科學(xué)技術(shù)大學(xué)學(xué)報,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)

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

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

A

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

2015-06-29；

2015-10-08

國家自然科學(xué)基金項目(11602250)；

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