平面桁架有限元分析matlab程序

平面桁架有限元分析matlab程序，分别使用乘大数法和对角元置1法对刚度矩阵进行处理，可以使用该程序直接完成普通桁架的有限元分析。支持输入相关参数，即可得到相应的分析结果。

%-----------------------------------------------------%
%该程序为平面桁架结构分析有限元MATLAB程序
%不能用于计算含倾斜支座的桁架结构
%输入的信息有节点数、节点信息、单元数、单元信息
%边界条件处理采用的方法是：对角元素乘大数法
%返回的信息主要有：
%   （1）每个单元的刚度矩阵
%   （2）整体刚度矩阵
%   （3）节点位移向量
%   （4）节点载荷向量
%   （5）每个单元的位移向量
%   （6）每个单元的内力向量
%   （7）每个单元的应力
%
%ver:2017/1/6
%-----------------------------------------------------%

clear
clc

np=3;%输入节点数
ne=3;%输入单元数

%node（节点信息）：节点坐标，约束情况，节点载荷
%矩阵规模：np x 6
node=[
     0 0 1 1 0 0;
    4 0 0 1 0 0;
    2 3 0 0 5000 -10000;
];

%element（单元信息）：左右节点ij，弹性模量E，截面积A
%矩阵规模：ne x 4
element=[
    1 2 210e9 0.0001;
    1 3 210e9 0.0001;
    2 3 210e9 0.0001;
]; 

kk=zeros（2*np2*np）;%分配整体刚度矩阵内存空间
f =zeros（2*np1）;%分配载荷内存空间
Elementinfo=zeros（ne3）;%分配空间存储单元sincosEA
elementstiffnesssave=zeros（44ne）;%分配空间存储每个单元的刚度矩阵

%首先生成单元刚度矩阵
%然后集成整体刚度矩阵
for lop=1:ne
    i=element（lop1）;
    j=element（lop2）;
    xi=node（i1）;
    yi=node（i2）;
    xj=node（j1）;
    yj=node（j2）;
    L=sqrt（（xj-xi）*（xj-xi）+（yj-yi）*（yj-yi））;
    s=（yj-yi）/L;
    c=（xj-xi）/L;
    ea=element（lop3）*element（lop4）/L;
    ek = ea*[c*cc*s-c*c-c*s;s*cs*s-s*c-s*s;-c*c-c*sc*cc*s;-c*s-s*ss*cs*s];
    ek    %输出每个单元的刚度矩阵
    elementstiffnesssave（::lop）=ek;
    kk（2*i-1:2*i 2*i-1:2*i）=kk（2*i-1:2*i 2*i-1:2*i）+ek（1:2 1:2）;
    kk（2*i-1:2*i 2*j-1:2*j）=kk（2*i-1:2*i 2*j-1:2*j）+ek（1:2 3:4）;
    kk（2*j-1:2*j 2*i-1:2*i）=kk（2*j-1:2*j 2*i-1:2*i）+ek（3:4 1:2）;
    kk（2*j-1:2*j 2*j-1:2*j）=kk（2*j-1:2*j 2*j-1:2*j）+ek（3:4 3:4）;
    kk   %输出集成每个单元刚度矩阵后的整体刚度矩阵
end

StiffnessSave=kk;%保存整体刚度矩阵

%边界条件处理，采用对角元乘大数法
for lop=1:np
    f（2*lop-11）=f（2*lop-11） + node（lop5）;
    f（2*lop  1）=f（2*lop  1） + node（lop6）;
    if node（lop3） >= 1
        kk（2*lop-1 2*lop-1） = 1e8*kk（2*lop-12*lop-1）;
    end
    if node（lop4） >= 1
        kk（2*lop 2*lop） = 1e8*kk（2*lop 2*lop）;
    end
end

u=kk\f;%高斯消元法求节点位移向量
Force = StiffnessSave * u;%求节点支反力

Elementforce=zeros（ne1）;%存储每个杆单元的内力
Elementstress=zeros（ne1）;%存储每个杆单元的应力
for lop = 1:ne
    i=element（lop1）;
    j=element（lop2）;
    E=element（lop3）;
    A=element（lop4）;
    xi=node（i1）;
    yi=node（i2）;
    xj=node（j1）;
    yj=node（j2）;
    L=sqrt（（xj-xi）*（xj-xi）+（yj-yi）*（yj-yi））;
    s=（yj-yi）/L;
    c=（xj-xi）/L;
    %求单元内力
    Elementforce（lop 1） =E*A/L*（（u（2*j-1）-u（2*i-1））*c+（u（2*j）-u（2*i））*s）; 
    %求单元应力
    Elementstress（lop 1） =E/L*（（u（2*j-1）-u（2*i-1））*c+（u（2*j）-u（2*i））*s）;
end

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

     文件       2871  2017-01-07 12:09  PlaneTrussfem.m

     文件       3780  2017-01-08 12:52  PlaneTrussfem_1.m

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

                 6651                    2