Euler Bernoulli beam with cubic spring and damper

Generate model
Dynamical system setup
Linear Modal Analysis
Forced response curves using SSMs
Verification: Collocation using coco
Verification of Isola

Additionally to external excitation on the last node parametric excitation is applied as well. The tip of the beam is subject to linear parametric excitation. The experiment this example is based on can be found in

Chen, C. C. & Yeh, M. K.: Parametric instability of a beam under electromagnetic excitation. Journal of Sound and Vibration 240,747–764, https://doi.org/10.1006/jsvi.2000.3255; A schematic depiction of the model is given by

Generate model

clear all
nElements = 5;
kappa = 50; % cubic spring
gamma = 0.01; % cubic damping
mu = 45;
[M,C,K,fnl,fext] = build_model_external(kappa, gamma, nElements,mu);

n = length(M);

Dynamical system setup

We consider the parametrically excited system

$\mathbf{M}\ddot{\mathbf{x}}+\mathbf{C}\dot{\mathbf{x}}+\mathbf{K}\mathbf{x}+\mathbf{f}(\mathbf{x},\dot{\mathbf{x}})=\epsilon\mathbf{g}(\mathbf{x},\dot{\mathbf{x}},\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{G}(\mathbf{z},\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{G}(\mathbf{z},\mathbf{\phi})=\left[\begin{array}{c}\mathbf{g}(\mathbf{x, \dot{x}},\phi)\\\mathbf{0}\end{array}\right]$

order = 2;
DS = DynamicalSystem(order);
set(DS,'M',M,'C',C,'K',K,'fnl',fnl);
set(DS.Options,'Emax',5,'Nmax',10,'notation','multiindex')

epsilon = 0.002;;
DS.add_forcing(fext,epsilon);

Linear Modal Analysis

% Analyse spectrum
[V,D,W_evec] = DS.linear_spectral_analysis();

% Choose Master subspace (perform resonance analysis)

% Set up SSM object
S = SSM(DS);
set(S.Options, 'reltol', 0.3,'notation','multiindex')

%Choose Master subspace
resModes = [1,2];
S.choose_E(resModes);

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 8.840017e-04
modal damping ratio for 2 mode is 5.488086e-03
modal damping ratio for 3 mode is 1.541080e-02
modal damping ratio for 4 mode is 3.044305e-02
modal damping ratio for 5 mode is 5.052761e-02

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

  -0.0001 + 0.0700i
  -0.0001 - 0.0700i
  -0.0024 + 0.4389i
  -0.0024 - 0.4389i
  -0.0190 + 1.2327i
  -0.0190 - 1.2327i
  -0.0741 + 2.4343i
  -0.0741 - 2.4343i
  -0.2042 + 4.0370i
  -0.2042 - 4.0370i

(near) outer resonance detected for the following combination of master eigenvalues
     6     0
     7     1
     8     2
     0     6
     1     7
     2     8

These are in resonance with the follwing eigenvalues of the slave subspace
  -0.2409 +43.8926i
  -0.2409 +43.8926i
  -0.2409 +43.8926i
  -0.2409 -43.8926i
  -0.2409 -43.8926i
  -0.2409 -43.8926i

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

These are in resonance with the follwing eigenvalues of the master subspace
  -0.0062 + 7.0006i
  -0.0062 + 7.0006i
  -0.0062 + 7.0006i
  -0.0062 + 7.0006i
  -0.0062 - 7.0006i
  -0.0062 - 7.0006i
  -0.0062 - 7.0006i
  -0.0062 - 7.0006i

sigma_in = 3300

Forced response curves using SSMs

Obtaining forced response curve in reduced-polar coordinate

order = 5; % SSM approximation order
outdof = [n-2 n-1]; % degree of freedom at which output is displayed

setup options

set(S.Options, 'reltol', 1,'IRtol',0.02,'notation', 'multiindex','contribNonAuto',true)
set(S.FRCOptions, 'nt', 2^7, 'nRho', 200, 'nPar', 200, 'nPsi', 100, 'rhoScale', 2 )
set(S.FRCOptions, 'method','level set')
set(S.FRCOptions, 'outdof',outdof)

choose frequency range

omega0 = imag(S.E.spectrum(1));
OmegaRange = omega0*[0.95 1.05];

Extract forced response curve

startFRCSSM = tic;
FRC_SSM = S.extract_FRC('freq',OmegaRange,order);
timings.FRCSSM = toc(startFRCSSM)

*****************************************
Calculating FRC using SSM with master subspace: [1  2]
(near) outer resonance detected for the following combination of master eigenvalues
     6     0
     7     0
     7     1
     8     1
     8     2
     0     6
     0     7
     1     7
     1     8
     2     8

These are in resonance with the follwing eigenvalues of the slave subspace
  -0.2409 +43.8926i
  -0.2409 +43.8926i
  -0.2409 +43.8926i
  -0.2409 +43.8926i
  -0.2409 +43.8926i
  -0.2409 -43.8926i
  -0.2409 -43.8926i
  -0.2409 -43.8926i
  -0.2409 -43.8926i
  -0.2409 -43.8926i

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

These are in resonance with the follwing eigenvalues of the master subspace
  -0.0062 + 7.0006i
  -0.0062 + 7.0006i
  -0.0062 + 7.0006i
  -0.0062 + 7.0006i
  -0.0062 - 7.0006i
  -0.0062 - 7.0006i
  -0.0062 - 7.0006i
  -0.0062 - 7.0006i

sigma_in = 3300
Due to (near) outer resonance, the exisitence of the manifold is questionable and the underlying computation may suffer.
Attempting manifold computation
Manifold computation time at order 2 = 00:00:00
Estimated memory usage at order  2 = 1.62E-02 MB
Manifold computation time at order 3 = 00:00:00
Estimated memory usage at order  3 = 2.19E-02 MB
Manifold computation time at order 4 = 00:00:00
Estimated memory usage at order  4 = 3.17E-02 MB
Manifold computation time at order 5 = 00:00:00
Estimated memory usage at order  5 = 4.43E-02 MB
gamma = 
   1.0e+05 *

  -0.0005 + 0.0074i
  -0.0286 - 1.2825i

Total time spent on FRC computation upto O(5) = 00:00:12

timings = 

  struct with fields:

    FRCSSM: 13.9056

Verification: Collocation using coco

Dankowicz, H., & Schilder, F. (2013). Recipes for Continuation, SIAM Philadelphia. https://doi.org/10.1137/1.9781611972573

nCycles = 10;

coco = cocoWrapper(DS, nCycles, outdof);
set(coco,'initialGuess','linear')
set(coco.Options, 'NAdapt', 1);
set(coco.Options,'NTST', 20,'PtMX',200);

startcoco = tic;
bd = coco.extract_FRC(OmegaRange);
timings.cocoFRC = toc(startcoco)

 Run='FRC0.002': Continue primary family of periodic orbits.

    STEP   DAMPING               NORMS              COMPUTATION TIMES
  IT SIT     GAMMA     ||d||     ||f||     ||U||   F(x)  DF(x)  SOLVE
   0                          4.78e-03  9.50e+00    0.0    0.0    0.0
   1   1  1.00e+00  1.88e-02  7.95e-07  9.50e+00    0.0    0.1    0.0
   2   1  1.00e+00  4.21e-06  3.36e-13  9.50e+00    0.0    0.1    0.0
   3   1  1.00e+00  1.68e-13  1.58e-13  9.50e+00    0.0    0.2    0.0

 STEP      TIME        ||U||  LABEL  TYPE         omega    po.period          eps         amp8         amp9
    0  00:00:00   9.5028e+00      1  EP      6.6506e+00   9.4476e-01   2.0000e-03   1.9201e-06   3.7842e-03
   10  00:00:08   9.8063e+00      2          6.8627e+00   9.1555e-01   2.0000e-03   5.3478e-06   1.0577e-02
   20  00:00:14   1.0851e+01      3          7.0895e+00   8.8626e-01   2.0000e-03   3.8253e-05   7.6031e-02
   30  00:00:20   1.2245e+01      4          7.2575e+00   8.6575e-01   2.0000e-03   5.6644e-05   1.1310e-01
   40  00:00:27   1.2671e+01      5          7.3329e+00   8.5685e-01   2.0000e-03   6.2462e-05   1.2496e-01
   44  00:00:32   1.2692e+01      6  FP      7.3389e+00   8.5614e-01   2.0000e-03   6.2451e-05   1.2500e-01
   44  00:00:32   1.2692e+01      7  SN      7.3389e+00   8.5614e-01   2.0000e-03   6.2451e-05   1.2500e-01
   50  00:00:36   1.2509e+01      8          7.3180e+00   8.5859e-01   2.0000e-03   5.9473e-05   1.1897e-01
   60  00:00:43   1.1514e+01      9          7.2060e+00   8.7194e-01   2.0000e-03   4.5680e-05   9.1122e-02
   70  00:00:50   1.0162e+01     10          7.0506e+00   8.9115e-01   2.0000e-03   1.4169e-05   2.8139e-02
   71  00:00:52   1.0114e+01     11  FP      7.0471e+00   8.9160e-01   2.0000e-03   1.1001e-05   2.1838e-02
   71  00:00:52   1.0113e+01     12  SN      7.0471e+00   8.9160e-01   2.0000e-03   1.0986e-05   2.1809e-02
   80  00:00:57   1.0092e+01     13          7.0660e+00   8.8922e-01   2.0000e-03   6.7889e-06   1.3475e-02
   90  00:01:02   1.0217e+01     14          7.1689e+00   8.7645e-01   2.0000e-03   3.0391e-06   6.0430e-03
   97  00:01:04   1.0466e+01     15  EP      7.3506e+00   8.5478e-01   2.0000e-03   1.5409e-06   3.0741e-03

timings = 

  struct with fields:

     FRCSSM: 13.9056
    cocoFRC: 66.9466

Verification of Isola

Using contiunation of an initial periodic orbit, the isola detected by SSMtheory can be verified with coco.

% Initial Condition on Isola obtained from SSM
load('FRC_IC.mat')

IC = FRC.Zic.';
OmegaIC = FRC.Omega;

nCycles = 10;

coco = cocoWrapper(DS, nCycles, outdof);
set(coco,'initialGuess','linear')
set(coco.Options, 'NAdapt', 1);
set(coco.Options,'NTST', 80,'PtMX',90,'bi_direct',false);

% Extract Isola
startcoco = tic;
bd_isola = coco.extract_FRC_fromIC(OmegaRange,IC,OmegaIC);
timings.cocoFRCisola = toc(startcoco)

 Run='FRC0.002': Continue primary family of periodic orbits.

    STEP   DAMPING               NORMS              COMPUTATION TIMES
  IT SIT     GAMMA     ||d||     ||f||     ||U||   F(x)  DF(x)  SOLVE
   0                          4.45e-06  1.19e+01    0.0    0.0    0.0
   1   1  1.00e+00  3.38e-03  2.12e-08  1.19e+01    0.0    0.2    0.1
   2   1  1.00e+00  8.64e-06  1.54e-12  1.19e+01    0.0    0.5    0.1
   3   1  1.00e+00  6.84e-11  8.65e-13  1.19e+01    0.0    0.7    0.2

 STEP      TIME        ||U||  LABEL  TYPE         omega    po.period          eps         amp8         amp9
    0  00:00:01   1.1871e+01      1  EP      7.0556e+00   8.9053e-01   2.0000e-03   2.3854e-05   4.7361e-02
    1  00:00:03   1.1939e+01      2  EP      7.0556e+00   8.9053e-01   2.0000e-03   2.3854e-05   4.7360e-02

    STEP   DAMPING               NORMS              COMPUTATION TIMES
  IT SIT     GAMMA     ||d||     ||f||     ||U||   F(x)  DF(x)  SOLVE
   0                          8.74e-13  1.19e+01    0.0    0.0    0.0

 STEP      TIME        ||U||  LABEL  TYPE         omega    po.period          eps         amp8         amp9
    0  00:00:02   1.1904e+01      1  EP      7.0556e+00   8.9053e-01   2.0000e-03   2.3854e-05   4.7361e-02
   10  00:00:47   1.3250e+01      2          7.1078e+00   8.8398e-01   2.0000e-03   3.5079e-05   6.9747e-02
   20  00:01:20   1.3004e+01      3          7.1467e+00   8.7917e-01   2.0000e-03   4.0403e-05   8.0420e-02
   27  00:01:33   1.2544e+01      4  SN      7.1512e+00   8.7862e-01   2.0000e-03   4.0652e-05   8.0933e-02
   27  00:01:34   1.2544e+01      5  FP      7.1512e+00   8.7862e-01   2.0000e-03   4.0652e-05   8.0933e-02
   30  00:01:37   1.2315e+01      6          7.1502e+00   8.7875e-01   2.0000e-03   4.0302e-05   8.0233e-02
   40  00:01:48   1.1434e+01      7          7.1333e+00   8.8082e-01   2.0000e-03   3.7139e-05   7.3915e-02
   50  00:01:56   1.0532e+01      8          7.0811e+00   8.8732e-01   2.0000e-03   2.7080e-05   5.3824e-02
   60  00:02:02   1.0301e+01      9          7.0512e+00   8.9108e-01   2.0000e-03   2.1007e-05   4.1708e-02
   66  00:02:07   1.0306e+01     10  FP      7.0504e+00   8.9118e-01   2.0000e-03   2.1268e-05   4.2218e-02
   66  00:02:07   1.0306e+01     11  SN      7.0504e+00   8.9118e-01   2.0000e-03   2.1268e-05   4.2218e-02
   70  00:02:10   1.0322e+01     12          7.0510e+00   8.9110e-01   2.0000e-03   2.1925e-05   4.3527e-02
   80  00:02:15   1.0530e+01     13          7.0721e+00   8.8845e-01   2.0000e-03   2.8290e-05   5.6198e-02
   90  00:02:21   1.0985e+01     14  EP      7.1340e+00   8.8074e-01   2.0000e-03   3.8887e-05   7.7372e-02

timings = 

  struct with fields:

          FRCSSM: 13.9056
         cocoFRC: 66.9466
    cocoFRCisola: 158.8883

Plot for paper

BBplotFRC(FRC_SSM,bd{1},order,outdof,bd_isola{1})