2012结构力学2大作业1.doc

2012结构力学2大作业1
Ø 编制原理：
本程序的编制原理是基于矩阵位移法的基本思想，首先把斜拉桥整个复杂的结构体系拆分为单个简单的单元，以结点位移为基本未知量，列出各个单元的单元刚度矩阵，再把单元刚度矩阵组装为整体刚度矩阵，同时根据各个结点所受合力算出每个结点的位移，再根据单元刚度方程算出各个单元的弯矩。
Ø 程序说明：
27
左侧主塔左侧7根拉杆，右侧5跟拉杆，右侧主塔两侧各6根拉杆，把整个体系拆分为51个单元，共30个结点，，采用后处理法，共90个位移编号，对于梁两端结点按照铰支座处理，列出总体刚度矩阵后，用置0置1法处理总刚度矩阵。根据公式算出，同时在根据算出弯矩。整体单元编号和结点编号如图：
L

此处无搭接

K=5

1
1
n=6

2H/3
H/3

L
L
L
2L
L
L
L

1
26
28
29
30

位移编号为：
结点编号
位移编号
1
（1,2,3）
2
(1,2,3)
……
……
30
(88,89,90)

结点、单元、位移编号
结构离散化
Ø 算法流程：

将梁、拉杆、主塔的单元刚度矩阵转换为整体坐标下的单元刚度矩阵
单元分析，建立拉杆、梁、主塔的单元刚度矩阵

确定结构中各个结点荷载的矩阵
对单元刚度矩阵进行组装，形成整体刚度矩阵

利用单元刚度矩阵和所求位移求解内力
建立刚度矩阵方程，求解位移

Ø 源代码：
function F=zy(H,L)
H=input('塔高：H=');
L=input('节点间距：L=');
EIc=input('主塔抗弯刚度EIc=');
EAc=input('主塔抗压刚度EAc=');
EIb=input('主梁抗弯刚度EIb=');
EAb=input('主梁抗压刚度EAb=');
EA0=input('拉杆抗压刚度EA0=');
q=input('均布荷载q='); %输入数据
if H<=0||L<=0||EIc<=0||EAc<=0||EIb<=0||EAb<=0
clear;
fprintf('输入错误，请重新输入！\n');
hitwh101210232;
else
k1 = [ EAb/L 0 0 -EAb/L 0 0
0 12*EIb/L^3 6*EIb/L^2 0 -12*EIb/L^3 6*EIb/L^2
0 6*EIb/L^2 4*EIb/L 0 -6*EIb/L^2 2*EIb/L
-EAb/L 0 0 EAb/L 0 0
0 -12*EIb/L^3 -6*EIb/L^2 0 12*EIb/L^3 -6*EIb/L^2
0 6*EIb/L^2 2*EIb/L 0 -6*EIb/L^2 4*EIb/L] ;%长L的梁单元刚度矩阵
k2= [ 1/2*EAb/L 0 0 -1/2*EAb/L 0 0
0 3/2*EIb/L^3 3/2*EIb/L^2 0 -3/2*EIb/L^3 3/2*EIb/L^2
0 3/2*EIb/L^2 2*EIb/L 0 -3/2*EIb/L^2 EIb/L
-1/2*EAb/L 0 0 1/2*EAb/L 0 0
0 -3/2*EIb/L^3 -3/2*EIb/L^2 0 3/2*EIb/L^3 -3/2*EIb/L^2
0 3/2*EIb/L^2 EIb/L 0 -3/2*EIb/L^2 2*EIb/L];%长2L梁的单元刚度矩阵

kt = [ EAc/H 0 0 -EAc/H 0 0
0 12*EIc/H^3 6*EIc/H^2 0 -12*EIc/H^3 6*EIc/H^2
0 6*EIc/H^2 4*EIc/H 0 -6*EIc/H^2 2*EIc/H
-EAc/H 0 0 EAc/H 0 0
0 -12*EIc/H^3 -6*EIc/H^2 0 12*EIc/H^3 -6*EIc/H^2
0 6*EIc/H^2 2*EIc/H 0 -6*EIc/H^2 4*EIc/H] ;%主塔单元刚度矩阵

for j=2:8
l(j+24,1)=sqrt((9*L-j*L)^2+4*H^2/9);
end
for j=9:13
l(j+24,1)=sqrt((j*L-8*L)^2+4*H^2/9);
end
for j=14:19
l(j+24,1)=sqrt((20*L-j*L)^2+4*H^2/9);
end
for j=20:25
l(j+24,1)=sqrt((j*L-19*L)^2+4*H^2/9);
end %钢索长度
for j=26:49
k0(:,:,j)=[1 0 0 -1 0 0
0 0 0 0 0 0
0 0 0 0 0 0
-1 0 0 1 0 0
0 0 0 0 0 0
0 0 0 0 0 0]*EA0/l(j,1); %钢索的单元刚度矩阵，L（j,1)为钢索的长度；
x(j,1)=asin(2*H/3/l(j,1)); %钢索与梁的夹角
T0(:,:,j)= [cos(x(j,1)) sin(x(j,1)) 0 0 0 0
-sin(x(j,1)) cos(x(j,1)) 0 0 0 0
0 0 1 0 0 0
0 0 0 cos(x(j,1)) sin(x(j,1)) 0
0 0 0 -sin(x(j,1)) cos(x(j,1)) 0
0 0 0 0 0 1] ; %钢索的坐标转换阵
end
Tt=[0 -1 0 0 0 0
1 0 0 0 0 0
0 0 1 0 0 0
0 0 0 0 -1 0
0 0 0 1 0 0
0 0 0 0 0 1];%主塔的坐标转换阵
f1=zeros(90,2);f2=zeros(90,2);
f1(1,1)=0; %杆端弯矩
f1(1,2)=-q*L^2/8;
f2(1,1)=3*q*L/8;%杆端剪力
f2(1,2)=5*q*L/8;
f1(25,1)=q*L^2/8;
f1(25,2)=0;
f2(25,1)=5*q*L/8;
f2(25,2)=3*q*L/8;
f1(8,1)=q*L^2/3;
f2(8,1:2)=q*L;
f1(13,1)=q*L^2/3;
f2(13,1:2)=q*L;
f1(19,1)=q*L^2/3;
f2(19,1:2)=q*L;
for j=2:7
f1(j,1)=1/12*q*L^2;
f2(j,1:2)=q*L/2;
end
for j=9:12
f1(j,1)=1/12*q*L^2;
f2(j,1:2)=q*L/2;
end
for j=14:18
f1(j,1)=1/12*q*L^2;
f2(j,1:2)=q*L/2;
end
for j=20:24
f1(j,1)=1/12*q*L^2;
f2(j,1:2)=q*L/2;
end
for j=2:24
f1(j,2)=-f1(j,1);

end
K=zeros(90);
K1=zeros(90);
K2=zeros(90);
K3=zeros(90); %与梁总体刚度矩阵同阶的零阵；
for m=1:25
if m~=8&&m~=13&&m~=19
k(:,:,m)=k1;
else k(:,:,m)=k2;
end
end
for m=1:25
a=zeros(90);
a(3*m-2:3*m+3,3*m-2:3*m+3)=k(:,:,m);
K1=K1+a;
end %梁刚度矩阵组装
e=1;
for j=26:37
e=e+1;
b=zeros(90);
k0(:,:,j)=T0(:,:,j)'*k0(:,:,j)*T0(:,:,j);
b(3*e-2:3*e,3*e-2:3*e)=k0(1:3,1:3,j);
b(3*e-2:3*e,79:81)=k0(1:3,4:6,j);
b(79:81,3*e-2:3*e)=k0(4:6,1:3,j);
b(79:81,79:81)=k0(4:6,4:6,j);
K2=K2+b;
end
for j=38:49
e=e+1;
b=zeros(90);
k0(:,:,j)=T0(:,:,j)'*k0(:,:,j)*T0(:,:,j);
b(3*e-2:3*e,3*e-2:3*e)=k0(1:3,1:3,j);
b(3*e-2:3*e,82:84)=k0(1:3,4:6,j);
b(82:84,3*e-2:3*e)=k0(4:6,1:3,j);
b(82:84,82:84)=k0(4:6,4:6,j);
K2=K2+b;
end %拉杆总刚度矩阵
kt=Tt'*kt*Tt;
b=zeros(90);
b(79:81,79:81)=kt(1:3,1:3);
b(79:81,85:87)=kt(1:3,4:6);
b(85:87,79:81)=kt(4:6,1:3);
b(85:87,85:87)=kt(4:6,4:6);
K3=K3+b;
b=zeros(90);
b(82:84,82:84)=kt(1:3,1:3);
b(82:84,88:90)=kt(1:3,4:6);
b(88:90,82:84)=kt(4:6,1:3);
b(88:90,88:90)=kt(4:6,4:6);
K3=K3+b; %主塔刚度矩阵
K=K1+K2+K3; %总体刚度矩阵
K(1:2,:)=0;
K(:,1:2)=0;
K(1,1)=1;
K(2,2)=1;
K(76:77,:)=0;
K(:,76:77)=0;
K(76,76)=1;
K(77,77)=1;
K(85:90,:)=0;
K(:,85:90)=0;
for j=85:90
K(j,j)=1;
end %置0置1
f=zeros(90,1);
for j=1:25
f(3*j,1)=-(f1(j,2)+f1(j+1,1));
f(3*j-1,1)=-(f2(j,2)+f2(j+1,1));
end %f为单元等效结点荷载矩阵
v=inv(K)*f;
disp('左边主塔塔顶水平位移值:');
disp(v(79,1));
z(1:3,1)=v(79:81,1);
z(4:6,1)=v(85:87,1);%z主塔位移矩阵
M=kt*z; %M为主塔结点力
disp('左边主塔塔底弯矩值:');
disp(M(6,1));
end
Ø 试算算例输入数据及计算结果
>> hitwh101210232
塔高：H=1
节点间距：L=0
主塔抗弯刚度EIc=1
主塔抗压刚度EAc=1
主梁抗弯刚度EIb=1
主梁抗压刚度EAb=1
拉杆抗压刚度EA0=1
均布荷载q=1
输入错误，请重新输入！
塔高：H=1
节点间距：L=1
主塔抗弯刚度EIc=1
主塔抗压刚度EAc=0
主梁抗弯刚度EIb=1
主梁抗压刚度EAb=1
拉杆抗压刚度EA0=1
均布荷载q=1
输入错误，请重新输入！
塔高：H=1
节点间距：L=1
主塔抗弯刚度EIc=1
主塔抗压刚度EAc=1
主梁抗弯刚度EIb=1
主梁抗压刚度EAb=1
拉杆抗压刚度EA0=0
均布荷载q=1
左边主塔塔顶水平位移值:
0

左边主塔塔底弯矩值:
0

>> hitwh101210232
塔高：H=1
节点间距：L=1
主塔抗弯刚度EIc=1
主塔抗压刚度EAc=1
主梁抗弯刚度EIb=1
主梁抗压刚度EAb=1
拉杆抗压刚度EA0=1
均布荷载q=1
左边主塔塔顶水平位移值:
-1.5748

左边主塔塔底弯矩值:
-4.7243

K(2,2)=1;
K(76:77,:)=0;
K(:,76:77)=0;
K(76,76)=1;
K(77,77)=1;
K(85:90,:)=0;
K(:,85:90)=0;
for j=85:90
K(j,j)=1;
end %置0置1
f=zeros(90,1);
for j=1:25
f(3*j,1)=-(f1(j,2)+f1(j+1,1));
f(3*j-1,1)=-(f2(j,2)+f2(j+1,1));
end %f为单元等效结点荷载矩阵
v=inv(K)*f;
disp('左边主塔塔顶水平位移值:');
disp(v(79,1));
z(1:3,1)=v(79:81,1);
z(4:6,1)=v(85:87,1);%z主塔位移矩阵
M=kt*z; %M为主塔结点力
disp('左边主塔塔底弯矩值:');
disp(M(6,1));
end
Ø 试算算例输入数据及计算结果
>> hitwh101210232
塔高：H=1
节点间距：L=0
主塔抗弯刚度EIc=1
主塔抗压刚度EAc=1
主梁抗弯刚度EIb=1
主梁抗压刚度EAb=1
拉杆抗压刚度EA0=1
均布荷载q=1
输入错误，请重新输入！
塔高：H=1
节点间距：L=1
主塔抗弯刚度EIc=1
主塔抗压刚度EAc=0
主梁抗弯刚度EIb=1
主梁抗压刚度EAb=1
拉杆抗压刚度EA0=1
均布荷载q=1
输入错误，请重新输入！
塔高：H=1
节点间距：L=1
主塔抗弯刚度EIc=1
主塔抗压刚度EAc=1
主梁抗弯刚度EIb=1
主梁抗压刚度EAb=1
拉杆抗压刚度EA0=0
均布荷载q=1
左边主塔塔顶水平位移值:
0

左边主塔塔底弯矩值:
0

>> hitwh101210232
塔高：H=1
节点间距：L=1
主塔抗弯刚度EIc=1
主塔抗压刚度EAc=1
主梁抗弯刚度EIb=1
主梁抗压刚度EAb=1
拉杆抗压刚度EA0=1
均布荷载q=1
左边主塔塔顶水平位移值:
-1.5748

左边主塔塔底弯矩值:
-4.7243

9