Modal Analysis - Modal Parameter Extraction

The Mode Indicator Functions are usually used on \(n\times p\) FRF matrix where \(n\) is a relatively large number of measurement DOFs and \(p\) is the number of excitation DOFs, typically 3 or 4.

In these methods, the frequency dependent FRF matrix is subjected to a singular value decomposition analysis which thus yields a small number (3 or 4) of singular values, these also being frequency dependent.

These methods are used to determine the number of modes present in a given frequency range, to identify repeated natural frequencies and to pre process the FRF data prior to modal analysis.

From the documentation of the modal software:

The MIF consist of the singular values of the Frequency response function matrix. The number of MIFs equals the number of excitations. By the powerful singular value decomposition, the real signal space is separated from the noise space. Therefore, the MIFs exhibit the modes effectively. A peak in the MIFs plot usually indicate the existence of a structural mode, and two peaks at the same frequency point means the existence of two repeated modes. Moreover, the magnitude of the MIFs implies the strength of the a mode.

We compute the Complex Mode Indicator Function. The result is shown on figure 1. The exact same curve is obtained when computed using the OROS software.

MIF = zeros(size(FRFs, 2), size(FRFs, 2), size(FRFs, 3));

for i = 1:length(freqs)
  [~,S,~] = svd(FRFs(:, :, i));
  MIF(:, :, i) = S'*S;
end

We can also compute the CMIF using the FRF matrix expressed in the same global frame. We compare the two CMIF on figure 2.

They do not indicate the same resonance frequencies, especially around 110Hz.

MIF_O = zeros(size(FRFs_O, 2), size(FRFs_O, 2), size(FRFs_O, 3));
for i = 1:length(freqs)
  [~,S,~] = svd(FRFs_O(:, :, i));
  MIF_O(:, :, i) = S'*S;
end

An alternative is the Composite Response Function \(HH(\omega)\) defined as the sum of all the measured FRF:

Instead, we choose here to use the sum of the norms of the measured FRFs:

The result is shown on figure 3.

figure;
hold on;
plot(freqs, squeeze(sum(sum(abs(FRFs)))), '-k');
plot(freqs, squeeze(sum(sum(abs(FRFs_O)))), '--k');
hold off;
xlabel('Frequency [Hz]'); ylabel('Amplitude');
xlim([1, 200]);

Modal identification are done within the Modal software of OROS. The manual for the software is available here.

Several modal parameter extraction methods are available. We choose to use the “broad band” method as it permits to identify the modal parameters using all the FRF curves at the same time. It takes into account the fact the the properties of all the individual curves are related by being from the same structure: all FRF plots on a given structure should indicate the same values for the natural frequencies and damping factor of each mode.

Such method also have the advantage of producing a unique and consistent model as direct output.

In order to apply this method, we select the frequency range of interest and we give an estimate of how many modes are present.

Then, it shows a stabilization charts, such as the one shown on figure 4, where we have to manually select which modes to take into account in the modal model.

We can then run the modal analysis, and the software will identify the modal damping and mode shapes at the selected frequency modes.

The obtained modal parameters are:

• resonance frequencies in Hertz
• modal damping ratio in percentage
• (complex) modes shapes for each measured DoF
• modal A and modal B which are parameters important for further normalization

They are all exported in a text file named modes.asc. Its first 20 lines as shown below.

Created by N-Modal
Estimator: bbfd
16-Jul-19 09:35:42


Mode 1
freq    = 11.86509Hz
damp    = 12.20318%
modal A =   4.13559e+003+6.22828e+003i
modal B =   4.98475e+005-2.49344e+005i
Mode matrix of local coordinate [DOF: Re IM]
     1X+:  -1.08546e-001   3.17769e-002
     1Y+:   2.32355e-001   1.71539e-003
     1Z+:  -1.98455e-002   5.79676e-003
     2X+:  -8.59861e-002   2.84442e-002
     2Y+:   2.40594e-001   0.00000e+000
     2Z+:  -7.65144e-003   4.22952e-003
     3Y+:   2.14274e-001   7.53040e-003
     3X-:   1.22352e-001  -2.90328e-002
     3Z+:  -2.12103e-002  -3.81368e-003

We split this big modes.asc file into sub text files using bash. The obtained files are described one table 1.

sed '/^\s*[0-9]*[XYZ][+-]:/!d' mat/modes_2019_07_16.asc > mat/mode_shapes.txt
sed '/freq/!d'                 mat/modes_2019_07_16.asc | sed 's/.* = \(.*\)Hz/\1/' > mat/mode_freqs.txt
sed '/damp/!d'                 mat/modes_2019_07_16.asc | sed 's/.* = \(.*\)\%/\1/' > mat/mode_damps.txt
sed '/modal A/!d'              mat/modes_2019_07_16.asc | sed 's/.* =\s\+\([-0-9.e]\++[0-9]\+\)\([-+0-9.e]\+\)i/\1 \2/' > mat/mode_modal_a.txt
sed '/modal B/!d'              mat/modes_2019_07_16.asc | sed 's/.* =\s\+\([-0-9.e]\++[0-9]\+\)\([-+0-9.e]\+\)i/\1 \2/' > mat/mode_modal_b.txt

Then we import the obtained .txt files on Matlab using readtable function.

shapes_m = readtable('mat/mode_shapes.txt',              'ReadVariableNames', false);  % [Sign / Real / Imag]
freqs_m  = table2array(readtable('mat/mode_freqs.txt',   'ReadVariableNames', false)); % in [Hz]
damps_m  = table2array(readtable('mat/mode_damps.txt',   'ReadVariableNames', false)); % in [%]
modal_a  = table2array(readtable('mat/mode_modal_a.txt', 'ReadVariableNames', false)); % [Real / Imag]
modal_b  = table2array(readtable('mat/mode_modal_b.txt', 'ReadVariableNames', false)); % [Real / Imag]

We guess the number of modes identified from the length of the imported data.

acc_n = 23; % Number of accelerometers
dir_n = 3; % Number of directions
dirs = 'XYZ';

mod_n = size(shapes_m,1)/acc_n/dir_n; % Number of modes

As the mode shapes are split into 3 parts (direction plus sign, real part and imaginary part), we aggregate them into one array of complex numbers.

T_sign = table2array(shapes_m(:, 1));
T_real = table2array(shapes_m(:, 2));
T_imag = table2array(shapes_m(:, 3));

mode_shapes = zeros(mod_n, dir_n, acc_n);

for mod_i = 1:mod_n
  for acc_i = 1:acc_n
    % Get the correct section of the signs
    T = T_sign(acc_n*dir_n*(mod_i-1)+1:acc_n*dir_n*mod_i);
    for dir_i = 1:dir_n
      % Get the line corresponding to the sensor
      i = find(contains(T, sprintf('%i%s',acc_i, dirs(dir_i))), 1, 'first')+acc_n*dir_n*(mod_i-1);
      mode_shapes(mod_i, dir_i, acc_i) = str2num([T_sign{i}(end-1), '1'])*complex(T_real(i),T_imag(i));
    end
  end
end

The obtained mode frequencies and damping are shown below.

It seems that the modal analysis software makes the assumption of viscous damping for the model with which it tries to fit the FRF measurements.

If we note \(N\) the number of modes identified, then there are \(2N\) eigenvalues and eigenvectors given by the software:

for \(r = 1, \dots, N\) where \(\omega_r\) is the natural frequency and \(\xi_r\) is the critical damping ratio for that mode.

We would like to arrange the obtained modal parameters into two modal matrices: \[ \Lambda = \begin{bmatrix} s_1 & & 0 \\ & \ddots & \\ 0 & & s_N \end{bmatrix}_{N \times N}; \quad \Psi = \begin{bmatrix} & & \\ \{\psi_1\} & \dots & \{\psi_N\} \\ & & \end{bmatrix}_{M \times N} \] \[ \{\psi_i\} = \begin{Bmatrix} \psi_{i, 1_x} & \psi_{i, 1_y} & \psi_{i, 1_z} & \psi_{i, 2_x} & \dots & \psi_{i, 23_z} \end{Bmatrix}^T \]

\(M\) is the number of DoF: here it is \(23 \times 3 = 69\). \(N\) is the number of mode

eigen_val_M = diag(2*pi*freqs_m.*(-damps_m/100 + j*sqrt(1 - (damps_m/100).^2)));
eigen_vec_M = reshape(mode_shapes, [mod_n, acc_n*dir_n]).';

Each eigen vector is normalized: \(\| \{\psi_i\} \|_2 = 1\)

However, the eigen values and eigen vectors appears as complex conjugates: \[ s_r, s_r^*, \{\psi\}_r, \{\psi\}_r^*, \quad r = 1, N \]

In the end, they are \(2N\) eigen values. We then build two extended eigen matrices as follow: \[ \mathcal{S} = \begin{bmatrix} s_1 & & & & & \\ & \ddots & & & 0 & \\ & & s_N & & & \\ & & & s_1^* & & \\ & 0 & & & \ddots & \\ & & & & & s_N^* \end{bmatrix}_{2N \times 2N}; \quad \Phi = \begin{bmatrix} & & & & &\\ \{\psi_1\} & \dots & \{\psi_N\} & \{\psi_1^*\} & \dots & \{\psi_N^*\} \\ & & & & & \end{bmatrix}_{M \times 2N} \]

eigen_val_ext_M = blkdiag(eigen_val_M, conj(eigen_val_M));
eigen_vec_ext_M = [eigen_vec_M, conj(eigen_vec_M)];

We also build the Modal A and Modal B matrices:

With \(a_i\) is the “Modal A” parameter linked to mode i.

modal_a_M = diag(complex(modal_a(:, 1), modal_a(:, 2)));
modal_b_M = diag(complex(modal_b(:, 1), modal_b(:, 2)));

modal_a_ext_M = blkdiag(modal_a_M, conj(modal_a_M));
modal_b_ext_M = blkdiag(modal_b_M, conj(modal_b_M));

“Modal A” and “modal B” are linked through the following formula: \[ B = - A \Lambda \]

From the modal model, we want to synthesize the Frequency Response Functions that has been used to build the modal model.

Let’s recall that:

• \(M\) is the number of measured DOFs (\(3 \times n_\text{acc}\))
• \(N\) is the number of modes identified

We then have that the FRF matrix \([H_{\text{syn}}]\) can be synthesize using the following formula:

An alternative formulation is: \[ H_{pq}(s_i) = \sum_{r=1}^N \frac{A_{pqr}}{s_i - \lambda_r} + \frac{A_{pqr}^*}{s_i - \lambda_r^*} \] with:

• \(A_{pqr} = \frac{\psi_{pr}\psi_{qr}}{M_{A_r}}\), \(M_{A_r}\) is called “Modal A”
• \(\psi_{pr}\): scaled modal coefficient for output DOF \(p\), mode \(r\)
• \(\lambda_r\): complex modal frequency

From the modal software documentation:

Modal A Scaling constant for a complex mode. It has the same properties as modal mass for normal modes (undamped or proportionally damped cases). Assuming

• \(\psi_{pr}\) = Modal coefficient for measured degree of freedom p and mode r
• \(\psi_{qr}\) = Modal coefficient for measured degree of freedom q and mode r
• \(A_{pqr}\) = Residue for measured degree of freedom p, measured degree of q and mode r
• \(M_{Ar}\) = Modal A of mode r

Then \[ A_{pqr} = \frac{\psi_{pr}\psi_{qr}}{M_{Ar}} \]

Modal B Scaling constant for a complex mode. It has the same properties as modal stiffness for normal modes (undamped or proportionally damped cases). Assuming

• \(M_{Ar}\) = Modal A of mode r
• \(\lambda_r\) = System pole of mode r

Then \[ M_{Br} = - \lambda_r M_{Ar} \]

Hsyn = zeros(acc_n*dir_n, acc_n*dir_n, length(freqs));

for i = 1:length(freqs)
  Hsyn(:, :, i) = eigen_vec_ext_M*((j*2*pi*freqs(i)).^2*inv(modal_a_ext_M)/(diag(j*2*pi*freqs(i) - diag(eigen_val_ext_M))))*eigen_vec_ext_M.';
end

Because the synthesize frequency response functions are representing the displacement response in \([m/N]\), we multiply each element of the FRF matrix by \((j \omega)^2\) in order to obtain the acceleration response in \([m/s^2/N]\).

for i = 1:size(Hsyn, 1)
  Hsyn(i, :, :) = squeeze(Hsyn(i, :, :)).*(j*2*pi*freqs).^2;
end

If we note \(\xi_i = \{ \xi_{i,x}\ \xi_{i,y}\ \xi_{i,z} \}^T\) the mode shape corresponding to the accelerometer \(i\).

From the figure above, we can write:

With

\(\xi\) and \(\Omega\) represent to mode shapes in translation and rotation of the solid expressed in the frame \(\{O\}\).

We can rearrange the equations in a matrix form, and then we obtain the mode shape of the solid in the wanted frame \(\{O\}\):

This inversion is equivalent to a mean square problem.

The obtained mode shapes matrix that gives the mode shapes of each solid bodies with respect to the fixed frame \(\{O\}\), mode_shapes_O, is an \(n \times p \times q\) with:

• \(n\) is the number of identified modes
• \(p\) is the number of DOFs for each solid body (6)
• \(q\) is the number of solid bodies

mode_shapes_O = zeros(mod_n, 6, length(solid_names));

for mod_i = 1:mod_n
  for solid_i = 1:length(solid_names)
    acc_i = solids.(solid_names{solid_i});

    Y = mode_shapes(mod_i, :, acc_i);
    Y = Y(:);

    A = zeros(3*length(acc_i), 6);
    for i = 1:length(acc_i)
      A(3*(i-1)+1:3*i, 1:3) = eye(3);

      A(3*(i-1)+1:3*i, 4:6) = [ 0              acc_pos(i, 3) -acc_pos(i, 2) ;
                               -acc_pos(i, 3)  0              acc_pos(i, 1) ;
                                acc_pos(i, 2) -acc_pos(i, 1)  0];
    end

    mode_shapes_O(mod_i, :, solid_i) = A\Y;
  end
end

We then rearrange the eigen vectors in another way: \[ \Psi_O = \begin{bmatrix} & & \\ \{\psi_1\} & \dots & \{\psi_N\} \\ & & \end{bmatrix}_{M \times N} \] with \[ \{\psi\}_r = \begin{Bmatrix} \psi_{1, x} & \psi_{1, y} & \psi_{1, z} & \psi_{1, \theta_x} & \psi_{1, \theta_y} & \psi_{1, \theta_z} & \psi_{2, x} & \dots & \psi_{6, \theta_z} \end{Bmatrix}^T \]

With \(M = 6 \times n_\text{solid}\) is the new number of DOFs and \(N\) is the number of modes.

eigen_vec_O = reshape(mode_shapes_O, [mod_n, 6*length(solid_names)]).';

eigen_vec_ext_O = [eigen_vec_O, conj(eigen_vec_O)];

First, we load the position of the Center of Mass of each solid body with respect to the point of interest.

model_com = table2array(readtable('mat/model_solidworks_com.txt', 'ReadVariableNames', false));

Then, we reshape the data.

model_com = reshape(model_com , [3, 6]);

Now we convert the mode shapes expressed in the DOF of the accelerometers to the DoF of each solid body (translations and rotations) with respect to their own CoM.

mode_shapes_CoM = zeros(mod_n, 6, length(solid_names));

for mod_i = 1:mod_n
  for solid_i = 1:length(solid_names)
    acc_i = solids.(solid_names{solid_i});

    Y = mode_shapes(mod_i, :, acc_i);
    Y = Y(:);

    A = zeros(3*length(acc_i), 6);
    for i = 1:length(acc_i)
      acc_pos_com = acc_pos(i, :).' - model_com(:, i);

      A(3*(i-1)+1:3*i, 1:3) = eye(3);
      A(3*(i-1)+1:3*i, 4:6) = [ 0               acc_pos_com(3) -acc_pos_com(2) ;
                               -acc_pos_com(3)  0               acc_pos_com(1) ;
                                acc_pos_com(2) -acc_pos_com(1)  0];
    end

    mode_shapes_CoM(mod_i, :, solid_i) = A\Y;
  end
end

We then rearrange the eigen vectors in another way: \[ \Psi_{\text{CoM}} = \begin{bmatrix} & & \\ \{\psi_1\} & \dots & \{\psi_N\} \\ & & \end{bmatrix}_{M \times N} \] with \[ \{\psi\}_r = \begin{Bmatrix} \psi_{1, x} & \psi_{1, y} & \psi_{1, z} & \psi_{1, \theta_x} & \psi_{1, \theta_y} & \psi_{1, \theta_z} & \psi_{2, x} & \dots & \psi_{6, \theta_z} \end{Bmatrix}^T \]

With \(M = 6 \times n_\text{solids}\) is the new number of DOFs and \(N\) is the number of modes.

Each eigen vector is normalized.

eigen_vec_CoM = reshape(mode_shapes_CoM, [mod_n, 6*length(solid_names)]).';
% eigen_vec_CoM = eigen_vec_CoM./vecnorm(eigen_vec_CoM);
eigen_vec_ext_CoM = [eigen_vec_CoM, conj(eigen_vec_CoM)];

sed '/Mass = /!d' mat/1_granite_bot.txt | sed 's/^\s*\t*Mass = //' | sed 's/\([0-9.]*\) kilograms\r/\1/' >  mat/model_solidworks_mass.txt
sed '/Mass = /!d' mat/2_granite_top.txt | sed 's/^\s*\t*Mass = //' | sed 's/\([0-9.]*\) kilograms\r/\1/' >> mat/model_solidworks_mass.txt
sed '/Mass = /!d' mat/3_ty.txt          | sed 's/^\s*\t*Mass = //' | sed 's/\([0-9.]*\) kilograms\r/\1/' >> mat/model_solidworks_mass.txt
sed '/Mass = /!d' mat/4_ry.txt          | sed 's/^\s*\t*Mass = //' | sed 's/\([0-9.]*\) kilograms\r/\1/' >> mat/model_solidworks_mass.txt
sed '/Mass = /!d' mat/5_rz.txt          | sed 's/^\s*\t*Mass = //' | sed 's/\([0-9.]*\) kilograms\r/\1/' >> mat/model_solidworks_mass.txt
sed '/Mass = /!d' mat/6_hexa.txt        | sed 's/^\s*\t*Mass = //' | sed 's/\([0-9.]*\) kilograms\r/\1/' >> mat/model_solidworks_mass.txt

sed '/I[xyz][xyz]/!d' mat/1_granite_bot.txt | sed 's/^\s*\t*//' | sed 's/I[xyz][xyz] = //g' | sed $'s/\r//' >  mat/model_solidworks_inertia.txt
sed '/I[xyz][xyz]/!d' mat/2_granite_top.txt | sed 's/^\s*\t*//' | sed 's/I[xyz][xyz] = //g' | sed $'s/\r//' >> mat/model_solidworks_inertia.txt
sed '/I[xyz][xyz]/!d' mat/3_ty.txt          | sed 's/^\s*\t*//' | sed 's/I[xyz][xyz] = //g' | sed $'s/\r//' >> mat/model_solidworks_inertia.txt
sed '/I[xyz][xyz]/!d' mat/4_ry.txt          | sed 's/^\s*\t*//' | sed 's/I[xyz][xyz] = //g' | sed $'s/\r//' >> mat/model_solidworks_inertia.txt
sed '/I[xyz][xyz]/!d' mat/5_rz.txt          | sed 's/^\s*\t*//' | sed 's/I[xyz][xyz] = //g' | sed $'s/\r//' >> mat/model_solidworks_inertia.txt
sed '/I[xyz][xyz]/!d' mat/6_hexa.txt        | sed 's/^\s*\t*//' | sed 's/I[xyz][xyz] = //g' | sed $'s/\r//' >> mat/model_solidworks_inertia.txt

model_mass = table2array(readtable('mat/model_solidworks_mass.txt', 'ReadVariableNames', false));
model_inertia = table2array(readtable('mat/model_solidworks_inertia.txt', 'ReadVariableNames', false));

\[ M = \begin{bmatrix} M_1 & & & & & \\ & M_2 & & & & \\ & & M_3 & & & \\ & & & M_4 & & \\ & & & & M_5 & \\ & & & & & M_6 \\ \end{bmatrix}, \text{ with } M_i = \begin{bmatrix} m_i & & & & & \\ & m_i & & & & \\ & & m_i & & & \\ & & & I_{i,xx} & I_{i,yx} & I_{i,zx} \\ & & & I_{i,xy} & I_{i,yy} & I_{i,zy} \\ & & & I_{i,xz} & I_{i,yz} & I_{i,zz} \\ \end{bmatrix} \]

M = zeros(6*6, 6*6);

for i = 1:6
  M((i-1)*6+1:i*6, (i-1)*6+1:i*6) = blkdiag(model_mass(i)*eye(3), model_inertia((i-1)*3+1:i*3, 1:3));
end

To do so, it seems that we need to know the mass matrix \([M]\). Then: \[ \{\phi\}_r = \frac{1}{\sqrt{m_r}} \{\psi\}_r \text{ where } m_r = \{\psi\}_r^T [M] \{\psi\}_r \]

Mass-normalized eigen vectors are very important for the synthesis and spatial model extraction.

Mass Matrix could be estimated from the SolidWorks model.

We compute the modal masses that will be used for normalization.

mr = zeros(size(eigen_vector_matrix_CoM, 2), 1);

for i = 1:length(mr)
  mr(i) = real(eigen_vector_matrix_CoM(:, i).' * M * eigen_vector_matrix_CoM(:, i));
end

figure;
plot(freqs_m, mr, 'ko');
xlabel('Mode Frequency [Hz]'); ylabel('Modal mass [kg]');

And finally, we compute the mass-normalized eigen vectors.

eigen_vector_mass_CoM = zeros(size(eigen_vector_matrix_CoM));

for i = 1:size(eigen_vector_matrix_CoM, 2)
  eigen_vector_mass_CoM(:, i) = 1/sqrt(mr(i)) * eigen_vector_matrix_CoM(:, i);
end

Verification that \[ [\Phi]^T [M] [\Phi] = [I] \]

eigen_vector_matrix_CoM.'*M*eigen_vector_matrix_CoM

eigen_vector_mass_CoM*M*eigen_vector_mass_CoM.'

Other test for normalized eigen vectors

eigen_vector_mass_CoM = (eigen_vector_matrix_CoM.'*diag(diag(M))*eigen_vector_matrix_CoM)^(-0.5) * eigen_vector_matrix_CoM';
eigen_vector_mass_CoM = eigen_vector_mass_CoM.';

In general, the form of response model with which we are concerned is an FRF matrix whose order is dictated by the number of coordinates \(n\) included in the test. Also, as explained, it is normal in practice to measured and to analyze just a subset of the FRF matrix but rather to measure the full FRF matrix. Usually one column or one row with a few additional elements are measured.

Thus, if we are to construct an acceptable response model, it will be necessary to synthesize those elements which have not been directly measured. However, in principle, this need present no major problem as it is possible to compute the full FRF matrix from a modal model using:

\(\{\Phi\}\) is a mass-normalized eigen vector.

FRF_matrix_CoM = zeros(size(eigen_vector_mass_CoM, 1), size(eigen_vector_mass_CoM, 1), length(freqs));

for i = 1:length(freqs)
  FRF_matrix_CoM(:, :, i) = eigen_vector_mass_CoM*inv(eigen_value_M^2 - (2*pi*freqs(i)*eye(size(eigen_value_M, 1)))^2)*eigen_vector_mass_CoM.';
end

exc_dir = 3;
meas_mass = 6;
meas_dir = 3;

figure;
hold on;
plot(freqs, abs(squeeze(FRF_matrix_CoM((meas_mass-1)*6 + meas_dir, 6*2+exc_dir, :))));
plot(freqs, abs(squeeze(FRFs_CoM((meas_mass-1)*6 + meas_dir, exc_dir, :))));
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Amplitude');
xlim([1, 200]);

FRF_matrix = zeros(size(eigen_vector_mass, 1), size(eigen_vector_mass, 1), length(freqs));

for i = 1:length(freqs)
  FRF_matrix(:, :, i) = eigen_vector_mass*inv(eigen_value_M - (freqs(i)*eye(size(eigen_value_M, 1)))^2)*eigen_vector_mass.';
end

figure;
hold on;
plot(freqs, abs(squeeze(FRFs_O(1, 1, :))));
plot(freqs, abs(squeeze(FRF_matrix(1, 1, :))));
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Amplitude');
xlim([1, 200]);