2022年平面三角形单元有限元程序设计.pdf

《2022年平面三角形单元有限元程序设计.pdf》由会员分享，可在线阅读，更多相关《2022年平面三角形单元有限元程序设计.pdf（13页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、平面三角形单元有限元程序设计P9 m9m一 、 题 目如 图 1 所 示 ,一 个 厚 度 均 匀 的 三 角 形 薄 板 ,在 顶 点 作 用 沿 板 厚 方 向 均 匀分 布 的 竖 向 载 荷 。已 知 :P=150N/m,E=200GPa,=0 、25,t=0 、1m, 忽 略自 重 。 试 计 算 薄 板 的 位 移 及 应 力 分 布 。要 求 : 1. 编 写 有 限 元计 算 机 程 序 ,计 算 节点 位 移 及 单元 应 力 。(划 分 三 角 形单 元 ,单 元 数不 得 少 于 30 个 ); 2. 采 用 有 限 元软 件 分 析 该 问 题 (有限 元 软 件 网格

2、 与 程 序 设 计 网 格 必须 一 致 ),详 细给 出 有 限 元 软 件 每 一步 的 操 作 过程 ,并 将 结 果 与 程 序计 算 结 果 进 行 对 比 (任 选 取 三 个 点 ,对 比 位 移 值 ); 3. 提 交 程 序 编 写 过 程 的 详 细 报 告 及 计 算 机 程 序 ; 4. 所 有 同 学 参 加 答 辩 ,并 演 示 有 限 元 计 算 程 序 。有 限 元 法 中 三 节 点 三 角 形 分 析 结 构 的 步 骤 如 下 :1)整 理 原 始 数 据 ,如 材 料 性 质 、 荷 载 条 件、 约 束条 件 等 ,离 散 结 构 并 进行 单 元

3、编 码 、 结 点 编 码 、 结 点 位 移 编 码 、 选 取 坐 标 系 。2)单 元 分 析 ,建 立 单 元 刚 度 矩 阵 。3)整 体 分 析 ,建 立 总 刚 矩 阵 。4)建 立 整 体 结 构 的 等 效 节 点 荷 载 与 总 荷 载 矩 阵5)边 界 条 件 处 理 。6)解 方 程 ,求 出 节 点 位 移 。7)求 出 各 单 元 的 单 元 应 力 。8)计 算 结 果 整 理 。精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 1 页，共 13 页 - - - - - - -

4、 - - - 平面三角形单元有限元程序设计一、 程序 设 计网 格 划 分如 图 ,将 薄 板 如 图 划 分 为 6 行 ,并 建 立 坐 标 系 ,则刚 度 矩 阵 的 集 成建 立 与 总 刚 度 矩 阵 等 维 数 的 空 矩 阵 ,已 变 单 元 刚 度 矩 阵 的 集 成 。由 单 元 分 析 已 知 节 点 、 单 元 的 排 布 规 律 ,继 而 通 过 循 环 计 算 求 得 每 个单 元 对 应 的 节 点 序 号 。通 过 循 环 逐 个 计 算 :(1) 每 个 单 元 对 应 2 种 单 元 刚 度 矩 阵 中 的 哪 一 种 ; (2) 该 单 元 对 应 总 刚

5、度 矩 阵 的 那 几 行 哪 几 列(3) 将 该 单 元 的 单 元 刚 度 矩 阵 加 入 总 刚 度 矩 阵的 对 应 行 列循 环 又 分 为 3 层 循 环 :(1) 最 外 层 :逐 行 计 算(2) 中 间 层 :该 行 逐 个 计 算(3) 最 里 层 :区 分 为 第奇 / 偶数 个 计 算XYPXYP节点编号单元编号精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 2 页，共 13 页 - - - - - - - - - - 平面三角形单元有限元程序设计单 元 刚 度 的 集 成 :2

6、15656665656266256561661eZeeeZeZeeeekkkKkkkkkk边 界 约 束 的 处 理 :划 0 置 1 法适 用 :这 种 方 法 适 用 于 边 界 节 点 位 移 分 量 为 已 知 (含 为 0) 的 各 种 约束 。做 法 :(1)将 总 刚 矩 阵 K 中 相 应 于 已 知 位 移 行 主 对 角 线 元 素 置 1,其 她 元素 改 为 零 ;同时 将 载 荷 列 阵 R 中 相 应 元 素 用 已 知 位 移 置 换 。这 样 ,由 该 方 程 求 得 的 此 位 移 值 一 定 等 于 已 知 量 。(2)将 K 中 已 知 位 移 相 应 的

7、 列 的 非 主 对 角 成 元 素 也 置 0,以 保 持 K的 对 称 性 。当 然 ,在 已 知 位 移 分 量 不 为 零 的 情 况 下 ,这 样 做 就 改 变 了 方 程 左 端的 数 值 ,为保 证 方 程 成 立 ,须 在 方 程 右 端 减 去 已 知 位 移 对 该 方 程 的 贡 献 已 知位 移 与 相 应 总 刚 元 素 的 乘 积 。若 约 束 为 零 位 移 约 束 时 ,此步 则 可省 去 。特 点 :精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 3 页，共 13 页

8、- - - - - - - - - - 平面三角形单元有限元程序设计(1)经 以 上 处 理 同 样 可 以 消 除 刚 性 位 移 (约 束 足 够 的 前 提 下 ),去 掉 未知 约 束 反 力 。(2)但 这 种 方 法 不 改 变 方 程 阶 数 ,利 于 存 贮 。(3)不 过 ,若 就 是 要 求 出 约 束 反 力 ,仍 要 重新 计 算 各 个 划 去 的 总 刚 元素 。程 序 如 下 : 变 量 说 明NNODE 单 元 节 点 数NPION 总 结 点 数NELEM 单 元 数NVFIX 受 约 束边 界 点 数FIXED 约 束 信 息 数 组NFORCE 节 点 力

9、 数FORCE 节 点 力 数 组COORD 结 构 节 点 坐 标 数 组LNODS 单 元 定 义 数 组YOUNG 弹 性 模 量POISS 泊 松 比THICK 厚 度B 单 元 应 变 矩 阵 (3*6) D 单 元 弹 性 矩 阵 (3*3) S 单 元 应 力 矩 阵 (3*6) 精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 4 页，共 13 页 - - - - - - - - - - 平面三角形单元有限元程序设计A 单元 面 积ESTIF 单 元 刚 度 矩阵ASTIF 总 体 刚 度

10、矩 阵ASLOD 总 体 荷 载 向 量ASDISP 节 点 位 移 向 量ELEDISP 单 元 节 点 位移 向 量STRESS 单 元 应 力%*% 初 始 化clear format short e % 设 定 输 出 类 型clear % 清 除 内 存 变 量NELEM=36 % 单 元 个 数 (单 元 编 码 总 数 ) NPION=28 % 结 点 个 数 (结 点 编 码 总 数 ) NVFIX=2 % 受 约束 边 界点 数NFORCE=1 % 结 点 荷 载 个 数YOUNG=2e11 % 弹 性 模 量POISS=0 、 25 % 泊 松 比THICK=0 、 1 %

11、 厚 度LNODS=1 2 3;2 4 5;2 5 3;3 5 6;4 7 8;4 8 5;5 8 9;5 9 6; 6 9 10;7 11 12;7 12 8;8 12 13;精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 5 页，共 13 页 - - - - - - - - - - 平面三角形单元有限元程序设计8 13 9;9 13 14;9 14 10;10 14 15; 11 16 17;11 17 12; 12 17 18; 12 18 13;13 18 19; 13 19 14;14 19 2

12、0;14 20 15;15 20 21;16 22 23;16 23 17;17 23 24;17 24 18;18 24 25;18 25 19;25 19 26;19 26 20;20 26 27;20 27 21;21 27 28 % 单 元 定 义 数 组(单 元 结 点 号 ) % 相 应 为 单 元 结 点 号 (编 码 )、 按逆 时 针 顺 序 输 入COORD=0 0;-0、 75 1 、 5;0 、 75 1 、 5;-1 、 5 3;0 3; 1、 5 3;-2 、 25 4 、 5;-0 、 75 4 、 5;0 、 75 4 、 5; 2、 25 4 、 5;-3 6

13、;-1、 5 6;0 6;1、 5 6;3 6; -3 、 75 7 、 5;-2 、 25 7 、 5; -0 、 75 7 、 5;0 、 75 7 、 5; 2、 25 7 、 5;3、 75 7 、 5;-4 、 5 9;-3 9; -1 、 5 9;0 9;1、 5 9;3 9;4 、 5 9 % 结 点 坐 标 数 组% 坐 标 :x,y 坐 标 (共NPOIN 组 ) FORCE=1 0 -15 % 结点 力 数组 (受 力 结 点 编 号 , x 方 向 ,y 方 向 ) FIXED=22 1 1;28 1 1 % 约 束 信 息 (约 束 点 ,x 约 束 ,y 约 束) %

14、有 约 束 为1, 无 约束 为0 %*% 生 成 单 元 刚 度 矩 阵 并 组 成 总 体 刚 度 矩 阵ASTIF=zeros(2*NPION,2*NPION); % 生 成 特 定 大 小 总 体 刚 度精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 6 页，共 13 页 - - - - - - - - - - 平面三角形单元有限元程序设计矩 阵 并 置0 %* for i=1:NELEM % 生 成 弹 性 矩 阵D D= 1 POISS 0; POISS 1 0; 0 0 (1-POISS)/

15、2*YOUNG/(1-POISS2) %*% 计 算 当 前 单 元 的 面 积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,(r

16、em(j+1),3)+1),1)+COORD(LNODS(i,(rem(j+2),3)+1),1); 精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 7 页，共 13 页 - - - - - - - - - - 平面三角形单元有限元程序设计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

17、*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 精品资料 - - - 欢迎下载 - - - - - - -

18、 - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 8 页，共 13 页 - - - - - - - - - - 平面三角形单元有限元程序设计%*% 将 约 束 信 息 加 入 总 体 刚 度 矩 阵 (对 角 元 素 改 一 法 ) 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 %*% 生 成

19、 单 元 刚 度 矩 阵 并 组 成 总 体 刚 度 矩 阵%* 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 精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - -

20、 - - - -第 9 页，共 13 页 - - - - - - - - - - 平面三角形单元有限元程序设计ASLOD(FORCE(i,1)*2-1):FORCE(i,1)*2)=FORCE(i,2:3); end %* % 求 解 内 力ASDISP=ASTIFASLOD % 计 算 节 点 位 移 向 量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); % 取 出 当 前 单 元 的 节 点 位 移 向 量

21、end i STRESS=D*B1(:, :, i)*ELEDISP % 求 内 力end (程 序 计 算 结 果 与 有 限 元 软 件 得 出 的 结 果 稍 有 偏 差 ,可 能 就 是 程 序 某 些地 方 数 据 输 入 时 出 了 问 题 ,还 在 寻 找 具 体 原 因 ) 二、 有 限元 软 件分 析设置材料参数精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 10 页，共 13 页 - - - - - - - - - - 平面三角形单元有限元程序设计精品资料 - - - 欢迎下载 - -

22、 - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 11 页，共 13 页 - - - - - - - - - - 平面三角形单元有限元程序设计网格划分边界约束添加载荷精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 12 页，共 13 页 - - - - - - - - - - 平面三角形单元有限元程序设计精品资料 - - - 欢迎下载 - - - - - - - - - - - 欢迎下载 名师归纳 - - - - - - - - - -第 13 页，共 13 页 - - - - - - - - - -