PDE Toolbox: Using different temperature-dependent thermal conductivities

I'm playing around with the PDE toolbox and solving a problem with temperature dependent thermal conductivities. I followed the example from solving-a-heat-transfer-problem-with-temperature-dependent-properties. Using R2018b, I set up a simple two region geometry with only a temperature gradient:

device = createpde('thermal','steadystate');
r1 = [3 4 0 0.5e-6 0.5e-6 0 0 0 0.5e-6 0.5e-6];
r2 = [3 4 0 0.5e-6 0.5e-6 0 0.5e-6 0.5e-6 1e-6 1e-6];
gdm = [r1;r2]';
g = decsg(gdm, 'R1+R2',['R1';'R2']');
geometryFromEdges(device,g);

thermalBC(device, 'Edge', 2, 'Temperature', 500);  % 500K Dirichlet 
thermalBC(device, 'Edge', 1, 'Temperature', 300);  % 300K Dirichlet

generateMesh(device,'Hmax',1e-7);

Now, using a single thermal conductivity for both regions I get the expected result for constant and temperature dependent thermal conductivity.

% Case 1: constant, same thermal conductivities
k1 = 100;
k2 = 100;

% Case 2: T-dep same thermal conductivities
% k1 = @(~,state) 100*(state.u/300).^-1.3;
% k2 = @(~,state) 100*(state.u/300).^-1.3;
% thermalIC(device,400);

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

soln = solve(device);
temp = soln.Temperature;

The problem I'm running into is when using different thermal conductivities for each region as in:

% Case 1: constant, same thermal conductivities
k1 = 10;
k2 = 100;

% Case 2: T-dep same thermal conductivities
k1 = @(~,state) 10*(state.u/300).^-1.3;
k2 = @(~,state) 100*(state.u/300).^-1.3;
thermalIC(device,400);

In this situation, the solution I'm getting for Case 2 is the exact same as the previous simulation where I used 100 for each of the temperature-dependent thermal conductivity values. Am I specifying a property incorrectly perhaps?

Edited to include plot with solution I expect:

Answer 1

[Ravi Kumar]

Both k1 and k2 yields zero when state.u = 0, this is causing solver to misinterpret them as same functions. You can get around this by adding a small noise to one of the functions:

k1 = @(~,state) 10*(state.u/300).^-1.3;
k2 = @(~,state) 100*(state.u/300).^-1.3+eps;

[hooked-on-phonons]

I tried adding an eps of .001, .1, 1, and 100 all giving the same result. The state.u variable should represent the temperature here, which is always between 300 and 500 because of the boundary conditions, correct?

I checked the thermal properties, and they're assigned as the correct functions though.

>> findThermalProperties(device.MaterialProperties,'Face',1)

ans = 

  ThermalMaterialAssignment with properties:

             RegionType: 'face'
               RegionID: 1
    ThermalConductivity: @(~,state)10*(state.u/300).^-1.3
            MassDensity: []
           SpecificHeat: []

>> findThermalProperties(device.MaterialProperties,'Face',2)

ans = 

  ThermalMaterialAssignment with properties:

             RegionType: 'face'
               RegionID: 2
    ThermalConductivity: @(~,state)100*(state.u/300).^-1.3+100
            MassDensity: []
           SpecificHeat: []

No matter what, I seem to get the solution attached.

Answer 2

[Ravi Kumar]

I did not realize this in my initial response, both of your k values goes to inf when u=0. Try the following following instead, put the following code in a M-file and run.

device = createpde('thermal','steadystate');
r1 = [3 4 0 0.5e-6 0.5e-6 0 0 0 0.5e-6 0.5e-6];
r2 = [3 4 0 0.5e-6 0.5e-6 0 0.5e-6 0.5e-6 1e-6 1e-6];
gdm = [r1;r2]';
g = decsg(gdm, 'R1+R2',['R1';'R2']');
geometryFromEdges(device,g);
thermalBC(device, 'Edge', 2, 'Temperature', 500);  % 500K Dirichlet 
thermalBC(device, 'Edge', 1, 'Temperature', 300);  % 300K Dirichlet
generateMesh(device,'Hmax',1e-7);
%Case 1: constant, same thermal conductivities
k1 = 100;
k2 = 100;
thermalProperties(device, 'ThermalConductivity', 100);
% Case 2: T-dep same thermal conductivities
thermalIC(device,400);
thermalProperties(device, 'ThermalConductivity', @kFunc);
device.SolverOptions.ReportStatistics='on'
soln = solve(device);
temp = soln.Temperature;


y = linspace(0,1E-6,100);
x = 2E-7*ones(size(y));
iTemp = interpolateTemperature(soln,x,y);
plot(y,iTemp)
hold on
plot([0,1E-6],[300,500])
legend('nonlinear','linear')


function k = kFunc(region,state)

if region.subdomain == 1
     k = 10*(state.u/300).^-1.3;
else
     k = 100*(state.u/300).^-1.3;
end
    
end
    

Answer 3

[hooked-on-phonons]

Doing this I'm still getting the same result as before.

This is exactly what I'd expect if the thermal conductivies are the same in the two regions. But since in kFunc the values are 10 in the first subdomain and 100 in the second, I'd expect the nonlinear case to have 2 distinct regions in the temperature plot. In the linear case, this would be two different slopes.

Answer 4

[hooked-on-phonons]

Edited original question with better plot of the solution I'm expecting.