MATLAB Examples

Contents

Study case: Clamped-free beam

The dynamic response of a 100 m high clamped-free steel beam is studied. Simulated time series are used, where the first three eigen-modes have been taken into account. More precisely, the acceleration response to white noise excitation is used here. The simulated data are generated with the central difference method [1], and the eigen-modes/eigen-frequencies are generated using [2]. A modal damping ratio of 0.005 for every modes is used. These parameters are considered as the "target modal parameters".

The purpose of this example is to use the AFDD to retrieve the eigen frequencies, the mode shapes, and the modal damping ratio from Ambient vibrations analysis. This is more challenging than the classical response to impulse load, but also more common for large civil engineering structures.

The most challenging task is the estimation of the modal damping ratio from ambient vibrations data (higher sensitivity to measurement noise). Prior to a systematic analysis, the parameters used in the AFDD method have to be "tuned", i.e. manually tested, so that we extract the proper modal parameters. This operation is rather common in Structural Health Monitoring.

We use here two approaches:

  1. Manual selection of the eigen-frequencies using the Matlab code introduced by [3], where a friently-user interface is proposed.
  2. Automatic peak-picking of the eigen-frequencies by using the Matlab code developed by [4], which worked much better than the function findpeaks.

[1] https://se.mathworks.com/matlabcentral/fileexchange/55557-modal-parameters-identification-from-ambient-vibrations--sdof-

[2] https://se.mathworks.com/matlabcentral/fileexchange/52075-eigen-value-calculation-of-a-continuous-beam--transverse-vibrations-

[3] https://se.mathworks.com/matlabcentral/fileexchange/50988-frequency-domain-decomposition--fdd-

[4] https://se.mathworks.com/matlabcentral/fileexchange/42927-pickpeaks-v-select-display-

clearvars;close all;clc;
load('beamData.mat')
fn = wn/(2*pi);
Nmodes = numel(fn);
fs = 1/median(diff(t));
rng(1)

Manual procedure

This procedure has the advantage to be easier if one has only one
sample. The peak-picking method is straightforward, and lead to accurate
results. However, if multiple sample have to be analysed, the manual
procedure becomes cumbersome. For more details, see [3]
tic
[phi_FDD,fn_FDD,zeta] = AFDD(Az,t,Nmodes,'PickingMethod','manual');
toc
% we plot the mode shapes
figure
for ii=1:size(phi_FDD),
    subplot(2,2,ii)
    hold on;box on;
    h1 = plot(linspace(0,1,size(phi_FDD,2)),phi_FDD(ii,:),'ro','linewidth',1.5);
    h2 = plot(linspace(0,1,size(phi,2)),-phi(ii,:),'k-','linewidth',1.5);
    xlabel('$y$ (a.u.)','interpreter','latex')
    ylabel(['$\phi_',num2str(ii),'$'],'interpreter','latex')
    if ii==1,
        legend('Measured','Target','location','SouthWest')
    end
end
% The theoretical and measured eigen-frequencies agrees well !
disp('left: target eigen frequencies. Right: Measured eigen-frequencies')
disp([fn(:),fn_FDD(1:Nmodes)'])
disp('left: target damping. Right: Measured damping')
disp([5e-3*ones(Nmodes,1),zeta(:),])
Peak selection procedure
a: Draw rectangles around peaks while holding left click
b: Press "Space" key to continue the peak selection
c: Press "any other key" if you have selected a peak by mistake and want to ignore it
Elapsed time is 14.391932 seconds.
left: target eigen frequencies. Right: Measured eigen-frequencies
    0.1671    0.1683
    1.0472    1.0483
    2.9323    2.9356
    5.7460    5.8362

left: target damping. Right: Measured damping
    0.0050    0.0059
    0.0050    0.0042
    0.0050    0.0049
    0.0050    0.0053

Automated procedure 1: minimalist example

Minimalist example with automated procedure for those who don't want to
read too much
[phi_FDD,fn_FDD,zeta] = AFDD(Az,t,Nmodes);

% plot the mode shapes
figure
for ii=1:size(phi_FDD),
    subplot(2,2,ii)
    hold on;box on;
    h1 = plot(linspace(0,1,size(phi_FDD,2)),phi_FDD(ii,:),'gd','linewidth',1.5);
    h2 = plot(linspace(0,1,size(phi,2)),-phi(ii,:),'k-','linewidth',1.5);
    xlabel('$y$ (a.u.)','interpreter','latex')
    ylabel(['$\phi_',num2str(ii),'$'],'interpreter','latex')
    if ii==1,
        legend('Measured','Target','location','SouthWest')
    end
end

% Comparison between the measured and target eigen freq. and mode shapes
disp('left: target eigen frequencies. Right: Measured eigen-frequencies')
disp([fn(:),fn_FDD(1:Nmodes)'])
disp('left: target damping. Right: Measured damping')
disp([5e-3*ones(Nmodes,1),zeta(:),])
left: target eigen frequencies. Right: Measured eigen-frequencies
    0.1671    0.1683
    1.0472    1.0483
    2.9323    2.9356
    5.7460    5.8362

left: target damping. Right: Measured damping
    0.0050    0.0093
    0.0050    0.0039
    0.0050    0.0053
    0.0050    0.0054

Automated procedure 2: 2-step analysis

The automated procedure presented here is based on 2 steps, but can be reduced to a single step (cf. above) if the signal to noise ratio is low enough:

  1. identification of the eigen-frequencies
  2. Identification of the mode shapes and the damping ratio

An automated selection of the peak is easier to be carried out by smoothing the spectral density matrix by using a relatively low value for M. However, using a low value for M leads to an overestimation of the modal damping ratio, unless long records are used. In the present example, the length the data is 20 minutes, which is too short to carry out the automated procedure in one step. Using a low value for M is also more time consuming. Therefore, I use a limited number of value of Az only.

% First step: determination of the eigen-frequencies
[~,fn_FDD,zeta] = AFDD(Az(1:5:end,:),t,Nmodes,'dataPlot',1);
% we show that the estimated zeta is ca. 10 x larger for the first 2 modes
% than expected.
% the theoritical and measured eigen frequencies agrees however well!
disp('left: target eigen frequencies. Right: Measured eigen-frequencies')
disp([fn(:),fn_FDD(1:Nmodes)'])

% Second step: determination of the modal damping ratio
% We use a high value for M and prescribed eigen-frequencies
% We use the option 'dataPlot' to plot intermediate figures, to illustrate
% the method, and to check the accuracy of the results.
[phi_FDD,fn_FDD,zeta] = AFDD(Az,t,Nmodes,'fn',fn_FDD,'M',8192);
% Plot the mode shapes
figure
for ii=1:size(phi_FDD),
    subplot(2,2,ii)
    hold on;box on;
    h1 = plot(linspace(0,1,size(phi_FDD,2)),phi_FDD(ii,:),'csq','linewidth',1.5);
    h2 = plot(linspace(0,1,size(phi,2)),-phi(ii,:),'k-','linewidth',1.5);
    xlabel('$y$ (a.u.)','interpreter','latex')
    ylabel(['$\phi_',num2str(ii),'$'],'interpreter','latex')
    if ii==1,
        legend('Measured','Target','location','SouthWest')
    end
end
disp('left: target damping. Right: Measured damping')
disp([5e-3*ones(Nmodes,1),zeta(:),])
left: target eigen frequencies. Right: Measured eigen-frequencies
    0.1671    0.1683
    1.0472    1.0483
    2.9323    2.9356
    5.7460    5.8362

left: target damping. Right: Measured damping
    0.0050    0.0050
    0.0050    0.0051
    0.0050    0.0053
    0.0050    0.0045

Case of user-defined boundaries for the selected peaks.

The boundaries for the selected peaks (lines 209 in the main function AFDD) may not be adapted if the eigen-frequencies range from low to high frequencies. For this reason, it is possible to manually give the upper boundaries (UB) and the lower boundaries(LB) as shown below for the first 4 eigen frequencies of the beam studied:

% lower boundary for the first four modes (Default: LB = 0.9*fn)
LB = [0.15,0.9,2.8,5.5];
% upper boundary boundary for the first four modes (Default: UB = 0.9*fn)
UB = [0.18,1.15,3.1,6.1];

% Visualization of the boundaries
clf;close all;
figure
plot(1:4,fn,'ko-',1:4,LB,'r-',1:4,UB,'b-')
ylabel('$f_n$ (Hz)','interpreter','latex')
legend('Measured eigen frequency','user-defined lower boundary','user-defined upper boundary','location','best')
xlabel('Mode number')
set(gca,'xtick',[1,2,3,4])

% Calculation of the modal parameters with user-defined UBs and LBs
[phi_FDD,fn_FDD,zeta] = AFDD(Az,t,Nmodes,'M',8192,'UB',UB,'LB',LB);

figure
for ii=1:size(phi_FDD),
    subplot(2,2,ii)
    hold on;box on;
    h1 = plot(linspace(0,1,size(phi_FDD,2)),phi_FDD(ii,:),'msq','linewidth',1.5);
    h2 = plot(linspace(0,1,size(phi,2)),-phi(ii,:),'k-','linewidth',1.5);
    xlabel('$y$ (a.u.)','interpreter','latex')
    ylabel(['$\phi_',num2str(ii),'$'],'interpreter','latex')
    if ii==1,
        legend('Measured','Target','location','SouthWest')
    end
end
disp('left: target damping. Right: Measured damping')
disp([5e-3*ones(Nmodes,1),zeta(:),])
left: target damping. Right: Measured damping
    0.0050    0.0057
    0.0050    0.0044
    0.0050    0.0056
    0.0050    0.0049

Conclusion

The function AFDD was here illustrated with four examples: One manual
procedure, and three automatic ones. We have seen that the modal damping
ratio is different for each case. It comes form the estimation of the
autocorrelation function in the time domain that requires that I use
random phases, since we lost the information about the phase angle of
each signal, at each frequency step during the calculation of the psd.
Alternative approaches may be more efficient, and are welcomed !