计算力学课程设计

计算力学课程设计说明书

班级 姓名 学号

日期2015年12月25日

计算力学课程设计

桁架的有限元计算方法

引言

有限元分析是求取复杂微分方程近似解的一种非常有效的工具，是现代数字化科技的一种重要的基础性原理。将它应用于工程技术中，可成为工程设计和分析的可靠工具。而在桁架结构中，运用有限元的方法，通过现代有限元分析软件如MATLAB，ANSYS等可轻易的求得各个杆件的内力等。例如下图的桁架结构，运用有限元法可十分清晰的了解各杆件的受力状况。

1

基本力学模型

如下图1所示

图1

2 有限元计算原理

首先，明确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);

1安徽理工大学理学院

计算力学课程设计

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); end end z=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)

2安徽理工大学理学院

计算力学课程设计

P=KK*q;

以上5个函数写入MATLAB的M文件中完成函数的定义，以备命令窗口的调用。 至此，基于MATLAB的杆单元有限元分析的程序设计的准备工作已经完成。

3 计算结果及分析

针对上面的具体模型运用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 单元3与单元4的刚度矩阵

（3） 建立整体刚度方程

由于该结构共有4个节点，因此，设置结构总的刚度矩阵为KK(8×8)，先对KK清零，然后四次调用函数Bar2D2Node _Assembly进行刚度矩阵的组装。相关主程序代码为：

KK=zeros(8,8);%由于该结构共有4个节点，因此，结构总的刚度矩阵为KK(8×8)，先对K清零，然后四次调用函数Bar2D2Node _Assembly进行刚度矩阵的组装。 KK=Bar2D2Node_Assembly (KK,k1,1,2); KK=Bar2D2Node_Assembly (KK,k2,2,3); KK=Bar2D2Node_Assembly (KK,k3,1,3); KK=Bar2D2Node_Assembly (KK,k4,4,3);

结果如下

4安徽理工大学理学院

计算力学课程设计

图4 整体刚度矩阵

（4）边界条件的处理及刚度方程的求解

由图可以看出，被约束的自由度有：节点1的x,y方向自由度，节点2的y方向自由度，4节点的x、y方向两个自由度。则活动自由度编号为3，5，6.活动自由度对应的节点载荷F3=20000N，F5=0N，F6=25000N，采用高斯消去法进行求解，对应的代码为： num=[3,5,6];%可活动的自由度编号 p=[20000;0;-25000]; u=Bar2D2Node_Disp(KK,num,p)

结果如下

由此可知u2=0.2712*10^-3,u3=0.0565*10^-3,v3=-0.2225*10^-3。

（5）支反力的计算

在得到整个结构的节点位移后，由原整体刚度方程就可以计算出对应的支反力。这部分对应的主程序的代码如下： q=zeros(8,1);

q(num)=u;%节点位移阵列 P=Bar2D2Node_Forces(KK,q)

结果如下

5安徽理工大学理学院

计算力学课程设计

由此可知支反力为FRX1=-15833N,FRY1=3125N,FRY2=21879N,FRX4=-4167N,FRY4=0N。

（6）单元应力的计算

先从整体位移列阵q中提取出单元的位移列阵，然后，调用计算单元应力的函数Bar2D2Node_Stress，就可以得到各个单元的应力分量。 u1=[q(1);q(2);q(3);q(4)]

stress1=Bar2D2Node_Stress(E,x1,y1,x2,y2,u1) u2=[q(3);q(4);q(5);q(6)]

stress2=Bar2D2Node_Stress(E,x2,y2,x3,y3,u2) u3=[q(1);q(2);q(5);q(6)]

stress3=Bar2D2Node_Stress(E,x1,y1,x3,y3,u3) u4=[q(7);q(8);q(5);q(6)]

stress4=Bar2D2Node_Stress(E,x4,y4,x3,y3,u4)

结果如下：

stress1 =200000000，stress2 =-2.1875e+008，stress3 =-5.2083e+007，stress4 =4.1667e+007

参考文献

[1] 曾攀.有限元基础教程.北京：高等教育出版社2009,7 [2] 王勖成.有限单元法.北京:清华大学出版社,2003,7.

附录（程序）

子程序1：

function k=Bar2D2Node_Stiffness(E,A,x1,y1,x2,y2) L=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));

6安徽理工大学理学院

计算力学课程设计

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：

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); end end z=KK; 子程序3：

function u=Bar2D2Node_Disp(KK,num,p) k=KK(num,num); u=k\\p; 子程序4：

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：

function P= Bar2D2Node_Forces(KK,q) P=KK*q; 主程序： E=2.95e11; A=0.0001; x1=0; y1=0; x2=0.4; y2=0; x3=0.4; y3=0.3;

7安徽理工大学理学院

计算力学课程设计

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) KK=zeros(8,8);

KK=Bar2D2Node_Assembly (KK,k1,1,2); KK=Bar2D2Node_Assembly (KK,k2,2,3); KK=Bar2D2Node_Assembly (KK,k3,1,3); KK=Bar2D2Node_Assembly (KK,k4,4,3) num=[3,5,6]; p=[20000;0;-25000]; u=Bar2D2Node_Disp(KK,num,p) q=zeros(8,1); q(num)=u;

P=Bar2D2Node_Forces(KK,q) u1=[q(1);q(2);q(3);q(4)];

stress1=Bar2D2Node_Stress(E,x1,y1,x2,y2,u1) u2=[q(3);q(4);q(5);q(6)];

stress2=Bar2D2Node_Stress(E,x2,y2,x3,y3,u2) u3=[q(1);q(2);q(5);q(6)];

stress3=Bar2D2Node_Stress(E,x1,y1,x3,y3,u3) u4=[q(7);q(8);q(5);q(6)];

stress4=Bar2D2Node_Stress(E,x4,y4,x3,y3,u4)

8安徽理工大学理学院

本文来源：https://www.bwwdw.com/article/qhwd.html