runge_kutta_sin_hard_rigidity_shock_response2007_6_22.m

在基础半正弦冲击作用下

%%Runge-Kutta法计算硬特性刚度半正弦冲击响应
%%定义硬特性刚度为双曲正切刚度，抗力＝2k*d1/pi*tg(pi*相对位移/2d1)
function program;
global F;
global oemiga_shock;;
%%global x_compel1; global v_compel1; global a_compel1; global delta_t; global omeiga_n; global i;
  %%   global abs_a;  global abs_v; global xx;   global vv; global aa; global zuni; global abs_x;
    %% global nzongnengliang; global nshineng; global ndongneng; global nnengliang; global nratio_nengliang; 
     %%global nratio_shineng;global nratio_dongneng;

d1=0.005;
F=10;           %半正弦冲击力幅值 

shockin_max=10;

freq_shock=50;
oemiga_shock=freq_shock*2*pi;
j=1;
shock_stepnum=1000;
abs_v(1)=0;

    for gama=(20:200)/10
   n=1; 
for cgma=(0:199)/200

    omeiga_n=oemiga_shock/gama;
    shock_time=pi/(oemiga_shock);
     delta_t=shock_time/shock_stepnum;


    
   xx(1)=0;
   vv(1)=0;
   aa(1)=0;   
   abs_a(1)=0;
   abs_x(1)=0;
   abs_v(1)=0;
   x_compel1=0;
   v_compel1=0;
   a_compel1=0;
   base_x(1)=0;
   base_v(1)=0;
   base_a(1)=0;
   abs_x1=0;
   abs_v1=0;
   abs_a1=0;
   beita=sqrt(abs(1-cgma^2));
   omeiga_d=beita*omeiga_n;     
   TT=round(((2*pi)/omeiga_d)/delta_t);        T_end=shock_stepnum+TT*3/5;
   if (T_end>20000)
       
            T_end=20000;
        end
    i=2;  
  %冲击输入===================================================================
for t=(1:shock_stepnum)*delta_t;  
    x_compel2=x_compel1+v_compel1*delta_t/2;    v_compel2=v_compel1+a_compel1*delta_t/2;
    a_compel2=-f(t+delta_t/2)-2*cgma*omeiga_n*v_compel2-(omeiga_n*omeiga_n)*2*d1/pi*tan(pi*x_compel2/(2*d1));
    

    %==================================================
    x_compel3=x_compel1+v_compel2*delta_t/2;    v_compel3=v_compel1+a_compel2*delta_t/2;
    a_compel3=-f(t+delta_t/2)-2*cgma*omeiga_n*v_compel3-(omeiga_n*omeiga_n)*2*d1/pi*tan(pi*x_compel3/(2*d1));
    
 
    
    %==================================================
    x_compel4=x_compel1+v_compel3*delta_t;    v_compel4=v_compel1+a_compel3*delta_t;
    a_compel4=-f(t+delta_t)-2*cgma*omeiga_n*v_compel4-(omeiga_n*omeiga_n)*2*d1/pi*tan(pi*x_compel4/(2*d1));
    

    
    %求出x2,v2,a2=====================================
    x_compel=x_compel1+(v_compel1+2*v_compel2+2*v_compel3+v_compel4)*delta_t/6;
    v_compel=v_compel1+(a_compel1+2*a_compel2+2*a_compel3+a_compel4)*delta_t/6;
    a_compel=-f(t+delta_t)-2*cgma*omeiga_n*v_compel-(omeiga_n*omeiga_n)*2*d1/pi*tan(pi*x_compel/(2*d1));  
    x_compel1=x_compel;    v_compel1=v_compel;    a_compel1=a_compel;   
    
    abs_a(i)=-2*cgma*omeiga_n*v_compel-(omeiga_n*omeiga_n)*2*d1/pi*tan(pi*x_compel/(2*d1));    
    %保存数据=================================================
    xx(i)=x_compel1;    vv(i)=v_compel1;    aa(i)=a_compel1;  
    i=i+1;
    
end;




 %冲击输入后==============================================================
for t=(shock_stepnum+1:T_end)*delta_t    
   
    x_compel2=x_compel1+v_compel1*delta_t/2;    v_compel2=v_compel1+a_compel1*delta_t/2;
    a_compel2=-2*cgma*omeiga_n*v_compel2-(omeiga_n*omeiga_n)*2*d1/pi*tan(pi*x_compel2/(2*d1));

    %==================================================
    x_compel3=x_compel1+v_compel2*delta_t/2;    v_compel3=v_compel1+a_compel2*delta_t/2;
    a_compel3=-2*cgma*omeiga_n*v_compel3-(omeiga_n*omeiga_n)*2*d1/pi*tan(pi*x_compel3/(2*d1));
    

    
    %==================================================
    x_compel4=x_compel1+v_compel3*delta_t;    v_compel4=v_compel1+a_compel3*delta_t;
    a_compel4=-2*cgma*omeiga_n*v_compel4-(omeiga_n*omeiga_n)*2*d1/pi*tan(pi*x_compel4/(2*d1));
    

    
    %求出x2,v2,a2=====================================
    x_compel=x_compel1+(v_compel1+2*v_compel2+2*v_compel3+v_compel4)*delta_t/6;
    v_compel=v_compel1+(a_compel1+2*a_compel2+2*a_compel3+a_compel4)*delta_t/6;
    a_compel=-2*cgma*omeiga_n*v_compel-(omeiga_n*omeiga_n)*2*d1/pi*tan(pi*x_compel/(2*d1));    
    x_compel1=x_compel;    v_compel1=v_compel;    a_compel1=a_compel;   
  
    abs_a(i)=-2*cgma*omeiga_n*v_compel-(omeiga_n*omeiga_n)*2*d1/pi*tan(pi*x_compel/(2*d1)); 
  
    %保存数据=================================================
    xx(i)=x_compel1;    vv(i)=v_compel1;    aa(i)=a_compel1;  



    i=i+1;
end;    
[xx_max,xx_max_i]=max(abs(xx));
[vv_max,vv_max_i]=max(abs(vv));
[abs_aa_max,abs_aa_max_i]=max(abs(abs_a));
[abs_vv_max,abs_vv_max_i]=max(abs(abs_v));
[abs_xx_max,abs_xx_max_i]=max(abs(abs_x));


sin_x_max(n)=xx_max;
sin_x_max1(j,n)=xx_max;
sin_x_max_i(n)=xx_max_i;
sin_x_max_i1(j,n)=xx_max_i;

sin_v_max(n)=vv_max;
sin_v_max1(j,n)=vv_max;
sin_v_max_i(n)=vv_max_i;
sin_v_max_i1(j,n)=vv_max_i;

sin_abs_a_max(n)=abs_aa_max;
sin_abs_a_max1(j,n)=abs_aa_max;
sin_abs_a_max_i(n)=abs_aa_max_i;
sin_abs_a_max_i1(j,n)=abs_aa_max_i;

base_v_max=2*F/oemiga_shock;
guiyi_sin_abs_a(n)=abs_aa_max/(base_v_max*omeiga_n);
guiyi_sin_abs_a1(j,n)=abs_aa_max/(base_v_max*omeiga_n);

guiyi_sin_x_max(n)=abs_aa_max*xx_max/(base_v_max^2);
guiyi_sin_x_max1(j,n)=abs_aa_max*xx_max/(base_v_max^2);

guiyi_sin_x_max_mine(n)=xx_max*omeiga_n/(base_v_max);
guiyi_sin_x_max_mine1(j,n)=xx_max*omeiga_n/(base_v_max);



 %j=1:i-1;
% a_max=num2str(a_max);
% x_max=num2str(x_max);
%ss=strcat('最大加速度=',a_max,'最大位移=',x_max);
% subplot(2,1,1),plot(j,abs_a),title(ss);
% subplot(2,1,2),plot(j,xx);
 xx(:)=0;
 vv(:)=0;
 aa(:)=0;
 abs_a(:)=0;

n=n+1;
end;
plot(guiyi_sin_abs_a);
hold on;
j=j+1;

end
%i=1:200;
%s=num2str(cgma);
%k=num2str(k);
%ss=strcat('加速度响应最大值随阻尼变化曲线（RUnge-Kutta法）');
%plot(i,a_max(i),'b'),title(ss),xlabel('阻尼(1/200)'),ylabel('加速度(m/s2)');
 

i=1:500;
s=num2str(cgma);
k=num2str(k);
a_max=num2str(aa_max);
x_max=num2str(xx_max);
v_max=num2str(vv_max);
ss=strcat('加速度 阻尼=',s,'刚度=',k,'a max=',a_max);
 %plot(i,abs_a(i),'b',i,xx(i),'r',i,vv(i),'g'),title(ss),xlabel('时间(1/10000)'),ylabel('加速度(m/s2)');
 subplot(3,1,1), plot(i,abs_a(i),'b'),title(ss),ylabel('加速度(m/s2)');
 ss=strcat('位移 阻尼=',s,'刚度=',k,'x max=',x_max);
 subplot(3,1,3), plot(i,xx(i),'r'),title(ss),xlabel('时间(1/10000s)'),ylabel('位移(m)');
  ss=strcat('速度 阻尼=',s,'刚度=',k,'v max=',v_max);
 subplot(3,1,2), plot(i,vv(i),'r'),title(ss),ylabel('速度(m/s)');

 y=20000;
 
 
 %=========================================================================
 
 function shock_in=f(t)
 global F;
global oemiga_shock;
%%%半正弦冲击
shock_in=F*sin(oemiga_shock*t);
%%%%%%%%%%衰减正弦冲击
%%shock_in=F*exp(-cgma_shock*oemiga_shock*t)*sin(sqrt(1-cgma_shock^2)*oemiga_shock*t);
%%%%%%%%%%%矩形冲击
%%shock_in=2*F/pi;
%%%%%%%%%%%%%%%%%%%%%%%%%%计算冲击响应%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
function shock_response(i,t,cgma)
     global x_compel1; global v_compel1; global a_compel1; global delta_t; global omeiga_n; global i;
     global abs_a;  global abs_v; global xx;   global vv; global aa; global zuni; global abs_x;
     global nzongnengliang; global nshineng; global ndongneng; global nnengliang; global nratio_nengliang;
     global nratio_shineng;global nratio_dongneng;
    %%%%%%%================================================
    x_compel2=x_compel1+v_compel1*delta_t/2;    v_compel2=v_compel1+a_compel1*delta_t/2;
    a_compel2=-f(t+delta_t/2)-2*cgma*omeiga_n*v_compel2-(omeiga_n*omeiga_n)*x_compel2;
 
    %==================================================
    x_compel3=x_compel1+v_compel2*delta_t/2;    v_compel3=v_compel1+a_compel2*delta_t/2;
    a_compel3=-f(t+delta_t/2)-2*cgma*omeiga_n*v_compel3-(omeiga_n*omeiga_n)*x_compel3;
    
    %==================================================
    x_compel4=x_compel1+v_compel3*delta_t;    v_compel4=v_compel1+a_compel3*delta_t;
    a_compel4=-f(t+delta_t)-2*cgma*omeiga_n*v_compel4-(omeiga_n*omeiga_n)*x_compel4;
    %求出x2,v2,a2=====================================
    x_compel=x_compel1+(v_compel1+2*v_compel2+2*v_compel3+v_compel4)*delta_t/6;
    v_compel=v_compel1+(a_compel1+2*a_compel2+2*a_compel3+a_compel4)*delta_t/6;
    a_compel=-f(t+delta_t)-2*cgma*omeiga_n*v_compel-(omeiga_n*omeiga_n)*x_compel;    
    x_compel1=x_compel;    v_compel1=v_compel;    a_compel1=a_compel;   
    %保存数据=================================================
    xx(i)=x_compel1;    vv(i)=v_compel1;    aa(i)=a_compel1;  
    zuni(i)=cgma;
    abs_a(i)=-(2*cgma*omeiga_n*vv(i)+omeiga_n^2*xx(i));
    abs_v(i)=abs_v(i-1)+abs_a(i)*delta_t;
    abs_x(i)=abs_x(i)+abs_v(i)*delta_t+1/2*abs_a(i)*delta_t^2;
        ndongneng(i)=1/2*(abs_v(i))^2;
    nshineng(i)=1/2*omeiga_n^2*(xx(i))^2;
    nnengliang(i)=nshineng(i)+ndongneng(i);
    nratio_shineng(i)=nshineng(i)/nnengliang(i);
     nratio_dongneng(i)=ndongneng(i)/nnengliang(i);
    nratio_nengliang(i)=nnengliang(i)/nzongnengliang;
    i=i+1;