Time-dependent internal heat source

I am writing a code using Matlab pde toolbox. I need to include an internal heat source that is time dependent... I am attaching the code. I get this error "Undefined operator '+' for input arguments of type 'function_handle'.

Error in testthermal (line 43)

specifyCoefficients(model,'m',0,'d',d,'c',c,'a',a,'f',f+q);"

What I am doing wrong? Please help (I have highlighted the heat source part in bold letters)

k = 14; % thermal conductivity of Ni, W/(m-K)

rho = 8020; % density of NiTi, kg/m^3

specificHeat = 450; % specific heat of NiTi, J/(kg-K)

thick = .0004; % foil thickness in meters

stefanBoltz = 5.670373e-8; % Stefan-Boltzmann constant, W/(m^2-K^4)

hCoeff = 1; % Convection coefficient, W/(m^2-K)

% The ambient temperature is assumed to be 25 degrees-Celsius.

ta = 298.15;

emiss = .16; % emissivity of sma (average A and M)

numberOfPDE = 1; %returns a PDE model object for a system of 1 equations

model = createpde(numberOfPDE);

width = 0.02;

height = 0.003;

gdm = [3 4 0 width width 0 0 0 height height]'; %Define the square by giving the 4 x-locations followed by the 4 y-locations of the corners.

g = decsg(gdm, 'S1', ('S1')');

geometryFromEdges(model,g);

figure;

pdegplot(model,'EdgeLabels','on');

axis([-.01 .03 -.0075 .01]);

title 'Geometry With Edge Labels Displayed';

%mesh

hmax = .0004; % element size

msh = generateMesh(model,'Hmax',hmax);

figure;

pdeplot(model);

axis equal

title 'Plate With Triangular Element Mesh'

xlabel 'X-coordinate, meters'

ylabel 'Y-coordinate, meters'

%%Specify the coefficients. The expressions for the coefficients required by PDE Toolbox can easily be identified by comparing the equation above with the scalar parabolic equation in the PDE Toolbox documentation.

endTime = 615;

tlist = 0:.1:endTime;

q = @(~,state)7200/endTime*state.t; %latent heat is a simple linear function of time

internalHeatSource(model,q)

a = @(~,state) 2*hCoeff + 2*emiss*stefanBoltz*state.u.^3;

f = 2*hCoeff*ta + 2*emiss*stefanBoltz*ta^4;

d = thick*rho*specificHeat;

specifyCoefficients(model,'m',0,'d',d,'c',c,'a',a,'f',f+q);

p = msh.Nodes;

numNodes = size(p,2);

u0(1:numNodes) = ta; % Set the initial temperature of all nodes to ambient.

applyBoundaryCondition(model,'dirichlet','Edge',2,'u',ta);

applyBoundaryCondition(model,'dirichlet','Edge',4,'u',ta);

setInitialConditions(model,0);

model.SolverOptions.RelativeTolerance = 1.0e-3;

model.SolverOptions.AbsoluteTolerance = 1.0e-4;

R = solvepde(model,tlist);

u = R.NodalSolution;

figure;

plot(tlist,u(3, :));

grid on

title 'Temperature Along the Top Edge of the Plate as a Function of Time'

xlabel 'Time, seconds'

ylabel 'Temperature, degrees-Kelvin'

figure;

pdeplot(model,'XYData',u(:,end),'Contour','on','ColorMap','jet');

title(sprintf('Temperature In The Plate, Transient Solution( %d seconds)\n', ...

tlist(1,end)));

xlabel 'X-coordinate, meters'

ylabel 'Y-coordinate, meters'

axis equal;

fprintf('\nTemperature at the top edge(t=%5.1f secs)=%5.1f degrees-K\n', ...

tlist(1,end), u(4,end));

Answer 1

[Ravi Kumar]

You cannot add function handle as you have done, f+q in:

specifyCoefficients(model,'m',0,'d',d,'c',c,'a',a,'f',f+q);

I would recommand using ThemalModel workflow, using model = createpde('thermal'). Take a look at the examples.