PDE Toolbox - How to specify coefficient for temperature-dependent properties ?

Hi all,

I simulate the heating of a Thermal Mass (TM) under a given heat flux. I am interested in the temperature profile inside the thermal mass in transient state. For that, I use the PDE toolbox of MATLAB.

The problem I have is that I didn't succeed to implement non-constant coefficients to the model. Especially, I would like to implement temperature-dependent properties such as the thermal conductivity of the thermal mass (denoted by k_TM in the code).

The following code works when the conductivity is a constant, but not when I try to make it temperature-dependent.

Constant Properties

k_TM=2.1; %thermal conductivity
density_TM=3000;
cp_TM=800; %specific heat

Creation of the Model and Geometry

numberOfPDE=1;
model=createpde(numberOfPDE);
TM = [2; 4; 0; 0.3; 0.3; 0; 0; 0; 0.5; 0.5];
gd = [TM];
sf = 'TM';
ns = char('TM');
ns = ns';
dl = decsg(gd,sf,ns);
geometryFromEdges(model,dl);
generateMesh(model,'Hmax',0.03);

figure
pdegplot(model,'EdgeLabels','on','FaceLabels','on')
xlim([-0.1 0.4])
ylim([-0.1 0.6])
axis equal;
figure
pdeplot(model)
xlim([-0.1 0.4])
ylim([-0.1 0.6])
axis equal;

Initial and boundary conditions

%Initial Conditions
u0=254;
setInitialConditions(model,u0);

%Boundary conditions
applyBoundaryCondition(model,'dirichlet','Edge',1,'u',254);

heat_gain=20e3; 
applyBoundaryCondition(model,'neumann','Edge',3,'g',heat_gain);

Specify Coefficient: Case1: k_TM is constant

%specify coefficients
specifyCoefficients(model,'m',0,'d',density_TM*cp_TM,'c',k_TM,'a',0,'f',0,'face',1)

Case2: k_TM is temperature dependent (it doesn't work) How shall I write it to ensure that it works ?

% k_TM=@(~,state) 6e-7*state.u*state.u-0.0028*state.u+3.3753;
% specifyCoefficients(model,'m',0,'d',density_TM*cp_TM,'c',k_TM,'a',0,'f',0,'face',1)

SIMULATION and PLOT

%SIMULATION
S=[];
t1=0;
tstep=60;
tend=3600;
time=t1:tstep:tend;
results = solvepde(model,time);
S = results.NodalSolution;
%query points to compute the TM block mean temperature
Xtm=[0 0.075 0.15 0.225 0.3 0 0.075 0.15 0.225 0.3 0 0.075 0.15 0.225 0.3 0 0.075 0.15 0.225 0.3 0 0.075 0.15 0.225 0.3 0 0.075 0.15 0.225 0.3]
Ytm=[0 0 0 0 0 0.1 0.1 0.1 0.1 0.1 0.2 0.2 0.2 0.2 0.2 0.3 0.3 0.3 0.3 0.3 0.4 0.4 0.4 0.4 0.4 0.5 0.5 0.5 0.5 0.5]

Tintrptm = interpolateSolution(results,Xtm,Ytm,1:length(time));
tm_temp=mean(Tintrptm);
TM_temp=tm_temp(:,end);

fig=figure(3);
pdeplot(model,'XYData',S(:,tend/tstep),'Contour','off','ColorMap','hot');
xlim([-0.1 0.4])
ylim([-0.1 0.6])
axis equal;

[Alan Weiss]

PDE Toolbox is changing rapidly. Please tell us your toolbox version (such as R2016b) so we can give you relevant help.

Alan Weiss

MATLAB mathematical toolbox documentation

[Patrick Fleith]

I am using R2017a. Thank you in advance :)

Answer 1

[Alan Weiss]

I believe that the error comes from a misunderstanding on your part: the state.u field is a row vector. So when you want to square or multiply elements in that vector, you need to use .* or .^2, not * or ^2.

% You wrote
k_TM=@(~,state) 6e-7*state.u*state.u-0.0028*state.u+3.3753;
% It should be
k_TM = @(~,state) 6e-7*state.u.*state.u - 0.0028*state.u + 3.3753;
% or
k_TM = @(~,state) 6e-7*state.u.^2 - 0.0028*state.u + 3.3753;

Alan Weiss

MATLAB mathematical toolbox documentation

[Patrick Fleith]

You are right. It works now ! I cannot believe I did that mistake. This is very useful for me, thank you Alan. Have a good day.