MATLAB Answers

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

22 views (last 30 days)
onamaewa
onamaewa on 17 Jun 2019
Commented: Bjorn Gustavsson on 17 Jun 2019
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 Comments

onamaewa
onamaewa on 17 Jun 2019
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 on 17 Jun 2019
Why bother with Laplace transform, why not solve it analytically? It is a a system of linear ODEs...

Sign in to comment.

Answers (0)

Community Treasure Hunt

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

Start Hunting!

Translated by