有限元的matlab编程

加下划线的为单元编号
完整ppt
5
形成等效荷载列阵
f=[0;0;0;a;0;0;0;a;0;0;0;a;0;0;0;a;0;0;0;a];%每个节点两个自由度，a 为之前输入的节点力
集成总刚：
获取单元两端节点坐标
xi = Node( Element( ie, 1 ), 1 ) ;%ie为单元号，以下相同 yi = Node( Element( ie, 1 ), 2 ) ; xj = Node( Element( ie, 2 ), 1 ) ; yj = Node( Element( ie, 2 ), 2 ) ;
K(2*i-1,2*i)=K(2*i-1,2*i)+k(3,4);
K(2*i,2*i-1)=K(2*i,2*i-1)+k(4,3);
K(2*i,2*i)=K(2*i,2*i)+k(4,4);
end
end
end
完整ppt
9
求解位移：
u=zeros(20);
根据约束情况修改总刚，采用对角元素置1法
for i=1:1:20
for i=1:1:10 %按节点的顺序循环
for j=1:1:21 %对于每个节点，再按单元的顺序循环
k=PlaneTrussElementStiffness(j);
if Element(j,1)==I %如果i节点为j单元的首节点
K(2*i-1,2*i-1)=K(2*i-1,2*i-1)+k(1,1);
sc 00 0 0 c -s 0 0 s c ];
计算单元刚度矩阵k
k = [ E*A/L 0 -E*A/L
0
0
0
0
0
-E*A/L 0 E*A/L
0
0
0
0
0];
T = TransformMatrix( ie ) ;
k = T*k*transpose(T) ;% transp完o整spept(T) 为T的转置矩阵2
K(2*m-1,2*n-1)=k(3,1);
K(2*m,2*n-1)=k(2,3);
K(2*m-1,2*n)=k(3,2);
K(2*m,2*n)=k(2,4);
K(2*m,2*n-1)=k(4,1);
K(2*m,2*n)=k(4,2);
完整ppt end
8
集成总刚的对角线元素（这里的元素指2*2的小矩阵）
有限元编程示例
完整ppt
1
例一：桁架
题目描述:
如下图所示的平面桁架，杆件长度、弹性模量、截 面积以及所受节点力P的大小可以自行定义。求节 点位移及杆件轴力。
完整ppt
2
解题思路：
• 建立模型 • 集成总刚 • 求解位移 • 求解杆件轴力 • 输出结果
完整ppt
3
建立模型：
模型相关参数输入
H=input('竖杆长度(m):'); L=input('水平杆长度(m):'); E=input('杆件弹性模量(Gpa):'); A=input('杆件截面积(m^2):'); a=input('节点力P(kN)：');
7
集成整体刚度矩阵K
K=zeros(20,20);%用来存储整体刚集成总刚的非对角线元素(这里的元素指2*2的小矩阵）
for ie=1:1:21 %按单元顺序进行循环
k=PlaneTrussElementStiffness(ie); %计算第ie个单元的单刚
计算杆件长度
L = ( (xj-xi)^2 + (yj-yi)^2 )^(1/2) ;
完整ppt
6
计算从局部坐标到整体坐标的坐标转换矩阵T
function T = TransformMatrix( ie )%ie为单元号
c = (xj-xi)/L ; s = (yj-yi)/L ; T=[ c -s 0 0
f=f*1e15;
u=K\f;
完整ppt
10
求解轴力：
获取单元两端的节点号
i = Element( ie, 1 ) ;%ie为单元号 j = Element( ie, 2 ) ;
获取单元两端的节点位移
uElement = zeros( 4, 1 ) ;
uElement( 1:2 ) = u( (i-1)*2+1:(i-1)*2+2 ) ;
m=Element(ie,1); %ie单元的首节点号
n=Element(ie,2); %ie单元的末节点号 m=Element(ie,2); %ie单元的末节点号
K(2*m-1,2*n-1)=k(1,3);
n=Element(ie,1); %ie单元的首节点号
K(2*m-1,2*n)=k(1,4);
Element=zeros(21,2); for i=1:2:7
Element(5/2*i-3/2,:)=[i,i+1]; Element(5/2*i-1/2,:)=[i,i+2]; Element(5/2*i+1/2,:)=[i,i+3]; end for i=2:2:8 Element(5*i/2-1,:)=[i,i+1]; Element(5*i/2,:)=[i,i+2]; end Element(21,:)=[9,10];
uElement( 3:4 ) = u( (j-1)*2+1:(j-1)*2+2 ) ;
K(1,i)=0; K(2,i)=0; K(18,i)=0;
K(i,1)=0;K(i,2)=0; K(i,18)=0;
end %自由度1、2、18被约束了,所在的行和列的其他元素都改为0
K(1,1)=1;%对角线元素置1
K(2,2)=1; K(18,18)=1;
求解
K=K*1e15;%乘以一个大数，减小计算误差
K(2*i-1,2*i)=K(2*i-1,2*i)+k(1,2);
K(2*i,2*i-1)=K(2*i,2*i-1)+k(2,1);
K(2*i,2*i)=K(2*i,2*i)+k(2,2);
end
if Element(j,2)==i %如果i节点为j单元的末节点
K(2*i-1,2*i-1)=K(2*i-1,2*i-1)+k(3,3);
定义节点坐标
Node = zeros(10,2) ;
x=-1*L; %L为横杆长度
for i=1:2:10
x=x+L;
Node(i,:)=[x 0];
end x=-1*L;
节点编号方式
for i=2:2:10
x=x+L;
Node(i,:)=[x H];%H为竖杆长度
end
完整ppt
4
定义单元，即储存单元两端的节点号