有限元基础教程（曾攀）习题3-13Matlab实现

---

习题3-13

有一个梁结构如图所示，试采用梁单元来分析该结构，并求取支反力。

解题步骤

1. 进行单元及节点编号，建立整体坐标系，如上图所示；并将均布荷载简化为等效节点荷载，如下图所示。
   
2. 输入单元信息

>>E=3E10;I=204E-8;A=7.65E-4;
>>L1=1;L2=0.9;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);
>>KK=zeros(9,9);
>>KK=Beam2D2Node.Assemble(KK,k1,1,2);
>>KK=Beam2D2Node.Assemble(KK,k2,3,2)

输出如下结果：

4. 处理位移边界条件及节点力，求解位移

>>k=KK([4:6],[4:6]);
>>Fp=[-62500;-52083;-93750;39062;-31250;13021];
>>u=k\Fp(4:6);
>>U=[0;0;0;u]

输出如下：
U=[0 0 0 -0.0014 -0.0086 -0.0041] (将列排成了行)

5. 计算支反力，Fr=KK×U-Fp

Fr=Beam2D2Node.Forces(KK,U)-Fp	

输出：
Fr =1.0e+05 *[0.5469 0.3906 1.3281 -0.0000 -0.0000 0] (将列排成了行)

函数文件

（所有函数都保存在同一个m文件中）
函数文件名：Beam2D2NodeFunc.m

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,-12*E*I/(L^3),6*E*I/(L^2);
    0,6*E*I/(L^2),4*E*I/L,0,-6*E*I/(L^2),2*E*I/L;
    -E*A/L,0,0,E*A/L,0,0;
    0,-12*E*I/(L^3),-6*E*I/(L^2),0,12*E*I/(L^3),-6*E*I/(L^2);
    0,6*E*I/(L^2),2*E*I/L,0,-6*E*I/(L^2),4*E*I/L
    ];
x=alpha*pi/180;
C=cos(x);
S=sin(x);
T=[
    C,S,0,0,0,0;
    -S,C,0,0,0,0;
    0,0,1,0,0,0;
    0,0,0,C,S,0;
    0,0,0,-S,C,0;
    0,0,0,0,0,1
    ];
k=T'*k*T;
end
function y=Beam2D2Node_Assemble(KK,k,i,j)
%以上函数进行单元刚度矩阵的组装
%输入单元刚度矩阵k和单元的节点编号i,j
%输出整体刚度矩阵KK
DOF(1)=3*i-2;
DOF(2)=3*i-1;
DOF(3)=3*i;
DOF(4)=3*j-2;
DOF(5)=3*j-1;
DOF(6)=3*j;
for n1=1:6
    for n2=1:6
        KK(DOF(n1),DOF(n2))=KK(DOF(n1),DOF(n2))+k(n1,n2);
    end
end
y=KK;
end
function forces=Beam2D2Node_Forces(k,u)
%以上函数计算单元的应力，输入单元刚度矩阵k、节点位移u
%输出单元节点力forces
forces=k*u;
end

参考文献

有限元基础教程（曾攀）