2D四杆桁架结构的有限元分析实例.docx
- 文档编号:29876832
- 上传时间:2023-08-03
- 格式:DOCX
- 页数:14
- 大小:154.45KB
2D四杆桁架结构的有限元分析实例.docx
《2D四杆桁架结构的有限元分析实例.docx》由会员分享,可在线阅读,更多相关《2D四杆桁架结构的有限元分析实例.docx(14页珍藏版)》请在冰豆网上搜索。
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:
2
forn2=1:
2
KK(DOF(n1),DOF(n2))=KK(DOF(n1),DOF(n2))+k(n1,n2);
end
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各单元的长度及轴线方向余弦
节点
x
y
单元
节点1
节点2
单元
l
xn
yn
1
0
0
①
1
2
①
400
1
0
2
400
0
②
3
2
②
300
0
-1
3
400
300
③
1
3
③
500
0.8
0.6
4
0
300
④
4
3
④
400
1
0
(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:
4
forn2=1:
4
KK(DOF(n1),DOF(n2))=KK(DOF(n1),DOF(n2))+k(n1,n2);
end
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.4;
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进行刚度矩阵的组装。
相关的计算流程如下。
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)
(4)边界条件的处理及刚度方程求解
由图可以看出,节点1的位移将为零,即u1=0,v1=0,节点2的位移v2=0,节点4的u4=0,v4=0。
节点载荷F3=10N。
采用高斯消去法进行求解,注意:
MATLAB中的反斜线符号“\”就是采用高斯消去法。
该结构的节点位移为:
而节点力为:
k=KK([3,5,6],[3,5,6])
p=[20000;0;-25000];
u=k\p
结果与前面通过数学推导得到的相同
(5)支反力的计算
在得到整个结构的节点位移后,由原整体刚度方程就可以计算出对应的支反力。
将整体的位移列阵q(采用国际标准单位)代回原整体刚度方程,计算出所有的节点力P,按上面的对应关系就可以找到对应的支反力。
相关的计算流程如下。
q=[000.000271200.0000565-0.000222500]'
P=KK*q
(6)各单元的应力计算
先从整体位移列阵q中提取出单元的位移列阵,然后,调用计算单元应力的函数Bar2D2Node_Stress,就可以得到各个单元的应力分量。
当然也可以调用上面的Bar2D2Node_Forces(E,A,x1,y1,x2,y2,alpha,u)函数来计算单元的集中力,然后除以面积求得单元应力。
相关的计算流程如下。
u1=[q
(1);q
(2);q(3);q(4)]
stress1=Bar2D2Node_Stress(E,x1,y1,x2,y2,alpha1,u1)
u2=[q(3);q(4);q(5);q(6)]
stress2=Bar2D2Node_Stress(E,x2,y2,x3,y3,alpha2,u2)
u3=[q
(1);q
(2);q(5);q(6)]
stress3=Bar2D2Node_Stress(E,x1,y1,x3,y3,alpha3,u3)
u4=[q(7);q(8);q(5);q(6)]
stress4=Bar2D2Node_Stress(E,x4,y4,x3,y3,alpha1,u4)
计算结果与前面通过数学推导得到的结果相同
【四杆桁架结构的有限元分析—ANSYS】
ANSYS是大型的通用有限元分析系统,ANSYS操作流程,包括基于图形界面的操作以及基于命令流的操作。
这样将使得以基于详细推导的典型例题与基于MATLAB的编程实现、以及与基于ANSYS的分析都完整地结合起来,可以更好的理解和使用有限元方法这一工具。
1基于ANSYS图形界面(GUI,graphicuserinterface)的菜单操作流程
2完整的命令流
以下为命令流语句;注意:
以“!
”打头的文字为注释内容,其后的文字和符号不起运行作用。
!
%%%%%%%%%%%%begin%%%%%%
/PREP7!
进入前处理
/PLOPTS,DATE,0!
设置不显示日期和时间
!
=====设置单元、材料,生成节点及单元
ET,1,LINK1!
选择单元类型
UIMP,1,EX,,,2.95e11,!
给出材料的弹性模量
R,1,1e-4,!
给出实常数(横截面积)
N,1,0,0,0,!
生成1号节点,坐标(0,0,0)
N,2,0.4,0,0,!
生成2号节点,坐标(0.4,0,0)
N,3,0.4,0.3,0,!
生成3号节点,坐标(0.4,0.3,0)
N,4,0,0.3,0,!
生成4号节点,坐标(0,0.3,0)
E,1,2!
生成1号单元(连接1号节点和2号节点)
E,2,3!
生成2号单元(连接2号节点和3号节点)
E,1,3!
生成3号单元(连接1号节点和3号节点)
E,4,3!
生成4号单元(连接4号节点和3号节点)
FINISH!
前处理结束
!
=====在求解模块中,施加位移约束、外力,进行求解
/SOLU!
进入求解状态(在该状态可以施加约束及外力)
D,1,ALL!
将1号节点的位移全部固定
D,2,UY,!
将2号节点的y方向位移固定
D,4,ALL!
将4号节点的位移全部固定
F,2,FX,20000,!
在2号节点处施加x方向的力(20000)
F,3,FY,-25000,!
在3号节点处施加y方向的力(-25000)
SOLVE!
进行求解
FINISH!
结束求解状态
!
=====进入一般的后处理模块
/POST1!
进入后处理
PLDISP,1!
显示变形状况
FINISH!
结束后处理
!
%%%%%%%%%%%%end%%%%%%
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 桁架 结构 有限元分析 实例