1D melting problem with Neumann's analytical solution

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

[Image Analyst]

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:

TUTORIAL: How to ask a question (on Answers) and get a fast answer

[Elisa Revello]

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

Answer 1

[Alan Stevens]

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]

I am trying to solve the Analytical solution for the Heat conduction phase change problem and I want to find the Temperature profiles both for liquid and solid phases, but I have some troubles related to the definition of the equations attached ("analytical equations for T(x,t)".png), which are function of two variables (x,t).

Particularly, I want to plot the T_l(x,t) and T_s(x,t) distributions as a function of x and I don't understand how to solve this last part of the attached code:

% 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:0.1:300;  % define x interval [mm]
y = 0:0.1:3600; % define t interval [s]
[X, Y] = meshgrid(x, y); % define mesh of x and y which will be shown in figure
surf(X, Y, T_l(X,Y));
figure;
contourf(X, Y, T_l(X,Y));

Thanks in advance for the help.

[Alan Stevens]

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]

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