PDETool temperature gradient-dependent internal heat generation

Hi,

I am currently working on implementing advection into a heat transfer problem. I have a prescribed velocity field and I want to transform the internal heat generation term to the following form: "-u*dT/dx - v*dT/dy" to address the advection portion. For that, I am imposing a function for the heat generation coefficient that depends on the derivatives of temperature at a specific location. However, when I run with a heat generation function (Heatsource) that includes "state.ux" or "state.uy", I get the following error:

"Dot indexing is not supported for variables of this type.

Error in heatsourcefunction (line 2)

heat=-(state.ux*y_vel+state.u*x_vel);"

Can anyone help me with this?

I could not find any question answering this problem, only addressing how to get a heat source dependent on location. However, I have tried changing the function Heatsource to be only dependent on the input "location", and the code runs succesfully. So, it seems like this is a different issue.

I have copied the example code below. In this code, I am only imposing the advection part in the subdomain "5".

Thanks,

Tadeu

function Tnest=Test_advection_example


Ic=700;
k0=0.2;
wspeed=1.38;
q1=50;%W

thermalmodelT = createpde('thermal','steady');

R1=[2;46;-1;-0.104789030466194;0.104789030466194;1;0.996040213058802;0.986926578126287;0.969523477663727;0.967940634753745;0.953654605292743;0.933658304681164;0.907578665717915;0.881447429010727;0.872890456587088;0.849919226015183;0.816099772602718;0.800322812121640;0.766481564083305;0.736857226024835;0.725801863401366;0.710526363731541;0.678896188331572;0.668552573056067;0.661806880421772;0.634160268808002;0.324859166248063;-0.324859166248063;-0.545581690922761;-0.561887496043782;-0.599192341904830;-0.629367284642012;-0.663403654417232;-0.674248777133431;-0.699572367663238;-0.711394965175357;-0.716586945532947;-0.734106103145198;-0.762861059912315;-0.800345857506185;-0.820123790995309;-0.834147953469543;-0.861743943723934;-0.894014431614409;-0.926306740272103;-0.961062587759775;-0.961312747639009;-0.965491972145012;0;0;0;0;0.103002254301223;0.206004508602445;0.309006762903668;0.412009017204890;0.515011271506113;0.618013525807336;0.721015780108558;0.824018034409781;0.927020288711003;1.03002254301223;1.13302479731345;1.23602705161467;1.33902930591589;1.44203156021712;1.54503381451834;1.64803606881956;1.75103832312078;1.85404057742201;1.95704283172323;2.06004508602445;2.16304734032567;2.13275626747832;2.03119644521745;1.92963662295658;1.82807680069570;1.72651697843483;1.62495715617396;1.52339733391309;1.42183751165221;1.32027768939134;1.21871786713047;1.11715804486960;1.01559822260872;0.914038400347852;0.812478578086980;0.710918755826107;0.609358933565235;0.507799111304362;0.406239289043490;0.304679466782617;0.203119644521745;0.101559822260872];
R2=[2;12;-0.104789030466194;0.104789030466194;0.104789030466194;0.953654605292743;0.933658304681164;0.104789030466194;0.104789030466194;-0.104789030466194;-0.104789030466194;-0.861743943723934;-0.894014431614409;-0.104789030466194;0;0;0.515011271506113;0.515011271506113;0.618013525807336;0.618013525807336;2.15278726271807;2.14301634508593;0.609358933565235;0.609358933565235;0.507799111304362;0.507799111304362];

R2new = [R2;zeros(length(R1) - length(R2),1)];
gm = [R1 R2new];
sf = 'R1+R2';

ns = char('R1','R2');
ns = ns';
g = decsg(gm,sf,ns);

geometryFromEdges(thermalmodelT,g);

mesh=generateMesh(thermalmodelT,'Hmax',0.05);

%FIND FACE ID OF PATHWAY
hcp_index=5;
for ii=1:length(g)
    surface_angle(ii)=atan((g(5,ii)-g(4,ii))/(g(3,ii)-g(2,ii))); %%angle from vertical
end


%CALCULATE ANGLE OF EACH SURFACE
surface_angle=zeros(1,length(g));

% % % % % % % % PDE COEFFICIENTS

%%%MATERIAL PROPERTIES
for ii=1:thermalmodelT.Geometry.NumFaces
    thermalProperties(thermalmodelT,'Face',ii,'ThermalConductivity',0.2);
end
thermalProperties(thermalmodelT,'Face',hcp_index,'ThermalConductivity',0.05);


%%%%%ADVECTION TERM FOR HCP

internalHeatSource(thermalmodelT,0);

x_vel=0;
y_vel=0.05;
Heatsource=@(location,state)heatsourcefunction(location,state,y_vel,x_vel);
internalHeatSource(thermalmodelT,Heatsource,'Face',hcp_index);



%%%%%% BOUNDARY CONDITIONS

Tamb=(21.38+33.46)/2;

for ii=1:length(g)
    if g(4,ii) == 0 && g(5,ii) == 0
        if abs(g(2,ii) + 1) < 1e-6 || abs(g(3,ii) - 1) < 1e-6
            thermalBC(thermalmodelT,'Edge',ii,'HeatFlux',0); %bottom faces
        else
            thermalBC(thermalmodelT,'Edge',ii,'HeatFlux',q1); %nest
        end
    else
        thermalBC(thermalmodelT,'Edge',ii,'ConvectionCoefficient',...
            12.5+(wspeed*max(0,cos(surface_angle(ii)+pi/2))),...
            'AmbientTemperature',Tamb,...
            'HeatFlux',max(0,-Ic*cos(surface_angle(ii)+pi/2)));
    end
end

% % % % % % % % SOLVING THERMAL PDE
R = solve(thermalmodelT);
T = R.Temperature;


Tsum=0;
count=0;
for ii=1:length(mesh.Nodes(2,:))
    if abs(mesh.Nodes(2,ii)) < 0.2 && abs(mesh.Nodes(1,ii)) < 0.1
        Tsum=Tsum+T(ii,end);
        count=count+1;
    end
end
Tnest=Tsum/count;
Tave=mean(T);
Tmin=min(T);
Tmax=max(T);

fig3=figure(3);
pdeplot(thermalmodelT,'XYData',T(:,end),'Contour','off','ColorMap','jet');
caxis([20 50])

end

Answer 1

[Ravi Kumar]

Hi Tadeu,

My apologies, this is a bug in detecting pseudo-nonlinearity. Solver does not pass solution and its derivative to custom function that defined heat source as function of temperature or its derivative. This is fixed in upcoming release. However, you can workaround this bug by specifying one of the thermal conductivities as a function. If you change:

thermalProperties(thermalmodelT,'Face',hcp_index,'ThermalConductivity',0.05);

to

thermalProperties(thermalmodelT,'Face',hcp_index,'ThermalConductivity',@(location,state) 0.05.*ones(size(location.x));

You should not get this error. However, please be aware our solver handles condiction dominant heat transfer. You are solving advection heat transfer, there might be some numerical difficulties even after you get past this bug.

Regards,

Ravi

[Tadeu Fagundes]

Thanks for the swift answer! This change solves my problem! I will keep in mind numerical problems that this might cause.