2DOF System - What is wrong with the implemented differential equations here?

조회 수: 1 (최근 30일)
onamaewa
onamaewa 2019년 6월 17일
댓글: Bjorn Gustavsson 2019년 6월 17일
For the following system, you should see frequency and velocity decrease as depth is increased. However, my script is showing the opposite trend. I changed my Ks2 to have a negative height as a quick fix, but I don't know if that makes my results mathematically sound as far as the differential equations are concerned. Do you see an issue here?
Screen Shot 2019-06-17 at 11.12.26 AM.png
RHOs = 1600; % Density (kg/m^3)
RADIUS = 0.1525; % Radius (m)
HEIGHT = 8; % Depth (in)
NAT_FREQ = 188; % Natural Frequency (Hz)
% Mass Parameters
Mm = 22.34;
Km = 3.11 * 10^7;
Rm = 2730;
Cs = 265; % Wave Speed (m/s)
ATTENUATION = 0.06; % Attenuation
% Calculations-------------------------------------------------------------
AREA = pi * RADIUS^2; % Surface Area (m^2)
HEIGHT = HEIGHT * 0.0254; % Depth (m)
CsP = Cs / (1 + ATTENUATION * 1i); % Longitudinal Wave Speed (m/s)
CsS = 0.5 * CsP; % Transverse Wave Speed
K = CsP^2 * RHOs; % Bulk Modulus (N/m^2) - Using Longitudinal Wave Speed
G = CsS^2 * RHOs; % Shear Modulus (N/m^2) - Using Shear Wave Speed
FREQ = 1:1000; FREQ = transpose(FREQ); % Frequency Data (1 - 1000 Hz)
LAMBDA = CsS ./ FREQ; % Wavelength from 1 - 1000 Hz (m)
COMP = K.^-1; % Compressibility (m^2/N)
% Textbook Formulas
Ms = RHOs * AREA * HEIGHT; % Mass (kg)
Ks1 = ((8 * pi) ./ LAMBDA) * G * RADIUS; % Shear Spring 1 (N/m^2)
Ks2 = (1 / COMP) * (AREA / -HEIGHT); % Spring (m^3/N) [changed height to negative sign to make correct trend]
Ks2 = real(Ks2); Rs2 = imag(Ks2); % Parameters
Ks1 = real(Ks1); Rs1 = imag(Ks1); % Parameters
PRESSURE = 1; FORCE = PRESSURE * AREA; % Converts Pressure (Pa) to Force (N)
NAT_OMEGA = 2 * pi * NAT_FREQ; % Converts Natural Frequency to Angular Frequency (rad/s)
vs = []; vm = []; % Initializes arrays for the for-loop.
for j = 1:length(FREQ)
OMEGA = 2 * pi * j;
if (HEIGHT == 0) % NO DEPTH
X = -1i * OMEGA * FORCE / (-Km + 1i * OMEGA * Rm + OMEGA^2 * Mm);
vs = [vs X];
else % WITH DEPTH
Rs1 = 0; Ks1 = 0;
A11 = Ks1 + Ks2 - 1i * OMEGA * (Rs1 + Rs2) - OMEGA^2 * Ms;
A12 = Ks2 + 1i * OMEGA * Rs2;
A21 = A12;
A22 = Km + Ks2 - 1i * OMEGA * (Rs2 + Rm) - OMEGA^2 * Mm;
A = [A11 A12; A21 A22];
X = A \ [-1i * OMEGA * FORCE; 0];
vs = [vs X(1)]; % Velocity
vm = [vm X(2)]; % Velocity
end
end
plot(1:1000, abs(vs)); xlim([0 500]);
title('Without Shear Parameters');
xlabel('Frequency (Hz)'); ylabel('Velocity (m/s/1Pa)');
% legend('0in', '2in', '4in', '8in', '16in')
hold on
  댓글 수: 3
onamaewa
onamaewa 2019년 6월 17일
편집: onamaewa 2019년 6월 17일
I'm using a laplace transform and then applying cramer's rule to solve the equations. A & X do those steps.
Bjorn Gustavsson
Bjorn Gustavsson 2019년 6월 17일
Why bother with Laplace transform, why not solve it analytically? It is a a system of linear ODEs...

댓글을 달려면 로그인하십시오.

답변 (0개)

카테고리

Help CenterFile Exchange에서 Programming에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by