1D melting problem with Neumann's analytical solution
12 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Elisa Revello
el 11 de En. de 2023
Comentada: Elisa Revello
el 5 de Feb. de 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 comentarios
Image Analyst
el 11 de En. de 2023
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:
Respuesta aceptada
Alan Stevens
el 12 de En. de 2023
Editada: Alan Stevens
el 12 de En. de 2023
- eps is a built-in constant so better not to use it as a variable.
- Your png files show it is xi not epsilon.
- Try the following
- 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)
disp(steps)
3 comentarios
Alan Stevens
el 17 de En. de 2023
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:')
disp(xi); % root of the implicit equation
% 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]:')
disp(t_melt);
fprintf('Melting position calculated [mm]:')
disp(X_melt);
fprintf('PCM mass calculated [kg]:')
disp(m_pcm);
% 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);
Más respuestas (0)
Ver también
Categorías
Más información sobre Quantum Chemistry en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!