Non-constant parameters for PDE (thermal solver)

14 views (last 30 days)
JClarcq
JClarcq on 19 Feb 2018
Answered: JClarcq on 20 Feb 2018
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

Sign in to answer this question.

Accepted Answer

JClarcq
JClarcq on 20 Feb 2018
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

More Answers (1)

JClarcq
JClarcq on 20 Feb 2018
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

Categories

Find more on Heat Transfer in Help Center and File Exchange

Tags

Products

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!