材料非线性问题的有限元
一、材料非线性问题的有限单元法
1．1 引言
以前各章所讨论的均是线性问题。线弹性力学基本方程的特点是
1．几何方程的应变和位移的关系是线性的。
2．物性方程的应力和应变的关系是线性的。
3．建立于变形前状态的平衡方程也是线性的。
但是在很多重要的实际问题中，上述线性关系不能保持。例如在结构的形状有不连续变化（如缺口、裂纹等）的部位存在应力集中，当外载荷到达一定数值时该部位首先进入塑性，这时在该部位线弹性的应力应变关系不再适用，虽然结构的其他大部分区域仍保持弹性。又如长期处于高温条件下工作的结构，将发生蠕变变形，即在载荷或应力保持不变的情况下，变形或应变仍随着时间的进展而继续增长，这也不是线弹性的物性方程所能描述的。上述现象都属于材料非线性范畴内所要研究的问题。工程实际中还存在另一类所谓几何非线性问题。例如板壳的大挠度问题，材料锻压成型过程的大应变问题等，这时需要采用非线性的应变和位移关系，平衡方程也必须建立于变形后的状态以考虑变形对平衡的影响。
由于非线性问题的复杂性，利用解析方法能够得到的解答是很有限的。随着有限单元法在线性分析中的成功应用，它在非线性分析中的应用也取得了很大的进展，已经获得了很多不同类型实际问题的求解方案。
材料非线性问题的处理相对比较简单，不需要重新列出整个问题的表达格式，只要将材料本构关系线性化，就可将线性问题的表达格式推广用于非线性分析。一般说，通过试探和迭代的过程求解一系列线性问题，如果在最后阶段，材料的状态参数被调整得满足材料的非线性本构关系，则最终得到问题的解答。几何非线性问题比较复杂，它涉及非线性的几何关系和依赖于变形的平衡方程等问题，因此，表达格式和线性问题相比，有很大的改变，这将在下一章专门讨论。这两类非线性问题的有限元格式都涉及求解非线性代数方程组，所以在本章开始对非线性代数方程组的求解作—一般性的讨论。这对下一章也是必要的准备。
正如在前面已指出的，材料非线性问题可以分为两类。一类是不依赖于时间的弹塑性问题，其特点是当载荷作用以后，材料变形立即发生，并且不再随时间而变化。另一类是依赖于时间的粘（弹、塑）性问题，其特点是载荷作用以后，材料不仅立即发生变形，而且变形随时间而继续变化，在载荷保持不变条件下，由于材料粘性而继续增长的变形称之为蠕变。另方面在变形保持不变条件下，由于材料粘性而使应力衰减称之为松弛。本章中重点讨论不依赖于时间的弹塑性问题，包括它的本构关系，有限元表达格式、求解步骤及数值方法。至于依赖于时间的蠕变问题，本章只简要地介绍其本构关系和求解步骤的基本特点，而不太多涉及其具体细节。
1.2 非线性方程组的解法
非线性问题有限元离散化的结果将得到下列形式的代数方程组
?
??
=+≡+≡=0)()()()(f a a K f a P a Q
a a K ψ （1．2．1）
其中f Q =-。该方程的具体形式通常取决于问题的性质和离散的方法。上式中参数a 代表未知函数的近似解。在以位移为未知量的有限元分析中，它是结点位移向量。
对于线性方程组0Ka f +=，由于K 是常数矩阵，可以没有困难地直接求解，但对于非线性方程组，由于K 依赖于未知量a 本身则不可能直接求解。以下将阐述借助于重复求解线性方程组以得到非线性方程组解答的一些常用方法。
1．2．1直接迭代法
对于方程（1.2.1）式
()0K a a f +=
假设有某个初始的试探解：
a a = （1．2．2）
代入上式的()K a 中，可以求得被改进了的一次近似解
f K a 101)(--= （1．2．3）
其中
00()K K a =
重复上述过程，可以得到n 次近似解
11()n n a K f --=- （1．2．4）
一直到误差的某种范数小于某个规定的容许小量r e ，即
1
n
n r e a a e -=-≤ （1．2．5）
上述迭代过程可以终止。
从（1.2.3）式和（1.2.4）式可以看到，要执行直接迭代法的计算，首先需要假设一个初始的试探解0a 。在材料非线性问题中，0
a 通常可以从先求解一线弹性问题得到。其次是直接迭代法的每次迭代需要计算和形成新的系数矩阵1
()n K a
-，并对它进行求逆计算。这里还隐含着K 可以显式地表示成a 的
函数，所以只适用于与变形历史无关的非线性问题，例如非线性弹性问题及可以利用形变理论分析的弹塑性问题。而对于依赖于变形历史的非线性问题，直接迭代法是不适用的，例如加载路径不断变化或涉及卸载及反复加载等必须利用增量理论分析的弹塑性问题。
关于直接迭代法的收敛性可以指出，当P （a ）－a 是凸的情况（当a 是标量，即系统为单自由度的，P —a 表示如图1．1（a ）），通常解是收敛的。但当P —a 是凹的情况（如图1·1（b ）所示），则解可能是发散的。
图1．1直接迭代法
1．2．2 Newton －Raphson 方法（简称N －R 方法）
如果方程（1.2.1）式的第n 次近似解n
a 已经得到，一般情况下（1.2.1）式不能精确地被满足，即
()0n a ψ≠。为得到进一步的近似解1n a +，可将1()n a ψ+表示成在n a 附近的仅保留线性项的Taylor 展开
式，即
0)()(1
=???? ??+≡+n
n
n n a da d a a
ψψψ （1．2．6）
且有
1n n n a a a +=+? （1．2．7）
式中/d da ψ是切线矩阵，即
()T d dP
K a da da
ψ≡≡ （1．2．8） 于是从（1．2．6）式可以得到
()()()f P K K a n n T n n
T
n +-=-=?--1
1
ψ （1．2．9） 其中
)(),(n n n T n
T a P P a K K ==
由于Taylor 展开式（1．2．6）式仅取线性项，所以1
n a +仍是近似解，应重复上述迭代求解过程直
至满足收敛要求。
N －R 方法的求解过程可以表示于图1．2。一般情况下，它具有良好的收敛性。当然像图1．2（b ）表示的那种发散情况也是可能存在的。
关于N －R 方法中的初始试探解0a ，可以简单地设00a =。这样一来，0
T K 在材料非线性问题中就
是弹性刚度矩阵。当然，从（1．2．9）式可以看到N －R 方法的每次迭代也需要重新形成和求逆一个新的切线矩阵n
T K 。
图 1．2 Newton －Raphson 方法
1．2．3修正的Newton －Raphson 方法（简称mN －R 方法）
为克服N －R 方法对于每次迭代需要重新形成并求逆一新的切线矩阵所带来的麻烦，常常可以采用一修正的方案、即mN －R 方法。其中切线矩阵总是采用它的初始值，即令，
0n T T K K = （1．2．10）．
因此（1·2·9）式可以修正为，
01
()()n n T a K P f -?=-+ （1．2．11），
这样一来，每次迭代求解的是一相同方程组。事实上，在用直接法求解此方程组时，系数矩阵只需要分解一次，每次迭代只进行一次回代即可。显然计算是比较经济的。虽然付出的代价是收敛速度较低，但总体上还是可能合算的。如和加速收敛的方法相结合，计算效率还可进一步改进。
另一种折衷方案是再迭代若干次（例如m 次）以后，更新T K 为m
T K ，再进行以后的迭代，在某些情况下，这种方案是很有效的。修正的N －R 方法的算法过程可表示如图1．3。
以上讨论的 N －R 法和 mN －R 法也隐含着K 可以显式地表示为a 的函数。而我们将讨论的弹塑性，蠕变等材料非线性问题，一般情况下由于应力依赖于变形的历史，这时将不能用形变理论，而必须用增量理论进行分析。在此情况下，不能将K 表示成a 的显式函数，因而也就不能直接用上述方法求解，而需要和以下讨论的增量方法，相结合进行求解。
图1.3 修正的 Newton-Raphson 方法
1．2．4增量法
为便于理解，假设方程（1．2．1）式表达的是结构应力分析问题，其中a 代表结构的位移，f -代表结构的载荷。所谓增量解法首先将载荷分为若干步：012,,f f f ，…，相应的位移也分为若干步：
012,,a a a ，…。每二步之间的增长量称为增量。增量解法的一般做法是假设第m 步载荷m f 和相应的位
移m a 已知，而后让载荷增加为1()m m m f f f +=+?，再求解1()m m m a a a +=+?。如果每步载荷增量m f ?足够小，则解的收敛性是可以保证的。同时，可以得到加载过程各个阶段的中间数值结果，便于研究结构位移和应力等随载荷变化的情况。
为了说明这种方法，将（1．2．1）式改写成如下形式：
0()()0a P a f ψλ=+= （1．2．12）
其中人是用以表示载荷变化的参数，上式对λ求导可以得到
000T dP da da
f K f da d d λλ
+≡+= （1．2．13）
从上式可以进一步得到
1
0()T da K a f d λ
-=- （1．2．14） 其中T K 即（1．2．8）式所定义的切线矩阵。
上式所提出的是一典型的常微分方程组问题，可以利用很多解法，最简单的是Euler 法，它可被表达为
11
10()()m m T m m T m m a a K a f K f λ--+-=-?=-?
（1．2．15）
其中：
m
m m m m m f f f -=?-=?++11λλλ
其他改进的积分方案（例如 Runge －Kutte 方法的各种预测校正）可以用来改进解的精度。和二阶 Runge －Kutte 方法等价的一种校正的 Euler 方法是可以采用的。即先按（1．2．15）式计算得到1m a +的
预测值，并表示为1m
a +'，再进一步计算1m a +；的改进值如下 1
1()m m T m m a a K f θ-++-=-? （1．2．16）
其中：
()()T m T m K K a θθ++≡
1(1)m m m a a a θθ++'=-+ (01)θ≤≤
利用上式计算得到的1m a +较利用（1．2．1）式得到的预测值1m
a +'将有所改进。 需要指出，无论是利用（1．2．15）式还是用（1．2．16）式计算1m a +或它的改进值，都是近似积分（1．2．14）式的结果，而未直接求解（1．2．12）式，因此所得到的1,m m a a +，…一般情况下是不能精确满足方程（1．2．12）式的，这将导致解的漂移。而且随着增量数目的增加，这种漂移现象将愈来愈严重。当系统为单自由度时，用Euler 法求解的增量方程（1．2．1）式以及解的漂移现象如图1．4所示。
为克服解的漂移现象，并改进其精度，可采用的方法之一，是从（1．2．16）式求得1m a +的改进值
图 1．4 Euler 法求解增量方程和解的漂移
以后，将它作新的预测值1m a +，仍用（1．2．16）式再计算新的改进值，为此继续迭代，直至方程（1．2．12）式在规定的误差范围内被满足为止。但每次迭代需要重新形成新的切线刚度()T m K a θ+。
现在更多采用的方法，是将N －R 方法或mN －R 方法用于每一增量步。如采用N －R 方法，在每一增量内进行迭代，则对于人的m ＋1次增量步的第n ＋1次迭代可以表示为
11111011011()()()n n n n n
m m m m m T m m P a f P a f K a ψλλ+++++++++≡+=++? （1．2．17）
由上式解出
))(()(011111f a P K a m n
m m n T n m ++-+++-=λ （1．2．18）
于是得到1m a +的第n ＋1次改进值
1111n n n
m m m a a a ++++=+? （1．2．19）
(1．2．17)式中1()n T m K +是1()T m K +的第n 次改进值。开始迭代时用0
1m m a a +=。连续地进行迭代，最后可以使得方程（1．2．12）式能够在规定误差范围内被满足。
从（1．2．17）式可见，当采用N －R 迭代时，每次迭代后也都要重新形成和分解1()n
T m K +，无疑工作量是很大的，因此通常采用mN －R 方法，这时0
11()()()n
T m T m T m K K K a ++==。
如果（1．2．18）式只求解一次，而不继续进行迭代，则有
01
1110()()m m T m m m a a K P f λ-+++?=?=-+ （1．2．20）
若进一步假设在上一增量步结束时，控制方程（1．2．12）式是精确满足的，即
00m m P f λ+=
则 1
10()m T m m a K f λ-+?=-? （1．2．21）
实际上，上式就是（1．2．15）式。而（1．2．20）式和上式相比不同之处在于它考虑了上一增量步中方程（ 1． 2． 12）式未精确满足的因素，将误差0m m P f λ+合并到0m f λ?中进行求解。（1．2．12）式在结构分析中实质上是平衡方程，所以（1. 2．18）式或（1．2．20）式称为考虑平衡校正的迭代算法。对于一个自由度的系统，将N －R 法或mN －R 法和增量法结合使用时，计算过程可以表示如图1．5。
图1．5 （a ）用N －R 法解增量方程 （b ）用mN －R 法解增量方程
1．2．5加速收敛的方法
由前面的讨论中已知，利用mN －R 方法求解非线性方程组时，可以避免每次迭代重新形成和求逆切线矩阵，但降低了收敛速度。特别是P －a 曲线突然趋于平坦的情况（对于结构分析问题，是结构趋于极限载荷或突然变软），收敛速度会很慢。为加速收敛速度，可以采用很多方法。这里介绍一种常用的、简单而有效的Aitken 加速法。
首先讨论单自由度系统，具体算法表示于图1．6。其中（a ）、（b ）分别为未采用和采用Aitken 加速的情况。
图1．6（a ）无Aitken 加速的mN －R 迭代 （b ）有Aitken 加速的mN －R 迭代
假设1m f +的初始试探解已知为0
1m m a a +=，利用修正的N －R 法进行迭代，求得第1，2次迭代后的改进解为
1111()(())n n
m T m m m a K P a f -+++?=-+ (0,1)n =
1111n n n
m m m a a a ++++=+? （1．2．22） 在求得面11m a +?以后，可以考虑寻求它的改进值1
1m a
+? 以加速收敛。Aitken 方法首先利用这两次迭代的不平衡差值来估计起始切线刚度()T m K 与局部割线刚度S K 的比值，从图1．6（b ）可见
001
111()()S m T m m m K a K a a +++?=?-?
所以
011
01
11()()
T m m S m m K a a K a a +++?==?-? （1．2．23） 然后以此比值来确定1
1m a
+? 即令 11
11()S m T m m K a K a ++?=?
这样就得到
11
11111()T m m m m S
K a a a a K +++?=
?=? （1．2．24） 从（1.2.23）式可知11a >,称1
a 为加速因子。并从（1．2．24）式得到1
1
11m m a
a ++?>? ，于是2
1m a +可以表示成
211111
11111m m m m m a a a a a a +++++=+?=+? （1．2．25）
从以上讨论可以推知,Aitken 加速收敛方法是每隔一次迭代进行一次加速,一般化地表示为
1111n n n n
m m m a a a a ++++=+? （1．2．26）
其中 1
1
11
11(0,2,...)(1,3,...)
n n m n n
m m n a a n a a -+-++=??=??=??-??
推广到N 个自由度的系统，Aitken 方法则以表示为
1111n n n n
m m m a a a α++++=+? （1．2．27）
其中n α是对角矩阵,它的元素n
i α是
1
,1
1,1
,11(0,2,...)(1,3,...)
n n i i m n n
i m i m n n αααα-+-++=?
?=??=??-?? （1．2．28）
从上式可见, 如果对于某个自由度i ,分母项 1
,1,1n n
i m i m αα-++?-?是很小的数值,这时 n
i α将是很大的
数值, 特别是当分母项趋于零, n
i α趋于无穷时,将使计算发生困难.为避免此情况,有提出了修正的Atiken 方法,即用一个标量代替(1.2.28)式中的对角矩阵.n
α这时
11
111
1111111(0,2,...)()(1,3,...)
()()
n n n T n m m m n n T n n m m m m n a a a a n a a a a --+++--++++=?
?=?-???=??-??-??
（1．2．29）
计算实践表明，当采用Aitken 法或修正的Aitken 法以后，收敛速度将有很大的改进。
应当指出，本节所讨论的上述算法只是目前常用的一些基本算法。实际上，由于有限元非线性分析是既费时又麻烦的工作，引起了研究工作者和实际分析工作者的广泛兴趣。旨在改进非线性分析的精度和效率的新的研究成果，不断出现在工程数值分析（International Journal for Numerical Methods In Engineering）、计算机与结构（Computers ＆Structures）等国际刊物，这里不—一列举。但是以上介绍的几种算法仍是当前发展阶段的最基本的方法。
另外，在众多的解法当中各有优缺点，很难说哪一种方法在任何情况下都比其它方法优越。但如果在计算程序中只准备编入一种算法，建议仍应采用具有N－R或mN－R迭代的增量法。只要增量步长足够小，收敛性是可保证的，一般情况下收敛速度也是令人满意的。
再有，在实际的计算中还可能遇到不少问题，例如增量步长的选择，结构刚度突然变化的处理等，将在以后结合具体问题进行讨论。
1．4弹塑性增量分析的有限元格式
1．4．1弹塑性问题的增量方程
由于材料和结构的弹塑性行为与加载以及变形的历史有关。在进行结构的弹塑性分析时，通常将载荷分成若干个增量，然后对于每一载荷增量，将弹塑性方程线性化，从而使弹塑性分析这一非线性问题分解为一系列线性问题。
假设对于时刻t 的载荷和位移条件的位移t
i u ，应变t
ij ε；和应力t
ij σ已经求得，当时间过渡到
t t +?(在静力分析且不考虑时间效应的情况t 和t t +?都只表示载荷的水平），载荷和位移条件有一增量，
即
)))t t
t
i i i t t t i i i t t t
i i i F F F V T T T u u u σ+?+?+??
=+??
?
=+???=+??
?
u (在内(在S 内(在S 内 （1．4．l ）
现在要求解t t +?时刻的位移，应变和应力
???
?????+=?+=?+=?+?+?+ij ij t ij t
t ij ij t
ij t t i
i t i t
t u u u σσσεεε （1．4．2）
它们应满足的方程和边界条件是 平衡方程
)(0,,内在V F F i i t
j ij j ij t
=?+++σσ （1．4．3）
应力和位移的关系
,,,,1
1()()()2
2
t
t t ij ij i j j i i j j i u u u u V εε+?=++?+?在内 （1．4．4）
应力和应变关系（为讨论方便起见，这里未考虑温度和蠕变的影响）
))((内在V t t t D kl ep
ijkl ij ?+≤≤?=?τεστ （1．4．5）
边界条件
()t
t
i i i i T T T T S σ+?=+?在上 （1．4．6）
()t
t
i i i i u u u u u S +?=+?在上 （1．4．7）
上式中
j
ij i j
ij t i t
n T n T σσ?=?= （1．4．8）
需要指出，在小变形的弹塑性分析中，除应力应变关系而外，其他方程和边界条件都是线性的。所以（1．4．3）～（l ．4. 8）式中除（1．4．5）式而外都未作进一步简化。如果,t
t
t
i ij ij u εσ和已精确地满足时刻t 的各个方程和边界条件，则可以从上列方程和边界条件中消去它们。现在仍保留它们是由于
按照数值求解的结果，它们不一定精确地满足方程和边界条件。这样做相当于进行一次迭代，可避免解的漂移。
至于应力应变关系表示成（1．4．5）式这是一种线性化处理。因为ij σ?应通过积分下列非线性关系
t t
t t
ep ij ij ijkl kl t
t
d D d σσε+?+??==?
?
（1．4．9）
得到，上式中ep
ijkl D 是,,p
ij ij σαε等的函数，而它们本身都是待求的未知量，所以将ij σ?表示为（1．4．5）式是一种线性化处理。如取
ep t ep
ijkl ijkl D D τ
=相当于最简单的Euler 法。在弹塑性有限元分析中称为起点切
线刚度法。当然还可有其他预测校正的方法以提高计算的精度和效率。
1．4．2增量有限元格式
首先建立增量形式的虚位移原理如下：如果t t ?+时刻的应力ij ij t
σσ?+和体积载荷i i t
F F ?+及边界载荷i i t
T T ?+满足平衡条件，则此力系在满足几何协调条件的虚位移)(i u ?δ［在V 内，
)(21)(,,i j j i ij u u ?+?=?δεδ；在u S 上，0)(=?i u δ］上的总虚功等干0。即
()()()()()()0
t
t
i i ij ij ij i V
V
t
i i i S dV F F u dV
T T u dS σ
σσδεδδ+??-+??-+??=?
?? （1．4．10）
将（1．4．5）式代入上式，则可得到
()()()()()()ep
i i ijkl kl ij i i V
V
S t
t
t
i i ij ij i i V
V
S D dV F u dV T u dS
dV F u dV T u dS
σ
σ
τ
εδεδδσδεδδ??-??-??=-?+?+??
?????
（1．4．11）
或表示成矩阵形式如下
????
???+?+?-=??-??-??σ
σ
δδσεδδδεεδτS t
T V
t
T v
t
T V
V
S T T ep T dS
T u dV F u dV dS
T u dV F u dV D )()
()()()()( (1．4． 12)
上式实际上就是增量形式的最小势能原理。它的左端和全量的最小势能原理的表达式在形式上完全相同。只是将全量改为增量。上式的右端是考虑,t
t
t
F T σ和可能不精确满足平衡而引入的校正项，也可理解为不平衡力势能（相差一负号）的变分。
基于增量形式虚位移原理有限元表达格式的建立步骤和一般全量形式的完全相同。首先将各单元内的位移增量表示成结点位移增量的插值形式
e u N a ?=? （1．4．13）
再利用几何关系，得到
e B a ε?=? （1．4．14）
将以上两式代入（1．4．12）式，并由虚位移的任意性，就得到有限元的系统平衡方程
ep K a Q τ
?=? （1．4．15）
其中,,ep K a Q τ
??分别是系统的弹塑性刚度矩阵，增量位移向量和不平衡力向量。它们分别由单元的各个对应量集成，即
,e
e
ep ep e e
t t t t t e e l i l i e e K K a a Q Q Q Q Q τ
τ+?+??=?=????=-=-?
?
∑∑∑∑ （1．4．16）
并且
e
e e
e T ep ep V t t
t t t t
e T
T l V S t
e T t i V K B D BdV
Q N FdV N TdS Q B dV
σ
τ
τσ+?+?+??
=??
=+???=?
???
? （1．4．17）
上式中
,t t
t l i Q Q +?分别代表外加载荷向量和内力向量，所以Q ?称不平衡力向量。如果,t t l i Q Q 满足平衡
的要求，则Q ?表示载荷增量向量。
从（1．4．15）式解出a ?以后，利用几何关系（1．4．14）式可以得到ε?，再按（1．4．9）式对本构关系进行积分可以得到σ?，并进而得到t t
t σσσ+?=+?。
需要指出，这时如将t t
σ+?代入（1．4．1）式右端，将发现
()e
e
t t
t t t t e t t e
l i l i e
e
t t
t t
T T T t t
V S V e
e
Q Q Q Q N
FdV N
TdV B
dV
σ
σ+?+?+?+?+?+?+?-=-=+-∑∑∑∑???
在一般情况下并不等于0，这表明此时求得的应力
t t
σ+?和外载荷t t
F +?及
t t
T +?尚未完全满足平衡条件，
也即仍存在不平衡力向量Q ?。需要通过迭代，以求得新的,t t
a εσσ+????和及，直至（1．4．15）
式的右端为0，即求得和外载荷完全满足平衡要求的应力状态为止。
从以上讨论可见，弹塑性增量有限元分析在将加载过程划分为若干增量步以后，对于每一增量步包含下列三个算法步骤：
1·线性化弹塑性本构关系（如（1．4．5）式），并形成增量有限元方程（如（1·4·1） 式）。
2·求解有限元方程（每个增量步或每次迭代的ep K τ
。都可能发生局部的变化人
3·积分本构方程（如（1..4·9）式）决定新的应力状态，检查平衡条件，并决定是否进行新的迭代。。 上述每一步骤的算法方案和数值方法，以及载荷增量步长的选择关系到整个求解过程。程的稳定性、精度和效率，是有限元数值方法研究的重要课题之一。以下将对其中的几个问题进行必要的讨论。
1·5数值方法中的几个问题
1．5．1非线性方程组的求解方案
利用工5． 2节所讨论的非线性方程组的一般解法于材料非线性有限元分析，根据具体问题的特点，可以组成不同的求解方案。通常采用的有下列几种：
1．无迭代的增量解法
此解法是对于每个载荷增量，求解方程（1．4．15）式。如进一步令t τ=，实际上就是（1．2．20）式所表述的Euler 算法。由于应力应变关系线性化带来的近似性，无疑地，如果要求得到足够精确的解答，必须采用足够小的载荷增量。另外，每一增量步中要重新形成和分解刚度矩阵，一般情况下计算效率是不高的。
2．具有变刚度迭代（N －R 迭代）的增量解法
此解法就是将（1．2．18）和（1．2．19）式所表述的算法用于弹塑性分析的情况。此时系统求解方程（1．4．1）式可以改写成
()()()t t
n n n ep K a Q +??=? （1．5．1）
其中n 是迭代次数，n ＝0，1，2，…，
()()
()()()()()()
(,,)e
e
t t
n T t t n ep ep V e
t t p n t t n t t n t t n ep ep n t t T t t n l V e K B D dV
D D Q Q B dV σαεσ+?+?+?+?+?+?+?+??=???=???=-???
∑?∑? （1．5．2） 且有
p
t p t
t t t t t t
t εε
αασσ===?+?+?+)
0(,
)0(,)0( （1．5．3）
当n=0时，方程（1．5．1）式就是t τ=的（1．4．15）式。具体的迭代步骤是
（1）利用（1．5．2）式计算
()
(),t t
n n ep K Q +??形成方程组（1．5．1）式。
（2）求解方程组（1．5．1）式，得到本次迭代的位移增量修正量()
n a
?
()()1()()n t t n n ep a K Q +?-?=? （1．5．4）
于是得到
(1)()()t t
n t t n n a a a +?++?=+? （1．5．5）
（3）计算各单元应变增量和应力增量修正量
()()n n B a ε?=? (1．5．6)
()()
n n ep D d εσ
ε'
??=?
（1·5·7）
其中()n ε
'
?是()
n ε
?中的弹塑性部分。关于()n ε
'
?的确定和（1．5．7）式的积分计算见以后1．5．3
节。从()
n σ
?可得
(1)()()t t
n t t n n σσσ+?++?=+? （1．5．8）
（4）根据收敛准则检验解是否满足收敛要求。如已满足，则认为此增量步内迭代已经收敛。对于每个增量步执行上述迭代，直至全部时间内的解被求得。
常用的收敛准则有 位移收敛准则：
a er a t D n ≤?)( （1． 5． 9）
平衡收敛准则： ()
(0)n F Q er Q ?≤? （1． 5．10）
能量收敛准则： ()()()(0)()(()n T
n n T E a Q er a Q ??≤?? （1． 5．11）
,,D F E er er er 是规定的容许误差。关于收敛准则的选择和容许误差的规定需要考虑具体问题的特点
和精度要求。
3．具有常刚度迭代（mN －R 迭代）的增量解法
解法步骤和前一种迭代解法相同，只是刚度矩阵保持某时刻的数值。如保持为每个增量步开始时刻t 的数值，称为起点切线刚度，这时
∑?===?+?+e
V ep t T ep t
ep t t n ep t
t e
BdV D B K K K )0()( （1． 5． 12）
其中 (,,)t p
t
t t
ep ep D D σαε=
如果刚度矩阵始终保持为弹性刚度矩阵,解法等效于一般的初应力法.这时
()t t
n T ep e e e
K K B D BdV +?==∑? （1． 5． 13）
对于后一种常刚度迭代，刚度矩阵在整个求解过程中只要形成和分解一次。对于前一种常刚度迭代
也不限制为起点切线刚度，例如先用前一增量步的刚度矩阵迭代l ～2次以后，再形成本增量步的刚度矩阵，以及可以规定经一定的迭代次数，仍不满足收敛准则的要求，则重新形成和分解新的刚度矩阵。
以上讨论了变刚度迭代和常刚度迭代的基本步骤，但是具体采用何种求解方案，仍应根据具体问题的特点，综合考虑精度和效率两方面因素。变刚度迭代具有良好的收敛性，允许采用较大的时间步长，但每次迭代都要重新形成和分解新的刚度矩阵。而采用常刚度迭代可以节省上述计算费用，缺点是收敛速度较慢，特别在接近载荷的极限状况时，因此经常需要同时采用加速迭代的措施。另方面在变刚度迭代中，时间步长常常也受到一些实际因素的限制，例如包含时间效应的问题（如蠕变问题、动力问题），如果解的方案是条件稳定的，时间步长必须小于某个规定的临界步长。再如对非线性本构关系线性化时，曾假定在一个时间步长内，应变向量方向不变，实际上该方向在时间步长内是变化的，对于该方向变化较快的结构和受力情况，上述假定将会带来较大的误差，因此时间步长必须限制在较小的范围内，这时采用变刚度迭代就不必要和不适合了。基于上述讨论，在一个通用的计算程序中应包含若干控制参数，以便根据具体问题的特点选择合理的求解方案。
1．5．2载荷增量步长的自动选择
在以上的讨论中都是假定载荷增量（或时间步长）是事先已经规定了的。这种事先将外加载荷分为若干个规定大小的增量步进行分析的方法是最简单的、也是常用的。但是在很多情况下，这种方法是不可行的。例如由理想弹塑性材料组成的结构，当载荷到达某一极限值时，结构将发生垮塌。表现在基于小变形理论的弹塑性分析中，当载荷到达极限值时，结构的位移将可无限增长。对应地，结构的载荷一位移曲线如图1.15所示。实际分析，由于极限载荷lim p 正是待求的未知量，如采用事先规定载荷增量步长的方法进行分析，当
lim t t
P P +?>时，将导致求解失败。因为在采用N －R 迭代求解时，将会遇到K 奇
异的情况（如图1.15（a ）所示）；而如采用mN －R 迭代求解，则总不能收敛（如图1.15（b ）所示）。为求得比较精确的极限载荷lim p ，此时可以采用的方法，是将载荷增量减小为二分之一，继续进行分析。如再遇上述求解失败情况，则再将载荷增量对半，再继续分析。这样逐步试探地逼近lim p 。无疑，这是一种麻烦而效率不高的方法。因此有必要研究载荷增量步长自动选择的方法。
图15·15规定载荷增量条件计算结构的极限载荷
（a ） N －R 迭代（b ） mN －R 迭代
在研究载荷增量步长自动选择的方法时，首先是假设载荷的分布模式是给定的，变化 的只是它的幅值。在此情况下，外载荷可表示成
00
00,,t
t
t t F pF T pT F pF T pT ?==?
??=??=??? （1． 5． 14）
相应地，结构的等效结点载荷向量也表示成
00,t
t
l l Q pQ Q pQ =?=? （1． 5． 15）
以上式中000,,F T Q 分别体积力、表面力和结点载荷的模式，p 是载荷幅值，或称载荷因子。载荷的分步实际就是()p t 的分步。
假定对于()t
p p t =的解已求得，现在要确定
()t t
t p p p t +?=+的大小。在载荷步长自动选择的求解
过程中，通常p ?不是一次给定，而是通过多次迭代不断修正以最后确定。现约定如下表示式：
(1)()()
(0,1,2,.)t t
n t t
n n p p p n r +?++?=
+?= (1．５．16)
且有
()0
r
n n p p =?=?∑ （1. 5．17）
其中n 是迭代次数，r 是迭代收敛的次数，()
n p ?是n 次迭代确定的p ?的修正值。
()t t
n p +?经1n -次迭
代修正后得到的
t t
p +?的数值。
现具体讨论几种常用的自动选择载荷步长的方法。 1．规定“本步刚度参数”的变化量以控制载荷增量。
此方法中，每一增量步的载荷因子增量p ?
仍是一次确定的，但是它是通过事先规定结构在本增量
步的刚度变化来控制的，这是Bergan 等人提出的一种自动进行载荷分布的方法。此法对计算结构的极限载荷特别有效。其基本思想是对于每一载荷分步，让结构刚度的变化保持差不多的大小。 （l ）“本步刚度参数”的概念
第i 增量步结构的总体刚度可用下式度量
)
()()()()*
(i T i i T i i p
Q
a Q Q S
????= （1． 5． 18） 初始（全弹性）结构总体刚度的度量是
(0)*T e e
p
T e e
Q Q S
a Q = （1． 5． 19） 其中e e Q a 和是载荷向量和按弹性分析得到的位移向量。 用无量纲参数)
(i p S 作为第i 步的结构刚度参数
()*()(0)*i p
i p
p
S S
S
=
（1． 5． 20）
()i p S 称为第i 步的“本步刚度参数”，它代表结构本身的刚度性质，与载荷增量的大小无关。当结构处于
完全弹性时，()1i p S =。随着载荷的增加，结构中的塑性区逐渐扩大，结构逐渐变软，()
i p S 也逐渐减小。当到达极限载荷时，()
0i p S =。
对于比例加载情况，如前所述，这时可记()
00,i e e i Q p Q Q
p Q =?=?，于是（1.5.20）式可以简化为
()
0()0
T
i i e p
i T e p a Q S
p a Q ?=? （1． 5． 21） 其中e p 是弹性极限载荷参数，01Q p =是时的结点载荷向量（载荷模式）。
利用本步刚度参数可以使步长调整得比较合理，并可减少总的增量步数，特别适合于计算理想弹塑性材料组成的结构的极限载荷等情况。
（2）增量步长的自动选择
以计算结构的极限载荷为例，利用“本步刚度参数”的规定变化量自动选择增量步长时，具体的算法步骤是
（i ）求弹性极限载荷参数e p 先施加任意的载荷0pQ ，假定结构为完全弹性求解，求出结构内的最大等效应力max σ
，
则
max
s e p p
σσ= （1． 5． 22） 其中0s σ是材料的初始屈服应力。
(ii)给定第一步载荷参数增量1p ?（例如取1/,e p p N N ?=的值可以事先给定），用11
e p p p =+?求解第一增量步。
（iii ）给定第二及以后各增量步的刚度参数变化的预测值p
S ? （它的大小决定步长的大小，例如可在0.05～0.2之间选择），并给定刚度的最小允许值min
p S （到达极限载荷时0=p S ，但在计算中如0=p S ，则结构刚度矩阵奇异，所以min
p S 应为接近于0的小的正数），则每步载荷参数的增量为
{
}),3,2(,~min )
1()2(m in
)1(1
???=--??=?----i S S S S S p p n p
n p p
i p p i i （1． 5． 23）
然后用1i i p p p -=+?；求解第i 载荷增量步。
在第i 增量步的解答求得以后，利用（1． 5． 20）式计算本步刚度参数()
i p S 和它的变化值
()()(1)i i i p p p S S S -?=-和预测值p
S ? 会有一定差别。所以此算法是在保持本步刚度参数变化值接近为常数条件下选择载荷增量的大小。算法的执行示意如图 1.16。从图可见，i p ?是在给定p
S
? 情况下，利用2
11/i i i p
p
p S
S
---?-线性外推得到。实际计算中，由于~p S p 曲线不一定如图示那样规则，在给定p
S ? 的情况下，可能使i p ?时而过大、时而过小，影响计算的执行。此时可以利用(3)
(2)(1)
,,i i i p p p S S S ---二次外推
得到i p ?；，实际计算表明，这样确定i p ?可避免上述现象，从而使最终计算结果有较大改进。
图15·16载荷的自动分步
2．规定某个结点的位移增量以确定载荷增量
增量迭代的有限元求解方程在1. 5 . 1节中已给出。以mN －R 迭代为例，它可以表示成
()
(
)
(0,1,2,)n n ep K a
Q n τ
?=?= （1． 5． 24）
在载荷增量步长自动控制的求解方法中,()
n Q
?可以表示成
()(1)()n t t
n t t n l i Q Q Q +?++??=
- （1． 5． 25）
其中：
(1)()()()()0()t t
n t t
n n t t n n l l l Q Q Q p p Q +?++?+?=
+?=+? （1． 5． 26）
()()e
t t
n T t t n i V e
Q B dV σ+?+?=∑? （1． 5． 27）
（1. 5．25）式还可以改写成
()()()()
()0n n n n n u l u Q Q Q Q p Q ?=?+?=?+? （1． 5． 28）
其中()n u Q ?是n 一1次迭代后得到的()t t
n σ+?和外载荷()0t t n p Q +?。构成的不平衡结点力向量，
即
()
()()0n t t
n t t n u i Q p Q Q +?+??=
- （1． 5． 29）
()n p ?是在n 次迭代中由某个规定的约束条件来确定的载荷因子增量p ?的第n 次修正量。在现在的方法
中，这约束条件就是某个结点的位移增量的大小。例如规定a ?中的某个分量g a ?
是给定的，此条件可表示为
()()
(0)0
(1,2,)
n T n g
l n a
b a
n ?=??=?=?=? （1． 5． 30）
其中l ?是g a ?的规定值，b 是除第g 个元素为1，其余元素为0的向量。具体迭代的算法步骤如下：
（1）计算对于结点载荷模式0Q 的位移模式0a 。
01
0Q K a ep -=τ （1． 5． 31）
（2）计算对于不平衡结点力向量()n u Q ?的位移增量修正值()
n u a ?和n 次迭代后位移增量修改值的全量()
n a
?
()1()
n n u ep u a K Q τ-?=? （1． 5． 32）
()()
()0
(1,2,)n n n u a a p a n ?=?+?= （1． 5． 33）
其中()
n p
?是待定的载荷因子增量的修正值。
（3）利用条件（1． 5． 30）式确定()
n p
?
()
()
()0(0)0(1,2,)T n T
n n T
g u
l
n a b a
b a
p b a n ?=??=?=?+?=?=?
（1． 5． 34）
从上式可得
()()
()0
n T n g u
n T a b a p b a ?-??=
（1． 5． 35）
这样就确定了()
n p
?,从而得到()
()()0,n n n l a
Q p Q ??=?等。
（4）计算()
()(1)
n n n u Q ε
σ+???等，并检验收敛准则的要求是否满足。如未满足，回到步骤（2）进行
新的一次选代，直至收敛准则的要求满足为止。本增量步的外载荷增量、位移增量、应力增量的结果是
()
()000
()()00,,r r
n n l l
n n r r
n n n n Q Q
p Q a a σσ====?
?=?=??
???
?=??=???∑∑∑∑ （1． 5． 36） 其中r 为迭代收敛时的次数。
关于每一个增量步某个指定结点位移增量l ?本身的选择，通常的方法是第1个增量步可由某个给定的载荷因子增量1p ?（例如令N p p e =?1，其中e p 是弹性极限载荷因子，N 可取5～10）通过求解得到1l ?。以后各增量步的i l ?；可由下式确定，
11
--?=
?i i i l r r l （1． 5． 37） 式中1i l -?；是前一增量步的规定位移增量，1i r -是前一次增量步迭代收敛的次数，0r 是优化的迭代次数，例如取0r =4～6。上式的倾向是使在严重非线性阶段、l ?缩短，而在接近线性阶段、l ?加长。
实际计算表明，采用规定某个结点的位移增量以确定载荷增量步长的方法，结合采用（1. 5．37）式调节位移增量的大小，计算结构的极限载荷取得了相当满意的结果。
1. 5．3弹塑性状态的决定和本构关系的积分
从上节的讨论中可以看到，增量形式有限元格式中()
n ep D Q
τ?和的确定是基于已经求得上一增量步
或迭代结束时的σt ，αt 和p
ε，因此为了进行下一增量步或迭代的计算，需要根据本步或迭代计算得到的位移增量，决定,,p
σαε???等，从而得到本步或迭代结束时的弹塑性状态，即,,p
σαε等。这一步骤称为状态决定，它和增量方程的建立（线性化步骤）以及方程组求解一起构成了非线性分析的基本算法过程。它不仅在计算中占相当大的工作量，而且对整个计算结果的精度有很大影响，因此应当给予足够的重视。
1．决定弹塑性状态的一般算法步骤
在每一增量步或每次选代，求得位移增量或其修正量u ?以后，决定新的弹塑性状态的基本步骤如下：
（1）利用几何关系计算应变增量（或其修正量）
B a ε?=?
（2）按弹性关系计算应力增量的预测值以及应力的预测值
e D σ
ε?=? （1． 5． 38） t t
t σ
σσ+?=+? （1． 5． 39） 其中t
σ是上一增量步（即i —1步）结束时的应力值。
(3)按单元内各个积分点计算D 的预测值 1）计算屈服函数值),,~(p
t
t t
t F εασ
?+，然后区分三种情况 （i ）若0),,~(
≤?+p t
t t
t F εασ
，则该积分点为弹性加载，或由塑性按弹性卸载，这时均有 σσ
?=? （1． 5． 40） （ii ）若0),,~(0),,~(<>?+p t t t p t t t t F F εασεασ
且，则该积分点为由弹性进入塑性的过渡情况，应由
0),,~(=?+p t t t m F εασ
σ （1． 5． 41） 计算比例因子m 。上式隐含着假设在增量过程中应变成比例的变化。计算m 是为确定应力到达屈服面的时刻。采用V ．Mises 屈服准则时，m 是下列二次方程的解
22100a m a m a ++= （1． 5． 42）
其中： ),,(,~
)(,~~2
1012p t t t T t t T F a S S a S S a εασα=?-=??=
（1． 5． 43） 对于各向同性硬化情况，0≡α。因为σt
总是在屈服曲面上或屈服曲面之内，所以常数020.0,a a m ≤>同时必须取正值，所以
12(/2m a a =-+ （1． 5． 44）
（iii ）若0),,(0),,~(
=>?+p t t t p
t
t t
t F F εασεασ
且，则该积分点为塑性继续加载，这时令m ＝0。 2）对于（ii ）、（iii ）两种情况，均有对应于弹塑性部分的应变增量
(1)m εε'?=-? （1． 5． 45）
3）计算弹塑性部分应力增量
(,,)p
ep D d εσσαεε'
?'?=?
（1． 5． 46）
一般情况下用数值积分方法进行此积分，在积分过程中可以同时得到,p
αε??。 4）计算本增步或迭代结束时刻的
,,
t t
p
t t
t t σαε+?+?+?
,,t t
p t p p
t t
t t t t m σσσσαααεεε
+?+?+?'=+?+?=+?=+? （1． 5． 47）
2．本构关系的积分
关于本构关系的积分（1. 5．46）式，对于其中的某些情况，已经获得了解析解或解析的近似解，例如Key 得到了理想弹塑性本构关系的解析解。运动硬化材料本构关系的解析解和各向同性硬化材料本构关系的以/(3)p
p
E G E +为小参数的渐近解，近年来都已获得。它们用于实际计算取得较好的结果，即以较小的工作量得到较精确的结果。但是现有的计算程序中仍都采用数值积分方法，其中采用最多的是切向预测径向返回的子增量法。另外，广义中点法近年来也受到较多的重视。现对这二种方法作一较详细的介绍。
（1）切向预测径向返回子增量法
所谓切向预测就是将Euler 方法用于（1. 5．46）式，得到应力增量的预测值，即
(,,)t t t p ep D E σ
σαε?=? （1． 5． 48） 进一步得到应力的预测值
σσσ
~~?+=?+t t
t （1． 5． 49） 同时还可以得到
2
3
p s ε
λσ?=? （1． 5． 50） 2
(1)()(Pr )3p t t E M s ager αλα?=-?- 法则 （1． 5． 51） 2
(1)()()3
p t t E M Zeigler αλσα?=-?- 法则 （1． 5． 52）
其中
200(/)(/)(/)(4/9)(,)()T e T p e s t p t p
s s s s s f D f D f E M M σε
λσσσσσεσσεσ?????=?
?
????+????==+-??????
（1． 5． 53）
并且有
??
?
???+=?+=?+?+ααα
εεε~~~t t
t p
p t t t p （1． 5． 54）
以上各式是对于混合硬化的一般情况给出的，对于各向同性硬化，M=1；对于运动硬化，M=0。
因为（1. 5.50）式所表达的算法是显式的Euler 方法，其中的(,,)t
p
t t
ep D σαε是起点切线刚度，所
以σ
? 是在加载曲面的切线方向。同时由于加载曲面是外凸的，因此t t
σ
+? 总是在加载曲面之外。但是屈服准则要求应力
t t
σ+?只能在加载曲面之上或者之内，所以常需再采用径向返回的方法以求得满足屈
服条件的t t σ+?和t t
α+?。具体做法是令
非线性有限元方法及实例分析
非线性有限元方法及实例分析 梁军 河海大学水利水电工程学院，南京（210098） 摘 要：对在地下工程稳定性分析中常用的非线性方程组的求解方法进行研究，讨论了非线性计算的迭代收敛准则，并利用非线性有限元方法分析了一个钢棒单轴拉伸的实例。 关键词：非线性有限元，方程组求解，实例分析 1引 言 有限单元法已成为一种强有力的数值解法来解决工程中遇到的大量问题，其应用范围从固体到流体，从静力到动力，从力学问题到非力学问题。有限元的线性分析已经设计工具被广泛采用。但对于绝大多数水利工程中遇到的实际问题如地下洞室等，将其作为非线性问题加以考虑更符合实际情况。根据产生非线性的原因，非线性问题主要有3种类型[1]： 1．材料非线性问题（简称材料非线性或物理非线性） 2．几何非线性问题 3．接触非线性问题（简称接触非线性或边界非线性） 2 非线性方程组的求解 在非线性力学中，无论是哪一类非线性问题，经过有限元离散后，它们都归结为求解一个非线性代数方程组[2]： ()()()00 021212211=… …==n n n n δδδψδδδψδδδψΛΛΛ （1.1） 其中n δδδ,,,21Λ是未知量，n ψψψ,,,21Λ是n δδδ,,,21Λ的非线性函数，引用矢量记 号 []T n δδδδΛ21＝ （1.2） []T n ψψψψΛ21= （1.3） 上述方程组（1.1）可表示为 ()0=δψ （1.4） 可以将它改写为 ()()()0=?≡?≡R K R F δδδδψ （1.5） 其中()δK 是一个的矩阵，其元素 是矢量的函数，n n ×ij k R 为已知矢量。在位移有限 元中，δ代表未知的结点位移，()δF 是等效结点力，R 为等效结点荷载，方程()0=δψ表示结点平衡方程。 在线弹性有限元中，线性方程组
浅淡钢筋混凝土结构的非线性有限元
价值工程 0引言 钢筋混凝土结构是目前使用最为广泛的一种结构形式。钢筋混凝土是由两种性质不同的材料组合而成的，材料性能非常复杂，特别是在其非线性阶段，混凝土和钢筋本身的各种非线性特性，都不 同程度地在这种组合材料中反映出来。 传统的分析和设计方法往往采用线弹性理论来分析其内力。随着有限元理论和计算机技术的进步，钢筋混凝土非线性有限元分析方法也得以迅速的发展并发挥出巨大的作用。 1钢筋混凝土有限元分析原理钢筋混凝土有限元分析，主要是研究钢筋混凝土结构的基本性能、设计方法和构造措施。结合钢筋混凝土的力学特性，采用有限元分析的一般原理，是有限元分析和钢筋混凝土力学特性两者的结合。 Ngo 和Scordelis 在早期进行的研究中， 把有限元方法用于钢筋混凝土结构分析，它包含了钢筋混凝土有限元分析的基本原理。可以具体阐述为如下几点： 1.1确定各单元的单元刚度矩阵， 它与一般的有限元方法基本相同，并组合成结构的整体刚度矩阵。随着荷载和作用的不断增加，可以得到钢筋混凝土结构自开始受荷到破坏的整个过程的位移、应变、应力、裂缝的形成和发展、钢筋和混凝土结合面的粘结滑移、钢筋的屈服和强化以及混凝土压碎破坏等大量有用的数据，为研究结构的性能和合理的设计方法提供可靠的依据。根据结构所受的荷载和约束，解出节点的未知位移，进而求出单元的应力。 1.2确定适用于各类单元的本构关系， 这种关系可以是线性的，也可以是非线性的。即应力应变关系，或结点力位移关系。 1.3通过设置联结单元， 模拟裂缝两侧的混凝土之间的咬合作用，以及钢筋和混凝土之间的粘结滑移关系。 1.4把钢筋混凝土结构分割成有限个小的结构单元。这些单元可以是钢筋和混凝土的组合单元或分离式单元。 2钢筋混凝土的非线性有限元分析 2.1混凝土的破坏准则混凝土的破坏准则就是描述混凝土破坏时其应力状态或应变状态满足的条件。 根据混凝土破坏准则的函数f （ξ，r ，θ，k 1，k 2，k 3，……，k n ）=0中包含参数的个数，破坏准则可以分为单参数破坏准则、两参数破坏准则等等。单参数破坏准则有最大拉应力准则、最大剪应力准则及八面体剪应力准则。两参数破坏准则有Mohr -Coulomb 准则和 Drucker-prager 准则。 单参数和双参数都是早期提出的破坏准则。单参数或双参数的破坏准则不能全面反映混凝土的破坏特性。多参数破坏准则是适用性更广泛的破坏准则。它克服了单参数和双参数的一些不足，一些多参数破坏准则已能较好地描述混凝土的破坏特性。其中比较有代表性的二维的破坏准则有Kupfer-Gerstle 准则、 Hsieh-Ting-Chen 准则、李～过准则等。三维破坏准则有：Ottosen 准 则、Willam-Warnke 准则、 过－王、江－周准则等。2.2混凝土的本构模型混凝土的本构关系就是指混凝土的应力状态和应变状态的关系。目前，混凝土的本构模型主要类型有：以弹性模型为基础的线弹性和非线弹性的本构关系；以经典塑性理论 为基础的理想弹塑性和弹塑性硬化本构模型；采用断裂理论和塑性 理论组合的塑性断裂理论,并考虑用应变空间建立的本构模型；以粘性材料本构关系发展起来的内时程理论描述的混凝土本构模型；用损伤理论和弹塑性损伤断裂混合建立的本构模型等。 线弹性模型是工程上一般材料所采用的关系模型，线弹性类本构模型也是最简单、最基本的材料本构模型。材料变形在加载和卸载时都沿同一直线变化，完全卸载后无残余变形。因而，应力和应变有确定的一一对应的关系。直线的斜率为材料的弹性模量。如果混 凝土在单向受拉、 单向受压或多轴应力作用下，其应力-应变之间关系为曲线而非直线时，从原则线弹性模性已不适用。但在一些特定的情况下仍可使用线弹性模型，这样作的好处就是给分析带来方便、 快捷。非线性本构模型是能够比较正确模拟混凝土材料性质的本构模型，主要有非线性弹性本构模型和弹塑性本构模型。如Kupfer-Gerstle 的各项同性的全量模型、Darwin 正交异性增量模型和Ottosen 模型等。非线性弹性本构的优点是能反映混凝土受力变形的主要特点；计算公式和参数值都来自试验数据的回归分析，在单调比例加载的情况下有很高的计算精度；模型的表达式简明、直 观，易于理解和应用。因而， 这种模型在工程中应用最广。但它也有的缺点：不能反映卸载和加载的区别，卸载后没有残余变形等，故不能应用于加、卸载循环和非比例加载等情况。 2.3钢筋与混凝土之间的关系模型钢筋混凝土中钢筋和混凝土之间存在粘结力、骨料咬合力和销栓作用等，如何正确模拟钢筋和混凝土之间的相互作用，关系到有限元分析结果能否正确反映结构真实受力状态的关键。 钢筋与混凝土界面的有限元分析模型，根据是否考虑钢筋与混凝土之间的粘结滑移及销栓作用，以及用什么方式模拟这种作用，有两种基本不同的联结模型，一种是钢筋和混凝土之间位移完全协调的联结模式，另一种是两者之间位移不协调的连接模型，即采用粘结单元的联结模型。位移完全协调的联结模式，又分为分离式、埋置式和组合式三种模型。这些模式都认为钢筋和混凝土之间即无相对滑移，也无相对错动，不需要粘结滑移及销栓作用的模拟。粘结单元的联结模采用在钢筋单元和混凝土单元之间，设置粘结单元模拟两者之间的粘结力及销栓作用。在混凝土与钢筋之间的粘结模拟方面，人们提出了各种不相同的粘结单元的模型，比如无厚度四节点或六节点粘结单元、双弹簧粘结单元、、斜弹簧单元粘和结斜杆单元等。而关于粘结～滑移关系方面，在分析初期采用的是线性关系，随后发展为非线性关系，提出多种τ～S 曲线的表达式。因为存在的影响因素比较多，而且问题相对复杂，所以目前尚且还没有相对完善的计算模式。 2.4裂缝的模拟混凝土受拉开裂后形成裂缝，在钢筋混凝土的有限单元法中，裂缝的模型很多，一般比较常用的是单元边界的的单独裂缝和单元内部的弥散裂缝以及断裂力学模型这三种模型。第一种方法把裂缝处理为单元边界，一旦出现新的裂缝就增加新的节点，重新划分单元，使裂缝总是处于单元和单元之间的边界。这种方法的缺点是计算工作繁琐，费机时。第二种方法使得在计算过程中裂缝自动形成和发展，即不必增加结点也不用重新划分单元，所以由计算机自动进行处理比较容易，因而得到了较为广泛的应用。 —————————————————————— —作者简介：范治华（1980-），男，河南永城人，助理工程师，研究方向为城建。 浅淡钢筋混凝土结构的非线性有限元分析 Nonlinear Finite Element Analysis of Reinforced Concrete Structure 范治华Fan Zhihua ；史玉侠Shi Yuxia （神火集团有限公司，永城476600） （Shenhuo Group Company ，Yongcheng 476600，China ）摘要：随着有限元理论和计算机技术的进步，钢筋混凝土非线性有限元分析方法必定会在理论实践和工程实施中得到更大程度的发展，发 挥更加强大的作用。 Abstract:With the progress of finite element theory and computer technology,nonlinear finite element analysis of reinforced concrete must be implemented in the theory and engineering practice and get the greater degree of development,and play a more powerful role. 关键词：钢筋混凝土结构；有限元；分析Key words:reinforced concrete structures ；finite element ；analysis 中图分类号：TU37 文献标识码：A 文章编号：1006-4311（2012）05-0086-02 ·86·
非线性有限元分析
轨道结构的非线性有限元分析 姜建华 练松良 摘 要 实际轨道结构受载时的力学行为,属于典型的非线性力学问题。钢轨垫层刚度、钢轨抗扭刚度和扣件扣压力的大小是影响轨距扩大的主要因素。根据非线性有限元接触理论,建立了能准确反映扣件、钢轨与垫层的拧紧接触,以及受载车轮与钢轨侧向滑动接触的力学计算模型;并研究计算了不同扣件压力下,由于受载车轮与钢轨侧向滑动接触引起的轨距扩大问题。 关键词 轮轨关系,扣件压力,非线性弹性力学,有限元分析 1 引言 实际工程中常见的非线性问题一般可以归纳为三类:材料非线性、几何非线性以及边界条件非线性。材料非线性问题是由于材料的非线性本构关系所引起的,例如材料的弹塑性变形,材料的屈服和硬化等;几何非线性问题是由于结构的位移或变形相当大,以至必须按照变形后的几何位置来建立平衡方程;边界条件非线性问题是指边界条件随位移变化所引起的非线性问题。通常情况下,我们所遇到的非线性问题多数是上述三类非线性问题的组合[1,2]。 实际轨道结构受载时的力学行为,属于典型的非线性力学问题。比如基于轮轨接触的材料非线性、几何非线性及边界条件非线性问题,以及扣件、钢轨、垫层三者间相互作用时所表现的边界条件非线性行为等。所以,机车车辆在轨道结构上行驶时引起的力学现象是相当复杂的。以往在研究轨道各部分应力应变分布规律时,通常采用连续弹性基础梁理论或连续点支承,偶尔简单考虑扣件的作用和弹性垫层的使用。不管用哪一种支承方式建立模型,都由于这样那样的假设而带有一定程度的近似性。所以,如何利用现代力学理论的最新成果以及日益发展的计算机技术,根据轨道结构的具体情况,建立更为完整更为准确的轨道结构计算模型,为轨道设计部门提供更加可靠的设计依据或研究思路,已十分必要。 本文提出了用非线性有限元理论研究轮轨系统和轨道结构的思路。作为算例之一,本文将根据非线性有限元理论,建立能准确反映扣件、钢轨与垫层的拧紧接触,以及受载车轮与钢轨侧向滑动接触的力学计算模型。 2 轨道结构的有限元接触模型 对于非线性问题,不管是材料非线性、几何非线性,还是边界条件非线性,总是最终归结为求解一组非线性平衡方程及其控制方程。例如用位移作为未知数进行有限元分析时,最后可得到一组平衡方程及其控制方程为 : 图1 轮轨系统的对称性模型简图 [K(u)]{u}={R}(1) (u)= (u)(2)其中:{u}为节点位移列阵;{R}为节点载荷列阵; [K(u)]为总体刚度矩阵; (u)为边界条件。它们 36 姜建华:同济大学工程力学系,副教授、博士,上海200092
第18章 接触问题有限元分析技术
第18章接触问题的有限元分析技术 第1节基本知识 接触问题是一种高度非线性行为，需要较大的计算资源，为了进行准确而有效的计算，理解问题的特性和建立合理的模型是很重要的。 接触问题存在两个较大的难点：其一，在求解问题之前，不知道接触区域，表面之间是接触或分开是未知的、突然变化的，这些随载荷、材料、边界条件和其它因素而定；其二，大多数的接触问题需要计算摩擦，有几种摩擦和模型可供挑选，它们都是非线性的，摩擦使问题的收敛性变得困难。 一、接触问题分类 接触问题分为两种基本类型：刚体─柔体的接触和半柔体─柔体的接触。在刚体─柔体的接触问题中，接触面的一个或多个被当作刚体，（与它接触的变形体相比，有大得多的刚度），一般情况下，一种软材料和一种硬材料接触时，问题可以被假定为刚体─柔体的接触，许多金属成形问题归为此类接触；另一类，柔体─柔体的接触，是一种更普遍的类型，在这种情况下，两个接触体都是变形体（有近似的刚度）。 ANSYS支持三种接触方式：点─点、点─面和平面─面。每种接触方式使用的接触单元适用于某类问题。 二、接触单元 为了给接触问题建模，首先必须认识到模型中的哪些部分可能会相互接触，如果相互作用的其中之一是一点，模型的对立应组元是一个节点。如果相互作用的其中之一是一个面，模型的对应组元是单元，例如梁单元，壳单元或实体单元。有限元模型通过指定的接触单元来识别可能的接触匹对，接触单元是覆盖在分析模型接触面之上的一层单元。下面分类详述ANSYS使用的接触单元和使用它们的过程。 1．点─点接触单元 点─点接触单元主要用于模拟点─点的接触行为，为了使用点─点的接触单元，需要预先知道接触位置，这类接触问题只能适用于接触面之间有较小相对滑动的情况（即使在几何非线性情况下）。 如果两个面上的节点一一对应，相对滑动又以忽略不计，两个面挠度（转动）保持小量，那么可以用点─点的接触单元来求解面─面的接触问题，过盈装配问题是一个用点─点的接触单元来模拟面─与的接触问题的典型例子。 2．点─面接触单元 点─面接触单元主要用于给点─面的接触行为建模，例如两根梁的相互接触。 如果通过一组节点来定义接触面，生成多个单元，那么可以通过点─面的接触单元来模拟面─面的接触问题，面即可以是刚性体也可以是柔性体，这类接触问题的一个典型例子是
非线性有限元读书报告
非线性有限元读书报告 1 对张量分析及其在结构非线性力学和有限元中应用的认识。 张量分析将张量中的每一项的计算过程（这些计算过程可能是相似的）用张量之间的计算方法代替，并统一，提供了一种更加方便的计算方法。张量标记中，指标不出现，因此张量标记的表示是独立于坐标系统的，并且可以用于笛卡尔坐标系，柱坐标系，曲线坐标等。如同矩阵一般，或者说张量就是矩阵的更一般的形式，而矩阵是二维的张量。此外，在张量标记中的方程非常容易记忆。 2 对各种应变、应力张量和材料本构关系的认识。 本段回答除本课程讲义外大部分参考文献[1] 2.1 应变张量 为了使刚体运动时，特别是刚体转动时，应变度量为0，工程应变张量被抛弃了。根据采用连续介质力学中的欧拉（Eulerian ）方法或拉格朗日方法（Lagrangian ），得到的应变分别为阿尔曼西（Almansi ）应变（或称欧拉应变） 张量（E ij ε）与格林（Green ）应变（或称拉格朗日应变）张量（G ij ε），它们与工 程应变（或称柯西小应变张量）（此处记为ij ε）均不相同。 三者的几何方程表达式为： 阿尔曼西（Almansi ）应变（或欧拉应变）张量： ()11,,,,11 ()22E ij ij ki kj i j j i k i k j F F u u u u εδ--= -=+-； 格林（Green ）应变（或拉格朗日应变）张量： (),,,,11 ()22G ij ki kj ij i j j i k i k j F F u u u u εδ= -=++； 工程应变（柯西小应变）张量： ,,1()2 ij i j j i u u ε=+。 其中u 为位移，ij F 为变形梯度矩阵，即i ij j x F X ??? =????? 另外，阿尔曼西（Almansi ）应变（或欧拉应变）张量（E ij ε）和格林（Green ）应变（或拉格朗日应变）张量（G ij ε）之间的转换关系为：
非线性有限元课程报告
1 非线性问题的求解方法 无论是哪一类非线性问题，经过有限元离散后,它们都归结为求解一个非线性代数方程组： ()0......,211=n δδδψ ()0......,212=n δδδψ …… ()0......,21=n n δδδψ 其中：1δ,2δ……n δ是未知量，1ψ,2ψ……n ψ是1δ,2δ……n δ的非线性函数，现引用矢量记号 []T n δδδδ (21) []T n ψψψψ (21) 上述方程组可表示为： ()0=δψ 还可以将它改写为： ()()()0=-=-=R K R F δδδδψ ()δK 是一个n n ?的矩阵，其元素ij k 是矢量δ的函数，R 为已知矢量。在位移有限元中，δ代表未知的结点位移，()δF 是等效结点力，R 为等效结点荷载，方程()0=δψ表示结点的平衡方程。 在线弹性有限元中，线性代数方程组0=-R K δ可以毫无困难地求解，但对非线性方程组()0=δψ则 不行。一般来说，难以求得其精确解，通常采用数值解法，把非线性问题转化为一系列线性问题。为了使这一系列线性解收敛于非线性解，曾经有过许多方法，但这些解法都有一定的局限性。某一解法对某一类非线性问题有效，但对另一类问题可能不合适。因而，根据问题性质正确选用求解方法成为非线性有限元的一个极重要的问题。 1.1 直接迭代法 目前求解非线性方程组的方法一般为线性化方法。若对总荷载进行线性化处理，则称为迭代法。 对非线性方程组： ()0=-R K δδ 1-1 设其初始的近似解为0 δδ=，由此确定近似的K 矩阵()0 δ K K =，根据式1-1可得出改进的近似解 ()R K 1 01-=δ。重复这一过程，以第i 次近似解求出第1+i 次近似解的迭代公式为： ()i K K i δ= ()R K i i 1 1 -+=δ 1-2 直到Δi i i δδ δ-=+1 变得充分小，即近似解收敛时，终止迭代。在迭代过程中，得到的近似解一般不会满足式1-1，即()()0≠-=R K i i i δδδψ。()δψ作为对平衡偏离的一种度量，称为失衡力。 对于一个单变量问题的非线性方程，δ-F 为凸曲线时，直接迭代法的计算过程如图1-1所示，可以看出()δK 就是过曲线上点()()δδF ,与原点的割线斜率。对于单变量问题，这一迭代过程是收敛的，但对多自由度情况，由于未知量通过矩阵耦合，迭代过程可能不收敛。同样可以得出δ-F 为凹曲线时的情况。
非线性有限元分析
非线性有限元分析 1 概述 在科学技术领域，对于许多力学问题和物理问题，人们已经得到了它们所应遵循的基本方程（常微分方程或偏微分方程）和相应的定解条件（边界条件）。但能够用解析方法求出精确解的只是少数方程性质比较简单，并且几何形状相当规则的问题。对于大多数工程实际问题，由于方程的某些特征的非线性性质，或由于求解区域的几何形状比较复杂，则不能得到解析的答案。这类问题的解决通常有两种途径。一是引入简化假设，将方程和几何边界简化为能够处理的情况，从而得到问题在简化状态下的解答。但是这种方法只是在有限的情况下是可行的，因为过多的简化可能导致误差很大甚至是错误的解答。因此人们多年来一直在致力于寻找和发展另一种求解途径和方法——数值解法。特别是五十多年来，随着电子计算机的飞速发展和广泛应用，数值分析方法已成为求解科学技术问题的主要工具。 已经发展的数值分析方法可以分为两大类。一类以有限差分法为代表，主要特点是直接求解基本方程和相应定解条件的近似解。其具体解法是将求解区域划分为网格，然后在网格的结点上用差分方程来近似微分方程，当采用较多结点时，近似解的精度可以得到改善。但是当用于求解几何形状复杂的问题时，有限差分法的精度将降低，甚至发生困难。 另一类数值分析方法是首先建立和原问题基本方程及相应定解条件相等效的积分提法，然后再建立近似解法并求解。如果原问题的方程具有某些特定的性质，则它的等效积分提法可以归结为某个泛函的变分，相应的近似解法实际上就是求解泛函的驻值问题。诸如里兹法，配点法，最小二乘法，伽辽金法，力矩法等都属于这一类方法。但此类方法也只能局限于几何形状规则的问题，原因在于它们都是在整个求解区域上假设近似函数，因此，对于几何形状复杂的问题，不可能建立合乎要求的近似函数。 1960年，R.W.CLOUGH发表了有限单元法的第一篇文献“The Finite Element Method in Plane Stress Analysis”，这同时也标志着有限单元法（FEM）的问世。有限单元法的基本思想是将连续的求解区域离散为一组有限个，且按一定方式相互联接在一起的单元的组合体。由于单元能按不同的联结方式进行组合，且单元本身又可以有不同形状，因此可以模型化几何形状复杂的求解域。并且可以利用在每一个单元假设的近似函数来分片地表示全求解域上待求的未知场函数，从而使一个连续的无限自由度问题变成离散的有限自由度问题。 现已证明，有限单元法是基于变分原理的里兹法的另一种形式，从而使里兹法分析的所有理论基础都适用于有限单元法，确认了有限单元法是处理连续介质问题的一种普遍方法。利用变分原理建立有限元方程和经典里兹法的主要区别是有限单元法假设的近似函数不是在全求解域而是在单元上规定的，而且事先不要求满足任何边界条件，因此可以用来处理很复杂的连续介质问题。 在短短四十余年的时间里，有限单元的分析方法已经迅速地发展为适合于使用各种类型计算机解决复杂工程问题的一种相当普及的方法。如今，有限元广泛地应用于各个学科门类，已经成为工程师和科研人员用于解决实际工程问题，进行科学研究不可或缺的有力工具。有限单元法的应用围已由弹性力学平面问题扩展到空间问题，板壳问题，由静力平衡问题扩展到稳定问题，动力问题和波动问题。分析的对象从弹性材料扩展到塑性，粘弹性，粘塑性和复合材料等，从固体
过盈配合应力的接触非线性有限元分析
过盈配合应力的接触非线性有限元分析 作者：许小强赵洪伦 摘要基于非线性有限元软件MARC，提出过盈配合应力的动态和静态两种有限元分析方法，并以铁道车辆某高速轮对组装的过盈装配为例进行了有限元仿真计算，比较了两种方法的计算结果，分析了过盈量、摩擦系数、形状误差对装配应力的影响，结果对于确定合理过盈量和改进加工工艺具有参考意义。 关键词过盈配合接触非线性接触应力 0引言 在机械工程实际中普遍采用过盈配合来传递扭矩和轴向力，例如轴承配合、轴瓦配合、铁道车辆的轮轴、制动盘等。它是利用过盈量产生半径方向的接触面压力，并依靠由该面压力产生的摩擦力来传递扭矩和轴向力。由于过盈配合两个相配合的接触面上不能粘贴应变片，因此难以对其应力状态进行测定，对整个组装过程的应力状态更难以进行跟踪研究，而且这种配合方式往往承受着交变载荷的作用，配合面间可能发生相对滑动，这一滑动是随着应力变化而变化的，因而配合面边缘的接触状态和应力状态也随着应力的交变而变化，表现出复杂的状态，因此一般只能凭经验确定采用的过盈量。从力学角度看，这类问题属于接触非线性问题，传统的弹性接触解法已难以处理，可采用光弹性模拟实验进行研究，但只能反映应力分布趋势。近年来，随着非线性理论的不断完善和计算机技术的飞速发展，利用非线性有限元法来分析这类问题已日趋成熟。 铁道车辆随着向高速、重载不断发展，对轮轴的安全性要求也越来越高。研究表明，轮轴配合部位的应力状态对车轴的疲劳强度具有重要的影响，因此对轮对配合部位的宏观接触应力状态进行研究将有助于指导轮对制造标准的制定、高速重载轮对的设计和加工工艺的改进，以提高轮对的抗疲劳性能。 本文利用著名非线性有限元软件MARC，针对过盈配合的压力压装法和温差组装法对这类问题提出动态和静态两种仿真计算方法，并以铁道车辆某高速轮对的配合为例进行了计算，对比了两种计算方法的结果，分析了过盈量、摩擦系数、形状误差等因素对装配应力的影响。
非线性有限元的有关著作和简要历史
非线性有限元的有关著作和简要历史 已经发表的一些成功的实验和专题文章，完全或者部分地对非线性有限元分析做出了贡献。仅论述非线性有限元的作者包括Oden（1972），Crisfield（1991），Kleiber（1998）和Zhong（1993）。特别值得注意的时Oden的书，因为它是固体和结构非线性有限元的先驱作者。最近的作者有Simo和Hughes（1998）、Bonet 和Wood（1997）。某些作者还部分的为非线性分析做出了贡献，他们是Belytschko 和Hunhes（1983），Zienkiewicz和Taylor（1991），Bathe（1996），以及Cook，Malkushe和Plesha（1989）。对于非线性有限元分析，他们的书提供了有益的入门指南。作为姐妹篇，线性有限元分析的论述也是有用的，内容最全面的是Hughes（1987）、Zienkiewicz和Taylor（1991）的著作。 非线性有限元方法有多种溯源。通过波音研究的工作和Turner，Clough，Martin和Yopp（1956）的著名文章，使线性有限元分析得以闻名，不久之后，在许多大学和研究所里，工程师们开始将方法扩展至非线性、小位移的静态问题。但是，它难以燃起早期有限元社会的激情和改变传统研究者们对于这些方法的鄙视。例如，因为考虑到没有科学的是实质，《Journal of Applied Mechanics》许多年都拒绝刊登关于有限元方法的文章。然而。对于许多必须涉及工程问题的工程师们，他们非常清楚有限元方法的前途，因为它提供了一种处理复杂形状真实问题的可能性。 在20世纪60年代，由于Ed Wilson发布了他的第一个程序，这种激情终于被点燃了。这些程序的第一代没有名字。在遍布世界的许多实验室里，通过改进和扩展这些早期在Berkeley开发的软件，工程师们扩展了新的用途，这些带来了对工程分析的巨大冲击和有限元软件的随之发展。在Berkeley开发的第二代线性程序称之为SAP（structural analysis program）。由Berkeley的工作发展起来的第一个非线性程序使NONSAP，它具有隐式积分进行平衡求解和瞬时问题求解的功能。 第一批非线性有限元方法文章的主要贡献者由Argyris（1965），Marcal 和King（1967）。不久，大批文章激增，而且软件随之诞生。当时在Brown大
求解温度场的非线性有限元方法
Ξ 求解温度场的非线性有限元方法 刘福来1, 杜瑞燕2 (1.东北大学信息科学与工程学院,辽宁沈阳 110004;2.河北青年干部管理学院教务处,河北石家庄 050031) 摘要:从G alerkin 有限元方法出发,对自由表面上的辐射换热的数学表达式不作线性化处理,而是把温 度场的求解问题转化为非线性代数方程组的求解问题,并且用Newton 迭代法计算了温度场. 关键词:温度场;有限元方法;Newton 迭代法 中图分类号:O 242.21 文献标识码:A 文章编号:100025854(2005)0120021204 由文献[1]知,求解二维待轧过程的温度场,就是要求下面微分方程和初值问题的解: 52T 5 x 2+52T 5y 2=1α5T 5t ;(1) -k 5T 5n =0,(x ,y )∈S 2; (2) -k 5T 5n =σεA (T 4-T 4 ∞),(x ,y )∈S 3; (3) T (x ,y ,0)=T 0(x ,y ). (4)其中:α=λ ρc 称为导温系数,λ,ρ和c 分别为热导系数、密度和比热;S 2为给出热流强度Q 的边界面; T ∞为环境温度;S 3为给出热损失的边界面.对轧制问题的温度场,常常考虑的几种边界面[1] 是:对称 面、自由表面和轧件与轧辊的接触面.在辐射面上,边界条件的数学表达式为σεA (T 4-T 4 ∞)(其中:σ为 Stefan 2Boltzmann 常数,ε为物体表面黑度,A 为辐射面积,T ∞为环境温度)是温度T 的4次幂,具有强 烈的非线性.以往在实际计算中有2种处理方法[2],一种是简化问题的物理模型,有时将表达式看成常 数,有时将边界条件转化成h r A (T -T ∞)(其中h r =σ ε(T 2+T 2∞)(T +T ∞)),在轧制问题中求解温度场时文献[1,3]都采用了这一方法;另一种是处理问题的数学方法,即用近似方法求解非线性的偏微分方程问题.例如,用数值分析的方法,文献[4]中利用了差分方法. 本文中,笔者从G alerkin 有限元法出发,对自由表面上辐射换热的数学表达式不作线性处理,而是直接对非线性代数方程组用Newton 迭代法计算温度场,以二维待轧过程温度场的有限元解析进行讨论.1 G alerkin 有限元方法简介 将待求解区域Ω剖分为E 个单元,每个单元4个节点.设N i 是形函数(i =1,2,3,4),用4节点线性等参单元,则单元内的温度为 T e =N 1T 1+N 2T 2+N 3T 3+N 4T 4={N }T {T}e , (5) 其中:{N }=(N 1,N 2,N 3,N 4)T ;{T}e =(T 1,T 2,T 3,T 4)T .设ω1,ω2,…,ωn 是一组基函数,用 G alerkin 方法求方程(1)～(4)的解,实际上是求c 1,c 2,…,c n ,使T n =c 1ω1+c 2ω2+…+c n ωn 满足 κ Ω ρc 5T n 5t -k 52T n 5x 2+ 52T n 5y 2 ωi d x d y =0,i =1,2,…,n. (6) 对式(6)应用Green 公式,有 Ξ收稿日期:2004 0105;修回日期:20040420 作者简介:刘福来(1975),男,河北省唐山市人,东北大学博士研究生. 第29卷第1期2005年 1月河北师范大学学报(自然科学版) Journal of Hebei Normal University (Natural Science Edition )Vol.29No.1Jan.2005
ABAQUS有限元接触分析的基本概念
ABAQUS有限元接触分析的基本概念2009-11-24 00:06:28 作者：jiangnanxue 来源：智造网—助力中国制造业创新—https://www.wendangku.net/doc/9f12098458.html, CAE(计算机辅助工程)是一门复杂的工程科学，涉及仿真技术、软件、产品设计和力学等众多领域。世界上几大CAE公司各自以其独到的技术占领着相应的市场。ABAQUS有限元分析软件拥有世界上最大的非线性力学用户群，是国际上公认的最先进的大型通用非线性有限元分析软件之一。它广泛应用于机械制造、石油化工、航空航天、汽车交通、土木工程、国防军工、水利水电、生物医学、电子工程、能源、地矿、造船以及日用家电等工业和科学研究领域。ABAQUS在技术、品质和可靠性等方面具有卓越的声誉，可以对工程中各种复杂的线性和非线性问题进行分析计算。 《ABAQUS有限元分析常见问题解答》以问答的形式，详细介绍了使用ABAQUS建模分析过程中的各种常见问题，并以实例的形式教给读者如何分析问题、查找错误原因和尝试解决办法，帮助读者提高解决问题的能力。 《ABAQUS有限元分析常见问题解答》一书由机械工业出版社出版。 16.1.1 点对面离散与面对面离散 【常见问题16-1】 在ABAQUS/Standard分析中定义接触时，可以选择点对面离散方法(node-to-surface-dis - cre-tization)和面对面离散方法(surface-to-surface discretization)，二者有何差别? 『解答』 在点对面离散方法中，从面(slave surface)上的每个节点与该节点在主面(master surface)上的投影点建立接触关系，每个接触条件都包含一个从面节点和它的投影点附近的一组主面节点。 使用点对面离散方法时，从面节点不会穿透(penetrate)主面，但是主面节点可以穿透从面。 面对面离散方法会为整个从面(而不是单个节点)建立接触条件，在接触分析过程中同时考虑主面和从面的形状变化。可能在某些节点上出现穿透现象，但是穿透的程度不会很严重。 在如图16-l和图16-2所示的实例中，比较了两种情况。
通用显式非线性有限元程序：LS-DYNA
通用显式非线性有限元程序：LS-DYNA LS-DYNA 是世界上最著名的通用显式非线性有限元分析程序，能够模拟真实世界的各种复杂问题，特别适合求解各种二维、三维非线性结构的碰撞、金属成型等非线性动力冲击问题，同时可以求解传热、流体及流固耦合问题。在工程应用领域被广泛认可为最佳的分析软件包。与实验的无数次对比证实了其计算的可靠性。 LS-DYNA 是功能齐全的几何非线性（大位移、大转动和大应变）、材料非线性（140多种材料动态模型）和接触非线性（50多种）软件。它以Lagrange 算法为主，兼有ALE 和Euler 算法；以显式求解为主，兼有隐式求解功能；以结构分析为主，兼有热分析、流体-结构耦合功能；以非线性动力分析为主，兼有静力分析功能（如动力分析前的预应力计算和薄板冲压成型后的回弹计算）;是通用的结构分析非线性有限元程序。 特色功能 ? 显式求解为主，兼有隐式算法，适合于求解高度非线性问题； ? 具有多种求解算法，以Lagrange 算法为主，兼有ALE、Euler 算法、SPH （Smoothed Particle Hydrodynamics）光顺质点流体动力算法和边界元法BEM（Boundary Element Method）； ? 具有160多种材料模型，是材料模型非常丰富的有限元软件； ? 具有50多种接触类型，是接触类型非常齐全的有限元软件； ? 极好的并行计算能力，包括分布式并行算法（MPP）和共享内存式并行（SMP）； ? 良好的自适应网格剖分技术，包括自适应网格细分和粗化； ? 行业化的专用功能：如针对汽车行业的安全带单元、滑环、预紧器、牵引器、传感器、加速计、气囊等。 客户价值 ? 拥有显式和隐式算法，各向异性材料模型，使得板成型、回弹、预应力计算等，可以连续求解； ? 多种控制选项和用户子程序使得用户在定义和分析问题时有很大的灵活性； ? MPP 版本大幅度减少计算时间，计算效率随计算机数目增多而显著提高； ? 与大多数的CAD/CAE 软件集成并有接口。 广州有道科技培训中心 h t t p ://w w w .020f e a .c o m
过盈配合应力的接触非线性有限元分析
过盈配合应力的接触非线性有限元分析 摘要基于非线性有限元软件MARC，提出过盈配合应力的动态和静态两种有限元分析方法，并以铁道车辆某高速轮对组装的过盈装配为例进行了有限元仿真计算，比较了两种方法的计算结果，分析了过盈量、摩擦系数、形状误差对装配应力的影响，结果对于确定合理过盈量和改进加工工艺具有参考意义。 关键词过盈配合接触非线性接触应力 0 引言 在机械工程实际中普遍采用过盈配合来传递扭矩和轴向力，例如轴承配合、轴瓦配合、铁道车辆的轮轴、制动盘等。它是利用过盈量产生半径方向的接触面压力，并依靠由该面压力产生的摩擦力来传递扭矩和轴向力。由于过盈配合两个相配合的接触面上不能粘贴应变片，因此难以对其应力状态进行测定，对整个组装过程的应力状态更难以进行跟踪研究，而且这种配合方式往往承受着交变载荷的作用，配合面间可能发生相对滑动，这一滑动是随着应力变化而变化的，因而配合面边缘的接触状态和应力状态也随着应力的交变而变化，表现出复杂的状态，因此一般只能凭经验确定采用的过盈量。从力学角度看，这类问题属于接触非线性问题，传统的弹性接触解法已难以处理，可采用光弹性模拟实验进行研究，但只能反映应力分布趋势。近年来，随着非线性理论的不断完善和计算机技术的飞速发展，利用非线性有限元法来分析这类问题已日趋成熟。 铁道车辆随着向高速、重载不断发展，对轮轴的安全性要求也越来越高。研究表明，轮轴配合部位的应力状态对车轴的疲劳强度具有重要的影响，因此对轮对配合部位的宏观接触应力状态进行研究将有助于指导轮对制造标准的制定、高速重载轮对的设计和加工工艺的改进，以提高轮对的抗疲劳性能。 本文利用著名非线性有限元软件MARC，针对过盈配合的压力压装法和温差组装法对这类问题提出动态和静态两种仿真计算方法，并以铁道车辆某高速轮对的配合为例进行了计算，对比了两种计算方法的结果，分析了过盈量、摩擦系数、形状误差等因素对装配应力的影响。 1 过盈装配接触非线性问题的求解方法 1．1 接触非线性问题的求解方法 过盈问题是接触问题的一种，属于边界条件高度非线性的复杂问题，其特点是在接触问题中某些边界条件不是在计算开始就可以给出，而是计算的结果，两接触体间的接触面积和压力分布随外载荷的变化而变化，同时还包括正确模拟接触面间的摩擦行为和可能存在的接触传热。用有限元法解接触问题以往常采用的物理模型是节点对模型，即将两接触物体的接触面划分成相同的网格，组成一一对应的节点对，并假定两接触体的接触力通过节点对传递，这种模型需预先知道接触发生的确切部位，以便施加边界单元，对于结构复杂问题和考虑摩擦的动态接触问题，点对模型将给结构离散和方程求解带来极大困难，从而难以解决。近年来提出的点面接触模型是把两接触体分为主动体和被动体，在分析时研究主动体的节点与被动体接触表面上相接触的自由度关系及变形的一致关系，从而确定接触边界条件，然后从边界变形协调的变分原理出发，建立整个接触系统的控制方程。这种模型能有效处理复杂接触表面和动态接触问题。
非线性有限元(河海教授-任青文)
第一章 绪论 1.1 引言 有限单元法已成为一种强有力的数值解法来解决工程中遇到的大量问题，其应用范围从固体到流体，从静力到动力，从力学问题到非力学问题。事实上，有限单元法已经成为在已知边界条件和初始条件下求解偏微分方程组的一般数值方法。 对于有限元的线性分析，我们已经比较熟悉，它已作为设计工具被广泛采用。但绝大多数实际问题属于非线性问题，根据产生非线性的原因，非线性问题主要有三种类型： 1． 材料非线性(简称材料非线性或物理非线性) 其特点是应力σ与应变ε之间为非线性关系，通常与加载历史有关，加载和卸载不是同一途径（图1.1），因而其物理方程 D εσ= 中的弹性矩阵D 是应变ε的函数。但材料非线性问题仍属于小变形问题，位移和应变是微量，其几何方程是线性的。土、岩石、混凝土等具有典型的材料非线性性质，所以，土坝、岩土地基的稳定性和加固，地下洞室和边坡的稳定性都应当按材料非线性问题处理。 2．几何非线性 几何非线性属于大变形问题，位移和应变或者它们中一个是有限量。可能会有三种情况:大位移(包括线位移和角位移)、小应变，小位移、大应变和大位移、大应变。此时反映应变和位移关系的几何方程是非线性方程，例如，正应变 z ε可表示为 ""+??+??+??]()()222+??= [(21x w x v x u x u x ε 剪应变xy γ表示为 ""+????+????+????+y w x w y v x v y u x u ??+??= y u x v xy γ 如果应力和应变之间的关系也是非线性的，就变成了更复杂的双重非线性问题。不过，在几何非线性问题中一般都认为应力在弹性范围内，σ与ε之间呈线性关系。工程中的实体结构和板壳结构都存在几何非线性问题，例如弹性薄壳的大挠度分析，压杆或板壳在弹性屈曲后的稳定性问题。 在采用有限元方法分析非线性问题时，以上两类都表现为结构的整体劲度矩阵K 不再是常量矩阵，而是结点位移δ的函数，还有一类问题是结点荷载R 与δ有关，这就是边界非线性问题，又称接触非线性。 3．接触非线性 由于接触体的变形和接触边界的摩擦作用，使得部分边界条件随加载过程而变化，且不可恢复。这种由边界条件的可变性和不可逆性产生的非线性问题，称为接触非线性。工程中有许多接触非线性的问题，如混凝土坝纵缝和横缝缝面的接触，面板堆石坝中钢筋混凝土面板与垫层之间的接触，岩体节理面或裂隙面的工作状态等。 目前，研究工程非线性问题比较有效的工具是非线性有限元方法。要使这一方法实用化，有两个问题必须解决。第一，由于非线性问题的数值计算工作量大大增加，需要有相当大计
非线性有限元分析(学习总结报告)
非线性有限元 博士研究生专业课课程报告
目录 第一章绪言 (1) 1.1 非固体力学非线性问题的分类[1] (1) 1.2 非线性问题的分析过程[1] (2) 1.3 非线性有限元分析的基本原理 (2) 1.4 钢筋混凝土非线性分析的特点、现状及趋势 (3) 第二章非线性方程组的数值解法 (4) 2.1逐步增量法[3,4,5] (4) 2.2迭代法[3,4,5] (6) 2.3收敛标准 (8) 2.3.1.位移收敛准则 (8) 2.3.2.不平衡力收敛准则 (8) 2.3.3.能量收敛准则 (9) 2.4结构负刚度的处理[4,5] (9) 第三章材料的本构关系 (13) 3.1 钢筋的本构关系 (13) 3.1.1 单向加载下的应力应变关系 (13) 3.1.2 反复加载下的应力应变关系 (14) 3.2 混凝土的本构关系 (14) 3.2.1 单向加载下的应力应变关系 (14) 3.2.2 重复加载下的应力应变关系 (14) 3.2.3 反复加载下的应力应变关系 (14) 3.3 恢复力模型的分类 (14) 3.4 恢复力的获得方法 (15) 第四章非线性有限元在结构倒塌反应中的应用 (17) 4.1 钢筋混凝土结构倒塌反应研究现状 (17) 4.2 钢筋混凝土的有限元模型 (17) 4.2.1分离式模型 (18) 4.2.2组合式模型 (19) 4.2.3整体式模型 (20) 4.3 倒塌反应中RC结构有限元分析方法的选择 (20) 4.3.1隐式有限单元法 (21) 4.3.2显式有限单元法 (22) 4.4 钢筋混凝土框架结构的倒塌反应分析 (22) 4.4.1基于隐式有限单元法的倒塌分析 (22) 4.4.2 基于显式有限单元法的倒塌分析 (23) 4.5显式有限法在倒塌反应分析中的问题 (24)
第10章(非线性有限元)分解
公式号、图号等 第十章 非线性动力有限元法 当机械结构受到较大的外载荷，或受到持续时间较短的冲击载荷作用时，结构会产生过大的变形, 以至于必须考虑结构几何大变形对结构整体刚度及固有频率的影响，即所谓的几何非线性影响。另外, 对于多数非线性动力学问题，还需要考虑材料非线性、接触非线性等方面的影响。 非线性动力学分析求解的基本方程有如下形式 0=-+P I u M (4.141) 式中，Ku u C I += 为粘性效应项，考虑阻尼、粘塑、粘弹等效应。P 为外部激励。 对于考虑各种非线性效应的动力学问题求解，需要对动力学方程进行直接时间积分。即非线性动力有限元分析具有如下特点：（1）问题分析过程需要考虑时间积分效应，不必做模态分析，不必提取固有频率；（2）采用直接积分方法求解非线性动力学方程，需要对时间作积分计算，因此计算量远远大于线性模态动力学方法；（3）非线性动力学分析中可以施加不同类型的载荷，包括结点力、非零位移、单元载荷；（4）在每个时间步上，进行质量、阻尼、及刚度的集成，采用完整矩阵，不涉及质量矩阵的近似；（5）可以同时考虑几何、材料和接触等多种非线性效应。 非线性动力有限元分析程序常采用隐式Hilber-Hughes-Taylor 法进行时间积分运算。这种方法适于模拟非线性结构的动态问题，对于冲击、地震等激发的结构动态响应以及一些由于塑性或粘性阻尼造成的能量耗散，隐式算法特别有效。隐式积分方法需要对刚度矩阵求逆计算，并通过多次迭代求解增量步平衡方程。隐式Hilber-Hughes-Taylor 时间积分算法为无条件稳定，对时间步长没有特别的限制。 采用子空间法也可以对动力学平衡方程作时间积分运算。子空间法是提取模态分析得到的各阶特征模态，并采用与线性模态动力学分析方法相近的分析方式进行求解。对于带有微小非线性效应的问题，如材料小范围进行入屈服、结点转角不大的情况，子空间法效率比进接积分法要高。 此外，非线性动力有限元分析还可以采用显式动态算法，如中心差分法。显式时间积分算法为有条件稳定，其临界稳定时间步长限制了时间步长的大小，与有限元模型最小单元尺寸、材料应力波速等有关。显式时间积分法适于模拟高速冲击、接触等问题。 上述方法的选择需要综合考虑计算量、分析问题的规模、单元限制等多方面因素，需要丰富的有限元模拟的理论、经验和实践知识。以下以几何非线性问题和材料非线性问题为例介绍非线性有限元法，其中粘弹粘塑性非线性材料问题的分析是典型的非线性动力有限元的求解思想。 9.1 几何非线性问题的有限元法 几何非线性问题一般是指物体经历大的刚体位移和转动，但固连于物体坐标系中的应变分量仍假设为小量, 即大位移小应变情况。