第三节刚度矩阵
第三节 刚度矩阵
——节点载荷与节点位移之间的关系
一、 单元刚度矩阵
1. 单元刚度矩阵
xj
单元e 是在节点力作用下处于平衡。节点i 的节点力为
{}T
i xi
yi R R R ??=?? （i , j , m 轮换）
则单元e 的节点力列阵为
{}
T
e
T
T T m
i j
T
xm ym xi yi xj yj R R R R R R R R R R ???
?
????
=
=
单元应力列阵为
{}
T
e
x y xy σσστ????
=
假定弹性体的所有节点都产生一虚位移，单元e 的三个节点的虚位移为
{}
*
*****
*e
T m
m i i j j
u v u v u v δ???
?
= 单元虚应变列阵为
{}
****T
x y xy εεεγ????
??
=
参照式（3-7），则单元虚应变为
{}
{}*
*
e
e
B εδ
????=
作用在弹性体上的外力在虚位移上所做的功为：
{}{}*
e
T
e R δ?
? ?
?
?
单元内的应力在虚应变上所做的功为：
{}{}*T
e tdxdy εσ?
?? ??
?
??
根据虚位移原理，可得单元的虚功方程
{}
{}{}
{}**e
T
T
e e R tdxdy δεσ?
????
? ???
??
=
??
或
{}{}
{}{}*
*
e
T
T
T
e e B R tdxdy δδσ?
?
????? ?
?
??
?
?
?
?
=??
故有
{}
{}e
T
B R tdxdy σ?
????
=
??
将式（3-10）代入，的
{}
{}{}e
e
e
T
T
D B D B R B B tdxdy
tdxdy δδ??
???
????????????
?????????=
=
????
（3-27）
简记为
{}{}e
e e
k R δ????
= （3-29）
--------上式表征单元节点力与节点位移之间的关系，称为单元刚度方程（单元平衡方程） 其中
T
e
D B B k tdxdy ?
?????
???????????
=
??
（3-28） e
k ????称之为单元刚度矩阵（简称为单刚）
，是66?矩阵。 如果单元的材料是均质的，矩阵D ????中的元素也是常量，且在三角形常应变的情况下，矩阵B ????中的元素也是常
数，当单元的厚度也是常数时，注意到
dxdy ?
=???
，于
是单元刚度矩阵可简化为
T
e
B D B t k ??????????
????
????
?= （3-30）
将单元刚度矩阵按节点号写成分块矩阵形式：
66
e
ii ij im
ji jj jm mm mi
mj k
k k k k k k k k k ???
??
????????
????
????
= （3-31）
其中任一子块[]rs k (r ，s=i ，j ，m ）是一个2×2子矩阵，
为
[][][][]T
r
s
rs k B D B t =? (r ，s=i ，j ，m ）
（1）对于平面应力问题
将[]B 和平面应力问题的弹性矩阵[]D 代入，得
T
rs r s k B D B t ???
?????????????=? ()
2
112
2
114122r 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 c
c c b b μμμμμμμ--?
?
++??
=??
---???
++???
?
(r ，s=i ，j ，m ） （3-32）
（2）对于平面应变问题
将[]B 和平面应变问题的弹性矩阵[]D 代入，得
()()()()()()()
1212211211121212121e 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 μ分别换成2
1E μ
- 和
1μ
μ
-）
2. 单元刚度矩阵的性质 （1）
e
k ??
??
的物理意义
式（3-29）可完整写为
13141516
111221222324252633343536
3132434445464142555152535456616263646i i j
j
m
m
e
U 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
??
????
????????????
??????????????????????????????
??????????????
??
=??????????????????
5
66i
i j j m m e
u v u v u k v ??
?????
?
???????
???
??????????
????????????????????
??
可见每个节点在x 和y 方向上有二个平衡方程，3个节点共有六个平衡方程。
单元刚度矩阵[]e
k 中的任一元素称为刚度系数，其物理意义为：
ij k -----当单元的第j 个节点有单位位移，而其它节点位移为
零时，需在单元第i 个节点位移方向上施加的节点力的大小。 例如， 23k 表示是第3个节点有水平（x ）方向单位位移（即
31u =）时，而其它节点位移分量均为零时，在第2个节点
所引起的铅垂（y ）方向的节点力。
（2）单元刚度矩阵只取决于单元形状，大小，方向和弹性常数，而与单元的位置无关。即e
k ????不随单元坐标平移而改变，这叫单元刚度的平移原理。
23
5
例如图示结构，有 []
[]
(1)
(3)
k k =
另外，可以证明 []
[]
(1)
(2)
B B
=-
则有 []
[]
(1)
(2)
k k =
即单元旋转180?
后，单元刚度矩阵相等。这是单元刚度旋转原理。 （3） 单元刚度矩阵是对称矩阵。
因为 []
[][]T
e
k B D B t ????=?
所以有
[]()[][][]()
T
T
e
T
t k B D B =
?
[]
[
][]T
T
T
T B D B t =?????[][][]T B D B t =?e k ????=
（4）单元刚度矩阵是奇异矩阵。
即 0e
k ????
=
因为 {}[]
{}
e
e e R k δ=
当节点位移已知时，节点力是唯一确定的，而{}e
R 已知时，
{}e
δ不能唯一确定，因为单元没有支承，可以产生任意的刚
体位移。
根据上述性质：对于上图结构，在节点力为零时，单元仍可产生刚体位移，即
0111213141516U k u k v k u k v k u k v i i i j j m m
==+++++此时 0u u u
u i j m
=== ，0
v v v v i j m === ，单元产生
刚体位移0u ，0v 为任意的。故有
()()011131501214160k k k u k k k v +++++=
由于,0
u v 的任意性，则
0111315k k k ++= ， 0121416
k k k ++=
从而得
0111315121416
k k k k k k +++++= 同理可得：单元刚度矩阵任意一行（列）元素之和为零。 （5）单元刚度矩阵主对角线上的元素恒为正。
即 0(1,2,,6
ii k i >=
二、 整体分析
假设弹性体被分成m 个单元和n 个节点，对每一个单元进行前面的运算，则得到m 组型如
{}{}
e e
e k R δ????
=
的方程。把这些方程集合起来，便可得到表征整体弹性体平衡的刚度方程：
{}{}
2121
22n n n n k R δ?????????
= （3-37） 式中 1. {}21n δ
?-------整体结构的节点位移列阵，是由各节点
位移按节点号码从小到大顺序排列组成的，即
{}1
2
T T T T n δδδ
δ
???
?=
其中 {}T
i i
i u v δ???
?
= （i=1, 2，…，n ）
P 3
P
例如图示结构有
{}123411223344T
T T T T T
u v u v u v u v δδδδδ????
?
?
????
==
2.
{}21n R ?-------整体结构的节点载荷列阵， 是由各
节点载荷按节点号码从小到大顺序排列组成的，即
{}
1
2
T T T
T n R R R R ????=
其中
{}T
i i x i y
R R R ?????
?
= （i=1, 2，…，n ） 例如图示结构有
{}114431
2
4032T x y x y T
T T T P P R R P R R R R R R R ????
??
=??
??
=
--
2. 22n n
k ???????
-------整体结构的刚度矩阵（总刚）
（1）22n n
k ?????
的组集（“对号入座”法）
22m
e n n
e
k k ???????????
=∑
例图示结构有 单元1
22
2421(1)424441121411ii ij im
ji jj jm mm mi mj k k k k k k k k
k k k k k k k k k k k ???
?
????
??????????
??
????
?????
?
==
单元2
444243(2)24
2223343233ii ij im
ji jj jm mm mi mj k k k k k k k k
k k k k k k k k k k k ???
?
????
??????????
??
????
???
?
?
?==
注：在单元刚度矩阵中，各子块的下标表示该子块在总刚中
的位置。
则总刚为
(1)(1)(1)
111214(1)(1)(2)(2)(1)(2)21
2222232424(2)(2)(2)
323334(1)(1)(2)(2)(1)(2)41
424243
44448800
k k k k k k k k k k k k k k k k k k k ????????
??????????
??
????
?++=++
（2）总刚的性质 ⅰ. 整体刚度矩阵的物理意义
[]K 中每一列元素的物理意义为：
欲使弹性体的某一个节点在坐标轴方向发生单位位移，而其它节点都保持为零状态，在各节点上所需要加施加的节点力。
由式可以看出，令节点1在坐标轴x 方向有单位位移，即11u =，而其余的节点位移为零时，即v 1=u 2=v 2=u 3=v 3=······= u 2n =v 2n =0，这样就可得到节点载荷列阵等于总刚[]k 的第一列元素组成的列阵，即
…………112211213141(21)1(2)1T
nx
ny x y x y T
n n R R R R R R K K K K K K ?
?
??
???
?
-??
??
=
ⅱ. 总刚k ????是对称矩阵 ⅲ. 总刚k ????是奇异矩阵
ⅳ. 总刚k ????
主对角线上的元素恒为正，即 0(1,2,,2
)ii k i n >= ⅴ. 总刚[]k 是一个稀疏矩阵。
若遵守一定的节点编号规则，则非零元素集中在主对角线附近呈带状分布。 单元越多，总刚[]k 越稀疏。
0,0/rs r s r s
k ???
??????
??????
????
≠==同属于一个单元的两个节点
号码
非零元素集中在主对角线两侧，在包括对角线元素在内的半个带形区域中，具有最多元素的数目称为最大半带宽（半带宽），用B 表示：
B=（max{单元节点号码的最大差值}+1）×节点自由度数
半带宽取决于节点号码的最大差值。半带宽越窄，计算机的存储量就越少。所以，在划分有限元网格进行节点编号时，要采用合理的编码方式，使同一单元中两节点的号码差尽可能地小，以便使半带宽小，节省存储空间，提高计算效率。而且还可以大幅度减少求解方程所需的运算次数，其效果对大型结构显得尤为突出。
10
1234
56
7
8
9
1234567
8910
（a ） （b ）
对于图（a ） B=[(7-1)+1]×2=14 对于图（b ） B=[(4-1)+1]×2=8
ⅳ. 总刚[]k 的存储方式
通常的有限元程序，一般都利用总刚的对称性和稀疏性的特点，在计算时采用：
半带宽存储-------------只存储上半带的元素 一维变带宽存储 分块一维变带宽存储
常用单元的刚度矩阵
r u r r u r =-+= πππεθ22)(2 由于各点在圆周方向上无位移，因而剪应变θr v 和r v θ均为 零。将应变写成向量的形式，则{}?? ?? ? ?????? ?????? ???????+??????=??????????????=r w z u z w r u r u rz z r γεεεεθ 根据上式，可推导出几何方程{}[]{})(e B ?ε= 其中几何矩阵[]????????? ?????????? ??= ij ji ki ik jk kj ji ik kj k j i ij kj jk z r z r z r r r r r z r N r z r N r z r N z z z B 000 0),(0),(0),(00021 3.弹性方程和弹性矩阵[D] 依照广义虎克定律，同样可以写出在轴对称中应力和应变之间的弹性方程，其形式为 [])(1 θσσσε+-= z r r u E [])(1 z r u E σσσεθθ+-= [])(1 θσσσε+-=r z z u E rz rz E r τμ)1(2+= 所以弹性方程为{}[]{}εσD = 式中应力矩阵{}{}T rz z r τσσσσθ=
弹性矩阵[]? ? ??????? ???? ?-----+=221000010101)21)(1(μμμμμμμμμμ μμE D 4.单元刚度矩阵[])(e k 与平面问题相同，仍用虚功原理来建立单元刚度矩阵，其积分式为 [][][][]dV B D B k V T e ?=)( 在柱面坐标系中，drdz dV π2= 将drdz dV π2=代入[][][][]dV B D B k V T e ?=)(，则[][][][]rdrdz B D B k T e ??=π2)( 即为轴对称问题求单元刚度矩阵的积分式。 与弹性力学平面问题的三角形单元不同，在轴对称问题中，几何矩阵[B]有的元素（如r z r N i ),(等）是坐标r 、z 的函 数，不是常量。因此，乘积[][][]B D B T 不能简单地从式 [][][][]rdrdz B D B k T e ??=π2)(的积分号中提出。如果对该乘积逐项求 积分，将是一个繁重的工作。一般采用近似的方法：用三角形形心的坐标值代替几何矩阵[B]的r 和z 的值。用[]B 表示在形心),(z r 处计算出的矩阵[B]。其中 3 ) (,3 ) (k j i k j i z z z z r r r r ++= ++= 只要单元尺寸不太大，经过这样处理引起的误差也不大。被积函数又成为常数，可以提出到积分号外面：
《结构力学习题集及答案》(下)-1a
第八章 矩阵位移法 一、判断题： 1、单元刚度矩阵反映了该单元杆端位移与杆端力之间的关系。 2、单元刚度矩阵均具有对称性和奇异性。 3、局部坐标系与整体坐标系之间的坐标变换矩阵T 是正交矩阵。 4、结构刚度矩阵反映了结构结点位移与荷载之间的关系。 5、结构刚度方程矩阵形式为：[]{}{}K P ?=，它是整个结构所应满足的变形条件。 6、图示结构用矩阵位移法计算时（计轴向变形）未知量数目为8个。 7、在直接刚度法的先处理法中，定位向量的物理意义是变形连续条件和位移边界条件。 8、等效结点荷载数值等于汇交于该结点所有固端力的代数和。 9、矩阵位移法中，等效结点荷载的“等效原则”是指与非结点荷载的结点位移相等。 10、矩阵位移法既能计算超静定结构，也能计算静定结构。 11、已知图示刚架各杆EI = 常数，当只考虑弯曲变形，且各杆单元类型相同时，采用先处理法进行结点位移编号，其正确编号是： (0,1,2) (0,0,0) (0,0,0) (0,1,3) (0,0,0)(1,2,0) (0,0,0)(0,0,3) (1,0,2) (0,0,0) (0,0,0)(1,0,3) (0,0,0) (0,1,2) (0,0,0)(0,3,4) A. B. C. D. 2134123412341234 ( )
二、计算题： 12、用先处理法计算图示结构刚度矩阵的元素133322,,K K K 。 12 3l l 4 l 5EI 2EI EA (0,0,0) (0,0,1) (0,2,3) (0,0,0) (0,2,4)(0,0,0) EI 13、用先处理法计算图示刚架结构刚度矩阵的元素153422,,K K K 。EI ，EA 均为常数。 l 14、计算图示结构整体刚度矩阵的元素665544,,K K K 。E 为常数。 l l 1 3 4 2 A , I A A /222A I , 2A 15、写出图示结构以子矩阵形式表达的结构原始刚度矩阵的子矩阵[][]K K 2224,。 [][]k k 1112 [][] k k 2122 [] k = i i i i i 单刚分块形式为 ： 16、已知平面桁架单元在整体坐标系中的单元刚度矩阵，计算图示桁架结构原始刚度矩阵[]K 中的元素,,7877K K EA ＝常数。,cos α=C ,sin α=S ,C C A ?=
一般单元在局部坐标系下的单元刚度矩阵
9.3 一般单元在局部坐标系下的单元刚度矩阵 1.杆端内力与位移关系回顾 （轴向）； ；（弯曲）； 2.公式推导（图1） 图1 杆件性质：长度l，截面面积A，截面惯性矩I，弹性模量E；杆端位移u、v、θ。 （1） （2）列成矩阵形式：
（3） 即：（4）局部坐标系下单元刚度矩阵： （5） 9.4 梁单元 1.简支梁 简支梁单元见图1。 图1 说明：（a）梁单元通常忽略轴向变形；（b）图10-3中；相应的力分量也应该为零；（c）依据刚度矩阵的物理意义，可以由一般单元的刚度矩阵生成梁单元矩阵。即去掉位移分量为零 的相应行和列。
即：单元刚度方程：单元刚度矩阵： （1） 2.悬臂梁等 思考：建立图2的单元刚度矩阵：（固定端位移为零；自由端有转角和竖向位移) 图2 图a：图b： 3.桁架 仅有轴向位移 9.5 单元刚度系数的物理意义 1.单元刚度系数的意义 一般地，第j 个杆端位移分量取单位值1，其它杆端位移为0 时所引起的第i个杆端力分量的值。
例：的物理意义：当第3个杆端位移分量时引起的第5个杆端力分量。 对称性 （反力互等定理） 3.奇异性（，不存在逆矩阵） 根据式可由杆端位移求解杆端力，且是唯一解。但由杆端力求杆端位移，可能无解，如有解也是非唯一解。 说明：已知6个杆端力分量，（a）无法保证力状态的合法性——可能造成无解；（b）无法确定杆的支承条件——可能造成非唯一解。 9.6 单元坐标转换矩阵的物理意义 1.问题的提出 单元刚度矩阵——单根杆；多根根组成的复杂结构呢？（图1）
图1 分析（a）从数学的角度理解整体坐标系（xy）与局部坐标系（）的区别； （b）力分量应向整体坐标系转换，图f给出了两种坐标系下力分量之间的数学关系： 。 同理： 2.公式推导 矩阵形式： （1）同理：（2）
7.4 单元刚度矩阵组装及整体分析
根据全结构的平衡方程可知，总体刚度矩阵是由单元刚度矩阵集合而成的 个结构的计算模型分成个单元，那么总体刚度矩阵可由各个单元的刚度矩阵组装而成，即是由每个单元的刚度矩阵的每个系数按其脚标编号“对号入座”叠加而成的 将总体坐标轴分别用表示，对某单元有 式中，和分别是局部坐标系和总体坐标系下的单元结点位移向量 是该单元在总体坐标系下的单元刚度矩阵
. 将单元结点的局部编号换成总体编号，
其中右上角的上标表示第单元所累加上的子矩阵 具有相同的下标，的那些子矩阵的累加总体刚度矩阵第行的非零子矩阵是由与结点相联系的那，从环绕点各单元移置而来的结点载荷为 式中表示对环绕结点的所有单元求和，环绕结点的各单元施加于结点的结点力为
.因此，结点的平衡方程可表示为 得到以结点位移表示的结点的平衡方程， 为整体刚度矩阵，为全部结点位移组成的向量，为全部结点载荷组成的向量式中，是总体坐标系下的结点载荷向量，为坐标转换阵 . 构是处于自由状态，在结点载荷的作用下，结构可以产生任意的刚体位移 的条件下，仍不能通过平衡方程惟一地解出结点位移.
约束的种类包括使某些自由度上位移为零，，或给定其位移值，还有给定支承刚 为了理解这个方法，我们把方程分块如下： 其中，假设是给定的结点位移；是无约束的（自由）结点位移因而是已知的结点力；
其中，不是奇异的，因而可以解方程（ 一旦知道了，求得未知结点力. 殊情况下，我们可以删除对应于的各行和各列（即删行删列法），故可把方程简写为 由于全部给定的结点位移通常都不能在位移向量的开始或终了，故分块法的编号方法是很麻烦因此，为了引入给定的边界条件，可以采用下述等价的方法 如果把给定为，则载荷向量 为结点自由度总数
最新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 结点平衡方程 我们首先用结构力学方法建立结点平衡方程.连续介质用有限元法离散以后，取出其中任意一个结点，从环绕点各单元移置而来的结点载荷为 式中表示对环绕结点的所有单元求和，环绕结点的各单元施加于结点的结点力为
1试用先处理法建立图示刚架的结构刚度矩阵设各杆EI为常
1．试用先处理法建立图示刚架的结构刚度矩阵[] K 。设各杆EI 为常数，长度为l 。(只考虑弯曲变形) 3 1 24 56① ② ③ ④ ⑤ 附： EA l EA l EI l EI l EI l EI l EI l EI l EI l EI l EA l EA l EI l EI l EI l EI l EI l EI l EI l EI l 00000126012606406200000126012606206432322232322--------?????????????????????????????? 2．求桁架单元⑨的杆端力{}F ⑨ ，已知{}[]?=u v u v u v u v u 223344556 T []=----(/).........Pl EA 133337684026667702782500076840116677027838333 T 。 l l l 3 ．已知桁架结点位移列阵(结构坐标系)为 [][]u v u v u v Pl EA 112233********* 0 0 T T =---+ /(.) 。 试求单元①的杆端力列阵。(局部坐标系) l P
4．图示结构，求等效结点荷载列阵{} P E 。(力矩以顺时针为正) l /2 l 0/2l 5．图a 所示结构，整体坐标见图b ，图中圆括号内数码为结点定位向量(力和位移均按水平、竖直、转动方 向顺序排列)，不考虑轴向变形。求等效结点荷载列阵{} P E 。 3m 36(b)(a) 附： EA l EA l EI l EI l EI l EI l EI l EI l EI l EI l EA l EA l EI l EI l EI l EI l EI l EI l EI l EI l 00000126012606406200000126012606206432322232322------- -????????????????????????????? ? 6．求图示刚架结构刚度矩阵的任意两个主元素。已知各杆I A //=110。 4m
第三节刚度矩阵(汇编)
第三节 刚度矩阵 ——节点载荷与节点位移之间的关系 一、 单元刚度矩阵 1. 单元刚度矩阵 xj 单元e 是在节点力作用下处于平衡。节点i 的节点力为 {}T i xi yi R R R ??=?? （i , j , m 轮换） 则单元e 的节点力列阵为 {} T e T T T m i j T xm ym xi yi xj yj R R R R R R R R R R ??? ? ???? = = 单元应力列阵为 {} T e x y xy σσστ???? =
假定弹性体的所有节点都产生一虚位移，单元e 的三个节点的虚位移为 {} * ***** *e T m m i i j j u v u v u v δ??? ? = 单元虚应变列阵为 {} ****T x y xy εεεγ???? ?? = 参照式（3-7），则单元虚应变为 {} {}**e e B εδ????= 作用在弹性体上的外力在虚位移上所做的功为： {}{}* e T e R δ?? ?? ? 单元内的应力在虚应变上所做的功为： {}{}* T e tdxdy εσ? ?? ?? ? ?? 根据虚位移原理，可得单元的虚功方程 {} {}{} {}**e T T e e R tdxdy δεσ? ???? ? ??? ?? = ?? 或 {}{} {}{}* * e T T T e e B R tdxdy δδσ? ? ????? ? ??? ? ? ? ? =??
故有 {} {}e T B R tdxdy σ? ???? = ?? 将式（3-10）代入，的 {} {}{}e e e T T D B D B R B B tdxdy tdxdy δδ?? ??? ???????????? ?????????== ???? （3-27） 简记为 {}{}e e e k R δ?? ?? = （3-29） --------上式表征单元节点力与节点位移之间的关系，称为单元刚度方程（单元平衡方程） 其中 T e D B B k tdxdy ? ????? ??????????? = ?? （3-28） e k ????称之为单元刚度矩阵（简称为单刚） ，是66?矩阵。 如果单元的材料是均质的，矩阵D ????中的元素也是常量，且在三角形常应变的情况下，矩阵B ????中的元素也是常
结构力学11.4 结构的原始刚度矩阵
§11-4结构的原始刚度矩阵 整体分析：即建立求解基本未知量的结构刚度方程。而位移法中求解的基本未知量是结点位移， 包括线位移与角位移。 考虑如上图所示刚架，有4个结点，3个单元，受结点荷载的作用。至于非结点荷载作用的情况，需要将其转化为等效的结点荷载，这将在后面的内容中进行专门介绍。这里，暂只考虑结点荷载作用的情况。 各单元的局部坐标系与整体坐标系如下图所示。 各单元的单元刚度矩阵，进行坐标变换后，得到整体坐标系下各单元的单元刚度矩阵如下 1221][22211211??????=①①①①①k k k k k ，2332][33322322??????=②②②②②k k k k k ，4 3][44433433??????=③③③③③k k k k k 34每个结点，有x 方向线位移、y 方向线位移与角位移3个位移分量。结构的结点位移列向量为 ??????????=1111}{?v u Δ，??????????=2222}{?v u Δ，??????????=3333}{?v u Δ，?? ????????=4444}{?v u Δ??? ???????????=4321}{ΔΔΔΔΔ
与结点位移列向量对应的结点外力（包括荷载和反力）列向量为 结构刚度方的建立 前面学习位移法时，已经知道，位移法方程，即结构的刚度方程，就是结点的平衡方程。所以，通过结点的平衡条件，可建立结构的刚度方程。 下面，以结点2为例，如图示。结点2的3个平衡方程为 ②①222x x x F F F +=，②①222y y y F F F +=，②①2 22M M M += 写成矩阵形式，有 ?? ????????+??????????=??????????②②②①①①222222222M F F M F F M F F y x y x y x ，即，} {}{}{222②①F F F +=单元杆端力，可用杆端位移表示为 }]{[}]{[}{2221212①①①①①δδk k F +=，} ]{[}]{[}{3232222②②②②②δδk k F +=得到， } ]{[}]){[]([}]{[}{323222221212Δk Δk k Δk F ②②①①+++=此即结点2的平衡方程。 同理，可列出结点1、3、4的平衡方程。共有4个结点，结果为
求整体刚度矩阵matlab程序
班级：机自11-13 学号：03111323 求整体刚度矩阵。 部分分析思路： 各单元信息 单元编号 — ① ② ③ ④ 整体编码 1,2,3 2,4,5 5,3,2 3,5,6 局部编码 i, j,m i, j,m i, j,m i, j,m 以整编码体表示 的单元刚度矩阵 K ⑴ K ⑴ K ⑴ K 11 K 12 K 13 K ⑴ K ⑴ K ⑴ K 21 K 22 K 23 K ⑴ K ⑴ K ⑴ K 31 K 32 K 33 K ⑵ K ⑵ K ⑵ K 22 K 24 K 25 K ⑵ K ⑵ K ⑵ K 42 K 44 K 45 K ⑵ K ⑵ K ⑵ K 52 K 54 K 55 K ⑶ K ⑶ K ⑶ K 55 K 53 K 52 K ⑶ K ⑶ K ⑶ K 35 K 33 K 32 K ⑶ K ⑶ K ⑶ K 25 K 23 K 22 K ⑷ K ⑷ K ⑷ K 33 K 35 K 36 K ⑷ K ⑷ K ⑷ K 53 K 55 K 56 K ⑷ K ⑷ K ⑷ K 63 K 65 K 66 整体刚度矩阵 姓名：吴佳侣 例：如图所示有限元模型，弹性模量为 E =1，厚度为t =1，为简化计算取 ， K 12 K 13 K 14 K 15 K 16 K 22 K 23 K 24 K 25 K 26 K 32 K 33 K 34 K 35 K 36 K 42 K 43 K 44 K 45 K 46 K 52 K 53 K 54 K 55 K 56 K 62 K 63 K 64 K 65 K 66
aA2/ 2 程序： %%%%%%%%%%%%%%求%整体刚度矩阵 clc clear all syms E u t a%%定义变量 %%%%%%%%%%%%%%%%%%%%%%写%出%单元刚度矩阵 ,单元面积为 %%%%%%%%%%%单%元 1 bi=0;ci=a; bj=-a;cj=-a; bm=a;cm=0; mianji=a A 2/2; B1=1/2/mianji*[bi 0 bj 0 bm 0 0 ci 0 cj 0 cm ci bi cj bj cm bm]; u=0;E=1;t=1; D=E/(1-uA2)*[1 u 0 u 1 0 0 0 (1-u)/2]; %%%弹性矩阵 D k1=transpose(B1)*D*B1*t*mianji; %%%%%%%%%%%单%元 2 bi=0;ci=a; bj=-a;cj=a; bm=a;cm=0; mianji=aA2/2; B2=1/2/mianji*[bi 0 bj 0 bm 0 0 ci 0 cj 0 cm ci bi cj bj cm bm]; u=0;E=1;t=1; D=E/(1-uA2)*[1 u 0 u 1 0 0 0 (1-u)/2]; %%%弹性矩阵 D K3=transpose(B2)*D*B2*t*mianji; %%%%%%%%%%%单%元 3 bi=0;ci=-a; bj=a;cj=-a; bm=-a;cm=0; mianji=aA2/2; B3=1/2/mianji*[bi 0 bj 0 bm 0 0 ci 0 cj 0 cm ci bi cj bj cm bm];
求整体刚度矩阵matlab程序
姓名：吴佳侣 班级：机自11-13 学号：03111323 例： 如图所示有限元模型，弹性模量为E 1= ，厚度为t 1= ，为简化计算取=0μ ， 求整体刚度矩阵。 部分分析思路： 各单元信息 整体刚度矩阵 1112 1314151621222324252631 32333435364142434445465152535455566162 63 64 65 66?????????????????? ? ?K K K K K K K K K K K K K K K K K K K =K K K K K K K K K K K K K K K K K K
程序： %%%%%%%%%%%%%%%求整体刚度矩阵 clc clear all syms E u t a%%定义变量%%%%%%%%%%%%%%%%%%%%%%%%写出单元刚度矩阵,单元面积为a^2/2 %%%%%%%%%%%%单元1 bi=0;ci=a; bj=-a;cj=-a; bm=a;cm=0; mianji=a^2/2; B1=1/2/mianji*[bi 0 bj 0 bm 0 0 ci 0 cj 0 cm ci bi cj bj cm bm]; u=0;E=1;t=1; D=E/(1-u^2)*[1 u 0 u 1 0 0 0 (1-u)/2]; %%%弹性矩阵D k1=transpose(B1)*D*B1*t*mianji; %%%%%%%%%%%%单元2 bi=0;ci=a; bj=-a;cj=a; bm=a;cm=0; mianji=a^2/2; B2=1/2/mianji*[bi 0 bj 0 bm 0 0 ci 0 cj 0 cm ci bi cj bj cm bm]; u=0;E=1;t=1; D=E/(1-u^2)*[1 u 0 u 1 0 0 0 (1-u)/2]; %%%弹性矩阵D K3=transpose(B2)*D*B2*t*mianji; %%%%%%%%%%%%单元3 bi=0;ci=-a; bj=a;cj=-a; bm=-a;cm=0; mianji=a^2/2; B3=1/2/mianji*[bi 0 bj 0 bm 0 0 ci 0 cj 0 cm ci bi cj bj cm bm];
ansys单元刚度矩阵的提取
看了这么久了都没人回，查了一些质料终于找到答案了，， 下面提供三种方法：方便与其他程序进行接口编程1. Which matrix you would like? element stiffness matrix or full stiffness matrix? element stiffness is within file.emat. full stiffness matrix is within file.full A simple way to dump the matrix is as follow: ------------------- /aux2 fileaux2,file,emat form,long dump,all ------------------- 2. 可以使用/DEBUG命令来得到。详细步骤参见下面的宏文件 finish /clear PI=3.1415926 w1=3 w2=10 w3=6 w4=1.2 r=.8 t=0.08 /PREP7 !* ET,1,SHELL63 R,1,t ET,2,MASS21 R,2,500,500,500,2000,2000,2000, !* UIMP,1,EX, , ,2e11 UIMP,1,NUXY, , ,0.3, UIMP,1,DAMP, , ,0.2, UIMP,1,DENS, , ,7800,
BLC4,0,0,w2,w1 ESIZE,1.5,0, AMESH,all NSEL,S,LOC,X,0.0 D,all, , , , , ,ALL, , , , , allsel,all SFA,all,1,PRES,12 FINISH /OUTPUT,cp,out,, ! 将输出信息送到cp.out文件 /debug,-1,,,1 ! 指定输出单元矩阵 /SOLU SOLVE finish /OUTPUT, TERM ! 将输出信息送到output windows中 ! 这时用编辑器打开cp.out文件，可以看到按单元写出的质量、刚度等矩阵 3. 其原理很简单，即使用ansys的超单元即可解决问题。定义超单元，然后列出超单元的刚度矩阵即可。 面是一个小例题，自可明白。 /prep7 k,1 k,2,3000 l,1,2 et,1,beam3 mp,ex,1,2e5 mp,prxy,1,0.3 r,1,5000,2e7,200 lesize,all,,,10 lmesh,all finish !----以上正常建立模型，不必施加约束和荷载 /solu antype,7 !substructuring分析类型 seopt
结构刚度矩阵的特点
由前面的讨论可知结构的刚度矩阵K是由单元刚度矩阵集合而成，它与单元刚度矩阵类同也具有明显的物理意义。有限元的求解方程(32)式是结构离散后每个结点的平衡方程。结构刚度矩阵K的任一元素K ij的物理意义是：结构第j个结点位移为单位值而其它结点位移皆为零时，需在第i个结点位移方向上施加的结点力的大小。与单元不同之处在于结构是单元的集合体，每个单元都对结构起一定的作用。由于单元刚度矩阵是对称和奇异的，由它们集成的结构刚度矩阵K也是对称和奇异的，也就是说结构至少需给出能限制刚体位移的约束条件才能消除K的奇异性，以便由(32)式求得结点位移。 连续体离散为有限个单元体，由图1可见，每个结点的相关单元只是围绕在该结点周围为数甚少的几个，一个结点通过相关单元与之发生关系的相关结点也只是它周围的少数几个，因此虽然总体单元数和结点数很多，结构刚度矩阵的阶数很高，但刚度系数中非零系数却很少，这就是刚度矩阵的大型和稀疏性。只要结点编号是合理的，这些稀疏的非零元素将集中在以主对角线为中心的一条带状区域内，即具有带状分布的特点。如图7所示。 综上所述，有限单元法最后建立的方程组的大型系数矩阵K具有以下性质:(1)对称性（2)奇异性;(3)稀疏性；(4)非零元素呈带状分布。由于方程组的大型，在求解方程时，除引入位移边界条件使奇异性消失外，其他特点都必须在解方程中予以充分的考虑和利用，以提高解题的效率。" 七、实施步骤与注意事项 利用上面讨论的三角形常应变单元解平面问题，其具体步骤可归纳如下: 1)将要计算的弹性体划分成三角形单元。对结点进行编号，列出结点坐标作为输入信息。 (2)对单元进行编号，列出单元三个结点的号码作为输入信息。 (3)计算载荷的等效结点力，把等效结点力作为输入信息。 (4)按照(6)式计算各单元的常数b i、c i、b j、c j、b m、c m，再按照(4)计算2A。 (5)按照(35)式计算各单元的刚度矩阵。 (6)形成整体刚度矩阵。 (7)处理约束及消除刚体位移。 (8)解线性方程组(32)式，求结点位移。 (9)按照(20)式计算应力矩阵，再按(18)式计算单元应力。根据需要计算主应力和主方向。 通常步骤(4)至(9)均由计算机来完成，而步骤(1)至(3)可以用手工完成，也可由计算机来完成。在实现以上各步骤时，为了达到一定的计算精度，节约计算机存储量，缩短计算机运行时间等目的，还需要注意下列事项。 1、利用对称性 在划分单元前要研究一下，计算对象是否有对称变形或反对称变形存在，从而确定是否需要取整个物体，还是取部分物体作为计算模型。例如图8a所示受纯弯曲的梁，它对于x,y轴都对称，而载荷对于y轴对称，对于x轴反对称。可见，应力和应变亦将具有同样的对称和反对称特性，所以我们只需计算1/4梁就行了。分离体如图8b所示。对于删去部分结构的影响可以这样考虑：对于处于y轴对称面内各结点的x方向位移和y方向分布力都应等于零，而对于处在x轴反对称面上的各结点的x方向位移和y方向分布力亦都应等于零。这些条件相当于安置如图8b中的约束。图中o点上安置y方向的约束是为了消除刚体位移而设置的。又例如在分析图9中所承受均匀压力的厚壁圆筒时，根据结构和载荷轴对称的性质，我们可以取出一个小扇形(图中阴影部分)进行计算。扇形的两侧边上应加上约束，以消除周向位移和沿径向的分布力。(其中CD边界上的约束和坐标方向不一致，将给计算带来一定麻烦。此问题如采用下一节的轴对称有限元进行分析，将方便得多。)
单元刚度矩阵MATLAB编程
《有限元法》实验报告 专业班级力学（实验）1601 姓名田诗豪 学号 10 提交日期
实验一（30分） 一、实验内容 编写一个计算平面3结点三角形单元的应变矩阵、应力矩阵和单元刚度矩阵的MATLAB 函数文件[B3,S3,K3] = ele_mat_tri3(xy3,mat)，其中：输入变量xy3为结点坐标数组，mat 为材料参数矩阵；输出变量B3为应变矩阵，S3为应力矩阵，K3为单元刚度矩阵。（要求给出3个不同算例进行验证，并绘制出单元形状和结点号） 二、程序代码 通用函数 function [B3,S3,K3] = ele_mat_tri3(xy3,mat) %生成平面3结点三角形单元的应变矩阵、应力矩阵和单元刚度矩阵的功能函数 %*********变量说明**************** %xy3------------------结点坐标数组 %mat------------------材料参数矩阵（弹性模量，泊松比，壁厚） %B3-------------------应变矩阵 %S3-------------------应力矩阵 %K3-------------------单元刚度矩阵 %********************************* xyh=[1,xy3(1,1),xy3(1,2);1,xy3(2,1),xy3(2,2);1,xy3(3,1),xy3(3,2)]; A=*det(xyh); A=abs(A); D=mat(1)/(1-mat(2)^2)*[1,mat(2),0;mat(2),1,0;0,0,(1-mat(2))/2]; b=zeros(1,3);c=zeros(1,3); %********************************* for i=1:3 if i==1 j=2; m=3; elseif i==2 j=3; m=1; else j=1; m=2; end b(i)=xy3(j,2)-xy3(m,2); c(i)=xy3(m,1)-xy3(j,1); end
1试用先处理法建立图示刚架的结构刚度矩阵设各杆EI为常
1．试用先处理法建立图示刚架的结构刚度矩阵[]K 。设各杆EI 为常数，长度为l 。(只考虑弯曲变形) 3 1 24 56① ② ③ ④ ⑤ 附： EA l EA l EI l EI l EI l EI l EI l EI l EI l EI l EA l EA l EI l EI l EI l EI l EI l EI l EI l EI l 00000126012606406200000126012606206432322232322--------?????????????????????????????? 2．求桁架单元⑨的杆端力{}F ⑨ ，已知{}[]?=u v u v u v u v u 223344556 T []=----(/).........Pl EA 133337684026667702782500076840116677027838333 T 。 l l l 3 ．已知桁架结点位移列阵(结构坐标系)为 [][]u v u v u v Pl EA 112233********* 0 0 T T =---+/ (.) 。 试求单元①的杆端力列阵。(局部坐标系) l P
4．图示结构，求等效结点荷载列阵{}P E 。(力矩以顺时针为正) l /2 l 0/2l 5．图a 所示结构，整体坐标见图b ，图中圆括号内数码为结点定位向量(力和位移均按水平、竖直、转动方 向顺序排列)，不考虑轴向变形。求等效结点荷载列阵{}P E 。 3m 1436(b)(a) 附： EA l EA l EI l EI l EI l EI l EI l EI l EI l EI l EA l EA l EI l EI l EI l EI l EI l EI l EI l EI l 00000126012606406200000126012606206432322232322------- -????????????????????????????? ? 6．求图示刚架结构刚度矩阵的任意两个主元素。已知各杆I A //=110。 4m