最新7.4-单元刚度矩阵组装及整体分析
最新7.4-单元刚度矩阵组装及整体分析
根据全结构的平衡方程可知,总体刚度矩阵是由单元刚度矩阵集合而成的个结构的计算模型分成个单元,那么总体刚度矩阵可由各个单元的刚度矩阵组装而成,即是由每个单元的刚度矩阵的每个系数按其脚标编号“对号入座”叠加而成的将总体坐标轴分别用表示,对某单元有式中,和分别是局部坐标系和总体坐标系下的单元结点位移向量是该单元在总体坐标系下的单元刚度矩阵.将单元结点的局部编号换成总体编号,其中右上角的上标表示第单元所累加上的子矩阵具有相同的下标,的那些子矩阵的累加总体刚度矩阵第行的非零子矩阵是由与结点相联系的那,从环绕点各单元移置而来的结点载荷为式中表示对环绕结点的所有单元求和,环绕结点的各单元施加于结点的结点力为.因此,结点的平衡方程可表示为得到以结点位移表示的结点的平衡方程,为整体刚度矩阵,为全部结点位移组成的向量,为全部结点载荷组成的向量式中,是总体坐标系下的结点载荷向量,为坐标转换阵.构是处于自由状态,在结点载荷的作用下,结构可以产生任意的刚体位移的条件下,仍不能通过平衡方程惟一地解出结点位移.约束的种类包括使某些自由度上位移为零,,或给定其位移值,还有给定支承刚为了理解这个方法,我们把方程分块如下:其中,假设是给定的结点位移;是无约束的(自由)结点位移因而是已知的结点力;其中,不是奇异的,因而可以解方程(一旦知道了,求得未知结点力.殊情况下,我们可以删除对应于的各行和各列(即删行删列法),故可把方程简写为由于全部给定的结点位移通常都不能在位移向量的开始或终了,故分块法的编号方法是很麻烦因此,为了引入给定的边界条件,可以采用下述等价的方法如果把给定为,则载荷向量为结点自由度总数中对应于的行和列为零,而对角线元素为)在载荷向量中引入规定的值,即对全部规定的结点位移均应反复运用上述过程(步骤(置大数法的思路是:在总体刚度矩阵中,把指定位移所对应的行和列的对角元素乘上一个很大的数,如,此行其他元素保持不变,同时把该行对应的载荷项也相应地用来代替,这里为指定位移,于是原平衡方程组变为除第行外,其他各行仍保持原来的平衡特性,而第个方程式展开为由于上式中的比其他项的系数大得多,求和后可略去其小量,则上式变为即.边上有,若结构的总体坐标系为为斜支座的局部坐标系(见图对于边界结点,须限定方向位移,为此,将边界结点的位移及载荷都变换到局部坐标轴系设轴与斜支座的轴夹角为,逆时针为正,其中,.)中第行左右两边前乘以上式的系数矩阵仍然是对称的,而且此方程中结点位沿轴表示,这样,限定方向的位移异性,解这个线性代数方程组可求出结位移.阶线性代数方程,需进行次消元行元素作为主元行,为主元,对第行元素()的消元公式为式中等的上角码(次消元后的系数矩阵和载荷阵分别记为及.式表时第我们把消元最后结果记为,为上当回代求解时,已经解得总体刚度平衡方程中,,是单位上三角矩阵,.记,则.由其中第一个方程解得,再由第二个方程解得,向上回代,可得,由得依此类推可求得.由平衡方程组解出位移后,从中分离出各单元的结点位移,再通过方程)等计算各单元的应变、应力和结点力等内力。
matlab组装整体刚度矩阵
matlab组装整体刚度矩阵摘要:一、引言- 介绍MATLAB 软件及其在工程领域的应用- 简述组装整体刚度矩阵的重要性二、MATLAB 软件介绍- 解释MATLAB 软件的基本功能- 说明MATLAB 软件在工程计算中的优势三、刚度矩阵的定义与作用- 定义刚度矩阵- 说明刚度矩阵在工程计算中的作用四、MATLAB 组装整体刚度矩阵实例- 给出一个具体的刚度矩阵组装实例- 详细解释实例中MATLAB 的使用方法五、总结- 总结整体刚度矩阵的组装过程- 强调MATLAB 软件在整体刚度矩阵组装中的重要性正文:一、引言MATLAB 是一种功能强大的数学软件,广泛应用于科学计算、数据分析、可视化等领域。
在工程领域中,MATLAB 软件同样具有广泛的应用,尤其是在进行有限元分析时,MATLAB 软件可以有效地帮助工程师完成复杂的计算任务。
组装整体刚度矩阵是有限元分析中的一个重要步骤,通过这一步骤,工程师可以更好地分析结构的力学性能。
二、MATLAB 软件介绍MATLAB 软件是由美国MathWorks 公司开发的一种数学软件,其全称为“Matrix Laboratory”。
MATLAB 软件具有强大的数值计算、数据分析、可视化等功能,可以方便地处理各种数学问题。
在工程领域中,MATLAB 软件被广泛应用于结构分析、控制系统、信号处理、图像处理等领域。
由于MATLAB 软件具有丰富的函数库和强大的计算能力,因此,工程师可以利用MATLAB 软件快速地完成各种复杂的计算任务。
三、刚度矩阵的定义与作用在有限元分析中,刚度矩阵是一个重要的概念。
刚度矩阵反映了结构在受力时的变形程度,可以帮助工程师更好地了解结构的力学性能。
刚度矩阵的定义如下:K = [Kx, Ky, Kz]其中,Kx、Ky 和Kz 分别表示x、y 和z 方向的刚度。
刚度矩阵在工程计算中的作用主要体现在以下几个方面:1.反映结构的刚度特性:刚度矩阵可以反映结构在受力时的变形程度,帮助工程师了解结构的刚度特性。
单元刚度矩阵
单元刚度矩阵单元刚度矩阵是在结构力学中的一个重要概念,它是一个矩阵,用来表示刚性和结构特性。
它可以用来描述广泛的结构,如桥梁,大型建筑物和其他复杂的构造。
它的研究有助于更好地理解结构的运动和反应。
它也可以用来预测和控制结构的变形和损坏,从而减少结构建设过程中可能发生的各种问题。
单元刚度矩阵是一个n x n等阶矩阵,其中n是一个复杂结构中的单元数量。
它代表了单元之间的约束关系,表明它们如何互相影响。
这也就是所谓的单元刚度矩阵。
每个矩阵元素代表了任意两个单元之间的受拉或受压力的数量,可以用来计算结构中每一个单元之间的刚度和约束。
单元刚度矩阵有几种不同的类型,其中一种是静态刚度矩阵,它用来表示复杂结构在静态荷载作用下的刚度。
它可以用来预测荷载作用下结构变形的情况,并作出相应的改善。
另外一种是有限元分析,它可以用来对复杂结构在动态荷载下的变形,受力,反应,以及可能发生的结构破坏作出分析。
单元刚度矩阵的计算方法有很多。
有些是利用有限元分析的方法来进行的,也有些是直接从节点和单元的计算和配置来得出的。
有些方法只需简单地求解结构中一组特定问题,而另一些方法则要求对结构中所有部件进行复杂的数值计算。
单元刚度矩阵的计算可以帮助从两个角度来改善设计:一方面,单元刚度矩阵可以帮助改善结构运动的性能,另一方面,它可以帮助减少结构上可能发生的变形以及提高结构的耐久性。
单元刚度矩阵的计算和研究非常重要,现代的结构力学和建筑设计工程正在用这个技术来设计新型的可靠性更高,耐久性更强的建筑结构。
基于单元刚度矩阵的计算和研究,科学家们可以更好地理解结构力学,并减少建筑物的再建设和变形,以及可能发生的损坏。
总之,单元刚度矩阵的研究和计算存在着很多的优势。
现代的结构力学和建筑设计都需要用到它,以便更好地分析和控制结构的变形和损坏。
它的研究也有助于开发更安全,更高效的建筑结构,有助于结构力学中的其他方面的研究。
单元刚度矩阵(整体坐标系)[详细]
§9-1 概述 §9-2 单元刚度矩阵(局部坐标系) §9-3 单元刚度矩阵(整体坐标系) §9-4 连续梁的整体刚度矩阵(先讲) §9-5 刚架的整体刚度矩阵 §9-6 结构整体结点荷载 §9-7 计算步骤和算例
▲ 竖向杆件坐标变换的简化技巧 §9-8 忽略轴向变形时刚架的整体分析 §9-9 桁架及组合结构的整体分析
31 2
32
1
x
0 1 0
1 0 0
0
T
0
01
0
1 0
0
1 0 0
0 0 1
②
②
5 6
4
(局部坐标)
4 6
5
y
(整体坐标)
整体坐标下的单元刚度矩阵:
k② T T k ② T
结点位移码
(
(1 2 3 0 0 0)
结点码
1
2
12 0 30 12 0 30 1
0
300
0
0
300
k ② T T k ② T
9.65
7.13
k
②
0.45 9.65
7.13 0.45
1 3
7.13 5.50 0.6 .137 5.50 0.6
0.45 0.6 5.0 0.45 0.6 2.5
9.65 7.13 0.45 9.65 7.13 0.45
7.13 5.50 0.6 7.13 5.50 0.6
0.0 0.69 2.08
0.0 0.69 2.08
0.0
2.08
4.17
0.0
2.08
8.33
①
①
k k
单元②: 15.0
【优】刚架的整体刚度矩阵(assembly最全PPT资料
2
2结)束各播杆放方请向点不“尽后相退同”,。要进行坐标变换;
点1 击左2键,一3 步步0播放。0 0
结1)点结位点移位列移阵分:{量Δ}增=[Δ加1到Δ三2 Δ个3;Δ4]T
1、结2点位移3分量0的统一0编码—0 —总码
结点力列阵:{F}=[F1 F2 F3 F4]T
2)各杆方向不尽相同,要进行坐标变换;
-30300 13020
300 00
+111000000 300 300 13200
030 00
55000 30 ××110044
点=[u击A左v键A,θ一A θ步C步]T播放。 13)除2了刚3结点,0还要考0 虑铰结4 点等其它情况。 1结束播2放请点3 “后0退”。0 0
00 31020 030 00 30102 0304
点击左键,一步步播放。结束播放请点“后退”。
重新播放请点 重新播放
0
2
3
1
A1
2
4
C
0 x(2)
(1)
结点位移列阵:{Δ}=[Δ1 Δ2 Δ3 Δ4]T
=[uA vA θA θC]T
(3)结点力列阵:{F}=[F1 F2
F3 F4]T (5)
(2)
(3)
(6) (4)
B
2
(1)
1
0
0
0y
2、单元定位向量 (4) (6)
(5)
1
{λ} = [1 2 3 0 0 4]T
2
{λ} =
[1 2 3 0 0 0]T
点击左键,一步步播放。结束播放请点“后退”。
3、单元集成过程
11 2 33 0 0 04
+13120200 00 -0330013200 00 030
abaqus 单元刚度矩阵
Abaqus 单元刚度矩阵——解析有限元分析中的基本工具引言(Introduction):有限元分析(Finite Element Analysis,FEA)是一种广泛应用于工程和科学领域的数值模拟方法。
它通过将复杂的连续体划分为简单的几何形状,并对每个几何单元进行数学建模,来近似求解实际问题。
在有限元分析中,单元刚度矩阵是一个重要的概念,它描述了单元的刚度特性,对于计算整体系统的行为非常有用。
本文将重点介绍Abaqus软件中的单元刚度矩阵。
首先,我们将简要回顾有限元分析的基本概念和步骤。
接着,我们将探讨单元刚度矩阵的定义和计算方法。
然后,我们将通过一个简单的示例案例来说明单元刚度矩阵的应用。
最后,我们将总结单元刚度矩阵在有限元分析中的重要性和应用前景。
有限元分析基础(Basics of Finite Element Analysis):有限元分析的基本步骤通常包括几何建模、网格剖分、物理特性分配、边界条件设置和结果解析等。
在进行数学建模时,连续体被分割成称为单元的小体积区域,每个单元内部的行为则通过数学公式进行描述。
这些单元通常是三角形、四边形、六面体等几何形状。
单元刚度矩阵的定义(Definition of Element Stiffness Matrix):单元刚度矩阵是描述单元在给定边界条件下的刚度特性的矩阵。
它由单元的几何属性、材料特性和积分算法决定。
在Abaqus软件中,单元刚度矩阵是通过数值积分方法计算得出的。
单元刚度矩阵计算方法(Calculation of Element Stiffness Matrix):单元刚度矩阵的计算涉及到单元的几何形状、材料特性、积分算法等因素。
不同类型的单元有不同的刚度计算方法,通常包括弹性理论和数值积分。
以Abaqus中的三角形单元为例,其刚度矩阵通常可以通过以下步骤计算:1.定义单元的几何属性,如节点坐标。
2.根据几何属性和材料特性,计算出单元的刚度矩阵表达式。
7.4 单元刚度矩阵组装及整体分析
7.4 单元刚度矩阵组装及整体分析7.4.1 单刚组装形成总刚根据全结构的平衡方程可知,总体刚度矩阵是由单元刚度矩阵集合而成的.如果一个结构的计算模型分成个单元,那么总体刚度矩阵可由各个单元的刚度矩阵组装而成,即[K]是由每个单元的刚度矩阵的每个系数按其脚标编号“对号入座”叠加而成的.这种叠加要求在同一总体坐标系下进行.如果各单元的刚度矩阵是在单元局部坐标下建立的,就必须要把它们转换到统一的结构(总体)坐标系.将总体坐标轴分别用表示,对某单元有式中,和分别是局部坐标系和总体坐标系下的单元结点位移向量;[T]为坐标转换阵,仅与两个坐标系的夹角有关,这样就有是该单元在总体坐标系下的单元刚度矩阵.以后如不特别强调,总体坐标系下的各种物理参数均不加顶上的横杠.下面就通过简单的例子来说明如何形成总体刚度矩阵.设有一个简单的平面结构,选取6个结点,划分为4个单元.单元及结点编号如图3-27所示.每个结点有两个自由度.总体刚度矩阵的组装过程可分为下面几步:图7-27(1)按单元局部编号顺序形成单元刚度矩阵.图7-27中所示的单元③,结点的局部编号顺序为.形成的单元刚度矩阵以子矩阵的形式给出是(2)将单元结点的局部编号换成总体编号,相应的把单元刚度矩阵中的子矩阵的下标也换成总体编号.对下图3-27所示单元③的刚度矩阵转换成总体编号后为(3)将转换后的单元刚度矩阵的各子矩阵,投放到总体刚度矩阵的对应位置上.单元③的各子矩阵投放后情况如下:(4)将所有的单元都执行上述的1,2,3步,便可得到总体刚度矩阵,如式(3-9).其中右上角的上标表示第单元所累加上的子矩阵.(3-9)(5)从式(3-9)可看出,总体刚度矩阵中的子矩阵AB是单元刚度矩阵的子矩阵转换成总体编号后具有相同的下标,的那些子矩阵的累加.总体刚度矩阵第行的非零子矩阵是由与结点相联系的那些单元的子矩阵向这行投放所构成的.7.4.2 结点平衡方程我们首先用结构力学方法建立结点平衡方程.连续介质用有限元法离散以后,取出其中任意一个结点,从环绕点各单元移置而来的结点载荷为式中表示对环绕结点的所有单元求和,环绕结点的各单元施加于结点的结点力为.因此,结点的平衡方程可表示为(3-10)以[K]代入平衡方程,得到以结点位移表示的结点的平衡方程,对于每个结点,都可列出平衡方程,于是得到整个结构的平衡方程组如下:式中,[K]为整体刚度矩阵,为全部结点位移组成的向量,为全部结点载荷组成的向量.当然,如果各点的载荷向量也是在单元局部坐标下建立的,在合成以前,也应把它们转换到统一的结构(总体)坐标系下,即式中,是总体坐标系下的结点载荷向量,为坐标转换阵.7.4.3 位移边界条件在有限元法对结构进行整体分析时,建立了整体刚度矩阵[K],也得到了结构的刚度平衡方程,即.结构刚度方程的求解相当于总刚[K]求逆的过程.但是,从数学上看,未经处理的总刚是对称、半正定的奇异矩阵,它的行列式值为零,不能立即求逆.从物理意义看,在进行整体分析时,结构是处于自由状态,在结点载荷的作用下,结构可以产生任意的刚体位移.所以,在已知结点载荷的条件下,仍不能通过平衡方程惟一地解出结点位移.为了使问题可解,必须对结构加以足够的位移约束,也就是应用位移边界条件.首先要通过施加适当的约束,消除结构的钢体位移,再根据问题要求设定其他已知位移.所以,处理位移边界条件在有限元分析步骤中十分重要.约束的种类包括使某些自由度上位移为零,,或给定其位移值,还有给定支承刚度等,本书涉及前两种.处理约束的方法,常用的有删行删列法、分块法、置大数法和置“1”法等,下面分别予以介绍.1、删行删列法若结构的某些结点位移值为零时(即与刚性支座连接点的位移),则可将总体刚度矩阵中相应的行列、删行删列划掉,然后将矩阵压缩即可求解.这种方法的优点是道理简单.如果删去的行列很多,则总体刚度矩阵的阶数可大大缩小.通常用人工计算时常采用该方法.若用计算机算题,在程序编制上必带来麻烦,因为刚度矩阵压缩以后,刚度矩阵中各元素的下标必全改变.因而一般计算机算题不太采用.2.分块法为了理解这个方法,我们把方程分块如下:(3-11)其中,假设是给定的结点位移;是无约束的(自由)结点位移.因而是已知的结点力;是未知的结点力.方程(3-11)可以写为即(3-12)和(3-13)其中,不是奇异的,因而可以解方程(3-12)得出(3-14)一旦知道了,就可以由方程(3-13)求得未知结点力.在全部给定的结点自由度都等于零的特殊情况下,我们可以删除对应于的各行和各列(即删行删列法),故可把方程简写为(3-15)3.置“1”法由于全部给定的结点位移通常都不能在位移向量的开始或终了,故分块法的编号方法是很麻烦的.因此,为了引入给定的边界条件,可以采用下述等价的方法.可以把方程(3-12)和(3-13)合在一起写为(3-16)在实际计算中,方程(3-16)所示的过程可以在不重新排列所述方程的情况下用下述分块的方法为进行.步骤(1)如果把给定为,则载荷向量P可以修改为为结点自由度总数.步骤(2)除对角线元素以外,使[K]中对应于的行和列为零,而对角线元素为1,即步骤(3)在载荷向量中引入规定的值,即对全部规定的结点位移均应反复运用上述过程(步骤(1)到(3)).应当指出,由于这个过程保持了方程的对称性,因此,[K]可以按带状存储,而且几乎不会增加编制程序的工作量.4.置大数法置大数法的思路是:在总体刚度矩阵中,把指定位移所对应的行和列的对角元素乘上一个很大的数,如,此行其他元素保持不变,同时把该行对应的载荷项也相应地用来代替,这里为指定位移,于是原平衡方程组变为除第行外,其他各行仍保持原来的平衡特性,而第个方程式展开为由于上式中的比其他项的系数大得多,求和后可略去其小量,则上式变为即.这样就用近似方程组代替原方程组,得到近似满足边界条件的解.当指定位移为零时,只要将对角元素乘上一个大数,而相应的载荷项经证明可以不置零.删行删列法适用于指定零位移点,而置大数法适用于给定位移(包括零位移).5.斜支座的处理对于简单的约束情况(如限定某些结点位移为零或取得给定数值),可以用前述置大数法处理.有的结构在直角坐标系内建立了位移方程组,但在某个斜边上受有法向约束.如图3-28所示正方形固支板,受均布横向载荷,对此,可利用对称性而只计算其1/8,如图中ABC部分,其中AC为固支边,按对称性,AB边上有,但在BC边上应限定绕BC的转用等于零.为处理此类斜边上的约束,须对斜边上的结点做坐标变换.若结构的总体坐标系为为斜支座的局部坐标系(见图3-29).对于边界结点,须限定方向位移,为此,将边界结点的位移及载荷都变换到局部坐标轴系.设轴与斜支座的轴夹角为,逆时针为正,图7-28 图7-29 则依据第二单中坐标转换关系有其中,.或写成(3-17)与位移关系相同有(3-18)将上两式带入结构刚度方程有(3-19)这样把位移到列阵中凡是斜支座的结点位移矢量都用局部坐标表示了.将式(3-19)中第行左右两边前乘以(3-20)由上式可见:凡是边界点的斜支座,在刚度方程中对应于斜支座的位移和载荷向量均可直接斜支座的局部坐标值,总刚度距阵中的相应行列需作相应的变换.上式的系数矩阵仍然是对称的,而且此方程中结点位沿轴表示,这样,限定方向的位移就很方便了.实际计算中,并不需要建立结构总的位移方程组后再进坐标变换.而可以在形成单元刚度矩阵和结点载荷之后,就对斜支座点进行坐标变换,把变换后的单元刚度矩阵和结点载荷叠加入总刚度矩阵和总载荷的相应位置,最后叠加形成的也就是方程组(3-20),即需要处理的结点,应该在单元计算中完成坐标变换后再叠加,当结构有不同的斜边约束时,都可以这样处理,只不过对不同边上的结点,应按不同的方向余弦矩阵变换就是了.7.4.4 总刚度平衡方程的求解应用有限元法,最终都是归结为解总体刚度平衡方程,它实际上是以总体刚度矩阵为系数矩阵的大型线性代数方程组.通过对结构施加位移边界条件,消除了结构的刚体位移,从而消除总体刚度矩阵的奇异性,解这个线性代数方程组可求出结位移.我们已知,总体刚度矩阵具有大型、对称、稀疏、带状分布、正定、主元占优势的特点,稀疏表示将对称消元法进一步改造,使之适合总刚的等带宽二维存储.(4)因子化法(三角分解)又称Cholesky分解,适合一维变带宽存储总刚.这上方法储效率高,计算速度快,应用较为普遍.此外,还有一种方法,叫做波前法.波前法实际上也是一种改进的高斯消去法.它建立一个称为“波前”的空间,各单元刚度系数依次进入波前.一旦与某自由度有关的所有单元的刚度系数全部装入,便可将相应的变量消去.经过消元的方程的系数随即退出波前,存放在计算机的外存中.这样就可腾出空间装入新的刚度系数.所以,波前法不需要生成完整的总刚,而是边组装边消元,“成熟”一个消去一个.消元完成后,全部系数都已存储在计算机的外存或缓冲区中.回代时将各方程的系数按“先出后入”的顺序调入内存求解.由此可见,这种方法是利用计算机充裕的外存资源,以多耗取机时来缓解内存不足的矛盾,以便适应较大规模的问题.随着计算机技术的发展,内存资源不断扩大,对具有稀疏、带状性质的有限元刚度方程,这种以时间换取空间的办法得不偿失.另一方面,波前法的阐述和程序设计比较复杂,且对多种单元并存的结构使用不便.所以,本书不拟介绍波前法.本书第九章将详细讨论适合整体存储总刚的高斯消去法和适合一维变带宽存储的因子化法以及有关的程序设计问题,以下仅列出这两种方法的梗概.1、高斯消去法高斯循序消去法的一般公式:对于n阶线性代数方程,需进行次消元.采用循序消去时,第m次消元以m-1次消元后的m行元素作为主元行,为主元,对第行元素()的消元公式为(3-21)式中等的上角码(m),表示该元素是经过第m次消元后得到的结果.同样,可以把经过m次消元后的系数矩阵和载荷阵分别记为及.式表时第m 次消元是在经m-1次消元的基础上进行的.消元过程中,主元及被消元素的位置可见图3-30(a).图中阴影部分已完成消元过程的元素,主元行以下的矩阵为待消部分.在进行第m次时,1-m行元素的消元过程已经完成,其中的元素就是消元最后得到的上三角阵中的元素. m行发下的元素消元过程尚未结束,连同m行元素在内构成一个待消的方阵.消元共需进行n-1次.消元完成后,即可回代求解.我们把消元最后结果记为,为上三角阵,回代公式可写作(3-22)回代过程自后向前进行.当回代求解时,已经解得.回代示意图见图3-30(b),阴影部分为已求得解答的部分.图7-30 高斯消去法2.三角分解法总体刚度平衡方程中,[K]是对称、正定矩阵,因而可做如下分解(3-23)其中,则是单位上三角矩阵,.代入整本结构平衡方程记,则.即由向下回代.由其中第一个方程解得,再由第二个方程解得,……,依此类推可求得{Y}.又由向上回代,可得,由得依此类推可求得.由上述过程可见,三角分解法求解线性代数方程组的关键是对系数矩阵进行三角分解.7.4.5 求解内力由平衡方程组解出位移后,从中分离出各单元的结点位移,再通过方程(3-3)、(3-4)和(3-6)等计算各单元的应变、应力和结点力等内力。
第2章5_用整体坐标表示单元刚度矩阵
0
K (2)
0.5
0
0.866 0
0
0
0
0 0.866 0.5 1 0
0 1
0
EA
0.5
0.866
0
0 2 0
0 0.866
0
0.5
0
0
0.5
0.866
0
0
0
0
0
0 0.5 0.866
0.75
0.433
0.75
0.433
0.433 0.25 0.433 0.25
C
2 y
y
E
A
CxCy
C
2 y
l
— 上式即为平面桁架单元整体坐标表示的单元刚度矩阵
[例2-1]平面桁架如图所示,各杆截面 EA均为常数。已知P1=15kN,P2= 20kN,试桁架各杆轴力。
1.对结点和单元编号如图示; 2. 列表表示各单元参数;
单元 ① ②
单元坐标 x轴方向
1→2
3→2
α
0
0
0.04 0.12
0
0.04 0.12
K
(1)
K
(3)
0
0.48
0 4
0.12 0
0.24 0
105
对 称
0.04 0.12
0.48
单元(2)的单元坐标和整体坐标不一致,必须经过以下变换
第一种方法: 直接代入公式:
2 1 2i 2 BCx l2 Cy
(e)
K
1 2i (B l2 )CxC y
0
6EI
2
l 2EI
l
EA l
0
0
EA l 0
7.4单元刚度矩阵组装及整体分析报告材料
7.4 单元刚度矩阵组装及整体分析7.4.1 单刚组装形成总刚根据全结构的平衡方程可知,总体刚度矩阵是由单元刚度矩阵集合而成的.如果一个结构的计算模型分成个单元,那么总体刚度矩阵可由各个单元的刚度矩阵组装而成,即[K]是由每个单元的刚度矩阵的每个系数按其脚标编号“对号入座”叠加而成的.这种叠加要求在同一总体坐标系下进行.如果各单元的刚度矩阵是在单元局部坐标下建立的,就必须要把它们转换到统一的结构(总体)坐标系.将总体坐标轴分别用表示,对某单元有式中,和分别是局部坐标系和总体坐标系下的单元结点位移向量;[T]为坐标转换阵,仅与两个坐标系的夹角有关,这样就有是该单元在总体坐标系下的单元刚度矩阵.以后如不特别强调,总体坐标系下的各种物理参数均不加顶上的横杠.下面就通过简单的例子来说明如何形成总体刚度矩阵.设有一个简单的平面结构,选取6个结点,划分为4个单元.单元及结点编号如图3-27所示.每个结点有两个自由度.总体刚度矩阵的组装过程可分为下面几步:图7-27(1)按单元局部编号顺序形成单元刚度矩阵.图7-27中所示的单元③,结点的局部编号顺序为.形成的单元刚度矩阵以子矩阵的形式给出是(2)将单元结点的局部编号换成总体编号,相应的把单元刚度矩阵中的子矩阵的下标也换成总体编号.对下图3-27所示单元③的刚度矩阵转换成总体编号后为(3)将转换后的单元刚度矩阵的各子矩阵,投放到总体刚度矩阵的对应位置上.单元③的各子矩阵投放后情况如下:(4)将所有的单元都执行上述的1,2,3步,便可得到总体刚度矩阵,如式(3-9).其中右上角的上标表示第单元所累加上的子矩阵.(3-9)(5)从式(3-9)可看出,总体刚度矩阵中的子矩阵AB是单元刚度矩阵的子矩阵转换成总体编号后具有相同的下标,的那些子矩阵的累加.总体刚度矩阵第行的非零子矩阵是由与结点相联系的那些单元的子矩阵向这行投放所构成的.7.4.2 结点平衡方程我们首先用结构力学方法建立结点平衡方程.连续介质用有限元法离散以后,取出其中任意一个结点,从环绕点各单元移置而来的结点载荷为式中表示对环绕结点的所有单元求和,环绕结点的各单元施加于结点的结点力为.因此,结点的平衡方程可表示为(3-10)以[K]代入平衡方程,得到以结点位移表示的结点的平衡方程,对于每个结点,都可列出平衡方程,于是得到整个结构的平衡方程组如下:式中,[K]为整体刚度矩阵,为全部结点位移组成的向量,为全部结点载荷组成的向量.当然,如果各点的载荷向量也是在单元局部坐标下建立的,在合成以前,也应把它们转换到统一的结构(总体)坐标系下,即式中,是总体坐标系下的结点载荷向量,为坐标转换阵.7.4.3 位移边界条件在有限元法对结构进行整体分析时,建立了整体刚度矩阵[K],也得到了结构的刚度平衡方程,即.结构刚度方程的求解相当于总刚[K]求逆的过程.但是,从数学上看,未经处理的总刚是对称、半正定的奇异矩阵,它的行列式值为零,不能立即求逆.从物理意义看,在进行整体分析时,结构是处于自由状态,在结点载荷的作用下,结构可以产生任意的刚体位移.所以,在已知结点载荷的条件下,仍不能通过平衡方程惟一地解出结点位移.为了使问题可解,必须对结构加以足够的位移约束,也就是应用位移边界条件.首先要通过施加适当的约束,消除结构的钢体位移,再根据问题要求设定其他已知位移.所以,处理位移边界条件在有限元分析步骤中十分重要.约束的种类包括使某些自由度上位移为零,,或给定其位移值,还有给定支承刚度等,本书涉及前两种.处理约束的方法,常用的有删行删列法、分块法、置大数法和置“1”法等,下面分别予以介绍.1、删行删列法若结构的某些结点位移值为零时(即与刚性支座连接点的位移),则可将总体刚度矩阵中相应的行列、删行删列划掉,然后将矩阵压缩即可求解.这种方法的优点是道理简单.如果删去的行列很多,则总体刚度矩阵的阶数可大大缩小.通常用人工计算时常采用该方法.若用计算机算题,在程序编制上必带来麻烦,因为刚度矩阵压缩以后,刚度矩阵中各元素的下标必全改变.因而一般计算机算题不太采用.2.分块法为了理解这个方法,我们把方程分块如下:(3-11)其中,假设是给定的结点位移;是无约束的(自由)结点位移.因而是已知的结点力;是未知的结点力.方程(3-11)可以写为即(3-12)和(3-13)其中,不是奇异的,因而可以解方程(3-12)得出(3-14)一旦知道了,就可以由方程(3-13)求得未知结点力.在全部给定的结点自由度都等于零的特殊情况下,我们可以删除对应于的各行和各列(即删行删列法),故可把方程简写为(3-15)3.置“1”法由于全部给定的结点位移通常都不能在位移向量的开始或终了,故分块法的编号方法是很麻烦的.因此,为了引入给定的边界条件,可以采用下述等价的方法.可以把方程(3-12)和(3-13)合在一起写为(3-16)在实际计算中,方程(3-16)所示的过程可以在不重新排列所述方程的情况下用下述分块的方法为进行.步骤(1)如果把给定为,则载荷向量P可以修改为为结点自由度总数.步骤(2)除对角线元素以外,使[K]中对应于的行和列为零,而对角线元素为1,即步骤(3)在载荷向量中引入规定的值,即对全部规定的结点位移均应反复运用上述过程(步骤(1)到(3)).应当指出,由于这个过程保持了方程的对称性,因此,[K]可以按带状存储,而且几乎不会增加编制程序的工作量.4.置大数法置大数法的思路是:在总体刚度矩阵中,把指定位移所对应的行和列的对角元素乘上一个很大的数,如,此行其他元素保持不变,同时把该行对应的载荷项也相应地用来代替,这里为指定位移,于是原平衡方程组变为除第行外,其他各行仍保持原来的平衡特性,而第个方程式展开为由于上式中的比其他项的系数大得多,求和后可略去其小量,则上式变为即.这样就用近似方程组代替原方程组,得到近似满足边界条件的解.当指定位移为零时,只要将对角元素乘上一个大数,而相应的载荷项经证明可以不置零.删行删列法适用于指定零位移点,而置大数法适用于给定位移(包括零位移).5.斜支座的处理对于简单的约束情况(如限定某些结点位移为零或取得给定数值),可以用前述置大数法处理.有的结构在直角坐标系内建立了位移方程组,但在某个斜边上受有法向约束.如图3-28所示正方形固支板,受均布横向载荷,对此,可利用对称性而只计算其1/8,如图中ABC部分,其中AC为固支边,按对称性,AB边上有,但在BC边上应限定绕BC的转用等于零.为处理此类斜边上的约束,须对斜边上的结点做坐标变换.若结构的总体坐标系为为斜支座的局部坐标系(见图3-29).对于边界结点,须限定方向位移,为此,将边界结点的位移及载荷都变换到局部坐标轴系.设轴与斜支座的轴夹角为,逆时针为正,图7-28 图7-29 则依据第二单中坐标转换关系有其中,.或写成(3-17)与位移关系相同有(3-18)将上两式带入结构刚度方程有(3-19)这样把位移到列阵中凡是斜支座的结点位移矢量都用局部坐标表示了.将式(3-19)中第行左右两边前乘以(3-20)由上式可见:凡是边界点的斜支座,在刚度方程中对应于斜支座的位移和载荷向量均可直接斜支座的局部坐标值,总刚度距阵中的相应行列需作相应的变换.上式的系数矩阵仍然是对称的,而且此方程中结点位沿轴表示,这样,限定方向的位移就很方便了.实际计算中,并不需要建立结构总的位移方程组后再进坐标变换.而可以在形成单元刚度矩阵和结点载荷之后,就对斜支座点进行坐标变换,把变换后的单元刚度矩阵和结点载荷叠加入总刚度矩阵和总载荷的相应位置,最后叠加形成的也就是方程组(3-20),即需要处理的结点,应该在单元计算中完成坐标变换后再叠加,当结构有不同的斜边约束时,都可以这样处理,只不过对不同边上的结点,应按不同的方向余弦矩阵变换就是了.7.4.4 总刚度平衡方程的求解应用有限元法,最终都是归结为解总体刚度平衡方程,它实际上是以总体刚度矩阵为系数矩阵的大型线性代数方程组.通过对结构施加位移边界条件,消除了结构的刚体位移,从而消除总体刚度矩阵的奇异性,解这个线性代数方程组可求出结位移.我们已知,总体刚度矩阵具有大型、对称、稀疏、带状分布、正定、主元占优势的特点,稀疏表示将对称消元法进一步改造,使之适合总刚的等带宽二维存储.(4)因子化法(三角分解)又称Cholesky分解,适合一维变带宽存储总刚.这上方法储效率高,计算速度快,应用较为普遍.此外,还有一种方法,叫做波前法.波前法实际上也是一种改进的高斯消去法.它建立一个称为“波前”的空间,各单元刚度系数依次进入波前.一旦与某自由度有关的所有单元的刚度系数全部装入,便可将相应的变量消去.经过消元的方程的系数随即退出波前,存放在计算机的外存中.这样就可腾出空间装入新的刚度系数.所以,波前法不需要生成完整的总刚,而是边组装边消元,“成熟”一个消去一个.消元完成后,全部系数都已存储在计算机的外存或缓冲区中.回代时将各方程的系数按“先出后入”的顺序调入内存求解.由此可见,这种方法是利用计算机充裕的外存资源,以多耗取机时来缓解内存不足的矛盾,以便适应较大规模的问题.随着计算机技术的发展,内存资源不断扩大,对具有稀疏、带状性质的有限元刚度方程,这种以时间换取空间的办法得不偿失.另一方面,波前法的阐述和程序设计比较复杂,且对多种单元并存的结构使用不便.所以,本书不拟介绍波前法.本书第九章将详细讨论适合整体存储总刚的高斯消去法和适合一维变带宽存储的因子化法以及有关的程序设计问题,以下仅列出这两种方法的梗概.1、高斯消去法高斯循序消去法的一般公式:对于n阶线性代数方程,需进行次消元.采用循序消去时,第m次消元以m-1次消元后的m行元素作为主元行,为主元,对第行元素()的消元公式为(3-21)式中等的上角码(m),表示该元素是经过第m次消元后得到的结果.同样,可以把经过m次消元后的系数矩阵和载荷阵分别记为及.式表时第m 次消元是在经m-1次消元的基础上进行的.消元过程中,主元及被消元素的位置可见图3-30(a).图中阴影部分已完成消元过程的元素,主元行以下的矩阵为待消部分.在进行第m次时,1-m行元素的消元过程已经完成,其中的元素就是消元最后得到的上三角阵中的元素. m行发下的元素消元过程尚未结束,连同m行元素在内构成一个待消的方阵.消元共需进行n-1次.消元完成后,即可回代求解.我们把消元最后结果记为,为上三角阵,回代公式可写作(3-22)回代过程自后向前进行.当回代求解时,已经解得.回代示意图见图3-30(b),阴影部分为已求得解答的部分.图7-30 高斯消去法2.三角分解法总体刚度平衡方程中,[K]是对称、正定矩阵,因而可做如下分解(3-23)其中,则是单位上三角矩阵,.代入整本结构平衡方程记,则.即由向下回代.由其中第一个方程解得,再由第二个方程解得,……,依此类推可求得{Y}.又由向上回代,可得,由得依此类推可求得.由上述过程可见,三角分解法求解线性代数方程组的关键是对系数矩阵进行三角分解.7.4.5 求解内力由平衡方程组解出位移后,从中分离出各单元的结点位移,再通过方程(3-3)、(3-4)和(3-6)等计算各单元的应变、应力和结点力等内力。
单元类型及单元刚度矩阵课件
面积单元的刚度矩阵可以通过解析方 法或数值方法计算得到。
它具有四个节点,每个节点具有三个 自由度:x、y和z方向的位移。
体积单元
体积单元是一种几何 形状,通常用于模拟 结构中的三维实体或 区域。
体积单元的刚度矩阵 可以通过解析方法或 数值方法计算得到。
它具有八个节点,每 个节点具有三个自由 度:x、y、z方向的 位移。
移。
线性单元的刚度矩阵可以通过解 析方法或数值方法计算得到。
角点单元
角点单元是一种特殊类型的线 性单元,通常用于模拟结构中 的角点或连接两个线性单元的 节点。
它具有三个自由度:x、y和z方 向的位移。
角点单元的刚度矩阵可以通过 解析方法或数值方法计算得到。
面积单元
面积单元是一种几何形状,通常用于 模拟结构中的平面区域或曲面上的小 区域。
单击此处添加正文,文字是您思想的提一一二三四五 六七八九一二三四五六七八九一二三四五六七八九文, 单击此处添加正文,文字是您思想的提炼,为了最终 呈现发布的良好效果单击此4*25}
通过稳定性分析,可以评估结构的承载安全性和预防 失稳的措施。
PART 04
单元类型选择与注意事项
选择依据
计算精度
根据模型精度要求选择合适的单 元类型,例如,对于复杂形状或 精细结构,应选择高阶单元以提
2023 WORK SUMMARY
单元类型及单元刚度 矩阵课件
REPORTING
CATALOGUE
• 单元类型介绍 • 单元刚度矩阵
PART 01
单元类型介绍
线性单元
线性单元是一种简单的几何形状, 通常用于模拟结构中的直线段或 平面区域。
它具有两个节点,每个节点具有 三个自由度:x、y和z方向的位
单元刚度矩阵组装及整体分析
单元刚度矩阵组装及整体分析在结构力学中,单元刚度矩阵的组装是进行有限元分析的重要步骤之一、通过将多个单元刚度矩阵组装成整体刚度矩阵,可以得到结构的整体刚度矩阵,并进行相应的整体分析。
下面将介绍单元刚度矩阵的组装方法以及整体分析的步骤。
1.单元刚度矩阵的组装方法:在有限元分析中,结构通常划分成多个单元,每个单元通过节点与相邻单元相连。
在单元分析中,首先需要建立单元刚度矩阵,然后将所有单元刚度矩阵组装成整体刚度矩阵。
单元刚度矩阵的建立通常通过离散化方法来进行,常见的方法有刚度法和能量法。
在刚度法中,通过对单元进行力学建模,并应用弹性力学原理,可以得到单元刚度矩阵。
在能量法中,通过考虑单元的应变能和变形能,可以得到单元刚度矩阵。
对于线性单元,其刚度矩阵可通过以下公式得到:[K]=∫(∑[B]T[D][B])dV其中,[K]为单元刚度矩阵,[B]为单元应变矩阵,[D]为单元弹性矩阵,dV为单元体积的微元。
在计算单元刚度矩阵时,通常会使用数值积分方法,如高斯积分,以提高计算精度。
2.整体分析的步骤:在得到所有单元的刚度矩阵后,需要将其组装成整体刚度矩阵并进行整体分析。
整体分析的步骤如下:(1)确定结构的边界条件:边界条件是指结构的位移或载荷边界条件,如固支条件、弹簧支承等。
在进行整体分析前,需要确定结构的边界条件。
(2)根据结构的几何特征和边界条件,建立结构的刚度矩阵方程:[K]{u}={F}其中,[K]为整体刚度矩阵,{u}为结构的位移向量,{F}为结构的载荷向量。
(3)施加约束条件:根据结构的边界条件,将约束施加到整体刚度矩阵方程中。
这可以通过修改整体刚度矩阵和载荷向量中的相应行和列来实现。
(4)解方程:通过求解经过约束的整体刚度矩阵方程,可以得到结构的位移。
(5)计算应力和应变:根据结构的位移和单元的形状,可以计算出每个单元的应力和应变。
这可以通过单元的位移-应变关系来实现。
(6)结果分析:通过计算得到的位移、应力和应变,可以对结构进行进一步的分析和评估,如变形分析、应力分析等。
单元刚度矩阵组装及整体分析
单元刚度矩阵组装及整体分析一、单元刚度矩阵的建立梁单元刚度矩阵:梁单元刚度矩阵是用来描述梁单元的弹性变形行为的矩阵。
常见的梁单元有线性梁单元和非线性梁单元。
对于线性梁单元来说,其刚度矩阵的计算可以通过在梁上进行数学推导得到。
三角形单元刚度矩阵:三角形单元刚度矩阵是用来描述三角形单元的弹性变形行为的矩阵。
常见的三角形单元有线性三角形单元和非线性三角形单元。
对于线性三角形单元来说,其刚度矩阵的计算可以通过在三角形单元上进行数学推导得到。
二、单元刚度矩阵的组装在结构的离散过程中,将整个结构划分为若干个单元,并按照一定的规则将单元刚度矩阵组装成整体刚度矩阵。
单元刚度矩阵的组装可以使用两种常见的方法:全局坐标法和局部坐标法。
全局坐标法:全局坐标法是一种将单元刚度矩阵直接组装到整体刚度矩阵中的方法。
在这种方法中,我们通过将单元的自由度与整体自由度进行对应,将单元的刚度矩阵的每个元素放入整体刚度矩阵的相应位置。
局部坐标法:局部坐标法是一种将单元刚度矩阵通过变换到整体坐标系后再进行组装的方法。
在这种方法中,我们首先将单元的自由度与局部坐标进行对应,然后将刚度矩阵变换到整体坐标系,最后再将变换后的刚度矩阵的每个元素放入整体刚度矩阵的相应位置。
三、整体分析在完成单元刚度矩阵的组装后,我们可以得到整体刚度矩阵。
整体刚度矩阵描述了整个结构的刚度特性,通过求解整体刚度矩阵与载荷之间的关系,可以得到结构的位移和应力分布。
对于线性弹性结构而言,整体分析可以通过直接求解线性方程组的方法进行。
我们可以根据边界条件和载荷信息,将整体刚度矩阵和载荷向量建立成一个线性方程组,然后通过数值方法(例如高斯消元法、LU分解法)求解出方程组的解,即得到结构的位移和应力分布。
对于非线性结构而言,整体分析可以采用迭代法进行。
在每一步迭代中,我们都需要更新刚度矩阵和载荷向量,然后再求解线性方程组,最终得到结构的位移和应力分布。
整体分析的目的是求解结构的位移和应力分布,进而评估结构的稳定性和安全性。
第三节刚度矩阵
第三节 刚度矩阵——节点载荷与节点位移之间的关系一、 单元刚度矩阵1. 单元刚度矩阵xj单元e 是在节点力作用下处于平衡。
节点i 的节点力为{}Ti xiyi R R R ⎡⎤=⎣⎦ (i , j , m 轮换)则单元e 的节点力列阵为{}TeTT T mi jTxm ym xi yi xj yj R R R R R R R R R R ⎡⎤⎣⎦⎡⎤⎣⎦==单元应力列阵为{}Tex y xy σσστ⎡⎤⎣⎦=假定弹性体的所有节点都产生一虚位移,单元e 的三个节点的虚位移为{}*******eT mm i i j ju v u v u v δ⎡⎤⎣⎦= 单元虚应变列阵为{}****Tx y xy εεεγ⎡⎤⎢⎥⎣⎦=参照式(3-7),则单元虚应变为{}{}**eeB εδ⎡⎤⎣⎦=作用在弹性体上的外力在虚位移上所做的功为:{}{}*eTe R δ⎛⎫ ⎪⎝⎭单元内的应力在虚应变上所做的功为:{}{}*Te tdxdy εσ∆⎛⎫ ⎪⎝⎭⎰⎰根据虚位移原理,可得单元的虚功方程{}{}{}{}**eTTe e R tdxdy δεσ∆⎛⎫⎛⎫⎪ ⎪⎝⎭⎝⎭=⎰⎰或{}{}{}{}**eTTTe e B R tdxdy δδσ∆⎛⎫⎛⎫⎡⎤ ⎪ ⎪⎣⎦⎝⎭⎝⎭=⎰⎰故有{}{}eTB R tdxdy σ∆⎡⎤⎣⎦=⎰⎰将式(3-10)代入,的{}{}{}eeeTTD B D B R B B tdxdytdxdy δδ∆∆⎡⎤⎡⎤⎡⎤⎣⎦⎣⎦⎣⎦⎡⎤⎡⎤⎡⎤⎣⎦⎣⎦⎣⎦==⎰⎰⎰⎰(3-27)简记为{}{}ee ek R δ⎡⎤⎣⎦= (3-29)--------上式表征单元节点力与节点位移之间的关系,称为单元刚度方程(单元平衡方程) 其中TeD B B k tdxdy ∆⎡⎤⎡⎤⎡⎤⎡⎤⎣⎦⎣⎦⎣⎦⎣⎦=⎰⎰(3-28) ek ⎡⎤⎣⎦称之为单元刚度矩阵(简称为单刚),是66⨯矩阵。
如果单元的材料是均质的,矩阵D ⎡⎤⎣⎦中的元素也是常量,且在三角形常应变的情况下,矩阵B ⎡⎤⎣⎦中的元素也是常数,当单元的厚度也是常数时,注意到dxdy ∆=∆⎰⎰,于是单元刚度矩阵可简化为TeB D B t k ⎡⎤⎡⎤⎡⎤⎡⎤⎣⎦⎣⎦⎣⎦⎢⎥⎣⎦∆= (3-30) 将单元刚度矩阵按节点号写成分块矩阵形式:66eii ij imji jj jm mm mimj kk k k k k k k k k ⨯⎡⎤⎢⎥⎡⎤⎢⎥⎢⎥⎢⎥⎣⎦⎢⎥⎢⎥⎣⎦= (3-31)其中任一子块[]rs k (r ,s=i ,j ,m )是一个2×2子矩阵,为[][][][]Trsrs k B D B t =∆ (r ,s=i ,j ,m )(1)对于平面应力问题将[]B 和平面应力问题的弹性矩阵[]D 代入,得Trs r s k B D B t ⎡⎤⎡⎤⎡⎤⎡⎤⎣⎦⎣⎦⎣⎦⎣⎦=∆ ()21122114122r s r s r s r s r s r s r s r s b b c c b c c b Et c b b cc c b b μμμμμμμ--⎡⎤++⎢⎥=⎢⎥---∆⎢⎥++⎢⎥⎣⎦(r ,s=i ,j ,m ) (3-32)(2)对于平面应变问题将[]B 和平面应变问题的弹性矩阵[]D 代入,得()()()()()()()12122112114112121212121e rs k b b c c b c c b r s r s r s r s E t c b b c c c b b r s r s r s r sμμμμμμμμμμμμμμμ⎡⎤⎢⎥⎢⎥⎡⎤⎢⎥⎣⎦⎢⎥⎢⎥⎣⎦--++----=+-∆--++--- (r,s=i ,j ,m ) (3-33)(注:是将式(3-32)中的,E μ分别换成21E μ- 和1μμ-)2. 单元刚度矩阵的性质 (1)ek ⎡⎤⎣⎦的物理意义式(3-29)可完整写为131415161112212223242526333435363132434445464142555152535456616263646i i jjmmeU k k k k k k V k k k k k k k k k k U k k k k k k k k V k k k k k k U k k k k k V⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦=⎧⎫⎪⎪⎪⎪⎪⎪⎨⎬⎪⎪⎪⎪⎪⎪⎩⎭566ii j j m m eu v u v u k v ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎡⎤⎢⎥⎢⎥⎢⎢⎥⎥⎣⎣⎦⎦⎧⎫⎪⎪⎪⎪⎪⎪⎨⎬⎪⎪⎪⎪⎪⎪⎩⎭可见每个节点在x 和y 方向上有二个平衡方程,3个节点共有六个平衡方程。
单元刚度矩阵组装及整体分析
7.4 单元刚度矩阵组装及整体分析7.4.1 单刚组装形成总刚根据全结构的平衡方程可知,总体刚度矩阵是由单元刚度矩阵集合而成的.如果一个结构的计算模型分成个单元,那么总体刚度矩阵可由各个单元的刚度矩阵组装而成,即[K]是由每个单元的刚度矩阵的每个系数按其脚标编号“对号入座”叠加而成的.这种叠加要求在同一总体坐标系下进行.如果各单元的刚度矩阵是在单元局部坐标下建立的,就必须要把它们转换到统一的结构(总体)坐标系.将总体坐标轴分别用表示,对某单元有式中,和分别是局部坐标系和总体坐标系下的单元结点位移向量;[T]为坐标转换阵,仅与两个坐标系的夹角有关,这样就有是该单元在总体坐标系下的单元刚度矩阵.以后如不特别强调,总体坐标系下的各种物理参数均不加顶上的横杠.下面就通过简单的例子来说明如何形成总体刚度矩阵.设有一个简单的平面结构,选取6个结点,划分为4个单元.单元及结点编号如图3-27所示.每个结点有两个自由度.总体刚度矩阵的组装过程可分为下面几步:图7-27(1)按单元局部编号顺序形成单元刚度矩阵.图7-27中所示的单元③,结点的局部编号顺序为.形成的单元刚度矩阵以子矩阵的形式给出是(2)将单元结点的局部编号换成总体编号,相应的把单元刚度矩阵中的子矩阵的下标也换成总体编号.对下图3-27所示单元③的刚度矩阵转换成总体编号后为(3)将转换后的单元刚度矩阵的各子矩阵,投放到总体刚度矩阵的对应位置上.单元③的各子矩阵投放后情况如下:(4)将所有的单元都执行上述的1,2,3步,便可得到总体刚度矩阵,如式(3-9).其中右上角的上标表示第单元所累加上的子矩阵.(3-9)(5)从式(3-9)可看出,总体刚度矩阵中的子矩阵AB是单元刚度矩阵的子矩阵转换成总体编号后具有相同的下标,的那些子矩阵的累加.总体刚度矩阵第行的非零子矩阵是由与结点相联系的那些单元的子矩阵向这行投放所构成的.7.4.2 结点平衡方程我们首先用结构力学方法建立结点平衡方程.连续介质用有限元法离散以后,取出其中任意一个结点,从环绕点各单元移置而来的结点载荷为式中表示对环绕结点的所有单元求和,环绕结点的各单元施加于结点的结点力为.因此,结点的平衡方程可表示为(3-10)以[K]代入平衡方程,得到以结点位移表示的结点的平衡方程,对于每个结点,都可列出平衡方程,于是得到整个结构的平衡方程组如下:式中,[K]为整体刚度矩阵,为全部结点位移组成的向量,为全部结点载荷组成的向量.当然,如果各点的载荷向量也是在单元局部坐标下建立的,在合成以前,也应把它们转换到统一的结构(总体)坐标系下,即式中,是总体坐标系下的结点载荷向量,为坐标转换阵.7.4.3 位移边界条件在有限元法对结构进行整体分析时,建立了整体刚度矩阵[K],也得到了结构的刚度平衡方程,即.结构刚度方程的求解相当于总刚[K]求逆的过程.但是,从数学上看,未经处理的总刚是对称、半正定的奇异矩阵,它的行列式值为零,不能立即求逆.从物理意义看,在进行整体分析时,结构是处于自由状态,在结点载荷的作用下,结构可以产生任意的刚体位移.所以,在已知结点载荷的条件下,仍不能通过平衡方程惟一地解出结点位移.为了使问题可解,必须对结构加以足够的位移约束,也就是应用位移边界条件.首先要通过施加适当的约束,消除结构的钢体位移,再根据问题要求设定其他已知位移.所以,处理位移边界条件在有限元分析步骤中十分重要.约束的种类包括使某些自由度上位移为零,,或给定其位移值,还有给定支承刚度等,本书涉及前两种.处理约束的方法,常用的有删行删列法、分块法、置大数法和置“1”法等,下面分别予以介绍.1、删行删列法若结构的某些结点位移值为零时(即与刚性支座连接点的位移),则可将总体刚度矩阵中相应的行列、删行删列划掉,然后将矩阵压缩即可求解.这种方法的优点是道理简单.如果删去的行列很多,则总体刚度矩阵的阶数可大大缩小.通常用人工计算时常采用该方法.若用计算机算题,在程序编制上必带来麻烦,因为刚度矩阵压缩以后,刚度矩阵中各元素的下标必全改变.因而一般计算机算题不太采用.2.分块法为了理解这个方法,我们把方程分块如下:(3-11)其中,假设是给定的结点位移;是无约束的(自由)结点位移.因而是已知的结点力;是未知的结点力.方程(3-11)可以写为即(3-12)和(3-13)其中,不是奇异的,因而可以解方程(3-12)得出(3-14)一旦知道了,就可以由方程(3-13)求得未知结点力.在全部给定的结点自由度都等于零的特殊情况下,我们可以删除对应于的各行和各列(即删行删列法),故可把方程简写为(3-15)3.置“1”法由于全部给定的结点位移通常都不能在位移向量的开始或终了,故分块法的编号方法是很麻烦的.因此,为了引入给定的边界条件,可以采用下述等价的方法.可以把方程(3-12)和(3-13)合在一起写为(3-16)在实际计算中,方程(3-16)所示的过程可以在不重新排列所述方程的情况下用下述分块的方法为进行.步骤(1)如果把给定为,则载荷向量P可以修改为为结点自由度总数.步骤(2)除对角线元素以外,使[K]中对应于的行和列为零,而对角线元素为1,即步骤(3)在载荷向量中引入规定的值,即对全部规定的结点位移均应反复运用上述过程(步骤(1)到(3)).应当指出,由于这个过程保持了方程的对称性,因此,[K]可以按带状存储,而且几乎不会增加编制程序的工作量.4.置大数法置大数法的思路是:在总体刚度矩阵中,把指定位移所对应的行和列的对角元素乘上一个很大的数,如,此行其他元素保持不变,同时把该行对应的载荷项也相应地用来代替,这里为指定位移,于是原平衡方程组变为除第行外,其他各行仍保持原来的平衡特性,而第个方程式展开为由于上式中的比其他项的系数大得多,求和后可略去其小量,则上式变为即.这样就用近似方程组代替原方程组,得到近似满足边界条件的解.当指定位移为零时,只要将对角元素乘上一个大数,而相应的载荷项经证明可以不置零.删行删列法适用于指定零位移点,而置大数法适用于给定位移(包括零位移).5.斜支座的处理对于简单的约束情况(如限定某些结点位移为零或取得给定数值),可以用前述置大数法处理.有的结构在直角坐标系内建立了位移方程组,但在某个斜边上受有法向约束.如图3-28所示正方形固支板,受均布横向载荷,对此,可利用对称性而只计算其1/8,如图中ABC部分,其中AC为固支边,按对称性,AB边上有,但在BC边上应限定绕BC的转用等于零.为处理此类斜边上的约束,须对斜边上的结点做坐标变换.若结构的总体坐标系为为斜支座的局部坐标系(见图3-29).对于边界结点,须限定方向位移,为此,将边界结点的位移及载荷都变换到局部坐标轴系.设轴与斜支座的轴夹角为,逆时针为正,图7-28 图7-29 则依据第二单中坐标转换关系有其中,.或写成(3-17)与位移关系相同有(3-18)将上两式带入结构刚度方程有(3-19)这样把位移到列阵中凡是斜支座的结点位移矢量都用局部坐标表示了.将式(3-19)中第行左右两边前乘以(3-20)由上式可见:凡是边界点的斜支座,在刚度方程中对应于斜支座的位移和载荷向量均可直接斜支座的局部坐标值,总刚度距阵中的相应行列需作相应的变换.上式的系数矩阵仍然是对称的,而且此方程中结点位沿轴表示,这样,限定方向的位移就很方便了.实际计算中,并不需要建立结构总的位移方程组后再进坐标变换.而可以在形成单元刚度矩阵和结点载荷之后,就对斜支座点进行坐标变换,把变换后的单元刚度矩阵和结点载荷叠加入总刚度矩阵和总载荷的相应位置,最后叠加形成的也就是方程组(3-20),即需要处理的结点,应该在单元计算中完成坐标变换后再叠加,当结构有不同的斜边约束时,都可以这样处理,只不过对不同边上的结点,应按不同的方向余弦矩阵变换就是了.7.4.4 总刚度平衡方程的求解应用有限元法,最终都是归结为解总体刚度平衡方程,它实际上是以总体刚度矩阵为系数矩阵的大型线性代数方程组.通过对结构施加位移边界条件,消除了结构的刚体位移,从而消除总体刚度矩阵的奇异性,解这个线性代数方程组可求出结位移.我们已知,总体刚度矩阵具有大型、对称、稀疏、带状分布、正定、主元占优势的特点,稀疏表示将对称消元法进一步改造,使之适合总刚的等带宽二维存储.(4)因子化法(三角分解)又称Cholesky分解,适合一维变带宽存储总刚.这上方法储效率高,计算速度快,应用较为普遍.此外,还有一种方法,叫做波前法.波前法实际上也是一种改进的高斯消去法.它建立一个称为“波前”的空间,各单元刚度系数依次进入波前.一旦与某自由度有关的所有单元的刚度系数全部装入,便可将相应的变量消去.经过消元的方程的系数随即退出波前,存放在计算机的外存中.这样就可腾出空间装入新的刚度系数.所以,波前法不需要生成完整的总刚,而是边组装边消元,“成熟”一个消去一个.消元完成后,全部系数都已存储在计算机的外存或缓冲区中.回代时将各方程的系数按“先出后入”的顺序调入内存求解.由此可见,这种方法是利用计算机充裕的外存资源,以多耗取机时来缓解内存不足的矛盾,以便适应较大规模的问题.随着计算机技术的发展,内存资源不断扩大,对具有稀疏、带状性质的有限元刚度方程,这种以时间换取空间的办法得不偿失.另一方面,波前法的阐述和程序设计比较复杂,且对多种单元并存的结构使用不便.所以,本书不拟介绍波前法.本书第九章将详细讨论适合整体存储总刚的高斯消去法和适合一维变带宽存储的因子化法以及有关的程序设计问题,以下仅列出这两种方法的梗概.1、高斯消去法高斯循序消去法的一般公式:对于n阶线性代数方程,需进行次消元.采用循序消去时,第m次消元以m-1次消元后的m行元素作为主元行,为主元,对第行元素()的消元公式为(3-21)式中等的上角码(m),表示该元素是经过第m次消元后得到的结果.同样,可以把经过m次消元后的系数矩阵和载荷阵分别记为及.式表时第m 次消元是在经m-1次消元的基础上进行的.消元过程中,主元及被消元素的位置可见图3-30(a).图中阴影部分已完成消元过程的元素,主元行以下的矩阵为待消部分.在进行第m次时,1-m行元素的消元过程已经完成,其中的元素就是消元最后得到的上三角阵中的元素. m行发下的元素消元过程尚未结束,连同m行元素在内构成一个待消的方阵.消元共需进行n-1次.消元完成后,即可回代求解.我们把消元最后结果记为,为上三角阵,回代公式可写作(3-22)回代过程自后向前进行.当回代求解时,已经解得.回代示意图见图3-30(b),阴影部分为已求得解答的部分.图7-30 高斯消去法2.三角分解法总体刚度平衡方程中,[K]是对称、正定矩阵,因而可做如下分解(3-23)其中,则是单位上三角矩阵,.代入整本结构平衡方程记,则.即由向下回代.由其中第一个方程解得,再由第二个方程解得,……,依此类推可求得{Y}.又由向上回代,可得,由得依此类推可求得.由上述过程可见,三角分解法求解线性代数方程组的关键是对系数矩阵进行三角分解.7.4.5 求解内力由平衡方程组解出位移后,从中分离出各单元的结点位移,再通过方程(3-3)、(3-4)和(3-6)等计算各单元的应变、应力和结点力等内力。
matlab组装整体刚度矩阵
matlab组装整体刚度矩阵【原创实用版】目录一、引言二、组装整体刚度矩阵的方法1.逐个提取单元刚度矩阵2.合并单元刚度矩阵3.保存整体刚度矩阵三、MATLAB 编程实例1.创建刚度矩阵2.组装整体刚度矩阵3.保存整体刚度矩阵四、结论正文在结构分析和计算中,我们常常需要组装整体刚度矩阵,以研究结构的整体刚度特性。
MATLAB 作为一种广泛应用于科学计算和工程分析的软件,为我们提供了方便的矩阵操作和处理功能。
本文将介绍如何使用MATLAB 编程组装整体刚度矩阵。
一、引言在结构力学中,刚度矩阵是描述结构刚度特性的重要工具。
在实际应用中,我们通常需要将多个单元刚度矩阵组装成一个整体刚度矩阵,以计算结构的整体刚度。
MATLAB 提供了丰富的矩阵操作函数,可以方便地进行刚度矩阵的组装和处理。
二、组装整体刚度矩阵的方法组装整体刚度矩阵的方法主要包括以下三个步骤:1.逐个提取单元刚度矩阵在 MATLAB 中,我们可以使用提取单元刚度矩阵的函数,如ktriangle2d3nodestiffness,来计算每个单元的刚度矩阵。
该函数接受弹性模量 e、泊松比 nu、厚度 t 以及节点坐标 xi、yi、xj、yj、xm、ym 等参数,并计算出单元刚度矩阵。
2.合并单元刚度矩阵在提取出各个单元的刚度矩阵后,我们需要将它们合并成一个整体刚度矩阵。
在 MATLAB 中,我们可以使用循环结构逐个提取单元刚度矩阵,并将它们合并成一个大的刚度矩阵。
3.保存整体刚度矩阵在 MATLAB 中,我们可以使用 save 函数将整体刚度矩阵保存到一个文件中,以便后续进行结构分析和计算。
三、MATLAB 编程实例下面是一个简单的 MATLAB 编程实例,演示如何组装整体刚度矩阵:1.创建刚度矩阵假设我们有一个由三个节点构成的平面三角形单元,其弹性模量e=200 GPa,泊松比 nu=0.3,厚度 t=1 cm。
刚度矩阵法获奖课件
• 基本概念 • 位移模式 • 刚度矩阵 • 座标变换 • 基本环节
§2.1基本概念
• 刚度矩阵法旳基本思绪: 将构造划分为有限个单元
外力——内力(应力)——应变——位移——节点位移
位移模式 几何方程 本构方程
平衡方程
θj θi
L+ΔL
i
j
L
• 单元:构成构造旳基本部分,是人为划 分,对于杆系,一般一种杆件为一种单 元
K e e f e f TF
K eTe TF e
T T K eTe F e Κ e T T K eT Κ ee F e
总刚度方程旳建立:
Κ ee F e
Κ F
Κ Κe
•单元刚度矩阵旳建立 •引入边界条件 •求解刚度方程组 •求解应力
x y
ml xx
ly my
lz mz
x y
z nx ny nz z
x x cos(x, x' ) y cos(x' , y) z cos(x' , z)
y x cos( y ' , x) y cos( y, y ' ) z cos( y ' , z)
*T
*T
BT DBdV f
V
*T
*T
BT DBdV f
V
*T
*T
BT DBdV f V
BT DBdV f
V
BT DBdV f
V
K f
K BT DBdV V
§2.4 坐标变换
Z Y’
Z’
Y X’ X
x y z T
x y z T
x x cos(x, x' ) y cos(x, y ' ) z cos(x, z ' )
整体分析及总体刚度矩阵的性质
整体刚度矩阵的特点
2、稀疏性。
1 1 2 3 4 5 6 7 8 9 10 2 3 4 5 6 7 8 9 10
一般,一个节点的相关结 点不会超过九个,如果网格中 有200个节点,则一行中非零 子块的个数与该行的子块总数 相比不大于9/200,即在5%以 下,如果网格的节点个数越多, 则刚度矩阵的稀疏性就越突出。 利用矩阵[K]的稀疏性, 可设法只存贮非零元素,从而 可大量地节省存贮容量。
整体分析
2、根据支承条件修改整体刚度矩阵。 建立整体刚度矩阵时,每个节点的位移当作未知量看待,没有
考虑具体的支承情况,因此进行整体分析时还要针对支承条件加以 处理。
在上图的结构中,支承条件共有四个,即在节点1、4、6的四 个支杆处相应位移已知为零:
u1 0,u4 0,v4 0,v6 0
局部码 总码
j1
1
[K jj ](1)
m1 , j 2 , i 3
2
[ K jm ](1)
i1 , m 3 , j 4
3
[ K ji ](1)
m2
4
i 2 , j3 , m 4
5
i4
6
j1
m1 j2 i3
1
[K mm ](1) 2
[K jj ]
[K ii ]
( 2)
( 3)
[K mi ](1) [K im ]( 3) [K ii ](1)
K K
集成包含搬家和迭加两个环节: e A、将单元刚度矩阵 K 中的子块搬家,得出单元的扩 e 大刚度矩阵 。 K e B、将各单元的扩大刚度矩阵 K 迭加,得出结构刚度 矩阵[K]。 T T n 2) R R1 Rn 为节点载荷向量, 1 为节点位移向量。
单元刚度矩阵PPT课件
m x
最新课件
29
性质1 形函数Ni在节点i上的值等于1,在其它节点上的 值等于0。对于本单元,有:
u ( x , y ) N i u i N j u j N m u m( x , y ) N ii N jj N m m
N i( x i,y i) 1 N i( x j,y j) 0 N i( x m ,y m ) 0
令
1
Ni 2A(ai bixciy) (i, j,m)
(2-18)
位移模式(2-16)可以简写为
u N iu i N ju j N m u m N ii N jj N m m (2-19)
式(2-19)中的Ni、Nj、Nm是坐标的函数,反应了单 元的位移形态,称为单元位移函数的形函数。数学上它反
u v
j j
u
m
vm
F 最新课件k
F xi
F
yi
F
F xj F yj
F F
xm ym
3
第二章 单元分析
——平面问题常应变单元
本章主要讲单元分析的一般理论、方法。但为了便 于理解,以平面问题常应变三角形单元为对象进行说明、 演引。必须指出:尽管说明、演引中具有明显的针对性 (平面问题三角形单元),但原理、方法和主要矩阵公 式都具有普遍性。
位移函数中包含了单元的常应变。
x u x,y y v,最x新y 课件 u y x v (a2, a6, a3+a158)
①、②、③、④单元的位移函数都是
u a 1 a 2 x a 3 y a 4 a 5 x a 6 y
可以看出: 位移函数在单元内是连续的; 位移函数在单元之间的边界上也连续吗?是。
ui uj
最新7.4-单元刚度矩阵组装及整体分析
根据全结构的平衡方程可知,总体刚度矩阵是由单元刚度矩阵集合而成的个结构的计算模型分成个单元,那么总体刚度矩阵可由各个单元的刚度矩阵组装而成,即是由每个单元的刚度矩阵的每个系数按其脚标编号“对号入座”叠加而成的将总体坐标轴分别用表示,对某单元有式中,和分别是局部坐标系和总体坐标系下的单元结点位移向量是该单元在总体坐标系下的单元刚度矩阵.将单元结点的局部编号换成总体编号,其中右上角的上标表示第单元所累加上的子矩阵具有相同的下标,的那些子矩阵的累加总体刚度矩阵第行的非零子矩阵是由与结点相联系的那,从环绕点各单元移置而来的结点载荷为式中表示对环绕结点的所有单元求和,环绕结点的各单元施加于结点的结点力为.因此,结点的平衡方程可表示为得到以结点位移表示的结点的平衡方程,为整体刚度矩阵,为全部结点位移组成的向量,为全部结点载荷组成的向量式中,是总体坐标系下的结点载荷向量,为坐标转换阵.构是处于自由状态,在结点载荷的作用下,结构可以产生任意的刚体位移的条件下,仍不能通过平衡方程惟一地解出结点位移.约束的种类包括使某些自由度上位移为零,,或给定其位移值,还有给定支承刚为了理解这个方法,我们把方程分块如下:其中,假设是给定的结点位移;是无约束的(自由)结点位移因而是已知的结点力;其中,不是奇异的,因而可以解方程(一旦知道了,求得未知结点力.殊情况下,我们可以删除对应于的各行和各列(即删行删列法),故可把方程简写为由于全部给定的结点位移通常都不能在位移向量的开始或终了,故分块法的编号方法是很麻烦因此,为了引入给定的边界条件,可以采用下述等价的方法如果把给定为,则载荷向量为结点自由度总数中对应于的行和列为零,而对角线元素为)在载荷向量中引入规定的值,即对全部规定的结点位移均应反复运用上述过程(步骤(置大数法的思路是:在总体刚度矩阵中,把指定位移所对应的行和列的对角元素乘上一个很大的数,如,此行其他元素保持不变,同时把该行对应的载荷项也相应地用来代替,这里为指定位移,于是原平衡方程组变为除第行外,其他各行仍保持原来的平衡特性,而第个方程式展开为由于上式中的比其他项的系数大得多,求和后可略去其小量,则上式变为即.边上有,若结构的总体坐标系为为斜支座的局部坐标系(见图对于边界结点,须限定方向位移,为此,将边界结点的位移及载荷都变换到局部坐标轴系设轴与斜支座的轴夹角为,逆时针为正,其中,.)中第行左右两边前乘以上式的系数矩阵仍然是对称的,而且此方程中结点位沿轴表示,这样,限定方向的位移异性,解这个线性代数方程组可求出结位移.阶线性代数方程,需进行次消元行元素作为主元行,为主元,对第行元素()的消元公式为式中等的上角码(次消元后的系数矩阵和载荷阵分别记为及.式表时第我们把消元最后结果记为,为上当回代求解时,已经解得总体刚度平衡方程中,,是单位上三角矩阵,.记,则.由其中第一个方程解得,再由第二个方程解得,向上回代,可得,由得依此类推可求得.由平衡方程组解出位移后,从中分离出各单元的结点位移,再通过方程)等计算各单元的应变、应力和结点力等内力。
平面分析整体刚度矩阵
F (4) 4
0
0 0
0 0 0
0 K (4)
33
K (4) 43
0 K (4)
34
K (4) 44
0 K (4)
35
K (4) 45
2 3
4
K
(4)
F (4) 5
0
0
K (4) 53
K (4) 54
K
(4) 55
5
式中:
Fi ( 4 )——④号单元中第i(i=3,4,5)节点所受力; K (—4)—④号单元的扩大刚度矩阵。
0
E t ,求整体刚度矩阵。
y
1
a
a
2①
3
③
②
④
x
4 a
5
6
a
7
5.1 整体刚度矩阵
2021/11/24
解:该模型中共有6个节点,4个单元,各单元的信息如表所示。
单元编号
各单元信息
①
②
a
y
1
2①
3
③
②
④
x
4 a
5
6
a
a
③
④
整体编码 局部编码
1、2、3 i、j、m
2、4、5 i、j、m
5、3、2 i、j、m
4
5.1 整体刚度矩阵
2021/11/24
不失一般性,仅考虑计算模型中有4个单元,如图所示。四个单元的整体节点位移列阵为
1T
T 2
T 3
T 4
T T 5
式中:
T i
, ui
vi
(i 1, 2,
, 5)
y
4
④
5
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
7.4 单元刚度矩阵组装及整体分析7.4.1 单刚组装形成总刚根据全结构的平衡方程可知,总体刚度矩阵是由单元刚度矩阵集合而成的.如果一个结构的计算模型分成个单元,那么总体刚度矩阵可由各个单元的刚度矩阵组装而成,即[K]是由每个单元的刚度矩阵的每个系数按其脚标编号“对号入座”叠加而成的.这种叠加要求在同一总体坐标系下进行.如果各单元的刚度矩阵是在单元局部坐标下建立的,就必须要把它们转换到统一的结构(总体)坐标系.将总体坐标轴分别用表示,对某单元有式中,和分别是局部坐标系和总体坐标系下的单元结点位移向量;[T]为坐标转换阵,仅与两个坐标系的夹角有关,这样就有是该单元在总体坐标系下的单元刚度矩阵.以后如不特别强调,总体坐标系下的各种物理参数均不加顶上的横杠.下面就通过简单的例子来说明如何形成总体刚度矩阵.设有一个简单的平面结构,选取6个结点,划分为4个单元.单元及结点编号如图3-27所示.每个结点有两个自由度.总体刚度矩阵的组装过程可分为下面几步:图7-27(1)按单元局部编号顺序形成单元刚度矩阵.图7-27中所示的单元③,结点的局部编号顺序为.形成的单元刚度矩阵以子矩阵的形式给出是(2)将单元结点的局部编号换成总体编号,相应的把单元刚度矩阵中的子矩阵的下标也换成总体编号.对下图3-27所示单元③的刚度矩阵转换成总体编号后为(3)将转换后的单元刚度矩阵的各子矩阵,投放到总体刚度矩阵的对应位置上.单元③的各子矩阵投放后情况如下:(4)将所有的单元都执行上述的1,2,3步,便可得到总体刚度矩阵,如式(3-9).其中右上角的上标表示第单元所累加上的子矩阵.(3-9)(5)从式(3-9)可看出,总体刚度矩阵中的子矩阵AB是单元刚度矩阵的子矩阵转换成总体编号后具有相同的下标,的那些子矩阵的累加.总体刚度矩阵第行的非零子矩阵是由与结点相联系的那些单元的子矩阵向这行投放所构成的.7.4.2 结点平衡方程我们首先用结构力学方法建立结点平衡方程.连续介质用有限元法离散以后,取出其中任意一个结点,从环绕点各单元移置而来的结点载荷为式中表示对环绕结点的所有单元求和,环绕结点的各单元施加于结点的结点力为.因此,结点的平衡方程可表示为(3-10)以[K]代入平衡方程,得到以结点位移表示的结点的平衡方程,对于每个结点,都可列出平衡方程,于是得到整个结构的平衡方程组如下:式中,[K]为整体刚度矩阵,为全部结点位移组成的向量,为全部结点载荷组成的向量.当然,如果各点的载荷向量也是在单元局部坐标下建立的,在合成以前,也应把它们转换到统一的结构(总体)坐标系下,即式中,是总体坐标系下的结点载荷向量,为坐标转换阵.7.4.3 位移边界条件在有限元法对结构进行整体分析时,建立了整体刚度矩阵[K],也得到了结构的刚度平衡方程,即.结构刚度方程的求解相当于总刚[K]求逆的过程.但是,从数学上看,未经处理的总刚是对称、半正定的奇异矩阵,它的行列式值为零,不能立即求逆.从物理意义看,在进行整体分析时,结构是处于自由状态,在结点载荷的作用下,结构可以产生任意的刚体位移.所以,在已知结点载荷的条件下,仍不能通过平衡方程惟一地解出结点位移.为了使问题可解,必须对结构加以足够的位移约束,也就是应用位移边界条件.首先要通过施加适当的约束,消除结构的钢体位移,再根据问题要求设定其他已知位移.所以,处理位移边界条件在有限元分析步骤中十分重要.约束的种类包括使某些自由度上位移为零,,或给定其位移值,还有给定支承刚度等,本书涉及前两种.处理约束的方法,常用的有删行删列法、分块法、置大数法和置“1”法等,下面分别予以介绍.1、删行删列法若结构的某些结点位移值为零时(即与刚性支座连接点的位移),则可将总体刚度矩阵中相应的行列、删行删列划掉,然后将矩阵压缩即可求解.这种方法的优点是道理简单.如果删去的行列很多,则总体刚度矩阵的阶数可大大缩小.通常用人工计算时常采用该方法.若用计算机算题,在程序编制上必带来麻烦,因为刚度矩阵压缩以后,刚度矩阵中各元素的下标必全改变.因而一般计算机算题不太采用.2.分块法为了理解这个方法,我们把方程分块如下:(3-11)其中,假设是给定的结点位移;是无约束的(自由)结点位移.因而是已知的结点力;是未知的结点力.方程(3-11)可以写为即(3-12)和(3-13)其中,不是奇异的,因而可以解方程(3-12)得出(3-14)一旦知道了,就可以由方程(3-13)求得未知结点力.在全部给定的结点自由度都等于零的特殊情况下,我们可以删除对应于的各行和各列(即删行删列法),故可把方程简写为(3-15)3.置“1”法由于全部给定的结点位移通常都不能在位移向量的开始或终了,故分块法的编号方法是很麻烦的.因此,为了引入给定的边界条件,可以采用下述等价的方法.可以把方程(3-12)和(3-13)合在一起写为(3-16)在实际计算中,方程(3-16)所示的过程可以在不重新排列所述方程的情况下用下述分块的方法为进行.步骤(1)如果把给定为,则载荷向量P可以修改为为结点自由度总数.步骤(2)除对角线元素以外,使[K]中对应于的行和列为零,而对角线元素为1,即步骤(3)在载荷向量中引入规定的值,即对全部规定的结点位移均应反复运用上述过程(步骤(1)到(3)).应当指出,由于这个过程保持了方程的对称性,因此,[K]可以按带状存储,而且几乎不会增加编制程序的工作量.4.置大数法置大数法的思路是:在总体刚度矩阵中,把指定位移所对应的行和列的对角元素乘上一个很大的数,如,此行其他元素保持不变,同时把该行对应的载荷项也相应地用来代替,这里为指定位移,于是原平衡方程组变为除第行外,其他各行仍保持原来的平衡特性,而第个方程式展开为由于上式中的比其他项的系数大得多,求和后可略去其小量,则上式变为即.这样就用近似方程组代替原方程组,得到近似满足边界条件的解.当指定位移为零时,只要将对角元素乘上一个大数,而相应的载荷项经证明可以不置零.删行删列法适用于指定零位移点,而置大数法适用于给定位移(包括零位移).5.斜支座的处理对于简单的约束情况(如限定某些结点位移为零或取得给定数值),可以用前述置大数法处理.有的结构在直角坐标系内建立了位移方程组,但在某个斜边上受有法向约束.如图3-28所示正方形固支板,受均布横向载荷,对此,可利用对称性而只计算其1/8,如图中ABC部分,其中AC为固支边,按对称性,AB边上有,但在BC边上应限定绕BC的转用等于零.为处理此类斜边上的约束,须对斜边上的结点做坐标变换.若结构的总体坐标系为为斜支座的局部坐标系(见图3-29).对于边界结点,须限定方向位移,为此,将边界结点的位移及载荷都变换到局部坐标轴系.设轴与斜支座的轴夹角为,逆时针为正,图7-28 图7-29 则依据第二单中坐标转换关系有其中,.或写成(3-17)与位移关系相同有(3-18)将上两式带入结构刚度方程有(3-19)这样把位移到列阵中凡是斜支座的结点位移矢量都用局部坐标表示了.将式(3-19)中第行左右两边前乘以(3-20)由上式可见:凡是边界点的斜支座,在刚度方程中对应于斜支座的位移和载荷向量均可直接斜支座的局部坐标值,总刚度距阵中的相应行列需作相应的变换.上式的系数矩阵仍然是对称的,而且此方程中结点位沿轴表示,这样,限定方向的位移就很方便了.实际计算中,并不需要建立结构总的位移方程组后再进坐标变换.而可以在形成单元刚度矩阵和结点载荷之后,就对斜支座点进行坐标变换,把变换后的单元刚度矩阵和结点载荷叠加入总刚度矩阵和总载荷的相应位置,最后叠加形成的也就是方程组(3-20),即需要处理的结点,应该在单元计算中完成坐标变换后再叠加,当结构有不同的斜边约束时,都可以这样处理,只不过对不同边上的结点,应按不同的方向余弦矩阵变换就是了.7.4.4 总刚度平衡方程的求解应用有限元法,最终都是归结为解总体刚度平衡方程,它实际上是以总体刚度矩阵为系数矩阵的大型线性代数方程组.通过对结构施加位移边界条件,消除了结构的刚体位移,从而消除总体刚度矩阵的奇异性,解这个线性代数方程组可求出结位移.我们已知,总体刚度矩阵具有大型、对称、稀疏、带状分布、正定、主元占优势的特点,稀疏表示刚度矩阵含有大量的零元素,带状表示非零元素集中在主对角线两侧.求解方程组应抓住上述特点,才能提高效率.首先,要为总刚度矩阵选择适当的存储方式,常用的有:(1)整体存储总刚.总刚的全部元素以二维数组形式放在计算机内存中,存储效率最低,适用于小型问题的分析.(2)等带宽二维存储总刚.总刚的下三角或上三角的带内元素取最大半带宽以二维数组形式存放在计算机内存中.其行数同整体总刚,列数等于最大半带宽.(3)一维变带宽存储总刚.将总刚下三角实际半带宽内元素逐行存放在一个一维数组内.其次,根据总刚度矩阵的存储方式,选择适当的求解刚度方程的方法:线性代数方程组的解法分直接解法(如高斯消去法、三角分解法)和迭代解法(如高斯-赛德尔迭代、超松弛迭代),本书主要讨论几种常用的直解法.(1)高斯消去法适合整体存储总刚.由于需要集合完整总刚,内存利用和计算效率都比较低.但高斯消去法原理和程序简单,作为初学便于理解.(2)对称消元法利用刚度矩阵对称,有每次消元的子阵均对称的性质,对高斯消去法稍加改进形成.这样就是只需组装总刚的上三角或下三角部分(3)带消元法将对称消元法进一步改造,使之适合总刚的等带宽二维存储.(4)因子化法(三角分解)又称Cholesky分解,适合一维变带宽存储总刚.这上方法储效率高,计算速度快,应用较为普遍.此外,还有一种方法,叫做波前法.波前法实际上也是一种改进的高斯消去法.它建立一个称为“波前”的空间,各单元刚度系数依次进入波前.一旦与某自由度有关的所有单元的刚度系数全部装入,便可将相应的变量消去.经过消元的方程的系数随即退出波前,存放在计算机的外存中.这样就可腾出空间装入新的刚度系数.所以,波前法不需要生成完整的总刚,而是边组装边消元,“成熟”一个消去一个.消元完成后,全部系数都已存储在计算机的外存或缓冲区中.回代时将各方程的系数按“先出后入”的顺序调入内存求解.由此可见,这种方法是利用计算机充裕的外存资源,以多耗取机时来缓解内存不足的矛盾,以便适应较大规模的问题.随着计算机技术的发展,内存资源不断扩大,对具有稀疏、带状性质的有限元刚度方程,这种以时间换取空间的办法得不偿失.另一方面,波前法的阐述和程序设计比较复杂,且对多种单元并存的结构使用不便.所以,本书不拟介绍波前法.本书第九章将详细讨论适合整体存储总刚的高斯消去法和适合一维变带宽存储的因子化法以及有关的程序设计问题,以下仅列出这两种方法的梗概.1、高斯消去法高斯循序消去法的一般公式:对于n阶线性代数方程,需进行次消元.采用循序消去时,第m次消元以m-1次消元后的m行元素作为主元行,为主元,对第行元素()的消元公式为(3-21)式中等的上角码(m),表示该元素是经过第m次消元后得到的结果.同样,可以把经过m次消元后的系数矩阵和载荷阵分别记为及.式表时第m 次消元是在经m-1次消元的基础上进行的.消元过程中,主元及被消元素的位置可见图3-30(a).图中阴影部分已完成消元过程的元素,主元行以下的矩阵为待消部分.在进行第m次时,1-m行元素的消元过程已经完成,其中的元素就是消元最后得到的上三角阵中的元素. m行发下的元素消元过程尚未结束,连同m行元素在内构成一个待消的方阵.消元共需进行n-1次.消元完成后,即可回代求解.我们把消元最后结果记为,为上三角阵,回代公式可写作(3-22)回代过程自后向前进行.当回代求解时,已经解得.回代示意图见图3-30(b),阴影部分为已求得解答的部分.图7-30 高斯消去法2.三角分解法总体刚度平衡方程中,[K]是对称、正定矩阵,因而可做如下分解(3-23)其中,则是单位上三角矩阵,.代入整本结构平衡方程记,则.即由向下回代.由其中第一个方程解得,再由第二个方程解得,……,依此类推可求得{Y}.又由向上回代,可得,由得依此类推可求得.由上述过程可见,三角分解法求解线性代数方程组的关键是对系数矩阵进行三角分解.7.4.5 求解内力由平衡方程组解出位移后,从中分离出各单元的结点位移,再通过方程(3-3)、(3-4)和(3-6)等计算各单元的应变、应力和结点力等内力。