隐式Newmark动力积分算法及matlab源码

1月前浏览2009

本文摘要（由AI生成）：

这篇文章主要介绍了显式动力积分中心差分法的基础算法和实现过程，并提供了相关的 MATLAB 代码。文章首先回顾了之前提供的显式动力积分中心差分法基础算法的 MATLAB 实现，然后介绍了时域动力学有限元分析中常用的计算方法，包括模态叠加法和直接积分法，并详细讲解了直接积分法中的隐式积分法和显式积分法的特点和优缺点。文章还介绍了 Newmark-beta 法的具体积分公式和计算过程，并通过 MATLAB 代码进行了演示。最后，文章通过对比不同步长下显式中心差分法和 Newmark-beta 法的计算结果，展示了显式积分法的稳定时间增量步与网格密度和质量的关系。

在之前的推文显式动力学中心差分法matlab源码中，本公众号提供了显式动力积分中心差分法基础算法的matlab实现。

正如推文所讲，在时域动力学有限元分析中，常见的计算方法包括模态叠加法和直接积分法，而直接积分法又分为隐式积分法和显式积分法。

其中，隐式积分法在计算的过程中需要迭代，常见的方式是Newton-Raphson迭代。隐式积分法常见的算法包括newmark-beta法，HHT-alpha法，wilson-theta法等。该方法的优点在于其一定程度上为无条件稳定，但是在复杂模型上可能会不收敛。

显式积分法通常以中心差分法为基础，在计算过程中不需要迭代，因此不存在不收敛的情况，但是其是有条件稳定的，即必须增量步较小时才能稳定。在实际中，显式积分法的稳定时间增量步主要依赖于单元特征长度和波速，而波速又与弹性模量和密度相关，在实践中显式积分的计算效率与网格密度和质量相关度较大。在知名的通用有限元abaqus中，其采用显式积分时会在.sta文件中显示具体哪些单元对应的稳定时间增量步较小。

隐式积分算法Newmark-beta法和HHT法得以广泛的应用。Newmark-beta法的具体积分公式如下：

当

取不同的值时，其可以得到各种不同的积分方法。例如，

时，该积分格式即为之前提到的中心差分法。

时该格式为线性加速度法。

在商软中，当我们选择Newmark-beta法，其通常指的是

时的积分方法，此时的积分方法也称为平均加速度法。

以下是Newmark-beta的计算过程：

本文给出上述计算过程的matlab代码：

clcclearn=2;k=zeros(n);m=zeros(n);Cdamp=zeros(n);k=[6 -2;-2 4];m=[2 0;0 1];disp0=zeros(n,1);velocity0=zeros(n,1);acceler0=[0;0];dt=0.28; alpha=0.25;delta=0.5;a0=1/(alpha*dt^2);a1=delta/(alpha*dt);a2=1/(alpha*dt);a3=1/(2*alpha)-1;a4=delta/alpha-1;a5=dt/2*(delta/alpha-2);a6=dt*(1 - delta);a7=delta*dt;kbar=k+a0*m+a1*Cdamp;disp_0_dt=disp0-dt*velocity0+a3*acceler0;Rt_pdt=[0;10];ndtime=12;u=zeros(ndtime+1,n);t=zeros(ndtime+1,1);t(1)=0;u(1,:)=disp0;for i=1:ndtimeRt_pdt_bar=Rt_pdt+m*(a0*disp0+a2*velocity0+a3*acceler0)+Cdamp*(a1*disp0+a4*velocity0+a5*acceler0);disp_t_pdt=kbar\Rt_pdt_bar;u(i+1,:)=disp_t_pdt;acceler_t_pdt=a0*(disp_t_pdt - disp0) -a2*velocity0-a3*acceler0;velocity_t=velocity0+a6*acceler0+a7*acceler_t_pdt;t(i+1)=i*dt;disp0=disp_t_pdt;velocity0=velocity_t;acceler0=acceler_t_pdt;endfigure(1),plot(t,u(:,1),"r-"),grid on;xlabel("time");ylabel("disp_1");hold on;%legend("Newmark");figure(2),plot(t,u(:,2),"B-"),grid on;xlabel("time");ylabel("disp_2");hold on;%legend("Newmark");