matlab-矩阵位移法编程-结构力学. doc
矩阵位移法编程大作业(091210211 )一、编制原理本程序原理根据结构力学矩阵位移法原理，以结构节点位移为基本未知量，将待分析结构分解为已知节点力节点力位移关系的单跨梁集，保留在结构中的基本未知位移二、程序说明本程序是计算3层11框架之间右侧节点位移和弯矩的程序，编译过程按照矩阵位移法的预处理法进行。 首先，以结构杆件交点为节点，共有36个节点和108个位移编号，根据梁、柱、斜杆的不同分别建立单元刚度矩阵，转换为整体坐标系下的刚度矩阵，然后将所有杆件的单元刚度矩阵合并为整体刚度矩阵然后，通过分析载荷，自己决定载荷矩阵，直接写入程序中。 这样，可以求出36个节点的108个位移，用各单元的单元刚度矩阵和得到的位移求出单元杆的内力。 离散化编号如下图3所示，算法流程采用杆单元刚度矩阵和求解位移求解内力方程，求解位移确定综合节点荷载矩阵预处理法将各杆单元刚度矩阵合并为整体刚度矩阵，得到梁、柱、 斜杆局部坐标系下单元刚度矩阵确定梁、柱、斜杆整体坐标系下刚度矩阵单元分析结构离散化编号输入输出结果4、源代码结构力学大作业3楼11框架间矩阵位移法编程的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； 012 * eib/l * l6 * eib/l * l0-12 * eib/l * l6 * eib/l * l； 0 * eib/l * l4 * eib/l0-6 * eib/l * L2 * eib/l； -EAb/L 0 0 EAb/L 0 0； 0-12 * eib/l * l-6 * eib/l * l012 * eib/l * l-6 * eib/l * l； 0 * eib/l * L2 * eib/l0-6 * eib/l * l4 * eib/l； 柱的单元刚度矩阵kc0EAc/h 0 0 -EAc/h 0 0； 012 * EIC/h * h * h6 * EIC/h * h0-12 * EIC/h * h6 * EIC/h * h； 0 * EIC/h * H4 * EIC/h0-6 * EIC/h * H2 * EIC/h； -EAc/h 0 0 EAc/h 0 0； 0-12 * EIC/h * h * h-6 * EIC/h * h012 * EIC/h * h * h-6 * EIC/h * h； 0 * EIC/h * H2 * EIC/h0-6 * EIC/h * H4 * EIC/h斜杆的机组刚度矩阵Hsqrth*hL*L； o0EaO/h0-EaO/h0； 012 * eio/h * h * h6 * eio/h * h0-12 * eio/h * h6 * eio/h * h； 0 * eio/h * H4 * eio/h0-6 * eio/h * H2 * eio/h； -EAo/H 0 0 EAo/H 0 0； 0-12 * eio/h * h * h-6 * eio/h * h012 * eio/h * h * h-6 * eio/h * h； 0 * eio/h * H2 * eio/h0-6 * eio/h * H4 * eio/h； 国安kbT1*kb0*T1； 整体坐标下梁的单元刚度矩阵kcT2*kc0*T2； 整体坐标下柱的单元刚性矩阵koT3*ko0*T3； 整体坐标斜棒单元刚性矩阵xZeros108、108； y zeros 108、108； ZZ eros 108、108； 定义108阶0矩阵k1Zeros108、108； k2 zeros 108，108； k3 zeros 108，108； k4zeros 108、108； k5zeros 108、108； K6 zeros 108，108； k7zeros 108、108； k8 zeros 108，108； K9 zeros 108，108； 将梁杆单元矩阵合并为整体刚性矩阵的循环语句forii 111 x3* ii-23 * ii3，3 * ii-23 * ii3kb； k1x； x zeros 108、108； end for ii 1323 y3 * ii-23 * ii3，3 * ii-23 * ii3kb； k1y； y zeros 108、108； end for ii 2535 z3 * ii-23 * ii3，3 * ii-23 * ii3kb； K1K1Z； ZZ eros 108、108； 将end柱棒单元矩阵合并为整体刚性矩阵的循环语句for jj1
36 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矩阵位移法编程大作业asjdjzg 091210211