Beam2D2Node曾攀有限元基础教程算例3.3.7（2）

 Beam2D2Node_main

%曾攀算例3.3.7（2）
clc
clear all
%1进行单元及节点编号，建立整体坐标系，并将均布荷载简化为等效节点荷载
%2输入单元信息
E=3E11;I=6.5E-7;A=6.8E-4;
L1=1.44;L2=0.96;alpha1=0;alpha2=90;

%3计算单元刚度矩阵并组装为整体刚度矩阵
format short
Beam2D2Node=Beam2D2NodeFunc;	%调用函数接口
k1=Beam2D2Node.Stiffness(E,I,A,L1,alpha1);
k2=Beam2D2Node.Stiffness(E,I,A,L2,alpha2);
%单元2和3一样，只是节点编号不一样
KK=zeros(12,12);%4个节点总自由度为4*3
KK=Beam2D2Node.Assemble(KK,k1,1,2);
KK=Beam2D2Node.Assemble(KK,k2,3,1)
KK=Beam2D2Node.Assemble(KK,k2,4,2)

%4处理位移边界条件及节点力，求解位移
k=KK([1:6],[1:6]);
Fp=[3000;-3000;-720;0;-3000;720;0;0;0;0;0;0];%前6个已知的外荷载
u=k\Fp(1:6);%前6个是已知的外荷载
U=[u;0;0;0;0;0;0]
%5计算支反力，Fr=KK×U-Fp
P=KK*U   %下一行等价于R=P-Fp
Fr=Beam2D2Node.Forces(KK,U)-Fp	

%求解单元应变

%求解单元应力

2.Beam2D2NodeFunc

function Beam2D2Node=Beam2D2NodeFunc
Beam2D2Node.Stiffness=@Beam2D2Node_Stiffness;
Beam2D2Node.Assemble=@Beam2D2Node_Assemble;
Beam2D2Node.Forces=@Beam2D2Node_Forces;
end
function k=Beam2D2Node_Stiffness(E,I,A,L,alpha)
%以上函数计算单元刚度矩阵
%输入弹性模量E、横截面积A、惯性矩I和长度L
%输出单元刚度矩阵k(6*6)
k=[
    E*A/L,0,0,-E*A/L,0,0;
    0,12*E*I/(L^3),6*E*I/(L^2),0,-1