PDETool temperature gradient-dependent internal heat generation

2 次查看（过去 30 天）
Tadeu Fagundes
Tadeu Fagundes 2020-2-12
评论： 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

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

采纳的回答

Ravi Kumar
Ravi Kumar 2020-2-12
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
  1 个评论
Tadeu Fagundes
Tadeu Fagundes 2020-2-12
Thanks for the swift answer! This change solves my problem! I will keep in mind numerical problems that this might cause.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Structural Mechanics 的更多信息

标签

产品


版本

R2019b

Community Treasure Hunt

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

Start Hunting!

Translated by