FEATool Multiphysics  v1.16.1 Finite Element Analysis Toolbox
ex_heattransfer2.m File Reference

## Description

EX_HEATTRANSFER2 1D Stationary heat transfer with radiation.

[ FEA, OUT ] = EX_HEATTRANSFER2( VARARGIN ) NAFEMS T2 benchmark example for heat transfer with radiation [1]. The left end of a 0.1 m rod is held at a temperature of 1000 K while the right end is radiating with an emissivity, em = 0.98, and Stefan-Bolzmann constant, sigma = 5.67e-8 Wm^2/K^4.

       +---------- L=0.1m ----------+ T(0.1)?
T=1000K                   q_n = em*sigma*(T_amb^4-T^4)


The rod is made of iron with density 7850 kg/m^3, heat capacity 460 J/kgK, and thermal conductivity 55.563 W/mK. The steady state temperature at the left end is sought when the surrounding ambient temperature is 300 K.

Reference
  [1] The Standard NAFEMS Benchmarks,
The National Agency for Finite Element Standards, UK, 1990.


Accepts the following property/value pairs.

Input       Value/{Default}        Description
-----------------------------------------------------------------------------------
hmax        scalar {0.02}          Grid cell size
sfun        string {sflag1}        Finite element shape function
solver      string fenics/{}       Use FEniCS or default solver
istat       scalar {1}/0           Use stationary (=1), or time dependent solver
iplot       scalar {1}/0           Plot solution (=1)
.
Output      Value/(Size)           Description
-----------------------------------------------------------------------------------
fea         struct                 Problem definition struct
out         struct                 Output struct


# Code listing

 cOptDef = { 'hmax',     0.02;
'sfun',     'sflag1';
'solver',   '';
'istat',    1;
'iplot',    1;
'fid',      1 };
[got,opt] = parseopt(cOptDef,varargin{:});

% Grid generation.
L        = 0.1;
nx       = round(L/opt.hmax);
fea.grid = linegrid( nx, 0, L );

% Problem definition.
fea.sdim  = { 'x' };                      % Space coordinate name.
fea = addphys( fea, @heattransfer );      % Add heat transfer physics mode.
fea.phys.ht.sfun = { opt.sfun };          % Set shape function.

% Equation coefficients.
fea.phys.ht.eqn.coef{1,end} = 7850;       % Density
fea.phys.ht.eqn.coef{2,end} =  460;       % Heat capacity.
fea.phys.ht.eqn.coef{3,end} =   55.563;   % Thermal conductivity.
fea.phys.ht.eqn.coef{6,end} = { 1000 };   % Initial temperature.

% Boundary conditions.
fea.phys.ht.bdr.sel = [ 1 4 ];
fea.phys.ht.bdr.coef{1,end} = { 1000 [] };
fea.phys.ht.bdr.coef{4,end}{2}{4} = '0.98*5.67e-8';
fea.phys.ht.bdr.coef{4,end}{2}{5} = 300;

% Parse physics modes and problem struct.
fea = parsephys(fea);
fea = parseprob(fea);

% Compute solution.
if( strcmp(opt.solver,'fenics') )
fea = fenics( fea, 'fid', opt.fid, ...
'tstep', 10, 'tmax', 1000, 'ischeme', 2*(~opt.istat) );
else
if( opt.istat )
fea.sol.u = solvestat( fea, 'fid', opt.fid, 'init', {'T0_ht'} );
else
[fea.sol.u, tlist] = solvetime( fea, 'fid', opt.fid, 'init', {'T0_ht'}, ...
'tmax', 1000, 'tstep', 10 );
end
end

% Postprocessing.
if( opt.iplot>0 )
postplot( fea, 'surfexpr', 'T', 'axequal', 0 )
title('Temperature')
xlabel('x')
ylabel('T')
end

% Error checking.
T_sol = evalexpr( 'T', 0.1, fea );
T_ref = 926.97;
out.err  = abs(T_sol-T_ref)/T_ref;
out.pass = out.err<6e-4;

if( nargout==0 )
clear fea out
end