考虑剪切变形的空间梁单元有限元分析

《有限元方法》
2023年课程大作业
*名：***
学号：XS********指导老师：**教授
*********大学
二〇二三年五月
目录
1作业内容与要求 (2)
2考虑剪切变形空间梁的势能泛函和梁单元位移插值函数，并通过泛函变分原理推导得到梁单元的刚度方程 (2)
3、给出matlab编写的有限元程序源代码以及对应的说明，包括 (5)
1) 有限元程序整体架构，计算的流程图； (5)
2) 离散化与编号； (8)
3) 数据的准备； (8)
4) 单元分析形成单元刚度矩阵； (11)
5) 组装整体刚度矩阵； (13)
6) 节点载荷的计算； (14)
7) 边界条件的处理及刚度方程求解； (17)
8) 应力和应变的计算；.......................................................... 错误!未定义书签。

9) 和无开口梁情况分析比较 (18)
10) 结论。

(20)
1作业内容与要求
如图1所示，某悬臂梁，有三个矩形开口，梁端承受一集中载荷F ，大小为60000N 。

试采用考虑剪切变形的2结点空间梁单元，用matlab 编写有限元程序，计算梁端A 点的线位移，并和无开口梁的结果进行比较。

悬臂梁材料的杨氏模量
，泊松比。

图1 矩形开口悬臂梁结构末端受集中力F 作用示意图
要求：
一、纸质报告A4单面打印，随笔试答卷一起上交； 二、报告格式要求同研究生学位论文；
2考虑剪切变形空间梁的势能泛函和梁单元位移插值函数，并通过泛函变分原理推导得到梁单元的刚度方程
首先把空间梁单元分解成一个杆单元、一个轴单元，和xoy ，xoz 平面内的两个平面梁，由此可以吧得到空间梁单元的能量泛函为：
2220002220
0000111
Π(,,,)d d d 222
111
d ()d ()d 222 d d l l l x y xoz xoz z xoy xoz l l l x xoy xoy l l z z j iz i z j
i
m
j GA u v w EI x x EI x k d GA du x EA x JG x k dx dx
q w x P w M qv x P θκγκθγθ=⎰+⎰+⎰+
⎰+⎰+⎰-⎰-∑+∑-⎰-∑00 d -d ym m n n l l x z l x x k x k x
k
n
l l
v M f u x F u m x M θθθ+∑⎰-∑-⎰-∑ (1)
首先给出如下的形函数（插值函数）
()()231233
22332
3
3
456
1322321
N N l N N l N N ξξξξξξξξξξξ=-+=-+=-=-+=-= (2)
其中
,01i
x x l
ξξ-=
≤≤ (3)
对于考虑剪切变形的梁单元将挠度分解为
b S w w w =+
(4)
则节点位移为
112212b e b b s e
s s w q w w q w θθ⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦⎡⎤=⎢⎥
⎣⎦
(5)
其中
1212
d d ,d d b b w w x x θθ⎛⎫⎛⎫== ⎪ ⎪⎝⎭⎝⎭
(6)
考虑将s w 采用2结点的Lagrange 插值，即线性插值； b
w 采用与不考虑剪切变形梁单元的w 相同的Hermite 插值，则有
11213242
5162
b b b
s
s
s w N w N N w N w N w N w θθ=+++=+ (7)
利用平衡方程将s
w 凝聚掉，可以得到最终平面梁单元的能量泛函为1Π()2
T
e q q K q Fq =- 其中
()()1122
23
21220[,,,]
12
61266(4)6(2)12
6126(1)6(2)6(4)d //e e T T T
j j b l k l l l b l l b l EI K l l b l l b l l b w l P N ql N P dN d M q l
w ξξξθθξ-⎡⎤⎢⎥+--⎢⎥=⎢⎥
---+⎢⎥
--+⎣⎦
=⎰+∑-∑⋅= (8)
对于杆单元
5162u N u N u =+
(9)
杆单元的能量泛函为1Π()2
T
e q q K q Fq =- 其中Ke 杆的刚度矩阵
()0
2111111[,]
d e T T
x j j x EA Ke l P N ql N P q u u ξξ-⎛⎫=
⎪-⎝⎭
=⎰+∑= (10)
同理可以得到轴的
5162x N N θθθ=+
(11)
能量泛函为1Π()2
T
e q q K q Fq =- 其中Ke 杆的刚度矩阵
()0
2111111[,]
d e T T
x j j x GJ Ke l
P N ql N P q θξθξ-⎛⎫=
⎪-⎝⎭
=⎰+∑= (12)
最后将这些进行组合可以得到总的能量泛函为 1
Π()2
T e q q K q Fq =-
(13)
其中
333
333
3
3
22333
3
23
126126(1)(1)(1)(1)126126(1)(1)(1)(1)6(4)6(2)(1)(1)(1)(1)6(4)(1)(1)z z z z y y y y y
y y y z z EA
EA l l
EI lEI EI lEI b l b l b l b l EI lEI EI lEI b l b l b l
b l GJ GJ
l l
lEI b l EI lEI b l EI b l b l b l b l lEI b l EI b l b --++++---
++++-+--
+++++++233
33
3
333
33
3
22333
(2)6(1)(1)126126(1)(1)(1)(1)126126(1)(1)(1)(1)((2)(4)6(1)(1)(1)6(y z z z z z y y y y y y z z b l EI lEI l b l b l EA EA l l
EI lEI EI lEI b l b l b l b l EI lEI EI lEI b l b l b l b l GJ
GJ l l
b l EI b l EI lEI b l b l b l lEI --++---++++-++++---+-
+++3
11111122222222
3
33
3((2)66(4)1 [,,,,,,,,,)(1)(1)(1,,)(1)z z z z x y z x y z b l E u I lEI lEI b l EI b l b l b l b l b w l q v u v w θθθθθθ⎡⎤⎢⎥⎢
⎥
⎢
⎥⎢⎥
⎢
⎥⎢
⎥⎢⎥⎢
⎥⎢
⎥⎢⎥⎢
⎥⎢⎥⎢⎥⎢⎥
⎢⎥
⎢⎥
⎢
⎥
⎢⎥⎢⎥⎢⎥
⎢⎥⎢⎥⎢⎥
⎢⎥⎢⎥⎢
⎥⎢⎥⎢
⎥⎢
⎥⎢⎥⎢⎥⎢
⎥⎢⎥-+-+++++⎥⎢⎥⎢⎣
⎦
=]
(14)
对能量泛函
1
Π()
2
T
e
q q K q Fq
=-取变分为零即
Π()0
e
q K q F
=-=(15)
即空间梁单元的刚度方程
3、给出matlab编写的有限元程序源代码以及对应的说明，包括
1) 有限元程序整体架构，计算的流程图；
程序的流程图如下图所示
global E mu G; %将泊松比和杨氏模量声明为全局变量
E = 2.0e11; %给杨氏模量赋予指定值
mu = 0.25; %给泊松比赋予指定值
G = E/(1+mu)/2; %有泊松比杨氏模量计算剪切模量
L=[15,30,15,30,15,30,15]*10^(-2);
A=[30*10,15*10,30*10,15*10,30*10,15*10,30*10]*10^(-4);
Iy=[10*(30^3),10*(30^3)-10*(15^3),10*(30^3),10*(30^3)…
-10*(15^3),10*(30^3),10*(30^3)-10*(15^3),10*(30^3)]/12*10^(-8);
2) 离散化与编号；
在本次作业中根据梁的截面变化将梁分成七个单元从左到右依次编号为①②③④…⑦,为了保证最后组集得到最终刚度矩阵的稀疏性，减小存储矩阵的半带宽，节点编号也是采用从左到右依次1,2,3,4，…8的编号方法，具体离散化和单元以及节点编号的示意图如下图所示：
图一、悬臂梁单元的离散以及单元节点编号示意图
这一部分程序体现在main—function主函数中。

3) 数据的准备；
数据准备部分包括梁单元的划分、单元编号、节点编号、单元长度、单元横截面积、单元的转动惯量、惯性矩Iy，Iz，每个单元的界面剪切修正系数k等，梁单元截面比较简单，所以单元编号、单元长度、转动惯量、惯性矩要么可以直接获取，要么可以通过简单的公式计算的得到，这些部分的数据准备直接体现在主函数中，但界面截切修正系数没有直接现成的公式，相对比较复杂，所以对单元的剪切截面修正系数进行推导计算，具体过程如下所示：
求界面的剪切修正系数
2
2
2A A S k dA I
b =⎰
(16)
2
2
225
322222144
12h b h S b y y y A bh I bh bh ⎡⎤⎛⎫⎛⎫⎛⎫=-+=-⎢⎥
⎪⎪ ⎪⎝⎭⎝⎭⎝⎭⎢⎥⎣⎦
==⎛⎫ ⎪⎝⎭
(17)
所以
2
222
2
2
5
2
2
2
422
2
452
h 52352h 52
144144144 416236 16235111 3628385166 5
h b h b h h b h k y dxdy bh
b b h h y y dy
bh h h y y h --
-
⎧⎫⎡⎤⎪⎪=-⎨⎬⎢⎥⎣⎦⎪⎪⎩⎭
⎡⎤=-+⎢⎥⎣⎦
⎡⎛⎫=+-+⎢ ⎪
⋅⎢⎝⎭⎣⎡⎤=-+
⎢⎥***⎣⎦=
⎰
⎰⎰ (18)
同理，对于中间去掉一块的矩形
其界面剪贴修正系数为
2
2223222
322222144
1()(2)h b h S b y y y A bh I h w h bw bh b hw w w ⎡⎤⎛⎫⎛
⎫⎛⎫=-+=-⎢⎥
⎪⎪ ⎪⎝⎭⎝-⎭⎝⎭⎢⎥⎣⎦
-==
⎛⎫ ++⎪
⎭
-⎝ (19)
w 222222222222
222222222w 222222
2
2
w
22144144
144144
1 ) 2*()((()44)b
h b h b b h b h y dxdy b k b h y dxd b b h w h hw w b h w h hw y b h b w y ----⎛⎫⎧⎫⎡⎤⎪⎪ ⎪-⎨⎬⎢⎥ ⎪⎣⎦⎪⎪⎩⎭ ⎪= ⎪ ⎪⎧⎫⎡⎤⎪⎪ ⎪+-⎨⎬⎢⎥ ⎪⎣⎦⎪⎪ ⎪⎩⎭
⎝⎭
⎧⎡⎤=-⎨⎦-++-++⎢⎥⎣⎰⎰⎰⎰⎰根据积分区域的对称性可以简化为
()
()2
2
422
24w 2
h 4222222352
w 2334552
2
222
144 2*416272
*() 3223572
() =
32638(53)2(())()()b b h
h b h w dxdy b h h y y dy h w h y y h h h h h w h w h w h w h hw w
w w h w h hw w -
⎫⎪⎪⎬⎪⎪⎩⎭
⎡⎤=-+⎢⎥⎣⎦
⎡⎤⎛⎫-=
+-+⎢
⎥ ⎪⋅⎢⎥⎝⎭⎣⎦
⎡⎤---⎢⎥+⎥-++-+
⨯⨯⨯⎢-++⎣++⎦
⎰⎰(20) 在此问题中，带入h=30，w=15可得 1.305102041k =
4) 单元分析形成单元刚度矩阵；
根据前面推导得到的空间梁的刚度方程形成了具体的刚度矩阵，这一部分工作主要是通过子程序Beam3D2Node_Stifness(A,J,I,b,L)完成的。

5) 组装整体刚度矩阵；
依照节点编号将梁单元的刚度矩阵组装为整体的刚度矩阵主要是通过子程序Beam3D2Node_Assembly_Stifness(K_temp,k,i,j)完成的。

6) 节点载荷的计算；
根据之前推导刚度方程的过程，可以用Beam3D2Node_Elem_Equivload(L,pz,py,px,Mx)子程序计算得到各个单元的等效载荷，并通过函数
Beam3D2Node_Assembly_Equi(F_equi,f_equi,i,j)根据各个单元的节点编号将各个单元的等
通过对整个梁进行受力分析，其中在处理外载荷F的时候利用了力的平移等效原则，将F移动到梁的正中间，此时会有一个附加的力矩M=F*r=60000*5/100N.m=3000N.m,以得到梁所受到的外力大小，梁的受力分析如下图所示：
根据受力分析情况，形成空间梁所受外载荷的载荷向量[Fx,Fy,Fz,Mx,My,Mz,…,0, 0,-F ，M ，0，0]，将外部载荷向量和整体等效载荷向量相加即可得到最终的载荷向量，这个部分的工作主要体现在主程序中。

7) 边界条件的处理及刚度方程求解；
为了不改变刚度矩阵的规模，悬臂梁的左端是固定端约束，采用置一法的方法处理边界条件，具体操作原理为
11
121111,212222221
2
0000100
n i n i i i i
n n nn n n ni i K K K u F K u K K K u F K u u u K K K u F K u -⎡⎤⎧⎫⎧⎫
⎢⎥⎪⎪⎪⎪-⎢⎥⎪⎪⎪⎪⎢⎥⎪⎪⎪⎪⎪⎪⎪⎪
=⎢⎥⎨⎬⎨⎬⎢⎥⎪⎪⎪⎪
⎢⎥⎪⎪⎪⎪⎢⎥⎪⎪⎪⎪-⎢⎥⎪⎪⎪⎪⎣⎦⎩⎭⎩⎭
在该问题中，固定端的广义位移都为零，只需要将刚度矩阵的相应行和列置为1和置零操作，，将相应的广义位移直接置零，但在具体操作的过程中发现无法处理未知力，求解困难，最终还是采用了划行划列的思想，把部分矩阵直接拿出来然后求解，具体的过程和该做要工作体现在主函数中，对处理过后得到梁的总的刚度方程，采用gauss 消元的方法，即在matlab 中利用集成后的节点载荷F ‘\’总的刚度矩阵K ，对各个节点的广义位移进行求解具体结果体现在主程序中。

8) 和无开口梁情况分析比较
为了直观体现有开口梁和无开口梁在外载荷作用下两者广义位移的变化，将两者的结果通过如下的对比图展示(由于没有给定材料密度，所以转动惯量无法计算，在此问题中，假设密度为1，绕y方向转动的角度没有实际意义，只用图来表示两者的变化趋势)
9) 结论。

u v wθθθ,通过利用有限元计算，得到端点处的广义位移如下图所示从左到右分别是,,,,,
x y z 其中长度单位是m，角度单位是弧度
从挠度角度来看，在外载荷和约束的相互作用下，无开口的梁在z方向的挠度为0.0596938cm，有开口梁挠度为0.0649792cm，开口梁在一定程度上破坏了梁的抗弯、抗扭性能，但从该数据来看，破坏的程度不是很大，但开口可以大大减轻梁的自身质量，所以在现代的工程中有也能看到具有开口的梁部件，特变是在大型工程中。

第20 页。