The following tutorial discusses multiphysics modeling of resistive
Joule heating simulated with the FEATool Multiphysics
MATLAB FEM toolbox. In the model the resulting current from an applied
electric potential will heat a thin spiral shaped Tungsten wire, such
as can be found in incandescent light bulbs. After the start up phase
the filament reaches an equilibrium temperature where the internal
heat generation is balanced by radiative heat loss through the
boundaries. Although the model could quite easily be implemented in
the FEATool GUI, since it features a one way multiphysics coupling it
will be shown how the dependent variables and equations can be
decoupled and implemented as a MATLAB m-script file to save
computational time.
Model parameters and constants
The first modeling step is to define constants for the equation coefficients and boundary conditions. It is assumed the material of the filament is Tungsten and is held at 20 ℃ at its connecting points in the socket, and is also subject to radiative heat loss with an ambient temperature of 80 ℃. An electric potential of 0.2 V is applied across the wire resulting in temperature rise due to the source term Q. Here, we label this term with the coefficient string qfunction which FEATool will attempt to parse as a general string expression. (Valid equation, boundary, and postprocessing string expressions may consist of combinations of numerical values, dependent variables and derivatives, here V and Vx etc., space coordinates x, y, and z, built-in constants and functions such as pi and sin, and custom functions). Here we will use the fact that as a last resort FEATool will try to parse and evaluate qfunction as a custom user-defined MATLAB m-function the construction of which will be discussed later.
V0 = 0.2; % Potential difference [V].
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].
Q = 'qfunction'; % Resistive heating source term function.
C = 5.670367e-8; % Stefan-Bolzmann radiation constant.
T0 = 20 + 273.15; % Initial temperature [K].
Tamb = 80 + 273.15; % Mean ambient temperature [K].
Grid generation
Since the geometry in this case is quite simple we can omit it and
directly proceed to manually create the grid through using the grid
primitive functions (fea.geom fields are only used by the
unstructured grid generation function gridgen). The computational
grid is thus here created by taking a circle
circgrid
, extruding, and
revolving it to create the filament spiral
gridrevolve
.
% Geometry dimensions in [m].
r_wire = 0.0005; % Wire radius.
l_spiral = 0.005; % Length of spiral.
r_spiral = 0.004; % Outer radius of spiral.
fea.sdim = { 'x' 'y' 'z' };
fea.grid = gridrotate( gridrevolve( ...
circgrid( 2, 3, r_wire ), ...
80, l_spiral, 3, r_spiral ), ...
-pi/2, 2 );
The resulting grid and boundaries can be visualized with the
plotgrid(fea)
and plotbdr(fea)
commands,
respectively.
Electrostatics
The model is quasi-static in the sense that the electric potential solution is static while the overall heat transfer process is examined as a function of time. Although we could solve this as a monolithically coupled problem (solving for T and V simultaneously together), the electric potential problem is assumed to be independent of the temperature (the electrical conductivity $\sigma$ is constant) and can thus be solved separately as a steady state problem by itself, that is
\[ \nabla\cdot(\sigma\nabla V) = 0 \]
with insulation boundary conditions, $\mathbf{n}\cdot(\sigma\nabla V)=0$, everywhere except for the two ends where the potential difference V0 is applied. In the MATLAB m-script of language used by FEATool this looks like the following
% Set up and solve electrostatics problem.
feaV = addphys( fea, @conductivemediadc );
feaV.phys.dc.eqn.coef{2,end} = {sigma}; % Electrical conductivity.
feaV.phys.dc.bdr.sel([1 98]) = 1;
feaV.phys.dc.bdr.coef{1,end}{1} = V0; % Potential difference.
feaV = parsephys( feaV );
feaV = parseprob( feaV );
feaV.sol.u = solvestat( feaV );
Note that here an alternative help variable for the FEM struct
belonging to the electric potential has been used, feaV
.
Heat transfer
The heat transfer problem for the temperature field is governed by the heat equation PDE
\[ \rho c_p\frac{\partial T}{\partial t} - \nabla\cdot (k\nabla T) = Q \]
where here the electric potential is also the source of heat generation with the term $Q = \sigma|\nabla V|^2$. Thus this equation features a non-linear one-way multiphysics coupling, V coupled to the temperature T. As for interaction with the surroundings the temperature is fixed at T = T0 at both ends while radiation flux boundary conditions $\mathbf{n}\cdot(k\nabla T)=C(T_{amb}^4-T^4)$ are applied everywhere else. The following code sets up this heat transfer problem
% Set up heat transfer problem.
fea = addphys( fea, @heattransfer );
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.
% Use generalized heat flux boundary condition (4) for all boundaries
% except at the ends where the temperature T0 is prescribed(1).
fea.phys.ht.bdr.sel([2:97]) = deal(4);
[fea.phys.ht.bdr.coef{4,end}{2:97}] = deal({0 0 0 C Tamb});
fea.phys.ht.bdr.sel([1 98]) = deal(1);
fea.phys.ht.bdr.coef{1,end}{ 1} = T0;
fea.phys.ht.bdr.coef{1,end}{98} = T0;
% Specify heat source term - qfunction.
fea.phys.ht.eqn.coef{7,end} = {'qfunction'};
% Add sigma and electric potential solution used by qfunction.
fea.expr = { 's_sigma' '' '' sigma ;
'u_V' '' '' feaV.so.u };
Custom external source term function
We cannot solve the temperature problem yet since we have not defined
the heat source term qfunction
. We do this in a separate
MATLAB m-file with the same name as the equation coefficient itself,
that is qfunction.m. Although the following code for qfunction
might look intimidating, all it does is to essentially take the
solution for V stored in fea.expr and computes the heat source
term as sigma*(Vx^2+Vy^2+Vz^2)
. Although we could easily
solve the fully monolithically coupled problem with the FEATool GUI,
we would then have to solve a much larger coupled system in each time
step. This decoupled approach we save both time and memory at the cost
of having to use a custom function.
function [ Q ] = qfunction()
% Build up a source term expression from the precomputed solution
% vector stored in prob.expr{ u_V }. Expected to be called from evalexpr0.
% Workaround to get input variable string names from evalexpr0.
persistent inputname0
if( isempty(inputname0) )
ds = dbstack;
fid = fopen( which(ds(2).file), 'r' );
sline = fgetl( fid );
sline = sline( find( sline== '(', 1 )+1:find( sline== ')', 1 )-1 );
inputname0 = regexp( sline(~isspace(sline)), ',', 'split' );
fclose( fid );
end
% Extract variable from the calling function, evalexpr0.
xi = evalin( 'caller', inputname0{2} );
ind_s = evalin( 'caller', inputname0{3} );
ind_c = evalin( 'caller', inputname0{4} );
prob = evalin( 'caller', inputname0{6} );
aJac = evalin( 'caller', inputname0{7} );
n_sdim = size( prob.grid.p, 1 );
n_vert = size( prob.grid.c, 1 );
i_dvar = 1; % Assume same dep. var. fem basis function for u_V as for i_dvar=1.
[~,~,~,sfun] = evalsfun( prob.sfun{i_dvar}, 0, n_sdim, n_vert ); % Get shape function root string.
if( isempty(aJac) ) % Calculate Jacobian if required.
store_aJTmp = strcmpi(sfun(end-1:end),'H3') && ( (n_sdim==2 && n_vert==4) || (n_sdim==3 && n_vert==8) );
[~,aJac] = tfjac( 1, prob.grid.p, prob.grid.c(:,ind_c), [], xi, aJac, store_aJTmp );
end
u_V = prob.expr{ find(strcmp(prob.expr,'u_V')), end }; % V solution vector.
% Evaluate derivatives of V (1=function value, 2=x-derivative, 3=y-derivative, 4=z-derivative).
Vx = evaldvar( sfun, 2, n_sdim, n_vert, xi, aJac, prob.eqn.dofm{i_dvar}(:,ind_c), u_V );
Vy = evaldvar( sfun, 3, n_sdim, n_vert, xi, aJac, prob.eqn.dofm{i_dvar}(:,ind_c), u_V );
Vz = evaldvar( sfun, 4, n_sdim, n_vert, xi, aJac, prob.eqn.dofm{i_dvar}(:,ind_c), u_V );
% Value or expression for sigma.
s_sigma = prob.expr{ find(strcmp(prob.expr,'s_sigma')), end };
% Calculate final expression and return values.
Q = evalexpr0( s_sigma, xi, ind_s, ind_c, [], prob, aJac ).*( Vx.^2 + Vy.^2 + Vz.^2 );
Solution and postprocessing
Finally we can solve the time-dependent heat transfer problem and plot
the solutions keeping in mind that the electric potential is stored in
the feaV
FEM struct and temperature separately in
fea
.
fea = parsephys( fea );
fea = parseprob( fea );
fea.sol.u = solvetime( fea, 'init', {T0}, 'tstep', 1, 'tmax', 20 );
figure
subplot(1,2,1)
postplot( feaV, 'surfexpr', 'V' )
subplot(1,2,2)
postplot( fea, 'surfexpr', 'T' )
After the solver has finished we can see that the maximum reached
temperature is about 700 K. The temporal evolution of the
temperature field can also be seen in the video linked
below. Technically, we could continue our study by also examining the
heat induced stresses and strains, coupling the temperature to the
linear elasticity physics mode. However since for this particular
application the failure point is likely to be due to temperature,
which is far lower than the melting temperature of Tungsten, it is
safe to assume that we do not need to examine the stresses.
The described resistive heating MATLAB model m-file is available from the FEATool examples directory as ex_resistive_heating1.