《有限元法》课程设计报告.doc

《有限元法》课程设计报告 姓 名： 学 号： 指导教师： 土木工程学院 年 月 日 （一）问题描述 试计算图（a）所示对称拱式三铰桁架的内力、支座反力和各结点的位移。各杆材料相同，E=30GPa，截面积A1=A2=A3=144cm2，A4=A5=A6=A7=180cm2。结构几何尺寸及结点荷载如图（a）所示。 有限元计算模型 ①利用结构的对称性，图（a）所示的结构可简化为图（b）所示的半边结构进行计算； ②确定结点、划分单元，建立整体坐标系与局部坐标系如图（b）所示； ③自由结点位移编码如图（c）所示。 基本理论 ①单元刚度矩阵 ②坐标变换矩阵 计算结果及分析 结点位移（整体坐标） 结点编号 dx dy 1 -8.9912749053E-15 36.245459714 2 -9.5392102578 41.847587165 3 -11.420285459 14.629135802 4 -9.2079563559 34.393883461 5 0 0 桁架变形图 杆端内力（局部坐标） 单元 编号 结点i 结点j N Q N Q 1 2248.4158383 0 -2248.4158383 0 2 2248.4158383 0 -2248.4158383 0 3 1506.5624167 0 -1506.5624167 0 4 3220 0 -3220 0 5 -2553.4904538 0 2553.4904538 0 6 4404.5917455 0 -4404.5917455 0 7 3434.6666667 0 -3434.6666667 0 源程序 %ti4.1 clear; nelem=7;%单元 nnode=5;%节点 E=30000000000;A(1)=144e4;A(2)=144e4;A(3)=144e4;A(5)=180e4;A(6)=180e4;A(7)=180e4;A(4)=180e4;%基本参数 x(1,:)=[0 0];x(2,:)=[150 35];x(3,:)=[300 70];x(4,:)=[150 160];x(5,:)=[300 300];%节点坐标 id(1,:)=[1 2];id(2,:)=[2 3];id(3,:)=[1 4];id(4,:)=[2 4];id(5,:)=[3 4];id(6,:)=[4 5];id(7,:)=[3 5];%单元的节点编码 for i=1:nelem L(i)=sqrt((x(id(i,1),1)-x(id(i,2),1))^2+(x(id(i,1),2)-x(id(i,2),2))^2);%单元长度 end %局部坐标的单元刚度 k=zeros(4,4); for i=1:nelem k(:,:,i)=E*A(i)/L(i)*[1 0 -1 0;0 0 0 0;-1 0 1 0;0 0 0 0]; end %整刚矩阵转变 for i=1:nelem lamd(:,:,i)=[(x(id(i,2),1)-x(id(i,1),1))/L(i),(x(id(i,2),2)-x(id(i,1),2))/L(i); (x(id(i,2),2)-x(id(i,1),2))/L(i),(x(id(i,2),1)-x(id(i,1),1))/L(i)];%方向余弦矩阵 T(:,:,i)=zeros(4,4); T(1:2,1:2,i)=lamd(:,:,i);T(3:4,3:4,i)=lamd(:,:,i);%坐标转换矩阵 kee(:,:,i)=T(:,:,i)*k(:,:,i)*T(:,:,i);%整体坐标的单元刚度 end K=zeros(2*nnode,2*nnode); for i=1:nelem Lo=zeros(4,2*nnode); Lo(1:2,(id(i,1)-1)*2+1:id(i,1)*2)=eye(2,2);%第i单元在n节点的两个分量位置 Lo(3:4,(id(i,2)-1)*2+1:id(i,2)*2)=eye(2,2);%单元定位矩阵 K=K+Lo*kee(:,:,i)*Lo;%组集整刚矩阵 end %荷载计算 R=zeros(2*nnode,1); R(2)=1610;R(4)=3220;R(6)=1610; %位移约束 K(1,1)=K(1,1)*1e15;K(9,9)=K(9,9)*1e15;K(10,10)=K(10,10)*1e15;%乘大数法 det=inv(K)*R; det