下载文档原格式
非线性有限元方法及实例分析梁军河海大学水利水电工程学院，南京（210098）摘 要：对在地下工程稳定性分析中常用的非线性方程组的求解方法进行研究，讨论了非线性计算的迭代收敛准则，并利用非线性有限元方法分析了一个钢棒单轴拉伸的实例。
关键词：非线性有限元，方程组求解，实例分析1引 言有限单元法已成为一种强有力的数值解法来解决工程中遇到的大量问题，其应用范围从固体到流体，从静力到动力，从力学问题到非力学问题。
有限元的线性分析已经设计工具被广泛采用。
但对于绝大多数水利工程中遇到的实际问题如地下洞室等，将其作为非线性问题加以考虑更符合实际情况。
根据产生非线性的原因，非线性问题主要有3种类型[1]：1．材料非线性问题（简称材料非线性或物理非线性） 2．几何非线性问题3．接触非线性问题（简称接触非线性或边界非线性）2 非线性方程组的求解在非线性力学中，无论是哪一类非线性问题，经过有限元离散后，它们都归结为求解一个非线性代数方程组[2]：()()()00021212211=……==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 ×ijk R 为已知矢量。
在位移有限元中，δ代表未知的结点位移，()δF 是等效结点力，R 为等效结点荷载，方程()0=δψ表示结点平衡方程。
在线弹性有限元中，线性方程组0＝－R K δ （1.6）可以毫无困难地求解，但对线性方程组()0=δψ则不行。
一般来说，难以求得其精确解，通常采用数值解法，把非线性问题转化为一系列线性问题。
第十章 巴拿赫(Banach)空间中的基本定理1. 设X 是赋范线性空间，12,,,k x x x 是X 中K 个线性无关向量，12,,,k ααα是一组数，证明：在X 上存在满足下列两条件：（1）(),1,2,,v v f x v k α==,(2) M f ≤ 的线性连续泛函f 的充要条件为：对任何数12,,,k t t t ，11kkv vv vv v t Mt xα==≤∑∑都成立。
证明 必要性。
若线性连续泛函f 满足（1）和（2），则1111()kkkkv vv v v vv vv v v v t f t x ft xMt xα=====≤≤∑∑∑∑充分性。
若对任意数12,,,k t t t ，有11kkv vv vv v t Mt xα==≤∑∑。
令0X 为12,,,k x x x 张成的线性子空间。
对任意01kv vv t xX =∈∑，定义上线性泛函：0011:()k kv v v v v v f f t x t α===∑∑。
因0111()k kkv v v v v v v v v f t x t Mt x α====≤∑∑∑，故0f是有界的，且0f M ≤。
由泛函延拓定理，存在X 上的线性连续泛函f ，使f 限制在0X 上就是0f 。
f 显然满足条件（1）和（2）。
证毕。
2．设X 是赋范线性空间，Z 是X 的线性子空间，0x X ∈，又0(,)0d x Z >，证明存在'f X ∈，满足条件： 1）当x Z ∈时，()0f x =； 2）00()(,)f x d x Z = ；3）1f = 。
证明 记0{,}M x y C y Z λλ=+∈∈。
在M 上定义泛函0f ：000()(,)f x y d x Z λλ+=，则以下三条件成立：1）当y Z ∈时，0()0f y =； 2）00()(,)f x d x Z =；3）0f 在M 上有界，且01Mf =。
其中3）可以这样证明：若0x y M λ+∈，则00000()(,)yf x y d x Z x x y λλλλλ+=≤+=+，所以01Mf ≤。
非线性有限元分析1 概述在科学技术领域内，对于许多力学问题和物理问题，人们已经得到了它们所应遵循的基本方程（常微分方程或偏微分方程）和相应的定解条件（边界条件）。
但能够用解析方法求出精确解的只是少数方程性质比较简单，并且几何形状相当规则的问题。
对于大多数工程实际问题，由于方程的某些特征的非线性性质，或由于求解区域的几何形状比较复杂，则不能得到解析的答案。
这类问题的解决通常有两种途径。
一是引入简化假设，将方程和几何边界简化为能够处理的情况，从而得到问题在简化状态下的解答。
但是这种方法只是在有限的情况下是可行的，因为过多的简化可能导致误差很大甚至是错误的解答。
因此人们多年来一直在致力于寻找和发展另一种求解途径和方法——数值解法。
特别是五十多年来，随着电子计算机的飞速发展和广泛应用，数值分析方法已成为求解科学技术问题的主要工具。
已经发展的数值分析方法可以分为两大类。
一类以有限差分法为代表，主要特点是直接求解基本方程和相应定解条件的近似解。
其具体解法是将求解区域划分为网格，然后在网格的结点上用差分方程来近似微分方程，当采用较多结点时，近似解的精度可以得到改善。
但是当用于求解几何形状复杂的问题时，有限差分法的精度将降低，甚至发生困难。
另一类数值分析方法是首先建立和原问题基本方程及相应定解条件相等效的积分提法，然后再建立近似解法并求解。
如果原问题的方程具有某些特定的性质，则它的等效积分提法可以归结为某个泛函的变分，相应的近似解法实际上就是求解泛函的驻值问题。
诸如里兹法，配点法，最小二乘法，伽辽金法，力矩法等都属于这一类方法。
但此类方法也只能局限于几何形状规则的问题，原因在于它们都是在整个求解区域上假设近似函数，因此，对于几何形状复杂的问题，不可能建立合乎要求的近似函数。
1960年，R.W.CLOUGH发表了有限单元法的第一篇文献“The Finite Element Method in Plane Stress Analysis”，这同时也标志着有限单元法（FEM）的问世。
（一）非线性方程的迭代解法1.非线性方程的一般形式：f(x)=02.非线性方程的分类：⎩⎨⎧=为其他函数。
超越方程，次代数多项式；为代数方程，)()(0)(x f n x f x f 3.方程的根：若存在常数s 使f(s)=0，则称s 是方程(4.1)的根，又称s 是函数f(x)的零点。
4.重根：若f(x)能分解为)()()(x s x x f m ϕ-= 则称s 是方程(4.1)的m 重根和f(x)的m 重零点。
当m=1时，s 称为方程(4.1)的单根和f(x)的单零点。
5.结论：(1)零点存在定理：设函数f(x)在闭区间[a,b]上连续，且f(a)•f(b)<0，那么在开区间（a,b ）内至少有一点ξ，使f(ξ)=0.(2)根的唯一性判别：一阶导数不变号且不为零(3)n 次代数方程在复数域上恰有n 个根(4)高于4次的代数方程没有求根公式6.方法：(1)搜索根方法：①作图法：②逐步搜索法：确定方程根的范围的步骤：步骤1 取含f(x)=0根的区间[a,b]，即f(a)•f(b)<0；步骤2 从a 开始，按某个预定的步长h ，不断地向右跨一步进行一次搜索， 即检查kh a x k +=上的函数)(k x f 值的符号。
若0)()(1<•-k k x f x f ，则可以确定一个有根区间],[1k k x x -.步骤3 继续向右搜索，直到找出[a,b]上的全部有根区间],[1k k x x -（k=1,2,…,n).(2)二分法①基本思想：含根区间逐次分半缩小，得到一个区间长度以1/2的比例减小的含根区间序列 {}k I ,在给定根的误差界时，利用长度趋于零的特点，可得到在某个区间中满足要求的近似根。
②迭代终止的条件ε<)(k x fε2<-k k a b或者ε<-≤-2k k k a b s x(3)简单迭代法及其收敛性)(0)(x x x f ϕ=⇔=,2,1,0),(1==+k x x k k ϕ迭代法是一种逐次逼近法,用某个固定公式反复校正根的近似值,使之逐 步精确化,最后得到满足精度要求的解。
第10章 非线性动力有限元法 (1)10.1 几何非线性问题的有限元法 (2)10.1.1 几何非线性问题的牛顿迭代法 ........................................................................... 2 10.1.2 典型单元的切线刚度矩阵 ................................................................................. 4 10.2 材料非线性问题的有限元法 (8)10.2.1 弹/粘塑性问题的基本表达式 .............................................................................. 8 10.2.2 粘塑性应变增量和应力增量 ............................................................................... 9 10.2.3 弹/粘塑性平衡方程 ............................................................................................ 10 10.3 材料非线性问题的动力有限元法 ................................................................................ 11 10.4 应用举例 (14)10.4.1 粘弹粘塑性动力有限元分析举例 ................................................................... 14 习题.. (15)第10章 非线性动力有限元法当机械结构受到较大的外载荷，或受到持续时间较短的冲击载荷作用时，结构会产生过大的变形, 以至于必须考虑结构几何大变形对结构整体刚度及固有频率的影响，即所谓的几何非线性影响。
另外, 对于多数非线性动力学问题，还需要考虑材料非线性、接触非线性等方面的影响。
非线性动力学分析求解的基本方程有如下形式0=-+P I uM (4.141) 式中，Ku uC I += 为粘性效应项，考虑阻尼、粘塑、粘弹等效应。
P 为外部激励。
对于考虑各种非线性效应的动力学问题求解，需要对动力学方程进行直接时间积分。
即非线性动力有限元分析具有如下特点：（1）问题分析过程需要考虑时间积分效应，不必做模态分析，不必提取固有频率；（2）采用直接积分方法求解非线性动力学方程，需要对时间作积分计算，因此计算量远远大于线性模态动力学方法；（3）非线性动力学分析中可以施加不同类型的载荷，包括结点力、非零位移、单元载荷；（4）在每个时间步上，进行质量、阻尼、及刚度的集成，采用完整矩阵，不涉及质量矩阵的近似；（5）可以同时考虑几何、材料和接触等多种非线性效应。
非线性动力有限元分析程序常采用隐式Hilber-Hughes-Taylor 法进行时间积分运算。
这种方法适于模拟非线性结构的动态问题，对于冲击、地震等激发的结构动态响应以及一些由于塑性或粘性阻尼造成的能量耗散，隐式算法特别有效。
隐式积分方法需要对刚度矩阵求逆计算，并通过多次迭代求解增量步平衡方程。
隐式Hilber-Hughes-Taylor 时间积分算法为无条件稳定，对时间步长没有特别的限制。
采用子空间法也可以对动力学平衡方程作时间积分运算。
子空间法是提取模态分析得到的各阶特征模态，并采用与线性模态动力学分析方法相近的分析方式进行求解。
对于带有微小非线性效应的问题，如材料小范围进行入屈服、结点转角不大的情况，子空间法效率比进接积分法要高。
此外，非线性动力有限元分析还可以采用显式动态算法，如中心差分法。
显式时间积分算法为有条件稳定，其临界稳定时间步长限制了时间步长的大小，与有限元模型最小单元尺寸、材料应力波速等有关。
显式时间积分法适于模拟高速冲击、接触等问题。
上述方法的选择需要综合考虑计算量、分析问题的规模、单元限制等多方面因素，需要丰富的有限元模拟的理论、经验和实践知识。
以下以几何非线性问题和材料非线性问题为例介绍非线性有限元法，其中粘弹粘塑性非线性材料问题的分析是典型的非线性动力有限元的求解思想。
10.1 几何非线性问题的有限元法几何非线性问题一般是指物体经历大的刚体位移和转动，但固连于物体坐标系中的应变分量仍假设为小量, 即大位移小应变情况。
10.1.1 几何非线性问题的牛顿迭代法由数值分析技术可知，求解非线性方程组的数值方法的常规方法是Newton-Raphson 法，即牛顿迭代法，这是一种近似线性化迭代求解方法。
对于非线性方程0)(=x ψ，具有一阶导数，在n x 点作一阶泰勒级数展开，它在n x 点的线性近似为d ()()()()d n n n x x x x xψψψ=+- (4.142) 因此，非线性方程0)(=x ψ在n x 附近似为线性方程：d ()()()0d n n n x x x xψψ+-= (4.143) 当d () 0d n xψ≠时，由上式求得n 步的修正项 1d ()/()d n n n X x xψ∆ψ+=- (4.144) Newton-Raphson 方法的迭代公式为11++∆+=n n n X x X (4.145)在几何非线性有限元法中，结构的刚度矩阵与其几何位置有关，平衡方程由变形后的位形描述，因此，结构的刚度矩阵是几何变形的函数。
设变形为δ, 结构的平衡方程式()0-=K δδR (4.146)为一个非线性方程组。
记非线性方程()0K =-=ψδδR (4.147)用Newton-Raphson 方法求()0=ψδ的根时，迭代公式分别为11n n n ++=+δδδ (4.148)其中,1n ∆+δ满足下式1()T n n n ∆+=-K δR K δδ (4.149)式中, T n K 称为切线刚度矩阵，表达式为d ()()d T n n =ψδK δ(4.150) 在每一个迭代步中，通过求解切线刚度矩阵T n K ，进而用1n ∆+δ进行迭代求解，称为Newton-Raphson 方法，又称切线刚度法。
牛顿法的收敛性是好的。
但是某些非线性问题中，使用牛顿法迭代时，若T n K 出现奇异或病态，则对T K 的求逆出现困难。
关于这一点也可以采用其它修正办法，如引入阻尼因子。
对于已经建立的有限元方程，设ψ表示内为和外力矢量的总和，有***d 0T T T VV =-=⎰δψεσδR (4.151)式中, R 为载荷列阵；*δ为虚位移；*ε为虚应变用应变的增量形式d d =εB δ代入上式，消去*δ项，可以得到非线性问题的一般平衡方程式为()d 0T VV =-=⎰ψδB σR (4.152)该式不论位移或应变的大小与否均成立。
在有限变形中，应变和位移之间的关系是非线性的，即B 矩阵是δ的非线性函数。
但是，近似地可将进行如下分解：0L =+B B B (4.153)式中, 0B 为线性应变分析的部分； L B 为由非线性变形引起的，与δ有关。
假定应力应变关系为线弹性，于是有00()=-+σD εεσ (4.154)式中 ][D 为材料的弹性矩阵； }{0ε为初应变列阵；}{0σ为初应力列阵对于式(4.152)的非线性平衡方程式，可用Newton-Raphson 方法进行迭代求解。
对该式微分，有d d d d d T T VVV V =+⎰⎰ψB σB σ (4.155)不考虑初应变和初应力的影响，得d d d ==σD εDB δ并且d d L =B B这样可得d d d d T L VV =+⎰ψB σK δ (4.156)这里0d T L VV ==+⎰K B DB K K (4.157)式中0K 为通常的小位移的线性刚度矩阵。
L K 矩阵则是由于大位移引起，它可以写成()00d T T T L L L L L VV =++⎰K B DB B DB B DB (4.158)式(4.156)又可记成：()0d d d L T σ=++=ψK K K δK δ (4.159)式中d d T L VV σ=⎰K B σ (4.160)式中，σK 是关于应力水平的对称矩阵，称之为初应力矩阵或几何刚度矩阵。
因此，用Newton-Raphson 方法迭代求解几何非线性问题的步骤为： (1) 用线弹性解作为1δ，即一次近似； (2) 通过定义1()δB 求出1σ，求出1ψ； (3) 确定切线刚度矩阵1T K ；(4)211/T ∆=-δψK ， 212∆=+δδδ；(5) 重复上述迭代步骤，直至n ψ足够小。
在这里，没有考虑载荷R 可能由于变形而发生的变化，即在这里假设了载荷不因变形而改变其大小和方向，否则是非保守力作用下的大变形问题，在此不做讨论。
10.1.2 典型单元的切线刚度矩阵求解具体的几何非线性问题时，必须计算单元的切线刚度矩阵。
对于一般空间问题，无论位移和应变大小，都可以利用应变的基本定义写出位移和应变的关系式。
用变形前的坐标),,(z y x 做为自变量，可以用位移w v u ,,定义如下大变形问题的应变分量表达式⎥⎦⎤⎢⎣⎡∂∂+∂∂+∂∂+∂∂=222)()()(21x w x v x u x u x ε2221()()()2y v u v w y yy y ε⎡⎤∂∂∂∂=+++⎢⎥∂∂∂∂⎣⎦⎥⎦⎤⎢⎣⎡∂∂+∂∂+∂∂+∂∂=222)()()(21z w z v z u z w z εzwy w z v y v z u y u y w z v xy ∂∂∂∂+∂∂∂∂+∂∂∂∂+∂∂+∂∂=γ (4.161) zw x w z v x v z u x u z u x w yz ∂∂∂∂+∂∂∂∂+∂∂∂∂+∂∂+∂∂=γ ywx w y v x v y u x u x v y u zx ∂∂∂∂+∂∂∂∂+∂∂∂∂+∂∂+∂∂=γ 对于微小位移情况，可以略去二次以上的偏导数项，得到小变形时的应变公式。
在有限变形中，假设应变仍为小量。
应变和位移之间的关系为：0L =+εεε (4.162)式中0ε为线性应变部分。
对于非线性部分，可以写成：000001122000T x T y x T z L y T T z y T T z zx T Ty x ⎡⎤⎢⎥⎢⎥⎧⎫⎢⎥⎪⎪⎢⎥==⎨⎬⎢⎥⎪⎪⎢⎥⎩⎭⎢⎥⎢⎥⎣⎦θθθθεθCθθθθθθθθ (4.163) 式中[][][]Tx Ty Tz u v w x x x u v w y y y u v w z zz∂∂∂=∂∂∂∂∂∂=∂∂∂∂∂∂=∂∂∂θθθ (4.164) 式中C 为96⨯矩阵。
