MATLAB Newton Raphson method with array - Stephan problem

I am trying to model the analytical solution for a heat conduction problem and I have some troubles with the implementation of Newton-Raphson method to solve a for cycle with array variables.

Particularly, I obtain errors in the following loop to find xi root of a trascendental equation, which appears in the definition of the solid-liquid interface (X_melt in the attached code).

Thanks in advance, if someone could help me.

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.05; % PCM length [m] (supposed)
T0 = -10+273; % initial temperature of the PCM [K] (assumed from Mansouri's paper)
Tamb = 20+273; % ambient temperature [K]

% input data for the IMD (geometrical properties)
D_imd = 254*(10^-3); % diamiter [m]
h_imd = 186*(10^-3); % height [m]

% Wall temperature definition (square wave)
t = linspace(0,1000);
t1 = linspace(10,1000,90); % time interval for T_wall definition
Twall = 393 + 100*square(100*t1);
Tw = [293*ones(1,10) Twall]; % wall temperature (array)

figure(1)
plot(t,Tw)
xlabel('t [s]')
ylabel('Wall Temperature [K]')
title('Wall Temperature as a square wave')


%% Stephan number definition
St_s = cp_s.*(T_melt-T0)./L; % Stephan number solid phase [-]
St_l = cp_l.*(Tw-T_melt)./L; % Stephan number liquid phase [-] (array)

% Melting time definition (Solomon correlation)
t_star = 0.11+(0.25./St_l); % adimensional time t* liquid phase [-] (array)
t_melt = t_star.*(R_pcm^2./alpha_l); % melting time [s] (array)

% domain discretization
X0 = 0; % initial position
XL = R_pcm; % final position [m]
x = linspace(X0, XL); % space discretization [m]

nu = sqrt(alpha_l./alpha_s); % parameter

xi = 0.01*ones(1,length(x)); % initialization for xi (root of the implicit equation)

% Trascendental equation for 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

%% Newton-Raphson method to find xi
tol = 1e-8; % initialization for tolerance 
xiold = 0.01; 
n = 0;

for z = 1:length(xi)
    err = 1; % initialization for the error
xi_step = xi(z);
    while err > tol
      xiold = xi_step;
      xi_step = xi_step - f(xi_step)./fd(xi_step); % xi solution
      err = abs(xi_step-xiold);
      n = n+1;
    end

xi(z) = xi_step;
clear err
clear xi_step
clear xiold
end
%% Position of the liquid-solid interface
Xf = @(t) 2*xi.*sqrt(alpha_l.*t);
X_melt = Xf(t_melt); % [m]
% Show results
fprintf('Solution value xi:')
disp(xi); 
fprintf('Melting time calculated [s]:')
disp(t_melt);
fprintf('Melting position calculated [m]:')
disp(X_melt);

% plot X_melt(t)
plot(t,real(X_melt),'Color','b');

[Torsten]

So you really want to replace the in-built MATLAB functions erf for the error function and erfc for the complementary error function with your own functions

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

?

Answer 1

[Askic V]

Is this what you wanted to achieve?

clear
clc
close all
a0 = -0.00401;
a1 = 1.18669;
a2 = -0.14559;
a3 = -0.33443;
a4 = 0.16069;
a5 = -0.02155;

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.05; % PCM length [m] (supposed)
T0 = -10+273; % initial temperature of the PCM [K] (assumed from Mansouri's paper)
Tamb = 20+273; % ambient temperature [K]
% input data for the IMD (geometrical properties)
D_imd = 254*(10^-3); % diamiter [m]
h_imd = 186*(10^-3); % height [m]
% Wall temperature definition (square wave)
t = linspace(0,1000);
t1 = linspace(10,1000,90); % time interval for T_wall definition
Twall = 393 + 100*square(100*t1);
Tw = [293*ones(1,10) Twall]; % wall temperature (array)
St_s = cp_s.*(T_melt-T0)./L; % Stephan number solid phase [-]
St_l = cp_l.*(Tw-T_melt)./L; % Stephan number liquid phase [-] (array)
% Melting time definition (Solomon correlation)
t_star = 0.11+(0.25./St_l); % adimensional time t* liquid phase [-] (array)
t_melt = t_star.*(R_pcm^2./alpha_l); % melting time [s] (array)
% domain discretization
X0 = 0; % initial position
XL = R_pcm; % final position [m]
x = linspace(X0, XL); % space discretization [m]
nu = sqrt(alpha_l./alpha_s); % parameter

syms xi % symbolic variable

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 = (St_l./(exp(xi.^2).*erf(xi)))-(St_s./((nu*exp(nu^2.*xi.^2).*erfc(nu.*xi))))-(xi.*sqrt(pi));
df = diff(f); % derivative of the function f

xi_sol = -1*ones(1,length(x));

tol = 1e-6;

err = 1; % initialization for the error

for i = 1:length(xi_sol)
    xi = xi_sol(i); % initial solution for each element in xi_sol
    while err > tol
        xi_n = xi - double(subs(f, xi))/double(subs(df,xi)); % xi solution
        err = abs(xi_n-xi);
        xi = xi_n;
    end
    
    xi_sol(i) = xi;
end

% Position of the liquid-solid interface
Xf = @(t) 2*xi_sol.*sqrt(alpha_l.*t);
X_melt = Xf(t_melt); % [m]
% Show results
%fprintf('Solution value xi:')
%disp(xi_sol); 
%fprintf('Melting time calculated [s]:')
%disp(t_melt);
%fprintf('Melting position calculated [m]:')
%disp(X_melt);
% plot X_melt(t)
plot(t,real(X_melt),'Color','b');