结构力学2大作业

2015结构力学2大作业1

编制原理：

本程序的编制原理是基于矩阵位移法的基本思想，首先把斜拉桥整个复杂的结构体系拆分为单个简单的单元，以结点位移为基本未知量，列出各个单元的单元刚度矩阵，再把单元刚度矩阵组装为整体刚度矩阵，同时根据各个结点所受合力算出每个结点的位移，再根据单元刚度方程算出各个单元的弯矩。

程序说明：

左侧主塔左侧7根拉杆，右侧5跟拉杆，右侧主塔两侧各3根拉杆，把整个体系拆分为39个单元，共24个结点，，采用后处理法，共72个位移编号，对于梁两端结点按照铰支座处理，列出总体刚度矩阵后，用置0置1法处理总刚度矩阵。根据公式k F算出 ，同时在根据 算出弯矩。整体单元编号和结点编号如图：

算法流程：

源代码：

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');

Hitwh130720206;

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+18,1)=sqrt((9*L-j*L)^2+4*H^2/9);

end

for j=9:13

l(j+18,1)=sqrt((j*L-8*L)^2+4*H^2/9);

end

for j=14:16

l(j+18,1)=sqrt((17*L-j*L)^2+4*H^2/9);

end

for j=17:19

l(j+18,1)=sqrt((j*L-16*L)^2+4*H^2/9);

end %钢索长度

for j=20:37

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(72,2);f2=zeros(72,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(19,1)=q*L^2/8;

f1(19,2)=0;

f2(19,1)=5*q*L/8;

f2(19,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(16,1)=q*L^2/3;

f2(16,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:15

f1(j,1)=1/12*q*L^2;

f2(j,1:2)=q*L/2;

end

for j=17:18

f1(j,1)=1/12*q*L^2;

f2(j,1:2)=q*L/2;

end

for j=2:18

f1(j,2)=-f1(j,1);

end

K=zeros(72);

K1=zeros(72);

K2=zeros(72);

K3=zeros(72); %与梁总体刚度矩阵同阶的零阵；

for m=1:19

if m~=8&&m~=13&&m~=16

k(:,:,m)=k1;

else k(:,:,m)=k2;

end

end

for m=1:19

a=zeros(72);

a(3*m-2:3*m+3,3*m-2:3*m+3)=k(:,:,m);

K1=K1+a;

end %梁刚度矩阵组装

e=1;

for j=20:31

e=e+1;

b=zeros(72);

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,61:63)=k0(1:3,4:6,j);

b(61:63,3*e-2:3*e)=k0(4:6,1:3,j);

b(61:63,61:63)=k0(4:6,4:6,j);

K2=K2+b;

end

for j=32:37

e=e+1;

b=zeros(72);

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,64:66)=k0(1:3,4:6,j);

b(64:66,3*e-2:3*e)=k0(4:6,1:3,j);

b(64:66,64:66)=k0(4:6,4:6,j);

K2=K2+b;

end %拉杆总刚度矩阵

kt=Tt'*kt*Tt;

b=zeros(72);

b(61:63,61:63)=kt(1:3,1:3);

b(61:63,67:69)=kt(1:3,4:6);

b(67:69,61:63)=kt(4:6,1:3);

b(67:69,67:69 )=kt(4:6,4:6);

K3=K3+b;

b=zeros(72);

b(64:66,64:66)=kt(1:3,1:3);

b(64:66,70:72)=kt(1:3,4:6);

b(70:72,64:66)=kt(4:6,1:3);

b(70:72,70:72)=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(58:59,:)=0;

K(:,58:59)=0;

K(58,58)=1;

K(59,59)=1;

K(67:72,:)=0;

K(:,67:72)=0;

for j=67:72

K(j,j)=1;

end %置0置1

f=zeros(72,1);

for j=1:19

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('左边主塔塔顶水平位移值:'); …… 此处隐藏：2537字，全部文档内容请下载后查看。喜欢就下载吧 ……