有限单元法及程序设计

有限单元法及程序设计


绪 论

1.力学分析方法:解析法,数值法
有限元法——实际结构形状和所受载荷比较复杂,大多用解析法很困难,因而数值法得到不断发展,随着电子计算机的进步,而发展起来的一种新兴的数值分析方法.

2基本步骤：
（1）结构离散化：将结构从集合上用线或面划分为有限个单元。
（2）单元分析： 导出单元的节点位移和结点力之间的关系（单元刚度矩阵）。
（3）整体分析： 将各单元组成的结构整体进行分析，导出征个结构点位移与结点力之间的关系。

3程序设计的步骤：
提出问题，拟定解决方案
构造数学模型
画出程序流程图
编写程序
编译调试程序
试算验证程序


4.根据国家标准（GB-1526-89）规定的程序流程图标准化符号及规定 ：

EMBED PBrush

图表示程序流程图的起点和终点；
图表示数据信息的输入和输出；
图表示数据进行系列运算之前要完成的数据预置；
图表示判断条件；
图表示各种处理功能，如数学运算方式等；
图表示流程的路径和指向。

杆件结构的有限单元法及程序设计

平面杆件单元的有限单元法
有限单元法的基本概念


基本思路：先分后合（先单元分析，再整体分析）
基本概念：整体号：节点端点号按自然数1，2，3，……（在整体坐标系xOy下）
局部号：每一个单元始末用i，j标记（在单元的局部坐标 EMBED Equation.3 系下，方向与整体坐标系一致）。
3．Fe= keδe 其中： ke = EMBED Equation.3 单元刚度矩阵，各元素为刚度系数
δe= EMBED Equation.3 单元杆段位移列阵
Fe= EMBED Equation.3 单元杆端力列阵

K EMBED Equation.3 =P （1-7）
K= EMBED Equation.3 整体刚度矩阵 EMBED Equation.3 = EMBED Equation.3 EMBED Equation.3 位移列阵
P= EMBED Equation.3 节点载荷列阵

3.有限元位移法分析连续梁需要考虑的问题
刚度集成法:
①将(1-3)K扩阶,扩大的元素为0,得到单元贡献矩阵
单元①: K①= EMBED Equation.3 单元② : K②= EMBED Equation.3
②将单元贡献矩阵想叠加,形成整体刚度矩阵
K= K①+ K②= EMBED Equation.3
两端支承条件的引入
先不考虑约束条件,得到整体刚度矩阵后,将其主对角线元素k EMBED Equation.3 改为1,第i行,第j列其余元素改为0,对应的载荷元素也改为0.
非结点荷载的处理
利用等效结点荷载进行分析:
各结点(包括两端结点)加约束,阻止结点转动,其约束力矩分别为交于该结点的各相关单元的固端力矩之和,顺时针为正.
去掉附加约束(相当在各结点施加外力荷载P EMBED Equation.3 ,其大小与

约束力矩相同,方向相反)
将两部分杆端弯矩叠加起来.


第二节 局部坐标系中的单元刚度矩阵
一般单元
设单元 eq \o\ac(○,e) 的弹性模量、截面惯性矩、截面积分别为E、I、A,杆长为l。单元的i、j端各有三个杆端力、 EMBED Equation.3 (即轴力、剪力和弯矩)和与其相应的三个杆端位移 EMBED Equation.3 ，如图1-7所示。图中 EMBED Equation.3 为单元局部坐标系，取i点位于坐标原点， EMBED Equation.3 轴与杆轴重合，规定由i到j为 EMBED Equation.3 轴的正方向，由 EMBED Equation.3 轴顺时针旋转90?为 EMBED Equation.3 轴正方向。力和位移的正方向如图1-7所示。
EMBED PBrush

在此单元中，单元杆端力列阵和杆端位移列阵分别为
EMBED Equation.3 单元杆端力列阵
EMBED Equation.3 杆端位移列阵

为了导出一般单元杆端力与杆端位移之间的关系，我们分别考虑以下两种情况。
首先分析两个杆端轴力 EMBED Equation.3 与轴向位移 EMBED Equation.3 的关系。根据胡克定律，有
EMBED Equation.3 （a）
其次考虑杆端弯矩 EMBED Equation.3 与杆端剪力 EMBED Equation.3 与杆端转角 EMBED Equation.3 和横向位移 EMBED Equation.3 的关系。根据结构力学位移法的转角位移方程，并按照本节规定的符号和正负号，可得
EMBED Equation.3 (b)

将(a)、(b)两式合在一起，并写成矩阵形式如下
EMBED Equation.3 EMBED Equation.3 (1-17)
上式可简写成 EMBED Equation.3 (1-18)
其中单元刚度矩阵为
EMBED Equation.3 （1-19）
单元刚度矩阵的性质
每个元素代表单位杆端位移引起的杆端力,任一元素 EMBED Equation.3 (r、s取1至6)的物理意义是第s个杆端位移分量等于1时,所引起的第r各杆端力分量值.
是对称矩阵,其元素 EMBED Equation.3 .
是奇异矩阵,它的元素行列式等于零,即 EMBED Equation.3 .
具有分快性质.

轴力单元:只考虑轴向杆端位移和杆端力的单元

单元刚度矩阵的坐标变换
上述单元刚度方程和单元刚度矩阵实在局部坐标系 EMBED Equation.3 中建立起来的，对于一般杆件结构，分析时所划分的各单元的局部坐标系显然不同。因此在研究结构平衡条件和变形连续条件时，必须选定一个统一的坐标系xOy,称为整体坐标系。同时，还必须把在局部坐标系中建立的单元刚度矩阵转换为整体坐标系下的单元刚度矩阵。
图1-8a)、图1-8b)分别表示单元 eq \o\ac(○,e) 在局部坐标系 EMBED Equation.3 和整体坐标系xOy种的杆端力分量。
为了导出整体坐标系中杆端力Xi、Yi、Mi和局部坐标系中 EMBED Equation.3 之间的关系，

将Xi、Yi分别向 EMBED Equation.3 轴上投影，可得
EMBED Equation.3 a)
式中，α表示由x轴到 EMBED Equation.3 轴之间的夹角，以顺时针为正。
EMBED PBrush 图 1-8
在两个坐标系中，力偶分量不变，即
EMBED Equation.3 b)
同理，对于单元 eq \o\ac(○,e) j端的杆端力可得
EMBED Equation.3 c)
将a)、b)、c)合起来，并用矩阵形式表示，可得
EMBED Equation.3 （1-24）

此式即为两种坐标系中单元杆端力的变换式，亦可简写为
EMBED Equation.3 （1-25）

式中： EMBED Equation.3 局部坐标系中的单元杆端力列阵
EMBED Equation.3 整体坐标系中的单元杆端力列阵
EMBED Equation.3 单元坐标变换矩阵(1-26)
T为正交矩阵,其逆矩阵等于其转置矩阵
即 EMBED Equation.3

单元未知量编码

为了便于编程计算，需要按一定规律对结点的位移分量编号。结构的节点位移有自由结点位移和支座结点位移（亦称支座结点位移）之分。自由结点位移是未知量。建立结构整体结构方程求解未知节点位移的方式有两种：“前处理法”和“后处理法”。
EMBED PBrush
用后处理法分析改结构时，设所有点位移都是未知量，则结点位移列阵为（参看图1-9））
EMBED Equation.3
Pix、Piy、Piθ分别代表作用在结点i（i=1，2，3，4）上的水平力、竖向力和力偶。规定，结点力Pix、Piy的正方向与整体坐标系x、y的正方向相同，Piθ以顺时针指向为正；结点位移的正方向与结点力的正方向一致。
在求出各单元刚度方程之后，根据结点平衡条件和位移连续条件，可建立整个结构的位移法方程
EMBED Equation.3 （1-37）
或 EMBED Equation.3 （1-38）
其中
EMBED Equation.3 （1-39）
为结构的整体刚度矩阵，或称为结构的原始刚度矩阵。
在建立整体刚度方程式（1-38）时，假定所有点位移都是未知量，相当整体结构无知座，因而在外力作用下，除了弹性变形外，还有可能发生刚体位移，此时，各结点位移不能唯一确定。这说明式（1-39）为奇异矩阵，不能求逆，故利用（1-38）不能求结点位移。
实际在图2-9a）所示的刚架中，结点1、4为固定点，因此结点位移是已知的，支承条件全为零。将该支承条件引入到整体刚度方程，得
EMBED Equation.3 （1-40）
可以分为两组方程，一组是
EMBED Equation.3 （1-41）
它可以求结未知结点位移Δ2、Δ3。令一组是
EMBED Equation.3 （1-42）
称为反力方程。利用式（1-41）求出结点位移Δ2、Δ3并代入上式后，便可计算未知的

支座反力。
对于一般杆件结构，都可以按上述步骤进行分析。无论结构有多少个结点位移分量，经过调整其排列顺序，总可以将它分为两组：一组包括所有的未知结点的位移分量，以ΔF表示；另一组为支座结点位移分量，以ΔR表示。相应的，将全部结点分为两组，与ΔF相应者为已知的结点力列阵，以PF表示；与ΔR相应者为支座结点力列阵，以PR表示。于是有
EMBED Equation.3
与以上分析方法相配合，将整体刚度矩阵K0中的各元素重新排列，则K0Δ0=P0可写成
EMBED Equation.3 （1-43）
展开上式得
EMBED Equation.3 式（1-44）、式（1-45）
式（1-44）为“修正的整体刚度方程”，它与式（1-38）的区别在于引进了支承条件。
后处理法: 由单元刚度矩阵形成整体刚度矩阵,建立刚度方程后在引入支承条件,进而求解结点的位置位移的方法.
前处理法: 仅对未知的自由结点位移分量编号,得到的结点位移列阵中不包含已知的约束结点位移分量.
图1-10所示的具有组合结点的刚架划分
为三个单元，其编号为①、②、③，各
杆之间的箭头表示局部坐标系的正方向，
刚架结构编号为1～5。下面考虑各单元
结点位移分量编号。采用“先处理法”
作如下规定：
仅对独立的位移分量按自然数编号，
称为位移号。若某些位移分量由于联结
条件或直接杆轴向刚性条件的限制彼此
相等，则编号相同。
在支座处，由于刚性约束而使某些位移分量为零时，此位移分量编号为零。
因此图1-10编号如下
杆数 单元编号 单元结点编号 单元位移分量编号 始端 末端 AB ① 1 2 000 123 BC ② 2 3 123 456 DC ③ 5 4 000 457 在计算程序中，单元两端结点号可采用二维数组JE(i,e)表示，称为“单元两端结点号数组”。
JE(1,e)=单元始端的结点号
JE(2,e)=单元末端的结点号
在本例中
JE(1，1) =1， JE(2，1) =2
JE(1，2) =2， JE(2，2) =3
JE(1，3) =5， JE(2，3) =4
任意结点位移分量的位移号可用二维数组JN(i,j)表示，称为“结点位移号数组”
JN（1，j）= 结点j沿x方向的位移号
JN（2，j）= 结点j沿y方向的位移号
JN（2，j）= 结点j角位移的位移号
在本例中，对第三结点而言，
JN（1，3）= 4
JN（2，3）= 5
JN（3，3）= 6
将单元 eq \o\ac(○,e) 始端及末端得位移号排成一行（始端在前），此数码为“单元定位数组”，利用它可方便的形成杆端位移与相应结点位移间的协调条件。它的展开式为
EMBED Equation.3
式中，d个元素m1～md分别是单元 eq \o\ac(○,e) 的两端位移分量所对应的位移号数值。在本例中，
EMBED Equation.3

平面结构的整体刚度矩阵

在进行了

单元分析得出单元刚度矩阵之后，需要进行整体分析。以“先处理法”为例，将离散单元重新组合成原结构，使其满足结构结点的位移连续条件和力的平衡条件，从而得到修正的结构刚度方程，即前面给出的式（1-44）
KFF△F=PF
式中：KFF称为修正的结构整体刚度矩阵；△F、PF分别为自由结点位移与自由结点荷载列阵。当已知计算对象为自由结点位移分量而不至引起误解时，式（1-44）也常称为整体刚度方程，KFF简称为整体刚度矩阵，△F、PF分别简称为结点位移列阵与荷载列阵。为了书写方便，下表常常略去。
有单元刚度矩阵集成整体刚度矩阵，通常采用“直接刚度法”，把计算步骤分为两步，首先求出各单元贡献矩阵，然后将它们叠加起来，得出整体刚度矩阵。然而在实际电算中，不便采用，原因是在计算中需要先将所有单元的贡献矩阵Ke都保存起来，Ke的阶数与整体刚度矩阵K相同，这就占用了大量存储容量，因此实际运算中，采用“边定位，便累加”
的方法。其原理没有变，而且结果相同。
EMBED PBrush
图1-10中单元①的单元矩阵为
EMBED Equation.3
式中单元刚度矩阵的上面和右侧标记了单元结点位移分量编号。因为整体刚度矩阵个元素是按位移分量编号排列的，按先处理法，单元刚度矩阵中对应于分量编号为零的元素不进入整体刚度矩阵，非零编号指明了其余各元素在整体刚度矩正中的行、列号。
所以
EMBED Equation.3 EMBED Equation.3 EMBED Equation.3
单元②的刚度矩阵为
EMBED Equation.3
各单元在K中的位置为
EMBED Equation.3
单元③的刚度矩阵为
EMBED Equation.3
各元素在K中的位置为
EMBED Equation.3 EMBED Equation.3 EMBED Equation.3
按以上定位方法，将三个单元的刚度矩阵有关元素一道整体刚度矩阵对应位置，得到






EMBED Equation.3
在世机电算时，采用
将K置零,这是K=0 EMBED Equation.3 ;
将k EQ EQ EQ \S(①) 中的相关元素,按照“对号入座”，累加到K；
将k EQ \S(②) 中的相关元素，继续按照“对号入座”，累加到K；
将k EQ \S(③) 中的相关元素，继续按照“对号入座”，累加到K，整体刚度矩阵最后完成。
求图1-12所示刚架的整体刚度矩阵K。
设各杆截面尺寸相同。
A=0.5m2 I=1/24 m2
E=3×104MPa
解
整理数据并进行编号。
EMBED Equation.3
求局部坐标系中单元刚度矩阵 EMBED Equation.3 。由于单元
①、②的尺寸完全相同，故有 EMBED Equation.3 ，可直接利用公式
（1-19）求得
EMBED Equation.3
（3）求整体坐标系中单元刚度矩阵 EMBED Equation.3 。按式（1-32）和

式（1-33），求得各单元在整体坐标系中的单元刚度矩阵，并将单元结点位移分量编号标记于上面与左侧。
EMBED Equation.3
EMBED Equation.3
形成整体刚度矩阵K。采用本节中介绍的方法建立结构整体刚度矩阵如下
EMBED Equation.3
EMBED PBrush


第六节 非结点荷载处理

为分析平面结构而建立的整体刚度方程，反映了结构的结点荷载与结点位移之间的关系。作用在结构上的荷载除了直接作用在结点上的荷载Pd之外，还有作用在杆件上的分布荷载等。这些非结点荷载应转换成等效结点荷载Pe。将Pd和Pe叠加，课的综合结点荷载（总结点荷载）Pc,其下标c通常可略去不写，即
P=Pd+Pe (1-46)
直接作用在结点上的荷载，可按其作用方位直接加入P之中，而等效结点荷载的计算步骤如下。
第一步：在局部坐标系下，求单元 eq \o\ac(○,e) 的固端力 EMBED Equation.3 。
EMBED Equation.3 （1-47）
式中，子向量 EMBED Equation.3 分别为单元 eq \o\ac(○,e) 在端点i、j的固端内力。几种非结点荷载作用下的单元固端力列于表1-2
表1-2 两端固定梁的固端力
简 图 剪 力 弯 矩 QAB QBA MAB MBA EMBED PBrush EMBED Equation.3 EMBED Equation.3 EMBED Equation.3 EMBED Equation.3 EMBED PBrush EMBED Equation.3 EMBED Equation.3 EMBED Equation.3 EMBED Equation.3 EMBED PBrush EMBED Equation.3 EMBED Equation.3 EMBED Equation.3 EMBED Equation.3
第二步：求单元 eq \o\ac(○,e) 的等效结点荷载 EMBED Equation.3 。
仿照局部坐标系与整体坐标系中单元杆端力的变换式（1-25）
EMBED Equation.3
固端内力在两种坐标系下的变换公式，可写成
EMBED Equation.3
EMBED Equation.3 （1-48）
因此，等效结点荷载列阵 EMBED Equation.3 可以由下式求出
EMBED Equation.3 （1-49）
将该式展开得
EMBED Equation.3 （1-50）
当α=0时， EMBED Equation.3
第三步：求整体结构的等效结点荷载 EMBED Equation.3
求得单元等效结点荷载 EMBED Equation.3 之后，利用单元结点位移分量编号，就可以将 EMBED Equation.3 中的各分量叠加到结构等效荷载列阵 EMBED Equation.3 中去。因为 EMBED Equation.3 中的各元素是按结点位移分量编号排列的， EMBED Equation.3 中的6个元素也与结点位移分量编号一一对应，所以也按对号入座方法，将其之一累加到 EMBED Equation.3 中相应的位置上去。当直接作用在结点上的荷载等于零时，即 EMBED Equation.3 时，由式（1-46）可知 EMBED Equation.3
在“后处理法”中，P、Pd和Pe应分别与自由

解点位移相对应，其表达式与式（1-46）相同。
EMBED PBrush











***[JimiSoft: Unregistered Software ONLY Convert Part Of File! Read Help To Know How To Register.]***