How to find Quadratic Equivalent Damping?

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.

[Mathieu NOE]

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')
xlabel('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 

damp_ratio_logdec = 0.0216

[Jake]

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!

Answer 1

[Mathieu NOE]

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

dzeta = 0.0305

pe =  dzeta*2*wn  %linear damping coefficient 

pe = 0.0959

%% 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)

p1 = 0.1101

p2 = p(1)

p2 = 0.0116

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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
        

[Mathieu NOE]

FYI, this is the publication I used