论坛 › 一般力学和振动理论 › 多自由度瞬态动力学强迫响应计算==Newmark&Runge-Kutta

mxlzhenzhu 发表于 2015-10-14 11:45

多自由度瞬态动力学强迫响应计算==Newmark&Runge-Kutta

本帖最后由 mxlzhenzhu 于 2015-10-14 11:47 编辑

%% 多自由度系统瞬态响应分析,New-mark Beta，适用于非比例阻尼，非线性刚度，非线性阻尼;
%% Inputs:
% K Stiff matrix
% C Damping matrix, Structural Damping will be transformed to viscous
% damping.
% M Mass matrix
% fi_set Force excitation DOFs.
% Force Force matrix for every i_set
% R_set Output DOFs, Disp, Velo,Acc output in cell format.

%% mxl.2015-5-24
% 单位制不做规定；默认自由度序号为从1到N
% 在本程序中没有考虑非线性刚度 K=K(t)，非线性阻尼C=C(t)这类问题，后续可以添加；
% 输入默认为：x(0)=0,x'(0)=1;t 表示时间；
% 强迫位移，强迫速度和强迫加速度功能没有考虑；
% 结构阻尼输入的时候，转换为等效粘性阻尼的功能还没有添加；
% 输出请求为位移，速度和加速度，不含应力；
clear
clc


dt=0.0001;
t=';% 延迟计算时间到15秒，可以看到明显的数值阻尼
method='Newmark';% Runge-Kutta,Runge-Kutta
% M=0.2533;
% K0=10;
% C=0.592;
m1=2e2;m2=5e3;
k1=2e6;k2=1.5e6;
c1=1000;c2=2000;

M=diag();
C=[
c1+c2,-c2;
-c2,c2;
];
K0=[
k1+k2,-k2;
-k2,k2;
];


fi_set=1;
Force=5*sin(pi*t/0.6);



N=size(M,1);

initial_disp=zeros(N,1);% You may change it.
initial_velo=;% You may change it.
initial_acc=[];% You may change it.

R_set=';
% %---------------------------------------------------------------------%
% (1)Input check
% %---------------------------------------------------------------------%

if max(fi_set)>N
error('fi_set is wrongly set.');
end

if numel(fi_set)-size(Force,2)
error('Input force must keep in consistence.');
end

if isempty(dt)
error('You must specify or calculate time step.');
end

if strcmp(method,'Newmark')
if isempty(initial_disp)||isempty(initial_velo)
error('You must specify the initial states.');
end
end
if max(R_set)>N||isempty(R_set)
error('Response DOFs are badly set.')
end



%-----初始化输出-----%
Nsteps=numel(t);
disp=zeros(Nsteps,numel(R_set));% 这是用来保存R_set 位移结果的；
velo=zeros(Nsteps,numel(R_set));% 这是用来保存R_set 速度结果的；
acc=zeros(Nsteps,numel(R_set));% 这是用来保存R_set 加速度结果的；
Fi=zeros(N,1);

update_disp=zeros(N,1);% 这是用来保存运算中全局位移结果的；
update_velo=zeros(N,1);% 这是用来保存运算中全局速度结果的；
update_acc=zeros(N,1);% 这是用来保存运算中全局加速度结果的；

Inv_M=inv(M);
% %---------------------------------------------------------------------%
% (2)New-Mark Beta
% gama=0.5;0.1667<=beta<=0.25,beta=0.25 is suggested.
% %---------------------------------------------------------------------%
gama=0.5;beta=0.25;
disp(['Gama= ',num2str(gama),'Beta= ',num2str(beta)]);
a2=1/beta/dt;a0=a2/dt;a1=a2*gama;
a3=1/2/beta-1;a4=gama/beta-1;a5=dt/2*(a4-1);

% if strcmp(method,'Newmark')
for loop=1:Nsteps
if loop>1


disp(loop,:)=initial_disp(R_set,1)';% 保存请求的计算结果
velo(loop,:)=initial_velo(R_set,1)';
acc(loop,:)=initial_acc(R_set,1)';
if loop==Nsteps
break;
end
%---we may update the K in the line in nonlinear problem----%
% K=Ki+a0*M+a1*C;

Fi=zeros(N,1);
for loopf=1:numel(fi_set)
Fi(fi_set(loopf),1)=Force(loop,loopf);
end
Fi=Fi+M*(a0*initial_disp+a2*initial_velo+a3*initial_acc)+...
C*(a1*initial_disp+a4*initial_velo+a5*initial_acc);


update_disp=K\Fi;
update_velo=a1*(update_disp-initial_disp)-a4*initial_velo-a5*initial_acc;
update_acc =a0*(update_disp-initial_disp)-a2*initial_velo-a3*initial_acc;


initial_disp=update_disp;% 将本次结果保存，下一次迭代时作为初值
initial_velo=update_velo;
initial_acc =update_acc;
else
for loopf=1:numel(fi_set)
Fi(fi_set(loopf),1)=Force(loop,loopf);
end

disp(loop,:)=initial_disp(R_set,1)';
velo(loop,:)=initial_velo(R_set,1)';
if isempty(initial_acc)
initial_acc=Inv_M*(Fi-C*initial_velo-K0*initial_disp);
acc(loop,:)=initial_acc(R_set,1)';
else
acc(loop,:)=initial_acc(R_set,1)';
end


K=K0+a0*M+a1*C;% 对于线性系统，以后K都不会变


Fi=Fi+M*(a0*initial_disp+a2*initial_velo+a3*initial_acc)+...
C*(a1*initial_disp+a4*initial_velo+a5*initial_acc);

update_disp=K\Fi;
update_velo=a1*(update_disp-initial_disp)-a4*initial_velo-a5*initial_acc;
update_acc =a0*(update_disp-initial_disp)-a2*initial_velo-a3*initial_acc;

initial_disp=update_disp;% 将本次结果保存，下一次迭代时作为初值
initial_velo=update_velo;
initial_acc =update_acc;
end



end % end of for
% end % end of if
subplot(2,1,1)
plot(t,disp)
hold on
plot(t,velo(:,1),'r')
hold off
legend('disp','velo')
title('Newmark')
% axis(gca,)
% %---------------------------------------------------------------------%
% (3)Runge-Kutta,注意R-K 可以应用于一般非线性的求解，因此很厉害的；
% 请参考盛宏玉书上308页内容。
% %---------------------------------------------------------------------%
% R-K 方法要求 输入刚度，质量，阻尼，载荷，以及非线性方程的形式；
% 初始条件只需给定初始位移和初始速度，初始加速度可以由此计算得到
disp=zeros(Nsteps,numel(R_set));% 这是用来保存R_set 位移结果的；
velo=zeros(Nsteps,numel(R_set));% 这是用来保存R_set 速度结果的；
acc=zeros(Nsteps,numel(R_set));% 这是用来保存R_set 加速度结果的；
% if strcmp(method,'Runge-Kutta')

% A=[
% zeros(size(M)),eye(N);
% -Inv_M*K,-inv_M*C;
% ];
initial_disp=zeros(N,1);% You may change it.
initial_velo=;% You may change it.
initial_acc=[];% You may change it.

for loop=1:Nsteps
if loop>1

if loop==Nsteps
break;
end

Fn1=zeros(N,1);
for loopf=1:numel(fi_set)
Fn1(fi_set(loopf),1)=Force(loop,loopf);
end


Kd1=Inv_M*(-K0*initial_disp-C*initial_velo+Fn1);% 这就是一般非线性的表达式，可以另外写为一个函数单独调用
Kd2=Inv_M*(-K0*(initial_disp+0.5*dt*initial_velo)-C*(initial_velo+0.5*dt*Kd1)+0.5*(Fn1+Fn0));
Kd3=Inv_M*(-K0*(initial_disp+0.5*dt*initial_velo+0.25*dt^2*Kd1)-C*(initial_velo+0.5*dt*Kd2)+0.5*(Fn0+Fn1));
Kd4=Inv_M*(-K0*(initial_disp+dt*initial_velo+0.5*dt^2*Kd2)-C*(initial_velo+dt*Kd3)+Fn1);

Fn0=Fn1;% 保存上一次的载荷，两个时间步之间的就用平均值

update_disp=initial_disp+dt*initial_velo+(dt^2)/6*(Kd1+Kd2+Kd3); % Kdi是第二组R-K法系数
update_velo=initial_velo+dt/6*(Kd1+2*Kd2+2*Kd3+Kd4);
update_acc =Kd1;


initial_disp=update_disp;% 将本次结果保存，下一次迭代时作为初值
initial_velo=update_velo;
initial_acc=update_acc;

disp(loop,:)=initial_disp(R_set,1)';% 保存输出
velo(loop,:)=initial_velo(R_set,1)';
acc(loop,:)=initial_acc(R_set,1)';

else


disp(loop,:)=initial_disp(R_set,1)';
velo(loop,:)=initial_velo(R_set,1)';
if isempty(initial_acc)
Fn0=zeros(N,1);
for loopf=1:numel(fi_set)
Fn0(fi_set(loopf),1)=Force(loop,loopf);
end


initial_acc=Inv_M*(Fn0-C*initial_velo-K0*initial_disp);
acc(loop,:)=initial_acc(R_set,1)';
else
acc(loop,:)=initial_acc(R_set,1)';
end


end



end

% end



subplot(2,1,2)
plot(t,disp)
hold on
plot(t,velo(:,1),'r')
hold off
legend('disp','velo')
title('Runge-Kutta')
% axis(gca,)
% %---------------------------------------------------------------------%
% (4)Post-Processing
% %---------------------------------------------------------------------%


% %---------------------------------------------------------------------%
% (5)
% %---------------------------------------------------------------------%

mxlzhenzhu 发表于 2015-10-14 11:52

瞬态动力学计算有很多算法，模态瞬态法是一种；另外，就是直接积分法，其中大多数存在稳定性和数值阻尼问题，而Newmark和Runge-Kutta算法用得比较多。在此分享这个程序。


论坛上的程序也多有人讨论和分享过，可是代码可读性差，移植性差，逻辑结构可能不怎么清楚，因此发一个自己的程序，自诩容易看懂得多。

oneonly 发表于 2015-10-18 17:55

老师您好，我的二阶动力学方程为（112个）（M C K F）都是时间t的函数


按照p308页下面的那个的那种不变化为一阶方程，而是直接用龙格库塔对其求解，结合上面的程序，我的编程如下：

clc
clear

tt=0:0.1:1;       %%这里我也取过0：0.01:1
for ij=1:length(tt)
   t=tt(ij)

dt=0.001;%%步长也换过0.01,0.0001

%%%以下程序就是求时间t时刻的 M K F C （是时变的）

M=。。。

K=。。。

F=。。。

C= 。。。
%%%

Inv_M=inv(M);

if ij==1
q=zeros(112,1);                      %%初始位移
q1d=zeros(112,1);                  %%初始速度
q2d=Inv_M*(QA-K*q-C*q1d);%%初始加速度
end
if ij>1
qq=q;
qq1d=q1d;
qq2d=q2d;   


Kd1=Inv_M*(-K*qq-C*qq1d+F);
Kd2=Inv_M*(-K*(qq+0.5*dt*qq1d)-C*(qq1d+0.5*dt*Kd1)+F);    %这里的F看公式是（tn+1/2*dt)
Kd3=Inv_M*(-K*(qq+0.5*dt*qq1d+0.25*dt^2*Kd1)-C*(qq1d+0.5*dt*Kd2)+F);
Kd4=Inv_M*(-K*(qq+dt*qq1d+0.5*dt^2*Kd2)-C*(qq1d+dt*Kd3)+F）;

update_disp=qq+dt*qq1d+(dt^2)/6*(Kd1+Kd2+Kd3);
update_velo=qq1d+dt/6*(Kd1+2*Kd2+2*Kd3+Kd4);
update_acc =Kd1;

q=update_disp;      %保存结果，下一次迭代时作为初值
q1d=update_velo;
q2d=update_acc;



    q107=q(107);%q的第107行是要求的X速度
    q108=q(108);%q的第108行是要求的y速度
    q109=q(109);%q的第109行是要求的z速度
    qqq1(ij)=q107;%%保存到qqq中
    qqq2(ij)=q108;%%
    qqq3(ij)=q109;%%

end%%%%%%%%%%%%%%%       我也试着放在这里
   %q107=q(107);%q的第107行是要求的X速度
    %q108=q(108);%q的第108行是要求的y速度
    %q109=q(109);%q的第109行是要求的z速度
    %qqq1(ij)=q107;%%保存到qqq中
    %qqq2(ij)=q108;%%
    %qqq3(ij)=q109;%%
%%%%%%%%%%%%%
   
end

figure(1)
plot(tt,qqq1,'b')   
title('x位移')
xlabel('时间t/s')
ylabel('位移/m')

figure(2)
plot(tt,qqq2,'m')
title('y位移')
xlabel('时间t/s')
ylabel('位移/m')




出来的数据（区间是分了十一份）


0        -0.0346758535292791        -15365335555923.0        -7.52481294950494e+27        -3.68121477935759e+42        -1.79891687691818e+57        -8.78113256011268e+71        -4.28158412948546e+86        -2.08529717616388e+101        -1.01446242889355e+116        -4.92950756226880e+130
能劳烦看一下吗，找了好长时间了，没有找到原因。
谢谢谢谢谢谢

oneonly 发表于 2015-10-18 19:41

{:{39}:}

oneonly 发表于 2015-10-18 19:46

这是我的变系数动力学方程：                                          http://forum.vibunion.com/thread-138097-1-1.html按照书中308页下边的那种龙格库塔的方法，结合您的帖子上面的程序，不做变换化为一阶方程，而直接用龙哥-库塔对其求解。我的编程思路大致如下： clcclear %%%%%%%%%%%%%%%%%%%%%%%             Ixy=0.81*10^-8; Ixz=0.56*10^-2;Iyx=0.81*10^-8; Iyy=0.66432;       Iyz=1.166*10^-7;Izx=0.56*10^-2; Izy=1.166*10^-7; Izz=0.64; 密度，质量 转动惯量等已知是的赋值 %%%%%%%%%%%%%%%%%%%%%%%%%%%% tt=0:0.01:1; %%%%%%这里也试过0.001,0.0001等for ij=1:length(tt)   t=tt(ij)dt=0.01;%%%%%%%%%步长也试过0.1,0.0001等 。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。 M= K= F= C= 。。。。。。。。。。。。。。。。。。。。。。 。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。%%%%%句号中间是求M K F C （都是关于t的函数）Inv_M=inv(M); if ij==1q=zeros(112,1);%%%%初始位移q1d=zeros(112,1);   %%%%初始速度q2d=Inv_M*(QA-K*q-C*q1d); %%%%初始加速度 endif ij>1qq=q;qq1d=q1d;qq2d=q2d; Kd1=Inv_M*(-K*qq-C*qq1d+F);%Kd2=Inv_M*(-K*(qq+0.5*dt*qq1d)-C*(qq1d+0.5*dt*Kd1)+F);Kd3=Inv_M*(-K*(qq+0.5*dt*qq1d+0.25*dt^2*Kd1)-C*(qq1d+0.5*dt*Kd2)+F);Kd4=Inv_M*(-K*(qq+dt*qq1d+0.5*dt^2*Kd2)-C*(qq1d+dt*Kd3)+F); update_disp=qq+dt*qq1d+(dt^2)/6*(Kd1+Kd2+Kd3);%update_velo=qq1d+dt/6*(Kd1+2*Kd2+2*Kd3+Kd4);update_acc =Kd1; q=update_disp; q1d=update_velo;q2d=update_acc; q107=q(107); %%%x位移    q108=q(108); %%%y位移    q109=q(109); %%%y位移   qqq1(ij)=q107;   qqq2(ij)=q108;   qqq3(ij)=q109; end    end figure(1)plot(tt,qqq1,'b')%%%x位移    title('λòÆ')xlabel('ê±¼ät/s')ylabel('λòÆ/m') figure(2)plot(tt,qqq2,'m')%%%y位移title('Ëù¶è')xlabel('ê±¼ät/s')ylabel('λòÆ/m') figure(3)plot(tt,qqq3,'m')%%%z位移title('Ëù¶è')xlabel('ê±¼ät/s')ylabel('λòÆ/m')

mxlzhenzhu 发表于 2015-10-20 16:32

oneonly 发表于 2015-10-18 19:46
这是我的变系数动力学方程：                                          http://forum.vibunion.com/th ...

这个程序就是专门为你发的，你自己再看看吧。

查看完整版本: 多自由度瞬态动力学强迫响应计算==Newmark&Runge-Kutta