continousbeam.m

结构力学连续梁节点弯矩求解

clc,clear all;
input('   **连续梁内力计算程序**');
NE=input('请输入单元数目：');
N=input('请输入结点转角未知量数目：');
JD=input('请输入单元定位向量矩阵：');
EI=input('请输入各单元抗弯系数数组：');
BL=input('请输入各单元长度数组：');
PJ=input('请输入结点荷载数组：');
FC=input('请输入集中荷载数组：');
FD=input('请输入均布荷载数组：');     %以上为用户通过输入为各组数据赋值

fprintf('输出原始数据');
fprintf('单元数目：');
NE
fprintf('结点转角未知量数目：');
N
fprintf('单元定位向量矩阵：');
JD
fprintf('单元长度数组：');
BL 
fprintf('单元截面刚度数组：');
EI
fprintf('单元集中力数组：');
FC
fprintf('单元均布荷载数组：');
FD
fprintf('单元结点荷载数组：');
PJ
fprintf('单元刚度系数数组');
EI                                 %以上为输出用户已输入的数据

FM=zeros(2,NE);                    %建立单元固端力矩矩阵
for i=1:N
        FM(:,i)=[-FC(i)*BL(i)/8-FD(i)*BL(i)^2/12;FC(i)*BL(i)/8+FD(i)*BL(i)^2/12];
end                                %计算单元固端力矩
P=zeros(N,1);
for i=1:NE
    for j=1:2
        if JD(j,i)>0
            P(JD(j,i),1)=P(JD(j,i),1)-FM(j,i);
        end
    end
end                                %计算等效结点荷载
F=zeros(N,1);
PJ=PJ';
F=PJ+P;                            %计算总结点荷载

KE=zeros(2,2,NE);
for i=1:NE
    KE(:,:,i)=[4*EI(i),2*EI(i);2*EI(i),4*EI(i)];
end                                %形成单元刚度矩阵
K=zeros(N,N);
for i=1:NE
    for m=1:2
        for n=1:2
            if JD(m,i)~=0&&JD(n,i)~=0
                K(JD(m,i),JD(n,i))= K(JD(m,i),JD(n,i))+KE(m,n,i);
            end
        end
    end
end                                %由单元刚度矩阵形成整体刚度矩阵
P=K\F;                             %计算结点转角位移

for i=1:NE
    for j=1:2
        if JD(j,i)~=0
            C(j,i)=P(JD(j,i));
        end
    end
end                                %计算各单元杆端转角位移
FJ=zeros(2,NE);
for i=1:NE
    FJ(:,i)=KE(:,:,i)*C(:,i)+FM(:,i);                               
end                                %计算单元杆端弯矩
input('结点位移：');
P
input('单元杆端弯矩');
FJ