FEATool Multiphysics  v1.16.4 Finite Element Analysis Toolbox
ex_planestress2.m File Reference

## Description

EX_PLANESTRESS2 NAFEMS benchmark challenge 1 plane stress example.

[ FEA, OUT ] = EX_PLANESTRESS2( VARARGIN ) NAFEMS benchmark example to calculate von Mieses stress thin plate under plane stress assumption with loads all around.

Accepts the following property/value pairs.

Input       Value/{Default}        Description
-----------------------------------------------------------------------------------
E           scalar {210e9}         Modulus of elasticity
nu          scalar {0.3}           Poissons ratio
igrid       scalar 1/{0}           Cell type (0=quadrilaterals, 1=triangles)
hmax        scalar {1/20}          Max grid cell size
sfun        string {sflag1}        Shape function for displacements
iplot       scalar 0/{1}           Plot solution (=1)
.
Output      Value/(Size)           Description
-----------------------------------------------------------------------------------
fea         struct                 Problem definition struct
out         struct                 Output struct


# Code listing

 cOptDef = { ...
'E',        210e9; ...
'nu',       0.3; ...
'igrid',    0; ...
'hmax',     1/20; ...
'sfun',     'sflag1'; ...
'iplot',    1; ...
'fid',      1 };
[got,opt] = parseopt(cOptDef,varargin{:});
fid       = opt.fid;

% Geometry definition.
gobj1 = gobj_rectangle( 0, 1, 0, 1, 'R1' );
fea.geom.objects = { gobj1 };
fea.sdim = { 'x' 'y' };

% Grid generation.
if( opt.igrid==1 )
fea.grid = gridgen(fea,'hmax',opt.hmax,'fid',fid);
else
nx = 1/opt.hmax;
fea.grid = rectgrid( nx, nx );
end

% Boundary conditions.
dtol = 0.1;
if( opt.igrid==1 )
lbdr = findbdr( fea, ['x<',num2str(dtol)] );     % Left boundary number.
rbdr = findbdr( fea, ['x>',num2str(1-dtol)] );   % Right boundary number.
tbdr = findbdr( fea, ['y>',num2str(1-dtol)] );   % Top boundary number.
bbdr = findbdr( fea, ['y<',num2str(dtol)] );     % Bottom boundary number.
else
lbdr = 4;
rbdr = 2;
tbdr = 3;
bbdr = 1;
end

% Add plane stress physics mode.
fea.phys.pss.eqn.coef{1,end} = { opt.nu };
fea.phys.pss.eqn.coef{2,end} = { opt.E  };
fea.phys.pss.sfun            = { opt.sfun opt.sfun };

% Set boundary condition types.
bctype = mat2cell( zeros(2,4), [1 1], [1 1 1 1] );
fea.phys.pss.bdr.coef{1,5}   = bctype;

bccoef = mat2cell( zeros(2,4), [1 1], [1 1 1 1] );
bccoef{1,tbdr} = '-x';
bccoef{2,tbdr} = 'x';
bccoef{1,bbdr} = '-(1-x)';
bccoef{2,bbdr} = '1-x';
bccoef{1,lbdr} = '1-y';
bccoef{2,lbdr} = '-(1-y)';
bccoef{1,rbdr} = 'y';
bccoef{2,rbdr} = '-y';
fea.phys.pss.bdr.coef{1,end} = bccoef;

% Parse and solve problem.
fea       = parsephys(fea);             % Check and parse physics modes.
fea       = parseprob(fea);             % Check and parse problem struct.
fea.sol.u = solvestat(fea,'fid',fid);   % Call to stationary solver.

% Postprocessing.
s_vm = fea.phys.pss.eqn.vars{1,2};
if ( opt.iplot>0 )
figure
postplot( fea, 'surfexpr', s_vm, 'isoexpr', s_vm )
title('von Mieses stress')
end

% Error checking.
dtol = sqrt(eps);
sm  = evalexpr( s_vm, [0.5;0.5], fea );
s01 = evalexpr( s_vm, [dtol;1],    fea );
s10 = evalexpr( s_vm, [1;dtol],    fea );
s00 = evalexpr( s_vm, [dtol;dtol], fea );
s11 = evalexpr( s_vm, [1;1],       fea );

sol = [ sm s01 s10 s00 s11 ];
ref = [  0   0   0   2   2 ];
out.err    = norm(sol-ref)/norm(ref);
out.pass   = out.err < 0.1;

if ( nargout==0 )
clear fea out
end