MATLAB有限元分析与应用

0 200
UU12



F1 F2

200 U3 F3
已知边界条件： U1 0, F2 0, F3 15
100 100 0
100 300 200
0 200
0 U 2



F1 0
f1 = -15.0000 15.0000
u2=U(2:3); f2=SpringElementForces(k2,u2);
f2 = -15.0000 15.0000
12
§3-1 弹簧元
5、实例计算分析应用
k1=SpringElementStiffness(100);
k2=SpringElementStiffness(200);
y = [E*A/L -E*A/L ; -E*A/L E*A/L];
2020/2/29
15
§3-2 线性杆元
2、MATLAB函数编写
2.2 整体刚度矩阵的形成
function y =LinearBarAssemble(K,k,i,j)
%LinearBarAssemble This function assembles the element stiffness

630000 U3 F3
已知边界条件： U1 0,U3 0.002, F2 10
420000 420000 0
420000 1050000 630000
0 0
630000

U2


F110
630000 0.002 F3
%SpringElementForces This function returns the element nodal force
%
vector given the element stiffness matrix k
%
and the element nodal displacement vector u.
%
stiffness matrix for a linear bar with
%
modulus of elasticity E, cross-sectional
%
area A, and length L. The size of the
%
element stiffness matrix is 2 x 2.
K= 100 -100 0 -100 100 0 0 00
K=
K = SpringAssemble(K,k2,2,3)
100 -100 0
2020/2/29
-100 300 -200
9
0 -200 200
§3-1 弹簧元
4、实例计算分析应用
步骤4：引入边界条件
100 100 0
100 300 200

200 U3 15
2020/2/29
10
§3-1 弹簧元
5、实例计算分析应用
步骤5：解方程
U=zeros(2,1); F=[0;15];
K(1,:)=[]; K(:,1)=[];
K = K(2:3,2:3);
300 200
200
200

UU23
3.1 单元刚度矩阵的形成
function y = SpringElementStiffness(k)
%SpringElementStiffness This function returns the element stiffness %matrix for a spring with stiffness k. %The size of the element stiffness matrix is 2 x 2.
19
§3-2 线性杆元 3、实例计算分析应用
步骤2：形成单元刚度矩阵
调用 function y = LinearBarElementStiffness(E,A,L)函数 k1=LinearBarElementStiffness(E,A,L1)
k2=LinearBarElementStiffness(E,A,L2)
第三章 MATLAB有限元分析与应用
§3-1 弹簧元 §3-2 线性杆元 §3-3 二次杆元 §3-4 平面桁架元 §3-5 空间桁架元
§3-6 梁元
2020/2/29
1
§3-1 弹簧元 1、有限元方法的步骤：
离散化域 形成单刚矩阵 集成整体刚度矩阵 引入边界条件 求解方程 后处理
2020/2/29
y = [k -k ; -k k];
2020/2/29
4
§3-1 弹簧元
3、MATLAB函数编写
3.2 整体刚度矩阵的形成
function y = SpringAssemble(K,k,i,j)
%SpringAssemble This function assembles the element stiffness
u2=U(2:3);
2020f/22/=29SpringElementForces(k2,u2)
13
§3-2 线性杆元
1、基本方程
线性杆元也是总体和局部坐标一致的一维有限单元，用线性函数描述
每个线性杆元有两个节点（node）
EA
单刚矩阵为：k


L

EA L

EA L

EA
y = k * u;
2020/2/29
6
§3-1 弹簧元
4、实例计算分析应用
如图所示二弹簧元结构，假定k1=100kN/m，k2=200kN/m，P=15kN。 求：系统的整体刚度矩阵；
节点2、3的位移； 节点1的支反力； 每个弹簧的内力
解：
步骤1：离散化域
2020/2/29
7
§3-1 弹簧元
4、实例计算分析应用
L
总刚矩阵： n n
结构方程： KU F
22
单元节点力： f ku
2020/2/29
14
§3-2 线性杆元
2、MATLAB函数编写
2.1 单元刚度矩阵的形成
function y = LinearBarElementStiffness(E,A,L)
%LinearBarElementStiffness This function returns the element
K = LinearBarAssemble (K,k2,2,3)
2020/2/29
21
§3-2 线性杆元
3、实例计算分析应用
步骤4：引入边界条件
420000 420000 0
420000 1050000 630000
0 630000
UU12



F1 F2
2020/2/29
20
§3-2 线性杆元
3、实例计算分析应用
步骤3：集成整体刚度矩阵
调用 function y = LinearBarAssemble(K,k,i,j)函数
n=3; K = zeros(n,n)
K= 000 000 000
K = LinearBarAssemble (K,k1,1,2)
%
matrix k and the element nodal
%
displacement vector u.
y = k * u;
2020/2/29
17
§3-2 线性杆元
2、MATLAB函数编写
2.4 节点应力计算
function y = LinearBarElementStresses(k, u, A)


0 15
U=K\F U=inv(K)*F
U=Leabharlann Baidu0.1500
2020/2/29
0.2250
11
§2-1 弹簧元
5、实例计算分析应用
步骤6：后处理
U=[0;U]
U= 0
0.1500 0.2250
F=K*U
F= -15.0000 0.0000 15.0000
2020/2/29
u1=U(1:2); f1=SpringElementForces(k1,u1);
%
k is assembled.
K(i,i) = K(i,i) + k(1,1);
K(i,j) = K(i,j) + k(1,2);
K(j,i) = K(j,i) + k(2,1);
K(j,j) = K(j,j) + k(2,2);
y 2020/2/29 = K;
16
§3-2 线性杆元
2、MATLAB函数编写
n=3;
K=zeros(n,n);
K=SpringAssemble(K,k1,1,2);
K=SpringAssemble(K,k2,2,3);
U=zeros(2,1);
F=[0;15];
K = K(2:3,2:3);
KK=K;
U=K\F
U=[0;U];
F=K*U;
u1=U(1:2);
f1=SpringElementForces(k1,u1)
步骤2：形成单元刚度矩阵
调用 function y = SpringElementStiffness(k)函数
k1=SpringElementStiffness(100);
k1 =
100 -100 -100 100
k2=SpringElementStiffness(200)；
k2 =
200 -200
%
matrix k of the linear bar with nodes i and j
%
into the global stiffness matrix K.
%
This function returns the global stiffness
%
matrix K after the element stiffness matrix
2020/2/29-200 200
8
§3-1 弹簧元
4、实例计算分析应用
步骤3：集成整体刚度矩阵
调用 function y = SpringAssemble(K,k,i,j)函数
n=3; K = zeros(n,n);
K= 000 000 000
K = SpringAssemble(K,k1,1,2)
%LinearBarElementStresses This function returns the element nodal
%
stress vector given the element stiffness
%
matrix k, the element nodal displacement
%
vector u, and the cross-sectional area A.
y = k * u/A;
2020/2/29
18
§3-2 线性杆元
3、实例计算分析应用
如图所示二线性杆元结构，假定E=210MPa，A=0.003m^2，P=10kN， 节点3的右位移为0.002m。
求：系统的整体刚度矩阵； 节点2的位移； 节点1、3的支反力； 每个杆件的应力
解：
步骤1：离散化域
2020/2/29
2
§3-1 弹簧元
2、基本方程
弹簧元是总体和局部坐标一致的一维有限单元 每个弹簧元有两个节点（node）
单刚矩阵为：
k

k k
k
k

总刚矩阵： n n
结构方程： KU F
22
单元节点力： f ku
2020/2/29
3
§3-1 弹簧元
3、MATLAB函数编写
%
matrix k of the spring with nodes i and j into the
%
global stiffness matrix K.
%
This function returns the global stiffness matrix K
%
after the element stiffness matrix k is assembled.
2020/2/29
22
§3-2 线性杆元 3、实例计算分析应用
步骤5：解方程
U=zeros(1,1); U3=0.002 F=[-10]; K = K(2,2)
2.3 节点载荷计算
function y = LinearBarElementForces(k,u)
%LinearBarElementForces This function returns the element nodal
%
force vector given the element stiffness
K(i,i) = K(i,i) + k(1,1);
K(i,j) = K(i,j) + k(1,2);
K(j,i) = K(j,i) + k(2,1);
K(j,j) = K(j,j) + k(2,2);
y = K;
2020/2/29
5
§3-1 弹簧元
3、MATLAB函数编写
3.3 节点载荷计算
function y = SpringElementForces(k,u)