FEATool Multiphysics
v1.17.1
Finite Element Analysis Toolbox
|
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
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 = addphys(fea,@planestress); 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; % Set loads on boundary. 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