Filter löschen
Filter löschen

PDETool temperature gradient-dependent internal heat generation

1 Ansicht (letzte 30 Tage)
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

Akzeptierte Antwort

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

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (0)

Produkte


Version

R2019b

Community Treasure Hunt

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

Start Hunting!

Translated by