필터 지우기
필터 지우기

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

조회 수: 4 (최근 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:
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 CenterFile Exchange에서 Ordinary Differential Equations에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by