0
  • 聊天消息
  • 系统消息
  • 评论与回复
登录后你可以
  • 下载海量资料
  • 学习在线课程
  • 观看技术视频
  • 写文章/发帖/加入社区
创作中心

完善资料让更多小伙伴认识你,还能领取20积分哦,立即完善>

3天内不再提示

基于Matlab有限元编程的变截面悬臂梁分析

8XCt_sim_ol 来源:仿真秀App 作者:SimPC 2022-09-08 11:11 次阅读

导读:大家好,我是SimPC博士,主要从事工程结构抗震及减隔震研究,玻璃成型热工设备流动及传热研究,玻璃材料力学性能研究。精通有限元等数值算法的实现,有限元软件二次开发,数据处理,偏微分方程求解,优化算法,GUI界面开发等。有多项科研成果,其中SCI论文4篇,EI3篇,专利2篇。



近日我注册并认证仿真秀专栏,将在仿真秀官网和App给平台用户带来Matlab有限元编程复杂函数拟合和matlab绘图相关内容。此外还会带来隔震建筑Abaqus建模仿真分析等内容。本次案例主要以受均布荷载和集中荷载的变截面悬臂梁为研究对象,通过matlab编制四节点和八节点四边形单元有限元程序来对悬臂梁进行受力分析。



一、问题概述


如图1-1 所示,某变截面悬臂梁长度为2m,截面面积由0.6m至0.2m线性变化,受作用在自由端节点的集中荷载2P=kN和竖直方向均布荷载q=1kN/m作用,按平面应力问题分析,求解自由端节点挠度。变截面悬臂梁采用C30混凝土,弹性模量为E= 4 3 10 MPa,泊松比为。编制四节点和八节点四边形单元有限元程序,最终得到梁的变形。a9ce4008-2e95-11ed-ba43-dac502259ad0.png


图1-1 变截面悬臂梁


二、求解思路


对于本问题采用基于MATLAB 编制有限元分析程序进行求解,其基本组成部分包括前处理模块、分析主程序模块和后处理模块。在前处理模块中,实现节点坐标输入、单元节点编号、网络划分以及边界条件输入等工作;在分析主程序模块中,求解整体刚度方程;在后处理模块中,实现结果显示、数据输出等工作。本文主要针对四节点四边形单元与八节点四边形单元理论和对应的计算程序进行讲解。

有限元法的基本步骤:
  • 几何域离散,获得标准化的单元;

  • 通过能量原理(虚功原理或最小势能原理,获得单元刚度方程;

  • 单元的集成(装配);

  • 处理位移边界条件;

  • 计算支反力;

  • 计算单元的其他物理量(应力应变)。
这几步中,最核心的内容是单元研究,具体包括:
  • 节点描述

  • 场描述

  • 单元刚度方程。
接下来的内容主要以单元的描述为核心内容结合matlab代码,为大家讲解本案例有限元matlab编程过程。



1、平面问题的平衡方程、几何方程、物理方程


平面问题的弹性力学基础理论是推导有限元方程的基础,所以先罗列出平面问题的平衡方程、几何方程、物理方程,具体如公式(1)-(3)所示。至于这些方程的推导过程大家可以参考任意弹性力学课本,都会进行详细的讲解。

a9f0193a-2e95-11ed-ba43-dac502259ad0.png


2、等参单元


在有限元方法中,若要离散边界为曲线或曲面的求解域,需要建立将形状规则的单元变换为边界为曲线或曲面的单元的方法,在有限元法中对应此问题所采用的变换方法是等参变换,即单元几何形状的变换和单元内长函数采用相同数目的节点及相同的插值函数进行变换。同样我们今天要讲的四边形单元也从其对应的等参单元的基础理论讲起。四边形单元可以由自然坐标系中的矩形单元映射而成,映射关系如图2-1所示。

aa0ae3be-2e95-11ed-ba43-dac502259ad0.png

图2-1 平面四节点矩形单元的映射关系

在自然坐标系下,矩形单元是规则化的,当自然坐标系中的单元取为双线性单元时(也即为四节点四边形单元),平面四节点矩形单元如图2-2所示,单元有4个节点,8个自由度。单元的形函数定义如下:aa232d02-2e95-11ed-ba43-dac502259ad0.png     



(4)
其中,aa373afe-2e95-11ed-ba43-dac502259ad0.jpgaa4b46ac-2e95-11ed-ba43-dac502259ad0.jpg为自然坐标系下的节点坐标值。单元从自然坐标系到物理坐标系的映射为

aa68605c-2e95-11ed-ba43-dac502259ad0.png

在进行映射变换时候,要求单元两个坐标系下的节点编号要对应。单元的节点变量用型函数进行插值,有aa7cef9a-2e95-11ed-ba43-dac502259ad0.png                   (7)
function N=ShapeFun(s,t)            N1=1/4*(1-s)*(1-t);N2=1/4*(1+s)*(1-t);N3=1/4*(1+s)*(1+t);N4=1/4*(1-s)*(1+t);N=[N1 0 N2 0 N3 0 N4 0;0 N1 0 N2 0 N3 0 N4];end

同理平面八节点矩形单元如图2-3所示,单元共有8个节点,16个自由度。单元的形函数定义如下:

aa8457a8-2e95-11ed-ba43-dac502259ad0.png  (8)


aa99b436-2e95-11ed-ba43-dac502259ad0.png       (9)

aaab5560-2e95-11ed-ba43-dac502259ad0.png   (10)

其中,aa373afe-2e95-11ed-ba43-dac502259ad0.jpgaa4b46ac-2e95-11ed-ba43-dac502259ad0.jpg为自然坐标系下的节点坐标值。单元从自然坐标系到物理坐标系的映射为

aadf56a8-2e95-11ed-ba43-dac502259ad0.png      (11)

aaed8980-2e95-11ed-ba43-dac502259ad0.png        (12)


在进行映射变换时候,要求单元两个坐标系下的节点编号要对应。单元的节点变量用型函数进行插值,有

ab017364-2e95-11ed-ba43-dac502259ad0.png   (13)
function N=ShapeFun(s,t)            %% 四边形八结点等参单元形函数矩阵 % 角点N1=1/4*(1-s)*(1+t)*(-s+t-1); N2=1/4*(1-s)*(1-t)*(-s-t-1); N3=1/4*(1+s)*(1-t)*(s-t-1); N4=1/4*(1+s)*(1+t)*(s+t-1); % 边中点 N5=1/2*(1-t^2)*(1-s); N7=1/2*(1-t^2)*(1+s); N6=1/2*(1-s^2)*(1-t); N8=1/2*(1-s^2)*(1+t); N=[N1 0 N2 0 N3 0 N4 0 N5 0 N6 0 N7 0 N8 0;0 N1 0 N2 0 N3 0 N4 0 N5 0 N6 0 N7 0 N8];

ab0b591a-2e95-11ed-ba43-dac502259ad0.png

图2-2 平面四节点矩形单元

ab27186c-2e95-11ed-ba43-dac502259ad0.png

图2-3 平面四节点矩形单元等参单元中除了完成如公式(5)(6)(10)(11)的坐标映射外,还需要完成坐标偏导数的映射和面积/体积的映射,因为在最终推导出的单元刚度矩阵表达式,即一个积分函数中会包含坐标的偏导项和坐标的面积积分项,如公式(x)所示,所以接下来我们研究坐标偏导项的映射关系。根据链式求导法则,形函数对自然坐标系的导数为


ab3f5c4c-2e95-11ed-ba43-dac502259ad0.png (14)


写成矩阵的形式就是


ab4b8db4-2e95-11ed-ba43-dac502259ad0.png     (15)



其中,J被称为Jacobi矩阵。反过来,形函数对物理坐标的导数为 ab5853aa-2e95-11ed-ba43-dac502259ad0.png         (16)



另外,对于二维平面单元还要完成面积的映射,为 ab71400e-2e95-11ed-ba43-dac502259ad0.png              (17)


可以看出Jacob矩阵在等参变化中扮演着至关重要的角色,Jacob矩阵具体的表达式如下所示, ab9749ca-2e95-11ed-ba43-dac502259ad0.png       (18)


公式18对应的八节点单元雅各比矩阵的求解代码为:
function J=Jacobi(ie,s,t,Elements,Nodes)            ENodes = Elements(ie,:);                 %获取单元结点 xe = Nodes(ENodes(:),:);                 %获取节点坐标 x1=xe(1,1);y1=xe(1,2); x2=xe(2,1);y2=xe(2,2); x3=xe(3,1);y3=xe(3,2); x4=xe(4,1);y4=xe(4,2); J=1/4*[-(1+t) -(1-t) 1-t 1+t;1-s -(1-s) -(1+s) 1+s]*[x1 y1;x2 y2;x3 y3;x4 y4];end
公式18对应的四节点单元雅各比矩阵的求解代码为:
function J=Jacobi(ie,kesi,yita,Elements,Nodes)            ENodes = Elements(ie,:);                 %获取单元结点 xe = Nodes(ENodes(:),:);                 %获取结点坐标 x1=xe(1,1);y1=xe(1,2); x2=xe(2,1);y2=xe(2,2); x3=xe(3,1);y3=xe(3,2); x4=xe(4,1);y4=xe(4,2); J=1/4*[-(1-yita),(1-yita),(1+yita),-(1+yita);-(1-kesi),-(1+kesi),(1+kesi),(1-kesi)]*[x1 y1;x2 y2;x3 y3;x4 y4];end
3、刚度矩阵的推导

为了求出上述平面四节点和八节点单元的单元刚度矩阵,需要借助能量原理(虚功原理、最小势能原理)进行推导,能量原理的推导过程大家可以参考任意一本有限元理论书籍,都会有详细的推导过程,这里就不做进一步推导讲解,直接给出物理坐标和几何坐标系下的刚度矩阵的公式

aba86c64-2e95-11ed-ba43-dac502259ad0.png   (19)

abb64028-2e95-11ed-ba43-dac502259ad0.png   (20)

其中B矩阵为应变矩阵,如下式;D矩阵为材料刚度矩阵,如公式(1)所示,是物理方程中表征应力应变关系的矩阵。从上述刚度矩阵的表达式可以看出,自然坐标和物理坐标间要完成坐标映射、偏导映射、面积隐射三个部分,具体映射公式已在上一节的等参单元讲解中详细给出。

abcd37e2-2e95-11ed-ba43-dac502259ad0.png        (21)

4、高斯积分

公式(20)中的单元刚度矩阵通过数值积分求得,本案例中的四节点和八节点四边形等参单元均采用2*2个积分点的高斯积分即可求得精确结果。高斯积分点的坐标具体如图所示。

abd6de14-2e95-11ed-ba43-dac502259ad0.png

4-1 Gauss积分点示意图

公式(20)写成数值积分的形式为

ac0adae8-2e95-11ed-ba43-dac502259ad0.png        (22)

对于8节点单元实现上述数值积分的代码如下所示:

r = [-sqrt(1/3) sqrt(1/3)];             % 2*2 高斯积分点 s = [r(1) r(1) r(2) r(2)]; t = [r(2) r(1) r(1) r(2)];              % 高斯积分点坐标for i=1:4         J = Jacobi(E_ID,s(i),t(i),Elements,Nodes);             % 雅可比矩阵         Nst = DiffShapeFun(s(i),t(i));        % 形函数关于自然坐标s,t求导         Nxy = zeros(8,2);         for j=1:8           Nxy(j,:) = (JNst(j,:)')';             % 形函数关于 x,y 求导=inv(J)*Nst         end         Bm = [Nxy(1,1) 0 Nxy(2,1) 0 Nxy(3,1) 0 Nxy(4,1) 0 Nxy(5,1) 0 Nxy(6,1) 0 Nxy(7,1) 0 Nxy(8,1) 0;             0 Nxy(1,2) 0 Nxy(2,2) 0 Nxy(3,2) 0 Nxy(4,2) 0 Nxy(5,2) 0 Nxy(6,2) 0 Nxy(7,2) 0 Nxy(8,2);             Nxy(1,2) Nxy(1,1) Nxy(2,2) Nxy(2,1) Nxy(3,2) Nxy(3,1) Nxy(4,2) Nxy(4,1) Nxy(5,2) Nxy(5,1) Nxy(6,2) Nxy(6,1) Nxy(7,2) Nxy(7,1) Nxy(8,2) Nxy(8,1)];         ke = ke+det(J)*Bm'*D*Bm*Width;  %数值积分  end

5、均布荷载的施加

在有限元中分布力要转为等效节点荷载,而确定等效节点荷载的方法也是通过能量原理推导得到

ac2b29d8-2e95-11ed-ba43-dac502259ad0.png       (22)

上式中,第一项代表体积力的等效荷载,第二项代表面积力的等效荷载,这个案例我们只考虑面力荷载。实现公式22的代码为

function Pe=UniLoad(ie,N_ID_p1,q0,Nodes,Elements)     k=-0.625e-3;                            % 均布荷载值 N/mms = [-sqrt(1/3) sqrt(1/3)];                 % 2*2 高斯积分点ENodes = N_ID_p1(ie,:);                    %获取单元结点号Pe=zeros(16,1);                          %生成临时单元节点力零列向量x1=Nodes(ENodes(1),1);x6=Nodes(ENodes(4),1);L16=abs(x6-x1);                          %单元长度for i=1:2                                 %用于高斯积分的求和循环    N_q=ShapeFun(s(i),1);                   % 4级子程序:ShapeFun(s(i),1)    q_x=q0;    Pe=Pe+N_q'*q_x*[0;L16/2];            endend

三、Matlab有限元编程精品课

网格划分及变形结果如图3-1所示。本案例的详细视频教程和对应的matlab源码,请关注我的仿真秀官网和APP精品课程Matlab有限元编程从入门到精通10讲》。

ac4745fa-2e95-11ed-ba43-dac502259ad0.png

图3-1 梁变形结果

审核编辑:汤梓红

声明:本文内容及配图由入驻作者撰写或者入驻合作网站授权转载。文章观点仅代表作者本人,不代表电子发烧友网立场。文章及其配图仅供工程师学习之用,如有内容侵权或者其他违规问题,请联系本站处理。 举报投诉
  • matlab
    +关注

    关注

    175

    文章

    2908

    浏览量

    228325
  • 编程
    +关注

    关注

    88

    文章

    3431

    浏览量

    92218

原文标题:教你Matlab有限元编程对悬臂梁进行受力分析-附源码及教程

文章出处:【微信号:sim_ol,微信公众号:模拟在线】欢迎添加关注!文章转载请注明出处。

收藏 人收藏

    评论

    相关推荐

    有限元分析

    `<p><font face="Verdana">有限元分析</font>
    发表于 06-17 10:50

    【PDF】matlab有限元法计算分析程序编写

    【PDF】matlab有限元法计算分析程序编写附件:
    发表于 02-28 11:04

    MATLAB有限元分析与应用

    有限元分析和应用。另外,本书还提供了大量免费资源。 第1章引言 1.1有限元方法的步骤 1.2用于有限元分析MATLAB函数 1.3MATLAB
    发表于 02-28 11:07

    悬臂梁冲击试验机原理

    悬臂梁冲击试验机冲击试验时,摆锤从垂直位置挂于机架扬臂上,把扬臂提升一扬角α,摆锤就获得了一定的位能。释放摆锤,让其自由落下,将放于支架上的样条冲断,向反向回升时,推动指针,从刻度盘读数读出冲断试样
    发表于 11-28 14:01

    【TL6748 DSP申请】检测悬臂梁传感器的振动频率并进行锁相,计算品质因数

    做的是基于MEMS的悬臂梁传感器测试气体的浓度,通过不同的气体在悬臂梁上的吸附质量不同,悬臂梁的谐振频率会发生变化,采集数据计算出频率的变化可以通过计算得出空气中所要检测的气体的浓度。之前我们
    发表于 09-09 16:59

    一种电磁型射频微机电系统开关的软磁悬臂梁制备工艺研究

    吴军民a, 汤学华b (上海电机学院a.汽车学院; b.机械学院,上海200245) 摘 要:研究了一种电磁型射频(RF)微机电系统(MEMS)开关的软磁悬臂梁的制备工艺。为了得到适合MEMS器件
    发表于 07-04 08:14

    悬臂梁式传感器的原理是什么?

    悬臂梁称重传感器的典型应用包括地秤,平台秤,料斗秤,吊车秤,飞机称重传统的杠杆系统规模的系统和转换,适用于固体、液体流动秤、人体秤、配料秤、案秤及精细化工配比秤。
    发表于 11-04 09:11

    如何有效的学习CAE有限元分析

      有限元分析在工程应用中正发挥着越来越重要的作用,尤其在复杂钢结构的创新设计中,离开有限元分析就很难使结构的强度、刚度以及稳定性等得到有效地保证。因而,为了更好地增加产品结构的可靠性,越来越多
    发表于 07-07 16:59

    求一种有限元分析中PCBA的简化建模方法

    本文基于某车型门窗控制器(DCM:Door Control Module)的PCBA提出一种有限元分析中PCBA的简化建模方法,并进行有限元仿真模态分析。通过仿真模态分析结果与试验模态
    发表于 04-19 06:20

    有限元法的原理

    1有限元法原理有限元法是以分原理为基础的一种数值计算方法。把所要求解的边值问题转化为相应的分问题,利用剖分、插值和离散化,把分问题转化
    发表于 09-06 08:08

    利用微电铸制作镍悬臂梁

    利用微电铸制作镍悬臂梁:为了设计出导电的Ni 微变截面悬臂梁,采用了基于成熟的电镀技术的微电铸工艺. 在分析悬臂梁的变
    发表于 12-29 23:52 7次下载

    静电式微开关硅悬臂梁的变形分析

    静电式微开关硅悬臂梁的变形分析:介绍了一种计算微开关的硅悬臂梁在电场力作用下的变形的方法, 由于作用在梁上的载荷随着梁的变形而变化,用积分法计算存在相当困难,以下提出
    发表于 01-01 11:35 21次下载

    测量仪悬臂梁拓扑优化

    针对高精度衡器载荷测量仪反力机构的悬臂梁质量较大的问题,对高精度衡器载荷测量仪的工作原理、悬臂梁加工工艺和拓扑优化方法进行了研究。利用ANSYS对优化前的反力机构进行了静力学分析,基于固体各向同性
    发表于 03-27 11:10 0次下载

    如何使用Matlab进行有限元分析和硬盘的选购与使用资料说明

    介绍了Matlab 语言特点,给出了Matlab 环境下实现有限元的步骤。并以一个具体实例说明如何使用Matlab 进行有限元分析。使用该方
    发表于 07-24 16:27 5次下载
    如何使用<b class='flag-5'>Matlab</b>进行<b class='flag-5'>有限元分析</b>和硬盘的选购与使用资料说明

    基于六面体单元热应力问题的Matlab有限元编程求解

    导读:上一篇《弹性地基梁matlab有限元编程,以双排桩支护结构计算为例》引起了Matlab有限元编程
    的头像 发表于 11-17 11:10 1818次阅读