matlab-矩阵位移法编程-结构力学.doc
矩阵位移法编程大作业(091210211)一、编制原理本程序的原理是基于结构力学矩阵位移法原理，以结构结点位移作基本未知量，将要分析的结构拆成已知节点力结点力位移关系的单跨梁集合，通过强令结构发生待定的基本未知位移，在各个单跨梁受力分析结果的基础上通过保证结构平衡建立位移法的线性方程组，从而求得基本未知量。二、程序说明本程序是计算3层11跨框架右侧结点的位移和弯矩的程序，编译过程是按照矩阵位移法的先处理法进行的。首先将结构杆件的交汇点作为结点，共有36个结点和108个位移编号，然后根据梁、柱、斜杆的不同分别建立单元刚度矩阵，然后转换为整体坐标系下的刚度矩阵，然后将所有杆件的单元刚度矩阵整合成为总体刚度矩阵，在进行整合时连续运用for函数，最终形成108阶的总体刚度矩阵。然后通过对荷载的分析自己确定出荷载矩阵，直接写进程序。这样就可以把36个结点的108个位移求得，然后再利用各个单元的单元刚度矩阵和所得的位移求得单元杆件的内力。离散化编号如下图三、算法流程利用杆件单元刚度矩阵和所求位移求内力建立方程，求解位移确定综合结点荷载矩阵先处理法把各个杆件的单元刚度矩阵整合成总体刚度矩阵建立梁、柱、斜杆的局部坐标系下的单元刚度矩阵确定梁、柱、斜杆在整体坐标系下的刚度矩阵单元分析结构离散化编号结束输出结果四、源代码结构力学大作业 3层11跨框架矩阵位移法编程 王贝 091210211h输入单层高h;L输入单跨度L;EIc输入柱子的抗弯刚度EIc;EAc输入柱子的抗压刚度EAc;EIb输入梁的抗弯刚度EIb;EAb输入梁的抗压刚度EAb;EIo输入斜杆的抗弯刚度EIo;EAo输入斜杆的抗压刚度EAo;q输入侧向均布荷载集度q;T11,0,0,0,0,0; 0,1,0,0,0,0; 0,0,1,0,0,0; 0,0,0,1,0,0; 0,0,0,0,1,0; 0,0,0,0,0,1;角度为0的转换矩阵T20,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;角度为90的转换矩阵xatanh/L;Tcosx,sinx,0,0,0,0; -sinx,cosx,0,0,0,0; 0,0,1,0,0,0; 0,0,0,cosx,sinx,0; 0,0,0,-sinx,cosx,0; 0,0,0,0,0,1;斜杆的转换矩阵T3T;梁的单元刚度矩阵kb0EAb/L 0 0 -EAb/L 0 0; 0 12*EIb/L*L*L 6*EIb/L*L 0 -12*EIb/L*L*L 6*EIb/L*L; 0 6*EIb/L*L 4*EIb/L 0 -6*EIb/L*L 2*EIb/L; -EAb/L 0 0 EAb/L 0 0; 0 -12*EIb/L*L*L -6*EIb/L*L 0 12*EIb/L*L*L -6*EIb/L*L; 0 6*EIb/L*L 2*EIb/L 0 -6*EIb/L*L 4*EIb/L;柱子的单元刚度矩阵kc0EAc/h 0 0 -EAc/h 0 0; 0 12*EIc/h*h*h 6*EIc/h*h 0 -12*EIc/h*h*h 6*EIc/h*h; 0 6*EIc/h*h 4*EIc/h 0 -6*EIc/h*h 2*EIc/h; -EAc/h 0 0 EAc/h 0 0; 0 -12*EIc/h*h*h -6*EIc/h*h 0 12*EIc/h*h*h -6*EIc/h*h; 0 6*EIc/h*h 2*EIc/h 0 -6*EIc/h*h 4*EIc/h斜杆的单元刚度矩阵Hsqrth*hL*L;ko0EAo/H 0 0 -EAo/H 0 0; 0 12*EIo/H*H*H 6*EIo/H*H 0 -12*EIo/H*H*H 6*EIo/H*H; 0 6*EIo/H*H 4*EIo/H 0 -6*EIo/H*H 2*EIo/H; -EAo/H 0 0 EAo/H 0 0; 0 -12*EIo/H*H*H -6*EIo/H*H 0 12*EIo/H*H*H -6*EIo/H*H; 0 6*EIo/H*H 2*EIo/H 0 -6*EIo/H*H 4*EIo/H;kbT1*kb0*T1;总体坐标下梁的单元刚度矩阵kcT2*kc0*T2;总体坐标下柱子的单元刚度矩阵koT3*ko0*T3;总体坐标斜杆的单元刚度矩阵Xzeros108,108;Yzeros108,108;Zzeros108,108;定义108阶0矩阵K1zeros108,108;K2zeros108,108;K3zeros108,108;K4zeros108,108;K5zeros108,108;K6zeros108,108;K7zeros108,108;K8zeros108,108;K9zeros108,108;把梁杆单元矩阵整合到总体刚度矩阵的循环语句for ii111 X3*ii-23*ii3,3*ii-23*ii3kb; K1K1X;Xzeros108,108;endfor ii1323 Y3*ii-23*ii3,3*ii-23*ii3kb; K1K1Y;Yzeros108,108;endfor ii2535 Z3*ii-23*ii3,3*ii-23*ii3kb; K1K1Z;Zzeros108,108;end把柱杆单元矩阵整合到总体刚度矩阵的循环语句for jj136 K23*jj-23*jj,3*jj-23*jjkc46,46;endfor jj124 K33*jj-23*jj,3*jj-23*jjkc13,13;endfor jj124 K43*jj-23*jj,3*jj343*jj36kc13,46;endfor jj124 K53*jj343*jj36,3*jj-23*jjkc46,13;end把斜杆杆单元矩阵整合到总体刚度矩阵的循环语句for gg31227 K63*gg-23*gg,3*gg-23*ggko46,46;endfor gg21214 K73*gg-23*gg,3*gg-23*ggko13,13;endfor gg21214 K83*gg-23*gg,3*gg373*gg39ko13,46;endfor gg21214 K93*gg373*gg39,3*gg-23*ggko46,13;endKK1K2K3K4K5K6K7K8K9;总体刚度矩阵Pzeros108,1;P1,1h*q;P37,1h*q;P73,1h*q/2;P75,1q*h*h/12;AKP;结构位移B1kb*A103108,1;B2kb*A6772,1;B3kb*A3136,1;D1zeros6,1;D113,1A7072,1;D146,1A106108,1;D2zeros6,1;D213,1A3436,1;D246,1A7072,1;C1kc*D1;C2kc*D2;C3kc46,46*A3436,1;M11,1B16,1;M21,1C16,1;M31,1B26,1;M41,1C26,1;M51,1C13,1;M61,1B36,1;M71,1C33,1;M81,1C23,1; for i13 m36*i-2;fprintf第d层最右侧节点的位移是dn,i,Am,1endfprintf第1层最右侧节点左侧杆的弯矩是fn,M6fprintf第1层最右侧节点下侧杆的弯矩是fn,M7fprintf第1层最右侧节点上侧杆的弯矩是fn,M8fprintf第2层最右侧节点左侧杆的弯矩是fn,M3fprintf第2层最右侧节点下侧杆的弯矩是fn,M4fprintf第2层最右侧节点上侧杆的弯矩是fn,M5fprintf第3层最右侧节点左侧杆的弯矩是fn,M1fprintf第3层最右侧节点下侧杆的弯矩是fn,M2五、试算算例输入数据输入单层高h1输入单跨度L1输入柱子的抗弯刚度EIc1输入柱子的抗压刚度EAc1输入梁的抗弯刚度EIb1输入梁的抗压刚度EAb1输入斜杆的抗弯刚度EIo1输入斜杆的抗压刚度EAo1输入侧向均布荷载集度q1计算结果第1层最右侧节点的位移是-6.219850e-003第2层最右侧节点的位移是-2.152659e-002第3层最右侧节点的位移是-4.131873e-002第1层最右侧节点左侧杆的弯矩是0.000729第1层最右侧节点下侧杆的弯矩是0.008642第1层最右侧节点上侧杆的弯矩是-0.009371第2层最右侧节点左侧杆的弯矩是0.000074第2层最右侧节点下侧杆的弯矩是0.004158第2层最右侧节点上侧杆的弯矩是-0.004232第3层最右侧节点左侧杆的弯矩是-0.000761第3层最右侧节点下侧杆的弯矩是0.000761矩阵位移法编程大作业王贝 091210211
