matlab求解移动力作用下简支梁变形问题？

%% 移动荷载模型

clear

clc

%% 基本参数定义

% Newmark-beta参数

dt=0.001;        % 积分步长

gama=1/2;        

beta=1/4;

% 车辆参数

P0=800;          %集中力荷载

V=10;            % 行车速度

%主梁参数

m=2738.5;  % 主梁单位长度质量

L=22;               % 主梁跨度

E=3.96e10;  % 弹性模量

I=158581e-6;  % 截面惯性矩

c=0.02;

% 计算参数e

T=L/V;             % 运行时间

n=5;               % 振型阶数

times=round(T/dt); % 积分步数

g=9.8;             % 重力加速度

%% 计算过程参数

% 积分常数

alpha0=1/beta/dt^2;

alpha1=gama/beta/dt;

alpha2=1/beta/dt;

alpha3=1/2/beta-1;

alpha4=gama/beta-1;

alpha5=dt/2*(gama/beta-2);

alpha6=dt*(1-gama);

alpha7=gama*dt;

%% 系统动力系数矩阵

K=zeros(n,n);           %刚度矩阵

M=zeros(n,n);           %质量矩阵

C=zeros(n,n);           %阻尼矩阵

F=zeros(n,1);           %荷载向量

Q=zeros(n,times);       %振型矩阵

q=zeros(n,1);           %振型

w=zeros(n,1);           %频率

kxi=zeros(n,1);         %阻尼比

pF=2*P0/m/L;

%% 位移、速度、加速度向量及结果保存

DP=zeros(n,1);          % 上一时间步位移向量

VP=zeros(n,1);          % 上一时间步速度向量

AP=zeros(n,1);          % 上一时间步加速度向量

DC=zeros(n,1);          % 当前时间步位移向量

VC=zeros(n,1);          % 当前时间步速度向量

AC=zeros(n,1);          % 当前时间步加速度向量

SD=zeros(n,times);      % 位移结果保存

SV=zeros(n,times);      % 速度结果保存

SA=zeros(n,times);      % 加速度结果保存

u=zeros(1,times);

v=zeros(1,times);

av=zeros(1,times);

FC=zeros(n,times);

%% 动力响应计算

% 广义位移、速度和加速度

t=dt:dt:T;

for i=1:times

    

for j=1:n

   w(j)=j^2*pi^2*sqrt(E*I/m/L^4);  %自振频率向量

   kxi(j)=0;                       %阻尼比向量

   q(j,i)=sin(j*pi*V*i*dt/L);        %振型向量

   F(j)=pF*q(j,i);                   %广义力向量

end

   FC(:,i)=F;

                  %储存振型结果

% 组装质量矩阵

for a=1:n

  M(a,a)=1;

  C(a,a)=2*kxi(a)*w(a);

  K(a,a)=w(a)^2;

end

    equK=K+alpha0*M+alpha1*C;                      % 等效刚度

    equF=F+M*(alpha0*DP+alpha2*VP+alpha3*AP)+C*(alpha1*DP+alpha4*VP+alpha5*AP);% 等效荷载

    DC=equK\equF;                                  % 当前时间步位移响应

    AC=alpha0*(DC-DP)-alpha2*VP-alpha3*AP;         % 当前时间步加速度响应

    VC=VP+alpha6*AP+alpha7*AC;                     % 当前时间步的速度响应    

    DP=DC;                                         % 将当前时间步的位移响应作为上一时间步的位移响应

    VP=VC;                                         % 将当前时间步的速度响应作为上一时间步的速度响应

    AP=AC;                                         % 将当前时间步的加速度响应作为上一时间步的加速度响应  

    % 结果保存

    SD(:,i)=DC;                                    

    SV(:,i)=VC;

    SA(:,i)=AC;

    % 振型叠加

    zSD=0;zSV=0;zSA=0;

    for a=1:n    

%     zSD1=SD(a,i)*sin(a*pi*V*i*0.001/L);

%     zSD=zSD+zSD1;

%     zSA1=SA(a,i)*sin(a*pi*V*i*0.001/L);

%     zSA=zSA+zSA1; 

%     zSV1=SV(a,i)*sin(a*pi*V*i*0.001/L);

%     zSV=zSV+zSV1;

     zSD1=SD(a,i)*sin(a*pi/2);

     zSD=zSD+zSD1;

     zSA1=SA(a,i)*sin(a*pi/2);

     zSA=zSA+zSA1;

     zSV1=SV(a,i)*sin(a*pi/2);

     zSV=zSV+zSV1;

    end

     u(:,i)=-zSD;

     v(:,i)=-zSV;

     av(:,i)=-zSA;

end

y=u;

plot(t,1e3*y)

title('位移随时间变化曲线');

xlabel('时间/s');

ylabel('跨中位置位移/mm');

% y=v;subplot(132)

% plot(t,y)

% title('速度随时间变化曲线');

% xlabel('t(s)');

% ylabel('v(m/s)')      %速度随时间变化曲线

SD_jifen=1e3*y;