Non-constant parameters for PDE (thermal solver)

Dear Matlab experts, I have difficulties to use non-constant parameters for PDE solver beside simple case as in radioactiveRod example. In the following code, I want to solve the thermal equation of a rod which has a temperature dependent function as heat source. PDE complains that it supports only 2 input and 1 output. I had read how to pass extra parameters in ananymous function but did not really understand how to fix this.

Here is a code to improve. The faulty function is at the end.

Thanks for your help.

    % Rod model is 2D axissymmetric as in example radioactiveRod
  close all;
  clear all;

% Constant Thermal properties.
k = 0.293;  % thermal conductivity, W/(m-degree C)
rho = 1200; % density, kg/m^3
cp = 1800;  % specific heat, W-s/(kg-degree C)
Emissivity = 0.9;
Boltzmann = 5.670367e-8; % W?m?2?K?4
% others input are
T_ambient=25; %°C
Load=2.5;

% Non-constant Material properties
% this is what I want to be able to interpolate for each temperature during
% the solving of thermalModelS
global Mat_Temp  % I am trying if global declaration would help
global Mat_Prop

Mat_Temp=[1;2;3;4;5;6;7;8;9;10;11;12;13;14;15;16;17;18;19;20;21;22;23;24;25;26;27;28;29;30;31;32;33;34;35;36;37;38;39;40;41;42;43;44;45;46;47;48;49;50;51;52;53;54;55;56;57;58;59;60;61;62;63;64;65;66;67;68;69;70;71;72;73;74;75;76;77;78;79;80;81;82;83;84;85;86;87;88;89;90;91;92;93;94;95;96;97;98;99;100;101;102;103;104;105;106;107;108;109;110;111;112;113;114;115;116;117;118;119;120;121;122;123;124;125;126;127;128;129;130;131;132;133;134;135;136;137;138;139;140;141;142;143;144;145;146;147;148;149;150;151;152;153;154;155;156;157;158;159;160;161;162;163;164;165;166;167;168;169;170;171;172;173;174;175;176;177;178;179;180;181;182;183;184;185;186;187;188;189;190;191;192;193;194;195;196;197;198];
Mat_Prop=[30507032.3781981;29370504.5893094;28254938.3766849;27148809.6904823;25921908.3835235;24723352.8423951;23610251.6473241;22754300.9289475;21911800.5538344;21082750.5219850;20573831.9568911;20175585.1086136;19780594.3859345;19388859.7888540;19056059.6901305;18734575.2768598;18415580.8414250;18099076.3838257;17782498.8922778;17464831.0731487;17149728.6860157;16837191.7308793;16527220.2077390;16148158.5991148;15811583.4819250;15477915.3724935;15145345.6139273;14845088.1332907;14575384.1924548;14306558.9323080;14011179.5967969;13768006.6694989;13587918.8310611;13409016.4135237;13231299.4168870;13054767.8411508;12879421.6863151;12714959.9814305;12584003.9417475;12453630.4853288;12323839.6121743;12214763.1034955;12139301.8799586;12064073.6895734;11989078.5323398;11914316.4082580;11837384.9693072;11749545.2261897;11560161.5337132;11371838.0172300;11184574.6767404;10998371.5122442;10865280.8089144;10772309.3648469;10679616.1467752;10587201.1546991;10495064.3886187;10380322.1052734;10254705.9628154;10129500.9765798;9999793.59395625;9844605.73478730;9680327.26311754;9516776.98220725;9354401.16192604;9199258.82696638;9059758.30613417;8930702.28421105;8825033.97706680;8719664.32807701;8615464.37054548;8525790.02263360;8436364.55642918;8347187.97193222;8257843.65139554;8165640.31803141;8073750.64808374;7947468.53604545;7811583.57038311;7676200.46636492;7541319.22399095;7411723.33286573;7277873.14618985;7141423.21462658;7005726.07164403;6876389.56859266;6783235.91929127;6694856.93983988;6607093.63451095;6519710.16215555;6432706.52277371;6346082.71636543;6235138.12578329;6107706.90929559;5982822.81021327;5859038.70319387;5735772.87110252;5626300.07650937;5541237.23592349;5459965.92449659;5379905.61850384;5300070.58351147;5220460.81951947;5141076.32652786;5065829.01038247;4993366.93321174;4921012.82889513;4848766.69743259;4776628.53882415;4707558.98356207;4647441.96560696;4589492.10770651;4531542.24980606;4473592.39190560;4413985.77138545;4346944.20027075;4280145.01590968;4213588.21830226;4146610.19485479;4076178.10336159;4008754.40675467;3943605.43706652;3878567.37373352;3813640.21675565;3752239.21944086;3711900.88718129;3671562.55492171;3631224.22266214;3590885.89040256;3550547.55814298;3510209.22588341;3469870.89362383;3429532.56136426;3389194.22910468;3321882.29009508;3250361.87646869;3179019.41555779;3107854.90736235;3036868.35188237;2966059.74911784;2892869.10942723;2843087.50679328;2810229.09327789;2777551.92220668;2745055.99357964;2716282.07794596;2687967.55572282;2678055.45524065;2679424.20587832;2680730.63626793;2678792.06724812;2676497.75129149;2674115.96981452;2670620.80931011;2667524.63649561;2665306.39485724;2663028.74605799;2658439.81892946;2650829.31001418;2643129.00938759;2633297.91359214;2615822.73979788;2595036.17058021;2574249.60136255;2553463.03214490;2534100.33665622;2521629.03327273;2509048.33493049;2496358.24162944;2478542.93897787;2459963.93231009;2441244.26685284;2425388.87642435;2426950.00298769;2445830.16476708;2450690.91147399;2453981.69594784;2456238.54444371;2467225.93989086;2479340.35913718;2480912.13441168;2480842.28211890;2470104.28915241;2455868.32797627;2440902.96862651;2438323.40876118;2435577.99452272;2431242.47701829;2423835.23663984;2415716.91075331;2407193.27186701;2396035.40795831;2383634.83372078;2367258.54102386;2338954.47682607;2302154.45384933;2270010.87851920;2237867.30318908];

% Create a thermal model for steady-state analysis.
thermalModelS = createpde('thermal','steadystate');
thermalModelS.StefanBoltzmannConstant = 5.670373E-8;

% The rod is defined as a 2-D strip model as shown below.
g = decsg([3 4 -0.015 0.015 0.015 -0.015 0 0 .0125 .0125]');

% Convert the geometry and append it to the thermal model.
geometryFromEdges(thermalModelS,g);

PDE Toolbox allows definition of the non-constant coefficients as function of spatial coordinates and solution.

kFunc = @(region,state) k*region.y;
cFunc = @(region,state) cp*region.y;

% This is where it differs from the built-in example radioactiveRod
% the heat source is generated temperature dependent
qFunc = @(region,state) heatgeneration(state.temperature)*region.y ;

% For a steady-state analysis, specify the thermal conductivity of the
% material.
thermalProperties(thermalModelS,'ThermalConductivity',kFunc);

% Specify internal heat source.
internalHeatSource(thermalModelS,qFunc,'Face',1);

% When defining boundary conditions below, it is necessary to know the
% edge numbers for the boundary edges of the geometry. A convenient
% way to obtain these edge numbers is to plot the geometry using
% |pdegplot| with option |edgeLabels| set to |'on'|.
figure
pdegplot(thermalModelS,'EdgeLabels','on');
axis equal
xlim([-0.1 0.1]);
title 'Rod Section Geometry With Edge Labels Displayed';

% Define the boundary conditions assuming constant temperature on all edges
% but the axis of symmetry.
thermalBC(thermalModelS,'Edge',[2 3 4],'Temperature',T_ambient);

Generate the mesh.

generateMesh(thermalModelS,'Hmax',0.001);
figure; 
pdeplot(thermalModelS); 
axis equal
title 'Rod Section With Triangular Element Mesh'

Solve the model and plot the result.

result = solve(thermalModelS);
T = result.Temperature;
figure; 
pdeplot(thermalModelS,'XYData',T,'Contour','on','ColorMap','jet'); 
axis equal
title 'Steady State Temperature';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function heat = heatgeneration(temperature)
    ncProp=NonConstantProp(temperature);
    heat =0.1571*ncProp;    
    % Nested function to pick up interpolated value of mat properties
    function R = NonConstantProp(~)        
        R = interp1(Mat_Temp,Mat_Prop, temperature);     
    end
end

Answer 1

[JClarcq]

This seems to be the right solution

% Set initial condition to ambient everywhere
thermalIC(thermalModelS,T_ambient);

% the heat source is generated temperature dependent
qFunc = @(region,state) heatgeneration(region,state,Mat_Temp,Mat_Prop);

% Specify internal heat source.
internalHeatSource(thermalModelS,qFunc,'Face',1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function heat = heatgeneration(region,state,MatTemp,MatProp)
N = 1;                  % Number of equations
nr = length(region.y);  % Number of columns
heat = zeros(N,nr);     % Allocate to avoid size issue
try
    ncProp=interp1(MatTemp,MatProp, state.u);
    tmp =0.1571.*ncProp;
    heat=tmp.*region.y;
    if isnan(heat)
        heat(N,nr)=0;
    end
    catch
        heat(N,nr)=0;
    end
end

Answer 2

[JClarcq]

I tried also this without success:

% Rod model is 2D axissymmetric as in example radioactiveRod
close all;
clear all;

% Constant Thermal properties.
k = 0.293;  % thermal conductivity, W/(m-degree C)
rho = 1200; % density, kg/m^3
cp = 1800;  % specific heat, W-s/(kg-degree C)
Emissivity = 0.9;
Boltzmann = 5.670367e-8; % W?m?2?K?4
% others input are
T_ambient=25; %°C
Load=2.5;

% Non-constant Material properties
% this is what I want to be able to interpolate for each temperature during
% the solving of thermalModelS
global Mat_Temp  % I am trying if global declaration would help
global Mat_Prop

Mat_Temp=[1;2;3;4;5;6;7;8;9;10;11;12;13;14;15;16;17;18;19;20;21;22;23;24;25;26;27;28;29;30;31;32;33;34;35;36;37;38;39;40;41;42;43;44;45;46;47;48;49;50;51;52;53;54;55;56;57;58;59;60;61;62;63;64;65;66;67;68;69;70;71;72;73;74;75;76;77;78;79;80;81;82;83;84;85;86;87;88;89;90;91;92;93;94;95;96;97;98;99;100;101;102;103;104;105;106;107;108;109;110;111;112;113;114;115;116;117;118;119;120;121;122;123;124;125;126;127;128;129;130;131;132;133;134;135;136;137;138;139;140;141;142;143;144;145;146;147;148;149;150;151;152;153;154;155;156;157;158;159;160;161;162;163;164;165;166;167;168;169;170;171;172;173;174;175;176;177;178;179;180;181;182;183;184;185;186;187;188;189;190;191;192;193;194;195;196;197;198];
Mat_Prop=[30507032.3781981;29370504.5893094;28254938.3766849;27148809.6904823;25921908.3835235;24723352.8423951;23610251.6473241;22754300.9289475;21911800.5538344;21082750.5219850;20573831.9568911;20175585.1086136;19780594.3859345;19388859.7888540;19056059.6901305;18734575.2768598;18415580.8414250;18099076.3838257;17782498.8922778;17464831.0731487;17149728.6860157;16837191.7308793;16527220.2077390;16148158.5991148;15811583.4819250;15477915.3724935;15145345.6139273;14845088.1332907;14575384.1924548;14306558.9323080;14011179.5967969;13768006.6694989;13587918.8310611;13409016.4135237;13231299.4168870;13054767.8411508;12879421.6863151;12714959.9814305;12584003.9417475;12453630.4853288;12323839.6121743;12214763.1034955;12139301.8799586;12064073.6895734;11989078.5323398;11914316.4082580;11837384.9693072;11749545.2261897;11560161.5337132;11371838.0172300;11184574.6767404;10998371.5122442;10865280.8089144;10772309.3648469;10679616.1467752;10587201.1546991;10495064.3886187;10380322.1052734;10254705.9628154;10129500.9765798;9999793.59395625;9844605.73478730;9680327.26311754;9516776.98220725;9354401.16192604;9199258.82696638;9059758.30613417;8930702.28421105;8825033.97706680;8719664.32807701;8615464.37054548;8525790.02263360;8436364.55642918;8347187.97193222;8257843.65139554;8165640.31803141;8073750.64808374;7947468.53604545;7811583.57038311;7676200.46636492;7541319.22399095;7411723.33286573;7277873.14618985;7141423.21462658;7005726.07164403;6876389.56859266;6783235.91929127;6694856.93983988;6607093.63451095;6519710.16215555;6432706.52277371;6346082.71636543;6235138.12578329;6107706.90929559;5982822.81021327;5859038.70319387;5735772.87110252;5626300.07650937;5541237.23592349;5459965.92449659;5379905.61850384;5300070.58351147;5220460.81951947;5141076.32652786;5065829.01038247;4993366.93321174;4921012.82889513;4848766.69743259;4776628.53882415;4707558.98356207;4647441.96560696;4589492.10770651;4531542.24980606;4473592.39190560;4413985.77138545;4346944.20027075;4280145.01590968;4213588.21830226;4146610.19485479;4076178.10336159;4008754.40675467;3943605.43706652;3878567.37373352;3813640.21675565;3752239.21944086;3711900.88718129;3671562.55492171;3631224.22266214;3590885.89040256;3550547.55814298;3510209.22588341;3469870.89362383;3429532.56136426;3389194.22910468;3321882.29009508;3250361.87646869;3179019.41555779;3107854.90736235;3036868.35188237;2966059.74911784;2892869.10942723;2843087.50679328;2810229.09327789;2777551.92220668;2745055.99357964;2716282.07794596;2687967.55572282;2678055.45524065;2679424.20587832;2680730.63626793;2678792.06724812;2676497.75129149;2674115.96981452;2670620.80931011;2667524.63649561;2665306.39485724;2663028.74605799;2658439.81892946;2650829.31001418;2643129.00938759;2633297.91359214;2615822.73979788;2595036.17058021;2574249.60136255;2553463.03214490;2534100.33665622;2521629.03327273;2509048.33493049;2496358.24162944;2478542.93897787;2459963.93231009;2441244.26685284;2425388.87642435;2426950.00298769;2445830.16476708;2450690.91147399;2453981.69594784;2456238.54444371;2467225.93989086;2479340.35913718;2480912.13441168;2480842.28211890;2470104.28915241;2455868.32797627;2440902.96862651;2438323.40876118;2435577.99452272;2431242.47701829;2423835.23663984;2415716.91075331;2407193.27186701;2396035.40795831;2383634.83372078;2367258.54102386;2338954.47682607;2302154.45384933;2270010.87851920;2237867.30318908];

% Create a thermal model for steady-state analysis.
thermalModelS = createpde('thermal','steadystate');
thermalModelS.StefanBoltzmannConstant = 5.670373E-8;

% The rod is defined as a 2-D strip model as shown below.
g = decsg([3 4 -0.015 0.015 0.015 -0.015 0 0 .0125 .0125]');

% Convert the geometry and append it to the thermal model.
geometryFromEdges(thermalModelS,g);

%InitialConditions(thermalModelS,T_ambient,'Face',1);

PDE Toolbox allows definition of the non-constant coefficients as function of spatial coordinates and solution.

kFunc = @(region,state) k*region.y;
cFunc = @(region,state) cp*region.y;

% This is where it differs from the built-in example radioactiveRod
% the heat source is generated temperature dependent
qFunc = @(region,state) heatgeneration(region,state);

% For a steady-state analysis, specify the thermal conductivity of the
% material.
thermalProperties(thermalModelS,'ThermalConductivity',kFunc);

% Specify internal heat source.
internalHeatSource(thermalModelS,qFunc,'Face',1);

% When defining boundary conditions below, it is necessary to know the
% edge numbers for the boundary edges of the geometry. A convenient
% way to obtain these edge numbers is to plot the geometry using
% |pdegplot| with option |edgeLabels| set to |'on'|.
figure
pdegplot(thermalModelS,'EdgeLabels','on');
axis equal
xlim([-0.1 0.1]);
title 'Rod Section Geometry With Edge Labels Displayed';

% Define the boundary conditions assuming constant temperature on all edges
% but the axis of symmetry.
thermalBC(thermalModelS,'Edge',[2 3 4],'Temperature',T_ambient);

Generate the mesh.

generateMesh(thermalModelS,'Hmax',0.001);
figure; 
pdeplot(thermalModelS); 
axis equal
title 'Rod Section With Triangular Element Mesh'

Solve the model and plot the result.

result = solve(thermalModelS);
T = result.Temperature;
figure; 
pdeplot(thermalModelS,'XYData',T,'Contour','on','ColorMap','jet'); 
axis equal
title 'Steady State Temperature';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function heat = heatgeneration(region,state)
    ncProp=interp1(Mat_Temp,Mat_Prop, state.u);
    tmp =0.1571.*ncProp;
    heat=tmp.*region.y;
end