Commit fba8a3b8 authored by Andreas Zilian's avatar Andreas Zilian
Browse files

Add example to demonstrate base excitation with multiple TimeFunctions.

parent 2b9f4a46
%% Example Kelvin-Voigt mass system
% Two masses connected by rheological devices, one mass constrained by Dirichlet
addpath([ '..' filesep '..' filesep '..' filesep ]);
%% Definition of sections
s1 = mafe.Section0D();
s1.k = 4.0; % [N/m]
s1.d = 0.1; % [Ns/m]
s2 = mafe.Section0D();
s2.m = 1.0; % [Ns^2/m]
%% Definition of time functions: Omega, factor_cosine, factor_sine
tfuns = mafe.TimeFunction.empty(0,0);
% define a spectrum of excitation
K = 80;
% chosen frequencies
Omegas = linspace(0,8,K);
% chosen amplitudes
% Amplit = logspace(-2,-6,K); % == 'push over'
Amplit = 0.5 * [Omegas == 0.5 + min(abs(Omegas-0.5))]; % pronounce single frequency sharply == 'harmonic'
% Amplit = 1.e-3 * Omegas./((Omegas-0.5).^2 + 0.005); % pronounce single frequency smoothly == 'harmonic'
% Amplit = 1.e-3 * Omegas./(((Omegas-1.0).*(Omegas-3.0)).^2 + 0.05); % pronounce certain frequencies == 'harmonic'
% set time functions
for ii = 1:K
tfuns(end+1) = mafe.TimeFunction( Omegas(ii), 0.0, Amplit(ii) ); % time function = A*sin(Omega*t)
end
%% Definition of nodes in the system
% Node 1
n1 = mafe.Node1D( [0.5] );
% Node 2
n2 = mafe.Node1D( [1.5] );
% Vector of nodes
nodes = [ n1, n2 ];
%% Definition of elements in the system
% Element 1
e1 = mafe.RheoKelvinVoigt1D( [n1, n2], s1 );
% Element 2
e2 = mafe.Point( n2, s2 );
% vector of elements
elems = [ e1, e2 ];
%% Definition of constraints in the system
const = mafe.Constraint.empty(0,0);
for tf = tfuns
const(end+1) = mafe.ConstraintNode( n1, mafe.DofType.Disp1, mafe.ConstraintType.Dirichlet, 1.0, tf ); % [m]
end
%% Setup of the finite element problem
fep = mafe.FeProblem( nodes, elems, const );
%% Select an analysis type: dynamic
ana = mafe.FeAnalysisDynamicFD( fep, tfuns );
%% Calculate response
ana.analyse();
%% Obtain amplification matrix
V = ana.Z(:,2:end);
%% --- Inspect results
%% Show
% - excitation spectrum of the prescibed base displacement u(t) at n1
% - response spectrum of the resulting displacement u(t) at n2
fig = figure(1); clf; subplot('Position',[0.10 0.10 0.80 0.85]); hold on;
plot(Omegas,[Amplit; abs(Amplit.*V)], '.-');
xlim([ 0 8]);
grid;
legend('|A_k| of u_{n1}(t)', '|A_k| of u_{n2}(t)');
xlabel('excitation frequency \Omega_k [rad/s]');
ylabel('amplitude of displacement u_k = A_k * sin(\Omega t) [m]');
%% Show relative amplifications due to forced vibration
fig = figure(2); clf; subplot('Position',[0.10 0.10 0.80 0.85]); hold on;
plot(Omegas,abs(V), '.');
xlim([ 0 8]);
ylim([-4 4]);
grid;
legend('u_{n1} = prescribed', 'u_{n2} = resulting');
xlabel('excitation frequency \Omega [rad/s]');
ylabel('amplification of natural coordinate (displacement) [-]');
% Show the system response (total solution for given initial conditions)
fig = figure(3); clf; subplot('Position',[0.10 0.10 0.80 0.85]); hold on;
% apply initial conditions
u0 = zeros(fep.ndofs,1);
v0 = zeros(fep.ndofs,1);
ana.applyInitialConditions( u0, v0 );
% evaluate and plot for given time instants
u_h = [];
u_p = [];
time = 0:0.1:pi*50; % first undamped eigenfrequency is omega = 2 [rad/s] -> T = 2pi/omega = 3.14 [s]
for tt = time
ana.putSolToDof(tt, 'homogeneous');
u_h(:,end+1) = [n1.lhs( find( n1.dof == mafe.DofType.Disp1 ) ), n2.lhs( find( n2.dof == mafe.DofType.Disp1 ) )];
ana.putSolToDof(tt, 'particular');
u_p(:,end+1) = [n1.lhs( find( n1.dof == mafe.DofType.Disp1 ) ), n2.lhs( find( n2.dof == mafe.DofType.Disp1 ) )];
end
plot(time,[u_p; u_h+u_p], '-');
xlim([ time(1) time(end)]);
grid;
legend('u_p(t) at n1 = prescribed', 'u_p(t) at n2 = resulting', 'u(t) at n1 = prescribed', 'u(t) at n2 = resulting');
xlabel('time [s]');
ylabel('natural coordinate (displacement) [m]');
% % Visualise system response (total solution for given initial conditions)
% fig = figure(4); clf; subplot('Position',[0.05 0.05 0.90 0.90]); axis equal; hold on;
% % apply initial conditions
% u0 = zeros(fep.ndofs,1);
% v0 = zeros(fep.ndofs,1);
% ana.applyInitialConditions( u0, v0 );
% % evaluate and plot for given time instants
% time = 0:0.1:pi*50; % first undamped eigenfrequency is omega = 2 [rad/s] -> T = 2pi/omega = 3.14 [s]
% for tt = time
% ana.putSolToDof(tt);
% fig = figure(4); clf; hold on;
% xlim([0 2.5])
% ylim([-0.25 0.25])
% fep.plotSystem('reference');
% fep.plotSystem('deformed');
% text(2.0,0.2,sprintf('t = %2.3f s', tt));
% drawnow;
% end
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment