1
$\begingroup$

I am trying to code the one-phase, one-dimensional Stefan problem using finite differences in Matlab, similarly to what has already been done in Mathematica (see https://mathematica.stackexchange.com/questions/58596/pde-with-stefan-conditions-a-k-a-variable-boundary).

The problem consists in finding $(u,s)$ satisfying the coupled PDE-ODE system

\begin{equation} \begin{cases} u_t - u_{xx} = 0 &\text{for} \ 0 < t < T, \ 0 < x < s(t),\\ u(0, t) = f(t) &\text{for} \ 0 < t < T, \\ u(s(t),t) = 0, &\text{for} \ 0 < t < T,\\ \dot{s}(t) = -u_x(s(t),t) &\text{for} \ 0 < t < T,\\ u(x,0)=u_0(x), \ s(0)=s_0, \end{cases} \end{equation} where $T, u_0, f, s_0$ are given (positive) data (to assure the positivity of the state $u$ in the set where it solves the heat equation). One assumes that $u=0$ elsewhere in space.

Conducting the change of variable $$\xi = \frac{x}{s(t)}, \ z(\xi,t) = u(x,t),$$ the governing equations hold in the fixed space domain $(0, 1)$: \begin{equation} \begin{cases} z_t - \frac{1}{s^2}z_{\xi\xi} - \frac{\xi \dot{s}}{s}z_\xi = 0 & (\xi, t) \in (0,1)\times(0,T), \\ \dot{s}(t) = -\frac{1}{s(t)}u_x(1,t) &t \in (0, T), \end{cases} \end{equation} now with a homogeneous Dirichlet boundary condition at $x=1$.

I used an explicit Euler scheme for discretizing $\dot{s}$, and an implicit scheme for the time derivative of the state $z$, and centered differences for its space derivatives. I chose the initial and boundary data $f \equiv 1, u_0 \equiv 0.1$ and $s_0 = 0.05$.

The problems I encountered are the following:

  • The solution seems to equal $0$ after a very short time, no matter the choice of the initial data;

  • The interface $s$ appears to have a non-smooth behavior when $u_0 > s_0$;

I would be very grateful if someone could help me out with the above issues. I am looking to obtain a similar output as the answer in Mathematica SE.

The code is below:

clc;

nx = 75; nT = 25;
dt = 1/nT; dx = 1/nx;
M1 = 1; M2 = 0;
s0 = 0.1;

U = zeros(nx+1,nT+1);
U(1,1) = s0; U(2:nx+1,1) = 0.05;

Id = eye(nx);
A = -2*Id + diag(ones(nx-1,1),-1) + diag(ones(nx-1,1),1); A = (1/dx^2)*A;
Ma = zeros(nx,1); Ma(1,1) = M1; Ma(nx,1) = M2; Ma = Ma/(dx^2);
G = diag(ones(1,nx-1),1) + diag(-ones(1,nx-1),-1); G = 1/(2*dx)*G;
Mg = zeros(nx,1); Ma(1,1) = -M1; Ma(nx,1) = M2; Mg = Mg/(2*dx);
x = linspace(0, 1, nx+2)';

for j=1:nT
    U(1,j+1) = (dt/U(1,j))*U(nx+1,j)/dx + U(1,j);
    ds = (U(1,j+1)-U(1,j))/dt;
    U(2:nx+1,j+1) = (U(1,j)^2*Id-dt*A-U(1,j)*ds*dt*x(2:nx+1).*G)\(dt*Ma+U(1,j)*ds+dt*x(2:nx+1).*Mg+U(2:nx+1,j));
end

s = U(1,1:nT+1);
tint = linspace(0,1,nT+2);
plot(tint(2:nT+2),s, 'r--');
hold on;

[tp, esp] = meshgrid(tint(2:nT+2), x(2:nx+1));
% Going back to the original variable
esp2 = esp*diag(s);
y = U(2:nx+1,:);
%contourf(tp,esp2,y);
surf(tp,esp2,y);
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.