Using time dependent (heat) source in PDE Toolbox

Question

My goal is to apply a time-dependent heat source when solving the heat transfer problem.

The partial differential equation for transient conduction heat transfer is:

and more information can be found here: Solving a Heat Transfer Problem With Temperature-Dependent Properties

All parameters are constants in my case, except the source term, f, needs to be changed along with time.

I'm following the example code here: Nonlinear Heat Transfer In a Thin Plate which gives a way to solve the transient problem, and I'm able to plot the heat data at each time point.

The problem when applying it to my case is that, in the example the source is a constant value, throughout the entire region and entire time, and related to radiation and convection (in my case they should be all zero), but I need to feed a time-dependent source (Joule heating by time-varying electric current). The source could have one of the following formats:

1. Analytical: a positive value such as 1 W/m^2, within a time window such as 0< t< 1 ns, and 0 otherwise.
2. Numerical: data is provided by a 1xN vector where N is the number of time points.

And the source is confined in a certain region, eg. 0< x <1 mm and 0< y<1 mm.

I have seen a similar question but unanswered: How to use a variable coefficient in PDE Toolbox to solve a parabolic equation (Matlab)

Is there a way to implement this with the PDE Toolbox? Writing code from scratch would be so complicated....

Answer 1

You can quite easily define and solve problems with time dependent and nonlinear PDE coefficients with the FEATool FEM Matlab Toolbox as shown here in the m-script code snippet below. Note that the heat source (sink) term f is scaled as f*(t>2500) which means that it will only be active after t=2500 (as the switch expression evaluates to either 0 if false or 1 if true).

% Coefficents and problem definition from https://www.mathworks.com/help/pde/examples/nonlinear-heat-transfer-in-a-thin-plate.html
k=400;rho=8960;specificHeat=386;thick=.01;stefanBoltz=5.670373e-8;hCoeff=1;ta=300;emiss=.5;

% Set up 2D fea struct with geometry and grid.
fea.sdim = {'x' 'y'};
fea.geom = { gobj_rectangle( 0, 1, 0, 1 ) };
fea.grid = rectgrid( 10 );

% Add heat transfer physics mode.
fea = addphys( fea, @heattransfer );

fea.phys.ht.eqn.coef{1,end}{1} = rho*thick;      % Density eqn coefficient.
fea.phys.ht.eqn.coef{2,end}{1} = specificHeat;   % C_p eqn coefficient.
fea.phys.ht.eqn.coef{3,end}{1} = k*thick;        % Thermal condictivity.
f = sprintf( '%g*( %g - T ) + %g*( %g^4 - T^4 )', ...
             2*hCoeff, ta, 2*emiss*stefanBoltz, ta );
fea.phys.ht.eqn.coef{6,end}{1} = ['(',f,')*(t>2500)'];   % Heat source term.

fea.phys.ht.bdr.sel(1) = 1;  % Set prescribed temperature for boundary 1.
fea.phys.ht.bdr.coef{1,end}{1} = 1000;  
fea.phys.ht.bdr.sel(2:4) = 3;   % Isolation BCs for boundaries 2-4.

% Check, parse, and solve fea problem.
fea = parsephys( fea );
fea = parseprob( fea );
[fea.sol.u,t] = solvetime( fea, 'tstep', 50, 'tmax', 5000, 'init', {ta} );

% Postprocessing and visualization.
for i=1:size(fea.sol.u,2)
  T_top(i) = evalexpr( 'T', [.5;1-sqrt(eps)], fea, i );
end

subplot(1,2,1)
postplot( fea, 'surfexpr', 'T', 'title', 'T @ t=5000' )
subplot(1,2,2)
plot( t, T_top, 'r-' )
xlabel( 't' )
ylabel( 'T(0.5,1) @ t=5000' )
grid on

In the solution here you can see that the temperature of the top edge rises linearly due to the heat diffusion up from the lower boundary until t=2500 where the heat sink is activated.

Matlab FEATool Nonlinear Time dependent Heat Transfer Solution

Regarding your second point with numerical source term, you could in this case create and call your own external function which tabulates and interpolates your data, in this case it would look something like

fea.phys.ht.eqn.coef{6,end}{1} = 'my_fun( t )';

where you have a Matlab function my_fun.m somewhere accessible on the Matlab paths of the form

function [ val ] = my_fun( t )
data  = [7 42 -100  0.1];   % Example data.
times = [0 10   99 5000];   % Time points.
val = interp1( data, times, t );   % Interpolate data.

Lastly, although the model is defined using m-script Matlab code here for sharing on StackOverflow you could just as easily use the Matlab FEA GUI, and even export your GUI model as m-script code if required.