kl800.com省心范文网

基本平面刚架MATLAB程序2013


% 平面刚架 MATLAB 程序 % 2003.9.16 2007.2.28 2008.4.1 2009.10 2011.10 2013.9 %************************************************* % 变量说明 % NPOIN NELEM NVFIX NFPOIN NFPRES % 总结点数,单元数, 约束个数, 受力结点数, 非结点力数 % COORD LNODS YOUNG % 结构节点坐标数组, 单元定义数组, 弹性模量 % FPOIN FPRES FORCE FIXED % 结点力数组,非结点力数组,总体荷载向量, 约束信息数组 % HK DISP % 总体刚度矩阵,结点位移向量 %************************************************** format short e %设定输出类型 clear %清除内存变量 FP1=fopen('6-6.txt','rt') %打开初始数据文件 %读入控制数据 NELEM=fscanf(FP1,'%d',1); %单元数 NPOIN=fscanf(FP1,'%d',1); %结点数 NVFIX=fscanf(FP1,'%d',1); %约束数 NFPOIN=fscanf(FP1,'%d',1); %作用荷载的结点个数 NFPRES=fscanf(FP1,'%d',1); %非结点荷载数 YOUNG=fscanf(FP1,'%f',1); %弹性模量 % 读取结构信息 LNODS=fscanf(FP1,'%f',[4,NELEM])' % 单元定义: 左、右结点号,面积,惯性矩(共计 NELEM 组) COORD=fscanf(FP1,'%f',[2,NPOIN])' % 坐标: x,y 坐标(共计 NPOIN 组) FPOIN=fscanf(FP1,'%f',[4,NFPOIN])' % 节点力(共计 NFPOIN 组) :受力结点号、X 方向力(向右正) , % Y 方向力(向上正) ,M 力偶(逆时针正) FPRES=fscanf(FP1,'%f',[4,NFPRES])' % 均布力(共计 % NFPRES 组) :单元号、荷载类型、荷载大小、距离左端长度 FIXED=fscanf(FP1,'%d',NVFIX)' % 约束信息:约束对应的位移编码(共计 NVFIX 组) %--------------------------------------------------------HK=zeros(3*NPOIN,3*NPOIN); % 张成总刚矩阵并清零 FORCE=zeros(3*NPOIN,1); % 张成总荷载向量并清零 %形成总刚 for i=1:NELEM % 对单元个数循环 % 生成局部单刚(局部坐标) 右手坐标系 EK=ele_EK(i,LNODS,COORD,YOUNG); T=zbzh(i,LNODS,COORD); % 坐标转换矩阵

EKT=T’*EK*T; % 生成整体单刚(整体坐标系) % 组成总刚 按 3*3 子块加入总刚中(共计 4 块) for j=1:2 %对行进行循环---按结点号循环 N1=LNODS(i,j)*3; % j 结点第 3 个位移的整体编码 for k=1:2 %对列进行循环---按结点号循环 N2=LNODS(i,k)*3; % k 结点第 3 个位移的整体编码 HK((N1-2):N1,(N2-2):N2)=HK((N1-2):N1,(N2-2):N2)... +EKT(j*3-2:j*3,k*3-2:k*3); End % 单刚 3 x 3 子块叠加到总刚中 end end % 由结点力与非结点力生成总荷载向量列阵 for i=1:NFPOIN % 对结点荷载个数进行循环 N1=FPOIN(i,1); % 作用荷载的结点号 N1=N1*3-3; % 该结点号对应第一个位移编码 - 1 for j=1:3 FORCE(N1+j)=FORCE(N1+j)+FPOIN(i,j+1);%取结点荷载 end end % 计算由非结点荷载引起的等效结点荷载 for i=1:NFPRES % 对非结点荷载个数进行循环 F0=ele_FPRES(i,FPRES,LNODS,COORD); %计算单元固端力 % 对单元局部杆端力要进行坐标转换 ele=FPRES(i,1); % 取荷载所在的单元号 T=zbzh(ele,LNODS,COORD); % 坐标转换矩阵 F0=T’*F0; NL=LNODS(ele,1); NR=LNODS(ele,2); %单元的左右结点号 % 将单元固端力变成等效结点荷载(注意固端力与等效结点荷载符号相反) FORCE((3*NL-2):3*NL)=FORCE((3*NL-2):3*NL)-F0(1:3); FORCE((3*NR-2):3*NR)=FORCE((3*NR-2):3*NR)-F0(4:6); end % 总刚、总荷载进行边界条件处理----划 0 置 1 for j=1:NVFIX % 对约束个数进行循环 N1=FIXED(j); HK(N1,1:3*NPOIN)=0; HK(N1,N1)=1; % HK(1:3*NPOIN,N1)=0; % 将零位移约束对应的行、列变成零,主元变成 1 FORCE(N1)=0; end %--------------------------------------------------------DISP=HK\FORCE % 方程求解,HK 先求逆再与力向量左乘 %--------------------------------------------------------% 求结构各个单元内力 EDISP=zeros(6,1); % 单元位移列向量清零

for i=1:NELEM % 对单元个数进行循环 for j=1:2 %对杆端循环 % i 单元左右端结点号*3 = 该结点的最后一个位移编码 N1=LNODS(i,j)*3; % 取一端的单元位移列向量 EDISP(3*j-2:3*j)=DISP(N1-2:N1); end % 生成局部单刚(局部坐标) 右手坐标系 EK=ele_EK(i,LNODS,COORD,YOUNG); T=zbzh(i,LNODS,COORD); % 坐标转换矩阵 FE=EK*T*EDISP; %计算局部坐标杆端力(由结点位移产生) for j=1:NFPRES if FPRES(j,1) == i %成立时,当前单元上有非结点荷载 F0=ele_FPRES(j,FPRES,LNODS,COORD);%单元固端力 FE=FE+F0; % 考虑由非结点荷载引起的杆端力 end end FE % 打印杆端力 end %------------------------------------------------------ele_FPRES.m %计算单元固端力函数( 正方向:X 向右 Y 向上 M 逆时针) % 入口参数:荷载序号,荷载信息,单元信息,结点坐标 % 出口参数:单元固端力——左右两端的轴力、剪力、弯矩 function F0=ele_FPRES(iFPRES,FPRES,LNODS,COORD) ele=FPRES(iFPRES,1); %取荷载所在的单元号 G=FPRES(iFPRES,3); %单元荷载大小 C=FPRES(iFPRES,4); %单元荷载与左端距离 NL=LNODS(ele,1); NR=LNODS(ele,2); %单元的左右结点号 dx=COORD(NR,1)-COORD(NL,1); % x 坐标差 dy=COORD(NR,2)-COORD(NL,2); % y 坐标差 L=sqrt(dx^2+dy^2); %单元长度 % 计算公式中一些常出现的项 D=L-C; C1=C/L; C2=C1*C1; C3=C1*C2; B1=D/L; B2=B1/L; F0=[0;0;0;0;0;0]; %单元固端力清零 switch FPRES(iFPRES,2) case 1 %均布荷载 F0(2)=-G*C*(2-2*C2+C3)/2.0; F0(3)=-G*C*C*(6-8*C1+3*C2)/12.0; F0(5)=-G*C-F0(2); F0(6)=G*C*C*C1*(4-3*C1)/12.0; case 2 %横向集中力 F0(2)=-G*B1*B2*(L+2*C);

F0(3)=-G*C*B1*B1; F0(5)=-G*C2*(L+2*D)/L; F0(6)=G*D*C2; case 3 %纵向集中力 F0(1)=-G*B1; F0(4)=-G*C1; end return ele_EK.m % 计算单元刚度矩阵函数 EK % 入口参数:单元号、单元信息数组、结点坐标、弹性模量 % 出口参数:局部单元刚度矩阵 EK function EK=ele_EK(i,LNODS,COORD,E) NL=LNODS(i,1); NR=LNODS(i,2); %左右结点号 dx=COORD(NR,1)-COORD(NL,1); % x 坐标差 dy=COORD(NR,2)-COORD(NL,2); % y 坐标差 L=sqrt(dx^2+dy^2); %单元长度 A=LNODS(i,3); I=LNODS(i,4); %面积;惯性矩 % 生成单刚(局部坐标) 右手坐标系 EK =[E*A/L 0 0 -E*A/L 0 0;... 0 12*E*I/L^3 6*E*I/L^2 0 -12*E*I/L^3 6*E*I/L^2;... 0 6*E*I/L^2 4*E*I/L 0 -6*E*I/L^2 2*E*I/L;... -E*A/L 0 0 E*A/L 0 0;... 0 -12*E*I/L^3 -6*E*I/L^2 0 12*E*I/L^3 -6*E*I/L^2;... 0 6*E*I/L^2 2*E*I/L 0 -6*E*I/L^2 4*E*I/L]; return %-------------------------------------------------------zbzh.m % 形成第 i 单元的坐标转换矩阵函数 T(6,6) % 入口参数:单元号,单元信息,结点坐标 % 出口参数:坐标转换矩阵(整体向局部投影) function T=zbzh(i,LNODS,COORD) NL=LNODS(i,1); %左结点号 NR=LNODS(i,2); %右结点号 dx=COORD(NR,1)-COORD(NL,1); % x 坐标差 dy=COORD(NR,2)-COORD(NL,2); % y 坐标差 L=sqrt(dx^2+dy^2); % 单元长度 c=dx/L; % cos a (与 x 轴夹角余弦) s=dy/L; % sin a T=[ c s 0 0 0 0;... -s c 0 0 0 0;... 0 0 1 0 0 0;... 0 0 0 c s 0;... 0 0 0 -s c 0;... 0 0 0 0 0 1]; return


赞助商链接

2014基本平面刚架MATLAB程序2014

2014基本平面刚架MATLAB程序2014_计算机软件及应用_IT/计算机_专业资料。% 平面刚架 MATLAB 程序 % 2003.9.16 2007.2.28 2008.4.1 2009.10 2011.10 2013....

基本平面刚架MATLAB程序2013

基本平面刚架MATLAB程序2013 - % 平面刚架 MATLAB 程序 % 2003.9.16 2007.2.28 2008.4.1 2009.10 2011.10 2013.9 %***...

平面刚架MATLAB程序

平面刚架MATLAB程序_计算机软件及应用_IT/计算机_专业资料 暂无评价|0人阅读|0次下载|举报文档 平面刚架MATLAB程序_计算机软件及应用_IT/计算机_专业资料。...

2016基本平面刚架MATLAB程序

2016基本平面刚架MATLAB程序_理学_高等教育_教育专区。专业课课件 % 平面刚架 MATLAB 程序 % 2003.9.16 2007.2.28 2008.4.1 2009.10 2011.10 2013.9 ...

2016基本平面刚架各种荷载MATLAB程序

2016基本平面刚架各种荷载MATLAB程序 - % 平面刚架 MATLAB 程序 % 2003.9.16 2007.2.28 2008.4.1 2009.10 2011.10 2013.9 2...

基于MATLAB的平面刚架静力分析

工程计算实践作业 基于 MATLAB平面刚架静力分析为了进一步理解有限元方法计算的过程, 本文根据矩阵位移法的基本原理应 用 MATLAB 编制计算程序对以平面刚架结构...

基本平面刚架MATLAB程序2011

基本平面刚架MATLAB程序2011 隐藏>> % 平面刚架 MATLAB 程序 % 2003.9.16 ...现在就登录>> 你可能喜欢©2013 Baidu 使用百度前必读 | 文库协议...

有限元钢架结构分析~手算+matlab+ansys模拟

(2)编制通用平面钢架分析有限元 Matlab 程序,并计算盖提,与手工结果进行比较; (3)利用 Ansys 求解,表格列出 B 点的位移,A、C 处支反力,绘制弯矩图、剪力图...

有限元平面钢架

(2)编制通用平面钢架分析有限元 Matlab 程序,并计算...4 第二章原理说明一、平面刚架的解题思路(步骤) ...有限元和计算机发展共同构成了现代计算力学的基础。...

有限元刚架结构matlab程序

有限元刚架结构matlab程序_机械/仪表_工程科技_专业资料。clc clear format compact format shortG jd=input('请输入节点数:'); dy=input('请输入单元数:');...