I use the pdepe.m command to solve the heat equation.
Everything seems to work, and I do get results.
However when I calculate the heat that entered the domain through the boundary, it is not equal to the energy that was required to heat the particle.

In other words, the particle temperature increases faster than possible, since a smaller amount of heat is fed into the particle through the boundary of the particle.

At the center of the domain (it is a sphere) the neumann boundary is dT/dx = 0
At the edge a heat flux is prescribed as:
keff*(dT/dx) = hconv*(T_edge - T_sur)

I have no idea what goes wrong, but I think it is an error in the boundary conditions.


The M-file looks like:
%model heat



function pdex

clear all
close all
clc

global r0 rho_0 T_0 T_sur_0  Nx Nt t_end keff 



r0= 0.0005;         % radius [m]                     
t_end= 0.2;         % time [sec]                     
Nt=100;             % Number of timesteps
Nx=100;             % Number of spacesteps
rho_0=200;          % initial char density [kg/m^3]
T_0= 1073;          % initial particle temperature [K]
T_sur_0=1273;       % temperature of the surroundings [K]

m = 0;                                  % A parameter corresponding to the symmetry of the problem
x = linspace(0,r0,Nx)    ;              % the first two numbers specify the range, the third number gives the number of intervals
t = linspace(0,t_end,Nt)    ; 

sol = pdepe(m,@pdexpde,@pdexic,@pdexbc,x,t);


u = sol;                % temperature

%% A surface plot of the temperature in space and time
figure(1);
surf(x,t,u);
title('Temperature as function of space and time');
xlabel('Distance [m]');
ylabel('Time [sec]');
zlabel('Temperature [K]')



function [c,f,s] = pdexpde(x,t,u,DuDx)

global rho_0 keff

keff = 0.1;                     % effective thermal conductivity of char. [W/(m*K)] 
Cp_char =1300;                  % heat capacity of char, assumed constant [J/(kg*K)]
rho_char_init =rho_0;              % inital char density [kg/m^3]
                             
%% system of equations:
c = [1*rho_char_init*Cp_char];
f = [(keff*(DuDx))];
s = [0];
                                                                                                                                                                                                                                                                                                        
% --------------------------------------------------------------------------
function u0 = pdexic(x)

global T_0 
u0 = [T_0];            % initial conditions of char density, oxygen mass fraction and temperature.                       


% --------------------------------------------------------------------------
function [pl,ql,pr,qr] = pdexbc(xl,ul,xr,ur,t)  

global T_sur_0

%boundary conditions
pl = [0];
ql = [1];
pr = [(1/(keff))*(hc*(ur-T_sur_0))];
qr = [1];

On 24 Mrz., 12:13, Niek <u...@compgroups.net/> wrote:
> I use the pdepe.m command to solve the heat equation.
> Everything seems to work, and I do get results.
> However when I calculate the heat that entered the domain through the bou=
ndary, it is not equal to the energy that was required to heat the particle=
..
>
> In other words, the particle temperature increases faster than possible, =
since a smaller amount of heat is fed into the particle through the boundar=
y of the particle.
>
> At the center of the domain (it is a sphere) the neumann boundary is dT/d=
x =3D 0
> At the edge a heat flux is prescribed as:
> keff*(dT/dx) =3D hconv*(T_edge - T_sur)
>
> I have no idea what goes wrong, but I think it is an error in the boundar=
y conditions.
>
> The M-file looks like:
> %model heat
>
> function pdex
>
> clear all
> close all
> clc
>
> global r0 rho_0 T_0 T_sur_0 =A0Nx Nt t_end keff
>
> r0=3D 0.0005; =A0 =A0 =A0 =A0 % radius [m] =A0 =A0 =A0 =A0 =A0 =A0 =A0 =
=A0 =A0 =A0
> t_end=3D 0.2; =A0 =A0 =A0 =A0 % time [sec] =A0 =A0 =A0 =A0 =A0 =A0 =A0 =
=A0 =A0 =A0
> Nt=3D100; =A0 =A0 =A0 =A0 =A0 =A0 % Number of timesteps
> Nx=3D100; =A0 =A0 =A0 =A0 =A0 =A0 % Number of spacesteps
> rho_0=3D200; =A0 =A0 =A0 =A0 =A0% initial char density [kg/m^3]
> T_0=3D 1073; =A0 =A0 =A0 =A0 =A0% initial particle temperature [K]
> T_sur_0=3D1273; =A0 =A0 =A0 % temperature of the surroundings [K]
>
> m =3D 0; =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =
=A0% A parameter corresponding to the symmetry of the problem

m =3D 2; (because you work in a spherical coordinate system)

> x =3D linspace(0,r0,Nx) =A0 =A0; =A0 =A0 =A0 =A0 =A0 =A0 =A0% the first t=
wo numbers specify the range, the third number gives the number of interval=
s
> t =3D linspace(0,t_end,Nt) =A0 =A0;
>
> sol =3D pdepe(m,@pdexpde,@pdexic,@pdexbc,x,t);
>
> u =3D sol; =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0% temperature
>
> %% A surface plot of the temperature in space and time
> figure(1);
> surf(x,t,u);
> title('Temperature as function of space and time');
> xlabel('Distance [m]');
> ylabel('Time [sec]');
> zlabel('Temperature [K]')
>
> function [c,f,s] =3D pdexpde(x,t,u,DuDx)
>
> global rho_0 keff
>
> keff =3D 0.1; =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 % effective thermal=
 conductivity of char. [W/(m*K)]
> Cp_char =3D1300; =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0% heat capacity of ch=
ar, assumed constant [J/(kg*K)]
> rho_char_init =3Drho_0; =A0 =A0 =A0 =A0 =A0 =A0 =A0% inital char density =
[kg/m^3]
>
> %% system of equations:
> c =3D [1*rho_char_init*Cp_char];
> f =3D [(keff*(DuDx))];
> s =3D [0];
>
> % -----------------------------------------------------------------------=
---
> function u0 =3D pdexic(x)
>
> global T_0
> u0 =3D [T_0]; =A0 =A0 =A0 =A0 =A0 =A0% initial conditions of char density=
, oxygen mass fraction and temperature. =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0=
 =A0 =A0
>
> % -----------------------------------------------------------------------=
---
> function [pl,ql,pr,qr] =3D pdexbc(xl,ul,xr,ur,t) =A0
>
> global T_sur_0
>
> %boundary conditions
> pl =3D [0];
> ql =3D [1];
> pr =3D [(1/(keff))*(hc*(ur-T_sur_0))];
> qr =3D [1];

pr =3D hc*(ur-T_sur_0);
qr =3D 1;

And don't forget to make hc available in this m-file
(it's not defined in your code above).

Best wishes
Torsten.

Hello Torsten,

Thanks for your answer! When removing the factor (1/keff) the results are correct. The energy that goes through the boundary into the particle due to convection is equal to the increase in internal energy.

But in that case the dimensions of my boundary conditions are not correct? How could you explain that?

at the edge:
dT/dx  = hconv * dT
In SI units:
[K/m] =# (W/(m^2 * K)) * K
so the dimensions are not correct, and I would say that the conductivity should be introduced on the left hand side which has dimensions [W/(m*K)], so in that case we get as dimensions [W/m^2]  However this gives errors in the result.

Do you have any idea on that?

Thanks

On 24 Mrz., 15:07, Niek <u...@compgroups.net/> wrote:
> Hello Torsten,
>
> Thanks for your answer! When removing the factor (1/keff) the results are=
 correct. The energy that goes through the boundary into the particle due t=
o convection is equal to the increase in internal energy.
>
> But in that case the dimensions of my boundary conditions are not correct=
? How could you explain that?
>
> at the edge:
> dT/dx =A0=3D hconv * dT
> In SI units:
> [K/m] =3D# (W/(m^2 * K)) * K
> so the dimensions are not correct, and I would say that the conductivity =
should be introduced on the left hand side which has dimensions [W/(m*K)], =
so in that case we get as dimensions [W/m^2] =A0However this gives errors i=
n the result.
>
> Do you have any idea on that?
>
> Thanks

The boundary condition in pdepe is of the form
p+f*q =3D 0.
Since you set f =3D k_eff*dT/dr in the routine where you define
the pde, setting q=3D1 and p=3Dhconv*(ur-T_u) gives
-k_eff*dT/dr =3D hconv*(T-T_u)
as boundary condition at r=3DR - as required.

Best wishes
Torsten.

I think that is the answer I was looking for,

I also have diffusion equations, so I go check my boundaries therea also!

Thanks!