Commit 3180efba authored by Andreas Zilian's avatar Andreas Zilian
Browse files

Add examples for stability analysis.

parent 67a6a523
%% Example Beam1
% A horizontal beam in the plane composed of a single beam finite element
% subject to axial load
addpath([ '..' filesep '..' filesep '..' filesep ]);
%% Definition of sections
s1 = mafe.Section1D();
s1.E = 2.0e11; % [N/m^2]
s1.Iy = 1.0e-7; % [m^4] = 0.01m * (0.05m)^3 / 12
s1.A = 1.0e-4; % [m^2] = 0.01m * 0.01m
%% Definition of nodes in the system
% Node 1
n1 = mafe.Node2D( [0.0, 0.0] );
% Node 2
n2 = mafe.Node2D( [1.0, 0.0] );
% Vector of nodes
nodes = [ n1, n2 ];
%% Definition of elements in the system
% Element 1
e1 = mafe.Member2D( [n1, n2], s1 );
% vector of elements
elems = [ e1 ];
%% Definition of constraints in the system
% Constraint 1 and 2 (fixed vertically and horizontally at left end)
c1 = mafe.ConstraintNode( n1, mafe.DofType.Disp1, mafe.ConstraintType.Dirichlet, 0.0 );
c2 = mafe.ConstraintNode( n1, mafe.DofType.Disp2, mafe.ConstraintType.Dirichlet, 0.0 );
c3 = mafe.ConstraintNode( n1, mafe.DofType.Rota3, mafe.ConstraintType.Dirichlet, 0.0 );
% Constraint 3 (fixed vertically and horizontally at right end)
c4 = mafe.ConstraintNode( n2, mafe.DofType.Disp1, mafe.ConstraintType.Neumann, -1.0 );
c5 = mafe.ConstraintNode( n2, mafe.DofType.Disp2, mafe.ConstraintType.Dirichlet, 0.0 );
% vector of constraints
const = [ c1, c2, c3, c4, c5 ];
%% Setup of the finite element problem
fep = mafe.FeProblem( nodes, elems, const );
%% Select an analysis type: stability
ana = mafe.FeAnalysisStability(fep);
%% Calculate response
ana.analyse();
%% Extract some data from stability analysis
ana.L
ana.X
%% Visualise system response
cla; clf;
fig = figure(1); subplot('Position',[0.05 0.05 0.90 0.90]); axis equal; hold on;
fep.plotSystem('reference');
ana.putModeToDof(1)
fep.plotSystem('deformed');
fep.printSystem();
%% Example Beam2
% A horizontal beam in the plane composed of multiple (N) beam finite element
% subject to axial load
addpath([ '..' filesep '..' filesep '..' filesep ]);
%% General properties
L = 1.0; % length of the cantilever
N = 10; % number of finite elements
%% Definition of sections
s1 = mafe.Section1D();
s1.E = 2.0e11; % [N/m^2]
s1.Iy = 1.0e-7; % [m^4] = 0.01m * (0.05m)^3 / 12
s1.A = 1.0e-4; % [m^2] = 0.01m * 0.01m
%% Definition of nodes in the system
nodes = mafe.Node.empty(0,0);
for ii = 1:N+1
% create node
nodes(end+1) = mafe.Node2D( [L/N*(ii-1), 0.0] );
end
%% Definition of elements in the system
elems = mafe.Element.empty(0,0);
for ii = 1:N
% create element
elems(end+1) = mafe.Member2D( [nodes(ii), nodes(ii+1)], s1 );
end
%% Definition of constraints in the system
const = mafe.Constraint.empty(0,0);
% Constraint: kinematic, clamped at the left
const(end+1) = mafe.ConstraintNode( nodes(1), mafe.DofType.Disp1, mafe.ConstraintType.Dirichlet, 0.0 );
const(end+1) = mafe.ConstraintNode( nodes(1), mafe.DofType.Disp2, mafe.ConstraintType.Dirichlet, 0.0 );
const(end+1) = mafe.ConstraintNode( nodes(1), mafe.DofType.Rota3, mafe.ConstraintType.Dirichlet, 0.0 );
% Constraint: kinematic, pinned at the right
const(end+1) = mafe.ConstraintNode( nodes(end), mafe.DofType.Disp2, mafe.ConstraintType.Dirichlet, 0.0 );
% Constraint: static, axial force at the right
const(end+1) = mafe.ConstraintNode( nodes(end), mafe.DofType.Disp1, mafe.ConstraintType.Neumann ,-1.0 );
%% Setup of the finite element problem
fep = mafe.FeProblem( nodes, elems, const );
%% Select an analysis type: stability
ana = mafe.FeAnalysisStability(fep);
%% Calculate response
ana.analyse();
%% Extract some data from stability analysis
%ana.L
%ana.X
%% Visualise system response
cla; clf;
fig = figure(1); subplot('Position',[0.05 0.05 0.90 0.90]); axis equal; hold on;
fep.plotSystem('reference');
ana.putModeToDof(1)
fep.plotSystem('deformed');
fep.printSystem();
%% Example Frame1
% A half-open frame structure with axial forces in the left colum and hinge
% at right support
addpath([ '..' filesep '..' filesep '..' filesep ]);
%% Definition of sections
s1 = mafe.Section1D();
s1.E = 2.0e11; % [N/m^2]
s1.Iy = 1.0e-7; % [m^4] = 0.01m * (0.05m)^3 / 12
s1.A = 1.0e-4; % [m^2] = 0.01m * 0.01m
%% Definition of member loads
l1 = mafe.EleLoad( mafe.DofType.Disp2, -1.0e5 ); % [N/m]
%% Definition of nodes in the system
% Node 1
n1 = mafe.Node2D( [0.0, 0.0] );
% Node 2
n2 = mafe.Node2D( [0.0, 1.0] );
% Node 3
n3 = mafe.Node2D( [1.0, 1.0] );
% Node 4
n4 = mafe.Node2D( [1.0, 0.0] );
% 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, l1 );
% Element 3
e3 = mafe.Member2D( [n3, n4], s1 );
% vector of elements
elems = [ e1, e2, e3 ];
%% Definition of constraints in the system
% Constraints: fully clamped at left support
c1 = mafe.ConstraintNode( n1, mafe.DofType.Disp1, mafe.ConstraintType.Dirichlet, 0.0 );
c2 = mafe.ConstraintNode( n1, mafe.DofType.Disp2, mafe.ConstraintType.Dirichlet, 0.0 );
c3 = mafe.ConstraintNode( n1, mafe.DofType.Rota3, mafe.ConstraintType.Dirichlet, 0.0 );
% Constraint: vertical force in left column, applied at top
c4 = mafe.ConstraintNode( n2, mafe.DofType.Disp2, mafe.ConstraintType.Neumann, -0.0 );
% Constarint: pin supported at the right
c5 = mafe.ConstraintNode( n4, mafe.DofType.Disp1, mafe.ConstraintType.Dirichlet, 0.0 );
c6 = mafe.ConstraintNode( n4, mafe.DofType.Disp2, mafe.ConstraintType.Dirichlet, 0.0 );
% vector of constraints
const = [ c1, c2, c3, c4, c5, c6 ];
%% Setup of the finite element problem
fep = mafe.FeProblem( nodes, elems, const );
%% Select an analysis type: stability
ana = mafe.FeAnalysisStability(fep);
%% Calculate response
ana.analyse();
%% Extract some data from stability analysis
ana.L % eigenvalues = buckling load factors
ana.X % eigenvectors = buckling shapes
%% Visualise system response
cla; clf;
fig = figure(1); subplot('Position',[0.05 0.05 0.90 0.90]); axis equal; hold on;
fep.plotSystem('reference');
%ana.putModeToDof(3)
fep.plotSystem('deformed');
ana.putModeToDof(1)
fep.plotSystem('deformed');
fep.printSystem();
%% Example Frame2
% Framework bridge in the plane
addpath([ '..' filesep '..' filesep '..' filesep ]);
%% Definition of sections
s1 = mafe.Section1D();
s1.E = 2.0e11; % [N/m^2]
s1.Iy = 1.7e-5; % [m^4] = 0.02m * (0.10m)^3 / 12
s1.A = 2.0e-3; % [m^2] = 0.02m * 0.10m
%% Definition of member loads
l1 = mafe.EleLoad( mafe.DofType.Disp2, -1.0e6 ); % [N/m]
%% Definition of nodes in the system
n01 = mafe.Node2D( [ 0.0, 0.0] );
n02 = mafe.Node2D( [ 2.0, 0.0] );
n03 = mafe.Node2D( [ 4.0, 0.0] );
n04 = mafe.Node2D( [ 6.0, 0.0] );
n05 = mafe.Node2D( [ 8.0, 0.0] );
n06 = mafe.Node2D( [10.0, 0.0] );
n07 = mafe.Node2D( [12.0, 0.0] );
n08 = mafe.Node2D( [ 2.0, 2.0] );
n09 = mafe.Node2D( [ 4.0, 2.0] );
n10 = mafe.Node2D( [ 6.0, 2.0] );
n11 = mafe.Node2D( [ 8.0, 2.0] );
n12 = mafe.Node2D( [10.0, 2.0] );
% Vector of nodes
nodes = [ n01, n02, n03, n04, n05, n06, n07, n08, n09, n10, n11, n12 ];
%% Definition of elements in the system
e01 = mafe.Member2D( [n01, n02], s1, l1 );
e02 = mafe.Member2D( [n02, n03], s1, l1 );
e03 = mafe.Member2D( [n03, n04], s1, l1 );
e04 = mafe.Member2D( [n04, n05], s1, l1 );
e05 = mafe.Member2D( [n05, n06], s1, l1 );
e06 = mafe.Member2D( [n06, n07], s1, l1 );
e07 = mafe.Member2D( [n08, n09], s1 );
e08 = mafe.Member2D( [n09, n10], s1 );
e09 = mafe.Member2D( [n10, n11], s1 );
e10 = mafe.Member2D( [n11, n12], s1 );
e11 = mafe.Member2D( [n02, n08], s1 );
e12 = mafe.Member2D( [n03, n09], s1 );
e13 = mafe.Member2D( [n04, n10], s1 );
e14 = mafe.Member2D( [n05, n11], s1 );
e15 = mafe.Member2D( [n06, n12], s1 );
e16 = mafe.Member2D( [n01, n08], s1 );
e17 = mafe.Member2D( [n02, n09], s1 );
e18 = mafe.Member2D( [n03, n10], s1 );
e19 = mafe.Member2D( [n05, n10], s1 );
e20 = mafe.Member2D( [n06, n11], s1 );
e21 = mafe.Member2D( [n07, n12], s1 );
% vector of elements
elems = [ e01, e02, e03, e04, e05, e06, e07, e08, e09, e10, ...
e11, e12, e13, e14, e15, e16, e17, e18, e19, e20, e21 ];
%% Definition of constraints in the system
% Constraint 1 and 2 (fixed vertically and horizontally at node 1)
c1 = mafe.ConstraintNode( n01, mafe.DofType.Disp1, mafe.ConstraintType.Dirichlet, 0.0 );
c2 = mafe.ConstraintNode( n01, mafe.DofType.Disp2, mafe.ConstraintType.Dirichlet, 0.0 );
% Constraint 3 and 4 (fixed vertically and horizontally at node 7)
c3 = mafe.ConstraintNode( n07, mafe.DofType.Disp1, mafe.ConstraintType.Dirichlet, 0.0 );
c4 = mafe.ConstraintNode( n07, mafe.DofType.Disp2, mafe.ConstraintType.Dirichlet, 0.0 );
% vector of constraints
const = [ c1, c2, c3, c4 ];
%% Setup of the finite element problem
fep = mafe.FeProblem( nodes, elems, const );
%% Select an analysis type: static
ana = mafe.FeAnalysisStability(fep);
%% Calculate response
ana.analyse();
%% Extract some data from stability analysis
ana.L(1) % eigenvalues = buckling load factors
ana.X; % eigenvectors = buckling shapes
%% Visualise system response
cla; clf;
fig = figure(1); subplot('Position',[0.05 0.05 0.90 0.90]); axis equal; hold on;
fep.plotSystem('reference');
fep.plotSystem('deformed');
ana.putModeToDof(1)
fep.plotSystem('deformed');
fep.printSystem();
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