航空动力学报
JOURNAL OF AEROSPACE POWER
1998年 第4期 No.4 1998科技期刊
粘塑性材料结构的有限元分析方法
杨晓光*　*　耿　瑞
北京航空航天大学
【摘要】　本文针对Walker粘塑性本构方程，通过发展一自适应变步长积分算法，用户子程序的方法把该本构方程结合进了国际著名的非线性有限元软件ABAQUS，并进行了算例验证。它包括了四种单元类型，可考虑材料变形的率相关，加载的循环、非等温、变温等情况。适用于发动机热端部件非线性结构分析。
　主题词：　粘塑性　本构方程　有限元
　分类号：　V231.9　O346
1　引　言
　　很多的研究及应用表明［1］：目前基于经典弹塑性和蠕变理论的非线性结构分析方法不足以正确设计、改进和分析航空航天及核动力设备关键热端零件的强度及疲劳寿命。其主要原因之一是结构承受各种复杂恶劣的循环加载和瞬态加载条件，这些条件使材料经受了与应力、温度和时间相关的并且是耦合的非线性行为。基于经典弹塑性及蠕变理论的本构方程，通常是把非弹性应变分为与时间相关的蠕变和与时间无关的塑性。但是大量的实验表明：在高温下其本质上是时间相关的，并且两者有交互作用。这表明非弹性变形是由一单一的机理在控制，宜用统一的方法，即把塑性及蠕变相联系起来的弹粘塑性本构方程来处理。
　　目前，这类统一型的本构理论在近20年得到了很大的发展，已有十多种方法。尤其是与有限元结构分析方法后，逐渐在实际的工程分析得到了应用。这些本构模型大都根据不可逆热力学及位错动力学的研究成果，恰当地选择单个或多个内变量及其演变方程用以描述非弹性变形的诸方面，如塑性、蠕变、松驰等。本文将根据Walker粘塑性本构理论，通过发展自适应显式积分方法，把该本构方程结合进通用的大型非线性有限元程序ABAQUS，并进行了该本构理论率相关，循环，热恢复及非等温等分析能力的验证。
2　Walker粘塑性本构理论
　　一般应变率可分为弹性应变率和非弹性应变率，并且弹性部分遵循Hooke定律。

式中：E为杨氏模量，ν为泊松比，δij为Kronecker符号。非弹性部分根据对材料硬化和恢复的不同处理将有不同的形式。
　　Walker型本构理论［2］是一种无屈服条件的理论。这就意味非弹性变形将在加载的所有阶段出现，虽然在弹性范围内它很小。其流动准则为：＝λ（Sij-Ωij）。式中λ是与加载历史有关的标量，Ωij称为背应力或“平衡应力”，它表示了运动硬化。各向同性硬化用标量K表示。非弹性应变率［2］：

这里：　　　　　　　　Sij＝σij-(1/3)σkkδij
背应力率：
其中：
式中：K1，K2，n1，n2，n3，n4，n5，n6，n7，m，材料常数，R为累积等效非弹性应变。
　　一般来说本构方程是温度相关的，也即方程材料常数是与温度相关的。这样为了遵循不可逆热力学定律，以考虑到非等温的情况，Ωij变为：

3　本构方程积分方法
　　目前统一型粘塑性本构方程总是以率形式给出，这就需要通过合适的积分方法把率形式的本构方程转化为增量形式的本构方程。通常采用的积分方法有显式的、隐式的和混合式的。这类率形式的统一型本构方程在数学上称之为“刚”的，也即自变量一个小的变化将引起因变量一个很大的变化。这通常是运动方程的高度非线性特点所决定的。因此在其积分方法上需特别注意。
　　Kumar［3］等人对这类方程的各种积分方法进行了系统比较研究。有时，一种相对比较简单的显式积分方法可能优于高阶的或隐式的积分方法。Cassenti［4］曾成功地将Euler积分方法应用于Walker方程。显式积分方法虽然是条件稳定的，但是在积分时恰当地控制积分步长，积分结果和效率是可以接受的。有时它优于高阶的Runge-Kutta方法。特别是它可以很好地与目前的有限元软件相结合。
　　下面给出一种自适应变步长的Euler前差显式积分方法（SAFE），该方法与细分步长技术相结合，将结合进ABAQUS程序。
　　如上所述的本构方程，通常可以用一阶的非线性偏微分方程表示，这样在某一给定时刻t时的应力、应变和其他内变量可以以一更紧凑的形式给出：dy/dt=F(y，t)。y可以是应力，非弹性应变，背应力等；F(y,t)是上述本构方程右端非线性项的缩写。这些偏微分方程通常是高度非线性的；积分时，为了保证其稳定性，需要较小的积分步长。当采用细分步长技术积分本构方程时，假定初始值y(t=0)=yA已知，这样从时刻t（此时各状态变量已知，即应力、非弹性应变、背应力已知）到t+△t时刻的积分可按下面的方法完成：首先积分步长△t细分为NSPLIT等份，再采用Euler显式积分公式：
yt+△t/NSPLIT=yt+(△t/NSPLIT)F(y,t)　　　　　(7)
上方程重复NSPLIT遍即可得到y在t+△t时的解。△t一般为有限元中整体时间载荷增量（步长）。采用上面Euler显式方法进行本构方程积分的成功与否，很重要的一点是要恰当的选择时间步长（△t/NSPLIT），时间步长不能太大以至破坏“时间步长内为y为常变化率”的假定。对粘塑性本构方程，这通常要求很小的时间步长。为此，将采用下面的方法来预报积分的成功与否，即对每一细分步长，采用应变度量的方法，确保积分成功（时间步长合适）。首先定义下面的应变度量（应变的改变）：ERROR=△R+，△J=(3/2)△Sij△Sij。经过数值积分试算：当ERROR不超过10-4，积分总是稳定的。因此，可以采用变步长的方法改进在时间步长△t内，在细分步长△t/NSPLIT上的积分。由（7）式可知，细分的时间步长取决于NSPLIT，可以有三种方法决定NSPLIT。(1)如果ERROR在指定的范围内，通常上限为10-4，下限为10-5，NSPLIT不需改变；(2)如果ERROR小于下限10-5，NSPLIT减半；(3)如果ERROR大于10-4，NSPLIT翻倍，并且在该细分时间步上的积分重新进行。
　　有限元时间载荷步长△t的确定，一般依靠ABAQUS中时间步长增长策略。通常，可以选择由ABAQUS提供的自动时间步长增长策略，也可人为的使用固定时间步长增长策略，或在ABAQUS提供的用户子程序上，通过变量PNEWDT直接控制。
4　与ABAQUS的结合
　　ABAQUS是目前国际上著名的非线性有限元分析软件，它提供了一个用户子程序UMAT，以帮助用户把自己的材料模型结合进ABAQUS程序。ABAQUS是采用Newton-Raphson方法求解非线性有限元平衡方程，对准静态加载情况，有：

　　　　　　　　　　　　　△uk+1=△uk+Ck　　　　　　　　　　　　　　　　　　　　(9)
式中：B是位移增量与应变增量转换矩阵，△ε=B△u。应力增量是应变增量的函数，因此也是位移增量的函数。方程（8）给出了第K次迭代对位移增量的修正。当应用Newton-Raphson迭代方法求解有限元平衡方程时，需给出方程（8）中Jacobian阵J=，也即为UMAT子程序中的DDSDDE矩阵，它定义为在t+△t时，第J个应变分量的微小增量对I个应力分量的扰动。通常DDSDDE矩阵是对称的，但是当采用细分技术时不能得到J的显式表达式，因此这里采用弹性刚度矩阵，这就象有限元方法中的“常刚度”法。ABAQUS将在每个单元的每个积分点调用UMAT。
　　在求解的过程中，需记录保存与求解过程相关的变量，如非弹性应变，背应力等，这可以储存在UMAT子程序的NSTATV数组中，并且可以定义储存的数量。对Walker本构方程，这里定义了14个与求解过程相关的变量，即6个非弹性应变，6个背应力，1个累积等效非弹性应变，1个NSPLIT变量。
　　当运行ABAQUS时，每当调用材料模式，都进入UMAT子程序。ABAQUS主程序在进入UMAT时，向UMAT传进应变增量△εI，温度增量△θ，时间步长△t及外载增量△Pi，同时也传入当前已知状态的应力，应变，温度及其它与求解过程相关的变量。在UMAT中，用户根据本构方程，需更新应力及其它相关的变量，并提供DDSDDE给ABAQUS主程序。对所有加载步（STEP）的第一个载荷增量，ABAQUS都向UMAT传入零应变增量。如果该步内考虑循环加载，则对每一步内的所有增量，ABAQUS都传入零应变增量。 
　　考虑到应用于航空发动机热端部件，在对UMAT编程时，主要包括了：平面应力、平面应变、轴对称和3D实体元。它可以通过ABAQUS的NDI和NSHEAR变量进行识别。
5　算例验证
　　为了验证模型及所编的UMAT程序，采用了一个单元(8节点三维元或轴对称元)的计算模型。ABAQUS为5.6版，计算机为SUN Ultra I，内存64 M。材料为Hastelloy-X，其Walker理论材料数据来自文献［2］，见表1。计算的结果见图1～图4。比较图1和图2，表明程序

图1　649℃，应变率为3.7×10-5，1.2×10-6循环的应力应变曲线(a)3.7×10-5(b)1.2×10-6
表1　Hastelloy-X Walker理论材料数据

温度／℃427.537.648.760.871.962.
杨氏模量／Psi26.0e624.0e624.0e622.6e618.6e613.2e6
泊松比ν0.3220.3280.3340.3390.3450.351
K1／Psi50931.75631.95631.251886.91505.59292.
K2／Psi000000
m0.0590.0590.0790.2440.1950.233
1/n1.161.161.161.161.161.16
n1／Psi000000
n2／Psi30.e76.e71.5e72.0e75.0e61.0e6
n38000.1000.781.1179.673.312.
n4000000
n5000000
n6／（Psi1-m*Sec-1）00008.89e-42.73e-3
n7000000
／Psi00-2000.-2000.-1434.-1200.

具有预测在同温度下，循环变形的率相关性的能力。由图1(a)和图2(a)，图1(b)和图2(b)，可看出程序可预测在同应变率下，循环变形对温度的相关性。图3表明了程序预测变加载的的情况。特别是由图4，程序可以很好的预测变温时的材料变形响应，并给出了考虑温度历史相关

图2　871℃，应变率为3.7×10-5，1.2×10-6循环的应力应变曲线(a)3.7×10-5(b)1.2×10-6

图　3760℃时，变应变率加载的应力应变曲线　　　图4　应变率为3.7×10-5，变温条件下的　　　（从3.0×10-6到 3.0×10-5到 3.0×10-4）　　　应力应变曲线（从650℃到850℃到650℃）
项（dT/dt）与不考虑温度历史相关项时，计算结果的差别。与文献［2］给出的实验结果比，计算结果的准确性是满意的，积分算法是成功的，计算效率是可以接收的。目前由于只进行了单一单元的模型验证，对于复杂工程结构，本构理论及所编程序的能力还需进一步验证。
参　考　文　献
1　Chan K S, et al.A Survey of Unified Constitutive Theories.N85-31531,1985
2　Walker K P.Research and Development Program for Non-Linear Structural Modeling with Advanced Time-Temperature Dependent Constitutive Relationships.NASA CR-165553,1981
3　Kumar V,et al.Numerical Integration of Some Stiff Constitutive Models of Inelastic Deformation. Trans.ASME J. of Eng. Mat. and Tech., 1980,102:91-96
4　Cassenti B N.Research and Development Program for the Development of Advauced Time-Temperature Dependent Constitutive Relationships.NASA CR-168191,1983,2:
（责任编辑　杨再荣）
1998年5月收稿；1998年6月收到修改稿。
*　*男　36岁　教授　北京航空航天大学405教研室　100083
