PDE Toolbox - How to specify coefficient for temperature-dependent properties ?
9 views (last 30 days)
Show older comments
Patrick Fleith
on 19 Jun 2017
Edited: Patrick Fleith
on 20 Jun 2017
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
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
Accepted Answer
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
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!