2000.4 哈尔滨建筑大学 王焕定教授制作 1
第三章 材料非线性
有限元分析
1 非线性弹性问题的有限单元法
2 弹塑性问题的有限单元法
2000.4 哈尔滨建筑大学 王焕定教授制作 2
1 非线性弹性问题的有限单元法
前提,材料处于弹性状态，但是应力 -应变
关系是非线性的。位移和应变是微小的。因此
),,(
2
1
ijjiij uu ???几何方程
klti j k lijklsi j k lij DD ???? dd ?? 或物理方程
象线性问题一样，设位移和应变分别为
ee BNu ??? ?? 则全量形式的应力为
es BD ??? )(?增量形式的应力为
eT BD ??? d)(d ?
2000.4 哈尔滨建筑大学 王焕定教授制作 3
同线性问题分析一样，可得单元刚度方程为
eeses kVBDBFFVB ????? )(d)(d
T
E
T ???? ??
进行先处理（定位向量）集成，可得
? ? ???? RPPUUKVB s EdT )(d?
与线性问题不同，上式是非线性的方程组，因
此要用第一章介绍的方法来求解。
1）切线刚度法 —— 牛顿法
?? ??? VBDBVBk TeT e d)(d,,,TT e ???? ???
集成
? ? ??? 0dT RVB ??非线性方程 用牛顿法求解时，切线刚度矩阵为 (这里认为 )
)(??? ?
2000.4 哈尔滨建筑大学 王焕定教授制作 4
经整体集成后，可得整体切线刚度矩阵，由
此可建立（自修正的）牛顿法迭代公式为
nnnnnTn UUURRKU ?? ???? ?? 11 )()(
式中 Rn是应力 σ n引起的结点力，因此
其中 σ n为第 n步位移对应的非线性单元应力。
讲义上列出切线刚度法分析的计算步骤，这
里不再赘述。 (P.22)
因为 R-Rn物理含义是 不平衡力,所以牛顿法
也可理解为按不平衡力修正位移，使不平衡力
足够小。
? ?? VBR nn dT ? 表示集成
2000.4 哈尔滨建筑大学 王焕定教授制作 5
5) 由式 （ 3,6） 进行迭代, 直到
满足精度要求 。
从此可得切线刚度法的计算步骤为：
0U 0 1U1）设 =,求线弹性解 ；
1U2）由 求各单元的应变、应力；
? dV)(TT BDB n?
eTk nTK
3）从 计算单元切线
刚度矩阵 并集装 ；
? dVT nB ? nR4）计算 并集装 ；
2000.4 哈尔滨建筑大学 王焕定教授制作 6
2）应力转移、初应力法 —— 修正牛顿法
nnnnTn UUURRKU ?? ???? ?? 110 )()(
为避免每次迭代形成切线矩阵并求解，以初
始切线矩阵（即线弹性的刚度矩阵）迭代，则
这相当于按弹性刚度分配不平衡力。迭代的
过程就是不断调整个单元的应力，使刚度弱的
单元不能承受的应力逐渐转移到刚度大的单元
或边界上，因此也称为,应力转移法,。它先
求位移修正值，然后求下一迭代步的位移。
因为初始切线刚度矩阵,故?? VDBBk
T d
T0
nnn
T RVBUK 0e
T0 d ?? ? ? ?
表示集成
2000.4 哈尔滨建筑大学 王焕定教授制作 7
式中
是第 n步非线性位移对应的 弹性应力 。 由此从
修正牛顿法迭代公式可得
nenn DBD ??? ??e
nnnT RRRUK ???? 010
因为 ?
?? VBR nn dT ?
非线性应力
所以若将 视作,初应力,，并记nn ?? ?
e
则 ? ?
? VσBR nnn )d-( eT ??
nnT RRUK ???? 10
表示集成
它是不断修改初应力，使趋于一常量（弹性应
力和真实应力之差）。因此也称 初应力法 。
2000.4 哈尔滨建筑大学 王焕定教授制作 8
4） 由 (3,8)求, 反复迭代, 直到 足够
小 。
1R?2U
从讲义式 （ 3,8）可得初应力法的计算步
骤为：
?? dVT0 DBBk e
0TK
1） 由 集装初始切线刚度矩
阵 ；
RUK ?10T2） 由 求线弹性的解；
1U ? )d V-( 11T ?? eB
nR?
3） 由 计算各单元的,并集
装 ；
2000.4 哈尔滨建筑大学 王焕定教授制作 9
又设第 n步单元非线性应力对应的 弹性应变 为
nn D ?? 1e ??
则非线性的应变可表为
nnn D ???? 10 )( ???
残余 (初 )应变
nnnnnn D 010e)( ??????? ????? ?
式中 也可视作,初应变,，由上式可得n
0?
因此单元刚度方程为
VDBFFVB nnn )d-)((d 0TET ???? ?? ???
有些问题的本构关系是用应力表示应变,即
)(??? ?
3）初应变法 —— 修正牛顿法
2000.4 哈尔滨建筑大学 王焕定教授制作 10
也即初应变可作为等效结点荷载考虑
? ?? VDBR nn d0T0 ?
nnT RRUK 010 ???
由此，象初应力法一样，可得迭代公式为
因为它是不断修改初应变，使趋于一常量（总
应变和弹性应变之差）。因此也称 初应变法 。
其求解步骤如 下页 所示。
表示集成
? ? ?? VDBR nnn d))((T0 ???
为了更好地掌握上述知识，讲义上举了一个
简单例子，用以说明切线刚度法和初应力法。
2000.4 哈尔滨建筑大学 王焕定教授制作 11
初应变法的计算步骤为：
?? dVT0 DBBk e 0TK1）由 集装初始切线矩阵 ；
RUK ?10T2）由 求线弹性的解；
nne D ?? 1??
nnen 0??? ??
3）计算单元线弹性应变 和单元初
应变 ；
? dV0T nDB ? nR?4） 由 集装 ;
1?nU n0?5）由（ 3,8）求 ;迭代直到 很小。
2000.4 哈尔滨建筑大学 王焕定教授制作 12
为了更好地理解，这一例子各种方法的示意
图如下所示。
2000.4 哈尔滨建筑大学 王焕定教授制作 13
2000.4 哈尔滨建筑大学 王焕定教授制作 14
2000.4 哈尔滨建筑大学 王焕定教授制作 15
2000.4 哈尔滨建筑大学 王焕定教授制作 16
2000.4 哈尔滨建筑大学 王焕定教授制作 17
2 弹塑性问题的有限单元法
材料应力应变非线性、岩土工程的开挖和施
工过程等对结果有重大影响的问题，都必须用
增量法来求解。
在增量荷载 ΔRm作用下，位移、应力、应变和
内变量等的增量分别为
mmmm kU,、,??
下一步迭代时的荷载水平为
设 m迭代步的结果已知，位移、应力、应变
和内变量等分别记作
mmmm kU ??????,、、
mmm RRR ???? 1
2000.4 哈尔滨建筑大学 王焕定教授制作 18
由于所讨论的是小变形问题，因此
或
11T1 d)( ??? ?? ? m
e
m
e
m
ee RVB ???
对弹塑性问题，本构关系为
??? dd)(d epp DDD ???
),,(
2
1
ijjiij uu ???几何方程
也即单元增量应变为 。象弹性问
题一样，第 m+1步单元刚度方程为
mem B ???? ?
11
E
T )()d( ?? ????? m
e
m
e
m
e
m
e RFFVB ???
0)(dT ???? ? memeeme RVB ?????
弹性矩阵
塑性矩阵
弹塑性矩阵
m
e
m
e
m
ee RVB ?? ? d)(
T ??? )0(?
精确解时
2000.4 哈尔滨建筑大学 王焕定教授制作 19
因此应力增量为
0d)( 11T1 ??? ??? ? ? mmm RVBU ??
由于在 ΔRm作用前后，单元高斯积分点的应力
状态可能是弹性，也可能是塑性。因此首先得
解决状态判别问题。
?
?
?
???
?
??? depD
基于此，由单元集成可得整体非线性方程为
0)(d)( T1 ???? ? ?? mmmm URVBU ?????
在第二章已讨论了应变空间的加、卸载的规
则，它们是
2000.4 哈尔滨建筑大学 王焕定教授制作 20
为进一步讨论状态判别，设在荷载 Rm作用下
应力和内变量对应一弹性状态，也即
其几何解释为,弹性应力增量 指向屈服面内侧
或相切时，反应是弹性的。否则是塑性加载，
反应是弹塑性的。
0),( ?? mmm kff ?
增加荷载 ΔRm后，转入塑性状态，即
ed,d,
ijkli j k l ijij fDfl ?? ?? ??
?
?
?
?
??
0 1
0 0)(
加载
卸载、中性变载
l
llH
弹性应力
0),( e1 ???? mmmm kff ???
弹性应力
2000.4 哈尔滨建筑大学 王焕定教授制作 21
也可由下式线性内插确定
因此可用下式确定中性变载点处弹、塑性部分
的比例因子 r
0)1( 1 ??? ?mm rffr
基于此，加荷载 ΔRm后，应力增量为
0),( e ?? mmm krf ???
)( 1 mmm fffr ??? ?
???
?
?
??
???
mm
mm
mm
m
mm
m r
r
m DDD
???
???
???
?
???
?
????? ddd epep
?
?
?
??
mm
mm r
m DrD
???
???
??? dep
弹性应力 转图
2000.4 哈尔滨建筑大学 王焕定教授制作 22
当 ΔRm足够小时 (Δεm足够小 )，上式可写为
mmm DrD ?????? p)1( ???
跳
转
式中 Dp是按 计算的塑性矩阵。
mmm kr,e??? ?
当 r=1时，反应是纯弹性的，可以是弹性到
弹性、塑性卸载到弹性或中性变载。
当 r=0时，应力应变是塑性的，是弹性到塑
性的加载。 0<r<1时，应力应变是弹塑性的。
转图
?
?
?
???
mm
mm r
m DDrD
???
???
??? )d( p ?
?
?
??
mm
mm r
m DD
???
???
??? dp
讲义上给出了由应变求应力的计算步骤，可
供编程时参考。 (P.31)
弹性应力
2000.4 哈尔滨建筑大学 王焕定教授制作 23
弹性弹性 塑性
fm
fm+1
弹性
转
18
有状态改变 无状态改变
弹性
塑性
由中性状态
改变到塑性
2000.4 哈尔滨建筑大学 王焕定教授制作 24
解决了由应变求应力的计算，下面再解决弹
塑性问题非线性方程组求解问题。
将其代入单元刚度方程，可得 memmmm BDD ?????? epep ??
首先，以 下的弹塑性矩阵 代替增
量应力应变关系中的,也即有“线性关系”mm
k,? mD
ep
epD
0)(d)( T1 ???? ?? memeememe RVB ???????
0d)( T1 ??? ?? mememe RVB ?????或
0)d( T ??? ? emme RVB ???
进行整体集成，可得（认为 ）0)( ?mU?
]0))(d( [ T ???? ? meemme URVB ????或
集成
2000.4 哈尔滨建筑大学 王焕定教授制作 25
由此，若记
上述方程可以用前面所介绍的初应力法等进行
求解。初应力法的迭代公式为
则切线刚度方程为
)]([1 mnmmnmmT URRUK ???? ????
? ?? )d( epT VBDBK mmT
0?? mmmT RUK ??
]0)( [ ??? mmmmT URUK ???或
表示集成
? ??? ? )d( T1 VBRR mmm ??
? ?? ))d-((,eT VBR nmn mnm ?????由此出发，讲义上给出了编程参考的步骤 (P.33)。
第三章 材料非线性
有限元分析
1 非线性弹性问题的有限单元法
2 弹塑性问题的有限单元法
2000.4 哈尔滨建筑大学 王焕定教授制作 2
1 非线性弹性问题的有限单元法
前提,材料处于弹性状态，但是应力 -应变
关系是非线性的。位移和应变是微小的。因此
),,(
2
1
ijjiij uu ???几何方程
klti j k lijklsi j k lij DD ???? dd ?? 或物理方程
象线性问题一样，设位移和应变分别为
ee BNu ??? ?? 则全量形式的应力为
es BD ??? )(?增量形式的应力为
eT BD ??? d)(d ?
2000.4 哈尔滨建筑大学 王焕定教授制作 3
同线性问题分析一样，可得单元刚度方程为
eeses kVBDBFFVB ????? )(d)(d
T
E
T ???? ??
进行先处理（定位向量）集成，可得
? ? ???? RPPUUKVB s EdT )(d?
与线性问题不同，上式是非线性的方程组，因
此要用第一章介绍的方法来求解。
1）切线刚度法 —— 牛顿法
?? ??? VBDBVBk TeT e d)(d,,,TT e ???? ???
集成
? ? ??? 0dT RVB ??非线性方程 用牛顿法求解时，切线刚度矩阵为 (这里认为 )
)(??? ?
2000.4 哈尔滨建筑大学 王焕定教授制作 4
经整体集成后，可得整体切线刚度矩阵，由
此可建立（自修正的）牛顿法迭代公式为
nnnnnTn UUURRKU ?? ???? ?? 11 )()(
式中 Rn是应力 σ n引起的结点力，因此
其中 σ n为第 n步位移对应的非线性单元应力。
讲义上列出切线刚度法分析的计算步骤，这
里不再赘述。 (P.22)
因为 R-Rn物理含义是 不平衡力,所以牛顿法
也可理解为按不平衡力修正位移，使不平衡力
足够小。
? ?? VBR nn dT ? 表示集成
2000.4 哈尔滨建筑大学 王焕定教授制作 5
5) 由式 （ 3,6） 进行迭代, 直到
满足精度要求 。
从此可得切线刚度法的计算步骤为：
0U 0 1U1）设 =,求线弹性解 ；
1U2）由 求各单元的应变、应力；
? dV)(TT BDB n?
eTk nTK
3）从 计算单元切线
刚度矩阵 并集装 ；
? dVT nB ? nR4）计算 并集装 ；
2000.4 哈尔滨建筑大学 王焕定教授制作 6
2）应力转移、初应力法 —— 修正牛顿法
nnnnTn UUURRKU ?? ???? ?? 110 )()(
为避免每次迭代形成切线矩阵并求解，以初
始切线矩阵（即线弹性的刚度矩阵）迭代，则
这相当于按弹性刚度分配不平衡力。迭代的
过程就是不断调整个单元的应力，使刚度弱的
单元不能承受的应力逐渐转移到刚度大的单元
或边界上，因此也称为,应力转移法,。它先
求位移修正值，然后求下一迭代步的位移。
因为初始切线刚度矩阵,故?? VDBBk
T d
T0
nnn
T RVBUK 0e
T0 d ?? ? ? ?
表示集成
2000.4 哈尔滨建筑大学 王焕定教授制作 7
式中
是第 n步非线性位移对应的 弹性应力 。 由此从
修正牛顿法迭代公式可得
nenn DBD ??? ??e
nnnT RRRUK ???? 010
因为 ?
?? VBR nn dT ?
非线性应力
所以若将 视作,初应力,，并记nn ?? ?
e
则 ? ?
? VσBR nnn )d-( eT ??
nnT RRUK ???? 10
表示集成
它是不断修改初应力，使趋于一常量（弹性应
力和真实应力之差）。因此也称 初应力法 。
2000.4 哈尔滨建筑大学 王焕定教授制作 8
4） 由 (3,8)求, 反复迭代, 直到 足够
小 。
1R?2U
从讲义式 （ 3,8）可得初应力法的计算步
骤为：
?? dVT0 DBBk e
0TK
1） 由 集装初始切线刚度矩
阵 ；
RUK ?10T2） 由 求线弹性的解；
1U ? )d V-( 11T ?? eB
nR?
3） 由 计算各单元的,并集
装 ；
2000.4 哈尔滨建筑大学 王焕定教授制作 9
又设第 n步单元非线性应力对应的 弹性应变 为
nn D ?? 1e ??
则非线性的应变可表为
nnn D ???? 10 )( ???
残余 (初 )应变
nnnnnn D 010e)( ??????? ????? ?
式中 也可视作,初应变,，由上式可得n
0?
因此单元刚度方程为
VDBFFVB nnn )d-)((d 0TET ???? ?? ???
有些问题的本构关系是用应力表示应变,即
)(??? ?
3）初应变法 —— 修正牛顿法
2000.4 哈尔滨建筑大学 王焕定教授制作 10
也即初应变可作为等效结点荷载考虑
? ?? VDBR nn d0T0 ?
nnT RRUK 010 ???
由此，象初应力法一样，可得迭代公式为
因为它是不断修改初应变，使趋于一常量（总
应变和弹性应变之差）。因此也称 初应变法 。
其求解步骤如 下页 所示。
表示集成
? ? ?? VDBR nnn d))((T0 ???
为了更好地掌握上述知识，讲义上举了一个
简单例子，用以说明切线刚度法和初应力法。
2000.4 哈尔滨建筑大学 王焕定教授制作 11
初应变法的计算步骤为：
?? dVT0 DBBk e 0TK1）由 集装初始切线矩阵 ；
RUK ?10T2）由 求线弹性的解；
nne D ?? 1??
nnen 0??? ??
3）计算单元线弹性应变 和单元初
应变 ；
? dV0T nDB ? nR?4） 由 集装 ;
1?nU n0?5）由（ 3,8）求 ;迭代直到 很小。
2000.4 哈尔滨建筑大学 王焕定教授制作 12
为了更好地理解，这一例子各种方法的示意
图如下所示。
2000.4 哈尔滨建筑大学 王焕定教授制作 13
2000.4 哈尔滨建筑大学 王焕定教授制作 14
2000.4 哈尔滨建筑大学 王焕定教授制作 15
2000.4 哈尔滨建筑大学 王焕定教授制作 16
2000.4 哈尔滨建筑大学 王焕定教授制作 17
2 弹塑性问题的有限单元法
材料应力应变非线性、岩土工程的开挖和施
工过程等对结果有重大影响的问题，都必须用
增量法来求解。
在增量荷载 ΔRm作用下，位移、应力、应变和
内变量等的增量分别为
mmmm kU,、,??
下一步迭代时的荷载水平为
设 m迭代步的结果已知，位移、应力、应变
和内变量等分别记作
mmmm kU ??????,、、
mmm RRR ???? 1
2000.4 哈尔滨建筑大学 王焕定教授制作 18
由于所讨论的是小变形问题，因此
或
11T1 d)( ??? ?? ? m
e
m
e
m
ee RVB ???
对弹塑性问题，本构关系为
??? dd)(d epp DDD ???
),,(
2
1
ijjiij uu ???几何方程
也即单元增量应变为 。象弹性问
题一样，第 m+1步单元刚度方程为
mem B ???? ?
11
E
T )()d( ?? ????? m
e
m
e
m
e
m
e RFFVB ???
0)(dT ???? ? memeeme RVB ?????
弹性矩阵
塑性矩阵
弹塑性矩阵
m
e
m
e
m
ee RVB ?? ? d)(
T ??? )0(?
精确解时
2000.4 哈尔滨建筑大学 王焕定教授制作 19
因此应力增量为
0d)( 11T1 ??? ??? ? ? mmm RVBU ??
由于在 ΔRm作用前后，单元高斯积分点的应力
状态可能是弹性，也可能是塑性。因此首先得
解决状态判别问题。
?
?
?
???
?
??? depD
基于此，由单元集成可得整体非线性方程为
0)(d)( T1 ???? ? ?? mmmm URVBU ?????
在第二章已讨论了应变空间的加、卸载的规
则，它们是
2000.4 哈尔滨建筑大学 王焕定教授制作 20
为进一步讨论状态判别，设在荷载 Rm作用下
应力和内变量对应一弹性状态，也即
其几何解释为,弹性应力增量 指向屈服面内侧
或相切时，反应是弹性的。否则是塑性加载，
反应是弹塑性的。
0),( ?? mmm kff ?
增加荷载 ΔRm后，转入塑性状态，即
ed,d,
ijkli j k l ijij fDfl ?? ?? ??
?
?
?
?
??
0 1
0 0)(
加载
卸载、中性变载
l
llH
弹性应力
0),( e1 ???? mmmm kff ???
弹性应力
2000.4 哈尔滨建筑大学 王焕定教授制作 21
也可由下式线性内插确定
因此可用下式确定中性变载点处弹、塑性部分
的比例因子 r
0)1( 1 ??? ?mm rffr
基于此，加荷载 ΔRm后，应力增量为
0),( e ?? mmm krf ???
)( 1 mmm fffr ??? ?
???
?
?
??
???
mm
mm
mm
m
mm
m r
r
m DDD
???
???
???
?
???
?
????? ddd epep
?
?
?
??
mm
mm r
m DrD
???
???
??? dep
弹性应力 转图
2000.4 哈尔滨建筑大学 王焕定教授制作 22
当 ΔRm足够小时 (Δεm足够小 )，上式可写为
mmm DrD ?????? p)1( ???
跳
转
式中 Dp是按 计算的塑性矩阵。
mmm kr,e??? ?
当 r=1时，反应是纯弹性的，可以是弹性到
弹性、塑性卸载到弹性或中性变载。
当 r=0时，应力应变是塑性的，是弹性到塑
性的加载。 0<r<1时，应力应变是弹塑性的。
转图
?
?
?
???
mm
mm r
m DDrD
???
???
??? )d( p ?
?
?
??
mm
mm r
m DD
???
???
??? dp
讲义上给出了由应变求应力的计算步骤，可
供编程时参考。 (P.31)
弹性应力
2000.4 哈尔滨建筑大学 王焕定教授制作 23
弹性弹性 塑性
fm
fm+1
弹性
转
18
有状态改变 无状态改变
弹性
塑性
由中性状态
改变到塑性
2000.4 哈尔滨建筑大学 王焕定教授制作 24
解决了由应变求应力的计算，下面再解决弹
塑性问题非线性方程组求解问题。
将其代入单元刚度方程，可得 memmmm BDD ?????? epep ??
首先，以 下的弹塑性矩阵 代替增
量应力应变关系中的,也即有“线性关系”mm
k,? mD
ep
epD
0)(d)( T1 ???? ?? memeememe RVB ???????
0d)( T1 ??? ?? mememe RVB ?????或
0)d( T ??? ? emme RVB ???
进行整体集成，可得（认为 ）0)( ?mU?
]0))(d( [ T ???? ? meemme URVB ????或
集成
2000.4 哈尔滨建筑大学 王焕定教授制作 25
由此，若记
上述方程可以用前面所介绍的初应力法等进行
求解。初应力法的迭代公式为
则切线刚度方程为
)]([1 mnmmnmmT URRUK ???? ????
? ?? )d( epT VBDBK mmT
0?? mmmT RUK ??
]0)( [ ??? mmmmT URUK ???或
表示集成
? ??? ? )d( T1 VBRR mmm ??
? ?? ))d-((,eT VBR nmn mnm ?????由此出发，讲义上给出了编程参考的步骤 (P.33)。
