首页/文章/ 详情

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

29天前浏览1238
   
       

问题引入

     

在数值仿真中,常常涉及到静力和动力的问题。

静力问题通常忽略物体的惯性力和阻尼效应。也就是说,系统中的力只与物体的静态形变和外部载荷有关,而不考虑动态过程中速度和加速度带来的影响,因此静力问题可以用方程AX=B表示,这种方程用线性代数里面的内容就可以求解,主要是A矩阵求逆的计算速度决定了方程的求解速度,这部分的资料很多,在此不多说,这个帖子重点讲解动力问题。

动力问题极为广泛,动力问题(Dynamic Problem)是指在力学和物理学中研究物体或结构在外力作用下随时间变化的运动、变形和应力响应。与静力问题不同,动力问题考虑了系统的惯性力和阻尼效应,这使得其分析需要解决随时间变化的运动方程。动力问题广泛应用于涉及冲击、振动、地震、波动等现象的工程和物理系统中。

动力问题的特点

考虑惯性力和加速度:在动力问题中,物体或结构的运动和变形不仅受外力和内力的影响,还受到惯性力的作用。根据牛顿第二定律,惯性力与物体的加速度成正比,因此在动力分析中必须考虑加速度的影响。

时间依赖性:动力问题的解通常是时间的函数,要求在整个时间域内求解物体或结构的运动状态。随着时间的推移,外力、位移、速度和应力等物理量会发生变化。

系统响应的类型

自由振动:不受外部力作用的振动,仅依赖于系统的初始状态和自然特性。

强迫振动:在外力作用下,系统的响应随外力变化而变化,通常伴随频率和振幅的变化。

冲击和瞬态响应:系统受到瞬时力(如冲击或爆炸)时的快速响应,通常表现为瞬态振动。

阻尼效应:动力系统中还可能存在阻尼,阻尼力通常与速度成正比,并会导致系统的运动逐渐衰减。动力问题中必须考虑阻尼效应对系统振动和响应的影响。

动力问题的基本方程:

动力问题的基本方程是牛顿第二定律或达朗贝尔原理的推广,称为运动方程(Equation of Motion)。一个简单的动力问题可以用以下形式的运动方程来描述:

   

这个方程说明了外力 F(t) 如何引起系统的位移、速度和加速度的时间响应。

这是一个二阶微分方程,求解这个方程的方法大致分为解析法数值解法

解析解当然是最精确的,但只有简单的动力系统或边界条件较为规则的问题,动力问题可以通过解析方法求解。例如:

单自由度系统的自由振动:对于简单的弹簧-质量系统,运动方程可以通过特征值问题解析求解,得到系统的固有频率和振动模态。

正弦激励下的强迫振动:当系统受到周期性外力(如正弦波)作用时,可以通过解析方法求解系统的稳态响应。

数值解法是对于复杂的几何形状、边界条件或载荷,解析解往往难以得到,因此通常采用数值方法来求解动力问题。常见的数值方法包括:

有限元法(FEM):将连续体离散为有限个单元,通过求解离散的动力方程来求解整个结构的动态响应。Abaqus 等软件广泛使用有限元方法进行复杂动力学分析。

有限差分法(FDM):将动力方程的时间和空间导数进行离散化,用数值方法逐步推进时间步长,求解各个时间点的响应。

显式和隐式时间积分法

显式方法:如中心差分法,适合处理瞬态问题和高频振动问题。

隐式方法:如 Newmark方法,适合求解稳定性要求较高的动力问题。

这里就重点介绍newmark法求解运动方程,包括固定位移边界条件的施加。

       

newmark理论及程序

     

newmark法 是一种常用的时间积分方法,用于求解结构动力学中的运动方程,尤其是在有限元分析中求解动力问题(如振动、冲击、地震响应等)。该方法由 Nathan M. 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);  % 有效刚度矩阵求逆,只进行一次
for i=2:kinc    fprintf("kinc:%d  time:%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;  for i=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;for i=1:ndofel/ndim    if U(ndim*i-(ndim-1),1)==4e20        index=index+1;    else        U(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),:);    end    endend
 

注意传进newmark函数的刚度矩阵等都是已经添加过边界条件的。

       

数值算例

     

设计悬臂梁受动荷载算例,悬臂梁一端固定,另一端受动荷载,统计悬臂端下边界的重点数据作对比,示意图为

   

荷载曲线为

   

计算总时长为10s,增量步长为0.01s,采用六面体网格离散,离散后的网格图为

   

采用C3D8单元,共计5120个六面体单元,节点总数为6273。

设置两种工况



工况一:采用abaqus计算全部的结果工况二:采用自编newmark程序计算
 

统计了悬臂端下边界中点的位移曲线如下图

   

速度曲线如下图

   

加速度曲线如下图

   

自编newmark程序求解运动方程,与abaqus的计算结果相对比,最终结果一致(也有可能是碰巧对准了),但是依然不影响在这过程中对newmark加深了理解。详细的数据和参数不在此罗列了,有兴趣的可以私信我获取全部数据

点击卡片 关注我们

     


来源:有限元先生
Abaqus振动非线性MATLAB理论爆炸
著作权归作者所有,欢迎分享,未经许可,不得转载
首次发布时间:2024-10-26
最近编辑:29天前
外太空土豆儿
硕士 我们穷极一生,究竟在追寻什么?
获赞 13粉丝 1文章 52课程 0
点赞
收藏
作者推荐

abaqus没有cae文件如何计算inp

采用abaqus进行有限元分析的时候,通常在abaqus进行建模,最终在job模块创建job提交计算。但有些工作场景中,需要采用第三方专业软件进行建模,然后导出inp文件,这时候我们只有inp文件,就无法按照正常流程打开cae进行计算,只能直接计算inp文件。对于以上工作场景,abaqus提供了多种方式运行inp文件,这里介绍三种。分别为1、把inp文件导入abaqus计算2、用bat脚本运行inp文件3、在windows的命令行窗口运行inp文件把inp文件导入abaqus计算这种方式是最推荐的,inp计算完之后可以直接查看odb文件,非常的方便。具体步骤如下,首先打开abaqus主界面打开主界面后,点击箭头1指示的部分,直奔2-job模块,弹出下面窗口点击3-jobmanager,弹出4-jobmanager窗口,然后点击5-create,弹出6-createjob窗口,然后再7-source中选择inputfile,最后点击continue,createjob窗口变成了下面的样子再点击上面图片红框框部分,弹出下面窗口8-selectinputfile这个窗口是我们的job文件所在的位置,选择9-job-1.inp,10箭头的部分会显示选择的inp文件,点击ok,窗口变成上面图片里的箭头11指示,要计算的inp文件已经被添加。点击continue,弹出下图12-editjob窗口到了这一步,如果自己的inp文件没有对应的子程序,就可以直接点击ok,如果有子程序,就在13-general中选中14添加子程序,弹出的窗口为15-selectusersubroutinefile在文件夹中选择已经准备好的子程序,这里选择的是编译之前的16-fortran源程序,也可以选择编译之后的obj文件,选中之后,17部分会显示选择的子程序,一路点击OK,回到下面的窗口,后面就可以正常提交作业检查或计算了。用bat脚本运行inp文件这种方法运用bat脚本调用abaqus主程序,首先在存放job文件的文件夹下面创建一个txt文本文件,然后将后缀修改为“bat”,这个bat文件的名字不做要求,英文就行,如下图中的“run.bat”,命名的时候会弹出修改文件后缀名不可用的提示,不用管,直接确定即可。注意:用记事本打开这个bat文件,要往里面写东西,如下图上面没有箭头指示的部分,必须要写上。如果不带子程序,则没有箭头1指示的内容,后面的2箭头指的是自己打开的核心数,依据自己的电脑核心数确定,箭头3指示的interactive表示允许用户与cmd窗口交互,abaqus主程序会向cmd窗口输出计算过程信息,第二行的箭头4表示计算完毕后cmd窗口不会自动退出,这个参数一般都要添加,方便查看cmd的报错信息,否则cmd窗口就会一闪而过。修改完毕,保存退出。运行bat文件只需要双击即可,自动调用abaqus主程序并弹出cmd窗口上图就是双击bat文件后的计算过程,cmd窗口显示了abaqus主程序计算的所有过程,特别注意下面红色框框的部分,提示我们输入任意键退出cmd窗口,如果bat文件中没有pause字样,cmd窗口在运行完之后就会直接退出,一闪而过,如果正常计算是没问题的,但是如果出错了就看不到模型的报错信息,所以一般都会加上pause。这种方法不如第一种方法方便,因为计算之后还是要打开abaqus进入odb模块查看结果。这种方法适合刚拿到inp文件之后的检查计算,如果报错了就不用打开abaqus了,直接前叉inp文件就行。在windows的命令行窗口运行inp文件这方式计算inp文件,需要首先打开windows的cmd窗口,然后用“cd”命令将工作路径切换到存放inp文件的位置,切换工作路径不敢走很简单,这里就不写了,我在另一个帖子里面详细写过,详细见cmd窗口切换工作路径。切换工作路径之后,把上文中bat文件中的内容直接复制,然后回车,会自动调用abaqus主程序进行计算,如下图这种方法与上一种基本相似,区别就是不用编辑创建一个bat文件,直接在cmd输入命令即可。这三种方法各自有各自的优缺点,适用于不同的人群,各自选择,适合计算计算场景的才是最好的。点击卡片关注我们来源:有限元先生

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