Solving an integral within a partial derivative

Hi everyone,

I am having trouble solving the following heat equation in Matlab

As you can see there will be an integral within the partial derivative.

When solving my equations work, but do not work over time. Thus i have a feeling I do something wrong with solving the integral. This is how I currently describe the problem. I did not fully iplement the higher order terms, as it does not work yet.

P = dimensionless term cf/kf

Z = 1/kf(w/W)

tmax = 1000 s

Any help would be greatly appreciated

m = 0;
sol = pdepe(m,@heatpde,@heatic,@heatbc,x,t);


function [c,f,s] = heatpde(x,t,u,dudx)
global B tmax Z P
c =P;
f = dudx;

y = @(u) (u.*(1-B/3*u+7*B^2/45*u.^2));
A = sqrt(2*B *integral(y,0,tmax))
s = -Z*u/A;
end

function u0 = heatic(x)

u0 = 0;
end

function [pl,ql,pr,qr] = heatbc(xl,ul,xr,ur,t)
pl = ul-1;
ql = 0;
pr = ur;
qr = 0;

Answer 1

[Torsten]

You have to solve a second ordinary differential equation for y:

dy/dt = 2*Ste * theta_f*(1-Ste/3 * theta_f + 7/45*Ste^2*theta_f^2)
y(t=0,x) = 0

Then y_LS = sqrt(y) in your equation for theta_f.

This can give problems because "pdepe" is designed to solve partial differential equations, not ordinary differential equations.

But since you prescribe theta_f at both ends of the spatial interval, you can compute y at these endpoints, too, and prescribe its value in "heatbc":

at x = 0: y = 2*Ste*integral(1-Ste/3+7/45*Ste^2) dt = 2*Ste*(1-Ste/3+7/45*Ste^2)*t

at x = L: y = 2*Ste*integral(0 dt) = 0

Is y_LS' differentiation with respect to t or with respect to x ? Same question for theta_f'.

[kjell revenberg]

Hi Torsten,

Thanksyou for the ultiple updates you have given. I am currently trying to solve it so I did not answer yet.

This is the overview, where theta_f is the dimensionless temperature of the fin. The fin Temperature changes over time and space as it goes further away from the center.

The current boundary conditions that can also be seen in my code are:

With theta_f being = (T-Tmelt)/(Tfin-Tmelt),

thus left side starts at T= Tfin

Rest of te fin starts at T = Tmelt

right side is adiabatic.

The probem is thus that the fin temperature is dependend on Yls, but also the other way around. Both thus depending on time and space.

This is the paper if youre interested: Theoretical modeling and optimization of fin-based enhancement of heat transfer into a phase change material - ScienceDirect

-----

I understand how you prescribed the problem, I only do not know how i can implement it in de boundary condition and then bring it back to the main function

[Torsten]

Since you set theta_f = 0 in your boundary condition part for "pdepe" instead of d(theta_f)/dx = 0 as in your description from above, I suggested to set

at x = 0: y = 2*Ste*integral(1-Ste/3+7/45*Ste^2) dt = 2*Ste*(1-Ste/3+7/45*Ste^2)*t

at x = L: y = 2*Ste*integral(0 dt) = 0

in the boundary condition part for y.

You can try y = 2*Ste*(1-Ste/3+7/45*Ste^2)*t at x = 0 and dy/dx = 0 at x = L instead :

function [pl,ql,pr,qr] = heatbc(xl,ul,xr,ur,t)
  Ste = ...;
  pl = [ul(1)-1;ul(2)-2*Ste*(1-Ste/3+7/45*Ste^2)*t];
  ql = [0;0];
  pr = [0;0];
  qr = [1;1];
end

[kjell revenberg]

Hi Torsten,

Thanks for the troubleshooting.

I agree that I did not implement the boundary condition right and understand why you want to put them in the boudnary condtions. However I still face the issue that this is now not correcly indexed. (Index exceeds the number of array elements)

[Torsten]

m = 0;
sol = pdepe(m,@heatpde,@heatic,@heatbc,x,t);
function [c,f,s] = heatpde(x,t,u,dudx)
  global B P
  theta_f = u(1);
  y_LS = sqrt(u(2));
  d_theta_f_dx = dudx(1);
  d_y_LS_dx = 1/(2*u(2))*dudx(2);
  c = [P;1];
  f = [dudx(1);0];
  s = [-1/(k_f*w/W)*(theta_f/y_LS+B*theta_f*(theta_f+2*d_theta_f_dx/d_y_LS_dx*y_LS)/(6*y_LS)-...
       B^2*theta_f^2*(40*(d_theta_f_dx/d_y_LS_dx)^2*theta_f^2+85*theta_f*d_theta_f_dx/d_y_LS_dx*y_LS+19*theta_f^2)/(360*y_LS));...
       2*B * theta_f*(1-B/3 * theta_f + 7/45*B^2*theta_f^2)]
end
function u0 = heatic(x)
 u0 = [0;0];
end
function [pl,ql,pr,qr] = heatbc(xl,ul,xr,ur,t)
  global B
  pl = [ul(1)-1;ul(2)-2*B*(1-B/3+7/45*B^2)*t];
  ql = [0;0];
  pr = [0;0];
  qr = [1;1];
end

[kjell revenberg]

Hi Torsten,

Sorry for the late response, it was due to my holiday.

I cleaned up the code, so it is more readable. This seems to be working, however the results are not there yet as the integration tolerances are not met. This is due to the fact that at point t =0 theta_f and Y_ls are both zero creating singularity.

This can be solved by adding the value for theta_f for t=0:t*

Where t* is significantly smaller then tmax.

I am just not sure how I can add this. Would you use a forloop?

Thanks in advance!

clc;close all;clear all;


%% Temperature input

TB = 64;                                % Fin base temperature;
TM = 44;                                % Temperature PCM



%% Fin input

k_f = 237;                              % Thermal conductivity fin
Cp_f = 2300;                            % Heat capacity 
rho_f = 2710;                           % Density

%% PCM properties

rho_p = 780;                            % Density
Cp_p = 2300;                            % Heat capacity
k_p = 0.15;                             % Conductivity
H = 244000;                             % Latent Heat
a_p = k_p/(rho_p*Cp_p);                 % Thermal diffusivity

%% Meshing input

w = 0.005;                              % 1/2 the fin thickness
W = 0.05;                               % 1/2 the pcm height

tmax = 1000;                            % Maximum timestep
L = 0.2;                                % Maximum length
x = 0:L/10:L;
t = 0:tmax/10:tmax;

%% Dimensionless numbers

STE = Cp_p*(TB-TM)/H;                  % Stefan number
Cp_d = (rho_f*Cp_f)/(rho_p*Cp_p);      % Dimensionless heat capacity
k_d = k_f/k_p;                         % Dimensionless conductivity
w_d = w/W;                             % Aspect ratio

t_d = a_p.*t/W^2;                      % Dimensionless time
x_d = x/W;                             % Dimensionless distance
tdmin = tmax/10000;

 %% Storing values in one place

C.STE =STE;
C.k_d = k_d;
C.w = w;
C.W = W;
C.Cp_d = Cp_d;



%% Solving the function

m = 0;
heat = @(x,t,u,dudx) heatpde(x,t,u,dudx,C);
bc = @(xl,ul,xr,ur,t) heatbc(xl,ul,xr,ur,t,C);
sol = pdepe(m,heat,@heatic,bc,x,t);

function [c,f,s] = heatpde(x_d,t_d,u,dudx,C)

k_d = C.k_d;
STE = C.STE;
w = C.w;
W = C.W;
Cp_d = C.Cp_d;

theta_f = u(1);
y_LS = sqrt(u(2));
dtheta_fdx = dudx(1);
d_y_LS_dx = 1/(2*u(2))*dudx(2);
c = [Cp_d/k_d;1];
f =[dudx(1);0];

s = [-1/(k_d*w/W)*(theta_f/y_LS+STE*theta_f*(theta_f+2*dtheta_fdx/d_y_LS_dx*y_LS)/(6*y_LS)-...
       STE^2*theta_f^2*(40*(dtheta_fdx/d_y_LS_dx)^2*theta_f^2+85*theta_f*dtheta_fdx/d_y_LS_dx*y_LS+19*theta_f^2)/(360*y_LS));...
       2*STE * theta_f*(1-STE/3 * theta_f + 7/45*STE^2*theta_f^2)];
end

function u0 = heatic(x_d)
 u0 = [0;0];
end
function [pl,ql,pr,qr] = heatbc(xl,ul,xr,ur,t_d,C)
 STE = C.STE;
 pl = [ul(1)-1;ul(2)-2*STE*(1-STE/3+7/45*STE^2)*t_d];
 ql = [0;0];
 pr = [0;0];
 qr = [1;1];
end

[Torsten]

Make a function in which you compute theta_f according to your infinite series for given t and x (e.g. by a while loop which stops if the absolute value of the n-th summand in the infinite series is smaller than a specified tolerance).

Call this function in "heatic" to give initial values to theta_f at t=t*. The initial values for y_ls can also be computed with this function by using "trapz" to evaluate the integral from 0 to t*.

[kjell revenberg]

Hi Torsten,

I calculated the function seperately as it added a lot of time to the simulation.

I used a t* of 10^-4

For theta_f it goes from 1 to 0.025 and at the third point it already zero.

Y_ls should be zero as it is neglected at this stage.

However I still get an error due to not meeting the integration tolerances. These values are similar to the paper so it should be possible, however I can't get it to work. I now changed u0(1) to a small avlue to see if i could find a point where it works, but it doesn''t

clc;close all;clear all;


%% Temperature input

TB = 64;                                % Fin base temperature;
TM = 44;                                % Temperature PCM



%% Fin input

k_f = 237;                              % Thermal conductivity fin
Cp_f = 2300;                            % Heat capacity 
rho_f = 2710;                           % Density

%% PCM properties

rho_p = 780;                            % Density
Cp_p = 2300;                            % Heat capacity
k_p = 0.15;                             % Conductivity
H = 244000;                             % Latent Heat
a_p = k_p/(rho_p*Cp_p);                 % Thermal diffusivity

%% Meshing input

w = 0.005;                              % 1/2 the fin thickness
W = 0.05;                               % 1/2 the pcm height

tmax = 5;                            % Maximum timestep
tmin = 0.1;
L = 0.2;                                % Maximum length
x = 0:L/10:L;
t = tmin:tmin/10^5:tmax;
ts = 0:tmin/10:tmin;

%% Dimensionless numbers

STE = Cp_p*(TB-TM)/H;                  % Stefan number
Cp_d = (rho_f*Cp_f)/(rho_p*Cp_p);      % Dimensionless heat capacity
k_d = k_f/k_p;                         % Dimensionless conductivity
w_d = w/W;                             % Aspect ratio

t_d = a_p.*t/W^2;                      % Dimensionless time
x_d = x/W;                             % Dimensionless distance
x_d = transpose(x_d);
ts_d =a_p.*ts/W^2;
ts_d2 =a_p.*tmin/W^2;

 %% Storing values in one place

C.STE =STE;
C.k_d = k_d;
C.w = w;
C.W = W;
C.Cp_d = Cp_d;



%% Solving the function


m = 0;
heat = @(x,t,u,dudx) heatpde(x_d,t_d,u,dudx,C);
bc = @(xl,ul,xr,ur,t_d) heatbc(xl,ul,xr,ur,t_d,C);
sol = pdepe(m,heat,@heatic,bc,x_d,t_d);

function [c,f,s] = heatpde(x_d,t_d,u,dudx,C)

k_d = C.k_d;
STE = C.STE;
w = C.w;
W = C.W;
Cp_d = C.Cp_d;

theta_f = u(1);
y_LS = sqrt(u(2));
dtheta_fdx = dudx(1);
d_y_LS_dx = 1/(2*u(2))*dudx(2);
c = [Cp_d/k_d;1];
f =[dudx(1);0];

s = [-1/(k_d*w/W)*(theta_f/y_LS+STE*theta_f*(theta_f+2*dtheta_fdx/d_y_LS_dx*y_LS)/(6*y_LS)-...
       STE^2*theta_f^2*(40*(dtheta_fdx/d_y_LS_dx)^2*y_LS^2+85*theta_f*dtheta_fdx/d_y_LS_dx*y_LS+19*theta_f^2)/(360*y_LS));...
       2*STE * theta_f*(1-STE/3 * theta_f + 7/45*STE^2*theta_f^2)];
end

function u0 = heatic(x_d)
u0 = [0.1;0];
end
function [pl,ql,pr,qr] = heatbc(xl,ul,xr,ur,t_d,C)
 STE = C.STE;
 pl = [ul(1)-1;ul(2)-2*STE*(1-STE/3+7/45*STE^2)*t_d];
 ql = [0;0];
 pr = [0;0];
 qr = [1;1];
end

[Torsten]

u0 = [0.1;0];
d_y_LS_dx = 1/(2*u(2))*dudx(2);

This will immediately result in a division by zero.

[kjell revenberg]

yes, that is right, but assuming Y_ls is almost zero. You could enter a really small value, which still gives the same error

[kjell revenberg]

Hi Torsten, I am still struggling. I found some mistakes and think that the boundary conditions for y are not correct, but should be the same as theta_f. However it is still not working. I would really appreciate it if you could look at it again:

%% Temperature input

TB = 64;                                % Fin base temperature;
TM = 44;                                % Temperature PCM



%% Fin input

k_f = 237;                              % Thermal conductivity fin
Cp_f = 2300;                            % Heat capacity 
rho_f = 2710;                           % Density

%% PCM properties

rho_p = 780;                            % Density
Cp_p = 2300;                            % Heat capacity
k_p = 0.15;                             % Conductivity
H = 244000;                             % Latent Heat
a_p = k_p/(rho_p*Cp_p);                 % Thermal diffusivity

%% Meshing input

w = 0.005;                              % 1/2 the fin thickness
W = 0.05;                               % 1/2 the pcm height

tmax = 1000;                            % Maximum timestep
tmin = 0.1;
L = 0.2;                                % Maximum length
x = 0:L/400:L;
t = tmin:tmin:tmax;
ts = 0:tmin/10:tmin;

%% Dimensionless numbers

STE = Cp_p*(TB-TM)/H;                  % Stefan number
Cp_d = (rho_f*Cp_f)/(rho_p*Cp_p);      % Dimensionless heat capacity
k_d = k_f/k_p;                         % Dimensionless conductivity
w_d = w/W;                             % Aspect ratio

t_d = a_p.*t/W^2;                      % Dimensionless time
x_d = x/W;                             % Dimensionless distance
x_d = transpose(x_d);
ts_d =a_p.*ts/W^2;
ts_d2 =a_p.*tmin/W^2;

 %% Storing values in one place

C.STE =STE;
C.k_d = k_d;
C.w = w;
C.W = W;
C.Cp_d = Cp_d;



%% Solving the function


m = 0;
heat = @(x_d,t_d,u,dudx) heatpde(x_d,t_d,u,dudx,C);
bc = @(xl,ul,xr,ur,t_d) heatbc(xl,ul,xr,ur,t_d,C);
options = odeset('RelTol', 1, 'AbsTol', 1);
sol = pdepe(m,heat,@heatic,bc,x_d,t_d,options);

function [c,f,s] = heatpde(x_d,t_d,u,dudx,C)

k_d = C.k_d;
STE = C.STE;
w = C.w;
W = C.W;
Cp_d = C.Cp_d;

theta_f = u(1);
y_LS = sqrt(u(2));
dtheta_fdx = dudx(1);
d_y_LS_dx = 1/(2*u(2))*dudx(2);
c = [Cp_d/k_d;1];
f =[dudx(1);0];

s = [-1/(k_d*w/W)*(theta_f/y_LS+STE*theta_f*(theta_f+2*dtheta_fdx/d_y_LS_dx*y_LS)/(6*y_LS)- ...
    STE^2*theta_f*(40*(dtheta_fdx/d_y_LS_dx)^2*y_LS^2+85*theta_f*dtheta_fdx/d_y_LS_dx*y_LS+19*theta_f^2)/(360*y_LS));...
       2*STE * theta_f*(1-STE/3 * theta_f + 7/45*STE^2*theta_f^2)];
end

function u0 = heatic(x_d)
u0 = [0.01;0.00001];
end

function [pl,ql,pr,qr] = heatbc(xl,ul,xr,ur,t_d,C)
 STE = C.STE;
 pl = [ul(1)-1;ul(2)-1];
 ql = [0;0];
 pr = [0;0];
 qr = [1;1];
end

[Torsten]

y(2) must always have initial value 0 for all x because t equals 0, and the integral from 0 to 0 of a function is 0. That's why I don't understand how the division by y(2) in your PDE can be dealt with. Maybe setting u0(2) = 0.00001 as you did is the best way to handle this problem.

If you set y(1) = 1 at x = 0, y(2) can by directly computed via the integral definition. It won't be y(2) = 1 - that's for sure. It will be t*(2*STE * 1*(1-STE/3 * 1 + 7/45*STE^2*1)).

[kjell revenberg]

HI Torsten,

I think I found a solution to part of the problem.

dtheta_fdx and dy_ls_dx are both 0 at some points. Thus creating singularity to the equation.

I solved this by adding a small epsilon to the problem, some of the results seem accurate while the PCM temperature is not right.

I wondered how you came to the formula for

dy_ls_dx as I think there lies the error.

Thanks in advance

[Torsten]

I wondered how you came to the formula for

dy_ls_dx as I think there lies the error.

Thanks in advance

Yes, it's sqrt(u(2)) instead of u(2) in the denominator:

y_ls = sqrt(u(2)) -> d(y_ls)/dx = 1/(2*sqrt(u(2))) * dudx(2)