有限元数值分析课程论文

《有限元数值分析》课程论文要求1．简述平面问题有限元法（自选单元类型）。2．写出计算程序框图。3．自编一份平面问题有限元弹性应力计算程序。4．以图1为算例，，分析计算结果。图11.简述平面问题有限元法1.1有限元法的基本思想：1.1.1假想把连续系统分割成数目有限的单元，单元之间只在数目有限的指定点处（称为节点）相互连接，构成一个单元集合体来代替原来的连续系统。在节点上引进等效载荷，代替实际作用于系统上的外载荷。1.1.2对每个单元由分块挖的思想，按一定的规则建立求解未知量与节点相互作用之间的关系。1.1.3把所有单元的这种特性关系按一定的条件集合起来，引入边界条件，构成一组以节点变量为未知量的代数方程组，求解之就得到有限个节点处的待求变量。1.2平面问题包括平面应力和平面应变问题1.2.1平面应力问题平面应力问题研究等厚度薄板状弹性体，厚度尺寸远远小于截面尺寸，受力方向沿板面方向，不沿厚度变化。1.2.2平面应变问题平面应变问题处理面内受力但垂直于平面方向上不产生变形的二维受力问题，沿截面方向的截面形状和大小相同且厚度尺寸远远大于截面尺寸，载荷垂直于厚度方向且沿厚度均匀分布。1.3结合所给模型说明平面问题的三角形单元求解1.3.1结构的离散化将所给实际模型简化，并划分为如图所示的2个三节点三角形单元，并对其进行编码，节点编号和坐标、单元编号以及相关节点如表1、表2所示。图2表1节点编码及其坐标节点编号1234x05005000y00200200表2单元编码及其相关节点单元编号12i14j21m33表3节点及节点载荷节点编号1234FxF1x0F3xF4xFyF1y0F3yF4y表4节点及节点位移节点编号1234u0u2u30v0v2v30节点载荷列阵(1)节点位移列阵：（2）1.3.2分别计算单元的刚度矩阵将每个单元局部编码，都可看做如图4的形式，每个单元都有其对应的i、j、m。图3单元刚度矩阵（t为板料厚度）（3）为几何矩阵，与节点坐标有关，其计算公式如下：（4）其中单元面积（5）为弹性矩阵，与材料的弹性模量E以及泊松比μ有关，本问题按照平面应力问题求解，平面应力问题的弹性矩阵的表达式如式（7）所示：（6）其中在所划分的模型中，这2个单元的面积是相同的。因为本问题是简单的平面应力问题，在求解几何矩阵时，没有用到位移插值函数，而是直接进行几何矩阵的求解。首先求出节点坐标的差和单元面积，得到几何矩阵，并将弹性模量E和泊松比μ带入式（6）得到弹性矩阵，最后根据公式（3）求得单元的刚度矩阵。这里有2个单元，因此将有2个单元刚度矩阵。1.3.3合成总体刚度矩阵当单元进行局部编码时（如图4），每个单元的刚度矩阵都有如下形式;（7）按单元的位移自由度所对应的位置进行组装可以得到整体刚度矩阵1.3.4引入边界条件求解刚度方程将所得到的总刚度矩阵、节点载荷列阵（1）式以及节点位移列阵（2）式代入整体刚度方程中，总体刚度方程如式（8）。（8）将位移边界条件和力边界位移条件代入总体刚度方程，对总体刚度矩阵进行降阶，划去的对应位移为0的行和列。总体刚度矩阵方程降阶后为线性方程组，采用牛顿迭代法求解，得到节点位移阵列，再将节点位移代入总体刚度矩阵方程中，可求得节点载荷F。1.3.5求解单元应变和单元应力将每个单元的位移和几何矩阵代入式（9），求得每个单元的应变eeeB（9）其中mmjjiievuvuvu将每个单元的应变代入式（10），求得每个单元的应力eD（10）2、计算程序框图3、平面应力问题有限元弹性应力计算程序#includemath.h#includestdafx.h#includeiostream#includeFYMatrixMN.husingnamespacestd;#defined11(E,V)E/(1-V)/(1+V)#defined12(E,V)E*V/(1-V)/(1+V)#defined33(E,V)E/2/(1+V)typedefFYNAMESPACE::FYMatrixMNdoubleFyMatrix;intmain(intargc,char*argv[]){cout程序中采用的单位是MPa、mm、Nendl;doubleE=210000,V=0.3,t=2;doublex1=0,y1=0,x2,y2=0,x3,y3,x4=0,y4;coutendlendl请按1、2、3、4节点号顺序依次输入节点坐标x1,y1,x2,y2,x3,y3,x4,y4endl;cinx1y1x2y2x3y3x4y4;//计算单元1的单元刚度矩阵,S11是面积doubleS11;FyMatrixdanyuan1(6,6);FyMatrixB1(3,6);FyMatrixD1(3,3);FyMatrixH1(3,6);doublex32,x21,x13;doubley12,y23,y31;x32=x3-x2;x21=x2-x1;x13=x1-x3;y12=y1-y2;y23=y2-y3;y31=y3-y1;S11=((x2*y3-x3*y2)-(x1*y3-x3*y1)+(x1*y2-x2*y1))/2.0;B1.AllZero();B1(0,0)=y23;B1(0,2)=y31;B1(0,4)=y12;B1(1,1)=x32;B1(1,3)=x13;B1(1,5)=x21;B1(2,0)=x32;B1(2,1)=y23;B1(2,2)=x13;B1(2,3)=y31;B1(2,4)=x21;B1(2,5)=y12;D1.AllZero();D1(0,0)=d11(E,V);D1(0,1)=d12(E,V);D1(1,0)=D1(0,1);D1(1,1)=D1(0,0);D1(2,2)=d33(E,V);H1=D1*B1*(0.5/S11);danyuan1=B1.Transpose()*D1*B1*t*(0.25/S11);//计算单元2的单元刚度矩阵，S22是面积doubleS22;FyMatrixdanyuan2(6,6);FyMatrixB2(3,6);FyMatrixD2(3,3);FyMatrixH2(3,6);doublex43,x31,x14;doubley13,y34,y41;x43=x4-x3;x31=x3-x1;x14=x1-x4;y13=y1-y3;y34=y3-y4;y41=y4-y1;S22=S11;B2.AllZero();B2(0,0)=y13;B2(0,2)=y34;B2(0,4)=y41;B2(1,1)=x31;B2(1,3)=x43;B2(1,5)=x14;B2(2,0)=x31;B2(2,1)=y13;B2(2,2)=x43;B2(2,3)=y34;B2(2,4)=x14;B2(2,5)=y41;D2=D1;H2=D2*B2*(0.5/S22);danyuan2=B2.Transpose()*D2*B2*t*(0.25/S22);//由单元刚度矩阵合成总体刚度矩阵FyMatrixzongti(8,8);FyMatrixK111(2,2),K121(2,2),K122(2,2),K131(2,2),K132(2,2),K133(2,2);FyMatrixK244(2,2),K214(2,2),K211(2,2),K234(2,2),K231(2,2),K233(2,2),K241(2,2),K243(2,2);//分解单元1的刚度矩阵K111(0,0)=danyuan1(0,0);K111(0,1)=danyuan1(0,1);K111(1,0)=danyuan1(1,0);K111(1,1)=danyuan1(1,1);K121(0,0)=danyuan1(2,0);K121(0,1)=danyuan1(2,1);K121(1,0)=danyuan1(3,0);K121(1,1)=danyuan1(3,1);K122(0,0)=danyuan1(2,2);K122(0,1)=danyuan1(2,3);K122(1,0)=danyuan1(3,2);K122(1,1)=danyuan1(3,3);K131(0,0)=danyuan1(4,0);K131(0,1)=danyuan1(4,1);K131(1,0)=danyuan1(5,0);K131(1,1)=danyuan1(5,1);K132(0,0)=danyuan1(4,2);K132(0,1)=danyuan1(4,3);K132(1,0)=danyuan1(5,2);K132(1,1)=danyuan1(5,3);K133(0,0)=danyuan1(4,4);K133(0,1)=danyuan1(4,5);K133(1,0)=danyuan1(5,4);K133(1,1)=danyuan1(5,5);//分解单元2的刚度矩阵K244(0,0)=danyuan2(0,0);K244(0,1)=danyuan2(0,1);K244(1,0)=danyuan2(1,0);K244(1,1)=danyuan2(1,1);K214(0,0)=danyuan2(2,0);K214(0,1)=danyuan2(2,1);K214(1,0)=danyuan2(3,0);K214(1,1)=danyuan2(3,1);K241=K214;K211(0,0)=danyuan2(2,2);K211(0,1)=danyuan2(2,3);K211(1,0)=danyuan2(3,2);K211(1,1)=danyuan2(3,3);K234(0,0)=danyuan2(4,0);K234(0,1)=danyuan2(4,1);K234(1,0)=danyuan2(5,0);K234(1,1)=danyuan2(5,1);K243=K234;K231(0,0)=danyuan2(4,2);K231(0,1)=danyuan2(4,3);K231(1,0)=danyuan2(5,2);K231(1,1)=danyuan2(5,3);K233(0,0)=danyuan2(4,4);K233(0,1)=danyuan2(4,5);K233(1,0)=danyuan2(5,4);K233(1,1)=danyuan2(5,5);//合成总体刚度矩阵zongti.AllZero();zongti(0,0)=K111(0,0)+K211(0,0);zongti(1,0)=K111(1,0)+K211(1,0);zongti(1,1)=K111(1,1)+K211(1,1);zongti(2,0)=K121(0,0);zongti(2,1)=K121(0,1);zongti(3,0)=K121(1,0);zongti(3,1)=K121(1,1);zongti(2,2)=K122(0,0);zongti(3,2)=K122(1,0);zongti(3,3)=K122(1,1);zongti(4,0)=K131(0,0)+K231(0,0);zongti(4,1)=K131(0,1)+K231(0,1);zongti(5,0)=K131(1,0)+K231(1,0);zongti(5,1)=K131(1,1)+K231(1,1);zongti(4,2)=K132(0,0);zongti(4,3)=K132(0,1);zongti(5,2)=K132(1,0);zongti(5,3)=K132(1,1);zongti(4,4)=K133(0,0)+K233(0,0);zongti(5,4)=K133(1,0)+K233(1,0);zongti(5,5)=K133(1,1)+K233(1,1);zongti(6,0)=K241(0,0);zongti(6,1)=K241(0,1);zongti(7,0)=K241(1,0);zongti(7,1)=