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

## Description

EX_PLANESTRESS5 Plane stress example for an elliptic membrane.

[ FEA, OUT ] = EX_PLANESTRESS5( VARARGIN ) Example to calculate displacements and stresses for an elliptic membrane with a hole in it. NAFEMS benchmark example LE1.

Accepts the following property/value pairs.

Input       Value/{Default}        Description
-----------------------------------------------------------------------------------
E           scalar {210e9}         Modulus of elasticity
nu          scalar {0.3}           Poissons ratio
hmax        scalar {0.1}           Max grid cell size
sfun        string {sflag2}        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; ...
'hmax',     0.1; ...
'sfun',     'sflag2'; ...
'iplot',    1; ...
'igeom',    1; ...
'tol',      0.05; ...
'fid',      1 };
[got,opt] = parseopt( cOptDef, varargin{:} );
fid       = opt.fid;

% Geometry definition.
gobj1 = gobj_ellipse( [0 0], 3.25, 2.75, 'E1' );
gobj2 = gobj_ellipse( [0 0], 2, 1, 'E2' );
gobj3 = gobj_rectangle( -3.25, 3.25, -2.75, 0, 'R1' );
gobj4 = gobj_rectangle( -3.25, 0, 0, 2.75, 'R2' );
fea.geom.objects = { gobj1 gobj2 gobj3 gobj4 };
if( opt.igeom==1 )
fea = geom_apply_formula( fea, 'E1-E2-R1-R2' );
else
fea = geom_apply_formula( fea, 'E1-E2' );
fea = geom_apply_formula( fea, 'CS1-R1' );
fea = geom_apply_formula( fea, 'CS2-R2' );
end
fea.sdim = { 'x' 'y' };

% Grid generation.
if( opt.igeom==1 )
fea.grid = gridgen( fea, 'hmax', opt.hmax, 'fid', fid, 'gridgen', 'gridgen2d' );
else
fea.grid = gridgen( fea, 'hmax', opt.hmax, 'fid', fid );
end
n_bdr = max(fea.grid.b(3,:));   % Number of boundaries.

% Boundary conditions.
dtol  = sqrt(eps);
lbdr  = findbdr( fea, ['x<=',num2str(dtol)] );   % Left boundary number.
lobdr = findbdr( fea, ['y<=',num2str(dtol)] );   % Lower boundary number.

% Problem definition.
E11 = opt.E/(1-opt.nu^2);
E12 = opt.nu*E11;
E22 = E11;
E33 = opt.E/(1+opt.nu)/2;

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 shape functions.

bctype = mat2cell( zeros(2,n_bdr), [1 1], ones(1,n_bdr) );
bctype{1,lbdr}  = 1;
bctype{2,lobdr} = 1;
fea.phys.pss.bdr.coef{1,5} = bctype;

dtol = 1e-3;
i_o = findbdr( fea, ['sqrt(x.^2+y.^2)>=',num2str(2.75-dtol)] );
bccoef = mat2cell( zeros(2,n_bdr), [1 1], ones(1,n_bdr) );
bccoef{1,i_o} = 'nx*10e6';
bccoef{2,i_o} = 'ny*10e6';
fea.phys.pss.bdr.coef{1,end} = bccoef;

% Parse and solve problem.
fea       = parsephys(fea);
fea       = parseprob(fea);
fea.sol.u = solvestat( fea, 'fid', fid );

% Postprocessing.
s_sy = [num2str(E12),'*ux+',num2str(E11),'*vy'];
if ( opt.iplot>0 )
figure
postplot( fea, 'surfexpr', s_sy, 'isoexpr', s_sy )
title('Stress, x-component')
end

% Error checking.
sy_D     = evalexpr( s_sy, [2;0]+sqrt(eps)*1e1, fea );
out.sy_D = sy_D;
out.err  = abs(sy_D - 92.7e6)/92.7e6;
out.pass = out.err <= opt.tol;

if ( nargout==0 )
clear fea out
end