1D melting problem with Neumann's analytical solution

40 views (last 30 days)
Problem: I am trying to model 1D heat transfer to analyze the melting process of a PCM (phase change material). Before solving the Stephan problem with initial and boundary conditions, I would like to implement the attached analytic solution code to calculate the volume of the PCM with the attached formula. I understand the derivation and math, but I have some problems with the implementation in MATLAB (the attached code is just a definition of the parameters needed). If someone expert in MATLAB could help me, I'd be very grateful.
Many thanks
  2 Comments
Image Analyst
Image Analyst on 11 Jan 2023
No one is going to type in all your code from an image. And we can't run an image.
If you have any more questions, then attach your data and code (m-file) to read it in with the paperclip icon after you read this:
Elisa Revello
Elisa Revello on 12 Jan 2023
Problem: model 1D heat transfer to analyze the melting process of a PCM (phase change material) - Analytic solution for Stephan problem.
How can I solve f (which has the attached expression) with Newton method?
Thanks in advance for the help.
eps0 = 0.1; % first approssimation for epsilon root
f = @(eps) (St_l/(exp(eps^2)*erf(eps)))-(St_s/((nu*exp(nu^2*eps^2))))-(eps*sqrt(pi));

Sign in to comment.

Accepted Answer

Alan Stevens
Alan Stevens on 12 Jan 2023
Edited: Alan Stevens on 12 Jan 2023
  1. eps is a built-in constant so better not to use it as a variable.
  2. Your png files show it is xi not epsilon.
  3. Try the following
  4. Don't forget erfc in the definition of f
% input data for the PCM selected
rho = 1370; % density [kg/m^3]
k_s = 0.830; % thermal conductivity solid phase [W/(m K)]
k_l = 0.660; % thermal conductivity liquid phase[W/(m K)]
cp_s = 1.69; % specific heat capacity solid phase [kJ/(kg K)]
cp_l = 1.96; % specific heat capacity liquid phase [kJ/(kg K)]
L = 227; % latent heat of the phase change [kJ/kg]
T_melt = 115+273; % melting temperature of the PCM [K]
alpha_s = k_s/(rho*cp_s); % thermal diffusivity solid phase [m^2/s]
alpha_l = k_l/(rho*cp_l); % thermal diffusivity liquid phase [m^2/s]
R_pcm = 0.5; % PCM radius [mm]
T_0 = -10+273; % initial temperature of the PCM [K]
% input data for the IMD (geometrical properties)
D_imd = 254; % diamiter [mm]
h_imd = 186; % height [mm]
T_wall = 135+273; % Temperature of the stator external surface [K]
% Stephan number definition
St_s = cp_s*(T_melt-T_0)/L; % Stephan number solid phase
St_l = cp_s*(T_wall-T_melt)/L; % Stephan number liquid phase
% Melting time definition
t_star = 0.11+(0.25/St_l); % adimensional time t* liquid phase [-]
t_melt = t_star*(R_pcm^2/alpha_l); % melting time [s]
t = linspace(0,10000); % time interval [s]
nu = sqrt(alpha_l/alpha_s); % parameter
f = @(xi) (St_l./(exp(xi.^2).*erf(xi)))-(St_s./((nu*exp(nu^2.*xi.^2).*erfc(nu*xi))))-(xi*sqrt(pi));
% Approximate numerical derivative (as I don't have the symbolic toolbox)
h = 1E-10;
fd = @(xi) (f(xi+h)-f(xi))/h;
% Newton-Raphson method
tol = 1e-8; err = 1; steps = 0;
xi = 0.1;
while err>tol
xiold = xi;
xi = xi - f(xi)/fd(xi);
err = abs(xi-xiold);
steps = steps+1;
end
disp(xi)
0.0943
disp(steps)
4
  3 Comments
Alan Stevens
Alan Stevens on 17 Jan 2023
Does this help?
%% 1) PCM Stephan problem: Analytical solution
% input data for the PCM selected (Urea-KCl)
rho = 1370; % density [kg/m^3]
k_s = 0.830; % thermal conductivity solid phase [W/(m K)]
k_l = 0.660; % thermal conductivity liquid phase[W/(m K)]
cp_s = 1.69; % specific heat capacity solid phase [kJ/(kg K)]
cp_l = 1.96; % specific heat capacity liquid phase [kJ/(kg K)]
L = 227; % latent heat of the phase change [kJ/kg]
T_melt = 115+273; % melting temperature of the PCM [K]
alpha_s = k_s/(rho*cp_s); % thermal diffusivity solid phase [m^2/s]
alpha_l = k_l/(rho*cp_l); % thermal diffusivity liquid phase [m^2/s]
R_pcm = 10; % PCM radius [mm] (supposed)
T_0 = -10+273; % initial temperature of the PCM [K] (assumed from Mansouri's paper)
% input data for the IMD (geometrical properties)
D_imd = 254; % diamiter [mm]
h_imd = 186; % height [mm]
T_wall = 135+273; % Temperature of the stator external surface [K]
% Stephan number definition
St_s = cp_s*(T_melt-T_0)/L; % Stephan number solid phase
St_l = cp_s*(T_wall-T_melt)/L; % Stephan number liquid phase
% Melting time definition
t_star = 0.11+(0.25/St_l); % adimensional time t* liquid phase [-]
t_melt = t_star*(R_pcm^2/alpha_l); % melting time [s]
t = linspace(0,10000,100); % time interval [s]
nu = sqrt(alpha_l/alpha_s); % parameter
x = linspace(0,300); % space interval
% Trascendental equation for xi
%% Newton-Raphson method to find xi
a0 = -0.00401;
a1 = 1.18669;
a2 = -0.14559;
a3 = -0.33443;
a4 = 0.16069;
a5 = -0.02155;
erf = @(xi) a0+a1*xi+a2*xi.^2+a3*xi.^3+a4*xi.^4+a5*xi.^5; % polynomial error function
erfc = @(xi) 1-(a0+a1*xi+a2*xi.^2+a3*xi.^3+a4*xi.^4+a5*xi.^5); % complementary error function
f = @(xi) (St_l./(exp(xi.^2).*erf(xi)))-(St_s./((nu*exp(nu^2.*xi.^2).*erfc(nu*xi))))-(xi*sqrt(pi));
% Approximate numerical derivative
h = 1E-10;
fd = @(xi) (f(xi+h)-f(xi))/h; % first derivative of f
tol = 1e-8; err = 1; n = 0;
xi = 0.01;
while err>tol
xiold = xi;
xi = xi - f(xi)/fd(xi);
err = abs(xi-xiold);
n = n+1;
end
fprintf('Solution value xi:')
Solution value xi:
disp(xi); % root of the implicit equation
0.0941
% Position of the liquid-solid interface
X = @(t) 2*xi*sqrt(alpha_l*t);
X_melt = X(t_melt);
% PCM volume formula
V = pi*((1/2*D_imd+X_melt)^2-(1/2*D_imd)^2)*h_imd; % volume [mm^3]
m_pcm = (rho*V)*(10^-9); % PCM mass [kg]
fprintf('Melting time calculated [s]:')
Melting time calculated [s]:
disp(t_melt);
7.2785e+05
fprintf('Melting position calculated [mm]:')
Melting position calculated [mm]:
disp(X_melt);
2.5174
fprintf('PCM mass calculated [kg]:')
PCM mass calculated [kg]:
disp(m_pcm);
0.5170
% T_l(x,t) and T_s(x,t) equations definition
T_l = @(x,t) T_wall+(T_melt-T_wall).*(erf(x./2.*sqrt(alpha_l.*t))./erf(xi));
T_s = @(x,t) T_0+(T_melt-T_0)*erfc((x./2.*sqrt(alpha_l.*t))./erfc(nu*xi));
% Plot Temperature distribution for T_l(x,t)
x = 0:100:3000; % define range and mesh of x and y which will be shown in figure
%y = 0:100:1e4; % define t interval (should be from 0 to 10000 [s])
[X, tm] = meshgrid(x, t);
surf(X, tm, T_l(X,tm));
figure;
contourf(X, tm, T_l(X,tm),20);
Elisa Revello
Elisa Revello on 5 Feb 2023
I am trying to solve the Analytical solution for the Transient heat transfer problem of a Phase Change Material and I want to evaluate the temperature profile through the PCM. I have some troubles related to the definition of the equations attached ("analytical equations for T(x,t)".png) in a for cycle.
In the following lines of MATLAB code I've just tryed to implement the temperature distribution, but I don't understand how obtain the correct overall T distribution along the spatial and temporal coordinates.
Many thanks in advance for the help.
%% Transient 1D Heat Conduction - PCM
% domain discretization
X0 = 0; % initial position
XL = R_pcm; % final position [m]
N = 100; % number of nodes
dx = XL/N; % space discretization
t = 1000; % final time [s]
dt = 1; % time discretization
M = t/dt; % number of time steps needed for the calculation
x = X0:dx:XL;
for j=1:M
T_l(1,j)=T_wall; % boundary condition for liquid phase T_l(0,t)=T_wall
for i=2:N
T_l(i,j) = T_wall+(T_melt-T_wall)*(erf(i./2.*sqrt(alpha_l.*j))./erf(xi));
if x(i)==X_melt
T_l(i,j) = T_melt; % boundary condition for liquid phase T_l(X_melt,t)=T_melt
end
end
end
T_s(i,1) = T_0; % initial condition for solid phase T_s(x,0)=T_0
for j=2:M
for i=1:N
T_s(i,j) = T_0+(T_melt-T_0)*erfc((i./2.*sqrt(alpha_s.*j))./erfc(nu*xi));
if x(i)==X_melt
T_s(i,j) = T_melt; % boundary condition for solid phase T_s(X_melt,t)=T_melt
elseif x(i)==XL
T_s(i,j) = T_0;
end
end
end
plot(t,T_l);
title('Temperature liquid phase vs Time');

Sign in to comment.

More Answers (0)

Categories

Find more on Quantum Chemistry in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!