SSMTool 2.2

examples/AxialMovingBeam/AxialMovingBeam.m

@@ -0,0 +1,271 @@
+%% 
+% We consider an axially moving beam under 1:3 internal resonance between the 
+% first two bending modes. The equation of motion is given by
+% 
+% $$\ddot{\mathbf{u}}+(\mathbf{C}+\mathbf{G})\dot{\mathbf{u}}+\mathbf{f}(\mathbf{u},\dot{\mathbf{u}})=\epsilon\Omega^2\mathbf{g}\cos\Omega 
+% t$$
+% 
+% where $\mathbf{G}^\top=-\mathbf{G}$ is a gyroscopic matrix. With viscoelastic 
+% material model, the system has nonlinear damping. In addition, the beam is subject 
+% to base excitation such that the forcing amplitude is a function of forcing 
+% frequency.
+% 
+% In this notebook, we will calculate the forced response curve of the system 
+% with/without the nonlinear damping taken into consideration to explore the effects 
+% of nonlinear damping
+%% Nonlinear Damping
+% Setup Dynamical System
+
+clear all;
+n = 10;
+[mass,damp,gyro,stiff,fnl,fext] = build_model(n,'nonlinear_damp');
+
+%% 
+% *Create model*
+%% 
+% * Given the assembled damping matrix $\mathbf{C}+\mathbf{G}$ is not Rayleigh 
+% damping anymore. One should set the corresponding Options in DS to be false.
+% * The BaseExcitation option in DS should be set true to account for the $\Omega$-dependent 
+% forcing amplitude.
+
+DS = DynamicalSystem();
+set(DS,'M',mass,'C',damp+gyro,'K',stiff,'fnl',fnl);
+set(DS.Options,'Emax',6,'Nmax',10,'notation','multiindex');
+set(DS.Options,'RayleighDamping',false,'BaseExcitation',true);
+
+% Forcing
+h = 1.5e-4; % h characterizes the vibration amplitude of base excitation. It plays the role of epsilon here
+kappas = [-1; 1];
+coeffs = [fext fext]/2;
+DS.add_forcing(coeffs, kappas, h);
+% Linear Modal analysis
+
+[V,D,W] = DS.linear_spectral_analysis();
+%% 
+% *Choose Master subspace*
+% 
+% Due to the 1:3 internal resonance, we take the first two complex conjugate 
+% pairs of modes as the spectral subspace to SSM. So we hvae resonant_modes = 
+% [1 2 3 4].
+
+S = SSM(DS);
+set(S.Options, 'reltol', 1,'notation','multiindex');
+resonant_modes = [1 2 3 4];
+order = 3;
+outdof = [1 2];
+% Primary resonance of the first mode
+% We consider the case that $\Omega\approx\omega_1$. Although the second mode 
+% is not excited externally ($f_2=\langle \phi_2,f^\mathrm{ext}\rangle=0$), the 
+% response of the second mode is nontrivial due to the modal interactions.
+
+freqrange = [0.98 1.04]*imag(D(1));
+set(S.FRCOptions, 'nCycle',500, 'initialSolver', 'fsolve');
+set(S.contOptions, 'PtMX', 300, 'h0', 0.1, 'h_max', 0.2, 'h_min', 1e-3);
+set(S.FRCOptions, 'coordinates', 'polar');
+%% 
+% We first compute the FRC with O(3,5,7) expansion of SSM and then check the 
+% convergence of FRC with increasing orders. 
+
+% O(3)
+start = tic;
+FRC_ND_O3 = S.SSM_isol2ep('isol-nd-3',resonant_modes, order, [1 3], 'freq', freqrange,outdof);
+timings.FRC_ND_O3 = toc(start);
+% O(5)
+sol = ep_read_solution('isol-nd-3.ep',1);
+start = tic;
+FRC_ND_O5 = S.SSM_isol2ep('isol-nd-5',resonant_modes, order+2, [1 3],...
+    'freq', freqrange,outdof,{sol.p,sol.x});
+timings.FRC_ND_O5 = toc(start);
+fig_ssm = gcf;
+% O(7)
+start = tic;
+FRC_ND_O7 = S.SSM_isol2ep('isol-nd-7',resonant_modes, order+4, [1 3],...
+    'freq', freqrange,outdof,{sol.p,sol.x});
+timings.FRC_ND_O7 = toc(start);
+%% 
+% Plot FRC at different orders in the same figure to observe the convergence 
+
+FRCs = {FRC_ND_O3,FRC_ND_O5,FRC_ND_O7};
+thm = struct();
+thm.SN = {'LineStyle', 'none', 'LineWidth', 2, ...
+  'Color', 'cyan', 'Marker', 'o', 'MarkerSize', 8, 'MarkerEdgeColor', ...
+  'cyan', 'MarkerFaceColor', 'white'};
+thm.HB = {'LineStyle', 'none', 'LineWidth', 2, ...
+  'Color', 'black', 'Marker', 's', 'MarkerSize', 8, 'MarkerEdgeColor', ...
+  'black', 'MarkerFaceColor', 'white'};
+color = {'r','k','b','m'};
+figure(20);
+ax1 = gca;
+for k=1:3
+    FRC = FRCs{k};
+    SNidx = FRC.SNidx;
+    HBidx = FRC.HBidx;
+    FRC.st = double(FRC.st);
+    FRC.st(HBidx) = nan;
+    FRC.st(SNidx) = nan;
+    % color
+    ST = cell(2,1);
+    ST{1} = {[color{k},'--'],'LineWidth',1.5}; % unstable
+    ST{2} = {[color{k},'-'],'LineWidth',1.5};  % stable
+    legs = ['SSM-$\mathcal{O}(',num2str(2*k+1),')$-unstable'];
+    legu = ['SSM-$\mathcal{O}(',num2str(2*k+1),')$-stable'];
+    hold(ax1,'on');
+    plot_stab_lines(FRC.om,FRC.Aout_frc(:,1),FRC.st,ST,legs,legu);
+    SNfig = plot(FRC.om(SNidx),FRC.Aout_frc(SNidx,1),thm.SN{:});
+    set(get(get(SNfig,'Annotation'),'LegendInformation'),...
+    'IconDisplayStyle','off');
+    HBfig = plot(FRC.om(HBidx),FRC.Aout_frc(HBidx,1),thm.HB{:});
+    set(get(get(HBfig,'Annotation'),'LegendInformation'),...
+    'IconDisplayStyle','off');   
+    xlabel('$\Omega$','Interpreter','latex'); 
+    ylabel('$||u_1||_{\infty}$','Interpreter','latex'); 
+    set(gca,'FontSize',14);
+    grid on; axis tight; 
+end
+% Validation using collocation method from COCO
+
+figure(fig_ssm); hold on
+nCycles = 500;
+coco = cocoWrapper(DS, nCycles, outdof);
+set(coco.Options, 'PtMX', 1000, 'NTST',20, 'dir_name', 'bd_nd');
+set(coco.Options, 'NAdapt', 0, 'h_max', 200, 'MaxRes', 1);
+coco.initialGuess = 'linear';
+start = tic;
+bd_nd = coco.extract_FRC(freqrange);
+timings.cocoFRCbd_nd = toc(start)
+%% Linear Damping
+% Setup Dynamical System
+
+n = 10;
+[mass,damp,gyro,stiff,fnl,fext] = build_model(n,'linear_damp');
+
+% Create model
+DS = DynamicalSystem();
+set(DS,'M',mass,'C',damp+gyro,'K',stiff,'fnl',fnl);
+set(DS.Options,'Emax',6,'Nmax',10,'notation','multiindex');
+set(DS.Options,'RayleighDamping',false,'BaseExcitation',true);
+DS.add_forcing(coeffs, kappas, h);
+% Linear Modal analysis
+
+[V,D,W] = DS.linear_spectral_analysis();
+%% 
+% *Choose Master subspace* 
+
+S = SSM(DS);
+set(S.Options, 'reltol', 1,'notation','multiindex');
+resonant_modes = [1 2 3 4];
+order = 5;
+outdof = [1 2];
+% Primary resonance of the first mode
+
+freqrange = [0.98 1.04]*imag(D(1));
+set(S.FRCOptions, 'nCycle',500, 'initialSolver', 'fsolve');
+set(S.contOptions, 'PtMX', 300, 'h0', 0.1, 'h_max', 0.2, 'h_min', 1e-3);
+set(S.FRCOptions, 'coordinates', 'polar');
+
+sol = ep_read_solution('isol-nd-3.ep',1);
+start = tic;
+FRC_LD_O5 = S.SSM_isol2ep('isol-ld-5',resonant_modes, order, [1 3],...
+    'freq', freqrange,outdof,{sol.p,sol.x});
+timings.FRC_LD_O5 = toc(start);
+% validation using collocation method from COCO
+
+nCycles = 500;
+coco = cocoWrapper(DS, nCycles, outdof);
+set(coco.Options, 'PtMX', 1000, 'NTST',20, 'dir_name', 'bd_ld');
+set(coco.Options, 'NAdapt', 0, 'h_max', 200, 'MaxRes', 1);
+coco.initialGuess = 'linear';
+start = tic;
+bd_ld = coco.extract_FRC(freqrange);
+timings.cocoFRCbd_ld = toc(start)
+%% Comparison of FRC with/without nonlinear damping
+
+% plot SSM results
+FRCs = {FRC_LD_O5,FRC_ND_O5};
+legs = {'LD-SSM-$\mathcal{O}(5)$-unstable','ND-SSM-$\mathcal{O}(5)$-unstable'};
+legu = {'LD-SSM-$\mathcal{O}(5)$-stable','ND-SSM-$\mathcal{O}(5)$-stable'};
+color = {'r','b'};
+fig25 = figure;
+fig26 = figure;
+for k=1:2
+    FRC = FRCs{k};
+    SNidx = FRC.SNidx;
+    HBidx = FRC.HBidx;
+    FRC.st = double(FRC.st);
+    FRC.st(HBidx) = nan;
+    FRC.st(SNidx) = nan;
+    % color
+    ST = cell(2,1);
+    ST{1} = {[color{k},'--'],'LineWidth',1.5}; % unstable
+    ST{2} = {[color{k},'-'],'LineWidth',1.5};  % stable
+    figure(fig25); hold on
+    plot_stab_lines(FRC.om,FRC.Aout_frc(:,1),FRC.st,ST,legs{k},legu{k});
+    SNfig = plot(FRC.om(SNidx),FRC.Aout_frc(SNidx,1),thm.SN{:});
+    set(get(get(SNfig,'Annotation'),'LegendInformation'),...
+    'IconDisplayStyle','off');  
+    xlabel('$\Omega$','Interpreter','latex'); 
+    ylabel('$||u_1||_{\infty}$','Interpreter','latex'); 
+    set(gca,'FontSize',14);
+    grid on, axis tight; 
+    legend boxoff;
+    figure(fig26); hold on
+    plot_stab_lines(FRC.om,FRC.Aout_frc(:,2),FRC.st,ST,legs{k},legu{k});
+    SNfig = plot(FRC.om(SNidx),FRC.Aout_frc(SNidx,2),thm.SN{:});
+    set(get(get(SNfig,'Annotation'),'LegendInformation'),...
+    'IconDisplayStyle','off'); 
+    xlabel('$\Omega$','Interpreter','latex'); 
+    ylabel('$||u_2||_{\infty}$','Interpreter','latex'); 
+    set(gca,'FontSize',14);
+    grid on; axis tight; 
+    legend boxoff
+end
+
+% load coco solution
+legs = {'LD-Collocation-unstable','ND-Collocation-unstable'};
+legu = {'LD-Collocation-stable','ND-Collocation-stable'};
+
+bd = bd_nd{:};
+ndom   = coco_bd_col(bd,'omega'); 
+ndamp1 = coco_bd_col(bd, 'amp1');
+ndamp2 = coco_bd_col(bd, 'amp2');
+ndst   = coco_bd_col(bd, 'eigs');
+ndst = all(abs(ndst)<1,1);
+logs  = [1 2 3 4:3:numel(ndst)-4 numel(ndst)-1 numel(ndst)];
+ndamp1 = ndamp1(logs);
+ndamp2 = ndamp2(logs);
+ndom = ndom(logs);
+ndst = ndst(logs);
+figure(fig25); hold on
+plot(ndom(ndst), ndamp1(ndst), 'ro', 'MarkerSize', 6, 'LineWidth', 1,...
+    'DisplayName', 'ND-Collocation-stable'); 
+plot(ndom(~ndst), ndamp1(~ndst), 'ms', 'MarkerSize', 6, 'LineWidth', 1,...
+    'DisplayName', 'ND-Collocation-unstable');
+figure(fig26); hold on
+plot(ndom(ndst), ndamp2(ndst), 'ro', 'MarkerSize', 6, 'LineWidth', 1,...
+    'DisplayName', 'ND-Collocation-stable'); 
+plot(ndom(~ndst), ndamp2(~ndst), 'ms', 'MarkerSize', 6, 'LineWidth', 1,...
+    'DisplayName', 'ND-Collocation-unstable');
+
+bd = bd_ld{:};
+ldom   = coco_bd_col(bd,'omega'); 
+ldamp1 = coco_bd_col(bd, 'amp1');
+ldamp2 = coco_bd_col(bd, 'amp2');
+ldst   = coco_bd_col(bd, 'eigs');
+ldst = all(abs(ldst)<1,1);
+logs  = [1 2 3 4:3:numel(ldst)-4 numel(ldst)-1 numel(ldst)];
+ldamp1 = ldamp1(logs);
+ldamp2 = ldamp2(logs);
+ldom = ldom(logs);
+ldst = ldst(logs);
+figure(fig25); hold on
+plot(ldom(ldst), ldamp1(ldst), 'kd', 'MarkerSize', 6, 'LineWidth', 1,...
+    'DisplayName', 'LD-Collocation-stable'); 
+plot(ldom(~ldst), ldamp1(~ldst), 'gv', 'MarkerSize', 6, 'LineWidth', 1,...
+    'DisplayName', 'LD-Collocation-unstable');
+figure(fig26); hold on
+plot(ldom(ldst), ldamp2(ldst), 'kd', 'MarkerSize', 6, 'LineWidth', 1,...
+    'DisplayName', 'LD-Collocation-stable'); 
+plot(ldom(~ldst), ldamp2(~ldst), 'gv', 'MarkerSize', 6, 'LineWidth', 1,...
+    'DisplayName', 'LD-Collocation-unstable');
+%%
+timings

examples/AxialMovingBeam/build_model.m

@@ -0,0 +1,52 @@
+function [mass,damp,gyro,stiff,fnl,fext] = build_model(n,varargin)
+
+A = 0.04*0.03; % m^2
+I = 0.04*0.03^3/12;
+rho = 7680;
+E = 30e9;
+% eta = 1e-5*E;
+eta = 1e-4*E;
+L = 1;
+P = 67.5e3;
+
+kf = sqrt(E*I/(P*L^2));
+k1 = sqrt(E*A/P);
+alpha = I*eta/(L^3*sqrt(rho*A*P));
+gamma = 0.5128;
+mass  = eye(n);
+damp  = zeros(n);
+gyro  = zeros(n);
+stiff = zeros(n);
+fext  = zeros(n,1);
+cubic_coeff = zeros(n,n,n,n);
+for i=1:n
+    damp(i,i)  = alpha*(i*pi)^4;
+    stiff(i,i) = kf^2*(i*pi)^4-(gamma^2-1)*(i*pi)^2;
+    fext(i)    = (1-(-1)^i)/(i*pi);
+    for j=1:n
+        if j~=i
+            gyro(i,j) = 4*gamma*i*j/(i^2-j^2)*(1-(-1)^(i+j));
+        end
+        cubic_coeff(i,j,j,i) = i^2*j^2;
+    end
+end
+
+disp('the first four eigenvalues for undamped system');
+lamd = eigs([zeros(n) eye(n);-mass\stiff -mass\(gyro)],4,'smallestabs')
+
+if numel(varargin)>0 && strcmp(varargin{1}, 'nonlinear_damp')
+    % visoelastic damping
+    tmp  = sptensor(cubic_coeff);
+    subs1 = tmp.subs;
+    subs2 = subs1;
+    subs2(:,3) = subs1(:,3)+n;
+    subs = [subs1; subs2];
+    vals = [0.25*k1^2*pi^4*tmp.vals; 0.5*alpha*k1^2/kf^2*pi^4*tmp.vals];
+    f3 = sptensor(subs, vals, [n,2*n,2*n,2*n]); 
+else
+    % viscoelastic damping but ignore nonlinear part
+    f3 = 0.25*k1^2*pi^4*sptensor(cubic_coeff);
+end
+fnl = {[],f3};
+
+end

examples/AxialMovingBeam/cal_parameters.m

@@ -0,0 +1,64 @@
+function [alpha, omega, nu] = cal_parameters(N,l)
+
+syms x
+% N = 2; % number of modes
+% l = 2;
+
+% construct modal functions
+phis   = [];
+dphis  = [];
+ddphis = [];
+alphal = zeros(N,1);
+alphal(1:4) = [3.927 7.069 10.210 13.352]';
+alphal(5:N) = (4*(5:N)'+1)*pi/4;
+alpha  = alphal/l;
+R = sin(alphal)./sinh(alphal);
+E = (0.5*l*(1-R.^2)+(R.^2.*sinh(2*alphal)-sin(2*alphal))/(4*alpha)).^(-0.5);
+for n=1:N
+   phin   = E(n)*(sin(alpha(n)*x)-R(n)*sinh(alpha(n)*x));
+   dphin  = diff(phin,x);
+   ddphin = diff(dphin,x);
+   phis   = [phis; phin];
+   dphis  = [dphis; dphin];
+   ddphis = [ddphis; ddphin];
+   normphi = int(phin^2, x, 0, l);
+   fprintf('L2 norm of phi%d is %d\n', n, normphi);
+end
+
+% compute nonlinear coeffficients
+alpha1 = zeros(N);
+alpha2 = zeros(N);
+alpha  = zeros(N,N,N,N);
+
+for i=1:N
+    for j=1:N
+        alpha1(i,j) = int(phis(i)*ddphis(j),x,0,l);
+        alpha2(i,j) = int(dphis(i)*dphis(j),x,0,l);
+    end
+end
+
+
+for n=1:N
+    for m=1:N
+        for p=1:N
+            for q=1:N
+                alpha(n,m,p,q) = alpha1(n,q)*alpha2(m,p);
+            end
+        end
+    end
+end
+
+% compute natural frequencies
+omega = zeros(N,1);
+for i=1:N
+    int1 = int(phis(i)^2,x,0,l);
+    int2 = int(diff(phis(i),x,4)*phis(i),x,0,l);
+    omega(i) = sqrt(int2/int1);
+end
+
+nu = 1/(2*l);
+
+end
+
+
+