1D melting problem with Neumann's analytical solution

11 次查看（过去 30 天）
Elisa Revello
Elisa Revello 2023-1-11
评论： Elisa Revello 2023-2-5
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 个评论
Image Analyst
Image Analyst 2023-1-11
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
Elisa Revello 2023-1-12
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));

请先登录，再进行评论。

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

采纳的回答

Alan Stevens
Alan Stevens 2023-1-12
编辑：Alan Stevens 2023-1-12
Ran in:
  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
  3 个评论
显示 1更早的评论
Alan Stevens
Alan Stevens 2023-1-17
Ran in:
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
Elisa Revello 2023-2-5
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');

请先登录，再进行评论。

更多回答（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