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

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

[Torsten]

For anyone who tries to give advice it would be helful if you could attach the set of equations you are trying to solve together with initial and boundary conditions.

Best wishes

Torsten.

[Jessica Hoffman]

I have attached them! Thank you, Torsten.

[Maria Lorena Richiusa]

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

Answer 1

[Bill Greene]

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]

It difficult for me please help me