平面桁架的有限元法

k =1, 2(单元节点)； kk =1, 2 (节点自由度)；
总节点力号： ss=2(jm[[e, k]]-1)+kk；
Do[Pk[[i]]=…., {i, 4}]; “生成单元节点力” ；
For[k=1, k<=2, k++,
For[kk=1, kk<=2, kk++,
rr =2(k-1)+kk；
平面桁架的有限元法
一、建立有限元模型
平面桁架是现成的有限元模型， 只有节点力
二、单元分析
位移假设，几何方程，物理方程， 单元刚度矩阵
三、总刚度矩阵、总载荷列阵、引入支撑
每个节点建立平衡方程 集成技术（叠加法、对号入座）
四、求解、结果分析，后处理
一、建立有限元模型
平面桁架是现成的有限元模型， 只有节点力
j
k11 k12 k13 k14 ui Rxi
i j
k21 kk3411
k22 k32 k42
k23 k33 k43
k24 k34 k44


vi
u j
v j


Ryi

Rxj
Ryj

i j
5、总刚度矩阵的排列规律
按节点号和节点自由度的顺序（从小到大）排列。
三、总刚度矩阵、总载荷列阵、引入支撑
建立节点号数组（用Mathematica)：
jm=Table[{0, 0}, {i, ne}]; “jm数组名，ne 总单元数;” jm={…，{i节点号， j节点号}，…};
建立单元长度数组
gc=Table[{0}, {i, ne}]; “gc数组名，ne 总单元数;” gc={L1, L2, … Lne};
{ f } [H s ]{a}
{ f得}21 [H s ]24[ Ab ]414{ e}41
a1
u v

1 0
x 0
0 1
0xaa32

a4
{
f
}21
[N
f
]24{
}e 41
节点位移与单元内位移 的关系 { f } [N f ]{ e} 其中 [N (x)] [H s ][ Ab ]1 为形函数。
kk =1, 2, (节 点自由度)
ui
{
e
}


vi
u j

v j
kk=1
k=1
kk=2
k=2
总节点位移向量 {d}的排列规律(从 小到大排序）：
ss=2(jm[[e, k]]-1)+kk; e=单元号
对每个单元 进行节点循 环和自由度 循环，从单 元地址到结 构地址搬家
(EA的两个字母连写作参数名)
{ f } [N f ]{ e}
4、单元刚度矩阵（用虚功原理）
给节点一组虚位移，单元内即产生相应的虚应变，实际节点力 （轴向力、横向力、弯矩）在虚位移上所做的功等于实际应力 （横向剪应力可忽略不计）在虚应变上产生的应变能。
{~ e }
{~} [BL ]{~e}
y y
ui ui cos wi sin
x
wi ui sin wi cos
o
{i} [t]{i}
x
wuii


cos sin
sin cos
wuii

三、总刚度矩阵、总载荷列阵、引入支撑
1、坐标变换
cos sin

ui 1 0 0 0a1

vi
u j
v j


0 1 0
0 b 0
1 0 1
0 0 b
aa32 a4

{ e} [ Ab ]{a}
解线性代数方程组，得 代入 {a} [ Ab ]1{ e}
写成矩阵形式： { }e [ Ab ]{a}
y
x
vj
uj
y
v uj
b
ui i vi
o
x
xi 0, x j b
ui 1 0 0 0a1

vi
u j
v j


0 1 0
0 b 0
1 0源自文库1
0 0 b
aa32 a4
For[k=1, k<=2, k++, For[kk=1, kk<=2, kk++, rr =2(k-1)+kk； ss=2(jm[[e, k]]-1)+kk；deg[[rr]]=d[[ss]]]];
del=Transpose[T[gj[[e]]]].deg; resultant=De.BL.del; se[[e, 1]]=e; se[[e, 2]]=resultant[[j,1]] ];
ss=2(jm[[e, k]]-1)+kk；
Rz[[ss]]=Rz[[ss]+Pk[[rr]]; “叠加法” ；
] ]];
5、总刚度矩阵的排列规律 （1）总刚度矩阵中的行和列均按节点号和节 点自由度的顺序（从小到大）排列；
（2）在总刚度矩阵中，与同一节点相连的单 元的刚度要叠加。
对一个单元而言：
i

([T
e ]{Re})rr

4、总节点载荷向量的排列规律
单元节点力号： rr =2(k-1)+kk；
用Mathematica编程： nj =ne+1总节点数
Pk=Table[0, {i, 4}]; “单元节点力向量” ； Rz=Table[0, {i, 2nj}]; “总节点力向量2nj行” ; For[e=1, e<=ne, e++,
dx
{Re}41 [Ke ]44{δ e}41
用Machematica推导单刚矩阵[Ke]：
用Machematica推导，所得结果以矩阵形式显示
三、总刚度矩阵、总载荷列阵、引入支撑
1、坐标变换
前面的推导是在局部坐标系中进行的，而总刚度
矩阵、总载荷列阵、引入支撑等要在整体坐标中进 行，因此，存在坐标变换问题。

d
ss


([T
e
]{
e
})rr

是单元到 结构的映 射(map)
3、总节点位移向量的排列规律
单元节点位移号：rr =2(k-1)+kk；k =1, 2 (单元节点)； kk =1, 2(节点自由度)
总节点位移号：ss=2(jm[[e, k]]-1)+kk；
例如，已解得结构的总节点位移向量{d}求各单元的应 力。{Se} [De ][BL ]{ e} 用Mathematica：
建立单元方向角数组
gj=Table[{0}, {i, ne}]; “gj数组名，ne 总单元数;” gj={af1, af2, … afne};
3、总节点位移向量的排列规律
jm={{1, 2}，{2, 3}，…，{i节点号， j节点号}，…};
单元节点位 移排列规律：
rr =2(k-1)+kk
k =1, 2 (单元 节点);
y
x P
桁架原型=有限元模型
二、单元分析
平面二力杆单元
y
x
vj
uj
y
j
b
ui i vi
o
x
注意：一个单元的节点数目（ 2）；
一个节点的自由度数目（2）；
有整体坐标系和局部坐标系之分
二、单元分析（先在局部坐标系中进行）
问题：当节点位移一定，相应的节点力为多大？
y
x
vj
uj
y
j
b
ui i vi
o
x
节点力与 节点位移 之间的比 例常数为 刚度系数
x v

y
EA
[De ]


0
{ } [De ]{}
0 EA
{
}


u
x v



u x


[
BL
]{
e
}
y
0
{Se} [De ][ BL ]{ e}
用Machematica推导应力矩阵：
t 0
[T ] 0
t

[t]1 [t]T [T ]1 [T ]T
[T ]{Re} [k ]{ e}
[T ][k]{ e} [k ][T ]{ e}
[k ] [T ][k][T ]T
三、总刚度矩阵、总载荷列阵、引入支撑
2、在总结构和整体坐标中讨论问题，所需要的信息： （1）节点号； （2）单元长度； （2）单元倾角(单元起点到终点的有向线段与整体横轴之夹角）； （4）节点力所在的节点号和作用方向； （5）支撑所在的节点号和作用方向。 3、总节点位移向量的排列规律
deg=Table[0, {i, 4}]; “整体坐标系中的单元节点位移向量” ； del=Table[0, {i, 4}]; “局部坐标系中的单元节点位移向量” ；del={ e} se=Table[{0, 0}, {i, ne}]; “单元内力数组 {e, N} ne行” ;
For[e=1, e<=ne, e++,
[t] sin
cos

{i} [t]{i}
y
x
vj
uj
y
v uj
b
ui i vi
o
x
{R e} [T ]{Re}

i j


t 0
0i
t


j

{Re} [k]{ e}
{R e} [k ]{ e}
{ e} [T ]{ e}
y
x
vj
uj
y
v uj
b
ui i vi
o
x
a1
u 1 v 0
x 0
0 1
0xaa32

a4
u a1 a2 x v a3 a4 x
在节点上满足：
ui a1 vi a3 u j a1 a2b vi a3 a4b
1 x 0 0 [Hs ] 0 0 1 x
y
x
vj
uj
y
v uj
b
ui i vi
o
x
用Machematica推导形函数：
在编程中只关心节点位移 { e} 的系数矩阵
2、应变（从位移求应变） 只有轴向应变
{
}


u v



u
x v
y
{Se} [De ][ BL ]{ e}
{~e}T {Re}
{~}T
l
{Se}dl

{~ e }T
b 0
[
BL
]T
[
De
][
BL
]{δ
e
}dx
Ni
{R
e
}


0 N
j

0
[Ke ]44
b
0[BL
]T42[
De
]22[BL
]24
4、总节点载荷向量的排列规律
按节点号和节点自由度的顺序（从小到大）排列。 设总节点力向量为{Rz}
单元节点位移号： rr =2(k-1)+kk； k =1, 2(单元节点)； kk =1, 2 (节点自由度)；
总节点位移号： ss=2(jm[[e, k]]-1)+kk；

R
zss


Rzss
（2）在总刚度矩阵中，与同一节点相连的单元 的刚度要叠加。用Mathematica编程：
单刚行号：r =2(i-1)+ii；i =1, 2(单元节点)； ii =1, 2 (节点自由度)； 单刚列号：s =2(j-1)+jj；j=1, 2(单元节点)； jj=1, 2 (节点自由度)； 总刚行号：rr=2(jm[[e, i]]-1)+ii； 总刚列号：ss=32(jm[[e, j]]-1)+jj；
N
u
A
Ei
A
u x
{S
e
}


N 0


Ei A

0
u
0 Ei A

x v

y
v 0 y
用内力表示的物理方程
N EA

0



0
N
{Se
}


0

u
0 EA
目标：建立单元刚度矩阵。 方法：结构力学的位移法加虚功原理，
皆用节点位移表示。
1、位移假设： 设各点的位移，
即位移场为
u a1 a2 x v a3 a4 x 其中 a1, a2 ,a4
为待定系数，由节点 坐标和节点位移确定。 将位移写成矩阵形式：
{ f } [H s ]{a}
单刚行号：r =2(i-1)+ii；i =1, 2(单元节点)； ii =1, 2(节点自由度)； 单刚列号：s =2(j-1)+jj；j=1, 2(单元节点)； jj=1, 2(节点自由度)； 总刚行号：rr=2(jm[[e, i]]-1)+ii； 总刚列号：ss=2(jm[[e, j]]-1)+jj； jm=Table[{0, 0}, {i, ne}]; “jm数组名，ne总单元数;” jm={…….，{i节点号，j节点号}，……};



u
x 0


[
BL
]{
e
}
u a1 a2 x v a3 a4 x
3、应力（从应变求应力，本构方程）
{
}


u v


Ei
uv


Ei
u x 0

4、内力（从应力求内力） 二力杆只有轴向力