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

9 views (last 30 days)
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;
  2 Comments
Alan Weiss
Alan Weiss on 19 Jun 2017
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

Sign in to comment.

Accepted Answer

Alan Weiss
Alan Weiss on 19 Jun 2017
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
  1 Comment
Patrick Fleith
Patrick Fleith on 20 Jun 2017
Edited: Patrick Fleith on 20 Jun 2017
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.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!