有限元-第1章

第二部分
有限单元法
第一章 有限单元法概述 § 1-1 引言
有限单元法作为固体力学的一种分析方法是在本世纪五十年代起源于航空工程中飞 机结构的矩阵分析方法。结构矩阵分析是把一个结构看成由许多元件互相连接而组合成的 集合体，通过对元件的受力分析，建立节点位移与节点力之间的关系，再将这些关系集合 起来形成结构方程组。根据选取的基本未知量是节点位移还是节点力，有位移法、力法和 混合法。结构矩阵分析的对象限于由杆、梁、受剪板等元件组成的结构。 1960 年 Clough R. W. 等人将这种处理问题的方法推广用来解弹性力学的平面应力问 题，并第一次采用“有限单元法”这个术语。应用有限元法对任意连续体进行分析时，首 先将连续体划分成有限个单元，并在每个单元上指定有限个节点，认为相邻单元在节点处 相互连接构成一组单元的集合体。用以模拟或逼近原来的连续体。然后，由对单元的分析 和集合，得到描述该离散结构的代数方程组。 常规的结构矩阵分析法是将每个元件的力与位移之间的关系精确推导出来。而将连续 体离散为单元的有限单元法，是选定场函数的节点值，例如取节点位移作为基本未知量， 对于每个单元根据分区近似的思想，在单元内假设近似的位移插值函数，利用弹性力学的 变分原理建立节点力与节点位移之间的关系，得到一组以节点位移为未知量的方程组。有 限单元法是一种近似的数值方法，显然，如果单元满足问题的收敛性要求，那么随着求解 区域内单元数目的增加，单元尺寸的缩小，近似解将收敛于精确解。 从有限单元法所依据的变分原理来看，早期的有限元法大都依据最小位能原理，以位 移作为基本位知量。后来，有依据最小余能原理的有限元法，以内力作为基本未知量。再 后来则有许多某种形式的广义变分原理，同时将位移和内力作为独立的基本未知量。还有 将这些变分原理结合起来应用，例如在每个单元内用最小余能原理而对整个系统用最小位 能原理求近似解，这就是所谓杂交应力有限元法的基本思想。基于最小位能原理的有限元 位移法是用得最广的一种方法，本教材在第四章介绍混合杂交有限元法的基本概念和基本 理论外，其余各章都采用有限元法的基本理论和方法。
[B] = [L][N ] —称单元的几何矩阵， [L] —微分算子矩阵，例如对于二维问题
∂ ∂x [L] = 0 ∂ ∂y 0 ∂ ∂y ∂ ∂x
2、利用物理方程由式（1-5）导出用单元节点位移表示单元应力的关系式 { σ }=[D] {ε}=[D] [B] {δe} 式中 { σ }——单元内任一点的应力列阵 [D]——与单元材料有关的弹模矩阵。对于二维平面应力问题
e e e
（1-7）
式中 K e 称单元刚度矩阵
[ ]
[K ] = ∫∫∫ [B] [D][B]dv
e T VE
（1-8）
由以上三部分内容导出了单元刚度矩阵的普遍表达式。如果计及单元的体力和考虑到 单元的边界有面力作用的情况，还可求得单元等效节点载荷。
四、建立结构的平衡方程 集合所有单元的平衡方程建立整个结构的平衡方程。这个集合过程包括两方面的内 容：一是将各单元刚度矩阵集合成整个物体的结构刚度矩阵；二是将作用于各单元的等效 节点载荷列阵集合成结构的节点载荷列阵。最常用的集合刚度矩阵的方法是直接刚度法， 将单元刚度矩阵中的各元素按照结构节点自由度排列顺序“对号入座”叠加成结构刚度矩 阵。下面用变分原理导出结构平衡方程的有限元列式。 一个结构的离散化计算模型由许多单元集合而成，结构的总位能等于所有单元的位能 之和，即
[ ]
{ }
[ ]
{ }
式中 [K ] = ∑ K e —— 结构刚度矩阵
{P} = ∑ {P e }——结构的节点载荷列阵
结构处于平衡，根据最小位能原理，位能的一阶变分为零， δΠ = 0 ，得
[K ]{∆} = {P}
这就是结构的平衡方程，也称结构的刚度方程。
（1-11）
上述单元刚度矩阵，单元节点载荷列阵都是相对于整体统一坐标的。若单元刚度矩阵 和单元载荷列阵是对单元局部坐标求得的，必须将它们转换到统一的整体坐标下才能进行 单元的集合。应用坐标变换矩阵 [λ ] ，有
T T e e T e


Ve
由于单元处于平衡，根据最小位能原理 δΠ e = 0 ，得
∂Π e ∂δ
{ }
e T
e T δ − se = 0 [ ] [ ][ ] = B D B dv ∫∫∫ V
{ } { }
得单元节点力与节点位移之间的关系
{s } = {K }{δ }
{P }—作用单元节点的集中载荷列阵。
e 0
Ve —单元体积。
se p ——单元的力边界。
将式（1-4） 、 （1-5） 、 （1-6）代入（1-9） ，并利用式（1-8）得 T T 1 Π = ∑ δ e K e δ e − ∑ δ e Pe 2 (1-10)
{ } [ ]{ }
{ }{ }
式中 P e = PVe + Pse + P0e
L K 1n u1 L K 2n u 2 L L M L K in u i L L M L K nn u n
为了看出该方程能给出所需的结果，取出该方程组中第 i 个方程式
10 30 u i* = K i1u1 + K i 2 u 2 + L + 10 30 u i + L + K in u n
Π = ∑Πe = ∑
1 T {u}T f dv + {u}T P ds + δ e { } { } − ε σ dv ∑ ∫∫∫ ∫∫ ∫∫∫ 2 Ve se p Ve
{}
{ } { } {P }
T —单元内位移函数列阵。
{f }—单元体力列阵； {P}—单元面力列阵。
{ } { } { } { }
e V
——单元的等效节点载荷列阵， PVe , Pse , P0e 分别是
{ }{ }{ }
e s
单元的体力、面力、节点上的集中力的等效节点载荷列阵，
{P } = ∫∫∫[N ] {f }dv
T Ve
{P } = ∫∫ [N ] {P}ds
T se p
(1-10a)
方程（1-10）表示了结构的总位能是各单元的应变能之和加上各单元等效节点载荷的 位能之和。将单元刚度矩阵 K e 和单元节点载荷列阵 P e 按结构节点位移列阵 {∆} 的自由 度数和排列顺序添零升阶，式（1-10）可进一步完成 1 1 T T T T Π = {∆} ∑ K e {∆} − {∆} ∑ P e = {∆} [K ]{∆} − {∆} {P} 2 2
P1 P2 M 30 * 10 u i M Pn
K 11 K 21 L = K i1 L K n1
K 12 K 22 L K i2 L K n2
L
K 1i
L K 2i L L L 10 30 L L L K ni
{δ }= [λ ]{δ } 或 {δ } = [λ ] {δ } [K ] = [λ ] [K ][λ ] （1-12） {P } = [λ ] {P } 式中 { δ }、 {K }、 [P ]是相对各单元局部坐标的单元节点位移列阵、单元刚度矩阵、
e e e T e e T e e T e e e e
1 E [D] = µ 1− µ 2 0
(1-6)
µ
1 0
0 0 1− µ 2
3、利用变分原理建立单元节点力与节点位移之间的关系，即单元的平衡关系。 位能泛函的表达式
Π = ∫∫∫ A(eij ) − f i u i dv − ∫∫ p i u i ds
§ 1-1 有限单元法进行结构分析的步骤及有限元列式 有限元位移法进行结构分析的过程概括起来可以分成以下几个步骤。 一、结构的离散化 所谓离散化，就是将连续体结构划分成有限个单元，并在单元的指定点设置节点，使 相邻单元的有关参数具有一定的连续性，构成单元的集合体，用来代替原来的结构。为了 使该离散的计算模型有效地逼近实际的连续体，就需要考虑合理地选择单元的形状和划分 网格的方案，确定单元和节点的数目、节点的位置和自由度等问题 。 1、单元的几何形状。选择单元的形状既要考虑单元的精度和计算的简便，又要考虑对边 界形状的适应性。如对二维问题，单元的形状有三角形、矩形、任意四边形……；对 于三维问题，有三棱锥的四面体、三棱柱和五面体、长方体……等。 2、单元的划分和节点的布置。节点的布置与单元的划分是相互联系的。通常集中载荷的 作用点、分布载荷集度的突变点、支承点、边界的转折点等都应取作节点，如果物体 的厚度有突变或者物体由不同材料组成，在布置节点时，不要把厚度不同或材料不同 的区域划在同一个单元中。至于节点的多少和分布的疏密、也即单元的大小，要根据 对计算精度的要求和计算机容量等因数综合考虑，单元划分得越细计算精度越高，但 也使计算时间及准备工作增加。因此，在保证足够精度下力求减少单元的数目。故在 划分单元时，对应力变化急剧的区域单元网格应划分得细些，应力变化平缓的区域单 元可划分粗些。单元节点的设置应使单元的几何形状能够容易地由节点位置坐标来描 述。一些在单元的端点或角点上设置节点，还可在角点之间的边上增设节点，这样可 描述单元边界是曲线的情况，用二次或更高次的曲线来描述单元的几何特征。 3、节点自由度数的确定，节点的自由度数应能正确地模拟单元的力学特征和易于构造单 元的近似解，保证相邻单元在公共边界上的连续性。 综合起来，结构离散化应符合几何近似和物理近似的原则。有限元的计算模型在几何 上应近似于真实结构。各单元的物理特性和受力特性应与真实结构在单元所在区域的情况 近似。 二、选择位移模式 在对典型单元进行特性分析时，为了能用节点位移表示单元内各点的位移、应变、应
单元节点载荷列阵。
五、 约束处理 由于结构刚度矩阵是奇异的、不能求逆。造成刚度矩阵奇异的原因是在建立刚度矩阵 时，解除了结构的外界约束而成了自由结构，使结构可产生刚体运动。因此必须排除结构 的刚体位移，使结构刚度矩阵成为非奇异的，才能求解出节点位移。一般情况下，引进边 界位移约束条件后可排除刚体运动，否则，还应适当给定某些节点的位移值。经引进边界 约束条件后，可适当地减少待求的节点位移的数目和方程的数目。 在这里介绍一种适合计算机实施的约束处理的作法，称“置大数法” 。设给定的节点 位移为 ui= u i* ,只需在 ui 所在行中， 将结构刚度矩阵的主对角元素 K ij 置入一个大数， 如 10 30 , 同时，将对应行的载荷项 Pi 用 10 30 u i* 代替,于是结构刚度方程成为
− V sp
[
]
若将单元视为一个独立的弹性力学问题，不计体力，单元只作用节点力 s e 的情况，对于 线形弹性体，单元的位能表示为
{ }
Πe =
1 {ε }T {σ }dv − δ e ∫∫∫ 2 VE
{ } {s }
T e
将式（1-5）和（1-6）代入上式，得
Πe =
1 e δ 2
{ } [B] [D][B]dv {δ }− {δ } {s } ∫∫∫
三、单元的力学特性分析。建立单元刚度矩阵，包括下面三部分内容： 1、利用几何方程由位移模式（1-4）导出单元的应变与单元节点位移的关系
{ε } = [L]{u} = [L][N ]{δ e } = [B]{δ e }
式中 {ε }—单元内的应变列阵。对于二维问题
（1-5）
{ε } = {ε x ε y γ xy }T
两边同除以 10 30 ,由于 10 30 ≥≥ K ij ，显然得到 ui= u i* 。 “置大数法”是处理边界约束的近似方法，但是当这个大数非常大时，将不影响计算 结构的精度。它的优点是保持了原有刚度矩阵的编排顺序，仅对其中某些项进行替换，这 在计算程序中是十分简便的。
力，必须对单元内的位移分布作出一定的假设，也就是在每个单元内假设位移是坐标的某 种简单函数，将单元内的位移场表示成单元节点位移的插值函数形式，称为位移模式。用 矩阵形式写成
{u} = [N ]{δ e }
式中 {u}—单元内位移函数列阵
（1-4）
{δ }—单元节点位移列阵
e
[N ]—称形状函数矩阵，它的元素是位置坐标的函数，称位移插值函数或形函数。