基于Matlab语言的按平面三角形单元划分的结构有限元程序设计模板

基于 Matlab 语言的按平面三角形单元 划分的结构有限元程序设计 专专 业 业 建筑与土木工程建筑与土木工程 班班 级 级 建工研建工研 12 212 2 姓姓 名 名 韩志强韩志强 学学 号 号 基于 Matlab 语言的按平面三角形单元划分 结构有限元程序设计 一 有限单元发及 Matlab 语言概述 1 有限单元法 随着现代工业 生产技术的发展 不断要求设计高质量 高水平的大型 复杂和精密的机械及工程结构 为此目的 人们必须预先通过有效的计算手段 确切的预测即将诞生的机械和工程结构 在未来工作时所发生的应力 应变和 位移因此 需要寻求一种简单而又精确的数值分析方法 有限单元法正是适应 这种要求而产生和发展起来的一种十分有效的数值计算方法 有限元法把一个复杂的结构分解成相对简单的 单元 各单元之间通过结 点相互连接 单元内的物理量由单元结点上的物理量按一定的假设内插得到 这样就把一个复杂结构从无限多个自由度简化为有限个单元组成的结构 我们 只要分析每个单元的力学特性 然后按照有限元法的规则把这些单元 拼装 成整体 就能够得到整体结构的力学特性 有限单元法基本步骤如下 1 结构离散 结构离散就是建立结构的有限元模型 又称为网格划分或单 元划分 即将结构离散为由有限个单元组成的有限元模型 在该步骤中 需要 根据结构的几何特性 载荷情况等确定单元体内任意一点的位移插值函数 2 单元分析 根据弹性力学的几何方程以及物理方程确定单元的刚度矩阵 3 整体分析 把各个单元按原来的结构重新连接起来 并在单元刚度矩阵 的基础上确定结构的总刚度矩阵 形成如下式所示的整体有限元线性方程 KF 式中 是载荷矩阵 是整体结构的刚度矩阵 是节点位移矩阵 F K 4 载荷移置 根据静力等效原理 将载荷移置到相应的节点上 形成节点 载荷矩阵 5 边界条件处理 对式 所示的有限元线性方程进行边界条件处理 6 求解线性方程 求解式 所示的有限元线性方程 得到节点的位移 在 该步骤中 若有限元模型的节点越多 则线性方程的数量就越多 随之有限元 分析的计算量也将越大 7 求解单元应力及应变 根据求出的节点位移求解单元的应力和应变 8 结果处理与显示 进入有限元分析的后处理部分 对计算出来的结果进 行加工处理 并以各种形式将计算结果显示出 2 Matlab 简介 在用有限元法进行结构分析时 将会遇到大量的数值计算 因而在实用上 是一定要借助于计算机和有限元程序 才能完成这些复杂而繁重的数值计算工 作 而 MatlabMatlab 是当今国际科学界最具影响力和活力的软件 它起源于矩阵运算 并已经发展成一种高度集成的计算机语言 它提供了强大的科学计算 灵活的 程序设计流程 高质量的图形可视化与界面设计 便捷的与其他程序和语言接 口的功能 MatlabMatlab 在各国高校与研究单位起着重大的作用 工欲善其事 必先利其器 如果有一种十分有效的工具能解决在教学与 研究中遇到的问题 那么 MatlabMatlab 语言正是这样的一种工具 它可以将使用者从 繁琐 无谓的底层编程中解放出来 把有限的宝贵时间更多地花在解决问题中 这样无疑会提高工作效率 目前 MatlabMatlab 已经成为国际上最流行的科学与工 程计算的软件工具 现在的 MatlabMatlab 已经不仅仅是一个 矩阵实验室 了 它已 经成为了一种具有广泛应用前景的全新的计算机高级编程语言了 有人称它为 第四代 计算机语言 它在国内外高校和研究部门正扮演着重要的角色 MatlabMatlab 语言的功能也越来越强大 不断适应新的要求提出新的解决方法 可以 预见 在科学运算 自动控制与科学绘图领域 MatlabMatlab 语言将长期保持其独一无 二的地位 为此 本例采用 MatlabMatlab 语言编程 以利用其快捷强大的矩阵数值计 算功能 二 问题描述 一矩形薄板 一边固定 承受 150kN 集中荷载 结构简图及按平面三角形 单元划分的有限元模型图如下所示 材料参数 弹性模量 泊松比 薄板厚度 28 102mKNE 2 0 mm2 在本例中 所取结构模型及数据主要用于程序设计理论分析 与工程实际无关 1 4 y F 150kN 4 3 2 123 56 2m 1m1m 1m 1m x 三 参数输入 单元个数 NELEM 4 节点个数 NNODE 6 受约束边界点数 NVFIX 2 节点荷载个数 NFORCE 1 弹性模量 YOUNG 2e8 泊松比 POISS 0 2 厚度 THICK 0 002 单元节点编码数组 LNODS 6 5 2 5 4 2 4 3 2 6 2 1 节点坐标数组 COORD 1 0 1 1 1 2 0 2 0 1 0 0 节点力数组 FORCE 4 0 150 约束信息数组 FIXED 1 1 6 1 1 1 以上数值数据为程序运行前输入的初始数据 存为 txt txt 文本格式 同 时必须放在 MatlabMatlab 工作目录下 路径不对程序不能自动读取指定初始文件 运 行出错 初始数据文本文件输入格式如下图 四 Matlab 语言程序源代码 1 程序中变量说明 NNODE 单元节点数 NPION 总结点数 NELEM 单元数 NVFIX 受约束边界点数 FIXED 约束信息数组 NFORCE 节点力数 FORCE 节点力数组 COORD 结构节点坐标数组 LNODS 单元定义数组 YOUNG 弹性模量 POISS 泊松比 THICK 厚度 B 单元应变矩阵 3 6 D 单元弹性矩阵 3 3 S 单元应力矩阵 3 6 A 单元面积 ESTIF 单元刚度矩阵 ASTIF 总体刚度矩阵 ASLOD 总体荷载向量 ASDISP 节点位移向量 ELEDISP 单元节点位移向量 STRESS 单元应力 FP1 数据文件指针 2 Matlab 语言程序代码 初始化及数据调用 format short e 设定输出类型 clear 清除内存变量 FP1 fopen txt txt rt 打开输入数据文件 读入控制数据 NELEM fscanf FP1 d 1 单元个数 单元编码总数 NPION fscanf FP1 d 1 结点个数 结点编码总数 NVFIX fscanf FP1 d 1 受约束边界点数 NFORCE fscanf FP1 d 1 结点荷载个数 YOUNG fscanf FP1 e 1 弹性模量 POISS fscanf FP1 f 1 泊松比 THICK fscanf FP1 f 1 厚度 LNODS fscanf FP1 d 3 NELEM 单元定义数组 单元结点号 相应为单元结点号 编码 按逆时针顺序输入 COORD fscanf FP1 f 2 NPION 结点坐标数组 坐标 x y 坐标 共 NPOIN 组 FORCE fscanf FP1 f 3 NFORCE 结点力数组 受力结点编号 x 方向 y 方向 FIXED fscanf FP1 d 3 NVFIX 约束信息 约束点 x 约束 y 约束 有约束为 1 无约束为 0 生成单元刚度矩阵并组成总体刚度矩阵 ASTIF zeros 2 NPION 2 NPION 生成特定大小总体刚度矩阵并置 0 for i 1 NELEM 生成弹性矩阵 D D 1 POISS 0 POISS 1 0 0 0 1 POISS 2 YOUNG 1 POISS 2 计算当前单元的面积 A A det 1 COORD LNODS i 1 1 COORD LNODS i 1 2 1 COORD LNODS i 2 1 COORD LNODS i 2 2 1 COORD LNODS i 3 1 COORD LNODS i 3 2 2 生成应变矩阵 B for j 0 2 b j 1 COORD LNODS i rem j 1 3 1 2 COORD LNODS i rem j 2 3 1 2 c j 1 COORD LNODS i rem j 1 3 1 1 COORD LNODS i rem j 2 3 1 1 end B b 1 0 b 2 0 b 3 0 0 c 1 0 c 2 0 c 3 c 1 b 1 c 2 b 2 c 3 b 3 2 A B1 i B 求应力矩阵 S D B S D B ESTIF B S THICK A 求解单元刚度矩阵 a LNODS i 临时向量 用来记录当前单元的节点编号 for j 1 3 for k 1 3 ASTIF a j 2 1 a j 2 a k 2 1 a k 2 ASTIF a j 2 1 a j 2 a k 2 1 a k 2 ESTIF j 2 1 j 2 k 2 1 k 2 根据节点编号对应关系将单元刚度分块叠加到总刚 度矩阵中 end end end 将约束信息加入总体刚度矩阵 对角元素改一法 for i 1 NVFIX if FIXED i 2 1 ASTIF FIXED i 1 2 1 0 一列为零 ASTIF FIXED i 1 2 1 0 一行为零 ASTIF FIXED i 1 2 1 FIXED i 1 2 1 1 对角元素为 1 end if FIXED i 3 1 ASTIF FIXED i 1 2 0 一列为零 ASTIF FIXED i 1 2 0 一行为零 ASTIF FIXED i 1 2 FIXED i 1 2 1 对角元素为 1 end end 生成荷载向量 ASLOD 1 2 NPION 0 总体荷载向量置零 for i 1 NFORCE ASLOD FORCE i 1 2 1 FORCE i 1 2 FORCE i 2 3 end 求解内力 ASDISP ASTIF ASLOD 计算节点位移向量 ELEDISP 1 6 0 当前单元节点位移向量 for i 1 NELEM for j 1 3 ELEDISP j 2 1 j 2 ASDISP LNODS i j 2 1 LNODS i j 2 取出当前单元的节点位移向量 end i STRESS D B1 i ELEDISP 求内力 end fclose FP1 关闭数据文件 五 程序运行结果 NELEM 4 NPION 6 NVFIX 2 NFORCE 1 YOUNG POISS 2 0000e 001 THICK 2 0000e 003 LNODS 1 2 6 2 3 4 2 4 5 2 5 6 COORD 0 0 1 0 2 0 2 1 1 1 0 1 FORCE 4 0 150 FIXED 1 1 1 6 1 1 D 2 0833e 008 4 1667e 007 0 4 1667e 007 2 0833e 008 0 0 0 8 3333e 007 A 5 0000e 001 D 2 0833e 008 4 1667e 007 0 4 1667e 007 2 0833e 008 0 0 0 8 3333e 007 A 5 0000e 001 D 2 0833e 008 4 1667e 007 0 4 1667e 007 2 0833e 008 0 0 0 8 3333e 007 A 5 0000e 001 D 2 0833e 008 4 1667e 007 0 4 1667e 007 2 0833e 008 0 0 0 8 3333e 007 A 5 0000e 001 ASDISP 0 0 8 0607e 004 1 5848e 003 1 0281e 003 4 4727e 003 1 1937e 003 4 6947e 003 8 6670e 004 1 8880e 003 0 0 说明 以上 12 个 ASDISP 输出结果数据表示节点位移向量 即各节点分别在 X Y 方向上产生的位移 i 1 STRESS 1 6793e 005 3 3586e 004 1 3207e 005 i 2 STRESS 5 5503e 004 5 5503e 004 5 5503e 004 i 3 STRESS 5 5503e 004 4 9526e 004 9 4497e 004 i 4 STRESS 1 6793e 005 2 7040e 004 1 7932e 004 说明 以上四组 STRESS 输出结果数据表示单元应力 SX SY SXY i 为单元号 六 ANSYS 建模比较分析 利用 ANSYSANSYS 有限元分析软件 完全按照前面 MatlabMatlab 程序设计的各项条件参 数以及单元划分方式建立 ANSYSANSYS 模型 其求解结果与以上程序计算结果比较 验证程序是否正确 ANSYS 模型节点单元如下图所示 ANSYS 模型求解变形图如下所示 ANSYS 求解节点位移结果列表显示如下 单位 m PRINT U NODAL SOLUTION PER NODE POST1 NODAL DEGREE OF FREEDOM LISTING THE FOLLOWING DEGREE OF FREEDOM RESULTS ARE IN THE GLOBAL COORDINATE SYSTEM NODE UX UY UZ USUM 1 0 0000 0 0000 0 0000 0 0000 2 0 80607E 03 0 15848E 02 0 0000 0 17780E 02 3 0 10281E 02 0 44727E 02 0 0000 0 45893E 02 4 0 11937E 02 0 46947E 02 0 0000 0 48441E 02 5 0 86670E 03 0 18880E 02 0 0000 0 20774E 02 6 0 0000 0 0000 0 0000 0 0000 MAXIMUM ABSOLUTE VALUES NODE 4 4 0 4 VALUE 0 11937E 02 0 46947E 02 0 0000 0 48441E 02 ANSYS 求解单元应力结果列表显示如下 单位 KN m2 PRINT S ELEMENT SOLUTION PER ELEMENT POST1 ELEMENT NODAL STRESS LISTING THE FOLLOWING X Y Z VALUES ARE IN GLOBAL COORDINATES ELEMENT 1 PLANE182 NODE SX SY SZ SXY SYZ 1 0 16793E 06 33586 0 0000 0 13207E 06 0 0000 2 0 16793E 06 33586 0 0000 0 13207E 06 0 0000 6 0 16793E 06 33586 0 0000 0 13207E 06 0 0000 ELEMENT 2 PLANE182 NODE SX SY SZ SXY SYZ 2 55503 55503 0 0000 55503 0 0000 3 55503 55503 0 0000 55503 0 0000 4 55503 55503 0 0000 55503 0 0000 ELEMENT 3 PLANE182 NODE SX SY SZ SXY SYZ 2 55503 49526 0 0000 94497 0 0000 4 55503 49526 0 0000 94497 0 0000 5 55503 49526 0 0000 94497 0 0000 ELEMENT 4 PLANE182 NODE SX SY SZ SXY SYZ 2 0 16793E 06 27040 0 0000 17932 0 0000 5 0 16793E 06 27040 0 0000 17932 0 0000 6 0 16793E 06 27040 0 0000 17932 0 0000 结论 通过比较 MatlabMatlab 语言设计程序运行结果和 ANSYSANSYS 建模分析结果可知 两种 方式算出的结果完全一致 说程序设计正确 所以本程序适用于按三角形单元 划分的平面结构有限元分析 format short e clear FP1 fopen LinearTriangleElement of Thin plate txt rt NELEM fscanf FP1 d 1 NPION fscanf FP1 d 1 NVFIX fscanf FP1 d 1 NFORCE fscanf FP1 d 1 YOUNG fscanf FP1 e 1 POISS fscanf FP1 f 1 THICK fscanf FP1 f 1 LNODS fscanf FP1 d 3 NELEM COORD fscanf FP1 f 2 NPION FORCE fscanf FP1 f 3 NFORCE FIXED fscanf FP1 d 3 NVFIX ASTIF zeros 2 NPION 2 NPION 生成特定大小总体刚度矩阵并置 0 for i 1 NELEM 生成弹性矩阵 D D YOUNG 1 POISS 2 1 POISS 0 POISS 1 0 0 0 1 POISS 2 计算当前单元的面积 A A det 1 COORD LNODS i 1 1 COORD LNODS i 1 2 1 COORD LNODS i 2 1 COORD LNODS i 2 2 1 COORD LNODS i 3 1 COORD LNODS i 3 2 2 生成应变矩阵 B for j 0 2 b j 1 COORD LNODS i rem j 1 3 1 2 COORD LNODS i rem j 2 3 1 2 c j 1 COORD LNODS i rem j 1 3 1 1 COORD LNODS i rem j 2 3 1 1 end B b 1 0 b 2 0 b 3 0 0 c 1 0 c 2 0 c 3 c 1 b 1 c 2 b 2 c 3 b 3 2 A B1 i B 求应力矩阵 S D B S D B ESTIF B S THICK A a LNODS i for j 1 3 for k 1 3 ASTIF a j 2 1 a j 2 a k 2 1 a k 2 ASTIF a j 2 1 a j 2 a k 2 1 a k 2 ESTIF j 2 1 j 2 k 2 1 k 2 end end end 将约束信息加入总体刚度矩阵 for i 1 NVFIX if FIXED i 2 1 ASTIF FIXED i 1 2 1 0 ASTIF FIXED i 1 2 1 0 ASTIF FIXED i 1 2 1 FIXED i 1 2 1 1 end if FIXED i 3 1 ASTIF FIXED i 1 2 0 ASTIF FIXED i 1 2 0 ASTIF FIXED i 1 2 FIXED i 1 2 1 end end 生成荷载向量 ASLOD 1 2 NPION 0 for i 1 NFORCE ASLOD FORCE i 1 2 1 FORCE i 1 2 FORCE i 2 3 end 求解内力 ASDISP ASTIF ASLOD ELEDISP 1 6 0 for i 1 NELEM for j 1 3 ELEDISP j 2 1 j 2 ASDISP LNODS i j 2 1 LNODS i j 2 end i STRESS D B1 i ELEDISP end fclose FP1 format short e clear FP1 fopen LinearTriangleElement of Thin plate txt rt NELEM fscanf FP1 d 1 NPION fscanf FP1 d 1 NVFIX fscanf FP1 d 1 NFORCE fscanf FP1 d 1 YOUNG fscanf FP1 e 1 POISS fscanf FP1 f 1 THICK fscanf FP1 f 1 LNODS fscanf FP1 d NELEM 3 COORD fscanf FP1 f NPION 2 FORCE fscanf FP1 f NFORCE 3 FIXED fscanf FP1 d NVFIX 3 ASTIF zeros 2 NPION 2 NPION 生成特定大小总体刚度矩阵并置 0 for i 1 NELEM 生成弹性矩阵 D D YOUNG 1 POISS 2 1 POISS 0 POISS 1 0 0 0 1 POISS 2 计算当前单元的面积 A A det 1 COORD LNODS