MATLAB: Solving 1D heat equation with constant heat flux (boundary condition)

boundary conditionheat equationheat fluxpdepe

Hello everyone, i am trying to solve the 1-dimensional heat equation under the boundary condition of a constant heat flux (unequal zero). The heat flux is on the left and on the right bound and is representing the heat input into the material through convective heat transfer. I am using pdepe to solve the heat equation and with dirichlet boundary conditions it is working. To represent the heat flux i was thinking of a Neumann boundary condition, but i can't figure out how to contribute the value of the heat flux into the boundary condition. For a heat flux of zero (isolation) the code would look like this (leaving out the declaration of the parameters):

function u0=IC(x)
u0=T_0;
end
function [pl,ql,pr,qr]=BC(xl,ul,xr,ur,t)
ql=1;
pl=0;
qr=1;
pr=0;
end
function [c,f,s]=DGL(x,t,u,DuDx)
c=1;
f=a*DuDx;
s=0;
end
% solution
x=linspace(0,L,10);
t=linspace(0,t_End,10);
sol=pdepe(m,@PDE,@IC,@BC,x,t);
u=sol(:,:,1);
uu=sol(end,:,1);
end

Now the heat flux being unequal zero, my idea was to modify pl and pr in dependence of the heat flux density, but i can't find a connection. Has anyone an idea?

Regards, Malte

Best Answer

What is the boundary condition you are trying to implement ?

lambda*dT/dn = alpha*(T-T_ext) (n: outer normal) ?

Best wishes

Torsten.

Related Solutions

MATLAB: Is the boundary condition for the heat transfer problem not satisfied using the PDEPE function in MATLAB

We have verified that there is a limitation in MATLAB 7.0 (R14)+ in the way that the PDEPE function handles boundary conditions with certain values.

PDEPE imposes Boundary Conditions in terms of flux "f()":

    p() + q()*f() = 0

At the left bound, "x=0". If the flux term is in a form such that it is zero for any specified boundary condition, the boundary condition will be ignored.

For example, "f() = (A+B*x)*du/dx" where "A" and "B" are constants. If "A = 0", flux is "(0+B*x)*du/dx".

You specify the boundary condition on the left bound to be "du/dx =0".

On the left bound, "B*(0)*du/dx = 0". This boundary condition is always satisfied regardless the value of "du/dx" and MATLAB will disregard the previously specified boundary condition.

Hence this particular combination:

     f() = x*du/dx
     0 <= x 
     Boundary Condition at 0,  du/dx = 0

gives an incorrect boundary result.

To work around this issue, shift the value of "x" from 0 to a very small number called delta that is close to 0. Define another point next to delta as follow:

x = [delta, 1.01*delta, linspace(2*delta,1000,100)];

The boundary condition on the left side should also be changed from:

pl = 0; 
ql = 1;
pr = 0; 
qr = 1;

to

pl = 0; 
ql = 1/(A+Bx*xl);
pr = 0; 
qr = 1;

MATLAB: Solving second order PDE

The reason that pdepe imposes a boundary condition of the flux equal zero at the

center is that this is required for the problem to be mathematically well-posed.

Imposing a prescribed temperature at the center would require that the flux go to

infinity there.

An easy way to understand this is to solve the problem with the left end a small distance

from the center and with a fine mesh. I have attached a short script below that shows this.

function matlabAnswers_6_27_2020
r0=1e-6;
x = linspace(r0,1,1000);
tf=1;
t = linspace(0,tf,40);
pdeFunc = @(x,t,u,DuDx) heatpde(x,t,u,DuDx);
icFunc = @(x) heatic(x);
bcFunc = @(xl,ul,xr,ur,t) heatbcDirichlet(xl,ul,xr,ur,t);
m=1;
sol = pdepe(m, pdeFunc,icFunc,bcFunc,x,t);
figure; plot(t, sol(:,end)); grid on; title 'Temperature at outer surface'
figure; plot(t, sol(:,1)); grid on; title 'Temperature at center'
figure; plot(x, sol(end,:)); grid; title 'Temperature at final time'
end
function [c,f,s] = heatpde(x,t,u,DuDx)
c = 1;
f = DuDx;
s = 0;
end
function u0 = heatic(x)
u0 = 0;
end
function [pl,ql,pr,qr] = heatbcDirichlet(xl,ul,xr,ur,t)
pl = ul-1;
ql = 0;
pr = 0;
qr = 1;
end

Related Question