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

2 次查看（过去 30 天）

onamaewa 2019-6-17

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/467500-2dof-system-what-is-wrong-with-the-implemented-differential-equations-here

评论： 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