FEATool Multiphysics  v1.16.6 Finite Element Analysis Toolbox
ex_planestrain2.m File Reference

## Description

EX_PLANESTRAIN2 Plane strain analysis of a pressure vessel.

[ FEA, OUT ] = EX_PLANESTRAIN2( VARARGIN ) Benchmark example for plane strain approximation of a pressure vessel (annular cross section with symmetry).

Reference

[1] Barber JR. Solid Mechanics and Its Applications, Intermediate Mechanics of Materials, pp. 461, vol. 175, 2nd ed. Springer, 2011.

Accepts the following property/value pairs.

Input       Value/{Default}        Description
-----------------------------------------------------------------------------------
hmax        scalar {0.02}          Grid 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 = { ...
'hmax',     0.02;
'sfun',     'sflag2';
'use_temp', true;
'iplot',    1;
'igrid',    1;
'tol',      0.1;
'fid',      1 };
[got,opt] = parseopt( cOptDef, varargin{:} );

E   = 2.1e11;
nu  = 0.3;

ri  = 0.1;
ro  = 0.8;
sri = -5e7;
sro = -2e7;

if( opt.use_temp )
Ti = 500;
To = 20;
K  = (Ti-To)/log(ro/ri);
al = 1.2e-5;
T  = @(r) K*log(ro./r) + To;
Texpr = [num2str(K),'*log(',num2str(ro),'/sqrt(x^2+y^2))+',num2str(To)];
int_rT = @(r) K*(1/2*r.^2.*log(ro./r) + r.^2/4) + r.^2/2*To;
else
Ti = 0;
T0 = 0;
K  = 0;
al = 0;
T  = @(r) 0;
Texpr = 0;
int_rT = @(r) 0;
end

% Reference solution.
C = @(r) E*al/(1-nu)./r.^2.*int_rT(r);
B = ( sro - sri + C(ro) - C(ri) )./(1./ro.^2 - 1./ri.^2);
A = sri + C(ri) - B./ri.^2;

ur  = @(r) al*(1+nu)./((1-nu)*r).*int_rT(r) + A*(1-2*nu)*(1+nu)*r/E - (1+nu)*B./(E*r);
sr  = @(r) -E*al/(1-nu)./r.^2.*int_rT(r) + A + B./r.^2;
sth = @(r)  E*al/(1-nu)./r.^2.*int_rT(r) - E*al*T(r)/(1-nu) + A - B./r.^2;
sz  = @(r) nu*(sr(r)+sth(r)) - E*al*(T(r)-To);   % Adding To here.

% Geometry and grid.
fea.sdim = { 'x' 'y' };   % Coordinate names.
gobj = gobj_ellipse( [0 0], [ro], [ro], 'E1' );
fea.geom.objects{1} = gobj;
gobj = gobj_ellipse( [0 0], [ri], [ri], 'E2' );
fea.geom.objects{2} = gobj;
fea.geom = geom_apply_formula( fea.geom, 'E1-E2' );
gobj = gobj_rectangle( [0], [1.2*ro], [0], [1.2*ro], 'R1' );
fea.geom.objects{2} = gobj;
fea.geom = geom_apply_formula( fea.geom, 'CS1&R1' );
fea.grid = gridgen( fea, 'hmax', opt.hmax, 'fid', opt.fid );

% Problem definition.
fea = addphys( fea, @planestrain );
fea.phys.psn.eqn.coef{1,end} = { nu };
fea.phys.psn.eqn.coef{2,end} = { E  };
fea.phys.psn.eqn.coef{6,end} = { al };
fea.phys.psn.eqn.coef{7,end} = { Texpr };
fea.phys.psn.sfun            = { opt.sfun opt.sfun };

% Boundary conditions.
dtol = 1e-3;
i_b = findbdr( fea, ['y<=',num2str(dtol)] );
i_l = findbdr( fea, ['x<=',num2str(dtol)] );
i_i = findbdr( fea, ['sqrt(x.^2+y.^2)<=',num2str(ri+dtol)] );
i_o = findbdr( fea, ['sqrt(x.^2+y.^2)>=',num2str(ro-dtol)] );

fea.phys.psn.bdr.sel = [ 1, 1, 1, 1 ];

c_setdn = cell(2,4);
[c_setdn{:}] = deal(0);
c_setdn{2,i_b} = 1;
c_setdn{1,i_l} = 1;

c_bdrcoef = cell(2,4);
[c_bdrcoef{:}] = deal(0);
c_bdrcoef{1,i_i} = [num2str(sri),'*nx'];
c_bdrcoef{1,i_o} = [num2str(sro),'*nx'];
c_bdrcoef{2,i_i} = [num2str(sri),'*ny'];
c_bdrcoef{2,i_o} = [num2str(sro),'*ny'];

fea.phys.psn.bdr.coef = { '', '', '', {}, c_setdn, [], c_bdrcoef };

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

% Error checking.
s_sx = fea.phys.psn.eqn.vars{5,end};
s_sy = fea.phys.psn.eqn.vars{6,end};
s_sz = fea.phys.psn.eqn.vars{7,end};

r = linspace(ri,ro,100);
p = [ r; zeros(size(r)) ];

u    = evalexpr( 'u', p, fea ).';
src  = evalexpr( s_sx, p, fea ).';
sthc = evalexpr( s_sy, p, fea ).';
szc  = evalexpr( s_sz, p, fea ).';
uref   = ur(r);
srref  = sr(r);
sthref = sth(r);
szref  = sz(r);
out.err = [ norm(u-uref)/norm(uref), norm(src-srref)/norm(srref), ...
norm(sthc-sthref)/norm(sthref), norm(szc-szref)/norm(szref) ];
p = p([2,1],:);
v    = evalexpr( 'v', p, fea ).';
src  = evalexpr( s_sy, p, fea ).';
sthc = evalexpr( s_sx, p, fea ).';
szc  = evalexpr( s_sz, p, fea ).';
out.err = [ out.err;
norm(v-uref)/norm(uref), norm(src-srref)/norm(srref), ...
norm(sthc-sthref)/norm(sthref), norm(szc-szref)/norm(szref) ];
out.pass = all( out.err(:) <= opt.tol(:) );

% Postprocessing.
if( opt.iplot>0 )
figure
postplot( fea, 'surfexpr', s_sx )
title( 'Stress in the x-direction' )

figure
subplot(2,2,1)
plot( r, uref, 'b-' )
hold on
plot( r, u, 'r.' )
plot( r, v, 'g.' )
xlabel('r')
ylabel('Displacement, u_r')
grid on

subplot(2,2,2)
plot( r, srref, 'b-' )
hold on
plot( r, src, 'r.' )
xlabel('r')
ylabel('Stress, \sigma_r')
grid on

subplot(2,2,3)
plot( r, sthref, 'b-' )
hold on
plot( r, sthc, 'r.' )
xlabel('r')
ylabel('Stress, \sigma_\theta')
grid on

subplot(2,2,4)
plot( r, szref, 'b-' )
hold on
plot( r, szc, 'r.' )
xlabel('r')
ylabel('Stress, \sigma_z')
grid on
end

if( nargout==0 )
clear fea out
end