|
FEATool Multiphysics
v1.17.5
Finite Element Analysis Toolbox
|
|
|
FEATool Multiphysics
v1.17.5
Finite Element Analysis Toolbox
|
EX_PLANESTRESS1 Example for plane stress for hole in plate.
[ FEA, OUT ] = EX_PLANESTRESS1( VARARGIN ) Example to calculate displacements and stresses for a hole in plate configuration under plane stress assumption.
Accepts the following property/value pairs.
Input Value/{Default} Description
-----------------------------------------------------------------------------------
E scalar {210e9} Modulus of elasticity
nu scalar {0.3} Poissons ratio
diam scalar {0.01} Diameter of hole
thick scalar {0.001} Plate thickness
force scalar {1000} Load force
sx_ref scalar {3e7} Reference stress in x-direction
hmax scalar {0.00125} Max grid cell size
sfun string {sflag1} Shape function for displacements
iphys scalar 0/{1} Use physics mode to define problem (=1)
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; ...
'diam', 0.01; ...
'thick', 0.001; ...
'force', 1000; ...
'sx_ref', 3e7; ...
'hmax', 0.00125; ...
'sfun', 'sflag1'; ...
'iphys', 1; ...
'iplot', 1; ...
'tol', 0.05; ...
'fid', 1 };
[got,opt] = parseopt(cOptDef,varargin{:});
fid = opt.fid;
% Model, geometry, and grid parameters.
h = 0.05; % Height of 1/4 rectangular domain.
l = 0.05; % Length of 1/4 rectangular domain.
xc = 0; % x-coordinate of hole center.
yc = 0; % y-coordinate of hole center.
area = 2*h*opt.thick; % Area on which load force is applied.
% Geometry definition.
gobj1 = gobj_rectangle( 0, l, 0, h, 'R1' );
gobj2 = gobj_circle( [xc,yc], opt.diam/2, 'C1' );
fea.geom.objects = { gobj1 gobj2 };
fea = geom_apply_formula( fea, 'R1-C1' );
fea.sdim = { 'x' 'y' }; % Coordinate names.
% Grid generation.
fea.grid = gridgen(fea,'hmax',opt.hmax,'fid',fid);
n_bdr = max(fea.grid.b(3,:)); % Number of boundaries.
% Boundary conditions.
dtol = opt.diam/10;
lbdr = findbdr( fea, ['x<=',num2str(dtol)] ); % Left boundary number.
rbdr = findbdr( fea, ['x>=',num2str(l-dtol)] ); % Right 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;
if ( opt.iphys==1 )
fea = addphys(fea,@planestress); % 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 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;
bccoef = mat2cell( zeros(2,n_bdr), [1 1], ones(1,n_bdr) );
bccoef{1,rbdr} = opt.force/area;
fea.phys.pss.bdr.coef{1,end} = bccoef;
fea = parsephys(fea); % Check and parse physics modes.
else
fea.dvar = { 'u' 'v' }; % Dependent variable names.
fea.sfun = { opt.sfun opt.sfun }; % Shape functions.
% Define equation system.
fea.eqn.a.form = { [2 3;2 3] [3 2;2 3]; ...
[3 2;2 3] [2 3;2 3] }; ...
fea.eqn.a.coef = { {E11 E33} {E12 E33}; ...
{E33 E12} {E33 E11} };
fea.eqn.f.form = { 1 1 };
fea.eqn.f.coef = { 0 0 };
% Define boundary conditions.
fea.bdr.d = cell(2,n_bdr);
fea.bdr.n = cell(2,n_bdr);
[fea.bdr.n{:}] = deal(0); % Assign zero to all Neumann boundaries.
fea.bdr.n{1,rbdr} = opt.force/area; % Set horizontal load force on right boundary.
fea.bdr.d{1,lbdr} = 0; % Set zero horizontal displacement on left boundary.
fea.bdr.d{2,lobdr} = 0; % Set zero vertical displacement on lower boundary.
end
% Parse and solve problem.
fea = parseprob(fea); % Check and parse problem struct.
fea.sol.u = solvestat(fea,'fid',fid); % Call to stationary solver.
% Postprocessing.
s_sx = [num2str(E11),'*ux+',num2str(E12),'*vy'];
if ( opt.iplot>0 )
figure
postplot(fea,'surfexpr',s_sx,'isoexpr',s_sx,'isomap','k')
title('Stress, x-component')
end
% Error checking.
[~,sx_max] = minmaxsubd( s_sx, fea );
out.sx_max = sx_max;
out.err = opt.sx_ref-sx_max;
out.pass = (sx_max>(opt.sx_ref*(1-opt.tol)))&&(sx_max<(opt.sx_ref*(1+opt.tol)));
if ( nargout==0 )
clear fea out
end