# Moving boundary problem with Stefan condition (modeling mass and heat transfer for sublimation front)

29 views (last 30 days)
Jessica Hoffman on 20 Aug 2018
Commented: CARLOS HOUNGBEME on 24 Aug 2021
Governing Equations and BCs attached.
Problem: I am trying to model 1D mass and heat transfer for sublimation with a porous,dried media (region I) through which gas flows and a frozen, solid section (region II), with a sublimation front at the interface. This imposes Stefan condition at the moving boundary. I understand the derivation and math, but I am new to MATLAB, and get so confused when I try to convert equations to it. The code I have attached here is garbage, but I wanted to show in case it better illustrates where my confusion lies.
I have performed coordinate transformation on the governing equations using a non-dimensionalized scaling factor, discretized in space using finite difference method, and tried to couple the equations for the 2 regions with the boundary equation (Stefan condition). I really could use some guidance from someone well-versed in MATLAB!
*I would use the following code as a function and then in command window use: [t,y]=ode15s(@solid_process,[0 100],[0 0 0]).*
.m code:
% this function sets up governing equations and boundary condition for solid
% process with mvoing boundary
% define function
function [systemofodes]=solid_process()
%--------params------------------
L = 8;
N = 51;
dx = L/N;
rho = 1;
del = 0.001;
kl = .0004;
alphal = 1e-8;
%-----------system of ODEs after Method of Lines-------------
dTk/dt = alphal*[(T(k+1)-2*T(k)+T(k-1))/(dx^2)];
dTi/dt = alphal*[(T(m)-2*T(i)+T(i-2))/(dx^2)] + alphal/kl*[(1-del)*rho*L/dx*ds/dt];
ds/dt = -kl/rho/L*[(T(i)-T(i-2))/2/dx + (1+del)*(T(i-2)-2*T(i-2)+T(I))/dx];
% set each of ODEs as item in vector
y(0) = dTk/dt; % IC: @t=0, y(0)=0
y(1) = dTi/dt; % IC: @t=0, y(1)=0
y(2) = ds/dt; % IC: @t=0, y(2)=0
%returns system of ODEs as array for solution in ode15s or other solver
systemofodes = [y];
end
Maria Lorena Richiusa on 1 May 2020
Hi Jessica,
have you managed to solve the Stefan problem? If yes, I would like to talk with you about it for knowing more.
I am running into it, trying to solve it by using the Goodman and Shea approach.
Many thanks,
Lorena

Bill Greene on 30 Aug 2018
I think you may be able to solve this problem using the pdepe function. The idea is that you model both phases as a single region with material properties that vary as a function of x. In particular, at the point in the region where the phase change is occurring, the latent heat associated with the phase change, is accounted for by adjusting the specific heat of the material. Voller refers to this as the "apparent heat capacity" method. See Voller Overview for more details.
CARLOS HOUNGBEME on 24 Aug 2021