Azzera filtri
Azzera filtri

How to find Quadratic Damping?

29 visualizzazioni (ultimi 30 giorni)
Jake
Jake il 17 Mar 2024
Modificato: Mathieu NOE il 17 Mag 2024 alle 12:51
Hi,
I have a time series data set from a motion decay test (sample attached).
load('sampleData.mat');
% whos
figure
plot(x, y);
xlabel('time [s]'); ylabel('amplitude'), grid minor
I want to find the quadratic equivalent damping using this data set, in the form of,
damping = A*amplitude + B*amplitude*|amplitude|
where, A and B are the coefficients that needs to be found. I tried the Logarithmic decrement approach, but I couldn't find the right peaks, plus I don't know how to use that approach to find two coefficients. Any help is very much appreciated.
  2 Commenti
Mathieu NOE
Mathieu NOE il 18 Mar 2024
Modificato: Mathieu NOE il 17 Mag 2024 alle 12:51
this is a starter - I will probably have a bit more time at the end of the week for the quadratic term
the code for linear damping (log decrement) is :
load('sampleData.mat');
ind = 1:20000; % I restrict to the first half of data to avoid the noise at the end of the record
[Ypk,Xpk,Wpk,Ppk] = findpeaks(y(ind),'MinPeakHeight',1,'MinPeakDistance',200);
plot(x,y,x(Xpk),Ypk,'dr')
label('time [s]'); ylabel('amplitude')
n_peaks=numel(Xpk);
%Vector length of interest
x_new = x(Xpk);
y_new = Ypk;
%---Calculate Logarithmic Decrement, undamped natural frequency, damped natural frequency, damping ratio
Log_Dec = zeros(length(n_peaks));
for nn = 1:n_peaks-1 %
Log_Dec(nn)= log(y_new(nn)/y_new(nn+1)); % computed with n = 1 periods apart
end
Mean_dec = mean(Log_Dec); %Calculate Average Logarithmic Decrement
%Damping
damp_ratio_logdec = 1/sqrt(1+((2*pi/(Mean_dec))^2)) %Assesses Damping Constant
Jake
Jake il 18 Mar 2024
Hi @Mathieu NOE, I see that by restricting the data, you have avoided the noise, as you also mentioned in a comment, which was the issue I had with logarithmic decrement approach. Thank you for this suggestion. I will study and also wait to see if the quadratic solution is possible. Thank you again!

Accedi per commentare.

Risposta accettata

Mathieu NOE
Mathieu NOE il 29 Mar 2024
hello again
with some delay, here now the code modified and extended
of course we already had the (equivalent) linear damping result , using the well know log decrement
dzeta = 0.0305 (normalized damping coefficient)
pe = 0.0959 (linear damping coefficient)
the second portion of the code allows to identify p1 and p2 for the quadratic damping coefficients
we get the following results
Nb that p1 is very close to pe from the linear damping above
p1 = 0.1101
p2 = 0.0116
code :
load('sampleData.mat');
y = detrend(y,'linear'); % remove some y offset (to be seen at the end of the record)
%% LINEAR DAMPING identification
ind = x<50; % I restrict to the first half of data to avoid the noise at the end of the record
[Yp,Xp,Wp,Pp] = findpeaks(y(ind),'NPeaks',30,'MinPeakHeight',1,'MinPeakDistance',200); % all positive and negative peaks
Xpp = x(Xp);
Ypp = y(Xp);
n_peaks=numel(Xpp);
figure(1),plot(x,y,Xpp,Ypp,'dr')
xlabel('time [s]'); ylabel('amplitude')
%---Calculate Logarithmic Decrement, undamped natural frequency, damped natural frequency, damping ratio
Log_Dec = zeros(length(n_peaks));
for nn = 1:n_peaks-1 %
Log_Dec(nn)= log(Ypp(nn)/Ypp(nn+1)); % computed with n = 1 periods apart
end
Mean_dec = mean(Log_Dec); %Calculate Average Logarithmic Decrement
wn = 2*pi/mean(diff(Xpp)); %Calculate Average Pulsation
%Damping
dzeta = 1/sqrt(1+((2*pi/(Mean_dec))^2)) % normalized damping coefficient
pe = dzeta*2*wn %linear damping coefficient
%% LINEAR + QUADRATIC DAMPING identification - Froude method
% ¨η + ( p1 + p2*|˙η|)*˙η + ω²*η = 0
for ck = 1:n_peaks-1
wn = 2*pi/(Xpp(ck+1) - Xpp(ck)); % computation of wn (from the peaks distance)
amplitude_mean = (Ypp(ck) + Ypp(ck+1))/2;
delta_amplitude = Ypp(ck) - Ypp(ck+1);
alpha(ck) = 8*wn*amplitude_mean/(3*pi);
pe(ck) = 2*wn*delta_amplitude/(pi*amplitude_mean);
end
% Plot pe vs alpha , first order polynomial fit :
% Fit a polynomial p of degree "degree" to the (x,y) data:
degree = 1;
p = polyfit(alpha,pe,degree);
% Evaluate the fitted polynomial p and plot:
pe_f = polyval(p,alpha);
eqn = poly_equation(flip(p)); % polynomial equation (string)
Rsquared = my_Rsquared_coeff(pe(:),pe_f(:)); % correlation coefficient
figure(2);plot(alpha,pe,'*',alpha,pe_f,'-')
xlabel('alpha'); ylabel('pe')
legend('data',eqn)
title(['Data fit , R² = ' num2str(Rsquared)]);
%linear damping coefficient pe = p1 + p2*alpha
p1 = p(2)
p2 = p(1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Rsquared = my_Rsquared_coeff(data,data_fit)
% R2 correlation coefficient computation
% The total sum of squares
sum_of_squares = sum((data-mean(data)).^2);
% The sum of squares of residuals, also called the residual sum of squares:
sum_of_squares_of_residuals = sum((data-data_fit).^2);
% definition of the coefficient of correlation is
Rsquared = 1 - sum_of_squares_of_residuals/sum_of_squares;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function eqn = poly_equation(a_hat)
eqn = " y = "+a_hat(1);
for i = 2:(length(a_hat))
if sign(a_hat(i))>0
str = " + ";
else
str = " ";
end
if i == 2
% eqn = eqn+" + "+a_hat(i)+"*x";
eqn = eqn+str+a_hat(i)+"*x";
else
% eqn = eqn+" + "+a_hat(i)+"*x^"+(i-1)+" ";
eqn = eqn+str+a_hat(i)+"*x^"+(i-1)+" ";
end
end
eqn = eqn+" ";
end

Più risposte (0)

Categorie

Scopri di più su Mathematics in Help Center e File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by