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

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

1、实例： 2D四杆桁架结构的有限元分析 学习有限元方法的一个最佳途径，就是在充分掌握基本概念的基础上亲自编写有限元分析程序，这就需要一个 良好的编程环境或平台。MATLAB软件就是这样一个平台，它以功能强大、编程逻辑直观、使用方便见长。将提 供有限元分析中主要单元完整的MATLAB 程序，并给出详细的说明。 1D杆单元的有限元分析MATLAB 程序 (Bar1D2Node) 最简单的线性杆单元的程序应该包括单元刚度矩阵、单元组装、单元应力等几个基本计算程序。下面给出编写 的线性杆单元的四个 MATLAB 函数。 Bar1D2Node _Stiffness(E,A,L) 该函数计算单元的刚度矩阵，

2、输入弹性模量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杆

3、单元的有限元分析的基本公式，写出具体实现以上每个函数的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; D

4、OF(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) %该函数计算单元节点力矢量 %输入单元刚度矩阵

5、 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

6、单元节点1 节点2 单元l xnyn 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 = u1v1u2v2u3

7、v3u4v4 T 节点力：P=R+F= Rx1Ry12 10 4 Ry20 2.5 10 4 Rx4Ry4 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) 支反力的计算 将节点位移的结果代入整体刚度方程中，可求出 2

8、D杆单元的有限元分析程序(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。 B

9、ar2D2Node_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 程序如下

10、。 % 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

11、*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; %- functi

12、on 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,

13、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所示结构，基于 MAT

14、LAB 平台求解该结构的节点位移、单元应力以及支反力。 解答：对该问题进行有限元分析的过程如下。 （1） 结构的离散化与编号 对该结构进行自然离散，节点编号和单元编号如图所示，有关节点和单元的信息见表1表3。 （2）计算各单元的刚度矩阵 (基于国际标准单位 ) 建立一个工作目录， 将所编制的用于平面桁架单元分析的四个MATLAB 函数放置于该工作目录中， 分别以各自 函数的名称给出文件名，即： Bar2D2Node_Stiffness， Bar2D2Node_Assembly， Bar2D2Node_Stress ， Bar2D2Node_Forces 。 然后启动 MATLAB ，将工作目录

15、设置到已建立的目录中，在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=B

16、ar2D2Node_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(88)，先对 KK 清零，然后四次调用函数 Bar2D2Node _Assembly

17、进行刚度矩阵的组装。相关的计算流程如下。 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

18、 中的反斜线符号 “ ” 就是采用高斯消去法。 该结构的节点位移为： 而节点力为： 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中提取出单元的位移列

19、阵，然后，调用计算单元应力的函数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

20、);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的分析都完整地结合起来， 可以更好的理解和使用有限

21、元方法这一工具。 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

22、.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 %