1.1 弹性力学与有限元
1.1.1 弹性力学中的基本概念
在本节中,我们主要介绍弹性力学中的基本概念,包括体力、面力、应力、应变、位移、主应力、相当应力和主应变。
1.体力
体力是一种外载荷,是随体积分布的力,如重力和惯性力。当材料具有磁性或者分布有非自由电荷,这时磁力和静电力也是体力。体力的单位为N/m3。体力在三个坐标轴上的投影为X、Y、Z,它们的总体可用体力列阵表示:。
2.面力
面力也是一种外载荷,是作用在物体表面的力,如接触力和流体压力,其单位为N/m2。面力在三个坐标轴上的投影为、、,它们的总体可用面力列阵表示:。集中力也是一种面力,它作用在物体表面,忽略其作用面积,认定只作用在一点,单位为N。
3.应力
应力是单位截面面积的内力(或内力的分布集度),是表示物体内某位置、沿某截面分布内力的大小和方向的物理量。物体由于外力或湿度、温度改变,其内部将产生内力。
在物体内的同一点,不同方向截面上的应力是不同的。过某点,各截面上应力的大小和方向的总和称为该点的应力状态。为了研究某点的应力状态,围绕该点取一个微小单元体,通常用与坐标面平行的平面截出微小的平行六面体,如图1-1所示,单元体三个方向上的尺寸都是非常小的,分别是dx、dy、dz。单元体每个面上的应力分解为一个正应力、两个切应力,分别与坐标轴平行。正应力的作用面用下标表示,如法线平行于x轴的面上正应力记作;切应力有两个下标,第一个下标为作用面,第二个下标为切应力方向,如表示其作用面垂直于x轴,方向平行于y轴。正应力与切应力正负号规定如下:正面上与坐标轴正方向一致为正,相反为负;负面上与坐标轴正方向一致为负,相反为正。正应力也简单记作拉为正、压为负。图1-1的单元体上,各面上的应力分量均为正。根据切应力互等定理,6个切应力有3组互等关系,即,,。
图1-1 单元体及其各面上的应力分量
物体中某点的应力状态完全可以由、、、、、6个应力分量确定,用应力阵列表示为。
4.应变
应变反映弹性体在外力作用下,其内部各部分发生变形的程度。变形可以归结为长度的改变与角度的改变,故有线应变与角应变。各线段单位长度的伸缩称为正应变或线应变;各线段之间直角的改变,用弧度表示,称为切应变或剪应变。正应变用表示,切应变用表示,用下标表示应变方向,如表示x方向上的线应变,表示x、y两方向的线段之间的直角改变量。正应变以伸长为正,缩短为负;切应变以直角变小为正,变大为负。正应变与切应变都是无量纲的量。由切应力互等和胡克定律,得到切应变也是两两互等的,即,,。
物体中某点的变形状态可用、、、、、这6个应变分量确定,用应变列阵表示为。
5.位移
位移是指物体受力过程中,物体上各点位置的变化量(位置的移动)。物体内一点(微元体)的位移由刚性位移和弹性位移两部分组成,刚性位移是由其他点的形变引起的位移,弹性位移是本身弹性变形产生的位移,与应变有着确定的几何关系。位移是一个矢量,用表示,它在空间直角坐标系中,三个坐标方向的位移分量用u、v、w表示。位移分量以沿坐标轴正方向为正,沿坐标轴负方向为负。位移及其分量的单位为m,用一个位移列阵表示为。
6.主应力
主应力指的是物体内某一点的微面积元上剪应力为零时的法向应力。这时,法向量n=(n1,n2,n3)的方向称为这一点的应力主方向。
7.相当应力
弹性体在外力作用下是否会破坏要通过应力来判断,有限元计算的直接应力结果是各点的6个应力分量,其中3个主应力是判断该点材料是否破坏的主要参数。对于不同的失效形式,适用不同的强度理论。根据这些强度理论求得某点的相当应力,用以判断该点强度是否足够。下面介绍几种常用的强度理论及其相当应力。
(1)第一强度理论(最大拉应力理论)。当材料发生断裂且受力弹性体内的某点有拉应力存在,即σ1大于零时,可以按照该点最大拉应力σ1是否小于许用应力来判断该点强度是否足够。第一强度理论的相当应力为
(1-1)
(2)第二强度理论(最大拉应变理论)。当材料发生断裂且受力弹性体内的某点没有拉应力存在,即σ1小于零时,可以按照该点最大拉应变是否小于许用值来判断该点强度是否足够。经过变换得到用主应力表示的第二强度理论的相当应力为
(1-2)
(3)第三强度理论(最大切应力理论)。当材料发生屈服,受力弹性体内的某点最大切应力大于某一定值时,根据第三强度理论,经过变换得到用主应力表示的相当应力为
(1-3)
(4)第四强度理论(最大形状改变比能理论)。当材料发生屈服,受力弹性体内的某点形状改变比能大于某一定值时,根据第四强度理论,经过变换得到用主应力表示的相当应力为
(1-4)
(5)莫尔强度理论。对于一些材料,如铸铁、混凝土等,它们的抗拉能力和抗压能力不同,当它们受到剪切作用、发生剪切破坏时,不仅与切应力大小有关,还与剪切面上的正应力有关,遵从莫尔强度理论,其相当应力为
(1-5)
根据这些强度理论,只要判断某点的相当应力是否小于相应材料的许用应力即可判断该点强度是否足够。
8.主应变
由单元体6个应变分量εx、εy、εz、γxy、γyz、γzx,可以求出过该点任意方向的线应变和任意两线段之间角度的改变,表达式如下。
(1-6)
(1-7)
其中,l、m、n为过物体内某点P的线段PN的方向余弦,l1、m1、n1为过P点与PN成θ角的线段PN1的方向余弦,θ'为物体受力变形后线段PN与PN1的夹角,如图1-2所示。
图1-2 过物体内某点P的线段PN和PN1
进一步分析还可知,物体内任意一点,一定存在3个相互垂直的应变主向,这3个方向的应变称为主应变。3个主应变中最大的一个就是该点的最大线应变,3个主应变中最小的一个就是该点的最小线应变。3个应变主方向与3个应力主方向是重合的,在线弹性范围内,主应力、主应变服从胡克定律(见后续的物理方程)。
1.1.2 弹性力学的基本方程
弹性力学研究弹性体受外力作用或由于温度变化等原因而发生的应力、应变和位移。弹性体占有三维空间,描述弹性体受力和变形的应力、应变、位移等物理量都是三维坐标的函数。
弹性力学基本方程的导出可从三方面分析:静力学方面,建立应力、体力和面力之间的关系;几何学方面,建立应变、位移和边界位移之间的关系;物理学方面,建立应变与应力之间的关系。通过分析分别得到平衡微分方程、几何方程和物理方程,统称为弹性力学基本方程。
1.平衡微分方程
围绕物体内任意一点,取如图1-1所示的一个微小平行六面体,它的3组面平行于3个坐标面,各边长度都是微量dx、dy、dz。外力作用下物体处于静力平衡状态,物体内任意一点也处于静力平衡状态,单元体各面上所受应力及单元体受到的体力满足平衡方程。3个力的平衡方程为
(1-8)
2.几何方程
如图1-3所示,可以列出平面问题中的几何方程:
(1-9)
同理可得空间问题的几何方程:
(1-10)
图1-3 平面应变与位移
3.物理方程
物理方程体现了应力与应变的关系,也称为本构方程,即
(1-11)
其中,为体积应变。
上述各方程均可用矩阵方程表示,如式(1-11)可用矩阵方程表示为
(1-12)
简写成σ=Dε。其中,D称为弹性矩阵,它完全由弹性常数E和泊松比μ决定。
4.边界条件
弹性力学基本方程共15个,由于平衡方程和几何方程都是微分方程,求解定解还需要边界条件。根据边界条件的不同,弹性力学问题分为位移边界问题、应力边界问题和混合边界问题。
在位移边界问题中,物体在全部边界上的位移是已知的,即
其中,、、在边界上是坐标的已知函数,这就是位移边界条件。
在应力边界问题中,物体在全部边界上的面力分量是已知的。根据面力分量和应力分量之间的关系,可以把面力已知的条件转换成应力方面的已知条件,这就是所谓的应力边界条件,即
(1-13)
其中,面力分量、、在边界上是坐标的已知函数,l、m、n为边界面外法线方向的方向余弦。
在混合边界问题中,物体的一部分边界具有已知位移,即具有位移边界条件,另一部分边界则具有已知面力,即具有应力边界条件。例如,图1-4(a)所示的固定铰支和可动铰支处为位移边界条件,DC边界上分布面力大小为q,其他边界上应力为零,为面力边界条件,整个问题为混合边界问题。图1-4(b)中同一边界存在两种边界条件:x方向位移u|x=a=0;y方向切应力τxy|x=a=0。
图1-4 混合边界问题
1.1.3 平面问题的基本理论
任何一个弹性体都是一个空间物体,任何一个实际的弹性力学问题都是空间问题。如果研究的弹性体具有某种特殊形状,并且所受的外力满足一定的条件,就可以把空间问题简化成平面问题。这样处理,分析和计算的工作量将大大减少,而所得的结果仍能满足工程精度要求。
1.平面应力问题
几何特征:一个方向尺寸比另外两个方向尺寸小得多,即t << a,t << b。
受力特征:假设有很薄的等厚平板,在板边上受有平行于板面且不沿厚度变化的面力,同时体力也平行于板面且不沿厚度变化,此类问题就是平面应力问题。
设薄板厚度为t,以薄板的中面(平分板厚的平面)为xy面,z轴垂直于中面,如图1-5所示,因为板面上不受力,所以有
图1-5 平面应力问题
因为板很薄,外力不沿厚度变化,故可以认为整个薄板的所有点都有
由切应力互等关系得,。这样只剩下平行于xy面的3个应力分量、、非零,所以称为平面应力问题。
因为板很薄,3个应力分量、3个应变分量和两个位移分量都可以认为不沿厚度变化,即它们只是x和y的函数,与z无关。
2.平面应变问题
几何特征:一个方向尺寸比另外两个方向尺寸大得多,且沿该方向截面尺寸和形状不变。
受力特征:设有很长的柱形体,在柱侧面上受有平行于横截面且不沿长度变化的面力,同时体力也平行于横截面且不沿长度变化,此类问题就是平面应变问题。
假想该柱体无限长,以任意一横截面为xy面,z轴垂直于xy面,如图1-6所示,则所有应力分量、应变分量和位移分量都不沿z方向变化,它们只是x和y的函数。此外,在这一情况下,由于柱体无限长,任意一横截面都可看作对称面,所有点都只会沿x和y方向移动,而不会有z方向位移,即w=0,εz=γzx=γzy=0,不为零的应变分量只有εx、εy、γxy,所以称为平面应变问题。
图1-6 平面应变问题
3.平面问题的基本方程
平面应力问题和平面应变问题都只有8个独立的未知量σx、σy、τxy、εx、εy、γxy、u、v,它们只是x和y的函数,因此统称平面问题。
平面问题的平衡微分方程为
(1-14)
平面问题中的几何方程为
(1-15)
平面应力问题中的物理方程为
(1-16)
写成矩阵形式为
(1-17)
记作σ=Dε。其中,D为弹性矩阵。
平面应变问题中的物理方程为
(1-18)
写成矩阵形式为
(1-19)
同样记作σ=Dε。其中,D为弹性矩阵。
平面问题的位移边界条件为
(1-20)
平面问题的应力边界条件为
(1-21)
1.1.4 弹性力学中的能量原理
对于弹性力学基本方程,只要给出边界条件,理论上完全可以解出空间问题15个未知量、平面问题8个未知量。这种问题在数学上称为微分方程的边值问题。通常有3种基本解法,即按应力求解、按位移求解和混合求解。按应力求解以应力分量为基本未知函数,先求应力分量,再求其他未知量,是超静定问题,需要补充变形协调条件。由于位移边界条件不能改用应力分量表达,按应力求解时,弹性力学问题只能包含应力边界条件。按位移求解则以位移分量为基本未知函数,此时应通过物理方程和几何方程将平衡微分方程改用位移分量表达。应力边界条件也可以用位移分量表达,按位移求解时,弹性力学问题可以包括位移边界条件和应力边界条件。混合法就是以一部分应力分量为基本未知量,再以一部分位移分量为基本未知量,既建立变形协调方程,又建立内力平衡方程,最后加以求解。不管用哪种方法,工程实际中提出的弹性力学问题,能求得解析解的极其有限,多数还要用数值方法求解。
弹性力学的变分解法属于能量法,是与微分方程边值问题完全等价的方法,将弹性力学问题归结为能量的极值问题。能量表达成位移分量的函数,位移本身又是坐标的函数,因此能量是函数的函数,称为泛函。变分法就是研究泛函的极值问题。
1.虚功原理
弹性体在外力作用下发生变形,外力对弹性体做功,若不考虑变形中的热量损失、弹性体的动能以及外界阻尼,则外力功将全部转换为储存于弹性体内的位能——应变能。把虚功原理应用于连续弹性体,则可表述为:弹性体在外力作用下处于平衡状态,外力在弹性体所能发生的任何一组虚位移上所做虚功的代数和等于弹性体所储存的虚应变能。
弹性体某位置处在外力作用下实际发生的位移分量u、v、w,既满足位移分量表达的平衡微分方程,又满足边界条件以及用位移分量表达的应力边界条件。假想这些位移分量发生了边界条件所允许的微小改变,即所谓虚位移或位移变分δu、δv、δw成为,,,则外力在虚位移上所做的虚功为
当弹性体发生虚位移后,虚应变能为应力在虚位移引起的虚应变上所做的虚功。
假定在发生虚位移的过程中,没有其他形式的能量损失,依据能量守恒定理,应变能的增加等于外力在虚位移上所做的功,即虚应变能等于外力虚功。这就是连续弹性体的虚功原理(或称虚位移原理)。
(1-22)
虚功原理(虚功方程)可具体表示为
(1-23)
也可以写成矩阵形式,即
(1-24)
其中,[δ*]为虚位移列阵,[Fb]为体力列阵,[FN]为面力列阵,[ε*]为虚应变列阵,[σ]为应力列阵。
2.最小势能原理
外力从位移状态退回到无位移的初始状态时所做的功,称为外力势能,记为E。弹性体在这个退回过程中,内部产生变形势能(应变能),记为U。
(1-25)
(1-26)
变形势能和外力势能的总和称为总势能,。从前面的虚功原理看到,位移状态d为真实位移状态的充分必要条件是:对应位移d的总势能一阶变分为零,即对应位移d的总势能取驻值。满足位移边界条件的所有位移中,实际发生的位移使弹性体的势能最小,这就是最小势能原理。
(1-27)
对比虚功原理与最小势能原理,可知二者是完全等价的,一个用功的形式表达,另一个以能的形式表达。通过运算,可以由它们导出平衡微分方程和应力边界条件。
1.1.5 弹性力学有限元法
弹性力学中的三大类变量为位移、应变和应力,三大类方程是平衡方程、几何方程和物理方程。一般求解弹性力学问题的方法包括以位移为基本未知量的位移法,以应力为基本未知量,通过假设应力函数进行求解的逆解法和半逆解法。一般来说,能够直接进行解析求解的弹性力学问题相当少,而绝大多数弹性力学问题的求解需要借助于数值解法。有限元法为偏微分方程(组)提供了有效的数值近似求解手段,是数值求解弹性力学问题的重要途径。目前,大多数商业有限元软件在对弹性力学问题的分析中采用了按位移求解的方法。
有限元分析的基本步骤主要包括前处理、求解和后处理三个阶段。前处理阶段是将问题的求解域离散成有限个结点和单元,以结点的某些物理量作为基本未知量(在弹性力学问题中,一般是结点的位移),对单元进行分析,构造描述单元物理属性的形函数,描述每个单元的解答并建立起单元刚度矩阵,组装单元,形成总体刚度矩阵,并施加载荷、边界条件和初值条件;求解阶段一般是求解大型稀疏线性方程组,得到各个结点的位移值;后处理阶段是在得到结点的位移值后,进一步计算应力、应变、主应力、相当应力等,例如考虑屈服的强度问题中所需的von Mises应力等。下面通过一个简单的实例向读者逐步说明有限元求解的各个具体步骤。
假设有一个如图1-7所示的变截面直杆,其一端承受10kN集中力,上端的横截面积为100mm2,下端的横截面积为50 mm2,杆的长度为1m,弹性模量为E=200GPa,利用有限元法分析沿杆长度方向上不同点变形的大小及其轴向应力(忽略杆件的自重)。
图1-7 变截面直杆
1.前处理阶段
(1)单元离散。离散化为有限个结点和单元,简单起见,我们将变截面杆离散化为4个单元、5个结点,如图1-8(a)所示(说明:理论上离散的结点和单元数越多,结果可能越精确,但带来的计算量也随之增大)。
(a) (b) (c)
图1-8 将杆离散化成单元
(2)材料力学描述。假定一个近似描述单元特性的解。对于横截面积为A,长度为l,弹性模量为E的等截面直杆,在受到轴向力F的作用下且材料处于线弹性阶段时,由材料力学可知:
(1-28)
(1-29)
其中,σ是杆件横截面的正应力,ε是杆件横截面的正应变,是杆件的伸长量,FN是横截面上的轴向内力。
类比线性弹簧的控制方程F=kx,将方程式(1-28)和式(1-29)合并为
(1-30)
取等效刚度。由于直杆在y方向上横截面积是有变化的,作为近似,可以将该杆件看作一系列受到轴向载荷作用的具有不同横截面积的阶梯杆,如图1-8(b)所示,亦可以看作由4根具有不同等效刚度的弹簧串联而成的模型,如图1-8(c)所示。
设第i个结点的位移是u,结点1的约束力为FR,对各个结点进行受力分析,如图1-9所示。建立平衡方程如下:
图1-9 结点受力分析
(1-31)
其中,k1、k2、k3、k4分别表示4个单元的等效刚度。
将式(1-31)重组,并改写成矩阵形式,分离出载荷和约束力,则有
(1-32)
为了不同时考虑未知的约束力和位移,我们将已知结点1的零位移,即用u1=0取代方程式(1-32)的第一行,将式(1-32)改写为
(1-33)
式(1-33)也可表示为
KU=F (1-34)
其中,K是系统的整体刚度矩阵,U是结点位移矩阵,F是载荷矩阵。式(1-33)中的5个方程包含5个未知量,可以求出全部的结点位移,然后可以利用式(1-32)求出约束力。
(3)单元刚度矩阵。前述方法直接求出了系统的整体刚度矩阵,若分析其中的一个单元,则会发现这种单元只有轴向力的作用,对任一结点只有一个轴向位移。考虑弹簧单元内力和位移的关系,根据图1-10所示第i个结点和第i+1个结点的单元内力fi和fi+1,有
(1-35)
图1-10 弹簧单元
列向量称为单元的结点力向量,矩阵称为单元刚度矩阵,列向量称为单元的结点位移列阵。这三者亦有如下关系:
(1-36)
(4)单元组装。将式(1-35)所描述的单元刚度方程应用到所有单元,并将它们组合成整体刚度矩阵K。以单元②为例,该单元连接结点2和结点3,因此。用k(2G)表示单元②进行扩展后的矩阵,用以进行单元刚度矩阵的组装,也能够清晰地看出单元②在整体刚度矩阵中的位置:
通过在单元刚度矩阵的右边列上结点位移,可以帮助观察该结点对临近单元的影响。同样可以写出单元③扩展后的矩阵:
将各个单元在整体刚度矩阵中的位置进行组合,即对扩展后的矩阵相加,得到最后的整体刚度矩阵K:
(1-37)
其中,n表示单元的总数,这里n=4。
因此,整体刚度矩阵可显式地写为
(1-38)
此式与式(1-32)的刚度矩阵一样。有了整体刚度矩阵以后,我们就可以施加边界条件,并进一步根据式(1-34)进行求解。
2.求解阶段:求解线性方程组
假定各单元的长度相等,并以单元的平均横截面积作为单元的计算截面面积,计算等效刚度。各单元属性见表1-1。
表1-1 单元属性
根据式(1-38)算出整体刚度矩阵,有
引入边界条件u1=0,F=10×103N,根据式(1-33),有
划去第一行和第一列,只需要求解4×4的矩阵方程:
解得:u2=0.133 333 mm,u3=0.287 179 mm,u4=0.468 998 mm,u5=0.691 220 mm。
3.后处理阶段:求出应力和约束力
每个单元的平均应力可按式(1-39)计算:
(1-39)
代入位移计算结果后,求得
无论在何处截断杆件,截面的内力都是FN=10 000 N,因此对于每个单元,可以利用已知材料力学公式进行验证。结果表明,通过位移信息计算得到的单元应力和采用材料力学公式计算得到的单元应力完全相同,因此就本问题而言,位移计算是正确的。同样也可将结果代入式(1-32),求出约束力FR=10 000 N=10kN,显然该结果与整体平衡条件相符。
说明:这里计算的应力是单元的平均应力,并不是指结点的应力,在有限元法中,经常把与结点相关联的单元的应力进行平均作为结点的应力值。
1.1.6 平面问题的单元构造
杆、梁单元一般可以根据其构件之间的连接进行“自然离散化”,而对于连续变形体,需要在对象的几何域上按照分析的需要,采用“人工布置”结点和划分单元的方式进行有限元建模。这种离散方法称为“逼近性离散”,如图1-11所示。
图1-11 逼近性离散
在单元的划分过程中,要求单元之间仅在结点处连接,外载荷将被等效作用到结点上,同时假设单元内部的几何和物理特性是均匀的,这样就把原先连续体的无限自由度问题转变成了近似的有限多个自由度问题。对于弹性力学平面问题而言,用一个全局的位移函数表示整个结构的变形,对于复杂结构几乎是不可能的。但是对于每个单元,采用其结点位移进行插值来构造比较简单的函数,近似表达单元的真实位移,这是简单易行的。当结构被划分成足够多的单元后,把各单元的位移连接起来,就可以近似地表示整个区域的真实位移。
在平面问题的有限元处理中,最常见的单元是3结点三角形平面单元,如图1-12所示。这种单元共有3个结点、6个自由度,对应的有6个结点力。单元的结点位移列阵可以表示为uT=[uix uiy ujx ujy ukx uky]T,单元的结点力列阵可以表示为f T=[fix fiy fjx fjy fkx fky]T。若单元承受分布载荷,可将其等效到结点上。
图1-12 3结点三角形平面单元
1.单元的位移模式
从单元的结点位移可以看出,x、y方向上的位移场u(x,y)、v(x,y)可分别由3个结点的x和y方向位移确定,因此分别假设单元中各个方向的位移模式为
(1-40)
其中,a1~a6是待定系数。单元的结点条件为
(1-41)
当已知各结点的位移时,上面的6个方程可以求解6个未知量,计算出a1~a6,并将其回代入式(1-40),经整理,可写成形函数和结点位移向量的乘积形式:
(1-42)
其中,形函数为
(1-43)
系数ai、bi、ci分别按式(1-44)计算:
(1-44)
假设A是单元的面积,有
(1-45)
为了保证面积为正,我们规定结点i、j、k的次序按逆时针排列,如图1-12所示。将式(1-42)写成矩阵形式,有
(1-46)
其中,N是形函数矩阵,I是二阶单位矩阵,ue是单元结点位移向量。对于形函数,同样具有单位分解的特性,即单元中任意点的形函数之和为1;形函数在其对应的结点上值为1,在其他点为0。例如,在结点i上,Ni=1;在结点j、k上,Ni=0。对于Nj和Nk也有同样的性质,这是由插值函数的基本性质所决定的。
在有限元法中,位移模式决定计算误差。载荷的移置以及应力矩阵、刚度矩阵的建立都依赖于位移模式。为了正确反映弹性体中的真实位移形态,位移模式需要满足完备性条件和连续性条件。完备性条件要求位移模式能够反映单元的刚体位移。每个单元的位移不仅包含本单元变形引起的位移,还包括研究对象移动、转动引起的和变形无关的刚体位移,因此位移模式中必须要反映单元的刚体位移。另外,完备性条件还要求位移模式能反映单元的常量应变,每个单元的应变包括与该单元中各点坐标位置无关的应变,即所有点都相同的常量应变。随着单元尺寸的变小,各个点的应变趋于相等,也就意味着单元的变形趋于均匀,常量应变就成为应变的主要部分,这就能保证单元划分逐步增加的情况下,结果趋于真实解。连续性条件则是要求相邻单元在它们的公共结点、公共边上具有相同的位移。
2.应力转换矩阵和单元刚度矩阵
有了单元位移模式,根据几何方程可求出应变:
(1-47)
把形函数表示的位移模式代入式(1-47),则有
(1-48)
其中,矩阵L是算子矩阵,而矩阵B是应变位移矩阵。式(1-43)表明形函数是关于x、y的一次函数,在求一次偏导以后,其结果就为常数,容易算出:
(1-49)
由式(1-44)可知,对于位置确定的单元,其面积A以及各个系数均只和结点的坐标位置有关,所以应变矩阵中的各个元素都是常量。可见,应变ɛ的各分量也是常量。所以,3结点三角形平面单元内部各点的应变均相等,它是常应变单元。
根据物理方程,考虑平面应力问题,有
(1-50)
其中,μ是材料的泊松比。
对于平面应变问题,只需要将弹性模量E替换为,将泊松比μ替换为即可。根据式(1-48),应力矩阵可进一步写为
(1-51)
其中,S矩阵即为应力位移矩阵。对于某种弹性模量和泊松比确定的材料,D是常量矩阵,那么S矩阵也是常量矩阵,所以在每一个单元中应力分量也是常量。相邻单元一般不具有相同的应力,在它们的公共边上应力具有突变。但是,随着单元的逐步减小,这种突变将急剧变小,并不妨碍有限元法的解答收敛于正确解答。
通过虚功原理,可以推出3结点三角形平面单元的单元刚度矩阵为
(1-52)
其中,t是单元的厚度,其他符号含义同前。式中的是2×2的矩阵,可通过式(1-53)计算:
(1-53)
对于平面应变问题,只需要将弹性模量E替换为,将泊松比μ替换为即可。
同样也有单元刚度方程,将其展开,有
(1-54)
单元刚度矩阵具有如下性质。
(1)刚性,单元刚度矩阵描述了单元的刚度特性,即描述了单元受到外力作用时的应变和应力的关系。单元刚度矩阵越大,能够承受越大的力和扭矩。
(2)对称性,即。
(3)奇异性,即。单元刚度矩阵的奇异性表明,给定了单元结点载荷列阵,不能得出单元位移列阵,因为即使它满足平衡条件,单元还可以有任意的刚体位移。
(4)主元素恒为正,即单元刚度矩阵的对角线上的元素恒为正。
其他类型的单元刚度矩阵也具有上述性质,对于线性位移模式的三角形单元,可以证明两个相似三角形单元具有相同的单元刚度矩阵;单元水平或竖向移动不会改变刚度矩阵的数值。这种性质使得将求解域划分成相似三角形的单元,只需要计算一次单元刚度矩阵即可。
3.等效结点载荷
对于作用在单元内部的体力和边界上的分布面力,需要将它们移置到结点上成为结点载荷。这种移置必须按照静力等效的原则来进行。所谓静力等效,是指原载荷与结点载荷在任何虚位移上的虚功相等。在一定的位移模式下,这种移置的结果是唯一的。对于线性位移模式的3结点三角形平面单元,静力等效意味着原载荷(体力、面力)与结点载荷向任意一点简化得到的主矢和主矩都相等。
设单元在坐标为(x,y)的任意一点处受到集中力P的作用,如图1-13所示,可以将P表示为,等效结点载荷用表示,根据静力等效,可导出
(1-55)
图1-13 三角形单元载荷移置
其中,N是形函数矩阵。将其展开,有
(1-56)
若单元上受到分布体力的作用,则单元结点载荷列阵:
(1-57)
其中,是单元所覆盖的求解域。
若上述单元的ij边上有分布面力作用,则
(1-58)
4.结构的整体分析和支配方程
有限元网格中任取一个结点i,如图1-14所示,该结点受到环绕该结点的单元对它的作用力(内力)fi,这些作用力与各单元的结点力等值反向。另外,该结点上还有从环绕该结点的那些单元上移置过来的等效结点荷载(外力)Ri。根据平衡关系,有
图1-14 典型结点
(1-59)
对于所有结点均可以建立这样的方程,若有n个结点,平面问题就有2n个这样的方程。代入结点力公式,可得关于结点位移u的2n个线性代数方程组,即整体的有限元支配方程:
(1-60)
其中,K是整体刚度矩阵,其拼装组合方式与杆、梁单元类似,u为整体结点位移列阵,F是整体载荷列阵。整体刚度矩阵有对称性、奇异性、稀疏性的特点。如果结点、单元合理编号,整体刚度矩阵就具有非零元素带状分布的特点,即非零元素集中在以主对角线为中心的一条带状区域。每行的第一个非零元素到主元素之间的元素个数称为半带宽。半带宽越小,计算机求解的效率越高。
5.三角形单元分析实例
悬臂深梁如图1-15(a)所示,已知右端面作用有均布拉力,其合力为P。若按照图1-15(b)的方式离散化为2个单元、4个结点,设泊松比μ=1/3,厚度为t,求结点位移。
图1-15 3结点三角形平面单元实例
考虑单元①,与结点1、2、3 关联,取i、j、k分别为1、2、3。求出各系数,。
单元面积A=1m2。
本实例属于平面应力问题,根据式(1-52)、式(1-53),算出单元刚度矩阵中的各个元素,可得
同理,单元②与结点 1、3、4 关联,取i、j、k分别为1、3、4。算出单元刚度矩阵,即
拼装为整体刚度矩阵,即
则载荷阵列为
其中,F1x、F1y、F4x、F4y实际上是作用在结点1和结点4上的约束力分量。
根据整体有限元支配方程Ku=F,有
考虑边界条件,划去支配方程的第1、2、7、8 行和第1、2、7、8列,处理后得到
解此矩阵方程,可得
6.4结点矩形单元和6结点三角形单元
3结点三角形平面单元是有限元法中最早提出的单元,其适应边界能力强,但由于其是常应力单元,因此单元内的应变和应力都是常量,精度相对较低。
图1-16所示的矩形单元也是最基本的单元,它有4个结点、8个自由度。考虑单元的结点位移,x方向的位移场可以由该方向上的4个结点位移来确定,而y方向的位移场可以由该方向上的4个结点位移来确定。其位移模式为
图1-16 4结点矩形单元
(1-61)
在单元边界上,当x或y是一个常量时,对于单元的每一条边,位移是线性变化的,所以也称为双线性单元。在单元内部,由于存在xy的乘积项,位移是非线性变化的。根据平面问题的几何方程,应变是位移的一阶导数,如;应力和应变是线性关系,所以4结点矩形单元的应变模式和应力模式是一次线性变化的,其单元内部的应力不再是常量。虽然相邻的矩形单元在公共边界处的应力也有差异,但这种差异较小。在整理应力结果时,采用绕结点平均法,即将环绕某一结点的各单元在该结点处的应力求平均,用来代表该结点处的应力,这种方法的表征性较好。但是,矩形单元也存在比较明显的缺陷,它不能适应曲线边界或斜线边界,也不便在不同部位采用不同大小的单元。为了弥补这种缺陷,我们可以混合使用矩形单元和三角形单元。
在三角形单元的3条边上各增设一个结点,这样每个单元就有6个结点、12个自由度,可以采用二次完全多项式的位移模式。由于应变是位移的一阶导数,而应力与应变是线性关系,因此单元中的应力按照线性变化,能够更好地反映弹性体中应力的变化。在结点数目大致相同的情况下,其精度远高于3结点三角形单元。但是6结点三角形单元对于非均匀性及曲线边界的适应性却不如简单三角形单元,而且由于一个结点的平衡方程涉及较多的结点位移,整体刚度矩阵的带宽较大。
限于篇幅,这里不对4结点矩形单元和6结点三角形单元的单元构造做具体的推导和描述,有兴趣的读者可参阅相关教材和参考书。
7.轴对称单元
在柱坐标下,轴对称问题中的一些非对称力学变量都为0,其三大类力学变量如下所示。
位移:
应变:
应力:
在轴对称问题中,以上这些变量都只和坐标r、z有关,与θ无关。对于轴对称问题的有限元离散,在每一个截面,它的单元形状与一般的平面问题相同,但需要注意的是,轴对称问题的单元都是环形单元,如图1-17所示。
图1-17 轴对称单元
8.等参数单元
三角形单元能够应用于曲折的几何边界,但其精度较低;矩形单元虽然精度较高,但其适应性较差,不便使用在曲线边界和非正交的直线边界。因此,需要采用一些具有斜边的四边形单元。任意四边形单元可以通过已有的4结点矩形单元进行坐标映射的方式来获得,我们将这种通过坐标映射的方式构造的单元称为参数单元。
在全局坐标系Oxy下的4结点四边形单元的坐标称为物理坐标;其映射母单元是4结点矩形单元,称为基准单元;采用坐标系,称为基准坐标,如图1-18所示。
图1-18 单元映射
设两个坐标系的坐标映射为
(1-62)
基准坐标中的一点对应于物理坐标中的一个相应点,4个角点的结点映射条件为
(1-63)
每一个方向有4个结点条件,可利用包含4个待定系数的多项式来建立映射关系,即
(1-64)
将式(1-63)代入式(1-64),可以求出待定系数。将各系数回代入式(1-64),并参照位移模式以形函数和结点位移乘积的模式进行改写,有
(1-65)
其中,称为几何形状插值函数,具体为
(1-66)
将位移模式表达为基于结点位移的显式表达,即
(1-67)
其中,位移插值函数N(x,y)也可根据式(1-64)改写为。如果几何形状插值函数与位移插值函数N的阶次相同,则这种单元称为等参元;若的阶次小于N,则称为亚参元。等参元和亚参元能保证单元的收敛。若的阶次大于N,则称为超参元,超参元不能保证单元的收敛。
限于篇幅,这里仅介绍参数单元的基本思路,构造参数单元的单元刚度矩阵、载荷列阵需要用到坐标系之间的偏导数变换,还要利用坐标系之间的雅可比矩阵,同时在计算的过程中需要用到勒让德-高斯求积。具体推导请读者参阅相关教材和参考书。
1.1.7 动力学分析简介
动力学分析是用来确定惯性和阻尼起作用时结构或构件动力学特性的响应技术,其特性一般包括如下几个方面:①振动特性,包括结构的振动方式与振动频率;②载荷随周期性变化的响应效应,主要是指施加周期性变化的载荷时结构的位移和应力的响应情况;③周期振动或者随机载荷的效应,主要是指结构受周期性载荷或者随机载荷时的变化规律。
结构动力学分析有如下类型。
1.模态分析
在设计工程结构或机器部件的振动特性时,设计人员需要进行模态分析,即确定承受动态载荷结构设计中的重要参数(固有频率和振型),同时也可以此作为瞬态动力学分析、谐响应分析等其他动力学分析的基础。
模态分析最终目标是识别系统的模态参数,为结构系统的振动特性分析、振动故障诊断和预报以及结构动力学特性的优化设计提供依据。
模态分析通常求解如下方程:
(1-68)
其中,K为刚度矩阵,M为质量矩阵,u为固有模态位移矢量,是角频率。
模态问题求解的方法有雅克比法、Givens法、Householder法、对分法、逆迭代法、QR法、子控件法、兰索斯法等。
2.谐响应分析
谐响应分析是用于确定线性结构在承受随时间按正弦(简谐)规律变化的载荷时稳态响应的一种技术,旨在计算结构在几种频率下的响应,并得到一些响应值对频率的曲线。该技术只计算结构的稳态受迫振动,不考虑结构激励开始时的瞬态振动。通过谐响应分析,设计人员能预测结构的持续动力特性,从而验证其设计是否能够克服疲劳、共振及其他受迫振动引起的有害因素。
3.瞬态动力学分析
瞬态动力学分析(也称时间历程分析)是用于确定承受任意的随时间变化载荷的动力学响应的一种求解问题的方法,可用于分析确定结构在稳态载荷、瞬态载荷和简谐载荷的任意组合作用下随时间变化的位移、应变、应力及力。
瞬态动力学分析通常求解如下方程:
(1-69)
其中,K为刚度矩阵,M为质量矩阵,u为位移矢量,C为阻尼系数矩阵(通常表示为比例阻尼,即质量矩阵和刚度矩阵的比例,)。
瞬态动力学问题的求解方法有中心差分法、纽马克法等。
1.1.8 有限元分析中的若干问题
有限元模型是有限元程序可以处理的对象,是对实际结构的合理模拟。有限元模型一方面要保证力学的完整性,另一方面要保证计算的有效性。也就是说,建立的有限元模型既要承载完整的力学信息,尽可能真实地反映实际情况,又要保证计算机可以快速计算。
在开始建立有限元模型之前,设计人员需要对要解决的问题有透彻的认识,理解问题的力学本质,弄清楚结构几何特征、所受载荷性质、结构材料特性,初步估计响应情况。此外,设计人员需要根据力学概念,分析判断研究对象属于哪一类性质的问题,是线性问题还是非线性问题,是静力问题还是动力问题。
线性问题是指受力与变形成正比关系,例如,做出连续性、均匀性、线弹性、各向同性、小变形等基本假设,最后得出的结构刚度为常量,这就属于线性问题。非线性问题是指受力与变形不成正比关系,引起非线性的因素主要有3种。一是材料本身的刚度随变形而变化,如一般低碳钢受力比较小时,受力与变形呈线弹性关系,受力较大发生屈服时刚度很小,进入强化阶段具有一定刚度但小于线弹性阶段,这类问题称为材料非线性问题。二是结构在载荷作用下发生过大变形,结构形状影响刚度,这类问题称为几何非线性问题。当物体变形的大小与物体某个几何尺寸可以相比拟时,应按大变形来处理;当应变量大于0.3时,应按大应变问题处理。大变形、大应变问题都属于几何非线性问题。三是状态非线性,如结构中两零件为接触关系,随着受力变形,接触位置、接触面积都可能发生变化,接触刚度当然也会发生变化,即状态改变影响刚度,这类问题称为状态非线性问题。
静力问题是指所有变量和关系式都与时间无关,只要任一变量或关系式与时间有关,即属于动力问题。作为动力问题,若要计算结构的固有特性,则属于模态分析;若要计算在随时间变化的载荷的作用下结构各结点随时间变化的位移、速度、加速度及应力等,则属于瞬态响应分析;若要计算在随时间按照正弦或余弦规律变化的载荷的作用下结构的稳态响应,则属于谐响应分析,分析目的在于得到结构在不同频率简谐载荷的作用下响应与频率的关系;若要计算结构在某载荷谱作用下的位移、应力等,则属于谱响应分析,载荷谱可以是位移、速度、加速度或力随频率变化的关系图,谱分析主要用于确定结构对随机载荷或随时间变化载荷(如地震、风载荷等)的动力响应情况,可代替瞬态响应分析。
建模前,设计人员还要根据物体的形状和受力情况,判断是一维问题、二维问题还是三维问题,有无对称性、反对称性或周期对称等可利用。这是一个把工程问题进行力学建模的过程。有了准确的力学模型,有限元建模就有了基础。有限元建模过程包括选择单元类型、确定单元的尺寸大小、保证网格划分质量、定义材料和单元特性、处理载荷和边界条件、确定计算方法和控制参数、求解、输出结果等。如果计算结果揭示了事物的内在规律,说明有限元模型和实际物理模型的力学性能是一致的,那么所建的模型就是个好模型。如果计算结果不符合实际情况或者计算进行不下去,那么可能出现单元类型不对、网格数量太少、材料模型错误、约束和载荷的施加方式不对、接触定义有问题、网格质量差、计算方法不对等问题,建模过程中的每个因素都可能造成计算结果错误或计算困难。针对工程问题,利用有限元分析方法,要想得到一个实用的计算结果,准确建模是关键。
1.有限元建模的准则
有限元建模的准则是根据工程分析的精度要求,建立合适的、能模拟实际结构的有限元模型。要使分析结果有足够的精度,所建立的有限元模型必须在能量上与原连续系统等价。有限元模型具体应满足下述准则。
(1)满足平衡条件。结构的整体和任意一单元在结点上都必须保持静力平衡。
(2)满足变形协调条件。交汇于一个结点上的各单元在受力变形后也必须保持交汇于同一结点。
(3)满足边界条件和材料的本构关系。边界条件包括整个结构的边界条件和单元间的边界条件。
(4)符合刚度等价原则。有限元模型的抗拉压、抗弯曲、抗扭转、抗剪切刚度应尽可能与原来结构等价。
(5)认真选取单元,包括单元类型、形状、阶次,使其能够很好地模拟几何形状、反映受力和变形情况。单元类型有杆单元、梁单元、平面单元、板单元、空间单元等,空间块体又分四面体块单元或六面体块单元,六面体块单元又分八结点六面体或二十结点六面体等。选取单元时应综合考虑结构的类型、形状特征、应力和变形特点、精度要求和硬件条件等因素。
(6)应根据结构特点、应力分布情况、单元的性质、精度要求及其计算量的大小等仔细划分计算网格。
(7)在几何上要尽可能地逼近真实的结构体,其中特别要注意曲线与曲面的逼近问题。
(8)仔细处理载荷模型,正确生成结点力,同时载荷的简化不应该跨越主要的受力构件。
(9)质量的堆积应该满足质心及惯性矩等效要求。
(10)超单元的划分尽可能单级化并使剩余结构最小。
2.边界条件的处理
对于基于位移模式的有限元法,在结构的边界上必须严格满足已知的位移约束条件。例如,弹性体某位置处有固定支撑,这些边界上的位移、转角等于零,如图1-19(a)所示,,图1-19(b)中,;或者弹性体某位置处位移或转角有已知值,如图1-19(c)中,,计算模型必须让它能实现这一点。
图1-19 各种边界条件示例
有时边界支撑不是沿坐标方向,称为斜支撑,如图1-19(d)所示,可以设定与整体坐标不一致的结点坐标方向来实现。还有的约束是单向的,例如,绳索只能承受拉力,光滑支撑面只能提供压力,这就需要按非线性对待。
当边界与另一弹性体相连,构成弹性边界时,我们可分两种情况处理。如果弹性体对边界点的支撑刚度已知,那么可将它的作用简化成弹簧,在此结点上加一弹簧单元,如图1-19(e)所示。如果弹性体对边界点的支撑刚度不清楚,那么可将此弹性体的一部分划出来和结构连在一起进行分析,所划区域的大小视其有影响的区域大小而定,如图1-19(f)所示。
如果整个结构存在刚体位移,就无法进行静力分析、动力分析。为此,我们必须根据实际结构的边界位移约束情况,对模型的某些结点施加约束,消除结构的刚体位移影响。对于平面问题,应消去结构的两个平移自由度、一个转动自由度;对于三维问题,须消去3个平移自由度、3个转动自由度。此外,要保证这些消除模型刚体位移的约束施加得当,如果不恰当,就会产生不真实的支反力,改变了原结构的受力状态和边界条件,从而得到错误的结果。例如,在图1-19(g)中根据对称性,C点两方向位移均为零,因此对C点施加约束是适当的。若把点A、B、D的两方向位移指定为零,则与实际情况不符。
3.连接条件的处理
复杂结构通常由杆、梁、板、壳及二维体、三维体等多种形式的构件组成。由于构件中各组件之间的自由度个数不匹配,因此在梁和二维体、板、壳和三维体的连接处必须妥善加以处理,否则模型会失真,得不到正确的计算结果。
在复杂结构中,我们还能遇到各种各样其他的连接关系,只要将这些连接关系彻底弄清楚,就能写出相应的位移约束关系式,这些关系式为构件间复杂的连接条件,需要在计算中使程序严格满足这些条件。
应当指出,在不少的实用结构有限元分析程序中,已为用户提供输入连接条件的接口,用户只需严格遵守用户使用规定,程序将自动处理自由度之间的用户所规定的位移约束条件。
1.1.9 减小解题规模的常用措施
有限元的计算时间和结点数的多少有很大关系。在保证计算精度的条件下,用户应采用各种手段减少结点数,以节约计算时间和成本。
1.对称性和反对称性
如果计算对象的结构具有对称性,就可利用这个特点减少参加计算的结点数。所谓结构的对称性,是指结构的几何形状和支撑条件对某轴(面)对称,同时截面和材料性质也对称于此轴(面)。也就是说,结构绕对称轴对折后,左、右两部分完全重合。
如果对称结构上有对称载荷作用,则变形和应力也是对称的。只需取一半的结构建模即可,对称轴上的结点给出对称边界条件,算完后还可以根据对称性扩展出另一半结果。这样解题规模可减小一半。
如果作用在对称结构上的载荷是反对称的,即将结构绕对称轴对折后,两载荷的作用点重合,载荷大小相同,但载荷方向相反。根据结构力学可知,在反对称载荷作用下,结构的位移及应力都将反对称于对称轴。
在常用的商业有限元软件中,只要用户给出对称或反对称条件,程序会自动加上相应的位移约束。
2.周期性条件
有些结构可以划分为若干形状完全相同的子结构,当任意一子结构绕对称轴旋转一定角度时,该子结构的形状将与其他子结构完全重合,这种结构称为循环对称结构或周期对称结构。工程中常见的风扇叶片、花键、螺旋桨、齿轮、法兰等都是周期对称结构。如果结构所受载荷和位移的约束也是周期对称的,且各子结构材料和物理特性也完全相同,则应力和变形关于同一轴周期对称。若所受载荷不是周期对称的,如齿轮、法兰,则不属于周期对称问题。对于周期对称问题,计算时可以只取一个子结构进行分析。注意:在取一个子结构时,应使应力集中区域在子结构内部而不在边界。
另外,还有一种周期对称结构可以看作由一个子结构沿某一方向多次重复得到,称为重复对称结构。如果结构所受载荷和约束同样满足重复对称条件,与循环对称类似,那么只需要模拟和分析一个子结构即可。
3.降维处理和几何简化
对于一个复杂的工程构件,设计人员可以根据其在几何学、力学或传热学上的特点,进行降维处理。对于一个三维物体,如果可以忽略某些几何上的细节或次要因素,就能按照二维问题来处理。例如螺纹连接结构中,由于螺纹升角很小,可认为螺纹牙的受力在周向是相同的,从而近似看作轴对称结构。一个二维问题,若能近似地看作一维问题,就尽量当作一维问题计算。维数降低,计算量将大大降低。例如,齿轮、连杆、径向轴承等许多零件的结构计算都可以近似作为平面问题。在复杂的结构计算中,应尽量减少按三维问题处理的部分。
此外,某些零件上会有许多小圆孔、小圆角、小凸台、浅沟槽等几何细节,细节的存在将影响网格的大小、数量及分布。因为在自动划分网格时,一段直线或曲线至少划分一个单元边,一个平面或曲面至少划分一个单元面,一个圆最少也要3个单元边来离散,所以细节将限制网格的大小。另外,单元由密到疏应该平缓过渡,这也会影响整个模型的网格数量和分布。但细节的取舍要遵循两条原则:一是细节处应力的大小,只要这些不是位于应力峰值区域中分析的要害部位,根据圣维南原理就可以将其忽略;二是与分析的内容也有关系,一般情况下,由于细节会影响应力的大小及分布,静应力、动应力计算中要注意细节的影响,而结构的固有频率和模态振型主要取决于质量分布和刚度,因此计算固有特性时就可以少考虑细节。
现代机械设计中进行力学计算的目的,往往在于求出结构最大承载能力或结构最薄弱区域,这些处理方法虽然会带来一定的误差,但一般都能满足工程上的设计要求,而计算成本却能大大降低。如果对个别部位分析后不能满意,则可将这块单独取出再作细致分析。
4.子结构技术
当计算的结构比较复杂时,整体刚度矩阵的阶数往往会很大,从而超出计算机容量,这时可以考虑一小块一小块地来计算,最后再将各子块边界结点归结在一起,这就是子结构分析法。
子结构方法还可以用在需要局部精确分析的场合,如应力集中位置、局部发生塑性变形需要进行非线性分析的地方、设计可能改变的局部等,设计人员可以只重复计算部分结构,节约计算时间和计算成本。
子结构分析法的基本思路:①几何分割;②子结构离散;③定义边界自由度;④凝聚内部自由度;⑤子结构集成;⑥求解整体模型;⑦回代。现有大型有限元程序一般包括子结构法的内容,用户根据需要调用即可。
5.线性近似化
在工程上,我们常常将一些呈微弱非线性的问题当作线性问题来处理,所得结果既能满足工程要求,又可降低成本。例如,许多混凝土结构(水坝、高层建筑、冷却塔、桥梁、大型机电设备地基等)实际上都是非线性结构,其非线性现象较弱,初步分析时,常看作线性结构。只有当分析其破坏形态时,才按非线性考虑。
6.多种载荷工况的合并处理
有时我们要对一个结构进行多种载荷工况的分析,如果每一种工况都作为一个新的问题重新分析一次,则计算量会很大,也没有必要。对于上述情况,用户可以将每一种载荷矢量合并成载荷矩阵R,一起进行求解。方程系数只需进行一次三角分解,计算量就会大大降低。对于线性问题,用户还可以先解出某些标准载荷模式、、下的解、、,若其他载荷模式可以写成这些载荷的线性组合,则它对应的解也是这些解的线性组合,即,其中a、b、c为线性组合系数。
7.结点编号的优化
有限元计算需要输入和存储大量信息,计算量也非常大,这就要求编程人员探索减少存储量、减少运算次数的方法。
有限元算法计算量大致与总体刚度带宽的平方成正比。为了减少存储量,我们可以根据总体刚度矩阵具有对称性、稀疏性、带状分布等特点,采用一种名为“半带存储”的技术,即只存储总体刚度矩阵中沿主对角线非零元素的一半的数据。对同一个模型,如果按不同的次序对各结点进行编号,那么得到的总体刚度矩阵形式是不同的,半带宽也不一样,相应存储量也就不同。