MATLAB Newton Raphson method with array - Stephan problem

2 次查看（过去 30 天）
Elisa Revello
Elisa Revello 2023-3-9
编辑： Torsten 2023-3-9
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 个评论
Torsten
Torsten 2023-3-9
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
?

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Askic V
Askic V 2023-3-9
编辑：Torsten 2023-3-9
Ran in:
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 个）

类别

在 Help Center 和 File Exchange 中查找有关 Thermal Analysis 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by