|
|
马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。
您需要 登录 才可以下载或查看，没有账号？我要加入
x
本帖最后由 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=[0:dt:10]';% 延迟计算时间到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([m1,m2]);
- 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=[0;1];% You may change it.
- initial_acc=[];% You may change it.
- R_set=[1]';
- % %---------------------------------------------------------------------%
- % (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,[0,max(t),-2e-5, 4e-5])
- % %---------------------------------------------------------------------%
- % (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=[0;1];% 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,[0,max(t),-2e-5, 4e-5])
- % %---------------------------------------------------------------------%
- % (4)Post-Processing
- % %---------------------------------------------------------------------%
- % %---------------------------------------------------------------------%
- % (5)
- % %---------------------------------------------------------------------%
|
|
发表于 2015-10-14 11:45
提升卡
变色卡
显身卡
楼主
