车桥耦合振动车体响应提取

Question

有谁懂车桥耦合吗，通过车桥耦合获得车体的加速度响应，为什么我这个代码实现不了，麻烦有懂的能教我怎么修改一下吗

Answer 1

%桥梁参数
L=1;
L1=3;
EI=1.65e9;
m1=2400;
g=9.8;
%车辆参数
v=3;
kv=500000;
mv=1000;
cv=0;
%%%%%%计算
mb=m1*L/420*[156,22*L,54,-13*L;22*L,4*L^2,13*L,-3*L^2;
    54,13*L,156,-22*L;-13*L,-3*L^2,-22*L,4*L^2];
kb=2*EI/L/L/L*[6,3*L,-6,3*L;3*L,2*L^2,-3*L,L^2;
    -6,-3*L,6,-3*L;3*L,L^2,-3*L,2*L^2];
%%%N=[1-3*(xc/L)^2+2*(xc/L)^3;xc*(1-(xc/L))^2;
 %%%  3*(xc/L)^2-2*(xc/L)^3;xc^2/L*(xc/L-1)];
K1=zeros(9,9);
    M1=zeros(9,9);
for i=1:2:5
   K1(i:i+3,i:i+3)=K1(i:i+3,i:i+3)+kb;
   M1(i:i+3,i:i+3)=M1(i:i+3,i:i+3)+mb;
end
%%%%组装第二部分
A=cell(1,100);
B=cell(1,100);
C=cell(1,100);
for i1=1:100
 K2=zeros(9,9);
 M2=zeros(9,9);
 F=zeros(9,1);
t=0.01;
x1=v*t*i1; %%%接触点位置
d=L1/3; %%%单元长度
n=ceil(x1/d); %%接触点所在单元
xc=x1-(n-1)*L;
%%%插值函数
N1=1-3*(xc/L)^2+2*(xc/L)^3;
N2=xc*(1-(xc/L))^2;
N3=3*(xc/L)^2-2*(xc/L)^3;
N4=(xc^2/L)*(xc/L-1);
N1d=-6*xc/(L.^2)+6*((xc)^2/(L)^3);
N2d=(1-(xc/L))^2-2*xc/L;
N3d=6*xc/(L^2)-6*((xc)^2/(L)^3);
N4d=(xc/L)*(3*(xc/L)-2);
%%叠加
kc=[v*cv*N1*N1d+kv*N1*N1,v*cv*N1*N2d+kv*N1*N2,v*cv*N1*N3d+kv*N1*N3,v*cv*N1*N4d+kv*N1*N4;
    v*cv*N2*N1d+kv*N2*N1,v*cv*N2*N2d+kv*N2*N2,v*cv*N2*N3d+kv*N2*N3,v*cv*N2*N4d+kv*N2*N4;
    v*cv*N3*N1d+kv*N3*N1,v*cv*N3*N2d+kv*N3*N2,v*cv*N3*N3d+kv*N3*N3,v*cv*N3*N4d+kv*N3*N4;
    v*cv*N4*N1d+kv*N4*N1,v*cv*N4*N2d+kv*N4*N2,v*cv*N4*N3d+kv*N4*N3,v*cv*N4*N4d+kv*N4*N4];

K2(2*n-1:2*n+2,2*n-1:2*n+2)=K2(2*n-1:2*n+2,2*n-1:2*n+2)+kc;
k1=[-v*cv*N1d-kv*N1,-v*cv*N2d-kv*N2,-v*cv*N3d-kv*N3,-v*cv*N4d-kv*N4];
K2(9,2*n-1:2*n+2)=K2(9,2*n-1:2*n+2)+k1;
k2=[-kv*N1;-kv*N2;-kv*N3;-kv*N4];
K2(2*n-1:2*n+2,9)=K2(2*n-1:2*n+2,9)+k2;
k3=[kv];
K2(9,9)=K2(9,9)+k3;
K=K1+K2;
mc=[mv];
M2(9,9)=M2(9,9)+mc;
M=M1+M2;
fc=[-mv*g*N1;-mv*g*N2;-mv*g*N3;-mv*g*N4];
F(2*n-1:2*n+2,1)=F(2*n-1:2*n+2,1)+fc;
A(1,i1)={K};
B(1,i1)={M};
C(1,i1)={F};
end
%%%%%%%%%%%%%%Newmark
gama=0.5;
beta=0.25;
dt=0.01;
a0=1/beta/dt/dt;
a1=gama/beta/dt;
a2=1/beta/dt;
a3=1/2/beta-1;
a4=gama/beta-1;
a5=dt/2*(gama/beta-2);
a6=dt*(1-gama);
a7=dt*gama;
d=zeros(9,101);
v1=zeros(9,101);
a=zeros(9,101);
tt=zeros(9,101);
for i2=1:100
    tt(1,i2+1)=i2*dt;
    if(i2*dt<=1)
    C{1,i2};
    end
ke=A{1,i2}+a0*B{1,i2};
fe=C{1,i2}+B{1,i2}*(a0*d(:,i2)+a2*v1(:,i2)+a3*a(:,i2));
  d(:,i2+1)=inv(ke)*fe;
  a(:,i2+1)=a0*(d(:,i2+1)-d(:,i2))-a2*v1(:,i2)-a3*a(:,i2);
  v1(:,i2+1)=v1(:,i2)+a6*a(:,i2)+a7*a(:,i2+1);
end
grid on 
hold on 
box on 
plot(tt(1,:),a(9,:));
xlabel('t/s'); 
ylabel('w/m');