Commit 3a43d9b9 authored by Andreas Zilian's avatar Andreas Zilian
Browse files

Provide step-bystep solution to tutorial1 in dynamics.

parent 58289564
%% Tutorial 1: Structural Dynamics
% A frame system with point mass and spring - static and dynamic analysis
%
% step 1 (task 2.1-2.2):
%
% definition of cross section and material properties
%
addpath([ '..' filesep '..' filesep '..' filesep ]);
% -----------------------------------------------------------DEFINE FE ITEMS---
%% 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]
% show the objects created above and the variables they hold
s1
s2
%% Tutorial 1: Structural Dynamics
% A frame system with point mass and spring - static and dynamic analysis
%
% step 2 (task 2.3):
%
% definition of nodes and beam/point elements
%
addpath([ '..' filesep '..' filesep '..' filesep ]);
% -----------------------------------------------------------DEFINE FE ITEMS---
%% 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 ];
%% Tutorial 1: Structural Dynamics
% A frame system with point mass and spring - static and dynamic analysis
%
% step 3 (task 2.4):
%
% definition of constraints (kinematic and static/dynamic)
%
addpath([ '..' filesep '..' filesep '..' filesep ]);
% -----------------------------------------------------------DEFINE FE ITEMS---
%% 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, 1000. ); % 1000 [N]
% vector of constraints
const = [ c1, c2, c3, c4 ];
%% Tutorial 1: Structural Dynamics
% A frame system with point mass and spring - static and dynamic analysis
%
% step 4 (task 3.1):
%
% create finite element problem (and plot the undeformed configuration)
%
addpath([ '..' filesep '..' filesep '..' filesep ]);
% -----------------------------------------------------------DEFINE FE ITEMS---
%% 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, 1000. ); % 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 );
% visualise the system
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');
% print the current state of the system
fep.printSystem();
%% Tutorial 1: Structural Dynamics
% A frame system with point mass and spring - static and dynamic analysis
%
% step 5 (task 3.2-3.6):
%
% create and perform a static finite element analysis (and plot result)
%
addpath([ '..' filesep '..' filesep '..' filesep ]);
% -----------------------------------------------------------DEFINE FE ITEMS---
%% 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, 1000. ); % 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: static
ana = mafe.FeAnalysisStatic(fep);
%% Calculate response
ana.analyse();
%% Visualise system response
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');
fep.plotSystem('deformed');
fep.printSystem();
%% Tutorial 1: Structural Dynamics
% A frame system with point mass and spring - static and dynamic analysis
%
% step 6 (task 4.1-4.5):
%
% create and perform a dynamic finite element analysis (and extract results)
%
addpath([ '..' filesep '..' filesep '..' filesep ]);
% -----------------------------------------------------------DEFINE FE ITEMS---
%% 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, 1000. ); % 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 );
%% 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
ana.putModeToDof(5, 'real'); % select the fifth mode (complex-conjugate eigenstates)
fep.plotSystem('deformed'); % plot the mode
%% Extract results for eigenvalues and eigenvectors
eigenvalues = ana.L([1 3 5])
eigenvectors = ana.X(:,[1 3 5])
% get the first three (distinct) eigenfrequencies
angular_frequency_in_rad_per_second = abs(imag(eigenvalues)) % omega
natural_frequency_in_hertz = abs(imag(eigenvalues)) / (2*pi) % f
%% Tutorial 1: Structural Dynamics
% A frame system with point mass and spring - static and dynamic analysis
%
% step 7 (task 5.1-5.3):
%
% provide initial conditions and plot system response at certain time instances
%
addpath([ '..' filesep '..' filesep '..' filesep ]);
% -----------------------------------------------------------DEFINE FE ITEMS---
%% 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, 1000. ); % 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 );
%% Calculate response
ana.analyse();
%% Visualise system response (overall 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
% % The same as above, but with a video being generated (*.avi format)
% video = VideoWriter( mfilename('fullpath') );
% open(video);
% 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');
% writeVideo(video, getframe);
% drawnow;
% end
% close(video);
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