MATLAB: Using a function handle for ThermalConductivity is not changing resulting temperature plot

differentialerrorfunctionhandleheatpartialPartial Differential Equation Toolboxpdeplotresulttemperaturewrong

When implementing the spatial dependency of thermal conductivity as follows:
k1=@(location,~) 1*location.x; 
k2=@(location,~) 20*location.x; 
thermalProperties(model,'Face',2,'ThermalConductivity',k1); 
thermalProperties(model,'Face',1,'ThermalConductivity',k2);
The results shows that the solver is only using 'k1' for both domains.
So you finally see the following plot, when using the code:
%%set geometry up and solve 
model = createpde('thermal','steadystate'); 
R1=[3 4 0 1 1 0 0 0 1 1]'; 
R2=[3 4 0.5 1 1 0.5 0.5 0.5 1 1]'; 
gm=[R1,R2]; 
sf='R1+R2'; 
ns=char('R1','R2') 
ns=ns'; 
g=decsg(gm,sf,ns); 
figure(1) 
pdegplot(g,'FaceLabels','on','EdgeLabels','on') 
%%create mesh 
geometryFromEdges(model,g) 
generateMesh(model,'Hmax',0.025); 
figure(2) 
pdemesh(model); 
%%set boundary conditions 
temperatureTop=@(location,~) 1*location.x 
thermalBC(model,'Edge',[1],'Temperature',0); 
thermalBC(model,'Edge',[7:8],'Temperature',temperatureTop); 
thermalBC(model,'Edge',[2],'HeatFlux',0); 
thermalBC(model,'Edge',[5:6],'HeatFlux',0); 
%Here you can find the function handles for the ThermalConductivity
k1=@(location,~) 1*location.x; 
k2=@(location,~) 20*location.x; 
thermalProperties(model,'Face',2,'ThermalConductivity',k2); 
thermalProperties(model,'Face',1,'ThermalConductivity',k1); 
%%solve model 
result = solve(model); 
temp = result.Temperature; 
[qx,qy] = result.evaluateHeatFlux; 
figure(3) 
pdeplot(model,'XYData',temp,'Contour','on'); 
set(gca,'box','on')

Best Answer

The Solver is misclassifying thermal conductivity as global assignment with k1=@(location,~) 1*location.x;. This is happening because at location.x = 0, which is the query point used to detect difference in conductivity, both k1 and k2 yield the same value of zero.

In the meantime a workaround is to add a small noise, like k2=@(location,~) 20*location.x + eps, this should give the correct results.

The code with the workaround looks like the following:

%% set geometry up and solve 
model = createpde('thermal','steadystate'); 
R1=[3 4 0 1 1 0 0 0 1 1]'; 
R2=[3 4 0.5 1 1 0.5 0.5 0.5 1 1]'; 
gm=[R1,R2]; 
sf='R1+R2'; 
ns=char('R1','R2') 
ns=ns'; 
g=decsg(gm,sf,ns); 
figure(1) 
pdegplot(g,'FaceLabels','on','EdgeLabels','on') 

%% create mesh 
geometryFromEdges(model,g) 
generateMesh(model,'Hmax',0.025); 
figure(2) 
pdemesh(model); 

%% set boundary conditions 
temperatureTop=@(location,~) 1*location.x 
thermalBC(model,'Edge',[1],'Temperature',0); 
thermalBC(model,'Edge',[7:8],'Temperature',temperatureTop); 
thermalBC(model,'Edge',[2],'HeatFlux',0); 
thermalBC(model,'Edge',[5:6],'HeatFlux',0); 

%Here you can find the function handles for the ThermalConductivity
k1=@(location,~) 1*location.x; 
k2=@(location,~) 20*location.x+eps; %here you see the fix with 'eps'

thermalProperties(model,'Face',2,'ThermalConductivity',k2); 
thermalProperties(model,'Face',1,'ThermalConductivity',k1); 

%% solve model 
result = solve(model); 
temp = result.Temperature; 

[qx,qy] = result.evaluateHeatFlux; 

figure(3) 
pdeplot(model,'XYData',temp,'Contour','on'); 
set(gca,'box','on')

This is the resulting plot:

Our development team is currently working on a fix for a future release. It seems to be that this will be fixed in one of the next releases after R2019a.

Related Question