matlab program in 有限元分析基础教程_曾攀.doc

有限元分析基础教程 曾攀 3.3.6 梁单元分析的 MATLAB 程序 【MATLAB 程序】3.3.6(1) 1D 梁单元的有限元分析程序(Beam1D2Node) 编写如图 3-12 所示局部坐标系中的梁单元的刚度矩阵、单元组装、单元应力的计算程序。 解答：编写的 2 节点梁单元的五个 MATLAB 函数如下。 Beam1D2Node_Stiffness(E,I,L) 该函数计算单元的刚度矩阵，输入弹性模量 E，横截面的惯性矩 I，梁单元的长度 L，输出单元刚度矩阵 k(4×4)。 Beam1D2Node _Assembly(KK,k,i,j) 该函数进行单元刚度矩阵的组装，输入单元刚度矩阵 k，单元的节点编号 i、j、m，输出整体刚度矩阵 KK。 Beam1D2Node_ Strain(x,L,y) 该函数计算单元的几何矩阵，输入所测点距梁单元左节点的水平距离 x，输入所测点以中性层为起点的 y 方向的坐标，梁单元的长度 L，输出单元几何形状函数矩阵 B(1×4）。 Beam1D2Node _Stress(E,B,u) 该函数计算单元内某点的应力，输入弹性模量 E，几何矩阵 B，节点位移列阵 u，输出单元的应力 stress Beam1D2Node_Deflection(x,L,u) 该函数计算单元内某点的挠度，输入所测点距梁单元左节点的水平距离 x，梁单元的长度 L，节点位移列 阵 u，输出该点的挠度 v。 基于第 3.3.2 节中的基本公式，可以编写出具体实现以上每个函数的 MATLAB 程序如下。 %%%%%%%%%% Beam1D2Node %%% begin %%%%%%%%%%%% function k =Beam1D2Node_Stiffness(E,I,L) %该函数计算单元的刚度矩阵 %输入弹性模量 E，横截面的惯性矩 I，梁单元的长度 L %输出单元刚度矩阵 k(4×4) %----------------------------------------- k = E*I/(L*L*L)*[12 6*L -12 6*L ; 6*L 4*L*L -6*L 2*L*L ; -12 -6*L 12 -6*L ; 6*L 2*L*L -6*L 4*L*L]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%% function z = Beam1D2Node_Assembly(KK,k,i,j) %该函数进行单元刚度矩阵的组装 %输入单元刚度矩阵 k，单元的节点编号 i、j、m %输出整体刚度矩阵 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 B= Beam1D2Node_Strain(x,L,y) %该函数计算单元的几何矩阵 %输入所测点距梁单元左节点的水平距离 x %输入所测点以中性层为起点的 y 方向的坐标，梁单元的长度 L %输出单元几何形状函数矩阵 B(1×4） %----------------------------------------- e=x/L; B1=(12*e-6)/(L*L); B2=(6*e-4)/L; B3=-(12*e-6)/(L*L); B4=(6*e-2)/L; B=-y*[B1,B2,B3,B4]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%% function stress= Beam1D2Node_Stress(E,B,u) %该函数计算单元内某点的应力 %输入弹性模量 E，几何矩阵 B，节点位移列阵 u %输出单元的应力 stress %----------------------------------------- stress = E*B*u; %%%%%%%%%%%%%%%%%%%%%%%%%%%%% function v=Beam1D2Node _Deflection(x,L,u) %该函数计算单元内某点的挠度 %输入所测点距梁单元左节点的水平距离 x %输入梁单元的长度 L，节点位移列阵 u %输出该点的挠度 v %----------------------------------------- e=x/L; N1=1-3*e*e+2*e*e*e; N2=L(e-2*e*e+e*e*e);