MATLAB Newton Raphson method with array - Stephan problem
5 views (last 30 days)
Show older comments
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');
1 Comment
Torsten
on 9 Mar 2023
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
?
Accepted Answer
Askic V
on 9 Mar 2023
Edited: Torsten
on 9 Mar 2023
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');
0 Comments
More Answers (0)
See Also
Categories
Find more on Thermal Analysis 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!