matlab编译平面有限元计算(附有完整代码)

完整代码下载链接
点击此处下载哦 下载后运行‘main.m’即可

问题描述
使用完成的代码,解决图1所示的平面应力问题。中心孔半径为A的均匀薄板承受单轴应力。a=0.5 in.,h=3 in.,w=6 in.,E=10(10)6 psi,泊松比=0.3。计算应力分布。

1. 模型绘制与网格划分
此问题第一步需要解决的就是网格的划分,模型绘制可以使用Matlab自带的矩形和圆工具,对于网格的划分,以一种较为简单的方式进行,这样也可以减低一些编程难度,效果如图所示,我全部使用三角形单元,对于中间的圆形孔洞,圆形的网格划分在程序上存在一定难度,我用方形空洞代替,在划分网格的同时,对于节点的编号和坐标将其记录在一个Node_information矩阵中,对于三角形单元的编号和包含的节点保存在Element_information中。(此部分对应的程序为Meshing.m)


2. 单元刚度矩阵的建立与整体刚度矩阵装配
得到单元信息和节点信息之后,下面的操作类似于桁架的有限元分析,需要建立每一个单元对应的刚度矩阵,对于三角形单元来说,课件中有详细的说明,其核心利用到的公式如下:
下面我们将公式转化为代码语言,D矩阵比较容易书写,B矩阵则要利用上步储存的Node_information和Element_information,通过矩阵的索引准确的找到每个三角形单元对应的节点坐标,再通过矩阵乘法则可以得到每个单元的刚度矩阵(66)。
为了利于后续的装配,我们将每个单元的刚度矩阵赋值到一个‘小装配矩阵’,此矩阵含义为把一个(6
6)矩阵按照公式1规则装配到一个(2*节点总数)的方阵中,此规则不同于桁架结构,因为其刚度矩阵是按照(u1 u2 u3 v1 v2 v3)的顺序来生成的,对应的匹配规则也要得到调整。(完成此部的程序文件对应为Single_triangular.m)

function [K,B,D]=Single_triangular(Num_node,Node_sum,E,Poisson)
%Node_sum 是3*3表示三角形三个节点坐标 前一个是编号 后两个是坐标
K=zeros(2*Num_node);
A=0.5*abs(det([1 Node_sum(1,2) Node_sum(1,3);...1 Node_sum(2,2) Node_sum(2,3);1 Node_sum(3,2) Node_sum(3,3)]));
B=(1/(2*A))*[Node_sum(2,3)-Node_sum(3,3) Node_sum(3,3)-Node_sum(1,3) Node_sum(1,3)-Node_sum(2,3) 0 0 0;...0 0 0 Node_sum(3,2)-Node_sum(2,2) Node_sum(1,2)-Node_sum(3,2) Node_sum(2,2)-Node_sum(1,2);...Node_sum(2,3)-Node_sum(3,3) Node_sum(3,3)-Node_sum(1,3) Node_sum(1,3)-Node_sum(2,3)...Node_sum(3,2)-Node_sum(2,2) Node_sum(1,2)-Node_sum(3,2) Node_sum(2,2)-Node_sum(1,2)];
D=(E/(1-Poisson^2))*[1 Poisson 0;Poisson 1 0;0 0 (1-Poisson)/2];
k=A*B'*D*B;
list=[Node_sum(1,1) Node_sum(2,1) Node_sum(3,1) ...Num_node+Node_sum(1,1) Num_node+Node_sum(2,1) Num_node+Node_sum(3,1)];
for i =1:6for j =1:6K(list(i),list(j))=k(i,j);end
end
end


由于得到了‘小装配矩阵’,接下来对每个三角形单元得到的‘小装配矩阵’作累加即可,这里只需要对Element_information遍历即可。完成此部的程序文件对应为Assemble.m)

function [KK,B,D]=Assemble(Node_information,Element_information,E,Poisson)
Num_node=size(Node_information,1);
Num_element=size(Element_information,1);
KK=zeros(2*Num_node);
B=[];
for i =1:Num_elementNode_sum=[Element_information(i,2) ...Node_information(Element_information(i,2),2) ...Node_information(Element_information(i,2),3);...Element_information(i,3) ...Node_information(Element_information(i,3),2) ...Node_information(Element_information(i,3),3);...Element_information(i,4) ...Node_information(Element_information(i,4),2) ...Node_information(Element_information(i,4),3)];[K,b,D]=Single_triangular(Num_node,Node_sum,E,Poisson);B=[B;b];KK=KK+K;
end
end

3. 外载转化至节点
在本题目中,板子两侧受到均布外载,我们需要将其转化到节点上,主要利用对形函数的积分上,其核心如公式2和所示,对于本题而言是均布荷载,简单的笔算可以得到节点的等效力,编程这里就直接给出。(对应程序文件为Load_plane.m)


4. 求解节点位移
经过上几步的准备工作,已经得到了刚度K矩阵和外载F矩阵,只需要对K求逆再与F相乘就可以得到结果了。但此时存在一个问题,对于刚度矩阵数量较为庞大,很容易出现矩阵接近奇异值的情况,Matlab处理此问题较为强大,提供了四种求逆的方法。如下表所示:

表1 求逆方法一览表

方法 介绍
Inv(K)*F 直接求逆,对于病态矩阵产生很大错误
Pinv(K)*F 伪逆,较为准确
K\F 高斯消主元,结果误差较大
Lsqminnorm(K,F) 最小范数,最为准确
————————————————————————

通过最小范数求得的结果较为准确,如表2展示的是四种方法得到的结果(前301个节点x方向位移)。由于计算出的位移非常小,下图是我将位移扩大1000倍展示的位移图,如图3,红色虚线为变形后图形。为了对比保证结果准确性,我们使用Ansys建模分析,得到图4,两者大致一样。(对应程序文件为Output_plane_displacement.m)

表2 不同算法结果

Lsqminnorm(K,F) Inv(K)*F Pinv(K)*F K\F
-6.1616E-04 -4.9019E-04 -6.1616E-04 -5.7071E-04
-6.1704E-04 -4.8256E-04 -6.1704E-04 -5.7158E-04
-6.1815E-04 -4.7112E-04 -6.1815E-04 -5.7268E-04
-6.1943E-04 -4.7684E-04 -6.1943E-04 -5.7395E-04
-6.2066E-04 -4.6921E-04 -6.2066E-04 -5.7517E-04
-6.2161E-04 -4.6539E-04 -6.2161E-04 -5.7611E-04
-6.2209E-04 -4.5013E-04 -6.2209E-04 -5.7657E-04
-6.2204E-04 -4.5013E-04 -6.2204E-04 -5.7651E-04
-6.2154E-04 -4.5013E-04 -6.2154E-04 -5.7601E-04
-6.2078E-04 -4.2343E-04 -6.2078E-04 -5.7524E-04
-6.1998E-04 -4.3869E-04 -6.1998E-04 -5.7443E-04
-6.1936E-04 -4.2343E-04 -6.1936E-04 -5.7380E-04
-6.1902E-04 -4.1962E-04 -6.1902E-04 -5.7345E-04
-5.6615E-04 -4.4250E-04 -5.6615E-04 -5.2070E-04
-5.6702E-04 -4.3297E-04 -5.6702E-04 -5.2156E-04
-5.6821E-04 -4.2725E-04 -5.6821E-04 -5.2274E-04
-5.6960E-04 -4.2343E-04 -5.6960E-04 -5.2411E-04
-5.7094E-04 -4.1962E-04 -5.7094E-04 -5.2545E-04
-5.7196E-04 -4.0817E-04 -5.7196E-04 -5.2646E-04
-5.7246E-04 -4.0817E-04 -5.7246E-04 -5.2695E-04
-5.7237E-04 -3.9673E-04 -5.7237E-04 -5.2685E-04
-5.7179E-04 -3.8910E-04 -5.7179E-04 -5.2626E-04
-5.7093E-04 -3.7384E-04 -5.7093E-04 -5.2538E-04
-5.7003E-04 -3.8147E-04 -5.7003E-04 -5.2448E-04
-5.6935E-04 -3.7384E-04 -5.6935E-04 -5.2379E-04
-5.6903E-04 -3.5858E-04 -5.6903E-04 -5.2346E-04
-5.1605E-04 -3.8528E-04 -5.1605E-04 -4.7059E-04
-5.1687E-04 -3.7766E-04 -5.1687E-04 -4.7141E-04
-5.1817E-04 -3.7384E-04 -5.1817E-04 -4.7270E-04
-5.1972E-04 -3.7575E-04 -5.1972E-04 -4.7423E-04
————————————————————————————

5. 应力分布求解
得到每个三角形位移之后,我需要求解每个单元的应力情况,这里采用公式进行计算,求出的应力就是三角形单元任意一点的应力,这里需要注意的是,相邻三角形单元的应力值是不连续的。为了将这种应力以图的形式直观展示出来,还需要做一些工作,Matlab中含有一个colorbar功能,我们选择jet形式,这和仿真软件得到的结果形式一致,其他的设置在代码中有详细解释,这里不再赘述,最后,得到的应力分布图如下图5、6所示。(对应程序文件为stress.m)通过图像可以看出,Ansys结果更为漂亮,因为其网格有加密,但是两者的变化趋势几乎一致。

function [Fx,Fy,Fxy]=stress(B,Node_information,Element_information,D,U)
Num_element=size(Element_information,1);
Num_node=size(Node_information,1);
Fx=[];
Fy=[];
Fxy=[];
for i = 1:Num_elementNode_num=Element_information(i,2:4);u=[U(Node_num);U(Node_num+Num_node)];f=D*B(3*i-2:3*i,:)*u;Fx=[Fx;f(1)];Fy=[Fy;f(2)];Fxy=[Fxy;f(3)];
endend


matlab计算 x y 方向应力


ansys 计算 x y 方向应力

6. 结果分析
凭我自己的直观想象,这个问题的模型对称,荷载对称,得到的应力结果和位移结果都应该是对称的,和我得到的结果显然一致,而且在孔洞处又很明显的应力集中的状态,而且随着网格的加密,其结果会往Ansys得到的结果逼近。所以我编成得到结果很好的和Ansys吻合在了一起。
但是限于自己的水平,代码有的地方比较冗杂繁复,而且计算速度明显低于成熟的商业软件Ansys,我也要继续努力,为了让我们自己的国家也有一套成熟的有限元软件而努力。

matlab编译平面有限元计算(附有完整代码)相关推荐

  1. 爬虫基础:python实现爬取无水印某瓜视频(附有完整代码,超详细)

    文章目录 一.前言 二.爬无水印的某瓜视频 1.分析网站 2.完整代码 三.总结 一.前言 爬虫真的很尴尬,稍微写点文章,分析网站什么的,就不给过,版权问题,哎,我会在边缘疯狂试探,一定要写详细点,让 ...

  2. 【通信】基于Matlab实现延时波束形成附完整代码

    1 内容介绍 现代社会发展要求通信系统功能越来越强,性能越来越高,构成越来越复杂;另一方面,要求通信系统技术研究和产品开发缩短周期,降低成本,提高水平.这样尖锐对立的两个方面的要求,只有通过使用强大的 ...

  3. 封装uniapp-uni-table组件,获取点击行事件,可传入自定义表头,传入后端数据,获取多选数据(其中行点击事件只有H5端可以用)附有完整代码

    功能如下:↓ 传值如下↓ data() {return {tableParam: {type: "selection", //是否可以多选loading: false,//加载状态 ...

  4. 基于MATLAB的求解线性方程组(附完整代码和例题)

    目录 前言 一. 直接求解:矩阵除法 例题1 例题2 例题3 二. 直接求解:判断求解 2.1 m=n且rank(A)=rank(C)=n 2.2 rank(A)=rank(C)=r<> ...

  5. 【影像配准】遥感影像配准结果输出tif影像(附有完整代码)

            遥感影像经过配准后会输出两幅影像:配准后的参考影像和待配准影像:         注:png 转 tif         因为两幅多时像影像在拍摄时即便是同源也无法保证拍摄到的景象完全 ...

  6. 高斯正反算—投影坐标转大地坐标、大地坐标转投影坐标(附有完整代码及测试结果)

            本文意在介绍高斯正反算的基本原理和代码实现!针对高斯正反算网上给出的方法很多,但是我试了之后发现多少都有些问题:或公式原理问题.或精度问题!         通过查找资料对其进行了总结 ...

  7. 入门C语言第三话:数组之实战篇——扫雷(进阶版——图形化界面,递归展开,播放音乐与音效,标记取消雷,记录雷的个数,鼠标点击,文末附有完整代码)

    文章目录 前言 每日鸡汤 基本思路 衔接基础班扫雷 准备阶段 正文 一.雷盘信息的存储 1.设置雷盘11*11与初始化 2.放置雷 3.放置雷周围的信息 二.图形化界面 1.创建与初始化窗口 2.加载 ...

  8. BP神经网络优化 | MATLAB基于飞蛾扑火算法优化BP神经网络(MFO-BP)的预测模型(完整代码在文末)

    飞蛾扑火( Moth-flame optimization algorithm,MFO) 是Seyedali Mirjalili等于2015年提出的一种新型智能优化算法.该算法具有并行优化能力强,全局 ...

  9. 【Matlab开发】MATLAB编译C/C++代码

    在使用MATLAB编译C/C++代码时,C/C++代码中要使用一个mexFunction函数,那么这个函数是如何定义,在编译时又是如何实现的呢?下面我将使用实例进行说明. 如一个简单的函数: doub ...

  10. MATLAB模拟陀螺仪的运动轨迹(附完整代码)

    本文讲述了陀螺仪运动轨迹的分析过程,并且通过MATLAB进行仿真. 顺时针旋转 MATLAB陀螺仪运动轨迹仿真(正转) 逆时针旋转 MATLAB陀螺仪运动轨迹仿真(反转) 第一步:建立坐标系 假设陀螺 ...

最新文章

  1. C语言的单链表求交点
  2. c 应用程序多语言版本,c – 在win32 API应用程序中实现全球化/多语言功能
  3. 澎思科技马原:AI安防竞争还未结束,落地进入后发优势时代 | MEET2020
  4. oracle中rownum和row_number()的区别
  5. syslog发送日志而docker容器接收不到的问题
  6. 分枝定界法解0/1背包问题
  7. 台达ms300变频器使用手册中文_台达变频器:满足未来驱动需求
  8. Jquery插件 bootstrap-datepicker 日期拾取器
  9. python爬取国内代理ip_python爬虫实战:爬取西刺代理的代理ip(二)
  10. Python 加性高斯白噪声 AWGN
  11. IOS开发把汉字转换成拼音的两种方法和返回拼音首字母
  12. 移动硬盘linux系统安装win7系统,超简单的移动硬盘安装系统win7教程
  13. 父债子偿有法可依吗?可法院却对这个案子说:不!
  14. 数据库架构设计——数据库选型
  15. 2023程序员找工作难?盘点目前IT各大热门行业,看看哪些更有前景
  16. python爬虫:获取菜鸟网站上url
  17. 2021最新上海互联网公司排名
  18. WebRTC-Failed to set remote answer sdp: The order of m-lines in answer doesn‘t match order in offer.
  19. python中int和input的区别_python中input()与raw_input()的区别分析
  20. 腾讯汤道生:2020年加大投入产业互联网生态建设

热门文章

  1. EasyDarwin —— ubuntu搭建rtsp服务,使用FFmpeg进行rtsp推拉流
  2. 微信公众号添加html,微信公众号如何在文章里添加超链接的方法教程
  3. Mac绿联USB转以太网无法连接解决方法
  4. 这些年,我工作上走过的路
  5. 前端工程化——Livereload和HMR、本地开发服务器
  6. 线性代数matlab的心得体会,关于线性代数心得体会
  7. cat3速度 rj45_综合布线当中,CAT8网线开始渐入佳境
  8. GIS | 坐标系统与地图投影
  9. 记录一些比较常用的简单jsp模板
  10. 如何在SVN创建分支版本