车桥耦合振动车体响应提取
有谁懂车桥耦合吗，通过车桥耦合获得车体的加速度响应，为什么我这个代码实现不了，麻烦有懂的能教我怎么修改一下吗
%桥梁参数
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');
“Devil组”引证GPT后的撰写：
你没有定义变量“xc”，导致在后面的代码中无法计算。
插值函数中的项与导数项之间存在错误的符号。正确的表达式如下：
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 = -6xc/L^2(1-xc/L);
N2d = (1 - xc/L)(1 - 3xc/L);
N3d = 6xc/L^2(1-xc/L);
N4d = xc/L*(4*xc/L-3);
在你的循环中，你计算了100个时间步长的加速度响应，但没有对它们进行任何处理。你应该将每个时间步长的加速度响应相加并求平均值，以得到车体的平均加速度响应。
在你的代码中，没有输出加速度响应的结果。
参考chatGPT的内容和自己的思路，以下是修改后的代码：
% 桥梁参数
bridgeLength = 1;
bridgeSpan = 3;
EI = 1.65e9;
m_b = 2400;
g = 9.8;
% 车辆参数
v = 3;
k_v = 500000;
m_v = 1000;
c_v = 0;
% 计算刚度和质量矩阵
M_b = m_b * bridgeLength / 420 * [156, 22*bridgeLength, 54, -13*bridgeLength;
22*bridgeLength, 4*bridgeLength^2, 13*bridgeLength, -3*bridgeLength^2;
54, 13*bridgeLength, 156, -22*bridgeLength;
-13*bridgeLength, -3*bridgeLength^2, -22*bridgeLength, 4*bridgeLength^2];
K_b = 2 * EI / bridgeLength^3 * [6, 3*bridgeLength, -6, 3*bridgeLength;
3*bridgeLength, 2*bridgeLength^2, -3*bridgeLength, bridgeLength^2;
-6, -3*bridgeLength, 6, -3*bridgeLength;
3*bridgeLength, bridgeLength^2, -3*bridgeLength, 2*bridgeLength^2];
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) + K_b;
M1(i:i+3, i:i+3) = M1(i:i+3, i:i+3) + M_b;
end
回答不易，还请能采纳！！！
参考gpt和自己的思路，你可以尝试将下面这段代码插入到你的代码中，这样就可以获得车体的加速度响应了：
# 计算车体加速度响应
a = np.zeros_like(v)
for i in range(1, len(v)):
a[i] = (v[i] - v[i-1]) / (t[i] - t[i-1])
# 绘制加速度响应曲线
plt.figure()
plt.plot(t, a, 'r-', linewidth=2)
plt.xlabel('时间 (s)')
plt.ylabel('加速度 (m/s^2)')
plt.title('车体加速度响应曲线')
plt.grid(True)
plt.show()
这段代码会计算出车体的加速度响应，并绘制出加速度随时间变化的曲线。你可以将这段代码插入到你的代码的末尾，就可以得到车体的加速度响应了。
以下答案由GPT-3.5大模型与博主波罗歌共同编写：
车桥耦合是指车辆运行时车身与车桥之间存在一定的耦合关系，而这种耦合关系会导致车桥系统的振动，进而影响车体的加速度响应。因此，提取车体的响应需要同时考虑车辆的动力学模型和物理学模型。
以下是一份python代码实现车桥耦合振动车体响应提取的示例，供参考：
import numpy as np
from scipy.integrate import odeint
# 定义参数（参考值）
m1 = 500 # 车体质量
m2 = 100 # 车轮质量
k1 = 8000 # 车体振动刚度
k2 = 250000 # 车轮振动刚度
c1 = 80 # 车体阻尼系数
c2 = 1000 # 车轮阻尼系数
L1 = 1.5 # 车体长度
L2 = 0.5 # 车轮位置
# 定义车体动力学方程
def veh_dyn(x, t):
x1, x2, x3, x4 = x # x1:车体位置；x2:车体速度；x3:车轮位置；x4:车轮速度
dxdt = [x2,
(-c1 * x2 - k1 * (x1 - x3) + c2 * (x4 - x2) + k2 * (x3 - x1)) / m1,
x4,
(-c2 * (x4 - x2) - k2 * (x3 - x1)) / m2]
return dxdt
# 设定时间范围
t = np.linspace(0, 10, 1000)
# 定义初始状态：车体位置为0，车体、车轮速度均为0.01
x0 = [0, 0.01, L2, 0.01]
# 求解ODE方程
x = odeint(veh_dyn, x0, t)
# 提取车体加速度响应并绘图
veh_acc = np.gradient(x[:, 1], t)
plt.plot(t, veh_acc)
plt.xlabel('Time (s)')
plt.ylabel('Acceleration (m/s^2)')
plt.title('Vehicle response to bridge-coupled vibration')
plt.show()
需要注意的是，由于车桥耦合振动的响应受到许多因素的影响，实际情况可能比这份示例代码复杂得多。因此，在实际应用中，需要遵循科学的方法论和数据分析技巧，通过合理的模型和算法来提升模拟和预测的准确度。
如果我的回答解决了您的问题，请采纳！
该回答引用ChatGPT
修正如下
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;
-13*L,-3*L^2,-22*L,4*L^2,13*L,-3*L^2,22*L,-L^2;
-13*L,3*L^2,-22*L,-4*L^2,13*L,3*L^2,-22*L,4*L^2];
kb = 2 * EI / L^3 * [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];
% 组装第一部分
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 = zeros(1, 100);
B = zeros(1, 100);
C = zeros(1, 100);
for il = 1:100
K2 = zeros(9,9);
M2 = zeros(9,9);
F = zeros(9,1);
t = 0.01;
xl = v * t * il;
d = L1 / 3;
n = floor(xl / d) + 1;
xc = xl - (n-1) * d;
N1 = 1 - 3 * (xc / L)^2 + 2 * (xc / L)^3;
N2 = xc * (1 - (xc / L))^2;
% 计算车桥耦合系统的刚度和质量矩阵
K2(1:4,1:4) = K1(1:4,1:4) + kv / mv * [N1, 0, -N1, 0;
0, N2, 0, -N2;
-N1, 0, N1, 0;
0, -N2, 0, N2];
M2(1:4,1:4) = M1(1:4,1:4) + mv / 4 * eye(4);
% 计算车体的载荷向量
F(1) = -kv / mv * N1 * xc;
F(2) = -kv / mv * N2 * xc;
% 计算车桥耦合系统的加速度响应
K = K1 + K2;
M = M1 + M2;
C = cv * eye(9);
A(:,il) = [zeros(4,1);1;zeros(4,1)];
B(:,il) = [zeros(4,1);1;zeros(4,1)];
C(:,il) = inv(M) * (F - K * A(:,il) - C * B(:,il));
end
% 绘制车体加速度随时间的变化曲线
t = 0:0.01:1;
acc = zeros(1, length(t));
for i = 1:length(t)
acc(i) = A(:,1)' * C(:,i);
end
plot(t, acc);
xlabel('Time (s)');
ylabel('Acceleration (m/s^2)');
title('Vehicle Acceleration Response');
- 这有个类似的问题, 你可以参考下: https://ask.csdn.net/questions/7705684