matlab求解移動力作用下簡支梁變形問題？
瀏覽：1744 回答：2
請教各位老師同學：一個簡支梁上作用了一個有速度的常量力，想要用MATLAB分析梁在不同速度下（不同工況，不是變速）的位移、速度、撓度等動力響應，還有內力等，是應該建立原始的動力學方程，輸入質量、阻尼之類的矩陣求解方程，還是輸入自己手動推導好的比如位移響應表達式直接計算呢？逐步積分法能用得到嗎？如圖所示的編程思路該如何實現呢？請大佬幫幫忙
第一次接觸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;