首页/文章/ 详情

基于matlab的二维杆单元开发

1月前浏览825

目录

  1. 概述
  2. 杆单元为什么要引入局部坐标
  3. 局部坐标系下的刚度矩阵
  4. 全局坐标系下的刚度矩阵
  5. 程序编写&算例
  6. 所有的计算代码

概述

   杆系结构是工程中应用较为广泛的结构体系,包含了平面或空间形式的梁、钢架和行桁架等等,这些结构的每个单独的杆件都可以视为一个单独的单元,单元与单元之间用节点连接,根据结构体系的不同,节点可以传递轴向力、切向力或者弯矩等等。
   帖子以二维平面杆单元为例,采用材料力学更直观的方法建立杆单元的刚度矩阵,详细介绍了具体的推导过程,最后以一个桁架为例讲解matlab程序并采用abaqus做对比计算,文末给出了详细的程序

杆单元为什么要引入局部坐标系

  我们在做一件事情的时候,一定要知道为什么做
  杆单元为什么不直接在全局坐标系中求解,非得人为的引入一个局部坐标系?我想第一次学习杆单元的人都会有这个疑问,然后很遗憾,单单用杆单元部分的知识并不能很好的回答这个问题。一个简单的答案是:为了统一单元、最终是为了简化计算
   我们知道,科学家最喜欢将特殊的东西泛化,然后给出一个简洁美丽的公式来解释广泛的事物。在科学计算中也是如此,就拿这次的杆单元来说,我们脑子里面浮现的是一个结构体系中横七竖八的段杆,就像是体育场上上空的架子,乱糟糟的,用数学的语言描述就是他们的空间位置都不相同,这在直觉上意味着极大的工作量。
   这时候,天马行空的想法随着工作量的增加会自然而然的从脑子里面冒出来(很多发明一开始都是天马行空的偷懒想法):如果我们能只求一次某个量,然后乘以一个系数把这个量映射到所有的杆单元上面就好了,这个想法就是引入局部坐标系的火苗,有了这个思想的指导,后面的数学家就开始干活了,最终的杆单元计算过程就是上面那句加粗的文字那样:首先在局部坐标系下计算刚度矩阵,然后乘以一个坐标转换矩阵,将据局部坐标系下的刚度矩阵映射到全局坐标系中参与组装。下面列出一些引入局部坐标系具体方便了哪些

  1. 简化计算    
       局部坐标系可以简化杆单元的几何和材料属性的描述,使得在进行刚度矩阵和质量矩阵的计算时更为方便。通过将杆单元的方向与局部坐标系对齐,可以使得计算过程中的数学表达更为直观。
  2. 方便载荷施加    
       在局部坐标系中,可以更方便地定义和施加载荷。比如,集中载荷或分布载荷可以直接在杆单元的局部坐标系中施加,避免了转换成全球坐标系的复杂性。
  3. 更好地处理变形和应力分析    
       杆单元在变形时会有不同方向的应力和应变,局部坐标系可以帮助明确杆单元的拉伸、压缩及剪切行为,进而提高计算精度。
  4. 支持多方向的连接    
       在复杂的结构中,杆单元可能以不同的方向连接。在局部坐标系中,每个单元的方向都是独立定义的,这使得处理复杂结构的连接关系更为灵活。
  5. 易于与其他单元类型结合    
       在实际工程应用中,结构通常由多种单元组成。引入局部坐标系可以更方便地将杆单元与其他类型的单元(如面单元或体单元)结合,确保它们之间的力学性质和边界条件协调一致。    
      我在这里只是从定性的角度给出了一些解释,要想彻底的从底层理解这一套工作思路,还需要学习等参元,等参元的引入把简化计算的思想发挥的淋漓尽致,在等参元中,科学家通过坐标变换解决了积分难题,

局部坐标系下的刚度矩阵

   材料力学给出了二维杆单元的受力特点

  1. 每个杆单元的节点只有两个平动自由度
  2. 二维杆单元的节点不传递弯矩
  3. 外荷载均为作用在节点上的集中力    
       基于此,我们从一个结构体系中取出一个杆件,并引入与杆轴向重合的局部坐标系,设杆单元的截面积为      ,长度为      ,两端的节点分别用      和      表示,并且从节点      到节点      的方向为局部有坐标轴      ,垂直于杆轴向为      ,具体的示意图为      根据材料力学的知识,对两端节点列方程
 
 

  整理成矩阵形式

 

  以上就是二维杆单元在局部坐标系中的刚度矩阵,下面推导局部坐标系与全局坐标系的转化矩阵,即全局坐标系的刚度矩阵。

全局坐标系下的刚度矩阵

  推导全局坐标系下的刚度矩阵的重点就是推导转换矩阵,即轴向量与全局坐标系之间的转换关系,如下图展示了轴向数据与全局坐标系数据   通过上面图片可知,整体坐标系下,杆单元两端节点的位移矩阵为

 

  局部坐标系下,杆单元两端节点的位移矩阵为

 

   根据三角函数关系,写出全局坐标系位移与局部坐标系位移之间的关系为

 
 

   整理成矩阵形式

 

  其中,转换矩阵    

 

  其中

 

  上面就是局部坐标系向全局坐标系转换数据的矩阵,局部和全局量尊遵循上述转换关系。对于全局坐标系下的刚度矩阵,则有

 

  结合上文中的局部坐标系下的刚度矩阵,代入转换矩阵有

 

  好了,现在有了全局坐标系下的杆单元刚度矩阵,下面就可以写程序了。

程序编写&算例

  这里直接给出一个算例,以算例的形式详细讲解程序。二维杆单元的弹性模量    ,截面面积    ,具体的网格和边界条件示意图为  上面的桁架一共15个杆单元,单元编号用斜体字表示;9个节点,节点用有背景色的方块表示,在节点1和6上施加竖直向的集中力,力的大小    ,节点9和4为的两个平动自由度被限制。
   下面开始介绍关键的代码,全部的代码在文末给出。
   首先是杆单元的材料属性,杆单元只需要弹性模量和截面积。

E=210e3; % 弹性模量 (GPa=1000 N/mm^2)
A=500; % 杆单元截面积 (mm^2)

   然后给出模型的信息。

node=9; % 模型节点总数
element=15; % 杆单元的数目
nodecord=[  1.5,  0.866025388;
  2.5,  0.866025388;
   2.,           0.;
   3.,           0.;
   1.,           0.;
  0.5,  0.866025388;
   0.,           0.;
 -0.5,  0.866025388;
  -1.,           0.];
xcord=nodecord(:,1); % 模型的x坐标
ycord=nodecord(:,2); % 模型的y坐标
elementcon=[1, 2   ; 
2, 3    ;
3, 4    ;
4, 2    ;
1, 5    ;
6, 1    ;
6, 7    ;
7, 5    ;
 3, 1   ;
 5, 3   ;
 5, 6   ;
 8, 6   ;
 7, 8   ;
 9, 7   ;
 8, 9   ];
Wd =2*node; % 模型的整体自由度

   下面计算杆单元刚度矩阵和施加节点力。计算杆单元的全部代码见文末

% 计算杆单元刚度矩阵
[S]=planetruss(E,A, Wd,element,elementcon,xcord,ycord);
% 施加荷载 (units: N)
F=zeros(Wd,1);
% 给第二个自由度施加45e3N的力,即给节点1的y方向施加45e3N的力
F(2)=-45e3; 
% 给第十二个自由度施加45e3N的力,即给节点6的y方向施加45e3N的力
F(12)=-45e3;

   下面给节点4和9施加固定边界条件。

%% 施加固定位移条件
gdof=1: Wd;
cdof=[7 8 17 18];

   求解方程,这里是静力计算,实质上的计算量就是矩阵求逆。

[U]=solve(Wd,gdof,cdof,S,F);
displacements=[gdof' U]

   下面是一些后处理的程序,以注释的形式给出解释说明。

%% 统计变形后的节点坐标
ux=1:2:2*node-1;
uy=2:2:2*node;
UX=U(ux);
UY=U(uy);
newcordx=xcord+UX;
newcordy=ycord+UY;
newcord=[newcordx,newcordy]
%% 计算杆单元的力
[p]=trussforce(E,A,element,elementcon,xcord,ycord,U);
A=1:size(p,1);
Elemental_Force=[A' p]
%% 绘制杆单元
SC=100; %缩放因子,便以观察变形
Snewcordx=xcord+SC*UX;
Snewcordy=ycord+SC*UY;
Snewcord=[Snewcordx,Snewcordy];
drawplanetruss(element,elementcon,nodecord,newcord,Snewcord);

  运行上面的matlab程序,并将变形前后的桁架绘制在同一个窗口中观察变形,图形为   可以看到,变形的形状符合直觉,为了进一步验证程序的计算效果,采用abaqus创建同样的模型并对比计算结果,abaqus的模型为   abaqus计算完毕后,统计所有节点的竖直向位移与matlab程序对比,绘制表格为


matlababaqus
-0.003214286-0.00321429
-0.001071429-0.00107143
-0.002285714-0.00228571
0-4.50E-32
-0.003428571-0.00342857
-0.003214286-0.00321429
-0.002285714-0.00228571
-0.001071429-0.00107143
0-4.50E-32
   可以看到,上面的计算结果一模一样,跟造数据了一样,事实上就算是我造数据了,你们也不知道,但是按我的习惯,我必然会把所的计算文件放在文末

所有的计算代码

   下面是matlab的计算主程序,主程序调用的函数我就不放在这里了,有需要的可以私聊我,私聊必发!

clc ; clear ; % 清空变量
%%
% 杆单元材料参数
E=210e3; % 弹性模量 (GPa=1000 N/mm^2)
A=500; %杆单元截面积 (mm^2)
node=9; % 模型节点总数
element=15; % 模型单元总数
nodecord=[  1.5,  0.866025388;
  2.5,  0.866025388;
   2.,           0.;
   3.,           0.;
   1.,           0.;
  0.5,  0.866025388;
   0.,           0.;
 -0.5,  0.866025388;
  -1.,           0.];
xcord=nodecord(:,1); ycord=nodecord(:,2);
elementcon=[1, 2   ; 
2, 3    ;
3, 4    ;
4, 2    ;
1, 5    ;
6, 1    ;
6, 7    ;
7, 5    ;
 3, 1   ;
 5, 3   ;
 5, 6   ;
 8, 6   ;
 7, 8   ;
 9, 7   ;
 8, 9   ];
Wd =2*node; % Whole DOFs
% 计算杆单元刚度矩阵
[S]=planetruss(E,A, Wd,element,elementcon,xcord,ycord);
% 施加节点力 (units: N)
F=zeros(Wd,1);
F(2)=-45e3;
F(12)=-45e3;
%% 施加固定位移边界条件
gdof=1: Wd;
cdof=[7 8 17 18];
%% 求解位移
[U]=solve(Wd,gdof,cdof,S,F);
displacements=[gdof' U]
%% 统计变形后的坐标
ux=1:2:2*node-1;
uy=2:2:2*node;
UX=U(ux);
UY=U(uy);
newcordx=xcord+UX;
newcordy=ycord+UY;
newcord=[newcordx,newcordy]
%% 计算单元力
[p]=trussforce(E,A,element,elementcon,xcord,ycord,U);
A=1:size(p,1);
Elemental_Force=[A'
 p]
%% 绘制变形前后的对比图
SC=100; % 缩放因子,便于观察
Snewcordx=xcord+SC*UX;
Snewcordy=ycord+SC*UY;
Snewcord=[Snewcordx,Snewcordy];
drawplanetruss(element,elementcon,nodecord,newcord,Snewcord);

   下面是abaqus的计算inp文件,这个inp文件拿过来就能直接算,算完的结果跟上面的matlab一样,曲线吻合的跟假的一样,非常能提升自己搞科研的自信心,建议保存起来,每当其他程序出bug的时候运行一下这个程序,以此欺骗和麻痹自己。

*Heading
** Job name: Job-1 Model name: Model-1
** Generated by: Abaqus/CAE 6.14-2
*Preprint, echo=NO, model=NO, history=NO, contact=NO
**
** PARTS
**
*Part, name=Part-1
*Node
      1,          1.5,  0.866025388
      2,          2.5,  0.866025388
      3,           2.,           0.
      4,           3.,           0.
      5,           1.,           0.
      6,          0.5,  0.866025388
      7,           0.,           0.
      8,         -0.5,  0.866025388
      9,          -1.,           0.
*Element, type=T2D2
1, 1, 2
2, 2, 3
3, 3, 4
4, 4, 2
5, 1, 5
6, 6, 1
7, 6, 7
8, 7, 5
 9, 3, 1
10, 5, 3
11, 5, 6
12, 8, 6
13, 7, 8
14, 9, 7
15, 8, 9
*Nset, nset=Set-1
 1, 5, 6, 7, 8
*Elset, elset=Set-1
  6,  7,  8, 11, 13
*Nset, nset=Set-6, generate
 1,  9,  1
*Elset, elset=Set-6, generate
  1,  15,   1
*Nset, nset=fixed
 4, 9
*Nset, nset=middle
 1, 6
** Section: Section-1
*Solid Section, elset=Set-6, material=Material-1
500.,
*End Part
**  
**
** ASSEMBLY
**
*Assembly, name=Assembly
**  
*Instance, name=Part-1-1, part=Part-1
*End Instance
**  
*End Assembly
** 
** MATERIALS
** 
*Material, name=Material-1
*Elastic
210000.,0.
** 
** BOUNDARY CONDITIONS
** 
** Name: fixed Type: Displacement/Rotation
*Boundary
Part-1-1.fixed, 1, 1
Part-1-1.fixed, 2, 2
** ----------------------------------------------------------------
** 
** STEP: Step-1
** 
*Step, name=Step-1, nlgeom=NO
*Static
1., 1., 1e-05, 1.
** 
** LOADS
** 
** Name: middle   Type: Concentrated force
*Cload
Part-1-1.middle, 2, -45000.
** 
** OUTPUT REQUESTS
** 
*Restart, write, frequency=0
** 
** FIELD OUTPUT: F-Output-1
** 
*Output, field, variable=PRESELECT
** 
** HISTORY OUTPUT: H-Output-1
** 
*Output, history, variable=PRESELECT
*End Step



来源:有限元先生
ACTAbaqusMATLABADSUG材料
著作权归作者所有,欢迎分享,未经许可,不得转载
首次发布时间:2024-11-08
最近编辑:1月前
外太空土豆儿
博士 我们穷极一生,究竟在追寻什么?
获赞 22粉丝 2文章 62课程 0
点赞
收藏
作者推荐

newmark法求解运动方程(附matlab代码)

问题引入在数值仿真中,常常涉及到静力和动力的问题。静力问题通常忽略物体的惯性力和阻尼效应。也就是说,系统中的力只与物体的静态形变和外部载荷有关,而不考虑动态过程中速度和加速度带来的影响,因此静力问题可以用方程AX=B表示,这种方程用线性代数里面的内容就可以求解,主要是A矩阵求逆的计算速度决定了方程的求解速度,这部分的资料很多,在此不多说,这个帖子重点讲解动力问题。动力问题极为广泛,动力问题(DynamicProblem)是指在力学和物理学中研究物体或结构在外力作用下随时间变化的运动、变形和应力响应。与静力问题不同,动力问题考虑了系统的惯性力和阻尼效应,这使得其分析需要解决随时间变化的运动方程。动力问题广泛应用于涉及冲击、振动、地震、波动等现象的工程和物理系统中。动力问题的特点考虑惯性力和加速度:在动力问题中,物体或结构的运动和变形不仅受外力和内力的影响,还受到惯性力的作用。根据牛顿第二定律,惯性力与物体的加速度成正比,因此在动力分析中必须考虑加速度的影响。时间依赖性:动力问题的解通常是时间的函数,要求在整个时间域内求解物体或结构的运动状态。随着时间的推移,外力、位移、速度和应力等物理量会发生变化。系统响应的类型:自由振动:不受外部力作用的振动,仅依赖于系统的初始状态和自然特性。强迫振动:在外力作用下,系统的响应随外力变化而变化,通常伴随频率和振幅的变化。冲击和瞬态响应:系统受到瞬时力(如冲击或爆炸)时的快速响应,通常表现为瞬态振动。阻尼效应:动力系统中还可能存在阻尼,阻尼力通常与速度成正比,并会导致系统的运动逐渐衰减。动力问题中必须考虑阻尼效应对系统振动和响应的影响。动力问题的基本方程:动力问题的基本方程是牛顿第二定律或达朗贝尔原理的推广,称为运动方程(EquationofMotion)。一个简单的动力问题可以用以下形式的运动方程来描述:这个方程说明了外力F(t)如何引起系统的位移、速度和加速度的时间响应。这是一个二阶微分方程,求解这个方程的方法大致分为解析法和数值解法。解析解当然是最精确的,但只有简单的动力系统或边界条件较为规则的问题,动力问题可以通过解析方法求解。例如:单自由度系统的自由振动:对于简单的弹簧-质量系统,运动方程可以通过特征值问题解析求解,得到系统的固有频率和振动模态。正弦激励下的强迫振动:当系统受到周期性外力(如正弦波)作用时,可以通过解析方法求解系统的稳态响应。数值解法是对于复杂的几何形状、边界条件或载荷,解析解往往难以得到,因此通常采用数值方法来求解动力问题。常见的数值方法包括:有限元法(FEM):将连续体离散为有限个单元,通过求解离散的动力方程来求解整个结构的动态响应。Abaqus等软件广泛使用有限元方法进行复杂动力学分析。有限差分法(FDM):将动力方程的时间和空间导数进行离散化,用数值方法逐步推进时间步长,求解各个时间点的响应。显式和隐式时间积分法:显式方法:如中心差分法,适合处理瞬态问题和高频振动问题。隐式方法:如Newmark方法,适合求解稳定性要求较高的动力问题。这里就重点介绍newmark法求解运动方程,包括固定位移边界条件的施加。newmark理论及程序newmark法是一种常用的时间积分方法,用于求解结构动力学中的运动方程,尤其是在有限元分析中求解动力问题(如振动、冲击、地震响应等)。该方法由NathanM.Newmark于1959年提出,它通过离散时间步长,将动力方程转化为一系列可以逐步求解的代数方程,从而计算出位移、速度和加速度的时间响应。newmark法的优点在于它既适用于求解线性问题,也适用于求解非线性问题,并且它在保持数值稳定性和精度之间提供了一种灵活的平衡。根据参数的选择,Newmark法可以实现不同的时间积分方案,如条件稳定和无条件稳定。newmark公式资料非常的多,基本随便拎起来一本涉及到动力学的书籍都会给出该方法的详细公式,如下图所示上面的资料来自于朱伯芳院士的《有限单元法原理与应用》第四版的第499页。朱伯芳院士在北京于2024年2月4日逝世,每次我想起朱院士序言中的内容,我都会对朱院士充满敬意。朱院士为我国大江大河治理和水电开发作出了卓著贡献,我在此借这个帖子沉痛悼念朱伯芳院士!深切缅怀朱伯芳院士!朱院士的书是我学习有限元过程中推公式用的最多的一本书,当时我是看了一些国内其他书籍,觉得完蛋了,学不会有限元了,感觉山重水复疑无路了,但是看到朱院士的这本书,顿时豁然开朗,柳暗花明。这里贴一下书的封面,有需要的可以私信我。书归正传,有了newmark公式之后,先不慌着写程序,我们通过自己的理论求解了刚度矩阵、质量矩阵和阻尼矩阵之后,还要添加边界条件,才能消除刚体位移,至于添加边界条件又是一个很大的领域,我涉及不多,这里重点讲解newmark,就以最简单的固定位移边界条件为例。添加固定边界的方法有乘大数法、划零置一法和划行划列法等等,我在这里采用的是行划列法,即假如模型的第7个节点是全固定的,就把刚度矩阵等的第3*7-2:3*7(三维问题中)的行和列都直接删除,然后将矩阵传进newmark法中进行求解,求解完成后,再将自由度还原,方便后处理,这里贴一下主要的matlab程序。function[U,V,A]=newmark(deltat,beta,kk,mm,cc,ndim,fbc,kinc,ndofel,nodeu)%结果初始化uu=zeros(ndofel-ndim*length(nodeu),kinc);%didpvv=zeros(ndofel-ndim*length(nodeu),kinc);%velaa=zeros(ndofel-ndim*length(nodeu),kinc);%acc%%基本参数计算alpha=0.25*(0.5+beta)*(0.5+beta);%newamrk计算参数a0=1/(alpha*deltat*deltat);a1=beta/(alpha*deltat);a2=1/(alpha*deltat);a3=1/2/alpha-1;a4=beta/alpha-1;a5=deltat/2*(beta/alpha-2);a6=deltat*(1-beta);a7=beta*deltat;%%迭代计算ek=a0*mm+a1*cc+kk;%有效刚度矩阵ekk=inv(ek);%有效刚度矩阵求逆,只进行一次fori=2:kincfprintf("kinc:%dtime:%f\n",i,deltat*(i-1));ft=fbc(:,i);u=uu(:,i-1);v=vv(:,i-1);a=aa(:,i-1);%计算i+1时刻的有效荷载ft=fbc(:,i);f=ft+mm*(a0.*u+a2.*v+a3.*a)+cc*(a1.*u+a4.*v+a5.*a);%计算i+1时刻的位移uu(:,i)=ekk*f;%计算i+1时刻的速度和加速度aa(:,i)=a0.*(uu(:,i)-u)-a2.*v-a3.*a;vv(:,i)=v+a6.*a+a7.*(aa(:,i));end%%还原自由度U=zeros(ndofel,kinc);V=U;A=U;fori=1:length(nodeu)node=nodeu(i);U(ndim*node-(ndim-1):ndim*node,:)=4e20;%在位移向量中标记固定边界的位置V(ndim*node-(ndim-1):ndim*node,:)=4e20;%在速度向量中标记固定边界的位置A(ndim*node-(ndim-1):ndim*node,:)=4e20;%在加速度向量中标记固定边界的位置end%还原自由度index=0;fori=1:ndofel/ndimifU(ndim*i-(ndim-1),1)==4e20index=index+1;elseU(ndim*i-(ndim-1):ndim*i,:)=uu(ndim*(i-index)-(ndim-1):ndim*(i-index),:);V(ndim*i-(ndim-1):ndim*i,:)=vv(ndim*(i-index)-(ndim-1):ndim*(i-index),:);A(ndim*i-(ndim-1):ndim*i,:)=aa(ndim*(i-index)-(ndim-1):ndim*(i-index),:);endendend注意传进newmark函数的刚度矩阵等都是已经添加过边界条件的。数值算例设计悬臂梁受动荷载算例,悬臂梁一端固定,另一端受动荷载,统计悬臂端下边界的重点数据作对比,示意图为荷载曲线为计算总时长为10s,增量步长为0.01s,采用六面体网格离散,离散后的网格图为采用C3D8单元,共计5120个六面体单元,节点总数为6273。设置两种工况工况一:采用abaqus计算全部的结果工况二:采用自编newmark程序计算统计了悬臂端下边界中点的位移曲线如下图速度曲线如下图加速度曲线如下图自编newmark程序求解运动方程,与abaqus的计算结果相对比,最终结果一致(也有可能是碰巧对准了),但是依然不影响在这过程中对newmark加深了理解。详细的数据和参数不在此罗列了,有兴趣的可以私信我获取全部数据。点击卡片关注我们来源:有限元先生

未登录
还没有评论
课程
培训
服务
行家
VIP会员 学习计划 福利任务
下载APP
联系我们
帮助与反馈