PDE Toolbox - How to specify coefficient for temperature-dependent properties ?
4 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Patrick Fleith
am 19 Jun. 2017
Bearbeitet: Patrick Fleith
am 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 Kommentare
Alan Weiss
am 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
Akzeptierte Antwort
Alan Weiss
am 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 Kommentar
Weitere Antworten (0)
Siehe auch
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!