有限元的matlab编程课题.pptx

有限元编程示例 题目描述: 如下图所示的平面桁架，杆件长度、弹性模量、截面积以及所受节点力P的大小可以自行定义。求节点位移及杆件轴力。 例一：桁架 解题思路： 建立模型 集成总刚 求解位移 求解杆件轴力 输出结果 建立模型： 定义节点坐标 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 模型相关参数输入 H=input(竖杆长度(m):); L=input(水平杆长度(m):); E=input(杆件弹性模量(Gpa):); A=input(杆件截面积(m^2):); a=input(节点力P(kN)：); 节点编号方式 定义单元，即储存单元两端的节点号 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]; 加下划线的为单元编号 集成总刚： 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 ) ; 获取单元两端节点坐标 L = ( (xj-xi)^2 + (yj-yi)^2 )^(1/2) ; 计算杆件长度 形成等效荷载列阵 f=[0;0;0;a;0;0;0;a;0;0;0;a;0;0;0;a;0;0;0;a];%每个节点两个自由度，a为之前输入的节点力 计算从局部坐标到整体坐标的坐标转换矩阵T function T = TransformMatrix( ie )%ie为单元号 c = (xj-xi)/L ; s = (yj-yi)/L ; T=[ c -s 0 0 s c 0 0 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) ;% transpose(T) 为T的转置矩阵2 集成整体刚度矩阵K for ie=1:1:21 %按单元顺序进行循环 k=PlaneTrussElementStiffness(ie); %计算第ie个单元的单刚 m=Element(ie,1); %ie单元的首节点号 n=Element(ie,2); %ie单元的末节点号 K(2*m-1,2*n-1)=k(1,3); K(2*m-1,2*n)=k(1,4); K(2*m,2*n-1)=k(2,3); K(2*m,2*n)=k(2,4); K=zeros(20,20);%用来存储整体刚度矩阵 集成总刚的非对角线元素(这里的元素指2*2的小矩阵） 在下面的集成中，将总刚看成10*10的矩阵，每个元素为2*2的小矩阵 m=Element(ie,2); %ie单元的末节点号 n=Element(ie,1); %ie单元的首节点号 K(2*m-1,2*n-1)=k(3,1); K(2*m-1,2*n)=k(3,2); K(2*m,2*n-1)=k(4,1); K(2*m,2*n)=k(4,2); end 集成总刚的对角线元素（这里的元素指2*2的小矩阵） for i=1:1:10 %按节点的顺序循环