2DOF Harmonic Oscillator - Is this right? (How would I do this with ODE45?)
5 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
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:
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
0 Kommentare
Antworten (1)
darova
am 7 Jun. 2019
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)
0 Kommentare
Siehe auch
Kategorien
Mehr zu Ordinary Differential Equations finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!