PDETool temperature gradient-dependent internal heat generation

Question

Tadeu Fagundes 2020-2-12

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/505133-pdetool-temperature-gradient-dependent-internal-heat-generation

评论： Tadeu Fagundes 2020-2-12

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Ravi Kumar 2020-2-12

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/505133-pdetool-temperature-gradient-dependent-internal-heat-generation#answer_415275

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