当前位置：文档之家› 常用单元的刚度矩阵

常用单元的刚度矩阵

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

kj ji ik kj k j i ij

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

)

(k j i k j i z z z z r r r r ++=

++=

只要单元尺寸不太大，经过这样处理引起的误差也不大。被积函数又成为常数，可以提出到积分号外面：

[][][][][][][]?==??r B D B rdrdz B D B k T

e ππ22)(式中?——三角形的面积。

由式[]

[][][][][][]?==??r B D B rdrdz B D B

k T

e ππ22)

(可以看出，两轴对

称的三角形单元，当形状、大小及方位完全相同而位置不同时，其刚度矩阵也不相同。距离主轴线越远的单元，其刚度越大。这与平面问题不一样。

二、等参数的刚度矩阵

对一些由曲线轮廓的复杂结构，如果采用直角边单元进行离散，由于用直线代替了曲线，除非网格划分得很细，否则不能获得较高的精度；对另一些应力随坐标急剧变化的结构，采用简单的常应力单元离散时，也必须划分成大量的微小单元，以保证足够的精度。为此引入一种高精度的单元——等参数单元。它既能简化复杂单元划分的工作，又能在满足同样精度的要求时，大大减少使用的单元数。目前流行的大程序中较常用，它成功地解决了许多二维和三维的弹性力学问题。

为导出等参数单元的刚度矩阵，首先要建立根据每个单元的形状确定的自然坐标系，然后将位移模式和形状函数都写成自然坐标的函数。

一个单元在自然坐标系内的点余元整体坐标系内的点成一一对应的关系。通过映射，可以将整体坐标系中的图形转化为自然坐标系中的相应徒刑。例如可以将整体坐标系中的一个任意四边形（实际单元）映射到自然坐标系中成为一

个正方形（基本单元）。同样也可以将任意四面体、六面体（包括直边和曲边的）分别映射成正四面体和正六面体。

这里只介绍较简单的一种平面问题的情况，将整体坐标系中的一个任意四边形映射成自然坐标系中的一个正方体，并导出单元刚度矩阵。其它种单元的映射，可依次原理进行。不再叙述。

1. 位移模式和形状函数

图4-2中的任意四边形单元上，作连接对边中点的直线，取其交点为原点，这两条直线分别为ξ和η轴，并令四条边上的ξ和η值分别为1±，建立一个新的坐标系，称之为该单元的自然坐标系。原坐标系XOY 称为整体坐标系。在整体坐标系中，自然坐标系非正交，它由任意四边形的形状所确定。

图4-19

如果将自然坐标系改画成直角坐标系，那么图4-19（a ）中的任意四边形单元就成为图4-19（b ）所示的正方形。

上述两个四边形的点（包括顶点）一一对应，即它们之间相互映射。因此，需要写出整体坐标X 、Y 和自然坐标ηξ、之间的坐标转换式，即

ξη

αηαξααξηαηαξαα87654321+++=+++=Y X *

四边形四个顶点的坐标值在XOY 坐标系中分别为()()()()44332211,,,,,,,Y X Y X Y X Y X ：在ηξo 坐标系中相应为

()()()()1,1,1,1,1,1,1,1----。将有关数据代入*中的第一式，则有

3214432134321243211,,αααααααααααααααα-+-=+++=--+=+--=X X X X

求解上述方程组得：

,44

321

4432134

321243211X X X X X X X X X X X X X X X X -+-=++--=-++-=

+++=

αααα

坐标变换方程*成为

()()()()[]432111114

X X X X X ξηηξξηηξξηηξξηηξ---+++++--+++--=

同理

()()()()[]432111114

Y Y Y Y Y ξηηξξηηξξηηξξηηξ---+++++--+++--=

当引入函数()ηξ,i N 后，坐标变换方程成为

()()i

i i i

i i Y N Y X N X ηξηξ,,41

1∑∑====

式中()()()ηηξξηξi i i N ++=114

变量ηξ、的正负号由相应节点的坐标值i i ηξ、决定。例如当i=4时，1,14

=-=ηξ，因此，()()()4

11,4

ηξηξ+-=N 。

下面再来研究函数()ηξ,i N 的特性。

对节点()11,1Y X ，相应的自然坐标值为（-1，-1）。从式

()()()ηηξξηξi i i N ++=

114

,中很容易看出，除N 1=1外，N 2=N 3=N 4=0。

对其余各节点也一样。总而言之，对节点i(i=1,2,3,4)，除N i =1外，其余三个N 值均为零。同时，不难看出

()()()()1,,,,4321=+++ηξηξηξηξN N N N ，即四个节点的

N i 函数之

和等于1。

函数()ηξ,i N 具备上章所介绍的形状函数应满足的条件，可作为本单元的形状函数。

采用()ηξ,i N 做形状函数，其位移模式为

()()i i i i i i v N v u N u ηξηξ,,,4

∑∑====

对比

()()i

i i i

i i Y N Y X N X ηξηξ,,4

1∑∑====和()()i i i i i i v N v u N u ηξηξ,,,4

∑∑====可以看出：

在这种实际单元（任意四边形）中，坐标变换式和位移模式不仅采用了相同的形状函数()ηξ,i N ，而且具有相同的数学模型。这种性质的实际单元称为等参数单元。

对用节点位移值u i （或v i 等）求单元内某一点位移量u(或v 等)的插值公式()i i i u N u ηξ,4

1∑==，只要将u(或v 等)换成X(或Y

等)，便成为利用节点值X i (或Y i 等)求相应点坐标X(或Y 等)的插值公式。相反也是这样。 2.几何矩阵[B]

由于几何矩阵[B]通过对位移求偏导数而得出，所以首先必须利用复合函数求导的规则得出下述公式

[]??????????????????=???????????????????????

????????????+

????????+????=??????????????????y u X u J u u Y Y u X X u Y Y u X X

u u u ηξηηξξηξ或写成式中[]????

??????????????=ηη

ξξ

Y X

J ，此式称为雅可比矩阵。

为了将几何矩阵[B]写成变量ηξ、的函数，必须将式

[]??

??????????????=???

???????????????y u X u J u u ηξ改写成 []??????????????????=??????????????????-ηξu u J y u X u 1，同理，[]???

??????????????=???????

???????????-ηξv v J y v X v 1 从表示单元内各点位移与其应变关系的几何方程可知：

{}T

Y v X

v Y

u v u X Y

Y X ?

????

?????????????

????

??=????????????

????????????????

=ξε0110100000010

将式[]??????????????????=??????????????????-ηξu u J y u X u 1和[]????

??????????????=??????????????????-ηξv v J y v X v 1合并，则

[][]T

v v u

J J Y v X

v Y

u X

u ?

???????????????????=???

???????????--ηξ

ηξ

110

对单元（e ），任意一点的位移u ，v 对自然坐标ηξ，的偏导数可利用上式求出，写成矩阵形式为：

[]

{}()e p T

N v v u u ψηξ

??????????? 式中{

}(){}

e v u v u v u v u 11

11=ψ

[][][][][][][][][][]

[][][][][][][]??

???=p

N N N N N N N N N 4321432100000000

对于i=1,2,3,4 []

i i ip N N N ?

???????=ηξ

将

[][]T

v v u

J J Y v X

v Y

u X

u ?

???????????????????=???

???????????--ηξ

ηξ

110

0和[]

{}()e p T

N v v u u ψηξ

???????????代

入

{}T

Y v X

v Y u

u v u X Y Y X ?

????

?????????????

??????=????????????

??????????????

=ξε0110100000010

0，则可得出表示

在整体坐标系中位移和应变关系的几何方程：{}()[](){}()e e e B ψε= 式中的几何矩阵[B]是自然坐标ηξ，的函数：

[][][][]p N J J B ??

??????

?????=--11

00011010000001

也可利用[]

i i ip N N

N ?

???????=ηξ求得的[]ip N 以及

()()i

i i i

i i Y N Y X N X ηξηξ,,4

1∑∑====和

[]????

??????????????=ηη

ξξY X

Y X J 求出[]J ，[][][][][]

{}T

p p p p Y Y Y Y X X X X N N N N J ??

???=43

43214321。

3.单元刚度矩阵[]()e k

设单元板厚为t ，根据虚功方程有：[]()[][][]tdA B D B k A

T e ??=，

此式中几何矩阵[B]和弹性矩阵[D]都已求出。因为几何矩阵[B]中的变量是自然坐标ηξ，，所以也要用自然坐标表示微分面积dA 。

在实际单元中任取一点p ，其整体坐标位X 、Y ，其相应的自然坐标为ηξ，。过p 点做ηξ，的等值线，同时做

ηηξξd d ++，的等值线，围成一小块微分面积dA ，如图4-20

（a ）所示。为便于分析，将四边形pqrs 放大，如图4-20（b ）所示。实际上，ηξd d ，取得很小，因此该四边形可视为

平行四边形。若相邻的两边用向量{}{

}b a ，表示，则两者的乘积等于该平行四边形的面积dA 。

图4-20

b a b a dA ?==θsin

若

b i b b j a i a a y x y x +=+=,，则

()()y

y x y x b a b a j b i b j a i a dA =

++=

为了求出y x y x b b a a ，，,的值，要先写出a 和b 两端节点p 、q 、s 的坐标值。点p ：()()ηξηξ,,,Y Y X X p p == 点q ：()()ηξξηξξ,,,d Y Y d X X q q +=+=

点s ：()()ηηξηηξd Y Y d X X s s

+=+=,,,

利用泰勒技术展开并略去高阶项，可得

()()()()η

ηηξηηξξξ

ηξηξξd X

X d X d X X d X ??+=+??+

=+,,,, 对()()ηηξηξξd Y d Y ++,,,，也可写出相应的展开式。利用式

()()()()η

ηηξηηξξξ

ηξηξξd X

X d X d X X d X ??+=+??+

=+,,,,可得：

ηξξηη

ξξd Y

Y Y b d X X X b d Y Y Y a d X X X a p s y p s x p q y p q x ??=-=??=-=??=-=??=

-=,,，

将此式代入式()()y

x x

y x y x b a b a j b i b j a i a dA =

++=得到：

ηξη

ξd d Y X Y X b a b a dA y

x x

????????==

，简写为ηξd d J dA = 单元刚度矩阵为：

[]

()

[][][][]ηξd d J t B D B k T

e ??--=

111

，这个积分可以采用“数值方法”

，用高斯求积分公式很方便的求出，在此不作介绍。例：求如图所示四边形的雅可比矩阵。

解：求雅可比矩阵可在整体坐标系中进行，也可以在实际单元的局部坐标系中进行。为便于计算，本例在局部坐标系中进行。

对单元（1）：

将四个节点的自然坐标值（-1，-1）、（1，-1）、（1，1）、（-1，1）代入下式：

()()()ηηξξηξi i i N ++=

114

,，再将所得到的()ηξ,i N 值及四个节点实际单元在局部坐标系的坐标值（-3，-2）、（3，-2）、（3，2）、（-3，2）代入下式计算：

()()i

i i i

i i Y N Y X N X ηξηξ,,4

1∑∑====，则ηξ2,3==Y X

雅克比矩阵为[]???????????

??????????????=2003ηηξξ

Y X

Y X J 。对单元（2）：

四个节点在局部坐标系中的坐标值分别为（-1，-3/4）、（1，-3/4）、（1，5/4）、（-1，1/4）。因此，

()()()()()()()()()()()(){}1111111111114

-+-++++-++---=

ηξηξηξηξX ()()()()()()()()?

???????? ??+-+??? ??+++??? ??--++??? ??---=

411145114311431141ηξηξηξηξY 因而雅可比矩阵为

[]()()??

??++=ξη30

41J 也可利用[][][][][]

{}T

p p p p Y Y Y Y X X X X N N N N J ??

????=43

4321

4321求雅可比矩阵，

其结果与上相同，同学们可自行验证。

图

常用单元的刚度矩阵

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 ++= ++= 只要单元尺寸不太大，经过这样处理引起的误差也不大。被积函数又成为常数，可以提出到积分号外面：

从虚功原理推导平面三角单元刚度矩阵

平面问题的三角形单元 ——从能量原理推导刚度矩阵一、虚功原理 1.1虚功如果使力作功的位移不是由于该力本身所引起，即作功的力与相应于力的位移彼此独立，二者无因果关系，这时力所作的功称为虚功。这个位移称为虚位移。 1.2虚位移虚位移指的是弹性体（或结构系）的附加的满足约束条件及连续条件的无限小可能位移。所谓虚位移的"虚"字表示它可以与真实的受力结构的变形而产生的真实位移无关，而可能由于其它原因（如温度变化，或其它外力系，或是其它干扰）造成的满足位移约束、连续条件的几何可能位移。对于虚位移要求是微小位移，即要求在产生虚位移过程中不改变原受力平衡体的力的作用方向与大小，亦即受力平衡体平衡状态不因产生虚位移而改变。 1.3虚功原理处于平衡状态的变形体发生虚位移后，全体外力在对应虚位移所作的外力虚功的等于内力在对应的虚应变上所做的内力虚功。对于一个单元的虚功原理的数学表达式为： {} {}{}{}**T T F d εσΩ?=Ω??? （1-1）

二、平面三角形单元相关矩阵 2.1平面三角形单元得几何和节点描述 3节点三角形单元如图5.1所示。3个节点得编号分别为i、j、m，各自得位置坐标为（）,（），（），各自节点在x方向和y方向的位移为（）,（）,（）。图2-1 2.2三角形单元的位移矩阵对于图5.1所示的平面3节点三角形单元，其位移矩阵为： (2-1)

2.3三角形单元的应变矩阵把位移函数u,v代入几何方程，写成矩阵的形式，则单元上任一点的应变为：（2-2）式（2-16）表示单元节点位移与单元应变的关系。令（2-3）则（2-4）矩阵称为应变矩阵。将其分块可写成：（2-5）式（2-5）表示应变矩阵为常数矩阵，再次证明三节点三角形单元为常应变单元。 2.4 三角形单元的应力矩阵由物理方程：解得：用矩阵表示：（2-6）令：

一般单元在局部坐标系下的单元刚度矩阵

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）

3c 虚功原理推导单元刚度矩阵

§3-3 虚功原理推导梁单元的（单元）刚度矩阵设在力P 的作用下，梁单元i-j 的两端点分别发生了线位移和角位移，用{}e δ来表示梁单元的端点位移（又称结点位移）： { }{}T e i i j j v v δθθ= 使梁单元发生结点位移{}e δ的单元结点力（杆端力）为： { }{}T e i i j j F F M F M = 根据材料力学，如果已知梁的两端点位移，则可求出等截面梁上任意一点的位移（挠度）。即梁上任意一点的位移v(x)可以用{}e δ表示出来，设二者的关系为： {}1234()()()()(){}{} i i T e j j v v x N x N x N x N x N v θδθ?? ???? ==???????? 又设由于某种其他原因，该梁发生了变形，引起梁单元○ e 两端点的位移为（用向量形式表示）： { } * * {}j e i i j v v δθθ= 梁中任意一点的位移为：

{}* ** 1234()()()()(){}{}i i T e j j v v x N x N x N x N x N v θδθ?????? ==???????? 相对于力P 引起的位移v(x),称v*(x)为虚位移计算梁单元○ e 的外力虚功和内力虚功对梁单元来说，两端点的力即是外力，则外力虚功为： **{}{}({}){}e T e e T e ex W F F δδ== 内力虚功 = 虚应变能 2*22* * 222in l l l d v dv dv d v W M d EI d EI dx dx dx dx dx θ??=== ??? ??? ∵ 2 22 2 2 312 422 222{''}{}{}[]{}T T e e e d N d N d N d N d v N B dx dx dx dx dx δδδ?? ===???? 22222** **312422 2 2 2 {''}{}{}[]{}T T e e e d N d N d N d N d v N B dx dx dx dx dx δδδ??===???? ∴ ****[]{}[]{}{}[][]{}{}[][]{}{}[]{} e e in l e T T e l e T T e l e T e e W EI B B dx B D B dx B D B dx k δδδδδδδδ====??? 式中： [][][]e T l k B D B dx =? 虚功原理：系统保持平衡状态的充要条件是外力虚功=内力虚功即： ex in W W = **{}{}{}[]{}e T e e T e e F k δ δδ= 而虚位移为任意、不为零，所以上式等价于：

第三节刚度矩阵(汇编)

第三节刚度矩阵 ——节点载荷与节点位移之间的关系一、单元刚度矩阵 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 ????中的元素也是常

单元刚度矩阵(等参元)MATLAB编程

《有限元法》实验报告专业班级力学（实验）1601 姓名田诗豪学号1603020210 提交日期2019.4.24

实验一（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=0.5*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 %********************************* B31=1/(2*A)*[b(1),0;0,c(1);c(1),b(1)];

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边界上的约束和坐标方向不一致，将给计算带来一定麻烦。此问题如采用下一节的轴对称有限元进行分析，将方便得多。)

文档之家

常用单元的刚度矩阵

常用单元的刚度矩阵

从虚功原理推导平面三角单元刚度矩阵

一般单元在局部坐标系下的单元刚度矩阵

3c 虚功原理推导单元刚度矩阵

最新7.4-单元刚度矩阵组装及整体分析

第三节刚度矩阵(汇编)

单元刚度矩阵(等参元)MATLAB编程

最新7.4-单元刚度矩阵组装及整体分析

ansys单元刚度矩阵的提取

结构刚度矩阵的特点