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

实例：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 % function k=Bar1D2Node_Stiffness(E, A, L) %该函数计算单元的刚度矩阵 %输入弹性模量E，横截面积A和长度L %输出单元刚度矩阵k(2×2) %- k=E*A/L -E*A/L; -E*A/L E*A/L; % function z=Bar1D2Node_Assembly(KK,k,i,j) %该函数进行单元刚度矩阵的组装 %输入单元刚度矩阵k，单元的节点编号i、j %输出整体刚度矩阵KK %- DOF(1)=i; DOF(2)=j; for n1=1:2 for n2=1:2 KK(DOF(n1), DOF(n2)= KK(DOF(n1), DOF(n2)+k(n1, n2); end end z=KK; %- function stress=Bar1D2Node_Stress(k, u, A) %该函数计算单元的应力 %输入单元刚度矩阵k, 单元的位移列阵u(2×1) %输入横截面积A计算单元应力矢量 %输出单元应力stress %- stress=k*u/A;%- % function forces=Bar1D2Node_Force(k, u) %该函数计算单元节点力矢量 %输入单元刚度矩阵k和单元的位移列阵u(2×1) %输出2×1的单元节点力分量forces %- forces=k*u; % Bar1D2Node % end %【四杆桁架结构的有限元分析数学推导】 如图所示的结构，各杆的弹性模量和横截面积都为E=29.54×10N/mm2, A=100mm 2，试求解该结构的节点位移、单元应力以及支反力。 图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 = u1 v1 u2 v2 u3 v3 u4 v4T 节点力： P=R+F= Rx1 Ry1 2×104 Ry2 0 2.5×104 Rx4 Ry4T其中(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 % function k=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*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; % function z = 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; 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; %- function stress= 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 -S C S*u; % function forces= 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 -S C S*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,角度alpha 1, alpha 2和alpha 3,然后分别针对单元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。节点载荷F310N。采用高斯消去法进行求解，注意：MATLAB中的反斜线符号“”就是采用高斯消去法。该结构的节点位移为：而节点力为：k=KK(3,5,6,3,5,6) p=20000;0;-25000; u=kp 结果与前面通过数学推导得到的相同（5）支反力的计算 在得到整个结构的节点位移后，由原整体刚度方程就可以计算出对应的支反力。将整体的位移列阵q（采用国际标准单位）代回原整体刚度方程，计算出所有的节点力P，按上面的对应关系就可以找到对应的支反力。相关的计算流程如下。 q=0 0 0.0002712 0 0.0000565 -0.0002225 0 0' 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, graphic user interface)的菜单操作流程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 %