2DOF Harmonic Oscillator - Is this right? (How would I do this with ODE45?)

13 次查看（过去 30 天）
onamaewa
onamaewa 2019-6-6
回答： darova 2019-6-7
I'm not entirely sure what's going on with this result. I don't think resonant frequency increases then decreases with depth.
I've attached my equations. If you know of an easier method for getting the same thing accomplished, please let me know.
I used this source as a reference for writing my script:
https://www.uml.edu/docs/Two-Degree-Vibration-Theory_tcm18-191227.pdf
clc; clear; close all;
for DEPTH = 0:0.01:0.05 % Buried depths
radius = 0.045; % Object radius (m)
SA = pi * radius^2; % SURFACE AREA (m^2)
NAT_FREQ = 220; Q = 20; % NATURAL FREQUENCY (Hz) | QUALITY FACTOR
RHO_SOIL = 1540; % (kg/m^3)
C_SOIL = 274; att_L = 0; att_T = 0;
C_SOIL_L = C_SOIL; C1 = C_SOIL_L / (1 + att_L * 1i); % (m/s)
C_SOIL_T = 120; C2 = C_SOIL_T / (1 + att_T * 1i); % (m/s)
mm = 12;
Km = mm * (2 * pi * NAT_FREQ)^2;
Rm = (2 * pi * NAT_FREQ * mm) / Q;
% Assigns values based on depths.
if DEPTH == 0
ms = 0; Ks2 = 0; Ks1 = 0;
else
ms = (1/3) * (RHO_SOIL * (DEPTH * (pi * radius^2)))/SA;
Ks2 = real(RHO_SOIL * C1^2 * SA / DEPTH);
Ks1 = real(RHO_SOIL * C2^2 * SA / DEPTH); % (Pa/m)
end
t = 1; P = 1; F = P * SA;
MASS = [ms 0; 0 mm];
SPRING = [(Ks1 + Ks2) -Ks2; -Ks2 (Km + Ks2)];
% THIS IS WHERE THE ERROR MOST LIKELY OCCURS.
for FREQ = 1:1000
OMEGA = 2 * pi * FREQ; OMEGA = transpose(OMEGA);
INPUT = F * exp(1i * OMEGA * t); FORCE = [INPUT; 0];
if (DEPTH == 0)
Rs2 = 0; Rs1 = 0;
else
Rs2 = ((-1 * imag(Ks1) / OMEGA))/SA;
Rs1 = -imag(Ks2) / OMEGA;
end
DAMPING = [(Rs1 + Rs2) -Rs2; -Rs2 (Rm + Rs2)];
Z = [(-OMEGA^2 * MASS(1,1) + 1i * OMEGA * DAMPING(1,1) + SPRING(1,1)) (-OMEGA^2 * MASS(1,2) + 1i * OMEGA * DAMPING(1,2) + SPRING(1,2));
(-OMEGA^2 * MASS(1,2) + 1i * OMEGA * DAMPING(1,2) + SPRING(1,2)) (-OMEGA^2 * MASS(2,2) + 1i * OMEGA * DAMPING(2,2) + SPRING(2,2))];
H = inv(Z);
xs(FREQ) = H(1,1) * FORCE(1) + H(1,2) * FORCE(2);
if DEPTH == 0
xm(FREQ) = H(2,2) * FORCE(1);
else
xm(FREQ) = H(2,1) * FORCE(1) + H(2,2) * FORCE(2);
end
vm(FREQ) = (FORCE(1) / SA)/Z(2,2);
vs(FREQ) = (FORCE(1) / SA)/Z(1,1);
end
xs = transpose(xs); vs = transpose(vs);
xm = transpose(xm); vm = transpose(vm);
if DEPTH == 0
plot(1:1000, abs(vm))
else
plot(1:1000, abs(vm))
end
xlabel('Frequency (Hz)'); ylabel('Velocity (m/s)');
legend('0m', '0.01m', '0.02m', '0.03m', '0.04m', '0.05m')
hold on
end

请先登录，再回答此问题。

回答（1 个）

darova
darova 2019-6-7
What about ode45?
% x(1) == xs
% x(2) == dxs
% x(3) == xm
% x(4) == dxm
% fun = @(t,x) [ Vs; d2xs; Vm; d2xm ];
fun = @(t,x) [x(2); (P0-Rs1*x(2)-Rs2*(x(2)-x(4))-Ks1*x(1)-Ks2*(x(1)-x(3)) )/ms; ...
x(4); ( 0- Rm*x(4)+Rs2*(x(2)-x(4)) -Km*x(3)+Ks2*(x(1)-x(3)) )/mm ];
[t,x] = ode45(fun,[0 5], [0 0 0 0]);
xs = x(:,1);
plot(t,xs)

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by