基于MATLAB的有限元法结构分析.docVIP

基于MATLAB的有限元法结构分析 　　摘要:本文结合平面梁单元实例，简单分析了用MATLAB进行有限元法分析结构的过程，有效证明了该方法的可行性和实用性。 　　关键词:MATLAB有限元法平面梁单元结构分析 　　中图分类号： O723 文献标识码： A 文章编号： 　　 　　1概述 　　MATLAB是一门计算机编程语言，专门以矩阵的方式来处理计算机数据。它将高性能的数值计算和可视化集成在一起，功能强大，适合多学科、多部门的要求，已广泛地应用于数值计算、图形处理、符号计算、数学建模、小波分析、系统辨识、实时控制和动态仿真等研究领域，是科学研究、工程设计和运算的得力助手。而在结构力学中常用的有限元法也可以用MATLAB软件快速准确的计算出结果。本文简单介绍了用MATLAB进行有限元分析的方法，并以较为常见的平面梁元为例近一步分析MATLAB的编程方法。 　　2用MATLAB进行有限元法分析过程概述 　　结构的离散化：将要分析的结构物分割成有限个单元体，并在单元体的指定点设置结点，使相邻单元的有关参数具有一定的连续性，并构成一个单元的集合体，以它代替原来的结构。[2] 　　写出单元刚度矩阵：用MATLAB写出域内每个单元的单元刚度矩阵。[1] 　　集成整体刚度矩阵：用直接刚度法写出整体刚度矩阵。 　　引入边界条件：引入支座、外加荷载、位移等边界条件。 　　求解方程：分解整体刚度矩阵并用高斯法求解方程组。 　　后处理：根据求得的结构节点位移矢量，求解支座反力、单元节点力和单元应力。 　　3MATLAB有限元分析编程方法 　　借助MATLAB的编程功能，可以写出三个M函数文件。 　　BeamNode_Stiffness(E,I,A,L)——计算单元的刚度矩阵，输入弹性模量E，横截面积A，惯性矩I，长度L，输出单元刚度矩阵k； 　　BeamNode_Assemble(KK,k,i,j)——进行单元刚度矩阵的组装，输入单元刚度矩阵k，单元的节点编号i、j，输出整体刚度矩阵KK； 　　BeamNode_Forces(k,u)——计算单元的应力，输入单元刚度矩阵k，节点位移u，输出单元节点力forces。 　　使用上述M函数时首先要给出参数(E,I,A,L,i,j)的值才可计算，而各M函数的具体编程代码可参考文献[2]。 　　4实例分析 　　如图所示的框架结构，其顶端受均布力作用，用有限元方法分析该结构的位移。结构中各个截面的参数都为：，，。 　　 　　 　　对该问题的解答过程如下： 　　4.1结构的离散化 　　将该结构离散为三个单元，有关节点和单元的信息见表1。 　　4.2写出单元刚度矩阵 　　首先在MATLAB环境下，输入弹性模量E、横截面积A、惯性矩I、长度L，然后针对单元1和单元2单元3，分别二次调用函数BeamNode_Stiffness，就可以得到单元的刚度矩阵k1(6×6)和k2(6×6)，且单元2和单元3的刚度矩阵相同。 　　 E=3E11; 　　 I=6.5E-7; 　　 A=6.8E-4; 　　 L1=1.44; 　　 L2=0.96; 　　 k1=BeamNode_Stiffness(E,I,A,L1); 　　 k2=BeamNode_Stiffness(E,I,A,L2); 　　4.3集成整体刚度矩阵 　　将单元2和单元3的刚度矩阵转换成整体坐标下的形式。由于该结构共有4个节点，则总共的自由度数为12，因此，结构总的刚度矩阵为KK(12×12)，对KK清零，然后两次调用函数BeamNode_Assemble进行刚度矩阵的组装。 　　 T=[0,1,0,0,0,0;-1,0,0,0,0,0;0,0,1,0,0,0;0,0,0,0,1,0;0,0,0,-1,0,0;0,0,0,0,0,1]; 　　 k3=T.*k2.*T; 　　 KK=zeros(12,12); 　　 KK=BeamNode_Assemble(KK,k1,1,2); 　　 KK=BeamNode_Assemble(KK,k3,3,1); 　　 KK=BeamNode_Assemble(KK,k3,4,2); 　　4.4引入边界条件并求解方程 　　该问题的位移边界条件为。因此，将针对节点1和节点2的位移进行求解，节点1和节点2的位移将对应KK矩阵中的前6行和前6列，则需从KK(12×12)中提出，置给k，然后生成对应的载荷列阵p，再采用高斯消去法进行求解。 　　 k=KK(1:6,1:6); 　　 p=[3000;-3000;-720;0;-3000;720]; 　　 u=k\p 　　4.5后处理 　　在得到整个结构的节点位移后，由原整体刚度方程就可以计算出对