显式动力学中心差分法matlab源码

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

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

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

以下是一个基础版本的中心差分的计算过程：

上述过程中需要对有效质量矩阵进行三角化，计算效率并不高。在实际中，往往可以通过对质量矩阵和阻尼矩阵C进行一些处理，使得左侧的系数矩阵为对角阵，从而使得计算过程更快。常见的方式是对于质量矩阵采用集中质量矩阵代替一致质量矩阵，此时质量矩阵为对角阵。但由于有效质量矩阵表达式中还包括阻尼矩阵，如果仅仅质量矩阵对角，有效质量矩阵并非对角阵，因此对阻尼矩阵也需要处理。

在显式积分中，常见的阻尼形式为瑞丽阻尼，即将阻尼矩阵表达为质量矩阵和刚度矩阵与对应系数的乘积之和。这种方式简化了阻尼矩阵的计算，但遗憾的是仍然不能形成对角阻尼矩阵。这时可以采用的常见方法有两种，一种是瑞丽阻尼的刚度矩阵的系数取0，阻尼矩阵仅为质量矩阵与对应的系数的乘积，此时如果质量矩阵为对角阵，则阻尼矩阵也为对角阵。另一种方式是修改原振动方程表达式将阻尼矩阵与速度乘积项取为前一次时间增量的阻尼矩阵与速度乘积，此时该项可以移至方程右侧，从而左侧仅有集中质量矩阵与加速度乘积项。

上述的具体实现细节，本文不作具体讨论，本文仅给出一个上面基本的中心差分计算的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;10];

dt=0.28;

a0=1/dt^2;a1=1/(2*dt);a2=2*a0;a3=1/a2;

disp_0_dt=disp0-dt*velocity0+a3*acceler0;

mbar=a0*m+a1*Cdamp;Rt0=[0;10];ndtime=12;u=zeros(ndtime+1,n);u(1,:)=disp0;dispt=disp0;disp_t_dt=disp_0_dt;t=zeros(ndtime+1,1);t(1)=0;for i=1:ndtimeRt=Rt0-(k-a2*m)*dispt-(a0*m-a1*Cdamp)*disp_t_dt;disp_t_pdt=mbar\Rt;u(i+1,:)=disp_t_pdt;

velocity_t=a1*(- disp_t_dt+disp_t_pdt);acceler_t=a0*(disp_t_dt-2*dispt+ disp_t_pdt);

disp_t_dt=dispt;dispt=disp_t_pdt;t(i+1)=i*dt;endfigure(1),plot(t,u(:,1)),grid on;hold on;xlabel("time");ylabel("disp_1");legend("Explicit")figure(2),plot(t,u(:,2)),grid on;hold on;xlabel("time");ylabel("disp_2");legend("Explicit")

最终的位移如图所示：

以上，即是本文的全部内容，感谢您的阅读！