2D四杆桁架结构的有限元分析实例学习资料.docx
《2D四杆桁架结构的有限元分析实例学习资料.docx》由会员分享，可在线阅读，更多相关《2D四杆桁架结构的有限元分析实例学习资料.docx（15页珍藏版）》请在冰豆网上搜索。
2D四杆桁架结构的有限元分析实例学习资料
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%%%%%%
升级会员 