Commit 40d1781a authored by Andreas Zilian's avatar Andreas Zilian
Browse files

Add examples for flat shell analysis.

parent 3f852dba
%% Example Flat Shell 1a
% A rectangular flat shell structure composed of 1 element
addpath([ '..' filesep '..' filesep '..' filesep ]);
%% Definition of sections
s1 = mafe.Section2D();
s1.E = 1.0; % [N/m^2]
s1.nue = 0.0; % [-]
s1.t = 1.0; % [m]
%% Definition of nodes in the system
% Node 1
n1 = mafe.Node2D( [0.0, 0.0] );
% Node 2
n2 = mafe.Node2D( [1.0, 0.0] );
% Node 3
n3 = mafe.Node2D( [1.0, 1.0] );
% Node 4
n4 = mafe.Node2D( [0.0, 1.0] );
% Vector of nodes
nodes = [ n1, n2, n3, n4 ];
%% Definition of elements in the system
% Element 1
e1 = mafe.FlatShell2D( [n1, n2, n3, n4], s1 );
% vector of elements
elems = [ e1 ];
%% Definition of constraints in the system
% Constraint 1 (fixed transversal displacement at left end, lower node)
c1 = mafe.ConstraintNode( n1, mafe.DofType.Disp3, mafe.ConstraintType.Dirichlet, 0.0 );
% Constraint 2 (fixed transversal displacement at left end, upper node)
c2 = mafe.ConstraintNode( n4, mafe.DofType.Disp3, mafe.ConstraintType.Dirichlet, 0.0 );
% Constraint 3 (fixed transversal displacement at left end, lower node)
c3 = mafe.ConstraintNode( n1, mafe.DofType.Rota1, mafe.ConstraintType.Dirichlet, 0.0 );
% Constraint 4 (fixed transversal displacement at left end, upper node)
c4 = mafe.ConstraintNode( n4, mafe.DofType.Rota1, mafe.ConstraintType.Dirichlet, 0.0 );
% Constraint 5 (moment at right end, lower node)
c5 = mafe.ConstraintNode( n2, mafe.DofType.Rota1, mafe.ConstraintType.Neumann, 0.1 ); % 0.1 [Nm]
% Constraint 6 (moment at right end, upper node)
c6 = mafe.ConstraintNode( n3, mafe.DofType.Rota1, mafe.ConstraintType.Neumann, 0.1 ); % 0.1 [Nm]
% Constraint 7 (fixed horizontal displacement at left end, lower node)
c7 = mafe.ConstraintNode( n1, mafe.DofType.Disp1, mafe.ConstraintType.Dirichlet, 0.0 );
% Constraint 8 (fixed vertical displacement at left end, lower node)
c8 = mafe.ConstraintNode( n1, mafe.DofType.Disp2, mafe.ConstraintType.Dirichlet, 0.0 );
% Constraint 9 (fixed horizontal displacement at left end, upper node)
c9 = mafe.ConstraintNode( n4, mafe.DofType.Disp1, mafe.ConstraintType.Dirichlet, 0.0 );
% Constraint 10 (horizontal force at right end, lower node)
c10 = mafe.ConstraintNode( n2, mafe.DofType.Disp1, mafe.ConstraintType.Neumann, -0.1 );
% Constraint 11 (horizontal force at right end, upper node)
c11 = mafe.ConstraintNode( n3, mafe.DofType.Disp1, mafe.ConstraintType.Neumann, -0.1 );
% vector of constraints
const = [ c1, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11 ];
%% 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;
fep.plotSystem('reference');
fep.plotSystem('deformed', 1.e0);
%fep.plotSystem('moment11'); colorbar;
%fep.plotSystem('shear13'); colorbar;
%fep.plotSystem('tensor');
fep.printSystem();
%% Example Plate Bending 2a
% 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)
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.85] );
% 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
%
testcase = 'bend';
%
% 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, -1e2] );
% bending effects
switch(testcase)
case 'bend'
% impose constant m11, support in Disp3
%
% left edge: supported
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 );
% right edge: distributed moment = 1.0 [Nm/m]
nn = nodes( [gl2.getGridNodes.id] ); % get nodes on GeoLine 2 (right)
const(end+1) = mafe.ConstraintEdge( nn, [mafe.DofType.Rota1, mafe.DofType.Rota2], mafe.ConstraintType.Neumann, [0.0, 1e-4] );
%
case 'shear'
% impose constant shear
%
% left edge: supported
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 );
% all nodes: fix rotations
nn = nodes( [ga1.getGridNodes.id] ); % get nodes on GeoArea 1
const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Rota1, mafe.ConstraintType.Dirichlet, 0.0 );
const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Rota2, mafe.ConstraintType.Dirichlet, 0.0 );
% right edge: vertical distributed force = 1.0 [N/m]
nn = nodes( [gl2.getGridNodes.id] ); % get nodes on GeoLine 2 (right)
const(end+1) = mafe.ConstraintEdge( nn, [mafe.DofType.Disp3, mafe.DofType.Disp3], mafe.ConstraintType.Neumann, [0.0, 1e-4] );
%
case 'twist'
% impose constant twist
%
% left edge: supported
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 );
% lower edge: supported
nn = nodes( [gl1.getGridNodes.id] ); % get nodes on GeoLine 1 (bottom)
const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Disp3, mafe.ConstraintType.Dirichlet, 0.0 );
% top right corner: single force = 1.0 [N]
nn = nodes( [gp3.getGridNodes.id] ); % get node at GeoPoint 3 (top right corner)
const(end+1) = mafe.ConstraintNode( nn, mafe.DofType.Disp3, mafe.ConstraintType.Neumann, 1e-4 );
%
otherwise
disp('Unknown test case! Constraints not set');
const = [];
end
% --------------------------------------------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(); 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.);
%fep.plotSystem('shear13'); colorbar();
%fep.plotSystem('tensor', 1.e-2);
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