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);figurepdegplot(model,'EdgeLabels','on','FaceLabels','on')xlim([-0.1 0.4])ylim([-0.1 0.6])axis equal;figurepdeplot(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;
Best Answer