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

Add examples for stability analysis.

parent 3180efba
%% Example Flat Shell Stability 1b
% example on how to use the integrated mesh generator for
% A two-dimensional quadratic plate structure composed of distorted elements
% -> extended patch test (using many distorted quadrilateral elements)
% subject to in-plane load, stability analysis
addpath([ '..' filesep '..' filesep '..' filesep ]);
%% Definition of sections
s1 = mafe.Section2D();
s1.E = 2.1e6; % [N/m^2]
s1.nue = 0.3; % [-]
s1.t = 1.0e-3; % [m]
% -----------------------------------------------------------MESH GENERATION---
%% Describe the goemetry of the problem for the mesh generator
% function that controls grid node distribution along geo lines
nf1 = mgen.MgNormFunction( mgen.MgNormFunctionType.Pow , [0.0, 0.5] );
% create all geo points
gp1 = mgen.MgGeoPoint( [ 0.0, 0.0] );
gp2 = mgen.MgGeoPoint( [ 1.0, 0.0] );
gp3 = mgen.MgGeoPoint( [ 1.0, 1.0] );
gp4 = mgen.MgGeoPoint( [ 0.0, 1.0] );
% create all geo lines
n = 12;
gl1 = mgen.MgGeoLine2Point( [gp2, gp1], n, nf1 ); % ask for 4 grid nodes along this line
gl2 = mgen.MgGeoLine2Point( [gp3, gp2], n, nf1 ); % ask for 6 grid nodes along this line
gl3 = mgen.MgGeoLine2Point( [gp4, gp3], n, nf1 ); % this will adapt automatically
gl4 = mgen.MgGeoLine2Point( [gp1, gp4], n, nf1 ); % this will adapt automatically
% create all geo areas
ga1 = mgen.MgGeoArea4Line( [gl1, gl2, gl3, gl4] );
% create geo object and generate the mesh
geo = mgen.MgGeo( [gp1, gp2, gp3, gp4], [gl1, gl2, gl3, gl4], [ga1] );
% print geo and mesh data info
%fprintf(geo.print());
% -----------------------------------------------------------DEFINE FE ITEMS---
%% Definition of nodes in the system
nodes = mafe.Node.empty(0,0);
for gn = geo.getAllGridNodes()
nodes(end+1) = mafe.Node2D( gn.coord );
end
%% Definition of elements in the system
elems = mafe.Element.empty(0,0);
for ge = geo.getAllGridElems()
elems(end+1) = mafe.FlatShell2D( nodes( [ ge.nodes.id ] ), s1 );
end
%% Definition of constraints in the system
% vector of constraints
const = mafe.Constraint.empty(0,0);
% membrane effects
nn = nodes( [gl4.getGridNodes.id] ); % get nodes on GeoLine 4 (left)
const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Disp1, mafe.ConstraintType.Dirichlet, 0.0 );
const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Disp2, mafe.ConstraintType.Dirichlet, 0.0 );
nn = nodes( [gl2.getGridNodes.id] ); % get nodes on GeoLine 2 (right)
const(end+1) = mafe.ConstraintEdge( nn, [mafe.DofType.Disp1, mafe.DofType.Disp2], mafe.ConstraintType.Neumann, [0.0, -1e-1] );
% bending effects
nn = nodes( [gl4.getGridNodes.id] ); % get nodes on GeoLine 4 (left)
const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Disp3, mafe.ConstraintType.Dirichlet, 0.0 );
const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Rota1, mafe.ConstraintType.Dirichlet, 0.0 );
% --------------------------------------------DEFINE FE PROBLEM AND ANALYSIS---
%% 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(); disp( sprintf('---> number dof = %5d', fep.ndofs) );
% -----------------------------------------------------------INSPECT RESULTS---
% %% 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', 1.);
ana.putModeToDof(1)
fep.plotSystem('deformed');
ana.putModeToDof(2)
fep.plotSystem('deformed');
ana.putModeToDof(3)
fep.plotSystem('deformed');
%
ana.L(1:3)
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