有限元平面应力及桁架机构matlab代码

包含一个桁架机构算例和平面应力算例，均较为基础，代码后大多有注释，可供初学者参考

%% 薄板算例（直接法——平面应力问题）
close all
clear all
clc

%% 基本参数定义
numele=4;
numnode=5;
elenode=3;
L=1;
t=1;%单元厚度
elenodecorre=[1 2 3;1 3 4;2 5 3;3 5 4];
E=206;%合金钢，单位Gpa
mu=0.25;%合金钢，0.24-0.28
q=20; %均布拉应力
%% 施加载荷前作图
%inicoornode=[0 0;1 0;1/2 1/2;0 1;1 1];
ininodevector=[0 0 1 0 1/2 1/2 0 1 1 1];
figure
axis off
axis equal
hold on
for inicoori=1:numnode
    plot（ininodevector（2*inicoori-1）ininodevector（2*inicoori）‘bo‘）;
end
for inicoorj=1:numele
    hold on
    line（[ininodevector（2*elenodecorre（inicoorj1）-1） ...
        ininodevector（2*elenodecorre（inicoorj2）-1）]...
        [ininodevector（elenodecorre（inicoorj1）*2） ...
        ininodevector（elenodecorre（inicoorj2）*2）]）;
    hold on
    line（[ininodevector（2*elenodecorre（inicoorj2）-1） ...
        ininodevector（2*elenodecorre（inicoorj3）-1）]...
        [ininodevector（elenodecorre（inicoorj2）*2） ...
        ininodevector（elenodecorre（inicoorj3）*2）]）;
    hold on
    line（[ininodevector（2*elenodecorre（inicoorj3）-1） ...
        ininodevector（2*elenodecorre（inicoorj1）-1）]...
        [ininodevector（elenodecorre（inicoorj3）*2） ...
        ininodevector（elenodecorre（inicoorj1）*2）]）;
end
%hold on

%%  
sysmatK=zeros（1010）;
u=zeros（81）;
delt=zeros（101）;
f=zeros（101）;


for elei=1:numele
    % initical coordinates of the nodes of each element
    x1=ininodevector（elenodecorre（elei1）*2-1）;
    y1=ininodevector（elenodecorre（elei1）*2）;
    x2=ininodevector（elenodecorre（elei2）*2-1）;
    y2=ininodevector（elenodecorre（elei2）*2）;
    x3=ininodevector（elenodecorre（elei3）*2-1）;
    y3=ininodevector（elenodecorre（elei3）*2）;    
    % interpolation function （shape function）
    area=0.5*（x2*y3+x1*y2+x3*y1-x1*y3-x2*y1-x3*y2）;
    
    % 对于每个节点进行计算
    %N1=0.5*（x2*y3-x3*y2+（y2-y3）*x+（x3-x2）*y）/area;
    %N2=0.5*（x3*y1-x1*y3+（y3-y1）*x+（x1-x3）*y）/area;
    %N3=0.5*（x1*y2-x2*y1+（y1-y2）*x+（x2-x1）*y）/area;
    b1=y2-y3;
    b2=y3-y1;
    b3=y1-y2;
    c1=x3-x2;
    c2=x1-x3;
    c3=x2-x1;
    B=[b1 0 b2 0 b3 0;0 c1 0 c2 0 c3;...
        ;c1 b1 c2 b2 c3 b3]/（2*area）;%由形函数偏导而来的单元应变矩阵
    D=[1 mu 0;mu 1 0;0 0 （1-mu）/2]*E/（1-mu^2）;%状态弹性矩阵
    K=t*B‘*D*B; %单元刚度矩阵
    A=area*B‘*t;%节点应力与单元边上应力之间的关系矩阵
    S=D*B;%单元应力矩阵（单元位移向量与应力向量的转化关系）
    inditemp=elenodecorre（elei:）;
    indicator=[2*inditemp（1）-1  2*inditemp（1） 2*inditemp（2）-1 ... 
        2*inditemp（2） 2*inditemp（3）-1  2*inditemp（3）];
    sysmatK（indicatorindicator）=K+sysmatK（indicatorindicator）; %系统刚度矩阵组装

end

f（3）=q*L*t;
f（9）=q*L*t;
soluK=sysmatK（[3 4 5 6 9 10][3 4 5 6 9 10]）;
soluf=f（[3 4 5 6 9 10]）;
soludelt=inv（soluK）*soluf;
delt（[3 4 5 6 9 10]）=soludelt;
delt
%% 施加载荷后作图
finnodevector=ininodevector‘+delt

hold on
for fincoori=1:numnode
    plot（finnodevector（2*fincoori-1）finnodevector（2*fincoori）‘ro‘）;
end
for fincoorj=1:numele
    hold on
    line（[finnodevector（2*elenodecorre（fincoorj1）-1） ...
        finnodevector（2*elenodecorre（fincoorj2）-1）]...

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     文件        2862  2019-01-05 09:48  有限元平面应力及桁架机构matlab代码\Truss.m
     文件        3830  2019-01-03 16:46  有限元平面应力及桁架机构matlab代码\thin_plate.m
     文件          79  2019-01-11 18:54  有限元平面应力及桁架机构matlab代码\文件说明.txt
     目录           0  2019-01-11 18:52  有限元平面应力及桁架机构matlab代码\