%平面桁架静力计算程序
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)
评论0