有限元（曾攀）B站杆单元构建及Matlab编程的代码

有限元（曾攀）杆单元构建及Matlab编程的代码

function k=Bar2D2Node_Stiffness(E,A,x1,y1,x2,y2,alpha)
%该函数计算单元的刚度矩阵
%输入弹性模量E，横截面积A
%输入第一个节点坐标（x1y1），第二个节点坐标（x2y2）角度 alpha（单位是度）
%输出单元刚度矩阵k（4X4）
%
L=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));
x=alpha*pi/180;
C=cos(x);
S=sin(x);
k=E*A/L*[C*C C*S -C*C -C*S;C*S S*S -C*S -S*S;
    -C*C -C*S C*C C*S;-C*S -S*S C*S S*S];
end

function z=Bar2D2Node_Assembbly(KK,k,i,j)
%该函数进行单元刚度矩阵的组装
%输入单元刚度矩阵k，单元的节点编号i、j
%输出整体刚度矩阵KK
%--
DOF(1)=2*i-1;
DOF(2)=2*i;
DOF(3)=2*j-1;
DOF(4)=2*j;
for n1=1:4
    for n2=1:4
        KK(DOF(n1),DOF(n2))=KK(DOF(n1),DOF(n2))+k(n1,n2);
    end
end
z=KK;

function stress=Bar2D2Node_Stress(E,x1,y1,x2,y2,alpha,u)
%该函数计算单元的应力
%输入弹性模量E，第一个节点坐标（x1y1），第二个节点坐标
%输入角度 alpha（单位是度）和单位节点位移矢量u
%