# Time-dependent internal heat source

9 views (last 30 days)
Giulia Ulpiani on 1 Jul 2019
Answered: Ravi Kumar on 2 Jul 2019
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));

Ravi Kumar on 2 Jul 2019
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.