FEATool Multiphysics  v1.16.4 Finite Element Analysis Toolbox
ex_resistive_heating1.m File Reference

Description

EX_RESISTIVE_HEATING1 Resistive heating of a tungsten filament.

[ FEA, OUT ] = EX_RESISTIVE_HEATING1( VARARGIN ) Example to calculate temperature generated by resistive heating of a tungsten spiral filament, such as for example in an incandescent bulb. Accepts the following property/value pairs.

Input       Value/{Default}        Description
-----------------------------------------------------------------------------------
ibc         scalar 1{2}            Boundary type (1=insulation, 2=radiation)
isolve      scalar 1{2}            Solve type (1=stat.coupled, 2=split quasistatic)
sfun        string {sflag1}        Shape function
iplot       scalar 0/{1}           Plot solution (=1)
.
Output      Value/(Size)           Description
-----------------------------------------------------------------------------------
fea         struct                 Problem definition struct
out         struct                 Output struct


Code listing

 cOptDef = { 'ibc',      2;
'isolve',   2;
'sfun',     'sflag1';
'iplot',    1;
'tol',      1e-2;
'fid',      1 };
[got,opt] = parseopt(cOptDef,varargin{:});
fid       = opt.fid;

% Set up geometry and grid.
r_wire   = 0.0005;
l_spiral = 0.005;
r_spiral = 0.004;

fea.sdim = { 'x' 'y' 'z' };
fea.grid = gridrotate( gridrevolve( circgrid( 2, 3, r_wire ), 80, l_spiral, 3, r_spiral ), -pi/2, 2 );
n_bdr = max(fea.grid.b(3,:));
ib2 = findbdr( fea, 'y<1e-3 & y>-1e-3' );
ib1 = setdiff(1:n_bdr,ib2);

% Problem constants and definitions.
V0    = 0.2;         % Voltage difference.
sigma = 1/52.8e-9;   % Electrical conductivity [1/Ohm/m].

rho   = 19.25e3;     % Density [kg/m3].
cp    = 133.9776;    % Specific heat capacity [J/kg/K].
k     = 173;         % Thermal conductivity [W/m/K].
if( opt.isolve==1 )
Q   = [num2str(sigma),'*(Vx^2+Vy^2+Vz^2)'];   % Resistive heating source term.
else
Q   = 'qfunction'; % Resistive heating source term function.
end

Const = 5.670367e-8;       % Constant for radiation BC (Stefan-Bolzmann constant).
T0    = 20+273.15;         % Initial temperature.
Tamb  = 80+273.15;         % Mean ambient temperature.

icub  = str2num(opt.sfun(end))+1; % Numerical quadrature order.
tmax  = 20;                 % Maximum time for temperature simulation.
dt    = tmax/20;            % Time step size.

% Set up and solve electrostatics problem.
feaV = fea;
feaV.phys.dc.sfun = { opt.sfun };

feaV.phys.dc.eqn.coef{2,end} = {sigma};    % Electrical conductivity.

feaV.phys.dc.bdr.sel(ib2) = 1;                 % Use Dirichlet / Electric potential BCs for end boundaries 57 and 58.
feaV.phys.dc.bdr.coef{1,end}{ib2(1)} = V0;     % Set electric potential on boundary 1 to V0 (boundary 57 is 0 by default).

if( opt.isolve==2 )
feaV = parsephys( feaV );
feaV = parseprob( feaV );
u_V  = solvestat( feaV, 'fid', opt.fid );
end

% Set up and solve heat transfer problem.
fea = addphys( fea, @heattransfer );
fea.phys.ht.sfun = { opt.sfun };

fea.phys.ht.eqn.coef{1,end} = {rho};       % Density.
fea.phys.ht.eqn.coef{2,end} = {cp};        % Heat capacity.
fea.phys.ht.eqn.coef{3,end} = {k};         % Thermal conductivity.
fea.phys.ht.eqn.coef{7,end} = {Q};         % Heat source term.

if( opt.ibc==1 )
fea.phys.ht.bdr.sel(ib1) = deal(3);     % Use insulation flux boundary condition for all boundaries except the ends.
else
fea.phys.ht.bdr.sel(ib1) = deal(4);     % Use generalized heat flux boundary condition for all boundaries except the ends.
[fea.phys.ht.bdr.coef{4,end}{ib1}] = deal({0 0 0 Const Tamb});
end
fea.phys.ht.bdr.sel(ib2) = deal(1);     % Prescribe fixed ambient temperature at the end points.
fea.phys.ht.bdr.coef{1,end}{ib2(1)} = T0;
fea.phys.ht.bdr.coef{1,end}{ib2(2)} = T0;

if( opt.isolve==2 )
% Add sigma and electric potential solution used by qfunction.
fea.expr = { 's_sigma' '' '' sigma ;
'u_V'     '' '' u_V   };

fea  = parsephys( fea );
fea  = parseprob( fea );
l_write_qfunction()
u_T  = solvetime( fea, 'icub', icub, 'init', {T0}, 'tstep', dt, 'tmax', tmax, 'fid', opt.fid );
delete( fullfile(tempdir(),'qfunction.m') );
clear qfunction;
end

% Merge solutions.
fea.phys.dc = feaV.phys.dc;
if( isfield(fea,'dvar') )
fea = rmfield( fea, 'dvar' );
fea = rmfield( fea, 'eqn'  );
end
fea = parsephys( fea );
fea = parseprob( fea );
if( opt.isolve==1 )
fea.sol.u = solvestat( fea, 'icub', icub, 'fid', opt.fid, 'nlrlx', 1-0.5*(opt.ibc==2) );
u_T = fea.sol.u(1:fea.eqn.ndof(1),:);
else
fea.sol.u = [ u_T; repmat(u_V,1,size(u_T,2)) ];
end

% Postprocessing.
if( opt.iplot>0 )
subplot(1,2,1)
postplot( fea, 'surfexpr', 'V' )
title( 'Electric potential, V')

subplot(1,2,2)
postplot( fea, 'surfexpr', 'T-273.15' )
title( 'Temperature, T [C]')
end

% Error checking.
if( opt.isolve==1 )
if( opt.ibc==1 )
T_max_ref = 840;
else
T_max_ref = 688;
end
else
if( opt.ibc==1 )
T_max_ref = 787;
else
T_max_ref = 680;
end
end
out.err  = abs( max(u_T(:)) - T_max_ref )/T_max_ref;
out.pass = out.err < opt.tol;

if( nargout==0 )
clear fea out
end

%-----------------------------