2DOF Harmonic Oscillator - Is this right? (How would I do this with ODE45?)
4 views (last 30 days)
Show older comments
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.mathworks.com/matlabcentral/answers/uploaded_files/223397/image.jpeg)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/223398/image.jpeg)
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 Comments
Answers (1)
darova
on 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 Comments
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!