Set a noncontinous 2d heat source in PDE toolbox

5 views (last 30 days)
hello,
i am trying to set a 2d heatsource in my simulation.
a bit more detail about the problem: i will solve a transient conduction problem of a rotation symmetric cylinder. The source terms are induced by radiation, and come from an external file, created by a different script. The source comes from a 2d matrix, that i want to map onto the PDE toolbox. however i cannot declare a mapping function for qFunc, as i get this error:
Error using pde.HeatSourceAssignment/checkCoefFcnHdlArgCounts (line 240)
Function specifying a heat source must accept two
input arguments and return one output argument.
Error in pde.HeatSourceAssignment/checkFcnHdlArgCounts (line 135)
self.checkCoefFcnHdlArgCounts(self.HeatSource, systemsize, ndims);
Error in pde.HeatSourceAssignment (line 62)
self.checkFcnHdlArgCounts(systemsize, numdims);
Error in pde.ThermalModel/internalHeatSource (line 99)
hs = pde.HeatSourceAssignment(hsContainer,argsToPass{:});
Error in Glass_melt_conduction_only (line 65)
internalHeatSource(thermalmodel,qFunc);
The problem seems to be, that the function is evaluated in the script, but the only thing that is passed on is @heatgeneration. I did it as follows:
...
kFunc = @(region,state) kFunc*region.x;
cFunc = @(region,state) SpecificHeat*region.x;
qFunc = @(region,state) heatgeneration(region)*region.x;
...
internalHeatSource(thermalmodel,qFunc);
...
%%funtion evaluation
function r = getglobal
global q
global Radius
global Height
global n_ring
global n_height
r.a =q;
r.b= Radius;
r.c= Height;
r.d= n_ring;
r.e= n_height;
end
function [result] = heatgeneration(region)
r=getglobal
x=region.x/r.b*r.d;
y=region.y/r.c*r.e;
q=r.a;
if mod(x,1)==0 %if region.x matches the segmentation
if mod(y,1)==0 %if region.y matches the segmentation
heat=q(y,x);
else %if region.y does not match the segmentation
y1=y-mod(y,1);
y2=y1+1;
q1=q(y1,x);
q2=q(y2,x);
heat=(q2-q1)/(y2-y1)*(y-y1) + q1;
end
else %if region.x does not match the segmentation
if mod(y,1)==0 %if region.y matches the segmentation
x1=x-mod(x,1);
x2=x1+1;
q1=q(y, x1);
q2=q(y, x2);
heat=(q2-q1)/(x2-x1)*(x-x1) + q1;
else %if region.y does not match the segmentation
y1=y-mod(y,1);
y2=y1+1;
x1=x-mod(x,1);
x2=x1+1;
q11=q(y1,x1);
q12=q(y2,x1);
q21=q(y1, x2);
q22=q(y2, x2);
q1=(q12-q11)/(y2-y1)*(y-y1) + q11;
q2=(q22-q21)/(y2-y1)*(y-y1) + q21;
heat=(q2-q1)/(x2-x1)*(x-x1) + q1;
end
end
result=heat;
end
i tried to use global variables, in order to account for the error message, and even tried to make the whole funtion heatgeneration=qfunc, that did not work either.

Answers (2)

Ravi Kumar
Ravi Kumar on 4 Oct 2018
As the message indicates, the function should take 2 input arguments (which your function does, no problem here). However, your function did not return one output argument. The scenario could be that your function definition did not return an output, it errored during execution. There are a couple of places where it might have errored:
1. In the definition: qFunc = @(region,state) heatgeneration(region)*region.x; the variable region.x would be of size 1 x N, where N can be 1 or a large number, meaning it is a rwo vector. The output of heatgeneration(region)*region.x should be of the same size as region.x, please check this first. That means heatgeneration(region) must return a scalar OR you have to use elementwise multiplication heatgeneration(region).*region.x if heatgeneration(region) returns a row vector of same size as that of region.x.
You can debug your heatgeneration(region) to make sure it works when region.x is a row vector. You can test this easily by creating your own region structure as:
region.x = rand(1,10); region.y = rand(1,10); outHeat = heatgeneration(region);
outHeat must be a scalar or a row vector of size 10, otherwise you have an error in the function.
Regards, Ravi
  4 Comments
Ravi Kumar
Ravi Kumar on 9 Oct 2018
I am missing something,. The working toy qFunc assigns size(region.x,1) heat flux values to size(region.x,1) points at which solver wants the heat flux value:
size(region.x,1) will be 1x1 for the initial few calls, where solver is determining the characteristics of the function, if you continue for few more call you will notice region.x will be of size 1xN, N is the order of your number of elements in the model. You may do this experiment by writing the full function as as mentioned at the end of my last post.
The requirement of qFunc is that it must return ONE the value of heat flux at each point point as defined by region.x and region.y vectors, hence heat flux output size should be 1xnumel(region.x).
Better yet, provide a complete reproduction code and sample data, then we can work from there.

Sign in to comment.


O Mueller
O Mueller on 10 Oct 2018
%%Solve glass melt heat conduction problem
clc;
close all;
clear all;
%%Input
Radius = 0.3; % [m]
Height = 0.0635; % [m]
n_rings = 10; % Number of rings
n_heights = 20; % Number of heights for which temperature is calculate (At least 2, bottom and top)
InitialTemperature = 300; % [K] input temperature
BottomTemperature = 300; % [K] input bottom temperature
SideTemperature = 300; % [K] input side temperature
TimeStep = 60; % [s] Timestep for one iteration
TransientTime = 3600; % [s] Time to calculate transient solution
SpecificHeat = 1464; % [J/kg/K] Specific heat capacity
MassDensity = 1202.5; % [kg/m^3] mass density of the glass melt
kFunc = 0.173; % [W/m/K]
MeshSize = 0.003; % [m] Size of the triangular mesh
HeatFlux = 867.5*ones(n_rings,1); % [W/m^2] Solar irradiance
%%End input section
%%Your input (now array with ones)
load('powerdensity');
q = PD2d;
%we declare now some global variables
setglobal(q,Radius,Height,n_rings,n_heights)
Create a transient thermal model and include the geometry.
thermalmodel = createpde('thermal','transient');
% Define rectangle with different input sources from the radiosity method
for j = 1:n_rings
xVector(j) = j*Radius/n_rings; % Define points on the glass melt surface
yVector(j) = zeros(length(n_rings));
end
g = decsg([2 3+n_rings xVector Radius 0 0 yVector -Height -Height 0]');
figure
pdegplot(g,'EdgeLabels','on','SubdomainLabels','on')
axis equal
% Convert the decsg format into a geometry object. Include the geometry in
% the model.
geometryFromEdges(thermalmodel,g);
% Define Thermal conductivity and specific heat depending on the radius of
% the rod
% data for k taken from pilon 2014 equation 25
kFunc = @(region,state) kFunc*region.x;
cFunc = @(region,state) SpecificHeat*region.x;
%qFunc = @heatgeneration;
% Specify thermal conductivity, mass density, and specific heat of the
% material.
thermalProperties(thermalmodel,'ThermalConductivity',kFunc,...
'MassDensity',MassDensity,...
'SpecificHeat',cFunc);
% Define internal heat source (Dein Input qFunc)
internalHeatSource(thermalmodel,qFunc);
% Define boundary conditions, fixed temperatures at the Bottom and side of
% the pot, the BC in the middle r=0 is by default adiabatic!
% thermalBC(thermalmodel,'Edge',1,'Temperature',BottomTemperature);
% thermalBC(thermalmodel,'Edge',2,'Temperature',SideTemperature);
% Create a mesh with elements no larger than MeshSize
msh = generateMesh(thermalmodel,'Hmax',MeshSize);
% Calculate the transient solution
tlist = 0:TimeStep:TransientTime;
thermalIC(thermalmodel,InitialTemperature);
R = solve(thermalmodel,tlist);
for j = 1:n_heights
tempdisk(j,:) = interpolateTemperature(R,xVector'-Radius/(2*n_rings),yVector'-(j-1)*Height/(n_heights-1),length(tlist):1:length(tlist));
Y(j,:) = yVector'-(j-1)*Height/(n_heights-1);
X(j,:) = xVector'-Radius/(2*n_rings);
end
save('Radius','X')
save('Height','Y')
save('tempdisk','tempdisk')
% Function to read heat out of array
function setglobal(q1,Radius1,Height1,n_rings1,n_heights1)
global q
global Rad
global Heig
global n_ri
global n_heig
q =q1;
Rad = Radius1;
Heig = Height1;
n_ri =n_rings1;
n_heig = n_heights1;
end
function r = getglobal
global q
global Rad
global Heig
global n_ri
global n_heig
r.a =q;
r.b= Rad;
r.c= Heig;
r.d= n_ri;
r.e= n_heig;
end
function [result] = qFunc(region, state) %heatgeneration(region, state)
r=getglobal
for p=1:size(region.x,2)
x=region.x(p)/r.b*r.d
y=region.y(p)/r.c*r.e
q=r.a;
y1=y-mod(y,1);
y2=y1+1;
x1=x-mod(x,1);
x2=x1+1;
if y1==0
y1=1;
end
if x1==0
x1=1;
end
if mod(x,1)==0 %if region.x matches the segmentation
if mod(y,1)==0 %if region.y matches the segmentation
heat(p)=q(y,x);
else %if region.y does not match the segmentation
q1=q(y1,x);
q2=q(y2,x);
heat(p)=(q2-q1)/(y2-y1)*(y-y1) + q1;
end
else %if region.x does not match the segmentation
if mod(y,1)==0 %if region.y matches the segmentation
q1=q(y, x1);
q2=q(y, x2);
heat(p)=(q2-q1)/(x2-x1)*(x-x1) + q1;
else %if region.y does not match the segmentation
q11=q(y1,x1);
q12=q(y2,x1);
q21=q(y1, x2);
q22=q(y2, x2);
q1=(q12-q11)/(y2-y1)*(y-y1) + q11;
q2=(q22-q21)/(y2-y1)*(y-y1) + q21;
heat(p)=(q2-q1)/(x2-x1)*(x-x1) + q1;
end
end
end
result=heat*region.x;
end
  1 Comment
Ravi Kumar
Ravi Kumar on 11 Oct 2018
Thanks. You have many coding errors in your qFunc. Before that first need to change the internalHeatSource specification as:
internalHeatSource(thermalmodel,@qFunc);
as qFunc is not a function_handle by itself in this version of your code.
Then insert a breakpoint on the first line of your qFunc and run the code. Step through to resolve all error. I tried to fix a few, but this is not possible without knowing what is your intent. Here are firs few error you will get:
1. heat(p) is assigned without initializing it. This can be fixed by assigning heat = zeros(size(region.x)), before its first use, beginning of the function.
2. Next you will get error because you try to access q(y,x) for x=y=0. Is did the same thing you have done for case x1=y1=0.
3. Then I got the error because q(y,x) is accessed with negative value of y.
Then I stopped. You need to continue debugging this function till the end.
Also at the end you need .* when you multiple heat and region.x as:
result=heat.*region.x;
You may need to consider a better approach to access your tabulated data. Please refer to the interpolation functions.
Regards, Ravi

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!