newmark法求解运动方程（附matlab代码）
问题引入
在数值仿真中，常常涉及到静力和动力的问题。
静力问题通常忽略物体的惯性力和阻尼效应。也就是说，系统中的力只与物体的静态形变和外部载荷有关，而不考虑动态过程中速度和加速度带来的影响，因此静力问题可以用方程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:kincfprintf("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/ndimif U(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加深了理解。详细的数据和参数不在此罗列了，有兴趣的可以私信我获取全部数据。
点击卡片 关注我们
