Finding a 2D SSM for a 3D finite element Plate

This is an example of how to reconstruct a slow 2D SSM of a mechanical system using all phase space variables measurements. In this example, we consider the perforated cover plate modelled in COMSOL Multiphysics. Details of this geometry can be found in: https://doi.org/10.1016/j.jsv.2016.12.037

generate model from COMSOL
Dynamical system setup
Linear Modal analysis and SSM setup
SSM setup
compute SSM coefficients and reduced dynamics
post-processing via SSM-based ROMs

generate model from COMSOL

clear all
close all
run ../../../install.m
[model, M, C, K, Null, Nullf, ud, outdof] = build_model();

Construction of M,C,K of the finite element model in COMSOL. We use the function provided in COMSOL to extract the internal stiffness force.

[fint,dfint] = get_fint(K,C,model,Null, Nullf, ud);

n = length(M);
disp(['Number of degrees of freedom = ' num2str(n)])
disp(['Phase space dimensionality = ' num2str(2*n)])

Number of degrees of freedom = 7405
Phase space dimensionality = 14810

Dynamical system setup

We consider the forced system

$\mathbf{M}\ddot{\mathbf{x}}+\mathbf{C}\dot{\mathbf{x}}+\mathbf{K}\mathbf{x}+\mathbf{f}(\mathbf{x},\dot{\mathbf{x}})=\epsilon\mathbf{f}^{ext}(\mathbf{\Omega}t),$

which can be written in the first-order form as

$\mathbf{B}\dot{\mathbf{z}}=\mathbf{Az}+\mathbf{F}(\mathbf{z})+\epsilon\mathbf{F}^{ext}(\mathbf{\phi}),\\\dot{\mathbf{\phi}} =\mathbf{\Omega}$

where

$\mathbf{z}=\left[\begin{array}{c}\mathbf{x}\\\dot{\mathbf{x}}\end{array}\right],\quad\mathbf{A}=\left[\begin{array}{cc}-\mathbf{K} & \mathbf{0}\\\mathbf{0} & \mathbf{M}\end{array}\right],\mathbf{B}=\left[\begin{array}{cc}\mathbf{C} & \mathbf{M}\\\mathbf{M} & \mathbf{0}\end{array}\right],\quad\quad\mathbf{F}(\mathbf{z})=\left[\begin{array}{c}\mathbf{-\mathbf{f}(\mathbf{x},\dot{\mathbf{x}})}\\\mathbf{0}\end{array}\right],\quad\mathbf{F}^{ext}(\mathbf{z},\mathbf{\phi})=\left[\begin{array}{c}\mathbf{f}^{ext}(\mathbf{\phi})\\\mathbf{0}\end{array}\right]$ .

DSorder = 2;
DS = DynamicalSystem(DSorder);
set(DS,'M',M,'C',C,'K',K);
set(DS.Options, 'Intrusion', 'none')

set(DS,'fnl_non',fint,'dfnl_non',dfint);
set(DS.Options,'Emax',5,'Nmax',10,'notation','multiindex')

We assume periodic forcing of the form

$\mathbf{f}^{ext}(\phi) = \mathbf{f}_0\cos(\phi)=\frac{\mathbf{f}_0}{2}e^{i\phi} + \frac{\mathbf{f}_0}{2}e^{-i\phi}$

Fourier coefficients of Forcing

f_0 = zeros(n,1);
f_0(outdof) = 1;
epsilon = 1.6;
masterModes = [1,2];
kappas = [-1; 1];
coeffs = [f_0 f_0]/2;

DS.add_forcing(coeffs, kappas,epsilon);

Linear Modal analysis and SSM setup

[V,D,W] = DS.linear_spectral_analysis();

Due to high-dimensionality, we compute only the first 5 eigenvalues with the smallest magnitude. These would also be used to compute the spectral quotients
Assuming a proportional damping hypthesis with symmetric matrices
modal damping ratio for 1 mode is 3.979552e-03
modal damping ratio for 2 mode is 6.463706e-03
modal damping ratio for 3 mode is 6.924137e-03
modal damping ratio for 4 mode is 1.008044e-02
modal damping ratio for 5 mode is 1.039666e-02
the left eigenvectors may be incorrect in case of asymmetry of matrices

 The first 10 nonzero eigenvalues are given as 
   1.0e+03 *

  -0.0051 + 1.2720i
  -0.0051 - 1.2720i
  -0.0134 + 2.0658i
  -0.0134 - 2.0658i
  -0.0153 + 2.2132i
  -0.0153 - 2.2132i
  -0.0325 + 3.2222i
  -0.0325 - 3.2222i
  -0.0345 + 3.3221i
  -0.0345 - 3.3221i

SSM setup

Choose Master subspace (perform resonance analysis)

S = SSM(DS);
S.choose_E(masterModes);
% %%
% setup options
set(S.Options, 'reltol', 1,'IRtol',0.1,'notation', 'multiindex','contribNonAuto',false)
set(S.FRCOptions, 'nt', 2^7, 'nRho', 300, 'nPar', 100, 'nPsi', 100, 'rhoScale', 2 )
set(S.FRCOptions, 'method', 'level set') %
set(S.FRCOptions, 'outdof',outdof)
set(DS.Options, 'outDOF',outdof)

% %%
% choose frequency range around the first natural frequency
omega0 = imag(S.E.spectrum(1));
omegaRange = omega0*[0.95 1.05];%
order = 7;

(near) outer resonance detected for the following combination of master eigenvalues
     2     0
     3     1
     4     2
     0     2
     1     3
     2     4
     2     0
     3     1
     4     2
     0     2
     1     3
     2     4
     3     0
     4     1
     0     3
     1     4
     3     0
     4     1
     0     3
     1     4

These are in resonance with the follwing eigenvalues of the slave subspace
   1.0e+03 *

  -0.0134 + 2.0658i
  -0.0134 + 2.0658i
  -0.0134 + 2.0658i
  -0.0134 - 2.0658i
  -0.0134 - 2.0658i
  -0.0134 - 2.0658i
  -0.0153 + 2.2132i
  -0.0153 + 2.2132i
  -0.0153 + 2.2132i
  -0.0153 - 2.2132i
  -0.0153 - 2.2132i
  -0.0153 - 2.2132i
  -0.0325 + 3.2222i
  -0.0325 + 3.2222i
  -0.0325 - 3.2222i
  -0.0325 - 3.2222i
  -0.0345 + 3.3221i
  -0.0345 + 3.3221i
  -0.0345 - 3.3221i
  -0.0345 - 3.3221i

sigma_out = 6
(near) inner resonance detected for the following combination of master eigenvalues
     2     1
     3     2
     1     2
     2     3

These are in resonance with the follwing eigenvalues of the master subspace
   1.0e+03 *

  -0.0051 + 1.2720i
  -0.0051 + 1.2720i
  -0.0051 - 1.2720i
  -0.0051 - 1.2720i

sigma_in = 6

compute SSM coefficients and reduced dynamics

[W_0, R_0] = S.compute_whisker(order);

save('SSMCoeff.mat','W_0','R_0')

Due to (near) outer resonance, the exisitence of the manifold is questionable and the underlying computation may suffer.
Attempting manifold computation
Warning: Integer operands are required for colon operator when used as index. 
Warning: Integer operands are required for colon operator when used as index. 
Manifold computation time at order 2 = 00:02:23
Estimated memory usage at order  2 = 1.84E+01 MB
Manifold computation time at order 3 = 00:05:05
Estimated memory usage at order  3 = 1.88E+01 MB
Manifold computation time at order 4 = 00:17:38
Estimated memory usage at order  4 = 2.06E+01 MB
Manifold computation time at order 5 = 00:31:43
Estimated memory usage at order  5 = 2.13E+01 MB
Manifold computation time at order 6 = 01:20:21
Estimated memory usage at order  6 = 2.29E+01 MB
Manifold computation time at order 7 = 02:02:37
Estimated memory usage at order  7 = 2.38E+01 MB

post-processing via SSM-based ROMs

cont_ep_run;

Due to high-dimensionality, we compute only the first 5 eigenvalues with the smallest magnitude. These would also be used to compute the spectral quotients
Assuming a proportional damping hypthesis with symmetric matrices
modal damping ratio for 1 mode is 3.979552e-03
modal damping ratio for 2 mode is 6.463706e-03
modal damping ratio for 3 mode is 6.924137e-03
modal damping ratio for 4 mode is 1.008044e-02
modal damping ratio for 5 mode is 1.039666e-02
the left eigenvectors may be incorrect in case of asymmetry of matrices

 The first 10 nonzero eigenvalues are given as 
   1.0e+03 *

  -0.0051 + 1.2720i
  -0.0051 - 1.2720i
  -0.0134 + 2.0658i
  -0.0134 - 2.0658i
  -0.0153 + 2.2132i
  -0.0153 - 2.2132i
  -0.0325 + 3.2222i
  -0.0325 - 3.2222i
  -0.0345 + 3.3221i
  -0.0345 - 3.3221i

No (near) outer resonances detected in the (truncated) spectrum
sigma_out = 6
(near) inner resonance detected for the following combination of master eigenvalues
     2     1
     3     2
     1     2
     2     3

These are in resonance with the follwing eigenvalues of the master subspace
   1.0e+03 *

  -0.0051 + 1.2720i
  -0.0051 + 1.2720i
  -0.0051 - 1.2720i
  -0.0051 - 1.2720i

sigma_in = 6

 Run='contep7.ep': Continue equilibria along primary branch.

    STEP   DAMPING               NORMS              COMPUTATION TIMES
  IT SIT     GAMMA     ||d||     ||f||     ||U||   F(x)  DF(x)  SOLVE
   0                          1.16e-02  1.80e+03    0.0    0.0    0.0
   1   1  1.00e+00  2.06e-02  1.33e-03  1.80e+03    0.0    0.0    0.0
   2   1  1.00e+00  4.16e-05  1.00e-08  1.80e+03    0.0    0.0    0.0
   3   1  1.00e+00  7.99e-11  5.33e-15  1.80e+03    0.0    0.0    0.0

 STEP      TIME        ||U||  LABEL  TYPE            om         rho1          th1          eps
    0  00:00:00   1.7989e+03      1  EP      1.2720e+03   6.7026e-02   3.5752e+00   1.6000e+00
   10  00:00:00   1.7863e+03      2          1.2631e+03   8.4096e-02   3.6970e+00   1.6000e+00
   20  00:00:00   1.7663e+03      3          1.2489e+03   1.0936e-01   3.8974e+00   1.6000e+00
   30  00:00:00   1.7463e+03      4          1.2348e+03   1.3221e-01   4.1199e+00   1.6000e+00
   40  00:00:01   1.7263e+03      5          1.2207e+03   1.5314e-01   4.4345e+00   1.6000e+00
   50  00:00:01   1.7190e+03      6          1.2155e+03   1.5901e-01   4.7596e+00   1.6000e+00
   57  00:00:01   1.7190e+03      7  SN      1.2155e+03   1.5895e-01   4.7668e+00   1.6000e+00
   57  00:00:01   1.7190e+03      8  FP      1.2155e+03   1.5895e-01   4.7671e+00   1.6000e+00
   60  00:00:01   1.7190e+03      9          1.2155e+03   1.5892e-01   4.7699e+00   1.6000e+00
   70  00:00:01   1.7190e+03     10          1.2155e+03   1.5879e-01   4.7832e+00   1.6000e+00
   80  00:00:01   1.7199e+03     11          1.2162e+03   1.5676e-01   4.8881e+00   1.6000e+00
   90  00:00:02   1.7383e+03     12          1.2291e+03   1.2975e-01   5.3321e+00   1.6000e+00
  100  00:00:02   1.7583e+03     13          1.2433e+03   9.5614e-02   5.6402e+00   1.6000e+00
  110  00:00:02   1.7690e+03     14          1.2508e+03   5.7677e-02   5.9133e+00   1.6000e+00
  119  00:00:02   1.7690e+03     15  FP      1.2508e+03   5.6254e-02   5.9228e+00   1.6000e+00
  119  00:00:02   1.7690e+03     16  SN      1.2508e+03   5.6254e-02   5.9228e+00   1.6000e+00
  120  00:00:02   1.7690e+03     17          1.2508e+03   5.6069e-02   5.9241e+00   1.6000e+00
  130  00:00:02   1.7690e+03     18          1.2508e+03   5.5055e-02   5.9308e+00   1.6000e+00
  140  00:00:03   1.7673e+03     19          1.2496e+03   4.4080e-02   6.0033e+00   1.6000e+00
  150  00:00:03   1.7555e+03     20          1.2413e+03   2.7627e-02   6.1092e+00   1.6000e+00
  160  00:00:03   1.7355e+03     21          1.2272e+03   1.8241e-02   6.1686e+00   1.6000e+00
  170  00:00:03   1.7155e+03     22          1.2131e+03   1.3760e-02   6.1969e+00   1.6000e+00
  174  00:00:03   1.7090e+03     23  EP      1.2084e+03   1.2741e-02   6.2033e+00   1.6000e+00

 STEP      TIME        ||U||  LABEL  TYPE            om         rho1          th1          eps
    0  00:00:03   1.7989e+03     24  EP      1.2720e+03   6.7026e-02   3.5752e+00   1.6000e+00
   10  00:00:03   1.8115e+03     25          1.2810e+03   5.0501e-02   3.4637e+00   1.6000e+00
   20  00:00:03   1.8315e+03     26          1.2951e+03   3.1126e-02   3.3379e+00   1.6000e+00
   30  00:00:04   1.8515e+03     27          1.3092e+03   2.0913e-02   3.2730e+00   1.6000e+00
   40  00:00:04   1.8715e+03     28          1.3234e+03   1.5477e-02   3.2387e+00   1.6000e+00
   49  00:00:04   1.8889e+03     29  EP      1.3356e+03   1.2585e-02   3.2205e+00   1.6000e+00