平面应力应变有限元分析程序设计课程设计.doc

平面应力/应变有限元分析程序设计 引言 有限元法是求解微分方程边值问题的一种通用数值方法，该方法是一种基于变分法（或变分里兹法）而发展起来的求解微分方程的数值计算方法，以计算机为手段，采用分片近似，进而逼近整体的研究思想求解物理问题。 摘要 有限元法的基本思想是将物体（即连续的求解域）离散成有限个且按一定方式相互联结在一起的单元的组合，来模拟或逼近原来的物体，从而将一个连续的无限自由度问题简化为离散的有限自由度问题求解的一种数值分析方法。物体被离散后，通过对其中各个单元进行单元分析，最终得到对整个物体的分析。网格划分中每一个小的块体称为单元。确定单元形状、单元之间相互联结的点称为节点。单元上节点处的结构内力为节点力，外力（有集中力、分布力等）为节点载荷。 基本思想：单元结点按右手法则顺序编号。 荷载类型：可计算结点荷载。 说明：主程序的作用是通过赋值语句、读取和写入文件、函数调用等完成算法的全过程，即实现程序流程图的程序表达。 1.1 程序流程图 根据题意将得到的程序运行所需的初始数据的input.txt导入。由MATLAB中的LinearTriangleElementStiffness计算单刚矩阵和LinearTriangleAssemble 依次加入总刚矩阵。根据问题将约束信息引入 计算得到载荷向量，由F=KU计算支反力，再次调用MATLAB中的LinearTriangleElementStresses 和 LinearTriangleElementPStresses计算得到到单元应力和主应力。详细说明见附录。 1.2 程序应用举例 【例1】 如图所示，受均匀分布载荷作用的薄平板结构。将平板离散化成两个线性三角单元，如图所示。假定E=210Gpa，，t=0.025m，w=3000kn/m。 解答 如图所示 单元连通性如下表所示 单元编号 节点i 节点j 节点m 1 1 3 4 2 1 2 3 由题意可得输入数据文件input.txt为： 4 2 2 4 210e6 0.3 0.025 1 3 4 1 2 3 0.0 0.0 0.5 0.0 0.5 0.25 0.0 0.25 2 9.375 0 3 9.375 0 1 1 1 2 0 0 3 0 0 4 1 1 说明： 第一行：读入程序控制信息 NPION=fscanf(FP1,%d,1) %结点个数（结点编码总数） NELEM=fscanf(FP1,%d,1) %单元个数（单元编码总数） NFORCE=fscanf(FP1,%d,1) ％结点荷载个数 NVFIX=fscanf(FP1,%d,1) ％受约束边界点数 YOUNG=fscanf(FP1,%e,1) ％弹性模量 POISS=fscanf(FP1,%f,1) ％泊松比 THICK=fscanf(FP1,%d,1) %厚度 第二、三行：读入单元连接信息： LNODS=fscanf(FP1,%d,[3,NELEM]); %单元定义数组，单元结点号，逆时针输入 第四、五、六、七行：读入结点坐标 COORD=fscanf(FP1,%f,[2，NPOIN]); %结点坐标值，第1列为x坐标值，第2列为y坐标值 第八、九行：读入结点荷载信息 FORCE=fscanf(FP1,%f,[3,NFORCE]); %结点号，X向结点荷载数值，Y向结点荷载数值（与坐标轴方向一致为正） 第十、十一、十二、十三行：读入零位移信息 FIXED=fscanf(FP1,%d,[3,inf]); %结点号，X向约束，Y向约束 下述两个例子思路相同 不再赘述 %----------------------------------------------------------------------------------------------------- 运行该程序得到输出文件output.txt如下所示 节点号 x位移 y位移 1 0+000 0+000 2 7-006 1-006 3 6-006 4-008 4 0+000 0+000 节点号 Fx Fy 1 -9+000