三种动态子结构方法的程序

给出了三种动态子结构方法的程序，对比了自由界面模态综合法，固定界面模态综合法，双协调自由界面模态综合法三种方法。-Shows three dynamic substructure method of procedure, compared to the free interface modal synthesis method, fixed-interface modal synthesis method, two-coordinate free interface modal synthesis method of three methods.

%%%%%%%%%%%%%双协调自由界面
clear;
b1=10;h1=10;b2=0.7;h2=0.7;
l1=50;A1=100;M1=4e5;EI1=1220;den1=M1/l1;I1=b1*（h1^3）/12;E1=EI1/I1;E1=7e10;
n1=50;%the number of the segments in length;
le1=l1/n1;
l2=5;A2=0.49;M2=3.1e4;EI2=218;den2=M2/l2;I2=b2*（h2^3）/12;E2=EI2/I2;E2=7e10;
n2=5;
le2 = l2 / n2;
%寻找单元的节点号
for i = 1:n1
    Element1（i:） = [i 1 + i];
end
for i = 1:n2
    Element2（i:） = [i 1 + i];
end
%单元节点的坐标号
for i = 1:n1
    for j = 1:2
        U1（（（j-1） * 3 + 1）:（j * 3）i） = [（Element1（ij） * 3 -2）;（Element1（ij） * 3 -1）;（Element1（ij） * 3）];
    end
end
for i = 1:n2
    for j = 1:2
        U2（（（j-1） * 3 + 1）:（j * 3）i） = [（Element2（ij） * 3 -2）;（Element2（ij） * 3 -1）;（Element2（ij） * 3）];
    end
end
Node1 = 3*（n1 + 1）;%总的位移坐标数
Node2 = 3*（n2 + 1）;
%边界点号
BN1 = [2 5]; BN2=[2 5];

syms x y z
X1 = 1 - （x/le1）;X2=x/le1;
g1=X1;g2=1-3*X2*X2+2*（X2^3）;g3=le1*（X2-2*（X2^2）+X2^3）;g4=X2;g5=3*（X2^2）-2*（X2^3）;g6=le1*（-X2^2+X2^3）;
N1=[g1 0 0 g4 0 0];N2=[0 g2 g3 0 g5 g6];N=[N1;N2];
NN1 = transpose（N） * N;

%单元质量阵
intx = int（NN1 x0 le1）;
Me1 = double（den1  * intx）;Me2 = double（den2  * intx）;


B1=diff（N1x）;B2=-y*diff（N2x2）;B=[B1;B2];
BB1=B‘*E1*B;
intx=int（BB1x0le1）;
inty=int（intxy-h1/2h1/2）;
intz=int（intyz0b1）;
Ke1 = double（intz）;
BB2=B‘*E2*B;
intx=int（BB2x0le2）;
inty=int（intxy-h2/2h2/2）;
intz=int（intyz0b2）;
Ke2 = double（intz）;



%组集（同一根梁内的组集）
MM1 = zeros（Node1Node1）;MM2 = zeros（Node2Node2）;
KK1 = zeros（Node1Node1）;KK2 = zeros（Node2Node2）;


for k = 1:n1
    for i = 1:6
        for j = 1:6
            MM1（U1（ik）U1（jk）） = MM1（U1（ik）U1（jk）） + Me1（ij）;
            KK1（U1（ik）U1（jk）） = KK1（U1（ik）U1（jk）） + Ke1（ij）;
        end
    end
end

for k = 1:n2
    for i = 1:6
        for j = 1:6
            MM2（U2（ik）U2（jk）） = MM2（U2（ik）U2（jk）） + Me2（ij）;
            KK2（U2（ik）U2（jk）） = KK2（U2（ik）U2（jk）） + Ke2（ij）;
        end
    end
end

M1 = MM1;M2=MM2;
K1 = KK1;K2=KK2;
e1=eye（length（MM1））;
e1=[e1;e1（4:5:）;e1（13:14:）];
e1（13:14:）=[];e1（4:5:）=[];
MM1=e1*MM1*inv（e1）;KK1=e1*KK1*inv（e1）;
M1 = MM1;
K1 = KK1;%153行列
B=zeros（length（MM1）1）;B（end-1）=1;B=e1*B;
C1=zeros（1length（MM1））;C1（end-1）=1;C1=C1*inv（e1）;
C2=zeros（1length（MM1））;


e2=eye（length（MM2））;
e2=[e2;e2（4:5:）;e2（13:14:）];
e2（13:14:）=[];e2（4:5:）=[];
MM2=e2*MM2*inv（e2）;KK2=e2*KK2*inv（e2）;
M2=MM2;
K2=KK2;%18行列


%计算频率和振型
disp（‘正在计算特征值……‘）;
[V1D1] = eig（K1M1）;%K1*V1=M1*V1*D1
[V2D2] = eig（K2M2）;
w1=eig（K1M1）.^（1/2）;w2=eig（K2M2）.^（1/2）;%所有特征值组成的向量
u1=w1/6.28;u2=w2/6.28;
%B=V1（:1:n1+1）‘*B;
%C1=C1*V1（:1:n1+1）;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%求剩余模态
nn1=length（M1）;nn2=length（K2）;%nn1=153   nn2=18
Bo1=zeros（nn1）;
Bo1（:5:nn1）=[];
Bo1（nn1-31）=1;
Bo1（nn1-22）=1;
Bo1（nn1-13）=1;
Bo1（nn14）=1;
Fk1=V1（:1:n1+1）;
Fd1=V1（:n1+2:nn1）;
Kd1=Fd1‘*K1*Fd1;
pasi1=Fd1*inv（Kd1）*Fd1‘*Bo1;
Bo2=zeros（nn2）;
Bo2（:5:nn2）=[];
Bo2（nn2-

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       4546  2011-07-26 14:05  ziyoujiemian.m

     文件       5421  2011-07-26 16:21  gaijinjiangjie.m

     文件       5631  2011-07-26 15:13  gudingjiemian.m

----------- ---------  ---------- -----  ----

                15598                    3