Solving Heat equation using pdepe does not give credible results

2 vues (au cours des 30 derniers jours)
Hugo Tessonneau
Hugo Tessonneau le 2 Nov 2022
Commenté : Torsten le 3 Nov 2022
Hello,
I am using pdepe to solve an heat equation, (the end goal being to simulate a fiber splicer and get the power input necessary). To create the interaction fiber to air (air heated by the graphite filament) I considered the following equation for my fiber (m=1, cylinder) :
function[c,f,s]=pdefun(x,t,u,DuDx)
ro=123.218;
C=281;
lambda=1.7;
alpha=ro*C/lambda;
c=alpha;
f=DuDx;
s=0;
end
And the boundary conditions to get heat from air to the fiber are :
function[pl,ql,pr,qr]=bcfun(xl,ul,xr,ur,t)
h=25;
R=2;
k=5;
Ts=300;%%fiber
Tf=350;%%air
pl=0;
ql=1;
pr=(h/k)*(Ts-Tf);
qr=1;
end
The initial conditions of the fiber being
function u0=icfun(x)
u0=300;
end
My problem is that when used at my fiber dimensions (270µm of diameter), I get a fiber temperature rising above the air temperature. I guess I forgot a term leading to stabilization of my sytem?
m=1;
xmesh=linspace(0,270E-6,270);
tspan=linspace(0,2,50);
sol=pdepe(m,@pdefun,@icfun,@bcfun,xmesh,tspan);
u=sol(:,:,1);
figure();
surf(xmesh,tspan,u);
The lack of visible variations on the x axe is only because of the size of my fiber (when I change xmesh I get the expected x variations).
I also get the same problem when I add a pvol term (s~=0 for pdepe model equation). Leading to absurdly low amount of power needed for the filament to generate high temperature (Code below).
%% m=0 because I chose to unroll it
function[c,f,s]=pdefunfilament(x,t,u,DuDx)%%equation for the filament
global PowerFil
ro=2.3E3;
C=720;
lambda=120;
alpha=ro*C/lambda;
c=alpha;
f=DuDx;
s=PowerFil/9E-9;%%(R*I^2)/dtau (dtau volume)
end
function u0=icfunfilament(x)%%initial conditions
u0=300;
end
function[pl,ql,pr,qr]=bcfunfilament(xl,ul,xr,ur,t)%%boundary conditions
pl=0;
ql=1;
pr=0;
qr=1;
end
For this kind of code I get results like this :
Does anyone has an idea of where my pdepe model did go wrong?
  2 commentaires
Walter Roberson
Walter Roberson le 2 Nov 2022
Modifié(e) : Walter Roberson le 2 Nov 2022
Pure black output from surf() typically means that edges are dense enough that they have taken over the display (since they are constant screen width no matter what magnification being used.)
m=1;
xmesh=linspace(0,270E-6,270);
tspan=linspace(0,2,50);
sol=pdepe(m,@pdefun,@icfun,@bcfun,xmesh,tspan);
u=sol(:,:,1);
figure();
surf(xmesh,tspan,u, 'edgecolor', 'none');
colorbar()
min(u(:)), max(u(:))
ans = 300
ans = 481.8633
function[c,f,s]=pdefun(x,t,u,DuDx)
ro=123.218;
C=281;
lambda=1.7;
alpha=ro*C/lambda;
c=alpha;
f=DuDx;
s=0;
end
function[pl,ql,pr,qr]=bcfun(xl,ul,xr,ur,t)
h=25;
R=2;
k=5;
Ts=300;%%fiber
Tf=350;%%air
pl=0;
ql=1;
pr=(h/k)*(Ts-Tf);
qr=1;
end
function u0=icfun(x)
u0=300;
end
Hugo Tessonneau
Hugo Tessonneau le 2 Nov 2022
Modifié(e) : Hugo Tessonneau le 2 Nov 2022
Thank you for your answer,
I get the part for the black surf, what I don't get is why I obtain temperature values above the temperature of the air which is heating the fiber by conducto-convection. That is impossible when considering thermodynamic.
Even if the time span I gave my algorithm is too wide in front of the system size, the temperatures should reach stable values after a certain time. Here they should never go beyond 350 K (because there is no other flux involved).

Connectez-vous pour commenter.

Réponse acceptée

Torsten
Torsten le 2 Nov 2022
Modifié(e) : Torsten le 2 Nov 2022
m=1;
xmesh=linspace(0,270E-6,270);
tspan=linspace(0,2,50);
sol=pdepe(m,@pdefun,@icfun,@bcfun,xmesh,tspan);
u=sol(:,:,1);
figure();
surf(xmesh,tspan,u, 'edgecolor', 'none');
function[c,f,s]=pdefun(x,t,u,DuDx)
rho = 123.218;
cp = 281;
lambda = 1.7;
c = rho*cp;
f = lambda*DuDx;
s = 0;
end
function [pl,ql,pr,qr] = bcfun(xl,ul,xr,ur,t)
alpha = 5;
Tf = 350;%%air
pl = 0;
ql = 1;
pr = alpha*(ur-Tf);
qr = 1;
end
function u0=icfun(x)
u0 = 300;
end
  4 commentaires
Hugo Tessonneau
Hugo Tessonneau le 3 Nov 2022
Yes sorry about that k==lambda in my case, I picked up formulas from different books and forgot to change the name of my variables.
Torsten
Torsten le 3 Nov 2022
But you defined k = 5 and lambda = 1.7 ...
But you know now about the boundary condition setup.

Connectez-vous pour commenter.

Plus de réponses (0)

Catégories

En savoir plus sur Thermal Analysis dans Help Center et File Exchange

Tags

Produits


Version

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by