Skip to content

Commit

Permalink
SSMTool 2.4
Browse files Browse the repository at this point in the history
  • Loading branch information
jain-shobhit committed Jul 18, 2023
1 parent 7df7d0b commit e369573
Show file tree
Hide file tree
Showing 318 changed files with 20,530 additions and 2,563 deletions.
41 changes: 27 additions & 14 deletions README.md
Original file line number Diff line number Diff line change
@@ -1,9 +1,10 @@
# SSMTool 2.3: Computation of invariant manifolds in high-dimensional mechanics problems
# SSMTool 2.4: Computation of invariant manifolds in high-dimensional mechanics problems
[![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.4614201.svg)](https://doi.org/10.5281/zenodo.4614201)

What's new in SSMTool 2.3
- Extension to constrained mechanical systems
- Demonstration to a fluid-structure interaction system with flow-induced gyroscopic and follower forces
What's new in SSMTool 2.4
- Computation of higher-order approximations to non-autonomous invariant manifolds,
- Applications to computation of stability diagrams and forced response curves in nonlinear mechanical systems under parametric resonance,
- Improvements in computation using the multi-index notation instead of the tensor notation.

This package computes invariant manifolds in high-dimensional dynamical systems using the *Parametrization Method* with special attention to the computation of Spectral Submanifolds (SSM) and forced response curves in finite element models.

Expand All @@ -22,6 +23,10 @@ How SSMs are extended to constrained mechanical systems are discussed in the fol
[4] **Li, M., Jain, S. & Haller, G. Model reduction for constrained mechanical systems via spectral submanifolds. To appear on Nonlinear Dyn (2023).
https://doi.org/10.48550/arXiv.2208.03119**

The treatment of systems subject to parametric resonance via higher-order approximations of nonautonomous SSMs is described in the following article:

[5] **Thurnher, T., Haller, G. & Jain, S. Nonautonomous spectral submanifolds for model reduction of nonlinear mechanical systems under parametric resonances. Preprint (2023)**

In this version, we demonstrate the computational methodology over the following small academic examples as well high-dimensional finite element problems using the FE package *YetAnotherFECode*

First-order examples:
Expand All @@ -31,22 +36,30 @@ First-order examples:
- Complex Dyn: example of 4D first-order dynamical system with complex coefficients benchmarked against SSMTool 1.0.

Second-order examples:
- Oscillator chain: two coupled oscillators with 1:2 internal resonance, three coupled oscillators with 1:1:1 internal resonance and n coupled oscillators without any internal resonances.
- Bernoulli beam: modeled using linear finite elements with localized nonlinearity in the form of a cubic spring with and without 1:3 **internal resonances** (IR) and demonstration of bifurcation to quasiperiodic response on a 3D torus.
- von Karman straight beam in 2D: geometrically nonlinear finite element model with and without internal resonance (1:3) and demonstration of bifurcation to quasiperiodic response on a 2D torus.
- von Karman plate in 3D: geometrically nonlinear finite element model of a square flat plate with demonstration of **parallel computing**, 1:1 internal resonance and bifurcation to quasiperiodic response on a 2D torus.
- von Karman shell-based shallow curved panel in 3D: geometrically nonlinear finite element model with and without 1:2 internal resonance.
- Oscillator chain: two coupled oscillators with 1:2 internal resonance, three coupled oscillators with 1:1:1 internal resonance and n coupled oscillators without any internal resonances. [1]
- Bernoulli beam: modeled using linear finite elements with localized nonlinearity in the form of a cubic spring with and without 1:3 **internal resonances** (IR) and demonstration of bifurcation to quasiperiodic response on a 3D torus. [3]
- von Karman straight beam in 2D: geometrically nonlinear finite element model with and without internal resonance (1:3) and demonstration of bifurcation to quasiperiodic response on a 2D torus. [1,2,3]
- von Karman plate in 3D: geometrically nonlinear finite element model of a square flat plate with demonstration of **parallel computing**, 1:1 internal resonance and bifurcation to quasiperiodic response on a 2D torus. [2,3]
- von Karman shell-based shallow curved panel in 3D: geometrically nonlinear finite element model with and without 1:2 internal resonance. [2]
- Prismatic beam: nonlinear beam PDE discretized using Galerkin method onto a given number of modes, comparison with the method of multiple scales, demonstration of 1:3 internal internal resonance
- AxialMovingBeam: an axially moving beam with gyroscopic and nonlinear damping forces
- PipeConveyingFluid: an **fluid-structure interaction** system with flow-induced gyroscopic and follower forces, demonstration of **asymmetric** damping and stiffness matrices and **global bifurcation**. More details can be found in https://doi.org/10.1016/j.ymssp.2022.109993
- TimoshenkoBeam: a cantilever Timoshenko beam carrying a lumped mass. This example demonstrates the effectiveness of SSM reduction for systems undergoing large deformations.
- NACA airfoil based aircraft wing model: shell-based nonlinear finite element model containing more than 100,000 degrees of freedom.
- AxialMovingBeam: an axially moving beam with gyroscopic and nonlinear damping forces [2]
- TimoshenkoBeam: a cantilever Timoshenko beam carrying a lumped mass. This example demonstrates the effectiveness of SSM reduction for systems undergoing large deformations. [2]
- NACA airfoil based aircraft wing model: shell-based nonlinear finite element model containing more than 100,000 degrees of freedom. [1]

Constrained mechanical systems
Constrained mechanical systems [4]
- 3D oscillator constrained in a surface
- Pendulums in Cartesian coordinates: single pendulum, pendulum-slider with 1:3 internal resonance, and a chain of pendulums
- Frequency divider: two geometrically nonlinear beams connected via a revolute joint

Computation of stability diagrams and forced response curves in mechanical systems under parametric resonance:
- 1-DOF Mathieu oscillator with periodically-varying linear and cubic spring
- Coupled system of two Mathieu oscillators with cubic nonlinearity [5]
- **Self-excited** 2-DOF-oscillator system with cubic nonlinearity, time-varying stiffness as well as external periodic forcing, where the forced response exhibits **multiple isolas**
- Self-excited 2-DOF-oscillator system with external and parametric excitation [5]
- **Parametric amplification** in coupled-oscillators due to nonlinear damping and time-varying stiffness and external periodic forcing with a phase lag relative to stiffness.
- Prismatic beam: nonlinear beam PDE discretized using Galerkin method onto a given number of modes, demonstration of axial stretching leading to parametric excitation of transverse degrees of freedom. [5]
- Bernoulli beam attached to a spring time-varying linear stiffness along with a nonlinear damper and spring: Stability diagrams as well as forced response curves (exhibiting isolas) are constructed over two-dimensional SSMs. [5]

This package uses the following external open-source packages:

1. Continuation core (coco) https://sourceforge.net/projects/cocotools/
Expand Down
Binary file modified examples/BenchamrkSSM1stOrder/demo.mlx
Binary file not shown.
Binary file modified examples/BernoulliBeam/BernoulliBeamBaseExcitation.mlx
Binary file not shown.
2 changes: 1 addition & 1 deletion examples/BernoulliBeam/build_model.m
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@
subs3 = [[nDOF-1,2*nDOF- 1,2*nDOF- 1,2*nDOF- 1];
[nDOF-1,nDOF-1,nDOF-1,nDOF-1]];
vals3 = [gamma;kappa];
f3 = sptensor(subs3, vals3, [2*nDOF,2*nDOF,2*nDOF,2*nDOF]);
f3 = sptensor(subs3, vals3, [nDOF,2*nDOF,2*nDOF,2*nDOF]);


%%
Expand Down
Binary file modified examples/BernoulliBeamIRs/demo.mlx
Binary file not shown.
Binary file modified examples/CharneyDeVore1stOrder/demo.mlx
Binary file not shown.
Binary file modified examples/ComplexDyn/demo.mlx
Binary file not shown.
Binary file modified examples/NACAWing/NACAWingWorkbook.mlx
Binary file not shown.
2 changes: 1 addition & 1 deletion examples/NACAWing/build_model.m
Original file line number Diff line number Diff line change
Expand Up @@ -57,7 +57,7 @@
V = MyAssembly.unconstrain_vector(V0);
mod = 1;
v1 = reshape(V(:,mod),6,[]);
PlotFieldonDeformedMesh(nodes,elements,v1(1:3,:).','factor',50)
PlotFieldonDeformedMesh(nodes,elements,v1(1:3,:).','factor',50);
title(['Mode ' num2str(mod) ', Frequency = ' num2str(omega(mod)/(2*pi)) ' Hz'] )

%% Damping matrix
Expand Down
7 changes: 5 additions & 2 deletions examples/OscillatorChain/build_model.m
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
function [M,C,K,fnl,f_ext] = build_model(n,m,c,k,kappa2,kappa3)
function [M,C,K,fnl,fext] = build_model(n,m,c,k,kappa2,kappa3)

[K,C,f2,f3] = assemble_global_coefficients(k,kappa2,kappa3,c,n);
%%
Expand All @@ -12,5 +12,8 @@

fnl = {f2, f3};

f_ext = ones(n,1);
f_0 = ones(n,1);
fext.data = forcing(n,f_0);
% f_ext = 0.1*ones(n,1)+8*(1:n)'/n; % ones(n,1) has no contribution to the second mode
end

37 changes: 37 additions & 0 deletions examples/OscillatorChain/demo_single_master_mode.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,37 @@
% %

n = 4;
m = 1;
k = 1;
c = 2;
kappa2 = 1;
kappa3 = 0.5;

[M,C,K,fnl,~] = build_model(n,m,c,k,kappa2,kappa3);
%% Dynamical system setup
DS = DynamicalSystem();
% set(DS,'M',M,'C',C,'K',K,'fnl',fnl); % second order works fine
A = [-K zeros(n);zeros(n) M];
B = [C M;M zeros(n)];
set(DS,'A',A,'B',B,'fnl',{-fnl{1},-fnl{2}});
% set(DS.Options,'Emax',5,'Nmax',10,'notation','multiindex')
set(DS.Options,'Emax',5,'Nmax',100,'notation','tensor')
%%
epsilon = 5e-3;
f_0 = ones(n,1);
%%
% Fourier coefficients of Forcing
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();
%%
S = SSM(DS);
% set(S.Options, 'reltol', 0.1,'notation','multiindex')
set(S.Options, 'reltol', 10,'notation','tensor')
masterModes = [7 8]; % single mode will fail
order = 5;
S.choose_E(masterModes);
[W_0,R_0] = S.compute_whisker(order);
12 changes: 12 additions & 0 deletions examples/OscillatorChain/forcing.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,12 @@
function [data] = forcing(n,f_0)
% Cosine type forcing of the oscillatory chain with forcing amplitude
% vector f_0

data(1).kappa = 1;
data(2).kappa = -1;
data(1).f_n_k(1).coeffs = f_0/2;
data(1).f_n_k(1).ind = zeros(1,n);
data(2).f_n_k(1).coeffs = f_0/2;
data(2).f_n_k(1).ind = zeros(1,n);

end
205 changes: 205 additions & 0 deletions examples/ParametricResonance/BernoulliBeamParametric/BBplotFRC.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,205 @@
function BBplotFRC(FRC,bd,order,outdof,varargin)
% FRC data from SSM
% bd data from coco
numOutdof = numel(outdof);


%% Prepare SSM data
Par = [FRC.Omega];
Aout = [FRC.Aout];
stab_ssm = [FRC.stability];

numPts = numel(stab_ssm);
Aout = reshape(Aout, [numOutdof, numPts]);
Aout = Aout';

%% Prepare COCO data
nSample = 1;

% extract results
omega = [];
stab = logical([]);
amp = cell(numOutdof,1);
omegai = coco_bd_col(bd, 'omega');
stabi = coco_bd_col(bd, 'eigs')';
stabi = abs(stabi);
stabi = all(stabi<1, 2);
omega = [omega, omegai];
stab = [stab; stabi];

ampNames = cell(1, numel(numOutdof));
for k = 1:numel(outdof)
ampNames{k} = strcat('amp',num2str(outdof(k)));
end

for j=1:numOutdof
ampj = coco_bd_col(bd, ampNames{j});
amp{j} = [amp{j}, ampj];
end

if ~isempty(varargin)
% Second Coco bd input
bd2 = varargin{1};

omega2 = [];
stab2 = logical([]);
amp2 = cell(numOutdof,1);
omegai2 = coco_bd_col(bd2, 'omega');
stabi2 = coco_bd_col(bd2, 'eigs')';
stabi2 = abs(stabi2);
stabi2 = all(stabi2<1, 2);
omega2 = [omega2, omegai2];
stab2 = [stab2; stabi2];

ampNames2 = cell(1, numel(numOutdof));
for k = 1:numel(outdof)
ampNames2{k} = strcat('amp',num2str(outdof(k)));
end

for j=1:numOutdof
ampj2 = coco_bd_col(bd2, ampNames2{j});
amp2{j} = [amp2{j}, ampj2];
end

end

%% Plotting parameters

%% Single Plot Setting

%axis_size = 20; % Set legend and Axis size
%xwidth = 500; % window size
%ywidth = 400;

%ssm_marker = 10*1.4;
%coco_marker = 6*1.4;


%% Double Plot Setting
% {
axis_size = 20; % Set legend and Axis size
xwidth = 600; % window size
ywidth = 500;

ssm_marker = 7*1.4;
coco_marker = 6*1.4;
%}

ustabc1 = [0.9290 0.6940 0.1250];
stabc1 = [0.4940 0.1840 0.5560];
ustabc2 = [0.4660 0.6740 0.1880];
stabc2 = [0.3010 0.7450 0.9330];
%% Start Plotting
for k=1:numOutdof
set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex');
figure(); hold on



% Plot Coco Data
om = omega(~stab); am = amp{k}(~stab);
plot(om(1:nSample:end), am(1:nSample:end), '*', 'MarkerSize',coco_marker,'color',ustabc1);

om = omega(stab); am = amp{k}(stab);
plot(om(1:nSample:end), am(1:nSample:end), '*', 'MarkerSize',coco_marker,'color',stabc1);


% Plot second set of Coco data
if ~isempty(varargin)
om2 = omega2(~stab2); am2 = amp2{k}(~stab2);
plot(om2(1:nSample:end), am2(1:nSample:end), '*', 'MarkerSize',coco_marker,'color',ustabc2);
% Plot Coco Data
om2 = omega2(stab2); am2 = amp2{k}(stab2);
plot(om2(1:nSample:end), am2(1:nSample:end), '*', 'MarkerSize',coco_marker,'color',stabc2);

%'color', [0.9294 0.6941 0.1255],
end


% Plot SSM Data
plot(Par(stab_ssm),Aout(stab_ssm,k),'bo','MarkerSize',ssm_marker);
plot(Par(~stab_ssm),Aout(~stab_ssm,k),'or','MarkerSize',ssm_marker);
add_labels('$\Omega$',strcat('$||z_{',num2str(outdof(k)),'}||_{\infty}$'))

% Add legend
if isempty(varargin)
add_legends(stab_ssm,order, get_cocolabel(stab))
else
add_legends(stab_ssm,order, get_cocolabel(stab,stab2))

end


%% Set axis limits by SSM axis
xmin = min(Par);
xmax = max(Par);
ymin = min(Aout(:,k));
ymax = max(Aout(:,k))*1.1;

xlim([xmin,xmax]);
ylim([ymin,ymax]);

%% Parameters for figure size, label size, legend size
set(gca,'fontsize',axis_size)

%% Size of window
set(gcf, 'Position', [100 100 xwidth ywidth])
end

end

function cclabel = get_cocolabel(stab,varargin)
cclabel = {'a'};


if all(stab)
cclabel{end+1} = 'collocation-stable' ;
elseif ~any(stab)
cclabel{end+1} = 'collocation-unstable';
else
cclabel{end+1} = 'collocation-unstable';
cclabel{end+1} = 'collocation-stable' ;

end

if ~isempty(varargin)
stab2 = varargin{1};
if all(stab2)
cclabel{end+1} = 'collocation-stable' ;
elseif ~any(stab2)
cclabel{end+1} = 'collocation-unstable';
else
cclabel{end+1} = 'collocation-unstable';
cclabel{end+1} = 'collocation-stable' ;

end


end

cclabel = cclabel(2:end);

end

function add_legends(stab,order,cclabel)

if all(stab)
legend(cclabel{:},strcat('SSM-$$\mathcal{O}(',num2str(order),')$$ - stable'),...
'Interpreter','latex','Location','east');
elseif ~any(stab)
legend(cclabel{:},strcat('SSM-$$\mathcal{O}(',num2str(order),')$$ - unstable'),...
'Interpreter','latex','Location','east');
else
legend(cclabel{:},strcat('SSM-$$\mathcal{O}(',num2str(order),')$$ - stable'),...
strcat('SSM-$$\mathcal{O}(',num2str(order),')$$ - unstable'),...
'Interpreter','latex','Location','east');
end

end

function add_labels(xlab,ylab)
xlabel(xlab,'Interpreter','latex');
ylabel(ylab,'Interpreter','latex');
set(gca,'FontSize',14);
grid on, axis tight; legend boxoff;
end
Loading

0 comments on commit e369573

Please sign in to comment.