2D四杆桁架结构的有限元分析实例.docx

资源描述

2D四杆桁架结构的有限元分析实例.docx

《2D四杆桁架结构的有限元分析实例.docx》由会员分享，可在线阅读，更多相关《2D四杆桁架结构的有限元分析实例.docx（15页珍藏版）》请在冰豆网上搜索。

2D四杆桁架结构的有限元分析实例.docx

2D四杆桁架结构的有限元分析实例

实例：

2D四杆桁架结构的有限元分析

学习有限元方法的一个最佳途径，就是在充分掌握基本概念的基础上亲自编写有限元分析程序，这就需要一个良好的编程环境或平台。

MATLAB软件就是这样一个平台，它以功能强大、编程逻辑直观、使用方便见长。

将提供有限元分析中主要单元完整的MATLAB程序，并给出详细的说明。

1D杆单元的有限元分析MATLAB程序（Bar1D2Node）

最简单的线性杆单元的程序应该包括单元刚度矩阵、单元组装、单元应力等几个基本计算程序。

下面给出编写的线性杆单元的四个MATLAB函数。

Bar1D2Node_Stiffness（E,A，L）

该函数计算单元的刚度矩阵，输入弹性模量E,横截面积A和长度L，输出单元刚度矩阵k（2×2）。

Bar1D2Node_Assembly（KK，k，i，j）

该函数进行单元刚度矩阵的组装，输入单元刚度矩阵k,单元的节点编号i、j，输出整体刚度矩阵KK。

Bar1D2Node_Stress（k，u，A）

该函数计算单元的应力，输入单元刚度矩阵k、单元的位移列阵u（2×1）以及横截面积A计算单元应力矢量,输出单元应力stress.

Bar1D2Node_Force（k,u）

该函数计算单元节点力矢量，输入单元刚度矩阵k和单元的位移列阵u（2×1），输出2×1的单元节点力矢量forces。

基于1D杆单元的有限元分析的基本公式，写出具体实现以上每个函数的MATLAB程序如下。

%％%％%%%％％%%Bar1D2Node%％begin％％％%%%%%%

functionk=Bar1D2Node_Stiffness（E，A，L）

%该函数计算单元的刚度矩阵

％输入弹性模量E，横截面积A和长度L

％输出单元刚度矩阵k（2×2）

%—-—-——---—-——---——————--—-—-——--——--——-

k=［E＊A/L-E＊A/L;-E＊A/LE＊A/L]；

％%%％%%％％%％%%%％％%％%%％％%％％%％

functionz=Bar1D2Node_Assembly（KK,k，i，j）

%该函数进行单元刚度矩阵的组装

％输入单元刚度矩阵k，单元的节点编号i、j

％输出整体刚度矩阵KK

%—-————----————--—-----—————--——---—

DOF

（1）=i;

DOF

（2）=j;

forn1=1:

forn2=1：

KK（DOF（n1）,DOF（n2））=KK（DOF（n1），DOF（n2））+k（n1，n2）；

end

z=KK;

％-—-----—----—----——————-——--—--——-----———---—--—--—----—----

functionstress=Bar1D2Node_Stress（k，u，A）

%该函数计算单元的应力

％输入单元刚度矩阵k，单元的位移列阵u（2×1）

％输入横截面积A计算单元应力矢量

%输出单元应力stress

%—-—-—--———--————-----——--——-———----

stress=k*u/A;

%-———-—---—————--—--——---———-—-———-—-—----——-——----—-—-----—

％％%%%％%%％％％％％％％%％％%%%％%%％

functionforces=Bar1D2Node_Force（k，u）

%该函数计算单元节点力矢量

%输入单元刚度矩阵k和单元的位移列阵u（2×1）

％输出2×1的单元节点力分量forces

%————-—-—-—---——-—--—---——————-—-——————-——

forces=k*u;

%%%％%％%%%％%Bar1D2Node%％end％%％％％%%%%

【四杆桁架结构的有限元分析—数学推导】

如图所示的结构，各杆的弹性模量和横截面积都为E=29。

54×10N/mm2，A=100mm2,试求解该结构的节点位移、单元应力以及支反力。

图1四杆桁架结构

解答：

对该问题进行有限元分析的过程如下.

（1）结构的离散化与编号

对该结构进行自然离散，节点编号和单元编号如图1所示，有关节点和单元的信息见表1-表3。

表1节点及坐标表2单元编号及对应节点表3各单元的长度及轴线方向余弦

节点

单元

节点1

节点2

单元

①

400

②

300

—1

400

300

③

500

0。

300

④

400

（2）各个单元的矩阵描述

由于所分析的结构包括有斜杆，所以必须在总体坐标下对节点位移进行表达，所推导的单元刚度矩阵也要进行变换，各单元经坐标变换后的刚度矩阵如下。

（3）建立整体刚度方程

将所得到的各个单元刚度矩阵按节点编号进行组装，可以形成整体刚度矩阵，同时将所有节点载荷也进行组装.

刚度矩阵：

K=K

（1）+K

（2）+K（3）+K（4）

节点位移：

q=［u1v1u2v2u3v3u4v4]T

节点力:

P=R+F=［Rx1Ry12×104Ry202。

5×104Rx4Ry4］T

其中（Rx1，Ry1）为节点1处沿x和y方向的支反力，Ry2为节点2处y方向的支反力，（Rx4，Ry4）为节点4处沿x和y方向的支反力。

整体刚度方程为

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

边界条件BC（u）为：

u1=v1=v2=u4=v4=0，代入整体刚度方程中，经化简后有

对该方程进行求解，有

则所有的节点位移为

（5）各单元应力的计算

其中T为坐标转换矩阵；同理，可求出其它单元的应力.

（6）支反力的计算

将节点位移的结果代入整体刚度方程中，可求出

2D杆单元的有限元分析程序（Bar2D2Node）

编写平面桁架单元的单元刚度矩阵、单元组装、单元应力的计算程序.编写的平面桁架单元的四个MATLAB函数如下。

Bar2D2Node_Stiffness（E，A，x1，y1,x2,y2,alpha）

该函数计算单元的刚度矩阵，输入弹性模量E，横截面积A，第一个节点坐标（x1,y1），第二个节点坐标（x2,y2）和角度alpha（单位是度），输出单元刚度矩阵k（4×4）。

Bar2D2Node_Assembly（KK,k，i,j）

该函数进行单元刚度矩阵的组装，输入单元刚度矩阵k，单元的节点编号i、j,输出整体刚度矩阵KK。

Bar2D2Node_Stress（E，x1，y1，x2，y2,alpha，u）

该函数计算单元的应力，输入弹性模量E，第一个节点坐标（x1,y1），第二个节点坐标（x2,y2），角度alpha（单位是度）和单位节点位移矢量u，返回单元应力标量。

Bar2D2Node_Forces（E,A，x1，y1，x2，y2,alpha，u）

该函数计算单元的应力，输入弹性模量E，横截面积A,第一个节点坐标（x1,y1），第二个节点坐标（x2,y2），角度alpha（单位是度）和单元节点位移矢量u,返回单元节点力。

基于2D杆单元的基本公式，可以编写出具体实现以上每个函数的MATLAB程序如下。

％%%%％%%％%％%Bar2D2Node％%begin％％％%%％％％%％％％％％

functionk=Bar2D2Node_Stiffness（E,A,x1，y1,x2，y2,alpha）

%该函数计算单元的刚度矩阵

%输入弹性模量E，横截面积A

%输入第一个节点坐标（x1,y1），第二个节点坐标（x2,y2），角度alpha（单位是度）

％输出单元刚度矩阵k（4×4）

％—--——-—-—---—-—-—-——-—--—-——-—----——-—------———--

L=sqrt（（x2—x1）＊（x2-x1）+（y2-y1）＊（y2—y1））;

x=alpha＊pi/180;

C=cos（x）;

S=sin（x）；

k=E*A/L*[C*CC*S—C＊C—C*S;

C＊SS*S—C*S—S＊S;

—C*C-C＊SC*CC*S；

-C＊S—S＊SC*SS＊S];

％%%％%％％%%%％%％％％％％%％%％％％％%%%％%%％

functionz=Bar2D2Node_Assembly（KK，k，i,j）

％该函数进行单元刚度矩阵的组装

％输入单元刚度矩阵k,单元的节点编号i、j

％输出整体刚度矩阵KK

％—--—-——-————-—-—--———--———————-————---——-—————--—--—————

DOF

（1）=2＊i-1；

DOF

（2）=2*i;

DOF（3）=2*j—1;

DOF（4）=2＊j；

forn1=1:

forn2=1:

KK（DOF（n1）,DOF（n2））=KK（DOF（n1）,DOF（n2））+k（n1，n2）;

end

z=KK；

%％％%％%％％%％%％％%％％％％%％%%%％％%％%％％％%%%—-

functionstress=Bar2D2Node_Stress（E，x1，y1,x2，y2,alpha,u）

％该函数计算单元的应力

%输入弹性模量E，第一个节点坐标（x1,y1）,第二个节点坐标（x2,y2）

%输入角度alpha（单位是度）和单位节点位移矢量u

％返回单元应力标量stress

％——--——-—-——————-—-----—--------———--——-—-—-—————

L=sqrt（（x2-x1）*（x2-x1）+（y2—y1）＊（y2-y1））；

x=alpha*pi/180;

C=cos（x）；

S=sin（x）；

stress=E/L*[—C—SCS]*u；

%％%%％%％％%％%％％％％%％%%%％%%%%％%％％％％％%

functionforces=Bar2D2Node_Forces（E，A,x1，y1,x2，y2，alpha，u）

%该函数计算单元的应力

％输入弹性模量E，横截面积A

％输入第一个节点坐标（x1,y1），第二个节点坐标（x2，y2），角度alpha（单位是度）

%输入单元节点位移矢量u

％返回单元节点力forces

%-———-—---—--——————--——--——-——---———-—--———-—-—-—-—-—-—--———--

L=sqrt（（x2-x1）*（x2—x1）+（y2-y1）＊（y2—y1））；

x=alpha＊pi/180;

C=cos（x）；

S=sin（x）；

forces=E*A/L＊［-C-SCS]＊u;

%％%%％%％%%%％Bar2D2Node％%end％％%%％%%％%%%％%％

【四杆桁架结构的有限元分析-MATLAB-（Bar2D2Node）】

仍就图1所示结构，基于MATLAB平台求解该结构的节点位移、单元应力以及支反力。

解答：

对该问题进行有限元分析的过程如下。

（1）结构的离散化与编号

对该结构进行自然离散,节点编号和单元编号如图所示,有关节点和单元的信息见表1—表3.

（2）计算各单元的刚度矩阵（基于国际标准单位）

建立一个工作目录，将所编制的用于平面桁架单元分析的四个MATLAB函数放置于该工作目录中,分别以各自函数的名称给出文件名,即：

Bar2D2Node_Stiffness，

Bar2D2Node_Assembly,

Bar2D2Node_Stress，

Bar2D2Node_Forces。

然后启动MATLAB，将工作目录设置到已建立的目录中，在MATLAB环境中,输入弹性模量E、横截面积A，各点坐标x1,y1，x2，y2，x3,y3，x4，y4，角度alpha1，alpha2和alpha3,然后分别针对单元1，2,3和4，调用四次Bar2D2Node_Stiffness,就可以得到单元的刚度矩阵.相关的计算流程如下。

E=2。

95e11；

A=0。

0001；

x1=0;

y1=0;

x2=0.4;

y2=0；

x3=0。

y3=0.3；

x4=0；

y4=0.3;

alpha1=0；

alpha2=90；

alpha3=atan（0。

75）＊180/pi；

k1=Bar2D2Node_Stiffness（E，A,x1,y1,x2，y2，alpha1）

k2=Bar2D2Node_Stiffness（E，A,x2,y2,x3,y3，alpha2）

k3=Bar2D2Node_Stiffness（E,A，x1,y1,x3，y3,alpha3）

k4=Bar2D2Node_Stiffness（E，A，x4,y4，x3，y3，alpha1）

（3）建立整体刚度方程

由于该结构共有4个节点，因此,设置结构总的刚度矩阵为KK（8×8），先对KK清零，然后四次调用函数Bar2D2Node_Assembly进行刚度矩阵的组装。