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

29 views (last 30 days)
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()
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];
Maria Lorena Richiusa
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,

Sign in to comment.

Answers (1)

Bill Greene
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.


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!