python结课设计力学方面_计算力学课程设计.docx

计算力学课程设计

计算力学课程设计说明书班级 姓名 学号 日期 2015年12月25日桁架的有限元计算方法引言 有限元分析是求取复杂微分方程近似解的一种非常有效的工具，是现代数字化科技的一种重要的基础性原理。将它应用于工程技术中，可成为工程设计和分析的可靠工具。而在桁架结构中，运用有限元的方法，通过现代有限元分析软件如MATLAB，ANSYS等可轻易的求得各个杆件的内力等。例如下图的桁架结构，运用有限元法可十分清晰的了解各杆件的受力状况。基本力学模型 如下图1所示 图1有限元计算原理 首先，明确Matlab程序要实现的5个重要模块分别为：单元刚度矩阵的求解、单元组装、节点位移的求解、单元应力的求解、节点力的求解。下面给出这5个模块的实现。(1)单元刚度矩阵求解 定义函数Bar2D2Node_Stiffness，该函数计算单元的刚度矩阵，输入弹性模量E，横截面积A，两个节点坐标输出单元刚度矩阵k(4X4)。具体代码如下：function k=Bar2D2Node_Stiffness(E,A,x1,y1,x2,y2)L=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));x=acos((x2-x1)/L);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];(2)单元组装 定义函数Bar2D2Node_Assembly，该函数进行单元刚度矩阵的组装，输入单元刚度矩阵k，单元的节点编号i、j。输出整体刚度矩阵KK，具体代码如下：function z=Bar2D2Node_Assembly(KK,k,i,j)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); endendz=KK;(3) 节点位移的求解 定义函数Bar2D2Node_Disp(KK,num,p)，该函数输入KK为总体刚度矩阵；num为活动自由度编号数组；p为活动自由度方向上的节点力；输出节点位移列阵。具体代码如下：function u=Bar2D2Node_Disp(KK,num,p)k=KK(num,num);u=k\p;(4)单元应力的求解 定义函数函数Bar2D2Node_Stress(E,x1,y1,x2,y2,u)，该函数计算单元的应力输入弹性模量E，第一个节点坐标(x1,y1)，第二个节点坐标(x2,y2)单元节点位移矢量u，返回单元应力标量stress 。具体代码如下：function stress=Bar2D2Node_Stress(E,x1,y1,x2,y2,u)L=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));x=acos((x2-x1)/L);C=cos(x);S=sin(x);stress=E/L*[-C -S C S]*u;(5)计算节点力 定义函数Bar2D2Node_Forces(KK,q)，该函数用于计算节点力，KK为刚度矩阵，q为节点位移阵列。function P= Bar2D2Node_Forces(KK,q)P=KK*q; 以上5个函数写入MATLAB的M文件中完成函数的定义，以备命令窗口的调用。 至此，基于MATLAB的杆单元有限元分析的程序设计的准备工作已经完成。计算结果及分析 针对上面的具体模型运用MATLAB实现有限元的计算：(1) 结构的离散化与编号 对该结构进行自然离散，节点编号和单元编号如上图所示(2)计算各单元的刚度矩阵(基于国际标准单位) 输入弹性模量E、横截面积A，各点坐标。然后分别针对单元1，2，3和4，调用4次Bar2D2Node_Stiffness，就可以得到单元的刚度矩阵。对应的主程序中代码：E=2.95e11;A=0.0001;x1=0;y1=0;x2=0.4;y2=0;x3=0.4;y3=0.3;x4=0;y4=0.3;k1=Bar2D2Node_Stiffness (E,A,x1,y1,x2,y2)k2=Bar2D2Node_Stiffness (E,A,x2,y2,x3,y3)k3=Bar2D2Node_Stiffness (E,A,x1,y1,x3,y3)k4=Bar2D2Node_Stiffness (E,A,x4,y4,x3,y3)计算结果如下 图2单元1与单元2的刚度矩阵 图3 单元3与单元4的刚度矩阵(3) 建立整体刚度方程由于该结构共有4个节点，因此，设置结构总的刚度矩阵为KK(8×8)，先