2DOF System - What is wrong with the implemented differential equations here?
조회 수: 1 (최근 30일)
이전 댓글 표시
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?
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
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 Center 및 File Exchange에서 Assembly에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!