刘明涛，李永池，胡秀章，章 杰

(1.中国工程物理研究院流体物理研究所,四川 绵阳 621999; 2.中国科学技术大学近代力学系,安徽 合肥 230027)

黏塑性本构计算的稳定性分析*

刘明涛1,2，李永池2，胡秀章2，章 杰2

(1.中国工程物理研究院流体物理研究所,四川 绵阳 621999; 2.中国科学技术大学近代力学系,安徽 合肥 230027)

提出本构方程计算方法的稳定性问题，针对黏塑性本构计算的显式精确算法的稳定性进行分析，发现该算法并非无条件稳定，使用小扰动方法给出了其计算稳定的必要条件，稳定性条件对数值计算中的时间步长提出限制要求。通过有限元算例验证了分析的正确性，计算结果也表明理论推导得到的稳定性公式能够准确预测满足计算稳定性条件要求的最大时间步长与各参数之间关系。

本构关系；算法稳定性；黏塑性本构；显式精确算法

一直以来，塑性和黏塑性材料的本构方程及其数值算法的研究是计算力学的核心问题。D.C.Drucker等[1-2]提出了著名的Drucker公设，该公设是经典塑性理论的基石。利用Drucker公设可直接推导出塑性流动的正交法则和屈服面的外凸性，但其缺点是只适用于稳定材料。A.C.Palmer[3]指出针对软化的非稳定材料也可以得到塑性流动的正交法则和屈服面的外凸性。李永池等[4]进一步发展了A.C.Palmer的思想，提出了广义Drucker公设，将稳定材料和非稳定材料统一在一个框架之内。

Drucker公设推导出的塑性流动法则是塑性本构计算的基石，其指出塑性应变必须沿着当前加载面的法线方向发展。迄今为止，有关冲击动力学问题的计算程序中，本构算法大部分采用传统的半径回归法[5-7]，这是人们最初针对理想塑性材料提出的一种本构更新算法。对于具有应变率效应的黏塑性材料，采用这种算法会带来较大的误差，尤其是当材料由弹性状态进入屈服状态时。针对此问题，李永池等[8]提出了一种新的显式本构计算方法，称为显式精确算法。基于广义Drucker公设，李永池等[8]理论推导指出当黏塑性材料进入屈服状态后，利用材料的实时状态量即可唯一确定塑性流变过程中的塑性流动因子。若在一阶精度下，利用上时刻的状态量即可求出现时刻的塑性流动因子，继而可计算现时刻的塑性变形，而后利用胡克定律可得现时刻的应力增量。

通常情况下，数值计算的稳定性是针对连续方程、动量方程和能量方程的差分格式而言。本文中探讨本构方程计算方法的稳定性问题。对李永池等[8]发展的显式精确算法进行稳定性分析，并通过对单个单元的有限元算例验证其正确性。

1 显式精确算法及其计算稳定性

显式精确算法的详细理论推导过程见文献[8]，文献[8]中是从一般性的热黏塑性本构关系出发推导的，具有一定的普适性。本文中对其推导过程进行一定的简化，从目前常见的Mises屈服准则下的黏塑性本构关系出发推导。

1.1显式精确算法

针对Mises类黏塑性本构关系，屈服函数可设为：

(1)

其中，Mises等效应力和等效塑性应变率计算公式分别为：

(2)

(3)

根据广义Drucker公设可知，塑性流动的正交法则为：

(4)

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

(5)

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

(6)

由黏塑性本构方程式(1)，可反解出等效塑性应变率为：

(7)

式(6)、(7)联立，得：

(8)

由式(8)可知，塑性流动因子可由材料当前的应力状态唯一确定。

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

(9)

根据材料的胡克定律，有：

(10)

(11)

1.2计算稳定性分析

式(11)是一个张量表达式，不利于计算的稳定性分析。下面对此进行简化，选择一种特殊情况进行计算稳定性理论分析。

(12)

(13)

(14)

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

(15)

由上式可计算得：

(16)

(17)

(18)

将函数g(s(n)+Δs)在s(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)

计算稳定性要求扰动的放大因子小于1，即|1-3Gdtg′(s(n))|<1，解得：

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

(21)

式(21)即为推导得出的 显式精确算法的计算稳定性条件，显式精确算法的计算稳定性对计算时间步长dt提出了限制要求。还需特别指出的是，式(21)是在式(12)、(13)的特殊情况下推导得到的，并不能保证本构计算绝对稳定。

若选择黏塑性材料的具体屈服准则形式为：

(22)

将式(22)代入式(21)，可得本构方程的计算稳定性要求为：

(23)

传统Courant稳定性条件为：

(24)

式中：L为单元特征长度，C为材料绝热声速，α为安全因数，通常取α=0.9。式(21)、(23)表明，本构计算的稳定性与材料的本构参数和实时塑性应变率密切相关，而传统Courant稳定性条件仅与材料的声速和单元的尺寸相关。

2 数值计算算例

四边形单元的运动过程设定为：在整个变形过程中4个节点在x方向均固定不动，同时节点1、2在y方向也固定不动，首先，节点3和4以速度-50 m/s沿着y方向匀速运动20 μs，四边形单元的y方向长度由20 mm压缩至19 mm；而后，节点3和节点4再以50 m/s匀速运动返回至初始位置，四边形单元的y方向长度由19 mm回复至初始时的20 mm。

四边形单元在上述变形过程中经历的加卸载路径较复杂，如图1所示。从状态1到状态2，经历了弹性加载、塑性加载；从状态2到状态3，经历了弹性卸载、反向弹性加载、反向塑性加载，共5个阶段。在第2阶段和第5阶段，材料发生了塑性应变，塑性应变率均约为1 700 s-1。根据显式精确算法的稳定性条件式(2)，可知时间步长dt需满足：0

2.1不同时间步长算例

2.1.1时间步长dt=0.3×10-7s时

计算结果如图2所示，黑色实线为等效应力历史曲线，蓝色点划线为等效塑性应变历史曲线。从图中可以看出，当选取的时间步长满足计算稳定性条件式(23)时，利用显式精确算法计算的结果能够准确描述四边形单元在整个变形过程中所经历的复杂加卸载过程。

(1)弹性加载阶段：随着变形的增加，等效应力匀速增加，没有发生塑性变形。

(2)塑性加载阶段：由于节点3和节点4以匀速运动，因此在整个加载过程中等效塑性应变率保持不变，因此根据本构关系式(22)，等效应力也保持不变，而图2中等效应力在该阶段呈现1个平台段，计算结果正确。

(3)弹性卸载阶段： 20 μs以后，节点3和节点4沿原路径匀速返回，材料进入弹性卸载阶段，该过程中不产生塑性应变增量，因此等效塑性应变呈现平台段，等效应力由材料的屈服点匀速降低至零，计算结果正确。

(4)反向弹性加载阶段：当等效应力降低为零后，重新开始匀速增大，材料进入反向弹性加载阶段，在此阶段等效塑性应变不增大，呈现平台段，计算结果正确。

(5)反向塑性加载：当材料由于反向加载再次进入屈服后，塑性变形重新开始累积，由于在该阶段节点3和节点4匀速运动，因此塑性应变率保持为恒定值，所以塑性应变线性增大，等效应力出现第2个平台段，计算结果正确。

2.1.2时间步长dt=1.0×10-7s时

该时间步长不满足显式精确算法计算得出的稳定性条件，计算结果如图3所示。可以看出，当材料屈服后，计算得到的等效应力出现了上下抖动现象，此时显式精确算法的计算不稳定，计算结果错误。

2.1.3时间步长dt=1.5×10-7s时

进一步增大了时间步长，其远不满足本构计算的稳定性条件。计算结果如图4所示。从图4可以看出，当材料屈服后，等效塑性应变历史曲线呈台阶式上升，四边形单元在弹性与塑性状态之间来回跳动，等效应力也变得极不稳定，计算结果与真实值相差巨大。但是当采用传统的近似算法(半径回归法)并仍取时间步长为dt= 1.5×10-7s时，数值计算结果稳定收敛,具体结果如图5所示。

综合上述算例可以看出，显式精确算法的确存在计算稳定性问题，取同样的时间步长dt=1.5×10-7s，当采用传统的近似算法(半径回归法)得到了正确的结果，而采用显式精确算法时结果失稳。当时间步长减小至满足本构计算稳定性条件式(21)时，显式精确算法给出的计算结果也稳定收敛。

2.2最大时间步长与材料参数和塑性应变率关系

进一步分析满足本构计算稳定的最大时间步长与材料参数和塑性应变率关系，分别研究各个参数与最大时间步长的关系。分为4组计算，每组只变化1个变量，其余参数的值取表1中的参数值，数值模拟得到满足本构计算稳定的最大时间步长，其随各参数的变化如图6～9所示。

可以看出，本构计算稳定性准则式(21)、(23)与模拟结果符合的较好，最大时间步长与各参量的依赖关系为：与应变率敏感因子β成正比、与静态条件下屈服强度Y*成正比、与剪切模量G成反比、与塑性应变率成反比。但需要特别指出的是，本构计算稳定性准则式(21)只是显式精确算法计算稳定的必要性条件。在实际工程计算过程中，为提高其可靠性，可以取一个安全系数。

3 结 论

首先对最常用的Mises类黏塑性材料，重新推导显式精确算法的计算公式和流程，然后通过理论推导得到显式精确算法的稳定性条件。通过数值算例，取不同的时间步长来验证对显式精确算法的稳定性分析。

数值模拟结果表明，当时间步长过大，不满足本构计算的稳定性条件时，计算得到的等效屈服应力出现了不稳定现象；而当时间步长满足稳定性条件时，计算结果准确地描述了材料的复杂变形过程为：弹性加载、塑性加载、弹性卸载、反向弹性加载和反向塑性加载。进一步的数值模拟结果表明：推导得到的稳定性条件可正确预测满足本构计算稳定的最大时间步长与各参数之间的关系。

[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] 李永池,唐之景,胡秀章.关于Drucker公设和塑性本构关系的进一步研究[J].中国科学技术大学学报,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] 李永池,谭福利,姚磊,等.含损伤材料的热粘塑性本构关系及其应用[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

(责任编辑 王易难)

Thenumericalstabilityoftheconstitutivecalculationonviscoplasticmaterials

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

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

中国工程物理研究院流体物理研究所发展基金项目(SFZ201401(04)02)

O345国标学科代码1301520

A

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

2015-06-29；

2015-10-08

国家自然科学基金项目(11602250)；

刘明涛(1986— )，男，博士，副研究员，liumingtao@caep.cn。