基于MATLAB平面桁架有限元分析研究_李远瑛

内力计算结果、 变形图和结构的动力特性 ， 并将计算结果与 ANSYS 有限元分析软件的结果进行比较 ． 结果表明， 二者的计算结果一致． 利用 MATLAB 可以编制结构有限元分析程序 ， 计算速度快， 数据准确， 并能为复杂结构有 限元分析提供参考． 关键词：MATLAB； 平面桁架； 有限元； ANSYS 中图分类号：TU311． 2 文献标识码： A 文章编号： 1006 － 642X（ 2012 ） 08 － 0029 － 05
图 3 基于 MATLAB 平面桁架静力有限元分析模拟流程图
theta10 ） ； % 生成整体坐标系下单元刚度矩阵 K = zeros（ 12 ， 12 ） ； K = PlaneTrussAssemble（ K ， k1 ， 1， 2） ； K = PlaneTrussAssemble（ K ， k2 ， 1， 3） ； K = PlaneTrussAssemble（ K ， k3 ， 2， 3） ； K = PlaneTrussAssemble（ K ， k4 ， 2， 4） ； K = PlaneTrussAssemble（ K ， k5 ， 3， 4） ； K = PlaneTrussAssemble（ K ， k6 ， 3， 5） ； K = PlaneTrussAssemble（ K ， k7 ， 4， 5） ； K = PlaneTrussAssemble（ K ， k8 ， 4， 6） ； K = PlaneTrussAssemble（ K ， k9 ， 5， 6） ； K = PlaneTrussAssemble（ K ， k10 ， 2， 5 ） ； % 生成总刚度 矩阵 K1 = ［ K （ 3 ： 11 ， 3 ： 11 ） ］ ； P =［ 0 ； － 200000 ； 0 ； 0 ； 0 ； － 200000 ； 0 ； 0 ； 0 ］ ； % 生成
始末端结 点编号 1 →2 1 →3 2 →3 2 →4 2 →5
K△ = P △为结构整体坐标下各结点位移向量 ． 1． 7 分析各杆件的单元轴力及应力［3］ F e = k e T△ e ； σ e = Fe A
（ 2） （ 3） （ 4） （ 5）
第 30 卷第 8 期
李远瑛， 张德生
基于 MATLAB 平面桁架有限元分析研究
2 ． 3 编制 MATLAB 平面桁架有限元静力分析程序 编制 MATLAB 平面桁架有限元静力分析程序． 程序中 PlaneTrussElementStiffness 为生成单元刚度矩 PlaneTrussAssemble 为集成结构总刚度 阵的子程序， 矩阵的子程序 clear； clc； E = 2 ． 0 E11 ； 杆件弹性模量 A = 2 ． 5 E － 3 ； % 杆件横截面积 L1 = 4 ； L4 = 4 ； L8 = 4 ； L6 = 4 ； L2 = 5 ； L5 = 5 ； L9 = 5 ； L3 = 3 ； L7 = 3 ； L10 = 5 ； % 计算各单元长度 theta1 = 0 ； theta2 = atan （ 3 / 4 ） * 180 / pi； theta3 = 90 ； theta4 = 0 ； theta5 = － atan（ 3 / 4 ） * 180 / pi； theta6 = 0 ； theta7 = 90 ；
e e
对于平面桁架单元， 通过坐标系旋转可以完成 从单元局部坐标系到结构整体坐标系的转换． 杆件 转换矩阵为 局部坐标系与整体坐标系的夹角为 α ，
Leabharlann Baidu
1 ． 3 整体坐标系下的杆件单元刚度矩阵 整体坐标系下桁架单元刚度矩阵与局部坐标
2
工程算例
某平面桁架如图 2 所示． 桁架的材料特性为：
e T e 即 系下的单元刚度矩阵关系为 k = T k T， 11 弹性模量 E = 2 ． 0 × 10 Pa： 所有单元的横截面积为
2 ． 5 × 10 － 3 m2 ． 在结点 2 和结点 4 施加竖向 200 kn 的集中力． 1． 4 整体坐标下结构总刚度矩阵 K 在结构的整体分析中涉及许多单元及联接它 们的结点， 为避免混淆， 必须对各单元、 结点及结点 位移（ 水平位移和竖向位移 ） 进行编号， 由单元刚度 矩 矩阵通过单元定位向量集成结构总刚度矩阵 K， 阵的阶数为结点总数的 2 倍． 单元定位向量由单元 两端的结点位移编号确定． 1． 5 建立结构的结点荷载向量 P
31
2 ． 2 编制程序流程图 编制 MATLAB 平面桁架静力有限元分析模拟 流程图， 见图 3 所示．
theta8 = 0 ； theta9 = － atan （ 3 / 4 ） * 180 / pi； theta10 = atan（ 3 / 4 ） * 180 / pi； % 计算各单元局部坐标与整体 坐标的夹角 k1 = PlaneTrussElementStiffness（ E ， A， L1 ， theta1 ） ； k2 = PlaneTrussElementStiffness（ E ， A， L2 ， theta2 ） ； k3 = PlaneTrussElementStiffness（ E ， A， L3 ， theta3 ） ； k4 = PlaneTrussElementStiffness（ E ， A， L4 ， theta4 ） ； k5 = PlaneTrussElementStiffness（ E ， A， L5 ， theta5 ） ； k6 = PlaneTrussElementStiffness（ E ， A， L6 ， theta6 ） ； k7 = PlaneTrussElementStiffness（ E ， A， L7 ， theta7 ） ； k8 = PlaneTrussElementStiffness（ E ， A， L8 ， theta8 ） ； k9 = PlaneTrussElementStiffness（ E ， A， L9 ， theta9 ） ； k10 = PlaneTrussElementStiffness （ E ， A， L10 ，
1
桁架结构有限元分析的基本原理
桁架结构在屋架、 桥梁等工程中应用很多． 在
桁架结构中， 各杆件的连接处近似为铰接， 杆件绕 结点可以自由转动， 结构中的各个杆件只受轴力，
收稿日期：2012 － 05 － 26 基金项目： 广 东 省 自 然 科 学 基 金 项 目 （ S2011010001679 ） ； 2011 年梅州市科技计划项目 作者简介： 李远瑛 （ 1973 － ） ， 女， 广东平远人， 副教授， 硕士， 主要研究方向： 结构理论研究与计算．
模态 1 2 3 4 5
Freque ncy（ hz） MATLAB 35 ． 108 76 ． 832 95 ． 892 184 ． 67 226 ． 57 ANSYS 35 ． 108 76 ． 832 95 ． 892 184 ． 67 226 ． 57
模态 6 7 8 9
2 ． 6 ANSYS 计算结果 为验证 MATLAB 编制该桁架有限元分析程序 使用 ANSYS 有限元分析软件分别对该结 的正确性， 构进行静力及模态分析计算． 在建立有限元分析模 使用 LINK1 单元， 结果如表 3 及图 5 、 图 6、 图 型中， 7 所示， 二者计算结果完全一致．
轴力 应力 单元 5 8 / 10 N / 10 N·m － 2 2 ． 666 7 1 ． 066 7 6 7 8 9 10
－ 1． 2 0． 7 1 ． 066 7
－ 3 ． 333 3 － 1 ． 333 3 1 ． 75 2 ． 333 3 0． 7 0 ． 933 33
图1
平面桁架杆件单元整体坐标与局部坐标
30 1． 2 坐标转换矩阵
嘉应学院学报（ 自然科学）
2012 年 8 月
其中： F e 为局部坐标下单元的轴力； k e 为局部坐标下单元的刚度矩阵； T 为坐标转换矩阵； △ 为整体坐标系下单元的结点位移向量， 由 单元定位向量确定． σ 为局部坐标下单元的应力．
展， 有限元分析得到了广泛应用． MATLAB 是效率 可以实现编程计算、 数据 高的工程计算软件工具， 图形化集成． 本文利用 MATLAB 软件， 利用有限元 分析的基本原理对平面桁架进行受力分析， 得到了 结构各杆件的内力、 应力、 结构变形及结构动力特 性， 并将计算结果同 ANSYS 有限元分析软件结果进 行比较． 结果表明， 利用 MATLAB 可以编制结构有 限元分析程序， 计算速度快， 数据准确， 并能为复杂 结构有限元分析提供参考．
单元 图2 平面桁架计算简图
2 ． 1 单元编号， 结点编号 各单元编号， 单元始末端结点编号， 见表 1 所示．
表1 各单元始末端结点编号 单元 （ 6） （ 7） （ 8） （ 9） （ 10 ） 始末端结 点编号 3 →5 4 →5 4 →6 5 →6 2 →5
1． 6 建立结构的刚度方程
（ 1）
［4 ］
．
32 结点荷载向量 U1 = K1 \ P ；
嘉应学院学报（ 自然科学）
2012 年 8 月
U =［ 0 ； 0 ； U1 （ 1 ： 9 ） ； 0 ］ ； % 生成整体坐标系下结点 位移向量 …… 2 ． 4 编制 MATLAB 平面桁架动力特性分析程序 结构动力特性分析方法有矢量迭代法、 子空间 迭代法． MATLAB 提供了计算矩阵特征值的专门函 数， 即 eig 和 eigs， 它们既可以用于求解标准特征值 也可以求解广义特征值问题． 其调用格式为： 问题， E = eig（ K ， M ） 等． 本文为研究桁架的动力特性， 编制了桁架的动 力特性 计 算 程 序， 得 到 了 结 构 前 9 阶 的 频 率， 见 表 3 ． 由于本文篇幅有限， 动力特性计算程序略． 2 ． 5 MATLAB 计算结果 表 2 为 MATLAB 编制程序的静力计算结果， 包 括桁架各杆件的轴力和应力． 为直观显示桁架受力 利用 MATLAB 中 line 命令编制了结构 之后的变形， 变形前后的示意图． 图 4 为桁架各结点位移变形示 意图（ 变形放大比例为 30 倍） ， 实线为变形前， 虚线 为变形后． 表 3 为 MATLAB 编制程序的动力特性计 算结果．
图 5 桁架各杆件轴力示意图 表2 单元 1 2 3 4 5 静力计算结果 轴力 / 10 5 N －3 1 ． 75 2 ． 666 7 应力 / 10 N·m － 2
8
图 4 桁架各结点变形示意图 表3 动力特性计算结果 Frequency（ hz） MATLAB 234 ． 66 242 ． 5 272 ． 92 313 ． 26 ANSYS 234 ． 66 242 ． 5 272 ． 92 313 ． 26
0
前言
有限元法是根据变分原理求解数学及物理问
属于一维问题， 各杆件内的位移线性变化， 杆件内 的位移分布函数是精确的
［2 ］
． 在进行结构离 散 化
题的数值计算方法
［1 ］
． 随着计算机技术的快速 发
时， 不需要对杆件再细分单元， 直接以杆件为分析 单元． 1． 1 局部坐标系下的杆件单元刚度矩阵 桁架中任意杆件 ij 为分析单元， 建立结构的整 体坐标系为 xoy， 局部坐标系为 xoy， 如图 1 所示． 一 维杆件单元在局部坐标中的单元刚度矩阵为 1 0 EA Ke = l －1 0 0 0 0 0 －1 0 1 0 0 0 0 0
第 30 卷
第8 期
嘉应学院学报（ 自然科学） JOURNAL OF JIAYING UNIVERSITY（ Natural Science）
Vol． 30
No． 08
2012 年 8 月
Aug． 2012
基于 MATLAB 平面桁架有限元分析研究
李远瑛， 张德生
（ 嘉应学院 土木工程学院 ， 广东 梅州 514015 ） 摘 要：基于平面桁架有限元分析的基本原理 ， 利用 MATLAB 软件对某平面桁架编制了计算程序 ， 得到了结构的