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
force       scalar {1e5}           Load force
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 = addphys(fea,@linearelasticity);
 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


%