triangleelementstructureparameters.m

在有限元分析中计算结构整体刚度矩阵、单元刚度矩阵、单元结点位移、应力应变等

%三角单元结构参数计算
%%%%%%%%%%%%%%%%输入参数
E=5*10^5;%弹性模量MPa
nu=0.3;%泊松比
t=1;%厚度cm
n=5;%结点个数
esn=4;%单元个数
ECN=[
    1   1   1   1;
    2   3   4   5;
    3   5   5   2
    ];%单元结点编号矩阵
CC=[
    0   5   -5   -5  5;
    0   5   5   -5  -5
    ];%结点坐标矩阵

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%控制程序
%%%%%%%%%结构整体刚度矩阵
K=zeros(2*n,2*n);
for e=1:1:esn
    ecn=ECN(:,e);%单元结点编号向量
    Ge=fGe_1(n,ecn);
    ECC=[CC(:,ecn(1)) CC(:,ecn(2)) CC(:,ecn(3))];%单元结点坐标阵
    Ke=fKe_1(E,nu,t,ECC);%单元刚度阵
    K=K+Ge'*Ke*Ge;
end
detK=det(K);

%%%%%%%%%结构结点载荷列阵
P=[-1000;0;1000;0;-1000;0;1000;0;0;0];%2n×1的载荷向量，单位N

%%%%%%%%%结构结点位移列阵
a=inv(K)*P;

%%%%%%%%%单元参数计算
e=2;%输入单元编号
ecn=ECN(:,e);%单元结点编号向量
ECC=[CC(:,ecn(1)) CC(:,ecn(2)) CC(:,ecn(3))];%单元结点坐标阵
Ae=fAe_1(ECC);%单元面积
Be=fBe_1(ECC);%单元应变矩阵
Ge=fGe_1(n,ecn);%单元结点自由度转换矩阵
Ke=fKe_1(E,nu,t,ECC);%单元刚度矩阵
ae=Ge*a;%单元结点位移列阵
epsilon=Be*ae;%单元应变矩阵
D=E/(1-nu^2)*[1 nu 0;nu 1 0;0 0 (1-nu)/2];
sigma=D*Be*ae;%单元应力矩阵