matlab.rar_MATLAB 桁架_桁架_桁架 matlab_桁架 matlab_桁架‘

%平面桁架静力计算程序 format short %数据显示设置 %从文件中读取数据 fid = fopen('data.txt','rt'); if fid == -1 disp('数据文件打开错误！'); return end %结构信息 [N,count] = fscanf(fid,'%d',1);%结点总数 [NM,count] = fscanf(fid,'%d',1);%单元总数 [NR,count] = fscanf(fid,'%d',1);%受约束的结点总数 DP = 10000; %一个大数 M = 2 * N; %结构自由度总数 %单元信息 [JD,count] = fscanf(fid,'%d',[2 NM]) %单元结点编号 [EA,count] = fscanf(fid,'%f',[2 NM]) %单元材料常数，第一行为弹性模量，第二行为橫截面面积 [XY,count] = fscanf(fid,'%f',[2 N]) %结点整体坐标 [LR,count] = fscanf(fid,'%d',[2 NR]) %约束结点的整体编号及与该结点相连接的单元总数 NLM = sum(LR(2,:)); %读入LM数据的总数 [LM,count] = fscanf(fid,'%d',NLM) %与NR个约束结点相连接的单元号 [DD,count] = fscanf(fid,'%f',[2 NR]) %约束结点的x和y方向已知位移，无约束则为较大的数(>10000) [P,count] = fscanf(fid,'%f',M) %存储沿整体坐标x和y方向的所有结点荷载值 fclose(fid); %统一单位 EA(1,:)=EA(1,:)*1e3; EA(2,:)=EA(2,:)*1e-4; %待求信息 U = zeros(1,M); %结点位移 FN = zeros(1,NM); %单元内力 RE = zeros(2,NR); %支座反力，第一行为x方向，第二行为y方向 EL = zeros(1,NM); %单元线刚度EA/L RF = zeros(2,NM); %单元方向正余弦值 %计算并存储单元参数 for IE = 1:NM %对单元IE NI = JD(1,IE); %单元的I端结点编号 NJ = JD(2,IE); %单元的J端结点编号 XI = XY(1,NI); %NI结点的x坐标值 YI = XY(2,NI); %NI结点的y坐标值 XJ = XY(1,NJ); %NJ结点的x坐标值 YJ = XY(2,NJ); %NJ结点的y坐标值 AL = sqrt((XJ-XI)^2+(YJ-YI)^2); %单元的长度 EL(IE) = EA(1,IE) * EA(2,IE) / AL; %单元的线刚度 RF(1,IE) = (YJ - YI) / AL; %与x轴夹角正弦值 RF(2,IE) = (XJ - XI) / AL; %与x轴夹角余弦值 end %整体刚度矩阵的集成 AK = zeros(M); %结构初始整体刚度矩阵 for IE = 1:NM %对单元IE SS = RF(1,IE) ^ 2; CC = RF(2,IE) ^ 2; CS = RF(1,IE) * RF(2,IE); EAL = EL(IE); NI = JD(1,IE); NJ = JD(2,IE); II = 2*NI; II1 = II-1; JJ = 2*NJ; JJ1 = JJ-1; AK(II1,II1) = AK(II1,II1)+EAL*CC; AK(II1,II) = AK(II1,II)+EAL*CS; AK(II,II1) = AK(II1,II); AK(II,II) = AK(II,II)+EAL*SS; AK(II1,JJ1) = -EAL*CC; AK(II1,JJ) = -EAL*CS; AK(II,JJ1) = -EAL*CS; AK(II,JJ) = -EAL*SS; AK(JJ1,II1) = AK(II1,JJ1); AK(JJ1,II) = AK(II,JJ1); AK(JJ,II1) = AK(II1,JJ); AK(JJ,II) = AK(II,JJ); AK(JJ1,JJ1) = AK(JJ1,JJ1)+EAL*CC; AK(JJ1,JJ) = AK(JJ1,JJ)+EAL*CS; AK(JJ,JJ1) = AK(JJ1,JJ); AK(JJ,JJ) = AK(JJ,JJ)+EAL*SS; end disp('整体刚度矩阵') disp(AK) %引入边界约束条件 for JE = 1:NR %采用置大数法 %对JE个约束结点 JN = LR(1,JE); %约束结点编号 DX = DD(1,JE); %x方向已知位移 DY = DD(2,JE); %y方向已知位移 %根据约束信息改变相应的刚度系数及荷载项 if DX < 10000 JN1 = 2*JN-1; if dot(AK(JN1,:),AK(JN1,:)) == 0 AK(JN1,JN1) = DP; else AK(JN1,JN1) = DP * AK(JN1,JN1); %刚度主元素自乘以系数DP end if DX == 0 P(JN1) = 0; else P(JN1) = AK(JN1,JN1)*DX; end end if DY < 10000 JN2 = 2*JN; if dot(AK(JN2,:),AK(JN2,:)) == 0 AK(JN2,JN2) = DP; else AK(JN2,JN2) = DP * AK(JN2,JN2); %刚度主元素自乘以系数DP end AK(JN2,JN2) = DP * AK(JN2,JN2); %刚度主元素自乘以系数DP if DY == 0 P(JN2) = 0; else P(JN2) = AK(JN2,JN2)*DY; end end end %求解方程AK * U = P U = AK \ P %左除 disp('结点位移(m)') disp(U) %计算单元内力 for IE = 1:NM NI = JD(1,IE); NJ = JD(2,IE); DIX = U(2*NI-1); DIY = U(2*NI); DJX = U(2*NJ-1); DJY = U(2*NJ); FN(IE) = -EL(IE)*(RF(2,IE)*(DIX-DJX)+RF(1,IE)*(DIY-DJY)); end disp('单元内力(kN)') disp(FN) %计算支座反力 for JE = 1:NR %对第JE约束结点 NI = LR(1,JE); %约束结点的整体编号，并假定为各连接单元的I端 ND = LR(2,JE); %得到与该结点相连接的单元总数 DXI = U(2*NI-1); DYI = U(2*NI); REX = 0; REY = 0; for IE = 1:ND LDE = LM(IE+sum(LR(2,1:JE-1))); %整体编号 NJ = JD(2,LDE); if NJ == NI NJ = JD(1,LDE); end DXJ = U(2*NJ-1); DYJ = U(2*NJ); DX = DXI - DXJ; DY = DYI - DYJ; REX = REX+EL(LDE)*(RF(2,LDE)^2*DX+RF(1,LDE)*RF(2,LDE)*DY); REY = REY+EL(LDE)*(RF(1,LDE)^2*DY+RF(1,LDE)*RF(2,LDE)*DX); end RE(1,JE) = REX-P(2*NI-1); RE(2,JE) = REY-P(2*NI); end disp('支座反力(kN)') disp(RE)