有限元分析作业1

《有限元分析》课程作业任课教师：徐亚兰学生姓名：林声佳学号：13040120093班级：1304012时间：2016-01-10一、问题描述及分析如图1所示为一矩形薄平板，在右端部受分布力P=10KN作用，材料常数为：弹性模量E=1×107Pa、泊松比μ=1/3，板的厚度为t=0.1m，试按平面应力问题计算各个节点位移及支座反力y13PP1m42x(a)问题描述(b)有限元分析模型图1右端部分受均布力作用的平面问题二、有限元建模及分析1、基于3节点三角形单元的有限元建模及分析（1）结构的离散化与编号对该结构进行离散，单位编号及节点编号如图（1-b）所示，即有两个3节点三角形单元。载荷F按静力等效原则向节点1，节点2移置等效。节点位移列阵q=[u1v1u2v2u3v3u4v4]T节点外载荷列阵F=[12Plt012Plt00000]T约束的支反力列阵R=F=[0000Rx3Ry3Rx4Ry4]T总的节点载荷列阵P=F+R=F=[12Plt012Plt0Rx3Ry3Rx4Ry4]T其中的Rx3、Ry3、Rx4、Ry4分别为节点3和节点4的两个方向的支反力（2）各单元的刚度矩阵及刚度方程a.单元的几何和节点描述1m②①图2如图2所示，一个单元体有6个节点位移自由度（DOF），将所有节点上的位移组成列阵q；同样，将所有节点上个力也组成列阵P，那么q=[u1v1u2v2u3v3]TP=[Px1Py1Px2Py2Px3Py3]T当两个单元取图a示中的局部编码（i,j,m)时，其各单元的位移场，应力场，应变场，势能，刚度矩阵完全相同。b.单元的位移场描述就如图2所示的平面3节点三角形单元，由于有3个节点，每一个节点有两个位移，因此共有6个节点位移，考虑到简单性、完备性、连续性及待定数的唯一确定性原则，分别选取单元中各个方向的位移模式为N(x,y)=123123000000NNNNNNNi=(a+bix+ciy),i=1,2,3其中a=x2y3-x3y2,b=y2-y3,c=-x2+x3上式中的符号（1，2，3）表示下标轮换，如1→2，2→3，3→1同时更换。u(x,y)=N（x,y)·qc.单元的应力场描述（x,y)=xxyyxy=00xyyx,(,)uxyvxy=[∂]u其中[∂]为几何方程的算子矩阵(operatormatrix)，即[∂]=00xyyx（x,y)=[∂]·N(x,y)·q=B(x,y)·qB=[B1B2B3]其中Bi=A21·00iiiibccb,i=1,2,3d.单元的应变场描述由弹性力学中平面问题的物理方程，将其写成矩阵形式2101011002xxxxyyyyxyxyuEuDuu其中平面应力问题的弹性系数矩阵D为D=2101011002uEuuu若为平面应变问题，则将上式中的系数(E,μ)换成平面应变问题的系数2,11E。DBqSqSDBe.单元的势能表达单元的势能的表达以上已将单元的三大基本变量(u,,εσ)用基于节点位移列阵来进行表达q来进行表达，将其带入单位势能的表达式，有12TTSdbdpudA12TTqKqpq其中K是单元刚度矩阵，即TTKBDBdBDBdAtt为平面问题的厚度，这时B为常系数矩阵，因此上式可以写成112233212223313233TKKKKBDBAtKKKKKK其中各个矩阵为1223441TrsrkkEtKBDBAtkk其中112rsrskbbcc212rsrskcbbc312srsrkcbbc412rssrkcbbbf.单元刚度阵及刚度方程将单元的势能式对节点位移取一阶极值,可得到单元的刚度方程qKqP（3）整体刚度阵及刚度方程a.单元贡献矩阵各个单元的描述当两个单元取图示中的局部编码(i,j,m)时，其单元刚度矩阵完全相同，即(1),(2)iiijimjijjjmmimjmmkkkKkkkkkk1002/312/301/32/302/31/302/34/304/32/392/30042/343212/34/32/37/34/32/31/32/344/313/3Etb.整体刚度阵及刚度方程按单元的位移自由度所对应的位置进行组装可以得到整体刚度矩阵，该组装过程可以写成(1)(2)KKK具体写出单元刚度矩阵的各个子块在总刚度矩阵的对应位置如下(2)(2)(2)(2)(2)(1)(2)(1)(1)(2)(2)(1)(2)(1)(1)88(1)(1)(1)mmmjmijmjjiijiijmjimjiijiijjjmmimjmmkkkkkkkkkKkkkkkkkkk带入整体刚度Kq=p中有，1122744221003333341322140033333427420013333321342140093333324742321003333321413200433333242740013333321241300433333uvuvEt333344440202xyxxFFuRvRuRvR3344120120xyxyPltPltRRRR（4)边界条件的处理及方程求解该问题的位移边界条件是u3=0,v3=0,u4=0,v4=0,将其带入上式中，划去已知节点位移对应的第5行至第8行（列），则11227442033334132493332427032033321324033uFvEtuvF120120PltPlt所以：节点位移u1=0.0020,v1=-0.0004,u2=0.0018,v2=-0.0000(5）结果分析支反力Rx3=-500.0000,Ry3=334.0757,Rx4=-500.0000,Ry4=-334.0757节点位移u1=0.0020,v1=-0.0004,u2=0.0018,v2=-0.00002、基于四节点四边形单元的有限元建模及分析3、（1）结构的离散化与编号3151mP4P1m642x(a)问题描述（b)有限元模型图3右端受均布力的薄板将结构离散为二个4节点矩形单元，单元编号及节点编号如图(3-b)所示，连接关系见表1，节点的几何坐标见表2，载荷F按静力等效原则向节点1、节点2移置。表1结构的单元连接关系单元号节点号1235641342表2节点的坐标节点节点坐标(m)xy123456212011100100节点位移列阵q=[u1v1u2v2u3v3u4v4u5v5u6v6]T①②节点外载荷列阵F=[12Plt012Plt00000000]T约束的支反力列阵R=F=[00000000Rx3Ry3Rx4Ry4]T总的节点载荷列阵P=F+R=F=[12Plt012Plt00000Rx3Ry3Rx4Ry4]T其中（Rx5,Ry5),（Rx6,Ry6)分别为节点5和节点6的两个方向的支反力。（2）局部坐标系下单元的描述a.单元的几何和节点描述图4平面四节点矩阵单元如图4所示平面4节点矩形单元，单元的节点位移共有8个自由(DOF)。节点的编号为1、2、3、4，各自的位置坐标为(xi,yi),i=1,2,3,4，各个节点的位移(分别沿x方向和y方向)为(ui,vi),i=1,2,3,4。b.单元的位移场描述若采取无量纲坐标xayb112233441,11,11,11,1将所有节点上的位移组成一个列阵q，记作；同样，将所有节点上的各个力也组成一个列阵p,那么qe=[u1v1u2v2u3v3u4v4]TP=[Px1Py1Px2Py2Px3Py3Px4Py4]Tc.单元的应力场描述单元位移场的表达从图2-b可以看出，节点条件共有8个，即x方向4个1234(,,,)uuuu，y方向4个1234(v,v,v,v)，因此，x和y方向的位移场可以各有4个待定系数，即取以多项式作为单元的位移场模式01230123u(x,y)av(x,y)axayaxybbxbybxy由节点条件，在x=xi，y=yi处，有u(x,y)u,1,2,3,4v(x,y)viiiiiii11212342123433440000u(x,y)(,)0000u(x,y)uvuNNNNvuxyNqNNNNuvuv其中1(1)(1),1,2,3,44IiiNid.单元的应变场描述由弹性力学平面问题的几何方程(矩阵形式)，有单元应变的表达Ɛ（x,y)=xxyyxy=00xyyx,(,)uxyvxy=[∂]u几何矩阵B（x,y)为B(x,y)=[]N=[B1B2B3B4]1,2,3,4iiiiiNxNBiyNNyxe.单元的势能表达单元势能表达式12TTqKqpq11212231323341424344eTkkksymKBDBtdAkkkkkkk1322441rskkEtkkkuab其中2212322411111323121211111,(,1,2,3,4)323rsrsrsrsrssrsrrsrsrsukbaukabuukabuukabrsf.单元刚度阵及刚度方程将单元的势能对节点位移取一阶极值，可得到单元的刚度方程KqP（3）整体刚度阵及刚度方程a.总体坐标系下刚度矩阵的计算（等参变换）112122(1),(2)31323341424344TkkksymKBDBtdAkkkkkkk1322441rskkEtkkkuab其中2212322411111323121211111,(,1,2,3,4)323rsrsrsrsrssrsrrsrsrsukbaukabuukabuukabrsb.整体刚度阵及刚度方程将针对节点1,2,3,4的位移进行求解，得到刚度方程K(1)；将针对节点3,4,5,6的位移进行求解，得到刚度方程K(2)(1)(2)KKK（4)边界条件的处理及方程求解节点5,6的两个方向的位移将为0，即55660,v0,0,v0uu(5）结果分析位移结果u1=0.0020,v1=-0.0002,u2=0.0020,v2=0.0002,u3=0.0010,v3=-