Solving an integral within a partial derivative

2023 12 21
1 个回答
更新时间：2024 1 9
7 次查看（30 天）

0 个投票

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;

请先登录，再回答此问题。

回答（1 个）

Torsten
Torsten 2023-12-21
编辑：Torsten 2023-12-21

0 个投票

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

13 个评论
显示 11更早的评论 隐藏 11更早的评论

kjell revenberg
kjell revenberg 2023-12-21
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
Torsten 2023-12-21
编辑：Torsten 2023-12-21
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
kjell revenberg 2023-12-21
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
Torsten 2023-12-21
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
kjell revenberg 2024-1-3
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
Torsten 2024-1-3
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
kjell revenberg 2024-1-5
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
Torsten 2024-1-5
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
kjell revenberg 2024-1-5
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
kjell revenberg 2024-1-8
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
Torsten 2024-1-8
编辑：Torsten 2024-1-8
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
kjell revenberg 2024-1-9
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
Torsten 2024-1-9
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)

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 Heat Transfer 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by