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

4 ビュー (過去 30 日間)
Patrick Fleith
Patrick Fleith 2017 年 6 月 19 日
編集済み: Patrick Fleith 2017 年 6 月 20 日
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 件のコメント
Alan Weiss
Alan Weiss 2017 年 6 月 19 日
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
Patrick Fleith 2017 年 6 月 19 日
I am using R2017a. Thank you in advance :)

サインインしてコメントする。

採用された回答

Alan Weiss
Alan Weiss 2017 年 6 月 19 日
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 件のコメント
Patrick Fleith
Patrick Fleith 2017 年 6 月 20 日
編集済み: Patrick Fleith 2017 年 6 月 20 日
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.

サインインしてコメントする。

その他の回答 (0 件)

Community Treasure Hunt

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

Start Hunting!

Translated by