Commit 97bebbd4 authored by Andreas Zilian's avatar Andreas Zilian
Browse files

Provide step-bystep solution to tutorial2 in dynamics.

parent 3a43d9b9
%% Tutorial 2: Structural Dynamics
% A frame system with point mass and spring - dynamic analysis (free and forced)
%
% step 1 (task 2.1):
%
% definition of time functions
%
addpath([ '..' filesep '..' filesep '..' filesep ]);
% -----------------------------------------------------------DEFINE FE ITEMS---
%% Definition of time functions (used to describe time-dependent values)
tF = mafe.TimeFunction( 2.5, 1.0, 0.0 ); % force: Omega = 2.5 [rad/s], a_c = 1.0, a_s = 0.0
tU = mafe.TimeFunction( 5.0, 0.0, 1.0 ); % base : Omega = 5.0 [rad/s], a_c = 0.0, a_s = 1.0
tfuns = [tF, tU];
%% Tutorial 2: Structural Dynamics
% A frame system with point mass and spring - dynamic analysis (free and forced)
%
% step 2 (task 3.1-3.4):
%
% solve the free vibration problem (no forcing) - inspection of modes
%
addpath([ '..' filesep '..' filesep '..' filesep ]);
% -----------------------------------------------------------DEFINE FE ITEMS---
%% Definition of time functions (used to describe time-dependent values)
tF = mafe.TimeFunction( 2.5, 1.0, 0.0 ); % force: Omega = 2.5 [rad/s], a_c = 1.0, a_s = 0.0
tU = mafe.TimeFunction( 5.0, 0.0, 1.0 ); % base : Omega = 5.0 [rad/s], a_c = 0.0, a_s = 1.0
tfuns = [tF, tU];
%% Definition of sections (containing geometrical data and material data)
s1 = mafe.Section1D(); % for the members
s1.A = 5.0e-4; % [m^2] = 0.01m * 0.05m
s1.Iy = 1.0e-7; % [m^4] = 0.01m * (0.05m)^3 / 12
s1.E = 2.0e11; % [N/m^2]
s1.rho = 8.0e+3; % [kg/m^3]
s2 = mafe.Section0D(); % for the point mass (m, but no k)
s2.m = 50.0; % [kg]
s2.k = 0.00; % [N]
s3 = mafe.Section0D(); % for the point spring (k, but no m)
s3.m = 0.00; % [kg]
s3.k = 1000; % [N]
%% Definition of nodes in the system
% Node 1
n1 = mafe.Node2D( [0.0, 0.0] );
% Node 2
n2 = mafe.Node2D( [0.0, 2.5] );
% Node 3
n3 = mafe.Node2D( [2.0, 2.5] );
% Node 4
n4 = mafe.Node2D( [4.0, 2.5] );
% Vector of nodes
nodes = [ n1, n2, n3, n4 ];
%% Definition of elements in the system
% Element 1
e1 = mafe.Member2D( [n1, n2], s1 );
% Element 2
e2 = mafe.Member2D( [n2, n3], s1 );
% Element 3
e3 = mafe.Member2D( [n3, n4], s1 );
% Element 4
e4 = mafe.Point( n3, s2, mafe.DofType.Disp1 );
% Element 5
e5 = mafe.Point( n3, s2, mafe.DofType.Disp2 );
% Element 6
e6 = mafe.Point( n4, s3, mafe.DofType.Disp1 );
% vector of elements
elems = [ e1, e2, e3, e4, e5, e6 ];
%% Definition of constraints in the system
% Constraint 1 and 2 (fixed vertically and horizontally at lower left base point)
c1 = mafe.ConstraintNode( n1, mafe.DofType.Disp1, mafe.ConstraintType.Dirichlet, 0.0 );
c2 = mafe.ConstraintNode( n1, mafe.DofType.Disp2, mafe.ConstraintType.Dirichlet, 0.0 );
% Constraint 3 (fixed vertically upper right end)
c3 = mafe.ConstraintNode( n4, mafe.DofType.Disp2, mafe.ConstraintType.Dirichlet, 0.0 );
% Constraint 4 (horizontal point force at the corner)
c4 = mafe.ConstraintNode( n2, mafe.DofType.Disp1, mafe.ConstraintType.Neumann, 0. ); % 0 [N]
% vector of constraints
const = [ c1, c2, c3, c4 ];
% --------------------------------------------DEFINE FE PROBLEM AND ANALYSIS---
%% Setup of the finite element problem
fep = mafe.FeProblem( nodes, elems, const );
%% Select an analysis type: dynamic
ana = mafe.FeAnalysisDynamicFD( fep );
%% Calculate response
ana.analyse();
%% Visualise system response (individual modes of vibration = eigenvectors)
cla; clf;
fig = figure(1); subplot('Position',[0.05 0.05 0.90 0.90]); axis equal; hold on;
xlim([-1 5])
ylim([-1 3])
fep.plotSystem('reference'); % plot the undeformed system
ana.putModeToDof(1, 'real'); % select the first mode (complex-conjugate eigenstates)
fep.plotSystem('deformed'); % plot the mode
ana.putModeToDof(3, 'real'); % select the third mode (complex-conjugate eigenstates)
fep.plotSystem('deformed'); % plot the mode
%% Extract results for eigenvalues and eigenvectors
eigenvalues = ana.L([1 3])
eigenvectors = ana.X(:,[1 3])
% get the first two (distinct) eigenfrequencies
angular_frequency_in_rad_per_second = abs(imag(eigenvalues)) % omega
natural_frequency_in_hertz = abs(imag(eigenvalues)) / (2*pi) % f
%% Tutorial 2: Structural Dynamics
% A frame system with point mass and spring - dynamic analysis (free and forced)
%
% step 3 (task 3.5-3.7):
%
% solve the free vibration problem (no forcing) - inspection of homogeneous sol
%
addpath([ '..' filesep '..' filesep '..' filesep ]);
% -----------------------------------------------------------DEFINE FE ITEMS---
%% Definition of time functions (used to describe time-dependent values)
tF = mafe.TimeFunction( 2.5, 1.0, 0.0 ); % force: Omega = 2.5 [rad/s], a_c = 1.0, a_s = 0.0
tU = mafe.TimeFunction( 5.0, 0.0, 1.0 ); % base : Omega = 5.0 [rad/s], a_c = 0.0, a_s = 1.0
tfuns = [tF, tU];
%% Definition of sections (containing geometrical data and material data)
s1 = mafe.Section1D(); % for the members
s1.A = 5.0e-4; % [m^2] = 0.01m * 0.05m
s1.Iy = 1.0e-7; % [m^4] = 0.01m * (0.05m)^3 / 12
s1.E = 2.0e11; % [N/m^2]
s1.rho = 8.0e+3; % [kg/m^3]
s2 = mafe.Section0D(); % for the point mass (m, but no k)
s2.m = 50.0; % [kg]
s2.k = 0.00; % [N]
s3 = mafe.Section0D(); % for the point spring (k, but no m)
s3.m = 0.00; % [kg]
s3.k = 1000; % [N]
%% Definition of nodes in the system
% Node 1
n1 = mafe.Node2D( [0.0, 0.0] );
% Node 2
n2 = mafe.Node2D( [0.0, 2.5] );
% Node 3
n3 = mafe.Node2D( [2.0, 2.5] );
% Node 4
n4 = mafe.Node2D( [4.0, 2.5] );
% Vector of nodes
nodes = [ n1, n2, n3, n4 ];
%% Definition of elements in the system
% Element 1
e1 = mafe.Member2D( [n1, n2], s1 );
% Element 2
e2 = mafe.Member2D( [n2, n3], s1 );
% Element 3
e3 = mafe.Member2D( [n3, n4], s1 );
% Element 4
e4 = mafe.Point( n3, s2, mafe.DofType.Disp1 );
% Element 5
e5 = mafe.Point( n3, s2, mafe.DofType.Disp2 );
% Element 6
e6 = mafe.Point( n4, s3, mafe.DofType.Disp1 );
% vector of elements
elems = [ e1, e2, e3, e4, e5, e6 ];
%% Definition of constraints in the system
% Constraint 1 and 2 (fixed vertically and horizontally at lower left base point)
c1 = mafe.ConstraintNode( n1, mafe.DofType.Disp1, mafe.ConstraintType.Dirichlet, 0.0 );
c2 = mafe.ConstraintNode( n1, mafe.DofType.Disp2, mafe.ConstraintType.Dirichlet, 0.0 );
% Constraint 3 (fixed vertically upper right end)
c3 = mafe.ConstraintNode( n4, mafe.DofType.Disp2, mafe.ConstraintType.Dirichlet, 0.0 );
% Constraint 4 (horizontal point force at the corner)
c4 = mafe.ConstraintNode( n2, mafe.DofType.Disp1, mafe.ConstraintType.Neumann, 0. ); % 0 [N]
% vector of constraints
const = [ c1, c2, c3, c4 ];
% --------------------------------------------DEFINE FE PROBLEM AND ANALYSIS---
%% Setup of the finite element problem
fep = mafe.FeProblem( nodes, elems, const );
%% Select an analysis type: dynamic
ana = mafe.FeAnalysisDynamicFD( fep );
%% Calculate response
ana.analyse();
%% Visualise system response (homogeneous solution for given initial conditions)
cla; clf;
fig = figure(1); subplot('Position',[0.05 0.05 0.90 0.90]); axis equal; hold on;
% apply initial conditions
modeshape = real(ana.X(:,1:2:fep.ndofs)); % get some selected (or all) modes
u0 = 0.00 * modeshape(:,1) + 1.00 * modeshape(:,2);
v0 = 0.00 * modeshape(:,1) + 0.00 * modeshape(:,2);
ana.applyInitialConditions( u0, v0 );
% evaluate and plot for given time instants
% first undamped eigenfrequency is omega = ~5.9 [rad/s] -> T = 2pi/omega = 1.1 [s]
for tt = 0:0.02:1.1*5
ana.putSolToDof(tt);
clf; subplot('Position',[0.05 0.05 0.90 0.90]); axis equal; hold on;
xlim([-1 5])
ylim([-1 3])
fep.plotSystem('reference');
fep.plotSystem('deformed');
drawnow;
end
%% Tutorial 2: Structural Dynamics
% A frame system with point mass and spring - dynamic analysis (free and forced)
%
% step 4 (task 4.1-4.6):
%
% solve the forced vibration problem (with zero ic) - inspection of particular sol
%
addpath([ '..' filesep '..' filesep '..' filesep ]);
% -----------------------------------------------------------DEFINE FE ITEMS---
%% Definition of time functions (used to describe time-dependent values)
tF = mafe.TimeFunction( 2.5, 1.0, 0.0 ); % force: Omega = 2.5 [rad/s], a_c = 1.0, a_s = 0.0
tU = mafe.TimeFunction( 5.0, 0.0, 1.0 ); % base : Omega = 5.0 [rad/s], a_c = 0.0, a_s = 1.0
tfuns = [tF, tU];
%% Definition of sections (containing geometrical data and material data)
s1 = mafe.Section1D(); % for the members
s1.A = 5.0e-4; % [m^2] = 0.01m * 0.05m
s1.Iy = 1.0e-7; % [m^4] = 0.01m * (0.05m)^3 / 12
s1.E = 2.0e11; % [N/m^2]
s1.rho = 8.0e+3; % [kg/m^3]
s2 = mafe.Section0D(); % for the point mass (m, but no k)
s2.m = 50.0; % [kg]
s2.k = 0.00; % [N]
s3 = mafe.Section0D(); % for the point spring (k, but no m)
s3.m = 0.00; % [kg]
s3.k = 1000; % [N]
%% Definition of nodes in the system
% Node 1
n1 = mafe.Node2D( [0.0, 0.0] );
% Node 2
n2 = mafe.Node2D( [0.0, 2.5] );
% Node 3
n3 = mafe.Node2D( [2.0, 2.5] );
% Node 4
n4 = mafe.Node2D( [4.0, 2.5] );
% Vector of nodes
nodes = [ n1, n2, n3, n4 ];
%% Definition of elements in the system
% Element 1
e1 = mafe.Member2D( [n1, n2], s1 );
% Element 2
e2 = mafe.Member2D( [n2, n3], s1 );
% Element 3
e3 = mafe.Member2D( [n3, n4], s1 );
% Element 4
e4 = mafe.Point( n3, s2, mafe.DofType.Disp1 );
% Element 5
e5 = mafe.Point( n3, s2, mafe.DofType.Disp2 );
% Element 6
e6 = mafe.Point( n4, s3, mafe.DofType.Disp1 );
% vector of elements
elems = [ e1, e2, e3, e4, e5, e6 ];
%% Definition of constraints in the system
% Constraint 1 and 2 (fixed vertically and horizontally at lower left base point)
c1 = mafe.ConstraintNode( n1, mafe.DofType.Disp1, mafe.ConstraintType.Dirichlet, 0.1, tU ); % 0.1 [m]
c2 = mafe.ConstraintNode( n1, mafe.DofType.Disp2, mafe.ConstraintType.Dirichlet, 0.0 );
% Constraint 3 (fixed vertically upper right end)
c3 = mafe.ConstraintNode( n4, mafe.DofType.Disp2, mafe.ConstraintType.Dirichlet, 0.0 );
% Constraint 4 (horizontal point force at the corner)
c4 = mafe.ConstraintNode( n2, mafe.DofType.Disp1, mafe.ConstraintType.Neumann, 1000., tF ); % 1000 [N]
% vector of constraints
const = [ c1, c2, c3, c4 ];
% --------------------------------------------DEFINE FE PROBLEM AND ANALYSIS---
%% Setup of the finite element problem
fep = mafe.FeProblem( nodes, elems, const );
%% Select an analysis type: dynamic
ana = mafe.FeAnalysisDynamicFD( fep, tfuns ); % analysis requires time functions!
%% Calculate response
ana.analyse();
%% Visualise system response (homogeneous solution for given initial conditions)
cla; clf;
fig = figure(1); subplot('Position',[0.05 0.05 0.90 0.90]); axis equal; hold on;
% apply (zero) initial conditions
u0 = zeros(fep.ndofs,1);
v0 = zeros(fep.ndofs,1);
ana.applyInitialConditions( u0, v0 );
% evaluate and plot for given time instants
% first undamped eigenfrequency is omega = ~5.9 [rad/s] -> T = 2pi/omega = 1.1 [s]
for tt = 0:0.02:1.1*5
ana.putSolToDof(tt);
clf; subplot('Position',[0.05 0.05 0.90 0.90]); axis equal; hold on;
xlim([-1 5])
ylim([-1 3])
fep.plotSystem('reference');
fep.plotSystem('deformed');
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