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

## Description

EX_LINEARELASTICITY4 Stress calculation of an I-beam attached to two brackets.

[ FEA, OUT ] = EX_LINEARELASTICITY4( VARARGIN ) Example to calculate displacements and stresses for an I-beam suppored by two brackets with circular holes.

Accepts the following property/value pairs.

Input       Value/{Default}        Description
-----------------------------------------------------------------------------------
E           scalar {200e9}         Modulus of elasticity
nu          scalar {0.3}           Poissons ratio
l           scalar {0.4}           Length of I-beam
ilev        scalar {2}             Grid regfinement level
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',        200e9; ...
'nu',       0.3; ...
'force',    1e5; ...
'l',        0.4; ...
'ilev',     2; ...
'sfun',     'sflag1'; ...
'iplot',    1; ...
'tol',      0.42; ...
'fid',      1 };
[got,opt] = parseopt(cOptDef,varargin{:});
fid       = opt.fid;

% Geometry definition.
fea.sdim = { 'x' 'y' 'z' };   % Coordinate names.

% Grid generation.
fea.grid = get_grid( opt.ilev );

% Problem definition.
fea.phys.el.eqn.coef{1,end} = { opt.nu };
fea.phys.el.eqn.coef{2,end} = { opt.E  };
fea.phys.el.sfun            = { opt.sfun opt.sfun opt.sfun };

% Boundary conditions.
dtol     = sqrt(eps);
fixbdr   = findbdr( fea, ['(sqrt(x.^2+z.^2)<=0.03+sqrt(eps))&(z>=-sqrt(eps))'] );
forcebdr = findbdr( fea, ['abs(y)>=0.2-sqrt(eps)'] );

% Fix boundaries (set zero Dirichlet BCs).
n_bdr  = max(fea.grid.b(3,:));        % Number of boundaries.
bctype = num2cell( zeros(3,n_bdr) );  % First set homogenous Neumann BCs everywhere.
[bctype{:,fixbdr}] = deal( 1 );       % Set Dirchlet BCs for right boundary.
fea.phys.el.bdr.coef{1,5} = bctype;

% Apply negative z-load to left boundary.
bccoef = num2cell( zeros(3,n_bdr) );
[bccoef{3,forcebdr}] = deal(-opt.force);
fea.phys.el.bdr.coef{1,end} = bccoef;

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

% Postprocessing.
if ( opt.iplot>0 )
DSCALE = 5000;

subplot(1,2,1)
postplot( fea, 'surfexpr', 'sqrt(u^2+v^2+w^2)', 'linestyle', 'none' )
title( 'Total displacement' )
view([30 20])

subplot(1,2,2)
dp = zeros(size(fea.grid.p));
for i=1:3
dp(i,:) = DSCALE*evalexpr( fea.dvar{i}, fea.grid.p, fea );
end
fea_disp.grid   = fea.grid;
fea_disp.grid.p = fea_disp.grid.p + dp;
plotgrid( fea_disp )
title(['Displacement plot (at ',num2str(DSCALE),' times scale)'])
view([30 20])

end

% Error check.
disp_max_ref = 6.204e-6;
xdisp = fea.sol.u(fea.eqn.dofm{1}(:));
ydisp = fea.sol.u(fea.eqn.dofm{2}(:)+fea.eqn.ndof(1));
zdisp = fea.sol.u(fea.eqn.dofm{3}(:)+sum(fea.eqn.ndof(1:2)));
disp  = sqrt(xdisp.^2+ydisp.^2+zdisp.^2);
disp_max = max(disp);

svm_max_ref = 4.410e6;
svm = evalexpr( fea.phys.el.eqn.vars{1,2}, fea.grid.p, fea );
svm_max = max(svm);

out.disp_max = disp_max;
out.svm_max  = svm_max;
out.err(1)   = abs(disp_max - disp_max_ref)/abs(disp_max_ref);
out.err(2)   = abs(svm_max - svm_max_ref)/abs(svm_max_ref);
out.pass     = all(out.err<opt.tol);

if ( nargout==0 )
clear fea out
end

%