4.1 材料非线性概述

许多与材料有关的参数可以使结构刚度在分析期间改变。塑性、非线性弹性、超弹性材料、混凝土材料的非线性应力—应变关系，可以使结构刚度在不同载荷水平下(以及在不同温度下)改变。蠕变、粘塑性和粘弹性可以引起与时间、率、温度和应力相关的非线性。膨胀可以引起作为温度、时间、中子流水平(或其他类似量)函数的应变。

ANSYS程序应可以考虑多种材料非线性特性：

1．率不相关塑性指材料中产生的不可恢复的即时应变。

2．率相关塑性也可称之为粘塑性，材料的塑性应变大小将是加载速度与时间的函数。 3．材料的蠕变行为也是率相关的，产生随时间变化的不可恢复应变，但蠕变的时间尺度要比率相关塑性大的多。

4．非线性弹性允许材料的非线性应力应变关系，但应变是可以恢复的。

5．超弹性材料应力应变关系由一个应变能密度势函数定义，用于模拟橡胶、泡沫类材料，变形是可以恢复的。

6．粘弹性是一种率相关的材料特性，这种材料应变中包含了弹性应变和粘性应变。 7．混凝土材料具有模拟断裂和压碎的能力。

8．膨胀是指材料在中子流作用下的体积扩大效应。 4.2 塑性分析

4.2.1 塑性理论简介

许多常用的工程材料，在应力水平低于比例极限时，应力—应变关系为线性的。超过这一极限后，应力—应变关系变成非线性，但却不一定是非弹性的。以不可恢复的应变为特征的塑性，则在应力超过屈服点后开始出现。由于屈服极限与比例极限相差很小，ANSYS程序在塑性分析中，假设这二个点相同，见 图4-1 。

图4-1 弹塑性应力-应变曲线

塑性是一种非保守的(不可逆的)，与路径相关的现象。换句话说，荷载施加的顺序，以及什么时候发生塑性响应，影响最终求解结果。如果用户预计在分析中会出现塑性响应，则应把荷载处理成一系列的小增量荷载步或时间步，以使模型尽可能附合荷载—响应路径。最大塑性应变是在输出(Jobname.OUT)文件的子步信息中打印的。

在一个子步中，如果执行了大量的平衡迭代，或得到大于15%的塑性应变增量，则塑性将激活自动时间步选项[ AUTOTS ](GUI：Main Menu>Solution> Sol\或 Main Menu>Solution>Unabridged Menu> Time /Frequenc>Time and Substps)。如果取了太大的时间步，则程序将二分时间步，并重新求解。

其他类型的非线性行为可以与塑性同时产生。实际上，大位移和大应变几何非线性经常伴随有塑性材料响应，如果用户预计在结构中存在大变形，则必须在分析中用 NLGEOM 命令激活这些效应(GUI:Main Menu>Solution>Sol\

Menu>Solution>Unabridged Menu>Analysis)。对于大应变分析，材料应力—应变特性必须按真实应力和对数应变输入。

在这一节中，我们将依次介绍塑性理论的三个主要方面：

· 屈服准则 · 流动准则 · 强化准则 4.2.1.1 屈服准则

对单向受拉试件，我们可以通过简单的比较轴向应力与材料的屈服应力来决定是否有塑性变形发生，然而，对于一般的应力状态，是否到达屈服点并不是明显的。

屈服准则是一个可以用来与单轴测试的屈服应力相比较的应力状态的标量表示。因此，知道了应力状态和屈服准则，程序就能确定是否有塑性应变产生。

在多轴应力状态下，屈服准则可以用下式来表示：

σ e =f({σ})=σ y

其中σ e 为等效应力，σ y 为屈服应力。

当等效应力超过材料的屈服应力时，将会发生塑性变形。 1． Von Mises屈服准则

Von Mises屈服准则是一个比较通用的屈服准则，尤其适用于金属材料。对于Von Mises屈服准则，其等效应力为：

其中σ 1 ，σ 2 ，σ 3 为三个主应力。

可以在主应力空间中画出Mises屈服准则，见 图 4-2 。

图 4-2 主应力空间中的Mises屈服面

在3D主应力空间中，Mises屈服面是一个以σ 1 =σ 2 =σ 3 为轴的圆柱面，在2D中，屈服面是一个椭圆，在屈服面内部的任何应力状态，都是弹性的，屈服面外部的任何应力状态都会引起屈服。注意：静水压应力状态（σ 1 =σ 2 =σ 3 ）不会导致屈服：屈服与静水压应力无关，而只与偏差应力有关，因此，σ 1 =180，σ 2 =σ 3 =0的应力状态比σ 1 =σ 2 =σ 3 =180的应力状态接近屈服。Mises屈服准则是一种除了土壤和脆性材料外典型使用的屈服准则，在土壤和脆性材料中，屈服应力是与静水压应力（侧限压力）有关的，侧限压力越高，发生屈服所需要的剪应力越大。 2． Hill屈服准则

以上介绍的Von Mises屈服准则是一个各向同性的屈服准则，而Hill屈服准则是各向异性的，可以考虑材料的弹性参数的各向异性和屈服强度的各向异性。它是Von Mises屈服准则的延伸。Hill屈服准则的等效应力被表达为：

其中六个材料常数 H , F , G , N , L , M 可以由试验结果取得，如下各式所示：

， ，

， ，

以上各式右端的六个常书为材料的屈服应力比率，如下：

， ，

，

其中

，

为各个方向的屈服应力，σ 0 为参考屈服应力。

Hill屈服面在主应力空间中如 图 4-3 。

图 4-3 主应力空间中的Hill屈服面

在3D主应力空间中，Hill屈服面是一个以σ 1 =σ 2 =σ 3 为轴的椭圆柱面，在2D中，屈服面也是一个椭圆，在屈服面内部的任何应力状态，都是弹性的，屈服面外部的任何应力状态都会引起屈服。

3． 广义Hill屈服准则

广义Hill屈服准则是对Hill屈服准则的进一步延伸，广义希尔准则不仅考虑了在材料的三个正交方向屈服强度的不同，而且还考虑了拉伸状态和压缩状态下屈服强度的不同。 广义Hill屈服准则的等效应力可以表示为：

其中：

,j=1 to 6

j=1 to 3

以上各式中，σ +j 和σ -j 分别是j方向的拉伸和压缩屈服强度，在此，压缩屈服应力被作为正值处理。对剪切屈服，σ +j =σ -j 。令M 11 =1，则有：

K=σ +x σ -x

单轴屈服强度值的选取需使得以下两式成立

广义Hill屈服准则使用等向强化准则，因此在应力空间中，屈服面是一个经过移动的椭圆柱面，其大小随塑性应变而胀缩，见 图4-4 。

图 4-4 主应力空间中的广义Hill屈服面

4． Drucker-Prager(DP)屈服准则

Drucker-prager屈服准则是对Mohr-Coulomb准则的近似，它修正了Von Mises 屈服准则，即在Von Mises表达式中包含一个附加项。其屈服面并不随着材料的逐渐屈服而改变，因此

没有强化准则, 塑性行为被假定为理想弹塑性，然而其屈服强度随着侧限压力(静水应力)的增加而相应增加，另外，这种材料考虑了由于屈服而引起的体积膨胀，但不考虑温度变化的影响。此材料选项适用于混凝土、岩石和土壤等颗粒状材料。 对DP材料，其屈服准则表达式为：

其中：

材料常数β和屈服强度σ y 的表达式如下：

，{S}偏差应力

其中：φ为内摩擦角，C为粘滞力。

对DP材料，当材料参数β,σ y 给定后，其屈服面为一圆锥面，此圆锥面是六角形的摩尔－库仑屈服面的外接锥面，如 图4-5 所示：

图4-5 Drucker-Prager 屈服面

4.2.1.2 流动准则

流动准则描述了发生屈服时，塑性应变的方向，也就是说，流动准则定义了单个塑性应变分量（

等）随着屈服是怎样发展的。流动准则由以下方程给出：

其中：

是塑性乘子（决定了塑性应变量），

为塑性势，是应力的函数（决定了塑

性应变方向）。

一般来说，流动方程是塑性应变在垂直于屈服面的方向发展的屈服准则中推导出来的，即

等于屈服函数，这种流动准则叫作关联流动准则，如果使用其它的流动准 则（从其它

不同的函数推导出来），则叫作不关联的流动准则。 4.2.1.3 强化准则

强化准则描述了初始屈服准则随着塑性应变的增加是怎样发展的。

一般来说，屈服面的变化是以前应变历史的函数，在ANSYS程序中，使用了三种强化准则:

等向强化 是指屈服面以材料中所作塑性功的大小为基础在尺寸上扩张。对Mises屈服准则来说，屈服面在所有方向均匀扩张。见 图 4 - 6 。

图4-6 等向强化时的屈服面变化图

由于等向强化，在受压方向的屈服应力等于受拉过程中所达到的最高应力。

随动强化 假定屈服面的大小保持不变而仅在屈服的方向上移动，当某个方向的屈服应力升高时，其相反方向的屈服应力应该降低。见 图 4-7 。

图4-7 随动强化时的屈服面变化图

在随动强化中，由于拉伸方向屈服应力的增加导致压缩方向屈服应力的降低，所以在对应的两个屈服应力之间总存一个

的差值，初始各向同性的材料在屈服后将不再是各向

同性的。

混合强化是等向强化和随动强化的结合，屈服面不仅在大小上扩张，而且还在屈服的方向上移动。见 图 4-8 。

图4-8 混合强化时的屈服面变化图

4.2.2 塑性材料选项

一些选项可用于描述塑性行为。用户也可以应用Ansys可编程特性在程序中加入其他特性。参见《ANSYS Guide to User Programmable Features》。 1 、双线性随动强化(BKIN)选项

该选项假设总应力范围等于屈服应力的二倍，以包括包辛格效应(见 图4-10 )。对于服从Von Misses屈服准则的一般小应变分析，建议用这一选项。不应该应用于大应变分析。用户可以用BKIN和HILL选项组合来模拟各向异性随动强化塑性。应力-应变-温度数据的示例如下。 图4-9 说明了这种材料选项的典型显示[ TBPLOT ]。

MPTEMP,1,0,500 ! Define temperatures for Young\MP,EX,1,12E6,-8E3 ! C0 and C1 terms for Young\TB,BKIN,1,2 ! Activate a data table TBTEMP,0.0 ! Temperature = 0.0

TBDATA,1,44E3,1.2E6 ! Yield = 44,000; Tangent modulus = 1.2E6 TBTEMP,500 ! Temperature = 500

TBDATA,1,29.33E3,0.8E6 ! Yield = 29,330; Tangent modulus = 0.8E6 TBLIST,BKIN,1 ! List the data table

/XRANGE,0,0.01 ! X-axis of TBPLOT to extend from varepsilon=0 to 0.01 TBPLOT,BKIN,1 ! Display the data table

上面的命令 MPTEMP,MP,TB,TBTEMP,TBDATA,TBLIST,/XRANGE 和 TBPLOT 参见《ANSYS Commands Reference》。

图4-9 (a)双线性随动强化;(b)多线性随动强化

图4-10 包辛格效应

2 、多线性随动强化(KINH 和 MKIN)选项

这一选项应用Besseling模型，也称为子层或覆盖模型，以包括包辛格效应。KINH要比MKIN更好，因为前者允许用户定义更多的应力-应变曲线(40：5)，并且每条曲线上允许定义更多的点(20：5)。对于这二个选项，如果用户定义了多于一条应力-应变曲线(对于温度相关特性)，则每条曲线应包含相同数目的点。其假设是，在不同应力-应变曲线上的相应点，代表一个特定的子层的温度相关屈服行为。这两个选项不建议用于大应变分析。 图4-9 也说明了典型的应力-应变曲线(MKIN)。应用KINH选项的典型应力-应变温度数据如下：

TB,KINH,1,2,3 ! Activate a data table TBTEMP,20.0 ! Temperature = 20.0

TBPT,,0.001,1.0 ! Stress = 0.001, Strain = 1.0 TBPT,,0.1012,1.2 ! Stress = 0.1012, Strain = 1.2 TBPT,,0.2013,1.3 ! Stress = 0.2013, Strain = 1.3 TBTEMP,40.0 ! Temperature = 40.0

TBPT,,0.008,0.9 ! Stress = 0.008, Strain = 0.9 TBPT,,0.09088,1.0 ! Stress = 0.09088, Strain = 1.0 TBPT,,0.12926,1.05 ! Stress = 0.12926, Strain = 1.05

在上面这个例子中，一条曲线的第3点定义第3子层的温度相关屈服行为。 应用MKIN选项的典型应力-应变温度数据如下：

MPTEMP,1,0,500 ! Define temperature-dependent EX, MP,EX,1,12E6,-8E3 ! as in BKIN example TB,MKIN,1,2 ! Activate a data table

TBTEMP,,STRAIN ! Next TBDATA values are strains TBDATA,1,3.67E-3,5E-3,7E-3,10E-3,15E-3 ! Strains for all temps TBTEMP,0.0 ! Temperature = 0.0

TBDATA,1,44E3,50E3,55E3,60E3,65E3 ! Stresses at temperature = 0.0 TBTEMP,500 ! Temperature = 500

TBDATA,1,29.33E3,37E3,40.3E3,43.7E3,47E3 ! Stresses at temperature = 500 /XRANGE,0,0.02 TBPLOT,MKIN,1

有关的命令 MPTEMP,MP,TB,TBPT,TBTEMP,TBDATA,/XRANGE 和 TBPLOT ,参见《ANSYS Commands Reference》。

3 、非线性随动强化(CHABOCHE)选项

这一选项应用Chaboche模型，这种模型是多分量非线性随动强化模型，允许用户迭加几种随动强化模型。见《ANSYS Theory Reference》。象BKIN和MKIN选项一样，用户可应用CHABOCHE选项来模拟单调强化和包辛格效应。这个选项还允许用户模拟材料的棘轮和调整(Shakedown)效应。把CHABOCHE选项与各向同性硬化模型选项BISO、MISO、NLISO

组合起来，可以进一步模拟周期强化或软化。这种模型有1+2n个常数(其中n为随动强化模型数目)，用 TB 命令的 NPTS 定义。见《ANSYS Theory Reference》。用户应用 TBTEMP 和 TBDATA 命令定义材料常数。这种模型适合于大应变分析。 下面是温度不相关及一个随动强化模型的典型数据：

TB,CHABOCHE,1 ! Activate CHABOCHE data table TBDATA,1,C1,C2,C3 ! Values for constants C1, C2, and C3

而下例则说明温度相关常数及2个随动强化模型(在２个温度点)的典型数据表：

TB,CHABOCHE,1,2,2 ! Activate CHABOCHE data table TBTEMP,100 ! Define first temperature

TBDATA,1,C11,C12,C13,C14,C15 ! Values for constants C11, C12, C13, ! C14, and C15 at first temperature TBTEMP,200 ! Define second temperature

TBDATA,1,C21,C22,C23,C24,C25 ! Values for constants C21, C22, C23,

! C24, and C25 at second temperature

有关命令 TB,TBTEMP,TBDATA 参见《ANSYS Commands Reference》。 4 、双线性各向同性强化(BISO)选项

这一选项与多线性各向同性强化MISO选项相似，只是用双线性曲线代替多线性曲线。其输入类似于双线性随动强化选项，只是现在 TB 命令要应用BISO标号。这一选项通常对大应变分析较佳。用户可以把这一选项与非线性随动强化(CHABOUCHE)选项组合，以定义材料的各向同性强化行为。用户也可以把这一选项与HILL选项组合来模拟各向异性塑性及各向同性强化，或与RATE选项组合以模拟率相关粘塑性。

5 、多线性各向同性强化(MISO)选项

这一选项应用von Mises屈服准则以及各向同性工作强化的假定。这个选项不建议用于周期荷载或高度非比例历史荷载的小应变分析。但可应用于大应变分析。MISO选项可包括20个不同温度曲线，每条曲线可以有最多100个不同的应力-应变点。在各条曲线上，应变点可以不同。用户可以把这个选项与非线性随动强化(CHABOCHE)选项组合，以模拟周期强化或弱化。用户还可以把MISO选项与HILL选项组合来模拟各异性塑性及各向同性强化，以及与RATE选项组合用以模拟率相关粘塑性。在MKIN选项的示例中的应力-应变-温度曲线可以作为多线性各向同性强化材料的输入，如下：

MPTEMP,1,0,500 ! Define temperature-dependent EX, MP,EX,1,12E6,-8E3 ! as in above example TB,MISO,1,2,5 ! Activate a data table TBTEMP,0.0 ! Temperature = 0.0

TBPT,DEFI,3.67E-3,29.33E3 ! Strain, stress at temperature = 0 TBPT,DEFI,5E-3,50E3 TBPT,DEFI,7E-3,55E3 TBPT,DEFI,10E-3,60E3 TBPT,DEFI,15E-3,65E3

TBTEMP,500 ! Temperature = 500

TBPT,DEFI,3.67E-3,29.33E3 ! Strain, stress at temperature = 500 TBPT,DEFI,5E-3,37E3 TBPT,DEFI,7E-3,40.3E3 TBPT,DEFI,10E-3,43.7E3 TBPT,DEFI,15E-3,47E3

/XRANGE,0,0.02 TBPLOT,MISO,1 有关命令 MPTEMP、MP、TB、TBTEMP、TBPT、/XRANGE 和 TBPLOT 见《ANSYS Commands Reference》。

6 、非线性各向同性强化(NLISO)选项

这一选项基于 Voce 强化准则，见《ANSYS Theory Reference》。NLISO选项是MISO选项的一个变种，即指数饱和强化项扩展到线性项，见 图4-11 。这一选项的优点是材料行为定义为特殊函数，其中四个材料常数通过 TBDATA 命令来定义。用户可以通过把材料拉伸应力-应变曲线适当地试配得到材料常数。与MISO选项不同的是，不需要注意如何恰当地定义成对的材料应力-应变点。但是，这一选项仅适用于拉伸曲线与 图4-11 所示相同者。这一选项适用于大应变分析并且用户可以把这一选项与非线性随动强化(CHABCHE)选项组合，用于定义材料的各向同性强化行为。用户也可以把这一选项与HILL选项组合，用于模拟各向异性塑性及各向同性强化；或与RATE组合，用于模拟率相关粘塑性。下面例子说明在二个温度点上的温度相关数据表：

TB,NLISO,1 ! Activate NLISO data table TBTEMP,100 ! Define first temperature

TBDATA,1,C11,C12,C13,C14 ! Values for constants C11, C12, C13, ! C14 at first temperature

TBTEMP,200 ! Define second temperature

TBDATA,1,C21,C22,C23,C24 ! Values for constants C21, C22, C23, ! C24 at second temperature

有关命令 TB、TBTEMP、TBDATA 参见《ANSYS Commands Reference》。

图4-11 NLISO应力-应变曲线

7 、各向异性(ANISO)选项

这一选项允许在材料x、y、z方向上有不同的双线性应力-应变行为，以及在受拉、受压、受剪时有不同的行为。这一选项适用于预先受到变形的金属(如轧制)。不推荐用于周期荷载或非比例荷历程，因为假设了工作强化。屈服应力和斜率不完全无关，见《ANSYS Theory Reference》。为了定义各向异性材料特性，应用 MP 命令(GUI：Main Menu>

Solution>Other>Change Mat Props)来定义弹性模量(EX、EY、EZ、NUXY、NUYZ、NUXZ)。然后，应用 TB 命令[ TB ，ANISO]和 TBDATA 命令定义屈服点和正切模量。参见《ANSYS Elements Reference》中的非线性应力-应变材料。 8 、HILL各向异性(HILL)选项

这一选项与其他选项组合，可模拟塑性、粘塑性、蠕变(应用HILL模型)，见§4.6。 Hill各向异性仅适用于下列单元：PLANE42、SOLID45、PLANE82、SOLID92、SOLID95、

LIMK180、SHELL181, PLANE182、PLANE183、SOLID185、SOLID186、SOLID187、BEAM188和BEAM189。

下面例子说明HILL选项与BISO选项的组合：

TB,HILL,1,2 ! Activate HILL data table for two temps. TBTEMP,100 ! Define first temperature as 100

TBDATA,1,1,1.0402,1.24897,1.07895,1,1 ! Values for Hill constants C1 to C6 TBTEMP,200 ! Define second temperature as 200

TBDATA,1,0.9,0.94,1.124,0.97,0.9,0.9 ! Values for Hill constants C1 to C6 TB,BISO,1,2 ! Activate BISO data table for two temps. TBTEMP,100 ! Define first temperature as 100

TBDATA,1,461.0,374.586 ! Values for BISO constants C1 and C2 TBTEMP,200 ! Define second temperature as 200

TBDATA,1,461.0,374.586 ! Values for BISO constants C1 and C2

9 、Drucker-Prager(DP)选项 这一选项用于颗粒状(摩擦)材料，如土、岩体、砼等，并利用圆锥面来近似Mohr-Coulomb定律。

MP,EX,1,5000 MP,NUXY,1,0.27 TB,DP,1

TBDATA,1,2.9,32,0 ! Cohesion = 2.9 (use consistent units), ! Angle of internal friction = 32 degrees, ! Dilatancy angle = 0 degrees

有关命令 MP、TB、TBDATA 见《ANSYS Commands Reference》。 4.2.3 怎样使用塑性

在这一节中，我们将介绍在程序中怎样使用塑性，重点介绍以下几个方面：

· 可用的ANSYS输入。 · ANSYS输出量。 · 使用塑性的一些原则。 · 加强收敛性的方法。 · 查看塑性分析的结果。 4.2.3.1 ANSYS输入

当使用TB命令选择塑性选项和输入所需常数时，应该考虑到：

· 常数应该是塑性选项所期望的形式，例如，我们总是需要应力和总的应变，而不是应力与塑性应变。 · 如果还在进行大应变分析，应力－应变曲线数据应该是真实应力－真实应变。

对双线性选项（BKIN，BISO），输入常数

和

可以按下述方法来决定，如果材

可

料没有明显的屈服应力，通常以产生0.2%的塑性应变所对应的应力作为屈服应力，而

以通过在分析中所预期的应变范围内来拟合实验曲线得到。 其它有用的载荷步选项： · 使用的子步数（使用的时间步长）。既然塑性是一种与路径相关的非线性，因此需要使用许多载荷增量来加载。 · 激活自动时间步长。

· 如果在分析所经历的应变范围内，应力－应变曲线是光滑的，使用预测器选项，这能够极大地降低塑性分析中的总的迭代数。 4.2.3.2 输出量

在塑性分析中，对每个节点都可以输出下列量：

EPPL－塑性应变分量

,

等

EPEQ－累加的等效塑性应变

SEPL－根据输入的应力－应变曲线估算出的对于EPEQ的等效应力 HPRES－静水压应力 PSV－塑性状态变量

PLWK－单位体积内累加的塑性功

上面所列节点的塑性输出量实际上是指离节点最近的那个积分点的值。

如果一个单元的所有积分点都是弹性的（EPEQ＝0），那么节点的弹性应变和应力从积分点外插得到，如果任一积分点是塑性的（EPEQ>0），那么节点的弹性应变和应力实际上是积分点的值，这是程序的缺省情况，但可以人为的改变它。 4.2.3.3 程序使用中的一些基本原则

下面的这些原则应该有助于执行一个精确的塑性分析。

1． 需要的塑性材料常数必须能够足以描述所经历的应力或应变范围内的材料特性。 2． 慢加载，应该保证在一个时间步内，最大的塑性应变增量小于5%，一般来说，如

果Fy是系统刚开始屈服时的载荷，那么在塑性范围内的载荷增量应近似为： · 0.05*Fy－ 对用面力或集中力加载的情况 · Fy－ 对用位移加载的情况

3． 当模拟类似梁或壳的几何体时，必须有足够的网格密度，为了能够足够地模拟弯曲

响应，在厚度方向必须至少有二个单元。

4． 除非那个区域的单元足够大，应该避免应力奇异，由于建模而导致的应力奇异有： · 单点加载或单点约束 · 凹角 · 模型之间采用单点连接 · 单点耦合或接触条件 5． 如果模型的大部分区域都保持在弹性区内，那么可以采用下列方法来降低计算时间： · 在弹性区内仅仅使用线性材料特性（ 不使用TB命令） · 在线性部分使用子结构 4.2.3.4 加强收敛性的方法

如果不收敛是由于数值计算导致的，可以采用下述方法来加强问题的收敛性：

1．使用小的时间步长

2．如果自适应下降因子是关闭的，打开它，相反，如果它是打开的 ，且割线刚度正在被连续地使用，那么关闭它。

3． 用线性搜索，特别是当大变形或大应变被激活时。

4． 预测器选项有助于加速缓慢收敛的问题，但也可能使其它的问题变得不稳定。 5．可以将缺省的牛顿－拉普森选项转换成修正的（MODI）或初始刚度（INIT）牛顿－拉普森选项，这两个选项比全牛顿－拉普森选项更稳定（ 需要更多的迭代），但这两个选项仅在小位移和小应变塑性分析中有效。 4.2.3.5 查看结果

1．感兴趣的输出项（例如应力，变形，支反力等）对加载历史的响应应该是光滑的，一个不光滑的曲线可能表明使用了太大的时间步长或太粗的网 格。

2．时间步长内的塑性应变增量应该小于5％，这个值在输出文件中以“Max plastic Strain Step”输出，也可以使用POST26来显示这个值（Main Menu:Time Hist Postpro> Define Variables）。

3．塑性应变等值线应该是光滑的，通过任一单元的梯度不应该太大。

4．画出某点的应力—应变图，应力是指输出量SEQV（Mises 等 效 应 力），总应变由累加的塑性应变EPEQ和弹性应变得来。 4.2.4 塑性分析实例－DP材料实例分析 4.2.4.1 问题描述

一根铁桩插入土壤中，铁桩上端受到垂直向下的载荷，试分析此时铁桩对土壤的影响。由于土壤区域无限大，我们只取一相对于铁桩来说足够大的半圆形区域作为研究对象，其与外界土壤的联系通过弹簧单元来模拟。

图4-12 问题描述图

4.2.4.2 问题详细说明 材料特性： Ex=5000 (泊松比)＝0.4 C(凝聚力)=10

(内摩擦角)=30 (膨胀角)=30

4.2.4.3 求解步骤（GUI方法） 步骤一：建立计算所需要的模型

在这一步中，建立计算分析所需要的模型，包括定义单元类型，创建结点和单元，施加载何并将数据库保存为“dp1.db”，在此对这一过程不再详述（用户可以从4.2.4.4节的命令流执行到这一步骤）。

步骤二：恢复数据库文件 “ dp1.db ” Utility menu: file>Resume from 步骤三：定义材料性质 1、选“Main Menu>Preprocessor>Material Props>Material Models”。出现“Define Material Model Behavior”对话框，选择Material Model Number 1。

2、在“Material Models Available”窗口，双击“Structural->Linear->Elastic-> Isotropic”。出现一个对话框。

3、对杨氏模量（EX）键入5000, 对泊松比（NUXY）键入0.4。 4、单击OK。

步骤四：定义DP数据表并输入相应值

1、在“Material Models Available”窗口，双击“Structural->Nonlinear->Inelastic->Non-metal Plasticity->Drucker-Prager，出现输入框。

2、根据所给材料特性，在数据表中输入相应值。

3、单击OK。并退出“Define Material Model Behavior”对话框。 步骤五：进入求解器

选择菜单路径Main Menu>Solution 步骤六：加载

根据所给条件，施加适当的约束和载荷。在此不作详述，参考命令流文件。 步骤七：定义分析类型和分析选项：

1、选择菜单路径Main Menu>Solution>-Analysis Type-New Analysis. 2、单击“Static”来选中它然后单击OK。

3、选择菜单路径Main Menu>Solution>Unabridged Menu>Analysis options。对话框出现。 4、单击Large deform effects option（大变形效应选项），使之为ON，然后单击OK。 步骤八：设置输出控制选项

1、选择菜单路径：Main Menu > Solution>Unabridged Menu >Load step opts-Output ctrls > DB/Resuls File：对话框出现 2、在“Item”中，选择“all”

3、对“FREQ”,选择“Every Substep” 4、单击OK

步骤九：设置载荷步选项

1、选择菜单路径Main Menu> Solution>Unabridged Menu>Load step opts-Time/Frequenc> Time and substps。对话框出现。

1、对time at end of load step（载荷步终止时间）键入1 2、对“DELTIM”(Time step size) 输入0.1

3、将“AUTOTS”(automatic time stepping)设置为ON 4、对 Minimum time step size 输入0.05 5、对 Maximum time step of size 输入0.2 步骤十：打开时间步长预测器 1、选择菜单路径Main menu> Solution>Unabridged Menu>Load step opts-Nonlinear> Predictor。 2、将predictor的状态设置为“ON”。 步骤十一：打开线性搜索

1、选择菜单路径Solution>Unabridged Menu>Load step opts-Nonlinear> Line Search。 2、将LNSRCH的状态设置为“ON”。 步骤十二：设置收敛性控制

1、选择菜单路径Main menu> Solution>Unabridged Menu>Load step opts-Nonlinear> Convergence Crit，对话框出现。 2、击Add,下一级对话框出现。 3、对“TOLER”，输入0.2 4、对“MINREF”，输入200 5、单击OK。

步骤十三：进行求解 步骤十四：进行后处理

4.2.4.4 求解步骤（命令流方法）

fini /cle l=60 l1=15 w=10 ri=50 ro=160 routn=200 pp=600 pp1=30 pex=5000 /prep7 et,1,42 et,2,14

keyopt,1,3,2 keyopt,2,3,2

r,1,3.1416*pex/19 r,2,3.1416*pex/19 mp,ex,2,2e5 mp,nuxy,2,0.3

rect,-w/2,w/2,-l1,l-l1 pcirc,,ro,-180,0 pcirc,,ri,-180,0 aovl,all

lsel,s,loc,y,-ro,-ri lesize,all,,,18 alls

lsel,s,loc,x,-ro,-ri lsel,a,loc,x,ri,ro lesize,all,,,10 alls

lsel,s,loc,y,(l-l1)/2 lesize,all,,,8 alls

lsel,s,loc,y,0

lsel,r,loc,x,-w/2,w/2 lcom,all alls

lsel,s,loc,x,0

lsel,r,loc,y,-l1,l-l1 lesize,all,,,6 alls

lsel,s,loc,y,-l1/2 lesize,all,,,6 alls

lsel,s,loc,x,-1*(ri+w/2)/2 lsel,a,loc,x,(ri+w/2)/2 lesize,all,,,16 alls

mshape,0,2d mshkey,2

asel,s,loc,y,-ro,-ri amesh,all alls

asel,s,loc,y,(l-l1)/2 aatt,2 amesh,all alls

asel,s,loc,y,-l1/2 aatt,2 amesh,all alls

lsel,s,loc,y,-l1/2 lsel,a,loc,y,-l1 lccat,all alls

asel,s,loc,y,-ri,-l1 amesh,all alls csys,1

n,1001,routn,-180 ngen,19,1,1001,,,,10 type,2 real,1 e,1,1001 *do,i,3,19 e,i,1000+i-1 *enddo e,2,1019

csys,1

nsel,s,loc,x,ro

ngen,2,1100,all,,,ro*(2*1.732-3)/3,30 type,2 real,2 *do,i,1,19 e,i,1100+i *enddo alls

save,dp1,db resume,dp1,db mp,ex,1,5000 mp,nuxy,1,0.4 tb,dp,1

tbdata,1,10,30,30 fini /solu

nsel,s,loc,x,routn d,all,all

nsel,s,loc,x,ro*(2*1.732)/3 d,all,all csys,0

nsel,s,loc,y,l-l1 sf,all,pres,pp alls time,1

deltime,0.1,0.05,0.2 autot,on

cnvtol,f,,0.2,,200 nlgeom,on pred,on lnsrch,on outpr,all,all outres,all,all solve

4.3 超弹性分析 4.3.1 超弹理论 4.3.1.1 超弹的定义

一般工程材料（例如金属）的应力状态由一条弹塑性响应曲线来描述，而超弹性材料存在一个弹性势能函数，该函数是一个应变或变形张量的标量函数，而该标量函数对应变分量的导数就是相应的应力分量。

上式中：[S]＝第二皮奥拉－克希霍夫应力张量 W＝单位体积的应变能函数 [E]＝拉格朗日应变张量

拉格朗日应变可以由下式表达：[E]＝1/2（[C]-I） 其中：[I]是单位矩阵，[C]是有柯西－格林应变张量

其中[F]是变形梯度张量，其表达式为：

x ：变形后的节点位置矢量 X ：初始的节点位置矢量

如果使用主拉伸方向作为变形梯度张量和柯西－格林变形张量的方向，则有：

其中： J=初始位置与最后位置的体积比 材料在第i个方向的拉伸率

在ANSYS程序中，我们假定超弹材料是各向同性的，在每个方向都有完全相同的材料特性，在这种情况下，我们既可以根据应变不变量写出应变能密度函数，也可以根据主拉伸率写出应变能密度函数。

应变不变量是一种与坐标系无关的应变表示法。使用它们就意味着材料被假定是各向同性的。Mooney－Rivlin和Blatz－Ko应变能密度函数都可以用应变不变量表示，应变不变量可以柯西－格林应变张量和主拉伸率表示出来：

一个根据应量不变量写出来的应变能密度函数如下：

为材料常数，上式是两个常数的Mooney－Rivlin应变能密度函数。

超弹材料可以承受十分大的弹性变形，百分之几百的应变是很普遍的，既然是纯弹性应变，因此超弹性材料的变形是保守行为，与加载路径无关。 4.3.1.2 不可压缩缩性

大多数超弹材料，特别是橡胶和橡胶类材料，都是几乎不可压缩的，泊松比接近于0.5，不可压缩材料在静水压力下不产生变形，几乎不可压缩材料的泊松比一般在0.48至0.5之间（不包含0.5），对这些材料，在单元公式中必须考虑不可压缩条件。在ANSYS程序中，不可压缩超弹单元修改了应变能密度函数，在单元中明确地包含了压力自由度。压力自由度

使不可压缩条件得到满足，而不降低求解速度。压力自由度是一种内部自由度，被凝聚在单元内部。

4.3.1.3 超弹单元

有三种单元适合于模拟超弹性材料：

不可压缩单元有HYPE56，58，74和158，这些单元适用于模拟橡胶材料。

可压缩单元有HYPER84和86，HYPER84既可以是4节点矩形也可以是8节点矩形单元，这种单元主要用来模拟泡沫材料。

18X族单元(除LIMK和BEAM单元外，包括SHELL181,

PLANE182,PLANE183,SOLID185,SOLID186,和 SOLID187)。18X族单元消除了体积锁定， 既适用于不可压材料，又适用于可压材料。参见《ANSYS Elements Reference》的“Mixed U-P Formulations”。 4.3.2 超弹材料选项

超弹性可用于分析橡胶类材料(elastomers)，这种材料可承受大应变和大位移，但体积改变极微(不可压缩)。这种分析需用到大应变理论[ NLGEOM ,ON]。 图4-13 是一个例子。

图4-13 超弹性结构

在ANSYS超弹性模型中，材料响应总是假设各向同性和等温性。由于这一假设，应变能势函数按应变不变量来表示。除非明确指出，超弹性材料还假设为几乎或完全不可压缩材料。材料热膨胀也假定为各向同性的。

ANSYS在模拟不可压缩或几乎不可压缩超弹性材料时，应变能势函数有几种选项。这些选项均适用于SHELL181,PLANE182, PLANE183, SOLID185, SOLID186, SOLID187 单元。可以通过 TB ,HYPER 命令的 TBOPT 参数进入这些选项。

其中一个选项，Mooney-Rivhlin 选项，也适用于 HYPER56, HYPER58, HYPER74,

HYPER158 单元，以及显式动力分析单元 PLANE162,SHELL163, SOLID164。这一选项可通过 TB ,MOONEY 命令进入。

4.3.2.1 Mooney-Rivlin超弹性选项(TB,HYPER)

请注意本小节论述应用Mooney-Rivlin 选项与单元 SHELL181, PLANE182, PLANE183, SOLID185, SOLID186, SOLID187。

如果你想应用Mooney-Rivlin选项于单元HYPER56, HYPER58, HYPER74, HYPER158, PLANE162, SHELL163,SOLID164，则参见§4.3.2.7。

Mooney-Rivlin选项( TB ,HYPER,,,,MOOMEY)是缺省项，允许用户通过 TB 命令的 NPTS 参数定义2，3，5或9个参数。例如，为了定义5参数模型，用户采用 TB ,HYPER,1,,5,MOONEY。

2参数Mooney-Rivlin选项，适用于应变大约为100%(拉)和30%(压)的情况。与其他选项相比，较高阶的Mooney-Rivlin选项，对于较大应变的求解，可得到较好的近似。

下例是3参数Mooney-Rivlin选项的输入实例：

TB,HYPER,1,,3,MOONEY !Activate 3 parameter Mooney-Rivlin data table TBDATA,1,0.163498 !Define c10 TBDATA,2,0.125076 !Define c01 TBDATA,3,0.014719 !Define c11

TBDATA,4,6.93063E-5 !Define incompressibility parameter !(as 2/K, K is the bulk modulus)

对于本选项所需要的材料常数的描述，见《ANSYS Elements Reference》。 4.3.2.2 Ogden选项

Ogden选项( TB ,HYPER,,,,OGDEN)允许用户通过 TB 命令的 NPTS 参数定义无限参数，例如，应用 TB ,HYPER,1,,3,OGDEN 定义3参数模型。

与其他选项相比，Ogden选项通常对大应变水平的求解提供最好的近似。可应用的应变水平可达到700%。较高阶的参数可提供更精确的解。但是这样也可能在拟合材料常数时引起数值困难，而且它要求在用户感兴趣的变形范围内要有足够的数据。 下面是2参数Ogden选项的输入列表:

TB,HYPER,1,,2,OGDEN !Activate 2 parameter Ogden data table TBDATA,1,0.326996 !Define μ1 TBDATA,2,2 !Define α1 TBDATA,3,-0.250152 !Define μ2 TBDATA,4,-2 !Define α2

TBDATA,5,6.93063E-5 !Define incompressibility parameter !(as 2/K, K is the bulk modulus)

!(Second incompressibility parameter d2 is zero)

对于这个选项所需要的材料常数的论述，请参见《ANSYS Elements Reference》。 4.3.2.3 Neo-Hookean超弹性选项

Neo-Hookean选项( TB ,HYPER,,,,NEO)代表应变能势能的最简单形式，可用于应变范围20-30%。下面是Neo-Hookean选项的一个输入列表示例：

TB,HYPER,1,,,NEO !Activate Neo-Hookean data table TBDATA,1,0.577148 !Define incompressibility parameter !(as 2/K, K is the bulk modulus)

对于这个选项所需要的材料常数的论述，请参见《ANSYS Elements Reference》。 4.3.2.4 多项式超弹性选项

多项式选项( TB ,HYPER,,,,POLY)允许用户通过 TB 命令的 NPTS 参数定义无限多个参数。例如应用 TB ,HYPER,1,,3,POLY 定义3参数模型。

与高阶Mooney-Rivlin选项相似，本选项对高应变水平可提供较好的近似。 在 NPTS =1，常数

=0，这一选项等价于Neo-Hookean选项(用户可参见§4.3.2.3)。在

NPTS =1时，本选项等价于2参数Mooney-Rivlin选项。在 NPTS =2时，来选项等价于5

参数Mooney-Rivlin选项。在 NPTS =3时，本选项等价于9参数Mooney-Rivlin选项(参见§4.3.2.1)。对于本选项所要求的材料常数，参见《ANSYS Elements Reference》。 4.3.2.5 Arruda-Boyce超弹性选项

Arruda-Boyce选项( TB ,HYPER,,,,BOYCE)可用于直到300%的应变水平。下面是本选项的一个例子:

TB,HYPER,1,,,BOYCE !Activate Arruda-Boyce data table

TBDATA,1,200.0 !Define initial shear modulus TBDATA,2,5.0 !Define limiting network stretch

TBDATA,3,0.001 !Define incompressibility parameter !(as 2/K, K is the bulk modulus)

对于本选项所要求的材料常数，参见《ANSYS Elements Reference》。 4.3.2.6 用户定义超弹性选项

用户定义选项( TB ,HYPER,,,,USER)允许用户应用子程序USERHYPER来定义应变能势对应变不变量的导数，参见《ANSYS Guide to User Programmable Features》。 4.3.2.7 Mooney-Rivlin超弹性选项(TB,MOONEY)

请注意这一选项适用于HYPER56,HYPER58,HYPER74,HYPER158,PLANE162, SHELL163,SOLID164 等单元。

如果要应用Mooney-Rivlin选项于SHELL181, PLANE182, PLANE183, SOLID185, SOLID186, SOLID187等单元，请参见§4.3.2.1 。

ANSYS的单元类型HYPER56,HYPER58,HYPER74,HYPER158，应用直到9个

Mooney-Rivlin弹性势能函数。如果用户已知2项、3项、5项或9项Mooney-Rivlin常数的值，则可以通过 TB 族命令直接输入，参见《ANSYS Theory Reference》中的Mooney-Rivlin函数。对于这些单元，用户也可以指定材料函数为用户可编程特性，参见《ANSYS Guide to User Programmable Features》。 下面是一个示例：

MP,NUXY,1,0.49999 ! NUXY should be almost equal to, but less than 0.5 TB,MOONEY,1,1 TBDATA,1,0.163498 TBDATA,2,0.125076 TBDATA,3,-0.0047583 TBDATA,4,0.014719 TBDATA,6,0.0003882

! (Constants 5, 7, 8, and 9 default to 0.0 in this example)

对于本选项所要求的材料常数，参见《ANSYS Elements Reference》。

对于任意给定的超弹性材料的Mooney-Rivlin常数，通常在公开文献中查不到。因而，用户可以用 * MOONEY 命令，从一组已知试验数据中自动生成Mooney-Rivlin常数。有时材料制造商可能提供所需的试验数据的一部分或全部，但用户可能发现还需要从试验取得更多的数据。

超弹性材料行为要比金属材料行为复杂得多。超弹性应力-应变关系，通常在拉伸、压缩和剪切变形中明显不同。因此，应用 * MOONEY 命令来生成一般应用的超弹性材料模型时，需要使用所有可能的三种变形模式，即拉伸、压缩和剪切。参见《ANSYS Theory Reference》中讨论的超弹性测试方法和相应的变形模式。 如果所得到的是一组不完全的数据(如只有单轴拉伸数据)，程序仍然可以确定出可用的超弹性材料特性。然而，在这种情况下，模型变形特征将只限于与测试特征完全一样的。换句话说，测试数据应该代表计算模型中所有变形模式和响应(应变)范围。

上面的劝告只是一个简单的说明。如果用户不知道某一个变形模式或应变范围，就不能准确地预测模型中具有这些变形或应变部分的行为。例如，如用户只知道单轴拉伸测试数据，就不能建立那种有很大剪切变形的模型；如只知道应变值为0%-100%之间的测试数据，就不要建立那种有150%应变的模型。如在计算完成之后用户发现所有测试数据不足以表征模型的响应，唯一的改正办法就是获取更多的测试数据。

用户可用 * MOONEY 命令来自动从实验数据中确定一组Mooney-Rivlin常数。ANSYS确定这些常数以一个数组的形式保存到数据库中。此外，程序还将这些常数以多个 TB 和 TBDATA 命令的格式把Mooney-Rivlin常数写到一个文本文件(Jobrame.TB)中。一旦形成这样的文件，用户就可以在将来的分析用于定义某些相同的Mooney-Rivlin常数。而毋需每次应用 * MOONEY 命令来生成这些常数。 确定和应用Mooney-Rivlin常数

计算并应用Mooney-Rivlin常数有如下五个步骤： 第一步：定义数组 命令：* DIM

GUI：Utility Menu>Parameters>Array Parameters>Define/Edit

用户在使用 * MOONEY 命令(GUI：Main Menu> Preprocessor>Material

Props>Mooney-Rivlin>Calc Constants)之前，必须定义数组，在绝大多数情况下，需至少定义六个不同的数组(用户可以给这些数组任意合法的参数名，但为了方便我们在这里应用特定的数组名，如STRAIN、SDTRESS等，用户可以用任何喜欢的有效参数名代替之)。它们是STRAIN, STRESS, CONST, CALC, SORTSN, 和SORTSS。

应变数组(STRAIN) ：这是从材料试验得到的工程应变的数组，分为三列： 第一列：单轴拉伸和/或压缩数据 第二列：等双轴拉伸和/或压缩数据

第三列：剪切数据(平面拉伸和/或压缩数据) 该数组的大小为 N ×3，其中， N 等于在三列测试中数据点个数的最大值。例如，如果从单轴拉伸/压缩数据得到20个数据点，从剪切试验得到10个数据点，则 N =20。即使只用了一种或两种测试，该数组的大小也必须为 N ×3。虽然以升序输入数据点是较好的，但不是必须的。

应力数组(STRESS) ：这是从材料试验得到的工程应力的数组。该数组的大小也是 N ×3，应力数据点的输入顺序必须与应变数据点的输入顺序完全一致。 常数数组(CONST) ：Mooney-Rivlin常数数组的大小为 M ×1， M 为所希望的常数的个数( M 必须为2，5，9三个数之一，如用其它数，则在应用 * MOONEY 命令时，会导致错误信息)。对该数组所作的定义同时就告诉了程序需要生成多少个Mooney-Rivlin常数，* MOONEY 命令自动读取该数组的大小并确定要生成多少个常数，并将其值写入该数组中。 确定Mooney-Rivlin常数的个数

概略地说，数据点的个数(即上面的 N )应至少为Mooney-Rivlin常数个数的两部。常数个数越多，曲线的统计量越与真实值相接近(即拟合得更好)，但曲线的形状可能会比常数个数少的曲线要差。鉴于此，用户可以按顺序尝试2项、5项、9项函数，并检查其所生成的应力-应变曲线以确定到底哪一个函数在综合曲线形状以及拟合质量两方面做得最好。

表4-1 建议的Mooney-Rivlin常数

应力-应变曲线中的点数 无拐点(即单个曲线) 一个拐点(即2条曲线) 2个拐点 建议的Mooney-Rivlin函数 2项 5项 9项

图4-14 典型的超弹性应力-应变曲线

输出应力数组(CALC) ：其大小为 N ×3， N 与前面的一样，该数组保存计算出来的应力值，这些应力值保存的顺序与其相应的应变保存值的顺序一致(后者以升序排列)。 排序应变数组(SORTSN) ：其大小为 N ×3，它保存经过排序的输入应变。 排序应变数组(SORTSN) ：其大小为 N ×3，它保存经过排序的输入应力。 例如，如果任一个类型的测试数据包含直到20个数据点，而想生成5项的Mooney-Rivlin常数，则可以应用下面的命令来定义所需要的数组(记住：用户可以代入任意有效的参数名)：

*DIM,STRAIN,,20,3 ! Dim. array (STRAIN) for 20 input strain-data points *DIM,STRESS,,20,3 ! Dim. array (STRESS) for input stress data (20 pts.) *DIM,CONST,,5,1 ! Dim. array (CONST) for 5-term M-R constants *DIM,CALC,,20,3 ! Dim. array (CALC) for sorted calculated stresses *DIM,SORTSN,,20,3 ! Dim. array (SORTSN) for sorted input strain data *DIM,SORTSS,,20,3 ! Dim. array (SORTSS) for sorted input stress data 参见 * DIM 命令的说明。 第二步：填充输入数据数组

当数组定义完成之后，就可以用 * SET 命令(GUI：Utility Menu>

Parameters>Array Parameters)将实验数据填入STRAIN数组和STRESS数组。请再一次记住，用户可以给这些数组任意有效的参数名；在这里所采用的参数名，只是为了讨论方便。 注意--*MOONEY 命令将所有输入应力和应变都解释为工程应力和工程应变。 这些数组大小都为 N *3，其每一列各自表示一类测试数据，顺序是：

第一列：单轴拉伸和/或单轴压缩； 第二列：等双轴拉伸和/或等双轴压缩； 第三列：剪切(平面拉伸或压缩)。

注意这并不是说，变形模式及其等同存在1：1的关系。第一个变形模式--单轴拉伸--与等效双轴压缩等同，但数组第一列包含从单轴拉伸和/或单轴压缩得到的数据。类似地，第二个变形模式--等双轴拉伸--与单轴压缩等同，但数组第二列包含从等双轴拉伸和/或等双轴压缩得到的数据。

表4-2 应力-应变输入数组中的数据位置 变形模态 单轴拉伸 等双轴拉伸 剪切 等效测试类型 单轴拉伸 等双轴压缩 等双轴拉伸 单轴压缩 平面拉伸 平面压缩 测试数据的在数组中的位置 第一列 第二列 第二列 第一列 第三列 第三列 如果只作了一种或两种测试，则须将未做的测试的相应的列置空。可用 图4-15 来说明。

图4-15 在应力和应变输入数组中的数据位置

现在考虑一个得到了单轴拉伸和剪切测试数据时的情况。在数组中存储应变和应力的命令可能象下面那样(当然，数组可以有任意的名字，而且在本例中用N1和N2表示的数据点数，可以是任意整数)：

! Uniaxial Tension Data

*SET,STRAIN(1,1), ... ! First 10 strain data points

*SET,STRAIN(11,1), ...! Strain data points 11 through N1 (if N1<21) *SET,STRESS(1,1), ... ! First 10 stress data points

*SET,STRESS(11,1), ...! Stress data points 11 through N1 ! Shear Data

*SET,STRAIN(1,3), ... ! Strain data points 1 through N2 (if N2<11) *SET,STRESS(1,3), ... ! Stress data points 1 through N2

参见 * SET 命令的说明。

第三步：计算Mooney-Rivlin常数 要自动生成Mooney-Rivlin常数，首先应执行 TB 命令，并使该命令的 Lab =MOONEY， TBOPT =1。然后，执行 * MOONEY 命令，将已生成好的数组名填入其中(用户可以给这些数组任意有效的参数名，在这里所采用的参数名，只是为了讨论方便)：

TB,MOONEY,MAT,NTEMP,,1

*MOONEY,STRAIN(1,1),STRESS(1,1),,CONST(1),CALC(1),SORTSN(1), SORTSS(1),Fname,Ext

程序自动计算出Mooney-Rivlin常数，将它们写入数据库和 CONST 数组(可以是任何有效的数组名)中，并还以 TB 和 TBDATA 命令的格式写入一个ASCII文件 Fname.Ext 中(缺省是Jobname.TB)。

单轴公式将用于第一列的数据，而等双轴公式将用于第二列，平面(纯剪)公式将用于第三列。

注意 --在 STARIN 和 STRESS 中输入的所有试验数据，将用于确定Mooney-Rivlin 超弹性材料常数。

第四步：估计Mooney-Rivlin常数的质量

在输出文件中(Jobname.OUT)，检查“ROOT—MEAN—SQUARE

ERROR(PERCENTAGE)”(均方根误差)和“COEFFICIENT OF DETERMINATION”(确定系数)两个输出信息，这两个值给出所计算出的应力-应变曲线与测试数据点相拟合的好坏程度的统计度量。以百分数表示的均方根误差(即2.5就表示2.5%应接近于零，确定系数接近于1.0(通常要大于0.99)。 另外，用户还应使用 * VEAL 和 * VPLOT 命令(GUI:Main Menu> Preprocessor>Material Props>Mooney-Rivlin>Evaluate Const 和 Utility Menu>Plot>Array Parameters)来以图形的方式显示输入的和计算出的应力-应变曲线，以直观地检查计算曲线与实验数据的匹配程度。在比较这些曲线时，应比较那些代表同一变形模式的数据。也就是说，所计算的单轴拉伸曲线形状(在 * EVAL 命令中的 EVPARM = 1 )，仅应当与单轴拉伸数据(在排序的STRAIN 和 STRESS数组第一列)比较。类似地，所计算的单轴压缩曲线形状，仅应当与单轴压缩数据比较；而所计算的剪切曲线形状，仅应当与剪切数据比较。

当用图形显示计算的应力-应变曲线时，用户可以将显示曲线扩展到那些没有实验数据的区域。这样就可以让用户对模型在超出实验数据区域以外的地方的响应有一个定性的认识。但是，应认识到当将显示曲线扩展到一个表示另一个不同的变形模式的区域时，在那个区域的显示就是没有意义的。例如，用户只能在正应变区域显示单轴拉伸曲线。总的来说，要得

到一个好的结果，所作的实验数据应能代表所分析模型的所有的变形模式和响应(应变)范围。

* MOONEY 命令自动把Mooney-Rivlin常数写到 CONST 数组中。因为 * EVAL 命令从 CONST 常数读入相同的常数，所以可以在同一个ANSYS 阶段中，在 * MOONEY 命令后，跟着用 * EVAL 命令。如已经有了Mooney-Rivlin常数(这时不必进行 * MOONEY 计算)，则必须在计算曲线前，首先定义CONST 数组 [* DIM ]并用Mooney-Rivlin常数填充该数组[* EVAL ]。可以方程容易地填充这一数组，通常给 CONST 数组1×2、1×5、1×9。用户也可以在Jobname.TB文件中，添加 * DIM 和数组填充命令，以方便操作。

要检查曲线的形状，首先还必须定义[* DIM ]两个表数组向量(* EVAL 命令中将这两个数组名定义为 XVAL 和 ECALC ，但可以应用任何有效的参数名)。这两个表数组的大小都是P维，此时P为用户想要在曲线中绘制的点的数目(通常应使用一个相对较大的P值，以使得所绘制的曲线尽可能光滑)。其次，定义变形模式、定义应变范围，用 * EVAL 命令将工程应变和计算出的工程应力数据填入数组中。最后，用 * VPLOT 命令来绘制计算出的应力-应变曲线。下面例子说明对单轴压缩变形模式的计算曲线绘图：

! Dimension strain and stress arrays for the calculated curve: ! (Any valid parameter names can be used) *DIM,XVAL,TABLE,1000 *DIM,ECALC,TABLE,1000

! Specify the mode of deformation (EVPARM), define the strain range

! (XMIN,XMAX), and use the M-R constants (CONST) to fill the strain (XVAL) ! and stress (ECALC) arrays with calculated data:

*EVAL,1,2,CONST(1),XMIN,XMAX,XVAL(1),ECALC(1) ! Label the graph axes:

/AXLAB,X,Engineering Strain /AXLAB,Y,Engineering Stress

! Plot the calculated uniaxial compression curve: *VPLOT,XVAL(1),ECALC(1)

参见 *DIM, *EVAL, /AXLAB, 和 *VPLOT 等命令的说明。

图4-16 典型的计算出的超弹性应力-应变曲线

第五步：使用Mooney-Rivlin常数

如果用户对曲线拟合的统计数值和整个曲线的形状都非常满意，则可使用生成的Mooney-Rivlin材料性质作后续的分析(* MOONEY 命令将在数据库中保存这些常数)。在将来的另外一些分析中，如果使用同样的材料模式，则只需用/ INPUT 命令读入文件“Jobname.TB”，即可将常数加载到新的数据库中。但别忘了定义材料泊松比的值

[ MP ,NUXY,...]。用户始终要记住的一点是，应对分析结果作仔细检查，以确定原始测试数据是否涵盖了模型的变形模式和最大应变。

含有超弹性单元的分析，有时对材料性质定义和载荷施加方式非常敏感。通常是某些Mooney-Rivlin常数会导致非常稳定的刚度矩阵而其它的常数则相反。因此，应当根据经验，并特别仔细地选取常数。

ANSYS内部提供超弹性材料(基于用户输入的Mooney-Rivlin 常数)稳定检查的功能。这些检查分两个级别：

在分析前进行第一次稳定检查。对6个典型的应力路径(单轴拉伸和压缩，等双轴拉

伸和压缩，平面拉伸和压缩)，延伸率范围0.1～10来进行检查。如果材料在这一范围不稳定，则出现一个提示，统计在材料变得不稳定时，名义应变的临界值，列出用户输入的Mooney-Rivlin 常数。如果材料在这一范围稳定，则不出现提示。下面的警告信息，列出材料1变得不稳定时的名义应变，然后列出用户输入的Mooney-Rivlin 常数：

*** WARNING *** CP= 1.110 TIME= 16:59:52 Material 1 can become unstable under certain loading.

The strain (nominal) limits where the material becomes unstable are: UNIAXIAL TENSION 0.645E+00 UNIAXIAL COMPRESSION -0.565E+00 EQUIBIAXIAL TENSION 0.516E+00 EQUIBIAXIAL COMPRESSION -0.220E+00 PLANAR TENSION 0.585E+00 PLANAR COMPRESSION -0.369E+00 Mooney-Rivlin constants of the hyperelastic material are: 0.170E+02, 0.000E+00, 0.150E+03 0.000E+00, 0.000E+00, 0.000E+00 0.000E+00, 0.000E+00, 0.000E+00 对于混合U-P公式的超弹性单元(HYPER56、HYPER58、HYPER74 和HYPER158)，

在设置 KEYOPT(8) = 1 时，也可以在ANSYS分析期间，执行稳定检查。对于每次平衡迭代，程序检查每个高斯点的稳定情况。如未通过稳定检查，则将在ANSYS输出窗口的求解历史阶段看到一个信息，报告该迭代上不稳定的高斯点总数。如问题通过检查，则不显示任何信息。下面的例子说明在分析期间，ANSYS检查出3个高斯点超过材料稳定极限：

DISP CONVERGENCE VALUE = 22.81 CRITERION= 0.5000

EQUIL ITER 1 COMPLETED. NEW TRIANG MATRIX. MAX 自由度 INC= 10.00

FORCE CONVERGENCE VALUE = 0.5018E+07 CRITERION= 0.1917E+05 >>> 3 Gauss points have exceeded the material stability limit

对于包含最少1个不稳定高斯点的单元，不稳定指标设置为1，并保存在结果文件中的SMISC记录中。在POST1中，可以绘出这个标识为STFLAG的SMISC

记录来显示不稳定区域。参见下列单元和不稳定指标STFLAG相应的SMISC表项：HYPER56 单元 ETABLE 和 ESOL 命令的Item and Sequence Numbers；HYPER74 单元 ETABLE 和 ESOL 命令的Item and Sequence Numbers ；HYPER158 单元 ETABLE 和 ESOL 命令的Item and Sequence Numbers。

用户应当清楚，虽然材料未通过稳定检查是收敛困难的一个原因，但并不说明一旦材料进入不稳定区域解就无效。材料稳定检查只是在收敛失效时，一个帮助用户诊断问题的简单工具。

对于许多接近不可压缩材料(泊松比大于0.49)，我们推荐应用混合U-P列式的超弹性单元(HYPER56、HYPER58、HYPER74 和 HYPER158)。

注意 --HYPER84和HYPER86号单元主要是用于模拟可压缩的、泡沫状的高弹体，它们使用Blatz-Ko函数来描述材料特性。将这些单元的KEYOPT(2)设置为1即可选取Blatz-Ko选项，然后用 MP 命令输入合适的 EX 和 NUXY 值，以定义初始的材料剪切模量。对这两种单元来说，不可压缩超弹性材料选项也是可用的，但只限于2项 Mooney-Rivlin 公式，通常应使用HYPER56、HYPER58、HYPER74、HYPER158(而不是 HYPER84 或 HYPER86)来计算所有的不可压缩超弹性材料。

用超弹性单元作分析时可能对载荷施加的快慢很敏感。在大多数情况下，应缓慢地施加载荷，以避免在收敛过程中，使单元发生过度变形。求解过程中所遇到的每一个问题，都可能是独特的，需要特殊考虑。在加载过程中的不同时间点有时会发生分叉解问题，亦即两个或多个不同的几何外形都具有相同的最小势能，此时可用具有二分的自动时间步长功能[ AUTOTS ，ON]来尽量避免之。 4.3.3 超弹分析实例 4.3.3.1 问题描述

一个由铁板和橡胶体所构成的物体（ 图4.17 ），其底面固定，顶面受一垂直向下的均匀位移，分析所需要提供的外力与顶面位移的变化关系。由于结构和载荷都是轴对称的，因此我们取一截面，对它进行轴对称分析。由于结构变形很大，橡胶体与铁板之间可能发生接触，因此在分析中应考虑到接触。 4.3.3.2 问题详细说明 铁板材料特性： Ex=2e5

(泊松比)＝0.3 橡胶体材料特性： (泊松比)＝0.499

单轴压缩时的实验数据：

应变 应力 应变 应力 -0.45 -256 0.0 0.0 -0.4 -128 0.1 1 -0.35 -0.3 -0.25 -0.2 -0.15 -0.1 -0.05 -64 0.2 1.5 -32 0.3 2.0 -16 0.4 2.9 0.5 3.6 -8 0.6 5 -4 0.7 7.5 -2 0.8 9.7 -1 0.9 17 单轴拉伸时的实验数据：

图4-17 问题描述图

4.3.3.3 分析过程（GUI方法） 步骤一：建立计算所需要的模型

在这一步中，建立计算分析所需要的模型，包括定义单元类型，创建结点和单元，并将数据库保存为“hypelastic.db”，在此对这一过程不再详细。 步骤二：恢复数据库文件 “ hypelastic.db ” Utility menu：file>Resume from 步骤三：定义材料性质 1、选“Main Menu>Preprocessor>Material Props>Material Models”。出现“Define Material Model Behavior”对话框，选择Material Model Number 1。

2、在“Material Models Available”窗口，双击“Structural->Linear->Elastic-> Isotropic”。出现一个对话框。

1、对泊松比（NUXY）键入0.499。 2、单击OK。

3、在“Define Material Model Behavior”对话框的菜单中，选择Material->New Model，弹出对话框，在Define Material ID中输入2，然后单击OK。 出现Material Model Number 2。 4、对弹性模量（EX）键入2e5 5、对泊松比（NUXY）键入0.3 6、单击ＯＫ。

步骤四：定义数组并输入相应的实验数据。 1、择菜单路径：utility menu>Parameters,在其下拉菜单中，将鼠标指向Array Parameters, Array

Parameters的下拉菜单出现。 2、Define/Edit,对话框出现。

3、在对话框中单击Add，另一对话框出现。 4、指定Parameter name 为strn 。

5、对No. of rows,cols,planes,依次输入19，3，1

6、单击Apply,结束此对话框的输入，另一相同的对话框出现。 7、在对话框中指定Parameter name 为strss 。

8、对No. of rows,cols,planes,依次输入19，3，1

9、单击OK,结束此对话框的输入，回到Array Parameters对话框。 10、选择数组strn，然后单击Edit,给数组赋值的对话框出现。

11、在数组的第一列依次输入单轴压缩与拉伸的实验应变数据，然后选择菜单路径

File>Apply/Quit,结束应变实验数据的输入。

12、在Array Parameters对话框中选择数组strss，然后单击Edit,给数组赋值的对话框出现。 13、在数组的第一列依次输入单轴压缩与拉伸的实验应力数据，然后选择菜单路径

File>Apply/Quit,结束应力实验数据的输入。 步骤五：计算MOONEY常数

1、选择菜单路径Main Menu>Preprocessor>Material Props>Mooney-Rivlin> Define Table对话

框出现。

2、在Material reference number（材料参考号）中，键入1, 在Number of tempratures中，键

入1，然后单击OK。

3、择菜单路径Main Menu>Preprocessor>Material props>Mooney-Rivlin> Calculate constant 。

对话框出现。

4、Number of Mooney-Rivlin const，选择5。

5、定strain data array 为strn, 指定stress data array 为strss。 6、单击OK。

步骤六：进入求解器

选择菜单路径Main Menu>Solution 步骤七：加载

根据所给条件，施加适当的约束和载荷。在此不作详述，参考命令流文件。 步骤八：定义分析类型和分析选项：

1、选择菜单路径Main Menu>Solution>-Analysis Type-New Analysis. 2、单击“Static”来选中它然后单击OK。

3、择菜单路径Main Menu>Solution>Unabridged Menu>Analysis Options。 Analysis Options对话框出现。

4、单击Large deform effects option(大变形效应选项）使之为ON， 然后单击OK。 步骤九：设置输出控制选项 1、选择菜单路径：Main Menu > Solution>Unabridged Menu >Load step opts-Output ctrls > Solu

Printout。对话框出现

2、在“Item”中，选择“all items” 3、对“FREQ”,选择“Every Substep” 4、单击OK

5、选择菜单路径：Main Menu > Solution>Unabridged Menu >Load step opts-Output ctrls >

DB/Resuls File：对话框出现 6、在“Item”中，选择“all items” 7、对“FREQ”,选择“Every Substep” 8、单击OK

步骤十：设置载荷步选项

1、选择菜单路径Main Menu> Solution>Unabridged Menu>Load step opts-Time/Frequenc>

Time and substps。对话框出现。

2、对time at end of load step（载荷步终止时间）键入1 3、对“DELTIM”(Time step size) 输入0.03

4、将“AUTOTS”(automatic time stepping)设置为ON 5、对 Mininim time step size 输入0.01 6、对 Maxinum time step of size 输入0.3 步骤十一：打开时间步长预测器

1、选择菜单路径Main Menu> Solution>Unabridged Menu>Load step opts-Nonlinear> Predictor。

2、将predictor的状态设置为“ON”。 步骤十二：打开线性搜索

1、选择菜单路径Main Menu> Solution>Unabridged Menu>Load step opts-Nonlinear> Line Search。

2、将LNSRCH的状态设置为“ON”。 步骤十三：设置收敛性控制

1、选择菜单路径Main Menu> Solution>Unabridged Menu>Load step opts-Nonlinear>

Convergence Crit，对话框出现。 2、单击Add,下一级对话框出现。 3、对“TOLER”，输入0.2 4、对“MINREF”，输入200 5、单击OK。

步骤十四：进行求解 步骤十五：进行后处理

4.3.3.4 分析过程（命令流方法）

fini /cle r=180 l1=185 l2=74 h1=6 h2=50 h3=182 r1=10 d=50 /prep7 et,1,56,,,1 et,2,42,,,1 et,3,48,,1 et,4,58 et,5,45 et,6,49,,1 keyopt,6,7,1 keyopt,3,7,1 r,1,10000,,0.5 rect,0,l1,0,h1 rect,0,l2,0,h2 cyl4,0,h3,r,-90,,0

aovlap,all

asel,s,loc,y,0,h1 asel,r,loc,x,0,l2 aadd,all alls

asel,s,loc,y,h1,h3 aadd,all alls

lsel,s,loc,x,0 lsel,r,loc,y,0,h1 lcom,all alls

lsel,s,loc,x,0 lsel,r,loc,y,h1,h3 lcom,all alls

lsel,s,loc,y,h1 lsel,r,loc,x,0,l2 lcom,all alls

lsel,s,loc,y,h1,h2 lsel,r,loc,x,l2

*get,line1,line,,num,max alls

lsel,s,loc,x,l2,l1 lsel,r,loc,y,h2,h3

*get,line2,line,,num,max alls

lfillt,line1,line2,r1 alls al,1,4,3

asel,s,loc,y,h1,h3 aadd,all alls

lsel,s,loc,x,l2,l1

lsel,r,loc,y,h1+0.1,h3-0.1 lcom,all alls

lsel,s,loc,x,l2,l1

lsel,r,loc,y,h1+0.1,h3-0.1 *get,line3,line,,num,max alls

kl,line3,0.12 lsel,s,loc,x,0

lsel,r,loc,y,h1,h3

*get,line4,line,,num,max kl,line4,0.4 kl,line4,0.7 alls lstr,5,7 lstr,13,8

asel,s,loc,y,h1,h3 lsel,s,,,3,4 asbl,all,all alls

asel,s,loc,y,h1,h3 aatt,1,1,1

asel,s,loc,y,0,h1 aatt,2,1,2 alls

lesize,3,,,10 lesize,8,,,12,2 lesize,7,,,6 lesize,12,,,1 lesize,19,,,8 mshkey,1 amesh,5 amesh,3 amesh,2 amesh,1 amesh,6 alls

lsel,s,,,19 nsll,s,1

cm,targ,node alls lsel,s,,,6 nsll,s,1

cm,cont1,node lsel,s,,,8 nsll,s,1

nsel,r,loc,y,h1,130 cm,cont2,node alls

cmsel,s,cont1 cmsel,a,cont2 !cmsel,a,cont3 cm,cont,node

alls type,3 mat,1 real,1

gcgen,cont,targ alls

save,hypelastic,db resume,hypelastic,db mp,nuxy,1,0.499 mp,ex,2,2e5 mp,nuxy,2,0.3 *dim,strn,,19,3 *dim,strss,,19,3 *dim,const,,5 *dim,calc,,19,3 *dim,sortss,,19,3 *dim,sortsn,,19,3 *dim,ecalc,table,100 *dim,xval,table,100

strn(1,1)=-0.45,-0.4,-0.35,-0.3,-0.25,-0.2,-0.15,-0.1,-0.05 strss(1,1)=-256,-128,-64,-32,-16,-8,-4,-2,-1 strn(10,1)=0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9 strss(10,1)=0.0,1,1.5,2.0,2.9,3.6,5,7.5,9.7,17 tb,mooney,,,,1

*mooney,strn(1,1),strss(1,1),,const(1),calc(1,1),sortsn(1,1),sortss(1,1) *eval,1,1,const(1),-0.2,0,xval(1),ecalc(1) *vplot,xval(1),ecalc(1) fini /solu alls

nsel,s,loc,y,0 d,all,all,0 alls

nsel,s,loc,x,0 d,all,ux,0 d,all,uz,0 alls

nsel,s,loc,y,h3 d,all,ux,0 d,all,uz,0 d,all,uy,-d alls

antype,static nlgeom,on

nropt,,,on outpr,all,all outres,all,all autots,on time,1

deltim,0.03,0.01,0.3 cnvtol,f,,0.02,2 lnsrch,on pred,on alls solve fini

4.4 蠕变分析 4.4.1 蠕变理论 4.4.1.1 定义

蠕变是率相关材料非线性，即在常荷载作用下，材料连续变形的特性。相反如果位移固定，反力或应力将随时间而变小，这种特性有时也称为应力松驰，见 图4-18a 。

图4-18 应力松弛和蠕变

蠕变的三个阶段如 图4-18b 所示。在初始蠕变阶段，应变率随时间而减小，这个阶段一般发生在一个相当短的时期。在第二期蠕变阶段，有一个常应变率，所以应变以常速率发展，在第三期蠕变阶段，应变率迅速增加直到材料失效。

由于第三期蠕变阶段所经历的时间很短，材料将失效，所以通常情况下，我们感兴趣的是初始蠕变和第二期蠕变。ANSYS程序中的蠕变行为用来模拟初始蠕变和第二期蠕变。蠕变系数可以是应力、应变、温度、时间或其它变量的函数。

在高温应力分析中(如核反应堆等)，蠕变分析非常重要。例如，假设在核反应堆中施加了预荷载，以保证与相邻部件保持接触而不松开。在高温下过了一段时间后，预荷载将降低(应力松驰)，可能使接触部件松开。对于一些材料如预应力砼，蠕变也可能十分重要。最重要的是要记住，蠕变是永久变形。 4.4.1.2 理论介绍 蠕变方程：我们通过一个方程来模拟蠕变行为，此方程描述了在实验中观测到的主要特征（特别是在一维的拉伸实验中）。这个方程以蠕应变率的方式表示出来，其形式如下：

上式中，A、B、C、D是从实验中得到的材料常数，常数本身也可能是应力，应变，时间或温度的函数，这种形式的方程被称为状态方程。

上式中，当常数D为负值时，蠕应变率随时间下降，材料处于初始蠕变阶段，当D为0时，蠕应变率为常值，材料处于第二期蠕变阶段。 对于2－D或3－D应力状态，使用VON Mises方程计算蠕应变率方程中所使用的标量等效应力和等效应变。

对蠕变方程积分时，我们使用经过修改的总应变，其表达式为:

经过修改的等效总应变为:

其等效应力由下式算出：

其中：G=剪切模量＝ 等效蠕应变增量

由程序给出的某一种公式进行计算，一般为正值，如果在数据表中

，则使用的是衰减的蠕应变率而不是常蠕变率，但这个选项一般不被推荐，因为在

初始蠕变所产生的应力为主的情况下，它可能会严重的低估蠕变值。如果 用修正的等效蠕应变增量来代替蠕应变增量。

，程序使

其中：e=2.718281828（自然对数的底数）

下面是计算积分点的蠕应变率与弹性应变比率的公式：

将本次迭代的所有单元的所有积分点的的最大值记为

，并且作为“CREEPRATIO”输出。

计算出等效蠕应变增量后，可将它转换成分量的形式，假设 Nc 是某个特定单元类型的应

变分量的个数。 如果

则有：

如果

，则有：

上式中，前三个为正应变分量，第四个是剪应变分量。 如果

，前四个分量与上式相同，另两个剪应变分量为：

接下来，可以按下式来计算弹性应变和总的蠕应变（以 X 方向的分量为例）：

为了从标量 来计算分量 , , ,程序使用相关流动准则：Prandtl-Reuss方程，

与塑性应变相同，蠕应变只有偏差分量（剪分量），没有由于蠕变引起的体积应变。

为了考虑应力随时间的变化，使用两种强化准则，时间强化和应变强化。我们以一简单拉伸试验来说明：刚开始时，杆被加载到应力为

，在时间 它被卸载到应力为

。

(a) 时间强化 (b) 应变强化

图4-19 典型的单轴蠕变曲线

时间强化假定蠕应变率仅仅依赖于蠕应变过程开始的时间。当应力从

材料的蠕变率由点A表示（相当于曲线向上移动）。 应变强化假定蠕应变率仅仅依赖于材料中的应变，当应力从

变到

变到

时，

时，材料的蠕变率

由点B表示（相当于曲线左移）。大多数实验数据与应变强化准则吻合得更好。 4.4.2 求解算法

ANSYS使用隐式和显式积分二种方法来进行蠕变分析，均可应用于静态和瞬态分析。隐式蠕变分析方法更强大、更快、更精确，一般推荐使用隐式蠕变分析。它可以处理温度相关蠕变常数，同时模拟蠕变与等向强化塑性模型。

对于需要很小时间步的情况，显式蠕变分析就非常有用。蠕变常数不能有温度相关性，而与其他塑性材料模型的耦合只能应用迭加法。

注意 --蠕变分析中的“隐式”和“显式”，与“显式动力分析”或“显式单元”没有任何关系。 隐式蠕变分析方法支持下列单元：PLANE42,SOLID45,PLANE82,SOLID92, SOLID95,LINK180，SHELL181，PLANE182，PLANE183，SOLID185，SOLID186，SOLID187，BEAM188 和 BEAM189。

显式蠕变分析方法支持下列单元：LINK1,PLANE2,LINK8,PIPE20,BEAM23,

BEAM24,PLANE42，SHELL43，SOLID45，SHELL51，PIPE60，SOLID62，SOLID65，PLANE82，SOLID92 和 SOLID95。

蠕变应变率可以是应力、应变、温度、电子流水平的函数。蠕变应变率方程已按初始蠕变、第二期蠕变和辐射引起的蠕变在ANSYS中建立。参见《ANSYS Elements Reference》中关于这些蠕变方程的讨论和输入方法。有一些方程需要特殊的单位。特别是，对于显式蠕变选项，蠕变方程中的温度应当基于绝对温度。 4.4.2.1 隐式蠕变方法

隐式蠕变方法的基本步骤包括应用 TB 命令( Lab =CREEP)，通过 TBOPT 值选择蠕变方程。TBOPT的输入值对应于特定的蠕变方程, ANSYS程序所提供的隐式蠕变方程如下: · TBOPT=1 所对应的蠕变方程(初始蠕变方程):

· TBOPT=2 所对应的蠕变方程(初始蠕变方程): · TBOPT=3 所对应的蠕变方程(初始蠕变方程): , · TBOPT=4 所对应的蠕变方程(初始蠕变方程):

· TBOPT=5 所对应的蠕变方程(初始蠕变方程):

,

,

· TBOPT=6 所对应的蠕变方程(初始蠕变方程):

· TBOPT=7 所对应的蠕变方程(初始蠕变方程):

· TBOPT=8 所对应的蠕变方程(初始蠕变方程): · TBOPT=9 所对应的蠕变方程(二期蠕变方程): · TBOPT=10 所对应的蠕变方程(二期蠕变方程):

· TBOPT=11 所对应的蠕变方程(初始蠕变+二期蠕变方程):

· TBOPT=12 所对应的蠕变方程(初始蠕变+二期蠕变方程):

· TBOPT=100 所对应的蠕变方程: 用户自定义的蠕变方程

在以上方程中:

=等效蠕应变

=等效蠕应变对时间的变化率

=等效应力

T=绝对温度, 程序内部温度偏移量(TOFFST)被加到所有的温 度上。

=通过TBDADA命令所输入的材料常数

t=子步的结束时间。

下例说明隐式蠕变分析方法。 TBOPT =2表示将应用初始蠕变方程于模型2。温度相关性通过 TBTEMP 命令来指定，与此方程有关的4个常数作为 TBDATA 命令的参数。

TB,CREEP,1,1,4,2 TBTEMP,100

TBDATA,1,C1,C2,C3,C4

用户也可以应用ANSYS的可编程特性，并设置 TBOPT =100 来输入其他蠕变表达式。可以用 TB 命令( Lab =STATE)来定义状态变量数。下例是如何定义5个状态变量的例子：

TB,STATE,1,,5

用户可以同时模拟蠕变[ TB, CREEP]和各向同性强化、双线性随动强化和HILL各向异性塑性来考察更复杂的材料行为。参阅《ANSYS Element Reference》中的《Material Model Combination》部分来了解可用的联合使用。另外参阅本书§4.6《Material Model Combination》中材料联合使用的输入命令。

为了执行隐式蠕变分析，用户必须应用求解 RATE 命令( Option =ON或1)。下面的例子说明一个时间强化蠕变分析，见 图4-20 。

图4-20 时间强化蠕变分析

用户在第1荷载步施加机械荷载，并把 RATE 命令设为 OFF，这样绕过(忽略)蠕变应变效应。由于在这一荷值步的时间间隔将影响其后的总时间，因此这一荷载步的时间间隔要充分小。例如，用户可指定时间值为1E-8秒。第2荷载步是蠕变分析。这时应把 RATE 命令设为ON。这里机械荷载保持为常数，而材料随时间增量而发生蠕变。

/SOLU !First load step, apply mechanical loading RATE,OFF !Creep analysis turned off

TIME,1.0E-8 !Time period set to a very small value ...

SOLV !Solve this load step

!Second load step, no further mechanical load RATE,ON !Creep analysis turned on

TIME,100 !Time period set to desired value ...

SOLV !Solve this load step

RATE命令仅对采用von Mises 和Hill势的隐式蠕变有效。

当采用von Mises势模拟隐式蠕变时，可以对如下单元运用RATE命令：LINK180, SHELL181, PLANE182, PLANE183, SOLID185, SOLID186, SOLID187, BEAM88, BEAM189。

当模拟各向异性蠕变时（ TB , CREEP 和 TB , HILL），可以对如下单元运用RATE命令：PLANE42, SOLID45, PLANE82, SOLID92, SOLID95, LINK180, SHELL181, PLANE182, PLANE183, SOLID185, SOLID186, SOLID187, BEAM88, BEAM189。

对于大多数材料，在早期阶段，蠕变应变率显著改变。因为这一原因，通常建议应用很小的初始时间步增量，然后应用求解命令 DELTIM 或 NSUBST 指定较大的最大增量时间步。对于隐式蠕变，用户可能需要在结果中仔细检验时间增量的影响，因为ANSYS缺省并不提供任何蠕变率的控制。用户可以应用 CRPLIM 或 CUTCONTROL ,CRPLIMIT 命令中的蠕变率控制选项来总是强迫采用一个蠕变极限比率。蠕变极限比率的推荐值是1～10。该比率可以随材料而变化，以便用户可以根据自己的经验来决定一个最佳值，从而获得需要的运行和精度。对于大型分析，建议首先在一个小模型中对时间增量收敛分析进行测试。

4.4.2.2 显式蠕变方法

显式蠕变方法求解蠕应变使用了欧拉朝前法，以时间步开始时的应力、应变为基础计算出蠕应变率，在每个时间步长内，蠕应变率被假定是常数，因此有：

使用这种方法所对应的蠕应变曲线如下：

图4-21 显式蠕变对应的蠕应变曲线

为了减小误差，需要小的时间步长，特别是在蠕应变率变化很大的区域。

蠕应变率越小，结果越精确。一个等于或小于0.1的蠕变率将产生相当精确的结果。如果步长太大，求解将变得不稳定，并且不收敛，稳定极限对应于0.25的蠕应变率。 对于在第二期蠕变阶段蠕变较小的情况，欧拉朝前法是一种十分有效的方法，对于蠕应变较大的情况，这种方法需要时间步长很小，对于高度非线性蠕变应变-时间曲线，在应用显式蠕变方法时，必须应用小的时间步。如果时间步小于1E-6，则不计算蠕变应变。若需自动调整时间步到合适的值，可以应用蠕变时间步优化选项[ AUTOTS 和 CRPLIM ]。

如果没有其它的非线性行为则：

l 只在每一个子步的开始修正蠕应变率，在此子步内将不再变化。 l 不执行牛顿－拉普森迭代，因此，求解精度依赖于时间步长。 如果有其它的非线性行为则有：

l 既然蠕应变影响应力分布（反过来应力分布又影响蠕变率），因此在每个时间步内都要

进行牛顿－拉普森迭代

l 缺省情况下，使用初始刚度的牛顿－拉普森选项（总是使用弹性刚度矩阵）。 l 既然时间步长很小，使用切线矩阵并不能使求解过程更有效。

显式蠕变方法的基本步骤，包括应用 TB 命令( Lab =CREEP)，选择蠕变方程(用

TBDATA 命令的参数加入适当的常数)。 TBOPT 为0或者为空白。下例是应用显式蠕变方法的输入。请注意，所有常数是作为 TBDATA 命令的参数加入的，而且无温度相关性。

TB,CREEP,1

TBDATA,1,C1,C2,C3,C4, ,C6

ANSYS程序中具有一个蠕变方程库供你选择，你可以根据材料特性选择相应的蠕变准则，也就是说，在ANSYS程序中，建立了某些特定材料的特定蠕变准则，特别是那些用于原子能工业方面的材料（例如304和316不锈钢，2 1/4 Cr-1锰低铝钢），另外，用户还可以建立自己的蠕变准则，应用程序的可编程特性加入用户蠕变表达式，见《ANSYS Guide to User Programmable Fealares》。

可通过数据表来输入蠕变材料常数：

Main Menu: Preprocessor>Material Props>Material Models 通过输入

的值进入到初始蠕变的计算，如果

将略过蠕变计算。通过输入

的＝

值进入到第二期蠕变阶段的计算，如果 将略过此阶段的计算，如果使用选项

9、10、11、13、14或15来进行初始蠕变的计算，由于在它们的公式中包括了第二期蠕变的计算，故第二期蠕变的计算被略过，通过输入 过输入

的值进入到辐射引起的蠕变计算，通

＝100，进入到用户自己定义蠕变准则的计算。下面我们详细说明一下ANSYS

程序所提供的蠕变准则。 ·

时的初始蠕变方程。

＝等效应变（以修改的总应变为基础）

＝等效应力

T＝绝对温度（所有给定温度加上TOFFST） t＝在子步结束的时间 e＝自然对数的底数

·

时的初始蠕变方程

·

时的初始蠕变方程。

·

的初始蠕变方程， 适用于退火304不锈钢：

可以使用几种不同的蠕变方程来计算 双指数蠕变方程

，输入

为了使用下面的双指数蠕变方程来计算

Psi(缺省)

S,r,

,G和H是温度和应力的函数。

－1100

Psi(缺省)

对于退火304不锈钢，当温度在800

时，上面的双指数方程是有效的，上式

中的前两项描述了初始蠕变，最后一项描述了第二期蠕变。 为了使用这个方程，输入一个非零的 用

（或

值，

，

，

。温度应该使

且TOFFST=460.0），如果温度低于有效范围，则不计算蠕变，时间以小时

为单位，应力使用Psi,有效的应力范围为6000-25000Psi。 使用公制单位的标准有理多项式蠕变方程（

＝1）

，输入

为了使用公制单位的标准有理多项式蠕变方程来计算

c＝初始蠕应变的几极值 p＝初始蠕变的时间因子

＝第二期蠕变阶段的蠕变率（最小蠕变率）

对于退火的304不锈钢，当温度在427

－704

时，上面的标准有理多项式蠕变方程

是有效的，上式中第一项描述了初始蠕变，最后一项描述了第二期蠕变。 为了使用上面的方程，输入

＝1.0，

=1.0,

=9.0和

＝0.0，温度必须用

为单

位（TOFFST＝273），如果温度低于有效范围，则不计算蠕变，时间以小时为单位，应力用Mpa。 通过选择 变强化，

的值来控制蠕变的强化准则， ＝2.0选择初始蠕变强化。

＝2.0）

＝2.0。此时温度应以

。

＝0.0,选择时间强化.

＝1.0选择总的蠕应

使用英制单位的有理多项式蠕变方程（

为了在英制单位下使用上面的标准有理多项式蠕变方程，输入

为单位（TOFFST＝460），应力为PSI，有效温度范围为800－1300 ·

＝10的初始蠕变方程。 适用于退火316不锈钢：

与退火304不锈钢相同，也可以使用几种不同的蠕变方程来计算

=0时的双指数蠕变方程。

使用与退火的304不钢相同的双指数蠕变方程（

＝9.0，

＝0.0），与退火304不锈

缺省为

。

钢所用方程不同的是有效应力范围为4000－30000PSI， 30000PSI。

使用公制单位的标准有理多项式蠕变方程（ 与

＝1.0）

缺省为4000PSI，

所描述的退火304不锈钢所用的标准有理多项式蠕变方程（

。

＝9.0，

＝1.0）

相同。不同的是：其有效温度范围是482－704 使用英制单位的有理多项式蠕变方程（ 与

＝2.0）

所描述的相同，所不同的是：温度必须以F为单位，TOFFST＝460，应力以PSI

。

为单位，且有效范围为0.0-24220PSI，温度范围为900－1300 ·

＝11的初始蠕变方程。

适用于退火2 1/4 Cr-1 Mo 低合金钢

可以使用几种不同的蠕变方程来计算 修改的有理多项式蠕变方程（

。

＝0.0）

，需输入

＝0.0：

为了使用下面修改的有理多项式蠕变方程来计算

A、B、

是温度和应力的函数。

范围内时，修改的有理多项式方程是有

对2 1/4 Cr-1 Mo低合含钢，当温度在700－1100

效的，第一项描述了初始蠕变，第二项描述了第二期蠕变。

为了使用上式须输入 以

＝1.0， ， =11.0和 =0.0温度必须以

为单位（或

为单位且TOFFST＝460.0），时间以小时为单位，应力以PSI为单位，有效应力范

围是1000－65000PSI。

公制单位的标准有理多项式蠕变方程（

＝1.0）

c＝初始蠕应变的极值 p＝初始蠕变时间因子

＝第二期蠕变阶段的蠕应变率

－593

，上式是有效的。

为单位，TOFFST

对退火的2 1/4 Cr-1 Mo 低合金钢，当温度371 为了使用上式，输入

＝1.0，

=1.0,

=11.0和 =0.0。温度以

＝273。时间用小时为单位，应力为Mpa，强化准则与 英制单位的有理多项式蠕变方程（

＝2.0）

=9.0所使用的相同。

与标准有理多项式蠕变方程相同，除了温度以 效温度范围为700－1100

＝12的初始蠕变方程

。

为单位，TOFFST＝460，应力为PSI，有

＝比例常数 M、N、K＝温度函数

：描述M、N、K函数的温度值的个数 ：第一个绝对温度值 ：第二个绝对温度值

：第

个绝对温度的值

+1：第一个M的值 ：第

个M的值

+1：第一个N的值 ：第

个N的值

+1：第一个K的值

??

系数随温度而变的幂函数蠕变准则与

，

＝1的方程相似，但

，

，

。温度必须以降序输入。

=13的初始蠕变方程

＝积累的蠕应变 A＝ B＝ C= 常数

=14的初始蠕变方程

应该输入0，

不应该输入0。

对于退火的316不锈钢，当温度为800 式，输入

=1.0和

=14.0。温度采用

-1300

（或用

，上面的蠕变方程有效。为了使用上

且TOFFST＝460，时间以小时为

单位，常数仅适用于英制单位，应力范围为0.0-45000PSI。

=15 的初始蠕变方程。

一般材料的有理多项式蠕变方程：

上面的有理多项式蠕变方程是标准有理多项式蠕变方程

的一般形式，对于等温情况，此方程变为标准方程，强化准则与

＝100的标准蠕变方程。 通过输入

＝100建立用户自已定义的蠕变方程。

和11.0(

=1.0和2.0

=9.0的情况相同。

＝0的第二期蠕变方程

＝等效应力 T＝绝对温度 t＝时间

＝1的等二期蠕变方程

＝5的由辐射引起的蠕变方程

=等效应力 =中子通量

对于冷加工的316不锈钢，当温度在700－1300

时，上式有效。

4.4.3 蠕变分析实例 4.4.3.1 问题描述

一块矩形板，其左端固定，而右端被拉伸至某一固定位置，然后保持在此位置不动。试分析板中应力随时间的变化。 4.4.3.2 问题详细说明 材料特性： Ex=2e5

(泊松比)＝0.3

C6=0的显式初始蠕变方程： C1=4.8e-23 C2=7 几何特性： L=100 H=10

图4-22 问题描述图

4.4.3.3 求解步骤（GUI方法） 步骤一：建立计算所需要的模型

在这一步中，建立计算分析所需要的模型，包括定义单元类型，创建结点和单元，并将数据库保存为“creep.db”，在此对这一过程不再详细。 步骤二：恢复数据库文件 “ creep.db ” utility menu>file>Resume from 步骤三：定义材料性质 1、选“Main Menu>Preprocessor>Material Props>Material Models”。出现“Define Material Model Behavior”对话框，选择Material Model Number 1。

2、在“Material Models Available”窗口，双击“Structural->Linear->Elastic-> Isotropic”。出现一个对话框。

3、对杨氏模量（EX）键入2e5 。 4、对泊松比（NUXY）键入0.3。 5、单击OK。

步骤四：定义creep数据表并输入相应值

1、在“Material Models Available”窗口，双击Structural->Nonlinear->Inelastic->Rate Dependent->Creep->Creep Only->Mises Potential->Explicit，出现一个对话框。 2、在对话框表格中的C1，C2位置输入相应值(C1=4.8e-23,C2=7)。 3、单击OK

4、退出“Define Material Model Behavior”对话框。 步骤五：进入求解器

选择菜单路径Main Menu>Solution 步骤六：加载

根据所给条件，施加适当的约束和载荷。在此不作详述，参考命令流文件。 步骤七：定义分析类型：

1、选择菜单路径Main Menu>Solution>-Analysis Type-New Analysis. 2、单击“Static”来选中它然后单击OK。 步骤八：设置输出控制选项

1、选择菜单路径：Main Menu >Solution>Unabridged Menu >Load step opts-Output ctrls> Solu printout。对话框出现

2、在“Item”中，选择“all items” 3、对“FREQ”,选择“Every Substep” 4、单击OK

5、选择菜单路径：Main Menu > Solution>Unabridged Menu >Load step opts-Output ctrls > DB/Results File。对话框出现。 6、在“Item”中，选择“all items” 7、对“FREQ”,选择“Every Substep” 8、单击OK

步骤九：设置载荷步选项

1、选择菜单路径Main Menu> Solution>Unabridged Menu>Load step opts-Time/Frequenc> Time and substps。对话框出现。

6、对“Time”（载荷步终止时间）键入10000 7、对“NSUBST”(子步数) 输入100 8、将“KBC”(加载方式)设置为Stepped 步骤十：进行求解 步骤十一：进行后处理

4.4.3.4 求解步骤（命令流方法）

fini /cle l=100 h=30 n1=10 n2=3

/prep7 et,1,42 rect,0,l,0,h lsel,s,loc,y,0 lesize,all,,,n1 lsel,s,loc,x,0 lesize,all,,,n2 mshkey,1 mshape,0,2d amesh,all save,creep,db resume,creep,db mp,ex,1,2e5 mp,prxy,0.3 TB,CREEP,1

TBDATA,1,4.8E-23,7 ! CREEP PROPERTIES /solu

nsel,s,loc,x,0 d,all,all nsel,s,loc,x,l d,all,ux,0.1 alls

BFUNIF,TEMP,900 ! UNIFORM TEMPERATURE TIME,10000 KBC,1

NSUBST,100

OUTPR,all,all OUTRES,all,all alls SOLVE fini

4.5 其它材料非线性 4.5.1 粘塑性

粘塑性是时间相关的塑性现象，其中塑性应变的发展，与荷载率有关。其主要应用是高温金属形成过程，如轧制和深拉，这时涉及大塑性应变和大位侈小弹性应变。见 图4-23 。典型的塑性应变十分大(例如50%或更大)，需要大应变理论[ NLGEOM ,ON]。

粘塑性应用VISCD106、VISCO107和VISCO108单元来模拟，其材料特性应用Anand模型，见《ANSYS Elements Reference》§4.2。

图4-23 轧钢中的粘塑性行为

粘塑性通过一组流动和演化方程，由一致塑性和蠕变定义。用约束方程来保留塑性区的体积。

率相关粘塑性( TB ,RATE)选项，允许用户在材料模式中引入应变率效应，以模拟材料的时间相关响应。在ANSYS中可采用二种材料选项，即Perzyna模型和Peirce模型，见《ANSYS Theory Reference》。与ANSYS中的其他率相关材料选项如蠕变或Anand模型相反，Perzyna和Peirce模型还包括了屈服面。塑性以及应变率强化效应只在塑性屈服后才有效。用户必须把这两个模型与各向同性塑性强化选项相结合( TB ,BISO，或 TB ,MISO，或 TB ,NLISO)来模拟粘塑性。此外，与HILL选项结合还可以模拟各向异性粘塑性。参阅《ANSYS Element Reference》中的《Material Model Combination》。也可参见本书§4.6《材料的联合使用》来了解有关材料联合使用的输入命令。对于各向同性强化，其目的是模拟材料的应变率强化，而不是软化。这个选项也适用于大应变分析，并且适合于下列单元：PLANE42,

SOLID45, PLANE82, SOLID92, SOLID95, LINK180, SHELL181, PLANE182, PLANE183, SOLID185, SOLID186, SOLID187, BEAM189和 BEAM189。这种材料模式选项的一些典型应用是金属形成和微电机系统(MEMS)分析。

下面的例子说明应用RATE选项与BISO选项结果(在2个温度点)：

TB,BISO,1 ! Activate BISO data table TBDATA,1,C1,C2 ! Values for constants C1, ! C2

TB,RATE,1,,,PERZYNA ! Activate RATE data table TBTEMP,1 ! Define first temperature TBDATA,1,C1,C2 ! Values for constants C1, ! C2

TBTEMP,2 ! Define second temperature TBDATA,1,C1,C2 ! Values for constants C1, ! C2

4.5.2 粘弹性

粘弹性类似于蠕变，但粘弹性材料的变形有弹性部分和粘性部分，弹性变形是瞬时发生的，粘性变形随时间而发展，在荷载移走后，弹性变形可以恢复。一种常见的粘弹性材料是玻璃。有些塑性也可认为是粘弹性的。一种粘弹性响应如 图4-24 所示。

图4-24 粘塑性行为(Maxwell模型)

粘弹性适用于单元VISCO88和VISCO89。用户必须应用 TB 族命令输入材料特性。关于如何用 TB 族命令输入粘弹性材料特性，参见《ANSYS Elements Reference》中的粘弹性材料常数、《ANSYS Theory Reference》§4.6。 4.5.3 混凝土材料

混凝土材料用来描述脆性材料的破坏行为，它可以考虑脆性材料的开裂或压碎。可以用 TB ，CONCR命令来定义混凝土材料模型，这种材料仅适用于钢筋混凝土单元SOLID65。

在多轴应力状态下，混凝土材料的破坏准则可以表示为：