§1.2 线性有限元的回顾
线性有限元的理论虽然不是本课程的重点，而线性有限元法是非线性有限元法的基础；非线性有限元法又是线性有限元法的发展。因为非线性问题的求解通常是采用分段线性化的思想，使其成为一系列的线性问题。因此有必要首先回顾线性有限元的一些基本内容。主要是线性有限元的基本理论和方法，以及当前应用最广的等参元理论。
固体力学的理论（材力、弹力、结构力学，无论是线性还是非线性）是建立在本构方程、几何（运动）方程以及平衡方程三方面方程的基础上的，由此导出的控制方程的解是满足上述充分、必要条件的唯一解，而且是反映了结构真实受载后的运动状态（变形）。在结构分析中许多情况下，本构和几何方程可以处理成线性，使控制方程线性化，求解也大大简化，同时可达到工程上的精度要求，这就是线性有限元的基本理论。 §1.2.1 线性本构关系（广义虎克定律）
影响结构材料性能有诸多因素（应力、应变、变形率、温度、湿度、时效等），而通常建立本构方程时仅考虑应力、应变两个物理参数，认为两者成线性关系的经典的弹性理论，即著名的虎克定律。
各向同性材料的Hooke 定律
ij ijkl kl D σε= 或 ij ijkl kl C εσ=
()1其中ijkl
D 和ijkl C 分别为弹性阵和柔度阵。
由剪应力互等定理，弹性阵()99ijkl D ?独立材料参数的个数由81个减少为21个。进而对于正交异性的材料参数为9个独立的参数，对于各向同性材料：
01121,
2,322(1)
,e
i j i j
i
j
E G G
E
d dS d i j νν
εσδ
-+=
+=
= ()2
仅有两个独立的材料常数，即E 和ν，其中E 为弹性模量，ν为泊松比。 §1.2.2线性几何方程（小变形情况）
线性（小变形）关系： ()(){}1
,,2
ij T
U U U i j j i ε+=?=
()3 位移边界条件：u S 边界上 U U = 其中U 为位移向量，
U 为边界u S 上的指定位移，()T
?为微分算子。
实际结构可根据它们的几何特点，将三维问题简化为二维问题，这里主要有：平面应力、平面应变和轴对称状态。
1）平面应力（薄壁结构）
外力仅作用在平面内，两表面无外力作用（离心力作用下圆盘）两表面应力分量（且在整个厚度方向上），yz σ=xz σ=zz σ=0，而xx σ、yy σ、xy σ沿厚度均匀分布。
按Hooke 定律，
()()
2410
xx yy xy
xx yy zz zz xz yz u x v y
v u x y ννεεεεεεεεε??????
??+ ???????+
?-??
==
=
=-===
得弹性阵
[]
()
102
/(1)105100
2
D E ν
ννν????????=-??????
??-???
?
2）平面应变（沿轴线的几何形状和外载无明显变化的物体），如水坝，挡土墙等很长的物体，在离两端一定距离处，可以假定任意横截面上的位移、应力、应变等力学量仅是x 、y 的函数，沿物体纵向无变化，进一步假设沿z 方向位移w 为常数或为0，则有
[]()()
(
)
()
,2101061121200
20,zz yy xz
z y x v u w y
z x E
D ννννννν
εεεσνσσ?????
+ ??????????
-??=
-??
+-??-???
?
==
=
=+
3）轴对称问题（广义的平面应变）
用轴坐标代替直角坐标，（例如厚壁筒，高压罐，炮筒）。
r u r
ε??=
，z w z
ε??=，u r
θε=
，2rz u w z
r ε????
+ ?????
= ()16
[]
()
()
1011101111211200
0210
111E D νννννννννννννννν
ν
?
?
??
--??????
---??=
-+-????-??????--??
()26
§1.2.3 平衡方程（虚功原理） 经典力学给出了结构的平衡方程：
,0ji j
i f σ
+= or ,0ij i j f σ+= ()17
边界条件 j i j i l f σ= or i e t t = 在S S σ=上
()27
ij σ为应力张量，i f 为体力，i t 为结构内力，e t 为外部表面力，直接由平衡方程出发导
出适合各类边界条件的有限元控制方程是困难的，但是从经典的虚功原理或加权余量原理出发，可得到广泛使用的有限元控制方程。
虚功原理认为：一个处于平衡的物体当发生满足位移边界条件的连续（微小）虚位移
i u δ时，其外力虚功等于物体中应力在虚应变ij δε上产生的变形能。其表达式为：
ij ij i i e i V
V S
dV f u dV t u dS σδεδδ=+
?
??
()8
其中δε为对应u δ（虚位移）的虚应变，而虚位移应满足位移—应变的几何关系
(),,12
ij i j
j i u
u ε=
+
和位移边界条件U U =。以位移作为基本未知量的方程是当今有限元的主
流称为位移法（位移模型），一个变形体总位能π为：π= U W -，式中U 是变形能，W 是外力功，按虚功原理 0U W δπδδ=-= ()9 按照变分运算规则，仅对位移取变分，而力和应力不变，因此
i j i j V
U d V δσδε=
?
， i i i i V
S
W f u dV t u d S δδδ=
+
?
?
()10
代入()9式即()8式。
按照上述能量原理，产生了多种求解结构位移的方法，如李兹法，伽辽金法等，有限元法与经典李兹法的一个重要区别是：有限元法不寻求在结构整个域上的连续位移函数，而是寻找各个子域（单元）上，满足域内及域边界上连续的位移函数。
设将物体分成若干个互不重迭的单元（包括三角单元，矩形单元等），对每个单元，设节点位移{}Q 未知量，在单元内部和边界上满足位移协调条件的插值函数N 。于是体内位移场 {}{}{}U N Q = ()11
代入几何关系()3可得 {}{}{}{}{}{}
T
N Q B Q ε=?= ()12 其中{}T
?为()3式中微分算子。将其代入方程()8则有
{}
[]
{}{}
[]
{}{}
[]
{}()13T
T
T
T
T
T
V V S q B dV q N f dV q N t dS
σδσδδ=+???进一步将
本构关系[]D 代入，上式成为：[][][][]K q F T =+ ()14 其中 []
[]
[][]T
V K B D B dV ∑=? 单元刚度矩阵
[]
[]{}T
V F N f dV ∑=? 作用于节点的当量体力 ()15
[][]{}T
s T N t dS σ
∑=
? 作用于节点的表面力
这里∑是指整个物体中所有单元求和，积分仅对每个单元域，这样在实际计算中对各单元，
可用单元自身坐标系（局部坐标），然后通过坐标转换成总体坐标后再进行叠加。
如果以结构的内力或对应应力为基本未知量，位移仅是因变量，那么构成的变形能为余能如图。
将虚功原理用于余能构成虚余能原理：在满足平衡条件及应力边界条件的所有可能的力与应力中，满足协调条件的状态使余能取最小值。其数学表达式为
c c c u w π=-，
其中c π为余能，c u 为余变形能，c w 为外载余能。
在有些情况下，如板壳结构，断裂分析以及接触分析，使用节点位移及应力（或者当量节点力）混合的方法或用节点位移和边界当量应力/应变杂交的方法。或许会更方便，且比单纯位移模式精度更高。这些模式的推导主要基于广义变分原理（或加权余能原理）胡海昌-鹫津久一郎变分原理是最早提出的，也是最基本的。 §1.2.4 有限元方法的三个基本步骤
1．离散化。把连续体分割为有限个只在结点上相连接的单元组合体。
2．单元分析。a.取位移插值函数；b.建立单元几何关系，本构关系；c.通过虚位移原理，
建立平衡方程，即得到[][][]{}{}e
T V e B D B dV q F =?，并令单元刚度矩阵为：
[][][][]e e T
V k B D B d V =? ()16
3．总体分析
集合所有的刚度矩阵和节点载荷向量，成为整个系统的平衡方程。由于离散化的原因，此时不再是弹性力学中连续体的无限个自由度的偏微方程组，而是一个有限自由度的线代数方程组；
[]{}{}e e
k q
F ∑∑= →kq F = ()17 通过求解，可得该系统各节点上的总位移向量{}q ，进而可求得各单元中的应变和应
力，从{}q 提取{}e
q ，得到单元中应力和应变；
{}[]{}{}[]{},e
e
e
e
B q D εσε== ()18 §1.3 等参单元
单元如同生物中的细胞，是有限元分析中最基本的部分，因此如何有效地建立单元刚度阵及计算它们是很重要的环节。在实际分析中，等参单元是非常重要的一类元素。艾恩斯（IRONS ）等人提出的等参单元理论，是对有限元的一大贡献。由于等参单元普遍具有精度高可靠性好的特点，在国内外当前流行的大型通用有限元程序中，品质好的等参元几乎成为单元库的主体，线性分析如此，非线性分析亦如此，因此有必要回顾一下它的基本理论。
§1.3.1 单元位移插值
在单元内选取适当的位移模式，将单元内部各点位移（场）用该单元节点上的位移表
示。 即：
{}
[]{}i i i
u u v N q w ????==∑??????
()1~i n = ()19 形函数应满足的条件：
a)单元内位移连续且保证单元间位移协调（conforming ）。 b)保证其刚体位移（常数项），保证单元作刚体运动时无应变。 c)包含任意线性项，满足常应变条件，保证其收敛性。
一旦位移模式选定之后，就可以按几何、物理、平衡三方面条件推导单元刚度矩阵，同时在一定的单元数目下，解的精度也就确定了。因此选择单元模式很关键。因此，要改善解的精度，也得从改进单元模式入手。 例如： 1. 一维两结点线性杆单元
设i j u N u N u =+ ()20
其中 11x N l
=-
2x N l
=
2. 常应变（常应力）三角形平面单元（二维）
定义： ?为ijk ?的面积。i ?，j ?，k ?分别为 pjk ?，pik ?，pji ?的面积。有
?
=?∑，
1i
L
=∑，且结点i 上， 1i L =，
0j k L L == 采用面积坐标：i i L ?=
?
，j j L ?=
?
，k k L ?=
?
。
?
内任一点的面积坐标为：（i L ，j L ，k L ），
设面内两个方向的位移为线性函数： 3
1
i i
i u L u ==
∑，
3
1
i i
i v L v ==
∑ ()21
以上直接选取的线性位移模式是对实际问题的最低阶逼近，精度受到很大限制，实际工程中需要的单元以形状、性能都比较复杂，下面特别要介绍的是等参单元。 §1.3.2 基于坐标变换的等参单元
前面介绍的杆和三角形平面单元，形状特别简单，可在物理坐标中直接找到合适的位移模式（插值函数），但又难于适应工程中的复杂形状和性能，为了解决这一矛盾，等参元利用了坐标变换。
1)坐标变换：
下面我们称实际具有曲线边界，形状复杂的单元，在物理空间,,x y z 坐标下定义为子单元；而假设通过坐标变换将子单元变换到某一参数空间的(),,ξηζ坐标系中形状简单的母单元。这样一方面子单元的形状，特征，荷载条件都来自于实际结构；另一方面，大量的计算工作是在母单元上进行的，由于它形状简单规律、计算方便。
举例：线性矩形单元
这也是一种常用的单元，它采用比三角形常应变单元有较高的位移模式，以便更好的反映弹性体中的位移和应力状态。
矩形单元1,2,3,4 面积 22a b ?，我们可以在xy 坐
标下（物理空间）推导单元。但如果采用坐标变换（参数坐标），引入一个局部坐标(),ξη可推导出非常简单的结果。则坐标变换为：
x x a ξ=+，0y y b η=+ ()22 其中：()0000,,,i j x y x y x y = ()1~4i =
为矩形单元形心坐标，与结点整体坐标有关。
这时的母单元各结点局部坐标为()1±。于是可在母单元上取位移模式：
4
1
i
i
i u N u
==
∑ 4
1
i
i
i v N v
==
∑ ()23
式中 ()()0011/4i N ξη=++ 其中 0i ξξξ=，0i ηηη=
，
i =
()24
上述的坐标变换函数和位移插值函数是不相同的。前者为线性函数（各边形状从直线到直线），而后者为二次函数。所以虽然引入了参数坐标，作了坐标变换，但不是等参变换。只有当两者的插值函数具有相同形式时，构造的单元是等参元。位移插值高于坐标插值的称为超参元，反之是亚参元。 2）等参变换及单元刚度阵的推导：
2-1) 以平面坐标系为例，将ξη参数平面上n 个结点的单元，映射到xy 物理平面上。
如：正方形→曲边单元 ，4
1
i
i i x N
x ==
∑， 4
1
i
i i y N
y ==
∑ ()25
(),i i N N ξη=是用局部坐标表示的形函数。,i i x y 是结点i 的整体坐标。 2-2)用同样的形函数，在母单元上（参数坐标下）()0ξη作位移插值
4
1
i i
i u N u
==
∑ 4
1
i
i
i v N v
==
∑ ()26
可以证明，将母单元上位移移植到子单元上时，能够满足完备性（刚体位移、常应变）和协调性。
2-3)在子单元上推导其单元刚度矩阵。
经过坐标变换的单元，由于物理坐标中单元已变换到参数坐标（局部坐标）系下，
所以单元刚度阵也将在局部坐标下计算。 由 e
T
k B DBdxdy =
?? 得
[]
[]()[][]()
11
11
,,e
T
k B D B J d d ξηξηξη--=
?? ()27
其中det J J = 3)高斯求积方法的应用：
精确快速计算结构刚度矩阵和等效结点力各元素，是有限元法中重要环节在等参元中
经常遇到的积分形式有：()111
111
,,f d d d ξηζξηζ
---?
?? 其中(),,f ξηζ往往很复杂，因此通
常采用数值积分，高斯积分法是其中精度较高的方法。原理是：在单元内选出某些点为积分点，计算出被积函数(),,f ξηζ在这些点上的函数值，然后用一些加权系数，乘上这些函数值，求出总和，作为近似积分值。
高斯积分公式：
()()111
1
11
111
,,,,n
m
l
i
j k i j k i j k f d d d H
H H f ξηζξηζ
ξηζ===---=
∑∑∑???式中 n —积
分点数目； i ξ；j η；k ζ—分别为i 点；j 点；k 点的局部坐标。i H ；j H ；k H —加权系数。
不同积分点数的局部坐标（参数坐标）和加权系数可查有关数据表格。在平面问题、轴对称问题及空间问题中，物体划分较多个等参单元时，使用22?和222??的高斯积分法则，通常能取得良好效果。其原因是位移插值函数中的不完全高次项，往往对计算精度起到不良影响，降阶积分的措施正好相当于把那些不完全的高次项对刚度阵的影响去除掉。
轨道结构的非线性有限元分析 姜建华 练松良 摘 要 实际轨道结构受载时的力学行为,属于典型的非线性力学问题。钢轨垫层刚度、钢轨抗扭刚度和扣件扣压力的大小是影响轨距扩大的主要因素。根据非线性有限元接触理论,建立了能准确反映扣件、钢轨与垫层的拧紧接触,以及受载车轮与钢轨侧向滑动接触的力学计算模型;并研究计算了不同扣件压力下,由于受载车轮与钢轨侧向滑动接触引起的轨距扩大问题。 关键词 轮轨关系,扣件压力,非线性弹性力学,有限元分析 1 引言 实际工程中常见的非线性问题一般可以归纳为三类:材料非线性、几何非线性以及边界条件非线性。材料非线性问题是由于材料的非线性本构关系所引起的,例如材料的弹塑性变形,材料的屈服和硬化等;几何非线性问题是由于结构的位移或变形相当大,以至必须按照变形后的几何位置来建立平衡方程;边界条件非线性问题是指边界条件随位移变化所引起的非线性问题。通常情况下,我们所遇到的非线性问题多数是上述三类非线性问题的组合[1,2]。 实际轨道结构受载时的力学行为,属于典型的非线性力学问题。比如基于轮轨接触的材料非线性、几何非线性及边界条件非线性问题,以及扣件、钢轨、垫层三者间相互作用时所表现的边界条件非线性行为等。所以,机车车辆在轨道结构上行驶时引起的力学现象是相当复杂的。以往在研究轨道各部分应力应变分布规律时,通常采用连续弹性基础梁理论或连续点支承,偶尔简单考虑扣件的作用和弹性垫层的使用。不管用哪一种支承方式建立模型,都由于这样那样的假设而带有一定程度的近似性。所以,如何利用现代力学理论的最新成果以及日益发展的计算机技术,根据轨道结构的具体情况,建立更为完整更为准确的轨道结构计算模型,为轨道设计部门提供更加可靠的设计依据或研究思路,已十分必要。 本文提出了用非线性有限元理论研究轮轨系统和轨道结构的思路。作为算例之一,本文将根据非线性有限元理论,建立能准确反映扣件、钢轨与垫层的拧紧接触,以及受载车轮与钢轨侧向滑动接触的力学计算模型。 2 轨道结构的有限元接触模型 对于非线性问题,不管是材料非线性、几何非线性,还是边界条件非线性,总是最终归结为求解一组非线性平衡方程及其控制方程。例如用位移作为未知数进行有限元分析时,最后可得到一组平衡方程及其控制方程为 : 图1 轮轨系统的对称性模型简图 [K(u)]{u}={R}(1) (u)= (u)(2)其中:{u}为节点位移列阵;{R}为节点载荷列阵; [K(u)]为总体刚度矩阵; (u)为边界条件。它们 36 姜建华:同济大学工程力学系,副教授、博士,上海200092
非线性有限元方法及实例分析 梁军 河海大学水利水电工程学院，南京（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=δψ表示结点平衡方程。 在线弹性有限元中，线性方程组
第9章非线性问题的有限单元法 9.1 非线性问题概述 前面章节讨论的都是线性问题，但在很多实际问题中，线弹性力学中的基本方程已不能满足，需要用非线性有限单元法。非线性问题的基本特征是变化的结构刚度，它可以分为三大类：材料非线性、几何非线性、状态非线性。 1. 材料非线性(塑性, 超弹性, 蠕变) 材料非线性指的是材料的物理定律是非线性的。它又可分为非线性弹性问题和非线性弹塑性问题两大类。例如在结构的形状有不连续变化（如缺口、裂纹等）的部位存在应力集中，当外载荷到达一定数值时该部位首先进入塑性，这时在该部位线弹性的应力应变关系不再适用，虽然结构的其他大部分区域仍保持弹性。 2. 几何非线性(大应变, 大挠度, 应力刚化) 几何非线性是有结构变形的大位移引起的。例如钓鱼杆，在轻微的垂向载荷作用下，会产生很大的变形。随着垂向载荷的增加，杆不断的弯曲，以至于动力臂明显减少，结构刚度增加。 3. 状态非线性(接触, 单元死活) 状态非线性是一种与状态相关的非线性行为。例如，只承受张力的电缆的松弛与张紧；轴承与轴承套的接触与脱开；冻土的冻结与融化。这些系统的刚度随着它们状态的变化而发生显著变化。 9.2 非线性有限元问题的求解方法 对于线性方程组，由于刚度方程是常数矩阵，可以直接求解，但对于非线性方程组，由于刚度方程是某个未知量的函数则不能直接求解。以下将简要介绍借助于重复求解线性方程组以得到非线性方程组解答的一些常用方法。 1.迭代法 迭代法与直接法不同，它不是求方程组的直接解，而是用某一近似值代人，逐步迭代，使近似值逐渐逼近，当达到允许的规定误差时，就取这些近似值为方程组的解。 与直接法相比，迭代法的计算程序较简单，但迭代法耗用的机时较直接法长。它不必存贮带宽以内的零元素，因此存贮量大大减少，且计算中舍入误差的积累也较小。以平面问题 为例，迭代法的存贮量一般只需直接法的14左右。在求解非线性方程组时，一般采用迭代 法。 2. 牛顿—拉斐逊方法 ANSYS程序的方程求解器计算一系列的联立线性方程来预测工程系统的响应。然而，非线性结构的行为不能直接用这样一系列的线性方程表示。需要一系列的带校正的线性近似来求解非线性问题。 一种近似的非线性救求解是将载荷分成一系列的载荷增量，即逐步递增载荷和平衡迭代。它可以在几个载荷步内或者在一个载荷步的几个子步内施加载荷增量。在每一个增量的
Ξ 求解温度场的非线性有限元方法 刘福来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
第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 u M (4.141) 式中，Ku u C I += 为粘性效应项，考虑阻尼、粘塑、粘弹等效应。P 为外部激励。 对于考虑各种非线性效应的动力学问题求解，需要对动力学方程进行直接时间积分。即非线性动力有限元分析具有如下特点：（1）问题分析过程需要考虑时间积分效应，不必做模态分析，不必提取固有频率；（2）采用直接积分方法求解非线性动力学方程，需要对时间作积分计算，因此计算量远远大于线性模态动力学方法；（3）非线性动力学分析中可以施加不同类型的载荷，包括结点力、非零位移、单元载荷；（4）在每个时间步上，进行质量、阻尼、及刚度的集成，采用完整矩阵，不涉及质量矩阵的近似；（5）可以同时考虑几何、材料和接触等多种非线性效应。 非线性动力有限元分析程序常采用隐式Hilber-Hughes-Taylor 法进行时间积分运算。这种方法适于模拟非线性结构的动态问题，对于冲击、地震等激发的结构动态响应以及一些由于塑性或粘性阻尼造成的能量耗散，隐式算法特别有效。隐式积分方法需要对刚度矩阵求逆计算，并通过多次迭代求解增量步平衡方程。隐式Hilber-Hughes-Taylor 时间积分算法为无条件稳定，对时间步长没有特别的限制。 采用子空间法也可以对动力学平衡方程作时间积分运算。子空间法是提取模态分析得到的各阶特征模态，并采用与线性模态动力学分析方法相近的分析方式进行求解。对于带有微小非线性效应的问题，如材料小范围进行入屈服、结点转角不大的情况，子空间法效率比进接积分法要高。
第八章 几何非线性问题的有限元法 引言 前面各章所讨论的问题都是在小变形假设的前提下进行的，即假定物体所发生的位移远小于物体自身的几何尺寸，应变远小于1。在此前提下，建立物体或微元体的平衡条件时可以不考虑物体的位置和形状（简称位形）的变化，因此在分析中不必区别变形前后位形的差别，且应变可用一阶无穷小的线性应变表达。 实际上，上述假设有时是不成立的。即使实际应变可能是小的，且不超过材料的弹性极限，但如果需要精确地确定位移，就必须考虑几何非线性，即平衡方程应该相对于变形后的位置得出，而几何关系应该计及二次项。例如平板大挠度理论中，由于考虑了中面内的薄膜应力，求得的挠度比小挠度理论的结果有显著的减低。再如在结构稳定性问题中，当载荷达到一定数值后，挠度比线性解答予示的结果更剧烈地增加，并且确实存在承载能力随继续变形而减低的现象。在冷却塔、薄壁结构及其它比较细长的结构中，几何非线性分析都显得十分重要。 几何非线性问题可以分为以下几种类型： （1）大位移小应变问题。一般工程结构所遇到的几何非线性问题大多属于这一类。例如高层建筑或高耸构筑物以及大跨度网壳等结构的分析常需要考虑到结构大位移的影响。 （2）大位移大应变问题，如金属压力加工中所遇到的问题就属于这一类型。 （3）结构的变形引起外载荷大小、方向或边界支承条件的变化等。 结构的平衡实际上是在结构发生变形之后达到的，对于几何非线性问题来说，平衡方程必须建立在结构变形之后的状态上。为了描述结构的变形需要设置一定的参考系统。一种做法是让单元的局部坐标系始终固定在结构发生变形之前的位置，以结构变形前的原始位形作为基本的参考位形，这种分析方法称作总体的拉格朗日（Ｌagrange ）列式法；另一种做法是让单元的局部坐标系跟随结构一起发生变位，分析过程中参考位形是不断被更新的，这种分析方法称作更新的拉格朗日列式法。 本章首先对几何非线性问题作一般性讨论，从中导出经典的线性屈曲问题的公式；然后建立平板大挠度问题和壳体的大位移（及大转动）分析的有限方法公式；接着还给出了大应变及大位移的一般公式，最后还详细讨论了杆系结构几何非线性问题的有关公式。在讨论中我们采用总体的拉格朗日列式法，但对杆系结构，为应用方便我们给出了两种列式法的公式。 & 一般性讨论 理论基础 无论是对于何种几何非线性问题，虚功原理总是成立的。由虚功原理，单元的虚功方程可以写成如下的形式 {}{}{}{}0=-???**v e eT e eT F dv δσε （） 其中{}F 为单元节点力向量，{}e *ε为单元的虚应变，{}e *δ为节点虚位移向量。 增量形式的应变一位移关系可表示为 {}[] {}e e d B d δε= （）
非线性有限元分析 1 概述 在科学技术领域，对于许多力学问题和物理问题，人们已经得到了它们所应遵循的基本方程（常微分方程或偏微分方程）和相应的定解条件（边界条件）。但能够用解析方法求出精确解的只是少数方程性质比较简单，并且几何形状相当规则的问题。对于大多数工程实际问题，由于方程的某些特征的非线性性质，或由于求解区域的几何形状比较复杂，则不能得到解析的答案。这类问题的解决通常有两种途径。一是引入简化假设，将方程和几何边界简化为能够处理的情况，从而得到问题在简化状态下的解答。但是这种方法只是在有限的情况下是可行的，因为过多的简化可能导致误差很大甚至是错误的解答。因此人们多年来一直在致力于寻找和发展另一种求解途径和方法——数值解法。特别是五十多年来，随着电子计算机的飞速发展和广泛应用，数值分析方法已成为求解科学技术问题的主要工具。 已经发展的数值分析方法可以分为两大类。一类以有限差分法为代表，主要特点是直接求解基本方程和相应定解条件的近似解。其具体解法是将求解区域划分为网格，然后在网格的结点上用差分方程来近似微分方程，当采用较多结点时，近似解的精度可以得到改善。但是当用于求解几何形状复杂的问题时，有限差分法的精度将降低，甚至发生困难。 另一类数值分析方法是首先建立和原问题基本方程及相应定解条件相等效的积分提法，然后再建立近似解法并求解。如果原问题的方程具有某些特定的性质，则它的等效积分提法可以归结为某个泛函的变分，相应的近似解法实际上就是求解泛函的驻值问题。诸如里兹法，配点法，最小二乘法，伽辽金法，力矩法等都属于这一类方法。但此类方法也只能局限于几何形状规则的问题，原因在于它们都是在整个求解区域上假设近似函数，因此，对于几何形状复杂的问题，不可能建立合乎要求的近似函数。 1960年，R.W.CLOUGH发表了有限单元法的第一篇文献“The Finite Element Method in Plane Stress Analysis”，这同时也标志着有限单元法（FEM）的问世。有限单元法的基本思想是将连续的求解区域离散为一组有限个，且按一定方式相互联接在一起的单元的组合体。由于单元能按不同的联结方式进行组合，且单元本身又可以有不同形状，因此可以模型化几何形状复杂的求解域。并且可以利用在每一个单元假设的近似函数来分片地表示全求解域上待求的未知场函数，从而使一个连续的无限自由度问题变成离散的有限自由度问题。 现已证明，有限单元法是基于变分原理的里兹法的另一种形式，从而使里兹法分析的所有理论基础都适用于有限单元法，确认了有限单元法是处理连续介质问题的一种普遍方法。利用变分原理建立有限元方程和经典里兹法的主要区别是有限单元法假设的近似函数不是在全求解域而是在单元上规定的，而且事先不要求满足任何边界条件，因此可以用来处理很复杂的连续介质问题。 在短短四十余年的时间里，有限单元的分析方法已经迅速地发展为适合于使用各种类型计算机解决复杂工程问题的一种相当普及的方法。如今，有限元广泛地应用于各个学科门类，已经成为工程师和科研人员用于解决实际工程问题，进行科学研究不可或缺的有力工具。有限单元法的应用围已由弹性力学平面问题扩展到空间问题，板壳问题，由静力平衡问题扩展到稳定问题，动力问题和波动问题。分析的对象从弹性材料扩展到塑性，粘弹性，粘塑性和复合材料等，从固体
第28卷第1期 V ol.28 No.1 工 程 力 学 2011年 1 月 Jan. 2011 ENGINEERING MECHANICS 82 ——————————————— 收稿日期：2009-06-19；修改日期：2010-03-11 基金项目：国家科技支撑计划项目(2006BA904B03) 作者简介：*周凌远(1968―)，男，四川成都人，副教授，工学博士，从事桥梁结构行为分析研究(E-mail: zhoulingyuan@https://www.doczj.com/doc/c217630355.html,)； 李 乔(1954―)，男，黑龙江铁力人，教授，工学博士，博导，西南交通大学土木工程学院院长，从事桥梁结构行为分析研究 文章编号：1000-4750(2011)01-0082-05 钢筋混凝土梁非线性有限元分析方法 * 周凌远，李 乔 (西南交通大学土木工程学院，成都 610031) 摘 要：针对钢筋混凝土结构有限元分析中，材料进入非线性阶段后，难以通过梁理论准确描述混凝土截面和钢筋应力状态的问题，提出了基于柔度法和分布式塑性理论的钢筋混凝土梁单元材料非线性方法——网格截面法。这种方法采用平面等参单元将梁单元网格化，由单元轴向积分点位置截面网格积分点的混凝土应力描述单元截面应力分布，同时考虑钢筋对刚度的贡献，并通过对截面网格材料的积分计算积分点位置的截面刚度矩阵，再利用力插值函数和能量原理得到梁单元的柔度矩阵，进而对柔度矩阵求逆计算单元刚度矩阵。通过算例验证该方法在钢筋混凝土承载力分析时的准确性。 关键词：有限元；钢筋混凝土梁；柔度法；网格截面；极限承载力 中图分类号：TU375.1; O241.82 文献标识码：A AN APPROACH OF NONLINEAR FINITE ELEMENT ANALYSIS OF REINFORCED CONCRETE BEAM * ZHOU Ling-yuan , LI Qiao (School of Civil Eng, Southwest Jiaotong University, Chengdu 610031, China) Abstract: A beam element with a meshed section based on distributed plasticity and flexibility theory is presented for the material nonlinear finite element analysis of a reinforced-concrete framed structure, the sections of a concrete beam element are discretized into the plane isotropic components in this formulation, the stress distribution on the sections is described with the stresses at quadrature points in the mesh, the stiffness matrices of the sections are calculated by integration of the stress-strain relations of the material on the meshes and the contribution of the stiffness by reinforcing steel is also counted, the flexibility matrix of the element is formed by integration of section flexibility matrices with force-interpolation functions, and then it is inverted to obtain the element stiffness matrix. Finally, a numerical example of the ultimate load capacity analysis of a reinforced concrete beam illustrates the accuracy of the formulation. Key words: finite element; reinforced concrete beam; flexibility method; meshed section; load capacity 钢筋混凝土结构的整体承载力问题一直为工程界所关注，材料非线性有限元方法是研究这类问题的有效手段，其分析模型主要包括集中塑性铰 法[1]和纤维模型法，1977年，Kang 提出了基于纤维模型的二维梁单元[2]，并运用于预应力混凝土框 架的分析，1993年Izzuddin B A 等提出了三次多项式插值的分布式塑性方法分析空间梁单元[3 ―4] ，通 过对沿梁轴方向两个积分点位置的截面划分监控区域，并假定每个监控区域内的法向应力均匀，得到单元的刚度矩阵和节点力，这样在同一个单元内
1. 轴对称问题的弹塑性分析 流程图 ： 节点号，单刚等各项参数 EN1 存储单元节点号， 局部坐标系转 换为全局坐标 N 打印错误 调用子函数 DEMATR 求[D] 调用子函数 BMATR 求 [B] 切线刚度阵 [EK]=[S][Q1]= · JD ·RN ·H(I1)H(J1) 返回各值 Y 读入单元号， B 矩阵位数，单刚位数，单元 开始 JD<0 [C]=[De ] [B] R=1 N [C]=[Dep][B]
解析解。厚壁筒受内压，采用Mises 屈服准则 经计算知，当t=（）时，材料处于弹塑性交界面。 弹性区为: 塑性区： 交界处有：， 最后解得残余应力为： (7a) 有限元网格信息图：(7b) (8a) (8b) (1) (2) (3) (4) (5) (6)
图1 有限元网格 输入数据文件内容（详细信息见附件）: DATA(1) NNODE MELEM IFU IFW IPF IPR NPP NRM HAC MSF NULOAD EXP NM(1-MELEM) NN NN(1-NNODE) R Z NFU(1-IFU) FU NFW(1-IFW) FW MPQ(1-IPF) NPQ*PQ NPRNRZ(1-IPR) PRNRZ E EMU SSS HH UNLOAD 对理想塑性材料厚壁筒，从初始状态开始，历经加载后完全卸载。这一过程中，厚壁筒内会产生残余应力。沿径向R的残余应力如图2-3 所示。
图 2 径向残余应力 -半径曲线 图 2-3 中分别给出了径向残余应力和切向残余应力随半径的变化， 比较。 从图中可以看出， 程序解和解析解在数值上能够很好的吻合， 大的地方 有少许偏差， 这验证了程序计算结果的正确性。 最大误差发生在径向残余应力达到 10 并且和解析解进行了 只是在径向残余应力最 -0.2 -0.4 -0.6 -0.8 -1 -1.2 -1.4 -1.6 12 14 16 Radius R 18 20 -5 -10 -15 图 3 切向残余应力 -半径曲线
公式号、图号等 第十章 非线性动力有限元法 当机械结构受到较大的外载荷，或受到持续时间较短的冲击载荷作用时，结构会产生过大的变形, 以至于必须考虑结构几何大变形对结构整体刚度及固有频率的影响，即所谓的几何非线性影响。另外, 对于多数非线性动力学问题，还需要考虑材料非线性、接触非线性等方面的影响。 非线性动力学分析求解的基本方程有如下形式 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 几何非线性问题的有限元法 几何非线性问题一般是指物体经历大的刚体位移和转动，但固连于物体坐标系中的应变分量仍假设为小量, 即大位移小应变情况。
非线性有限元 博士研究生专业课课程报告
目录 第一章绪言 (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)
1 §1.2 线性有限元的回顾 线性有限元的理论虽然不是本课程的重点，而线性有限元法是非线性有限元法的基础；非线性有限元法又是线性有限元法的发展。因为非线性问题的求解通常是采用分段线性化的思想，使其成为一系列的线性问题。因此有必要首先回顾线性有限元的一些基本内容。主要是线性有限元的基本理论和方法，以及当前应用最广的等参元理论。 固体力学的理论（材力、弹力、结构力学，无论是线性还是非线性）是建立在本构方程、几何（运动）方程以及平衡方程三方面方程的基础上的，由此导出的控制方程的解是满足上述充分、必要条件的唯一解，而且是反映了结构真实受载后的运动状态（变形）。在结构分析中许多情况下，本构和几何方程可以处理成线性，使控制方程线性化，求解也大大简化，同时可达到工程上的精度要求，这就是线性有限元的基本理论。 §1.2.1 线性本构关系（广义虎克定律） 影响结构材料性能有诸多因素（应力、应变、变形率、温度、湿度、时效等），而通常建立本构方程时仅考虑应力、应变两个物理参数，认为两者成线性关系的经典的弹性理论，即著名的虎克定律。 各向同性材料的Hooke 定律 ij ijkl kl D σε= 或 ij ijkl kl C εσ= ()1其中ijkl D 和ijkl C 分别为弹性阵和柔度阵。 由剪应力互等定理，弹性阵()99ijkl D ?独立材料参数的个数由81个减少为21个。进而对于正交异性的材料参数为9个独立的参数，对于各向同性材料： 01121, 2,322(1) ,e i j i j i j E G G E d dS d i j νν εσδ -+= += = ()2 仅有两个独立的材料常数，即E 和ν，其中E 为弹性模量，ν为泊松比。 §1.2.2线性几何方程（小变形情况） 线性（小变形）关系： ()(){}1 ,,2 ij T U U U i j j i ε+=?= ()3 位移边界条件：u S 边界上 U U = 其中U 为位移向量， U 为边界u S 上的指定位移，()T ?为微分算子。 实际结构可根据它们的几何特点，将三维问题简化为二维问题，这里主要有：平面应力、平面应变和轴对称状态。 1）平面应力（薄壁结构） 外力仅作用在平面内，两表面无外力作用（离心力作用下圆盘）两表面应力分量（且在整个厚度方向上），yz σ=xz σ=zz σ=0，而xx σ、yy σ、xy σ沿厚度均匀分布。 按Hooke 定律，
几何非线性大作业荷载增量法 和弧长法程序设计 系（所）：建筑工程系 学号：1432055 姓名：焦联洪 培养层次：专业硕士 指导老师：吴明儿 2015年6月19日
一、几何非线性大作业（ Newton-Raphson法） 用荷载增量法（Newton-Raphson法）编写几何非线性程序： （1）用平面梁单元，可分析平面杆系 （2）算例：悬臂端作用弯矩。悬臂梁最终变形形成周长为悬臂梁长度的圆。 1.1 Newton-Raphson算法基本思想 图1.1 Newton-Raphson算法基本思想 1.2 悬臂梁参数 基本参数：L=2m, D=0.03m, A=7.069E-4m2, I=3.976E-08m4 ,E=2.0E11N/m2
图1.2 悬臂梁单元信息 将悬臂梁分成10个单元，如图1.2所示 2.1 MATLAB输入信息 材料信息单元信息 约束信息（0为约束，1为放松）荷载信息（FX,FY,M）
节点信息 2.2 求解过程 梁弯成圆形：理论弯矩M=EIY"=24981.944N.m ,直径为0.642m 运用ABAQUS和MATLAB进行求解对比： 图1.3 加载图 图1.4 ABAQUS变形图
图1.5 MATLAB变形曲线 ABAQUS和MATLAB变形对比，最终在理论荷载作用下都弯成了一个圆，其直径为0.64716m，与理论值相对比值为：（0.64716-0.642）/0.642=0.00804.非常接近。 2.3 加载点荷载位移曲线 图1.5 加载点Y方向的荷载位移曲线
加载点的最大竖向位移分别为1.4525m和1.45246m，相对比值（1.4525-1.45246）/1.45246=2.75395E-05。完全相同，说明MATLAB的计算结果很好。
基于非线性有限元法的弹簧刚度分析 摘要本文以铁路车辆三大件式转向架用螺旋弹簧为研究对象。传统的弹簧的垂向和横向刚度分析一般采用经验公式来计算，这在线弹性范围不会存在问题。而实际工作中，弹簧运动过程往往存在大的变形，属于非线性的范畴，所以本文要研究其在非线性范围内有限元计算结果和传统经验公式的对比，以便于指导设计研究。 关键词非线性；有限元；弹簧；横向刚度；垂向刚度 前言 螺旋弹簧在铁路车辆三大件式转向架中起着垂向支撑和减震的双重作用，是三大件式转向架必不可少的组成部分之一，其刚度的大小和匹配关系着整个转向架的动力学性能，因此对弹簧刚度的研究有着非常重要的意义。 弹簧刚度作为弹簧的主要参数之一，在以往的设计中往往是按照经验公式对其轴向刚度和横向刚度进行计算，在线性阶段该方法也许不会有什么问题，可是当弹簧变形到一定程度的时候会出现弹簧自接触的问题，即弹簧由于变形而发生了自身的一部分与另一部分接触，此时的弹簧参数已经由类线性参数变成了非线性参数，而按照经验公式则无法判断何时弹簧进入非线性，所以弹簧的设计仅仅依靠经验公式会存在一定的风险。由于有限元软件的普及[1]，本文将使用有限元的方法对弹簧刚度进行分析，从而更进一步提高刚度计算的精度。 1 研究对象 本次分析使用的模型為某型转向架上的一种弹簧，该弹簧所用材料为60SiMnAT，有效圈数为5.5圈，线径24mm，中径115mm，剪切模量为78.5GPa，自由高252mm。其材料属性如下表。 2 研究方法 按照刚度的定义，即结构抵抗变形的能力，也就是产生单位位移所需要的力，其单位为N/mm。在进行弹簧横向刚度和轴向刚度的分析时，弹簧的两个端面与接触面之间做刚性接触处理，并假定在整个过程中上下支撑面保持平行，对弹簧进行强迫位移分析，并取得每一个位移值对应的支反力，从而求得其刚度曲线。分析采用UGNX软件，分析假想图如下。 3 结果及分析 根据分析结果可以得到如下轴向支反力与轴向位移关系图、轴向刚度与轴向位移关系图、横向刚度与轴向位移关系图等。
第四章材料非线性有限元法 以上三章分别研究了线性弹性有限元法，材料非线性本构方程和非线性方程组解法，本章就可以研究材料非线性有限元法了。 在材料非线性基本方程中，除第二章所述的本构方程外，与线性弹性一样，而非线性有限元法又归结为一系列线性弹性问题。因此，只要在第一章中改用第二章的本构方程，就可建立材料非线性有限元法的基本内容。 §4-1 非线性弹性有限元法 第二章提到，非线性弹性本构方程与形变理论弹塑性本构方程在形式上相同，所以与第二章一样，这里也按塑性力学形变理论，研究非线性弹性有限元法，以便把二者统一起来。 1．非线性弹性基本方程为了便于以后直接引用，这里列出全量形式的非线性弹性（或形变理论弹塑性）基本方程，并用矩阵表示。 几何方程： (1.14) 本构方程： =[D ] (2.13)
平衡方程： (在 内) (1.20) 边界条件： （在A 上）（1.22） （在A 上） (1.23) 虚功方程： (1.28) 位能变分方程： =0 ( 1.31) 其中 (1.32)
(4.1) 2．非线性方程组的建立由于虚功方程本身不涉及材料性质，所以第一章由虚功方程得到的单元平衡方程（1.48）式和总体平衡方程（1.109）式完全适用于非线性弹性（或形变理论弹塑性）问题。可见，只要把非线性弹性（或形变理论弹塑性）本构方程代入单元或总体的平衡方程，就可以建立非线性方程组。 （1）割线刚度方程仿照线性弹性有限元法，把（1.36）式代入（2.13）式后，再把（2.13）代入（1.48）式便得单元割线刚度方程，即 (4.2) 其中单元割线刚度矩阵 (4.3) 而割线本构矩阵[ ] ，如（2.14）式所示。 仿照（1.113）式的推导，同样可得总体割线刚度方程 即
1. 轴对称问题的弹塑性分析 开始 读入单元号，B 矩阵位数，单刚位数，单元节点号，单刚等各项参数 EN1存储单元节点号，局部坐标系转 换为全局坐标 JD<0 打印错误 调用子函数DEMATR 求[D] 调用子函数BMATR 求[B] [S]== 切线刚度阵[EK]=[S][Q1]=·JD ·RN ·H(I1)H(J1) 结束 N Y Y 流程图： R=1 [C]=[De][B] N Y [C]=[Dep][B] 返回各值
解析解。厚壁筒受内压，采用Mises屈服准则 (1) 经计算知，当t=（）时，材料处于弹塑性交界面。 弹性区为: (2) (3) 塑性区： (4) (5) 交界处有：， (6) 最后解得残余应力为： (7a) (7b) (8a) (8b) 有限元网格信息图：
图1 有限元网格 输入数据文件内容(详细信息见附件): DATA(1) NNODE MELEM IFU IFW IPF IPR NPP NRM HAC MSF NULOAD EXP NM(1-MELEM) NN NN(1-NNODE) R Z NFU(1-IFU) FU NFW(1-IFW) FW MPQ(1-IPF) NPQ*PQ NPRNRZ(1-IPR) PRNRZ E EMU SSS HH UNLOAD 对理想塑性材料厚壁筒，从初始状态开始，历经加载后完全卸载。这一过程中，厚壁筒内会产生残余应力。沿径向R的残余应力如图2-3所示。
-1.6 -1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 10 12 14 16 18 20 R a d i a l s t r e s s Radius R 图2 径向残余应力-半径曲线 -15 -10 -5 5 10 12 14 16 18 20 T a n g e n t i a l s t r e s s Radius R 图3 切向残余应力-半径曲线 图2-3中分别给出了径向残余应力和切向残余应力随半径的变化，并且和解析解进行了比较。从图中可以看出，程序解和解析解在数值上能够很好的吻合，只是在径向残余应力最大的地方有少许偏差，这验证了程序计算结果的正确性。最大误差发生在径向残余应力达到
§1.2 线性有限元的回顾 线性有限元的理论虽然不是本课程的重点，而线性有限元法是非线性有限元法的基础；非线性有限元法又是线性有限元法的发展。因为非线性问题的求解通常是采用分段线性化的思想，使其成为一系列的线性问题。因此有必要首先回顾线性有限元的一些基本内容。主要是线性有限元的基本理论和方法，以及当前应用最广的等参元理论。 固体力学的理论（材力、弹力、结构力学，无论是线性还是非线性）是建立在本构方程、几何（运动）方程以及平衡方程三方面方程的基础上的，由此导出的控制方程的解是满足上述充分、必要条件的唯一解，而且是反映了结构真实受载后的运动状态（变形）。在结构分析中许多情况下，本构和几何方程可以处理成线性，使控制方程线性化，求解也大大简化，同时可达到工程上的精度要求，这就是线性有限元的基本理论。 §1.2.1 线性本构关系（广义虎克定律） 影响结构材料性能有诸多因素（应力、应变、变形率、温度、湿度、时效等），而通常建立本构方程时仅考虑应力、应变两个物理参数，认为两者成线性关系的经典的弹性理论，即著名的虎克定律。 各向同性材料的Hooke 定律 ij ijkl kl D σε= 或 ij ijkl kl C εσ= ()1其中ijkl D 和ijkl C 分别为弹性阵和柔度阵。 由剪应力互等定理，弹性阵()99ijkl D ?独立材料参数的个数由81个减少为21个。进而对于正交异性的材料参数为9个独立的参数，对于各向同性材料： 01121, 2,322(1) ,e i j i j i j E G G E d dS d i j νν εσδ -+= += = ()2 仅有两个独立的材料常数，即E 和ν，其中E 为弹性模量，ν为泊松比。 §1.2.2线性几何方程（小变形情况） 线性（小变形）关系： ()(){}1 ,,2 ij T U U U i j j i ε+=?= ()3 位移边界条件：u S 边界上 U U = 其中U 为位移向量， U 为边界u S 上的指定位移，()T ?为微分算子。 实际结构可根据它们的几何特点，将三维问题简化为二维问题，这里主要有：平面应力、平面应变和轴对称状态。 1）平面应力（薄壁结构） 外力仅作用在平面内，两表面无外力作用（离心力作用下圆盘）两表面应力分量（且在整个厚度方向上），yz σ=xz σ=zz σ=0，而xx σ、yy σ、xy σ沿厚度均匀分布。 按Hooke 定律，