有限元的matlab编程.ppt

有限元编程示例 题目描述 如下图所示的平面桁架 杆件长度 弹性模量 截面积以及所受节点力P的大小可以自行定义 求节点位移及杆件轴力 例一 桁架 解题思路 建立模型集成总刚求解位移求解杆件轴力输出结果 建立模型 定义节点坐标 Node zeros 10 2 x 1 L L为横杆长度fori 1 2 10 x x L Node i x0 endx 1 L fori 2 2 10 x x L Node i xH H为竖杆长度end 模型相关参数输入 H input 竖杆长度 m L input 水平杆长度 m E input 杆件弹性模量 Gpa A input 杆件截面积 m 2 a input 节点力P kN 节点编号方式 定义单元 即储存单元两端的节点号 Element zeros 21 2 fori 1 2 7Element 5 2 i 3 2 i i 1 Element 5 2 i 1 2 i i 2 Element 5 2 i 1 2 i i 3 endfori 2 2 8Element 5 i 2 1 i i 1 Element 5 i 2 i i 2 endElement 21 9 10 加下划线的为单元编号 集成总刚 xi Node Element ie 1 1 ie为单元号 以下相同yi Node Element ie 1 2 xj Node Element ie 2 1 yj Node Element ie 2 2 获取单元两端节点坐标 L xj xi 2 yj yi 2 1 2 计算杆件长度 形成等效荷载列阵 f 0 0 0 a 0 0 0 a 0 0 0 a 0 0 0 a 0 0 0 a 每个节点两个自由度 a为之前输入的节点力 计算从局部坐标到整体坐标的坐标转换矩阵T functionT TransformMatrix ie ie为单元号 c xj xi L s yj yi L T c s00sc0000c s00sc 计算单元刚度矩阵k k E A L0 E A L00000 E A L0E A L00000 T TransformMatrix ie k T k transpose T transpose T 为T的转置矩阵2 集成整体刚度矩阵K forie 1 1 21 按单元顺序进行循环k PlaneTrussElementStiffness ie 计算第ie个单元的单刚m Element ie 1 ie单元的首节点号n Element ie 2 ie单元的末节点号K 2 m 1 2 n 1 k 1 3 K 2 m 1 2 n k 1 4 K 2 m 2 n 1 k 2 3 K 2 m 2 n k 2 4 K zeros 20 20 用来存储整体刚度矩阵 集成总刚的非对角线元素 这里的元素指2 2的小矩阵 在下面的集成中 将总刚看成10 10的矩阵 每个元素为2 2的小矩阵 m Element ie 2 ie单元的末节点号n Element ie 1 ie单元的首节点号K 2 m 1 2 n 1 k 3 1 K 2 m 1 2 n k 3 2 K 2 m 2 n 1 k 4 1 K 2 m 2 n k 4 2 end 集成总刚的对角线元素 这里的元素指2 2的小矩阵 fori 1 1 10 按节点的顺序循环forj 1 1 21 对于每个节点 再按单元的顺序循环k PlaneTrussElementStiffness j ifElement j 1 I 如果i节点为j单元的首节点K 2 i 1 2 i 1 K 2 i 1 2 i 1 k 1 1 K 2 i 1 2 i K 2 i 1 2 i k 1 2 K 2 i 2 i 1 K 2 i 2 i 1 k 2 1 K 2 i 2 i K 2 i 2 i k 2 2 endifElement j 2 i 如果i节点为j单元的末节点K 2 i 1 2 i 1 K 2 i 1 2 i 1 k 3 3 K 2 i 1 2 i K 2 i 1 2 i k 3 4 K 2 i 2 i 1 K 2 i 2 i 1 k 4 3 K 2 i 2 i K 2 i 2 i k 4 4 endendend 求解位移 u zeros 20 根据约束情况修改总刚 采用对角元素置1法 fori 1 1 20K 1 i 0 K 2 i 0 K 18 i 0 K i 1 0 K i 2 0 K i 18 0 end 自由度1 2 18被约束了 所在的行和列的其他元素都改为0 K K 1e15 乘以一个大数 减小计算误差f f 1e15 u K f 求解 K 1 1 1 对角线元素置1K 2 2 1 K 18 18 1 求解轴力 获取单元两端的节点号 i Element ie 1 ie为单元号j Element ie 2 获取单元两端的节点位移 uElement zeros 4 1 uElement 1 2 u i 1 2 1 i 1 2 2 uElement 3 4 u j 1 2 1 j 1 2 2 k PlaneTrussElementStiffness ie nodef k uElement 整体坐标下的节点力T TransformMatrix ie 计算坐标转换矩阵nodef transpose T nodef 从整体坐标转换到局部坐标 计算单元的节点力 输出求解结果 输出位移 fprintf 节点位移 n fori 1 1 10disp 节点号 num2str i x方向位移 num2str u 2 i 1 1 y方向位移 num2str u 2 i 1 end 输出节点力 fprintf n n节点力 n forie 1 1 21nodef NodeForce ie disp 单元号 num2str ie 节点号 num2str Element ie 1 节点号 num2str Element ie 2 轴力 num2str nodef 1 end 例二 网架 思路分析 网架是由多根杆件按照一定的网格形式通过节点连结而成的空间结构 构成网架的基本单元有三角锥 三棱体 正方体 截头四角锥等 鉴于网架的形式较多 本程序提供一种通用的网架输入方法 但录入较为繁琐 同时提供一种正放四角锥网架的简易输入方法作为典型 考虑几何非线性 本程序采用荷载分级加载的方式考虑网架的几何非线性 将总荷载分成1000份分步施加 求解各荷载步下的节点位移 修改网架相应节点坐标以及刚度矩阵 依次迭代求出网架的总位移 本程序的网架位移求解函数附在主程序后面 主程序运行时调用该函数 几点说明 e input 选择网架类型 0代表自由定义网架 1代表四角锥网架 网架类型的选择 网架类型的选择 用户自定义网架 网架信息的录入 包括节点 单元 截面 弹性模量等 ife 0 选择自定义网架 Node input 定义节点编号及对应坐标 按 1x1y1z1 2x2y2z2 输入 形成节点储存矩阵 Men input 定义单元与节点的关系 按 1node1node2 2node3node4 输入 node1 node2 依次类推 形成单元储存矩阵 Msum length Men 查找网架录入的单元数 Cont1 input 定义单元实常数 若所有杆件截面面积和弹性模量不变 则输入0 否则输入1 定义单元属性的输入方式 ifCont1 0 AE1 input 请输入统一的截面面积与弹性模量 按 AE 输入 AE zeros Msum 3 AE 1 1 Msum AE 2 AE1 1 1 AE 3 AE1 1 2 else AE input 请输入相应单元的截面面积与弹性模量 按 1 A1E1 2 A2E2 输入 end P input 定义节点荷载 按 node1P1 node2P2 输入 网架荷载输入 BC input 定义边界约束 按 node1ConxConyConz node2ConxConyConz 输入 Con代表x y z方向约束 取0为约束 取1无约束 网架边界条件 end 单元属性相同 单元属性不同 荷载及边界条件 正放四角锥网架定义 ife 1 hu input 输入网架上层节点行数 定义网架上层节点的行数lu input 输入网架上层节点列数 定义网架上层节点的列数 dis xu input 输入网架上层节点列间距 定义网架上层的行间距dis yu input 输入网架上层节点行间距 定义网架上层的列间距 hd hu 1 网架下层节点的行数ld lu 1 网架下层节点的列数 dis xd dis xu 网架下层的行间距dis yd dis yu 网架下层的行间距 dis z input 输入网架上下层间距 网架上下层间距 定义网架上层节点 定义网架下层节点 定义网架高度 fori 1 huforj 1 luNode i 1 lu j 2 j 1 dis xu Node i 1 lu j 3 i 1 dis yu Node i 1 lu j 4 dis z endend fori 1 hdforj 1 ldNode i 1 ld j hu lu 2 j 1 0 5 dis xd Node i 1 ld j hu lu 3 i 1 0 5 dis yd Node i 1 ld j hu lu 4 0 endend 网架上层节点编号与对应坐标 网架下层节点编号与对应坐标 Nsum length Node 查询网架的节点数 fori 1 Nsum 将节点编号录入节点矩阵Node i 1 i end fori 1 huforj 1 lu 1Men i 1 lu 1 j 2 i 1 lu j Men i 1 lu 1 j 3 i 1 lu j 1 endend fori 1 luforj 1 hu 1Men i 1 hu 1 j lu 1 hu 2 j 1 lu i Men i 1 hu 1 j lu 1 hu 3 j lu i endend 节点编号的录入 网架上层横向单元的拓扑 网架上层纵向单元的拓扑 fori 1 hdforj 1 ld 1Men i 1 ld 1 lu 1 hu hu 1 lu j 2 i 1 ld j hu lu Men i 1 ld 1 lu 1 hu hu 1 lu j 3 i 1 ld j hu lu 1 endend fori 1 ldforj 1 hd 1Men i 1 hd 1 ld 1 hd lu 1 hu hu 1 lu j 2 j 1 ld i hu lu Men i 1 hd 1 ld 1 hd lu 1 hu hu 1 lu j 3 j ld i hu lu endend 网架下层纵向单元的拓扑 网架下层横向单元的拓扑 网架腹杆单元的拓扑 fori 1 hdforj 1 ld Men i 1 ld j 1 4 hu 1 lu lu 1 hu hd 1 ld ld 1 hd 1 2 i 1 lu j Men i 1 ld j 1 4 hu 1 lu lu 1 hu hd 1 ld ld 1 hd 1 3 i 1 ld hu lu j Men i 1 ld j 1 4 hu 1 lu lu 1 hu hd 1 ld ld 1 hd 2 2 i 1 lu j 1 Men i 1 ld j 1 4 hu 1 lu lu 1 hu hd 1 ld ld 1 hd 2 3 i 1 ld j hu lu Men i 1 ld j 1 4 hu 1 lu lu 1 hu hd 1 ld ld 1 hd 3 2 i lu j Men i 1 ld j 1 4 hu 1 lu lu 1 hu hd 1 ld ld 1 hd 3 3 i 1 ld j hu lu Men i 1 ld j 1 4 hu 1 lu lu 1 hu hd 1 ld ld 1 hd 4 2 i lu j 1 Men i 1 ld j 1 4 hu 1 lu lu 1 hu hd 1 ld ld 1 hd 4 3 i 1 ld j hu lu endend 腹杆N 腹杆N 1 腹杆N 2 腹杆N 3 单元编号录入单元储存矩阵 Msum length Men 查询网架单元数 fori 1 Msum 将单元编号录入单元储存矩阵Men i 1 i end 定义截面属性 E 2 1e11 默认材料为steel A1 input 请输入网架上层单元的截面面积 默认网架上层单元截面尺寸相同A2 input 请输入网架下层单元的截面面积 默认网架下层单元截面尺寸相同A3 input 请输入网架腹杆单元的截面面积 默认网架腹杆单元截面尺寸相同 AE zeros Msum 3 定义单元属性矩阵 m1 hu 1 lu lu 1 hu 上层单元截止编号m2 m1 hd 1 ld ld 1 hd 下层单元截止编号 AE 1 m1 2 A1 将上层单元尺寸录入AE矩阵AE m1 1 m2 2 A2 将下层单元尺寸录入AE矩阵AE m2 1 Msum 2 A3 将腹杆单元尺寸录入AE矩阵 AE 1 1 Msum 将单元编号录入AE矩阵 AE 3 E 将材料弹性模量录入AE矩阵 定义荷载 cont2 input 定义节点荷载 若网架上层节点力与下层节点力均布 则输入0 否则输入1 ifcont2 0P1 input 请输入网架上层节点荷载 P2 input 请输入网架下层节点荷载 m3 hu lu P 1 Nsum 1 1 Nsum P 1 m3 2 P1 P m3 1 Nsum 2 P2 elseP input 定义节点荷载 按 node1P1 node2P2 输入 end 定义边界条件 cont3 input 定义边界约束 若网架上层周边节点全约束 则输入0 若下层周边节点全约束 输入1 否则输入2 ifcont3 0n1 2 hu lu 2 BC zeros n1 4 BC 1 lu 2 1 2 lu 1 BC lu 1 2 lu 4 1 lu hu 1 2 lu hu 1 BC 2 lu 3 2 lu 4 hu 1 1 lu lu hu 1 1 BC 2 lu 3 hu n1 1 lu lu hu lu elseifcont3 1n1 2 hd ld 2 BC zeros n1 4 BC 1 ld 2 1 2 ld 1 BC ld 1 2 ld 4 1 ld hd 1 2 ld hd 1 BC 2 ld 3 2 ld 4 hd 1 1 ld ld hd 1 1 BC 2 ld 3 hd n1 1 ld ld hd ld fori 1 n1BC i 1 BC i 1 hu lu end elseBC input 定义边界约束 按 node1ConxConyConz node2ConxConyConz 输入 Con代表x y z方向约束 取0为约束 取1无约束 endend Nsum length Node Msum length Men Psum length P BCsum length BC 提取各矩阵的行数 考虑几何非线性分析网架 fori 1 Psum 将力分为1000份P i 2 P i 2 1000 end U zeros 3 Nsum 1 总位移矩阵 fori 1 1000 u L1 Kz grid Node Men AE P BC Nsum Msum Psum BCsum forj 1 Nsum 根据节点位移修改网架的节点坐标Node j 2 Node j 2 u 3 j 2 1 Node j 3 Node j 3 u 3 j 1 1 Node j 4 Node j 4 u 3 j 1 endU U u 每次迭代位移的叠加end 迭代法修正刚度矩阵和网架位移 求解网架杆件的应力应变 L0 zeros Msum 1 所有根杆的最初长度 fori 1 Msum 单元两端的节点编号p Men i 2 q Men i 3 X1 Node p 2 单元端节点的坐标Y1 Node p 3 Z1 Node p 4 X2 Node q 2 Y2 Node q 3 Z2 Node q 4 L0 i 1 sqrt X2 X1 2 Y2 Y1 2 Z2 Z1 2 网架杆件的初始长度end Lt L1 L0 所有杆件长度的增量 strain zeros Msum 1 定义应变矩阵stress zeros Msum 1 定义应力矩阵 fori 1 MsumE AE i 3 strain i 1 Lt i 1 L0 i 1 第i根杆件应变stress i 1 E strain i 1 第i根杆件应力end disp U 输出网架节点位移disp stress 输出网架杆件应力 网架节点位移求解函数 function u L1 Kz grid Node Men AE P BC Nsum Msum Psum BCsum Kz zeros 3 Nsum 3 Nsum 定义刚度矩阵L zeros Msum 1 fori 1 Msum 单元两端的节点编号p Men i 2 q Men i 3 A AE i 2 E AE i 3 单元两端节点坐标X1 Node p 2 Y1 Node p 3 Z1 Node p 4 X2 Node q 2 Y2 Node q 3 Z2 Node q 4 单元长度L i 1 sqrt X2 X1 2 Y2 Y1 2 Z2 Z1 2 单元的方向余弦l X2 X1 L i 1 m Y2 Y1 L i 1 n Z2 Z1 L i 1 定义转化矩阵t sqrt l 2 n 2 ift 0r 0m0 m00 001 elser lmn l m tt m n t n t0l t endO zeros 3 3 T rO Or 整体坐标下的单刚矩阵ke E A L i 1 100 100 000000 00