Solving Heat equation using pdepe does not give credible results

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?

[Walter Roberson]

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]

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).

Answer 1

[Torsten]

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

[Hugo]

Thank you, so the reason was the fact that I kept my boundary condition constant despite the evolution of the fiber temperature.

[Torsten]

Yes. And be careful with the setting of "f" in "pdefun".

The boundary condition reads

p + q*f =0.

You set f = DuDx and qr = 1, so the heat conductivity of the fiber does not explicitly appear in your boundary condition although it should. But maybe "lambda" was hidden in "h" or "k" - I don't know.

In the code above, the boundary condition at r=R is set to

-lambda*dT/dr = alpha*(T-Tf)

[Hugo]

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]

But you defined k = 5 and lambda = 1.7 ...

But you know now about the boundary condition setup.