Cannot apply multiple, region-dependent thermal properties in PDE Toolbox

I am trying to solve a heat transfer problem on a 2D-geometry with multiple domains. I want to assign different thermal properties to the different faces of the geometry using the function 'thermalproperties'. The problem is axis-symmetric around the y-axis, and the suggested way to implement this is to multiply the material properties with 'r', or region.y/region.x, as explained here: https://nl.mathworks.com/help/pde/ug/heat-distribution-in-a-circular-cylindrical-rod.html.

This is how I tried applying the above (k1, k2, c1, c2 are constants):

kfun1 = @(region, state) k1*region.x;
kfun2 = @(region, state) k2*region.x;
cfun1 = @(region, state) c1*region.x;
cfun2 = @(region, state) c2*region.x;

thermalProperties(thermalmodel,'ThermalConductivity',kfun1, 'MassDensity',p1, 'SpecificHeat',cfun2,'Face',1);   
thermalProperties(thermalmodel,'ThermalConductivity',kfun2, 'MassDensity',p2, 'SpecificHeat',cfun1,'Face',2);  

The problem I run into is that when I use the above method, instead of assigning the thermal properties to the corresponding faces, the last listed properties are applied to all faces. In the above example, kfun2 and cfun2 are applied on the entire geometry. And when the thermal properties for face 1 are listed last, kfun1 and cfun1 are applied everywhere. Same goes for more than 2 faces, the last listed properties are used. This does not happen when the factor *region.x is left out, implying that the multiplication with the x-coordinate causes the properties to be applied across the geometry. Is there any way to cirumvent this problem?

Answer 1

[Ravi Kumar]

This is a bug where solver fails to distinguish the functions as seperate based on some queries. Workaround is to define variation within a single function using region.subdomain to return different values.

Change the two thermalProperties call to a single one:

thermalProperties(thermalmodel,'ThermalConductivity',@kfun, 'MassDensity',@mfunc, 'SpecificHeat',@cfun);   

Where, functions are defined to give differnt values of properites:

function k = kfun(region,state)
k1 = 1;
k2 = 2;

if region.subdomain == 1
     k = k1*region.x;
else
     k = k2*region.x;
end
    
end
   
function c = cfun(region,state)
c1 = 10;
c2 = 20;
if region.subdomain == 1
     c = c1*region.x;
else
     c = c2*region.x;
end
    
end

   
function p = mfun(region,state)
p1 = 100;
p2 = 200;
if region.subdomain == 1
     p = p1*ones(size(region.x));
else
     p = p2*ones(size(region.x));
end
    
end

Regards,

Ravi

[Burak Demirbas]

Dear Ravi,

Thank you for the detailed response. I ran the code you posted. Unfortunately, the second value was chosen (k2, c2, p2) everytime for the entire geometry again. I then changed

function k = kfun(region,state)
k1 = 1;
k2 = 2;
if region.subdomain == 1
     k = k1*region.x;
else
     k = k2*region.x;
end
   

to

function k = kfun(region,state)
k1 = 1;
k2 = 2;
if region.subdomain == 2
     k = k2*region.x;
else
     k = k1*region.x;
end
   

Which does work, the values appear to be properly assigned now. But I don't understand why this should make a difference?