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.