# 1D melting problem with Neumann's analytical solution

27 views (last 30 days)
Elisa Revello on 11 Jan 2023
Commented: Elisa Revello on 5 Feb 2023
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 CommentsShowHide 1 older comment
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));

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
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');