Programming Homework 06

Submission Deadline: 22.07.2026 - 23:55
  • The homework solution has to be handed in as a group solution via Moodle.
  • The notebooks must run without errors in our environment on the RWTH-Cluster. See here for more information on how to get started.
Participant list

Please add your names and student ID (Matrikelnummer) here:

  1. Name, ID
  2. Name, ID
  3. Name, ID
  4. Name, ID

Learning goals

  1. Get a feeling for a fundamental 1D phase-change problem: The Stefan problem.
  2. Understand important components of an analytical solution for two different freezing / melting configurations
  3. Use numerical root finding methods to solve nonlinear equations
  4. Understand the role of latent heat and the Stefan number in the phase-change problem
  1. Make sure that you update the library folder in your CMM directory. Download link: library.tar.gz
  2. Make sure that you update your dolfinx_latest.sif container image.

You can use the following command on the RWTH cluster/Linux:

    wget https://mbd.pages.rwth-aachen.de/courses/cmm/content/exercises/notebooks/homework06.out.ipynb
    rm -rf library
    wget https://mbd.pages.rwth-aachen.de/courses/cmm/content/exercises/notebooks/library.tar.gz
    tar -xvzf library.tar.gz
    rm library.tar.gz
    rm dolfinx_latest.sif
    apptainer pull oras://registry.git.rwth-aachen.de/mbd/courses/containers/dolfinx:latest

Todo’s

Make sure that all tests in the notebook run without errors. The following tasks need to be completed:

Imports

import os
import sys
import numpy as np

from mpi4py import MPI
import matplotlib as mpl 

# Colorblind friendly color cycle for matplotlib:
mpl.rcParams['axes.prop_cycle'] = mpl.cycler(color=['#377eb8', '#ff7f00', '#4daf4a',
                  '#f781bf', '#a65628', '#984ea3',
                  '#999999', '#e41a1c', '#dede00'])


CMM_DIR = os.getcwd()
os.environ["PYTHONPATH"] = f"{CMM_DIR}:{os.environ.get('PYTHONPATH', '')}"

# overload print function to print only on rank 0
def print_rank0(*args, **kwargs):
    if MPI.COMM_WORLD.rank == 0:
        print(*args, **kwargs)


# Check if we are running interactively (e.g. Jupyter notebook with visual frontend)
# or as a script: https://stackoverflow.com/a/64523765
interactive_mode_is_on = hasattr(sys, 'ps1')

Background: The Stefan problem

In this homework we implement analytical solutions to the Stefan problem, which we introduced in the lecture. The Stefan problem is a fundamental model for 1D solid-liquid phase change: a moving phase-change interface (PCI) separates the solid from the liquid phase, and the interface position \(X_m(t)\) evolves in time as latent heat is absorbed or released at the interface.

We analyze the analytical solution for two different configurations:

  1. The single-phase Stefan problem (P1): a solid, initially at melting temperature, is molten from a heated wall. Only the liquid phase carries a nontrivial temperature profile.
  2. The full two-phase Stefan problem (P2): a liquid, initially above melting temperature, is frozen from a cooled wall. The temperature varies in both phases.

The similarity solution

In both configurations, the position of the phase-change interface follows a similarity solution of the form

\[ X_m(t) = 2\,\sqrt{\alpha\,t}\,\lambda \quad , \]

where \(\alpha\) is a thermal diffusivity and \(\lambda\) is a dimensionless similarity parameter that links the time \(t\) to the position of the phase-change interface \(X_m(t)\) at that time. Once \(\lambda\) is known, we can express the rest of the analytical solution for the temperature profile in terms of it. Finding \(\lambda\) requires solving a nonlinear equation numerically — this is where the numerical root finding of this homework comes in.

A central role is played by the Stefan number

\[ \textit{Ste} = \frac{c_{p,l}\,\left(T_w - T_m\right)}{L} \quad , \]

which relates a characteristic amount of sensible heat (\(c_{p,l}\,\left(T_w - T_m\right)\)) to the latent heat \(L\) of the phase change. For the single-phase problem, the Stefan number is the only parameter that determines \(\lambda\); for the full two-phase problem, the material parameters of both phases enter as well.

P1 — The single-phase Stefan problem

We start with the setup from the lecture:

single phase stefan setup

In this setup, the solid, initially at melting temperature \(T_m\), is suddenly heated from the left boundary at \(x=0\) to a temperature \(T_{\text{wall}} > T_m\).

We are first interested in finding the position of the phase-change interface (PCI) \(X_m(t)\) as a function of time \(t\). For a given interface location we can then find an analytical expression for the temperature in the liquid phase.

“Single-phase” refers to the fact that the solid remains at a constant temperature throughout and we hence only have to worry about the heat transfer in the liquid (the flux on the solid side is zero).

In the following, we implement the analytical solution for the single-phase Stefan problem.

1. Material parameters

We start by defining liquid material parameters. Note that no solid-phase parameters are defined — indeed, the solid phase does not play a role in the analytical solution!

class LiquidMaterialParameters():
    """ 
    A small dataclass to hold the material parameters for the liquid phase.
    """
    specific_heat_capacity = 4200  # J / (kg * K)
    density = 1000  # kg / m³
    thermal_conductivity = 0.6  # W/(m * K)
    latent_heat = 333700  # J / kg
    thermal_diffusivity = thermal_conductivity / (density * specific_heat_capacity)

liquid_parameters = LiquidMaterialParameters()

2. The SinglePhaseStefanProblem class

We then set up a class for our SinglePhaseStefanProblem, which holds the material parameters and also implements functions to return the interface position \(X_m(t)\) as well as the analytical temperature solution \(T(x, t)\).

The interface position is computed from the similarity parameter \(\lambda\)

\[ X_m(t) = 2\,\sqrt{\alpha_l\,t}\,\lambda \]

so we will first need to tackle the problem of finding \(\lambda\).

class SinglePhaseStefanProblem():
    """Implements the analytical solution to the 1D single-phase Stefan problem."""

    def __init__(self, T_wall, T_m, liquid_material_params: LiquidMaterialParameters):
        cp_l = liquid_material_params.specific_heat_capacity
        hm = liquid_material_params.latent_heat
        k_l = liquid_material_params.thermal_conductivity
        rho_l = liquid_material_params.density

        self.Ste = cp_l * (T_wall - T_m) / hm 
        """ Stefan number """

        self.lam = None 
        """Similarity parameter lambda, corresponds to lambda from below (12.5) in slides"""

        self.T_wall = T_wall
        r"""Liquid / Wall temperature (temperature at `x=0`)"""
        self.T_m = T_m
        r"""Melting temperature (equals temperature at `x=\infty`)"""
        self.alpha_l = k_l / (rho_l * cp_l)


    def interface_location(self, t) -> float:
        r"""Compute phase-change interface (PCI) location :math:`x_{\text{PCI}}(t)` at
        time `t`.

        Evaluate Eqn. below (12.5) in slides.

        :param t: time
        :return: PCI location at time `t`
        """
        return float(2 * np.sqrt(self.alpha_l) * self.lam * np.sqrt(t))


    def temperature_analytical(self, T_liquid, x, t):
        """ 
        Return the temperature at position `x` and time `t` based on the analytical solution
        of the Stefan problem.  

        The function takes a numpy array `x` of positions and a scalar `t` for time.
        For each position `x`, it checks if it is less than the interface location at time `t`,
        and returns the temperature of the liquid phase if it is, or the melting temperature
        if it is not.

        :param T_liquid: Temperature of the liquid phase, given as a Python function of `x` and `t` 
                         that returns the liquid temperature at the given position and time.
        :param x: position in the domain (in dimensions)
        :param t: time (in dimensions)
        """
        interfaceLoc = self.interface_location(t)

        # Support both scalar and numpy array x
        x = np.asarray(x)
        interfaceLoc = np.asarray(interfaceLoc)
        # For each x, if x < interfaceLoc, use T_liquid(x, t), else use self.T_m
        mask = x < interfaceLoc
        result = np.where(mask, T_liquid(self, x, t), self.T_m)
        return result

We set up a case where the wall temperature is \(T_{\text{wall}} = 300\,\mathrm{K}\) and the melting temperature is \(T_m = 273\,\mathrm{K}\):

stefan_problem = SinglePhaseStefanProblem(
    T_wall=300,  # K
    T_m=273,  # K
    liquid_material_params=liquid_parameters
)

3. Finding the similarity parameter \(\lambda\)

A core realization of the derivation given in the lecture is that the analytical solution has a similarity parameter \(\lambda\) that links the time \(t\) to the position of the phase change interface \(X_m(t)\) at that time:

\[ \lambda = \frac{X_m(t)}{2\,\sqrt{\alpha_l\,t}} \]

As discussed in the background section, for the single-phase problem \(\lambda\) depends solely on the Stefan number \(\textit{Ste}\). The determining equation for \(\lambda\) is Eqn. (12.6) in the lecture script, which reads

\[ \textit{Ste} \, \frac{\text{exp}\left(-\lambda^2\right)}{\sqrt{\pi}\,\text{erf}\left(\lambda\right)} = \lambda \quad . \]

This equation can only be solved numerically and we can do so by applying a root-finding algorithm to find the \(\lambda\) that satisfies \(G(\lambda) = 0\), where

\[ G(\lambda) = \textit{Ste} \, \frac{\text{exp}\left(-\lambda^2\right)}{\sqrt{\pi}\,\text{erf}\left(\lambda\right)} - \lambda \overset{!}{=} 0 \quad . \]

STUDENT TODO — Task 1: Find the similarity parameter \(\lambda\)

Complete the find_similarity_parameter_lambda() function below.

  1. Code the function \(G(\lambda)\) as a Python lambda function.
  2. Use the scipy.optimize.fsolve function to find the root of \(G(\lambda)\), i.e. the value of \(\lambda\) that satisfies \(G(\lambda) = 0\), iteratively. You can find help on how to use fsolve here.

The find_similarity_parameter_lambda function should return the value of \(\lambda\) that satisfies \(G(\lambda) = 0\).

from scipy.optimize import fsolve
import numpy as np
from scipy.special import erf

def find_similarity_parameter_lambda(stefan_problem: SinglePhaseStefanProblem):
    """
    Find the similarity parameter lambda by solving the nonlinear equation (12.6)
    for lambda.
    """
    Ste = stefan_problem.Ste
    rootpi = np.sqrt(np.pi)

    eqn_12_6_script = lambda lam: 0 # YOUR CODE HERE
    lambda_sim = 0.0 # YOUR CODE HERE

    return lambda_sim

our_lambda = find_similarity_parameter_lambda(stefan_problem)
stefan_problem.lam = our_lambda
Verification

The cell below checks the computed similarity parameter \(\lambda\) against the expected value for the single-phase Stefan problem.

def test_lambda_solution():
    assert np.isclose(0.3914744836716143, our_lambda, rtol=1e-5), \
        f"Expected lambda to be 0.3914744836716143, but got {our_lambda}"
    
    print_rank0("TEST PASSED: similarity parameter is correct for the single phase Stefan problem.")
    
test_lambda_solution()
TEST PASSED: similarity parameter is correct for the single phase Stefan problem.

With that we can compute the interface location at a given time t from the interface_location(t) function of the SinglePhaseStefanProblem class. The interface location is given by the equation:

\[ X_m(t) = 2\,\sqrt{\alpha_l \, t}\,\lambda \]

We will later explain a bit what the \(\sqrt{\alpha_l \, t}\) term means. For now, we can just use it to compute the interface location at a given time t:

X_m = stefan_problem.interface_location(t)

4. Analytical solution in the liquid phase

The analytical solution for the 1D single-phase Stefan problem depends on the location \(x\) and time \(t\):

\[ \begin{cases} T(x, t) = T_{\text{liquid}}(x,t)= T_w - (T_w - T_m) \, \frac{\text{erf}\left(\frac{x}{2\,\sqrt{\alpha_l\,t}}\right)}{\text{erf}\left(\frac{X_m(t)}{2\,\sqrt{\alpha_l\,t}}\right)}, & x < X_m(t) \\ T(x, t) = T_m, & x \geq X_m(t) \end{cases} \]

STUDENT TODO — Task 2: Analytical solution in the liquid phase

Implement the first part of the analytical solution, the liquid temperature profile, in the function T_liquid(stefan_problem, x, t) template below. The function should return the analytical solution for the temperature in the liquid phase at a given position x and time t.

Note: the function T_liquid should accept an array of sample points x given as a numpy array and return the temperature at these points as a numpy array of the same shape.

Hint: We can get the interface location from the stefan_problem.interface_location(t) function.

from scipy.special import erf
from numpy import sqrt

def T_liquid(stefan_problem: SinglePhaseStefanProblem, x, t):
        alpha_l = stefan_problem.alpha_l
        T_wall = stefan_problem.T_wall
        T_m = stefan_problem.T_m
        # ...

        result = x # YOUR CODE HERE
        return result
Verification

The cell below evaluates T_liquid at three sample points (the wall, the midpoint, and the interface) and compares against the expected values.

def test_temperature_liquid():
    """
    Test the temperature function for the liquid phase.
    """
    t_eval = 7000000  # seconds
    interface_location = stefan_problem.interface_location(t_eval)
    x_eval = np.linspace(0, interface_location, 3) 

    check_result = np.array([stefan_problem.T_wall,  285.98631577, stefan_problem.T_m])
    assert np.allclose(
        T_liquid(stefan_problem, x_eval, t_eval), check_result, rtol=1e-5
    ), f"Expected {check_result}, but got {T_liquid(stefan_problem, x_eval, t_eval)}"

    print_rank0("TEST PASSED: temperature function for the liquid phase seems correct.")
    

test_temperature_liquid()
TEST PASSED: temperature function for the liquid phase seems correct.

5. The thermal diffusion timescale

We now have an analytical solution for the temperature field at a given time \(t\). We now try to find a sensible time scale at which to evaluate the problem for a given size \(L\) of the domain. In our case the governing PDE is given by the heat equation in the liquid:

\[ \frac{\partial T}{\partial t} - \alpha_l \Delta T = 0 \]

We can nondimensionalize this by introducing a temperature scale \(T_f\), a length scale \(L\) and a time scale \(t_f\), such that \[ \tilde{T} = \frac{T}{T_f}, \quad \tilde{x} = \frac{x}{L}, \quad \tilde{t} = \frac{t}{t_f} \] are the nondimensional temperature, position and time respectively.

After canceling out \(T_f\) and rearranging, the nondimensionalized PDE then reads

\[ \underbrace{\left(\frac{L^2}{\alpha_l \, t_f}\right) \frac{\partial \tilde{T}}{\partial \tilde{t}}}_{\text{transient term}} - \underbrace{\tilde{\Delta} \tilde{T}}_{\text{diffusive term}} = 0 \]

Aha! From that we see that the two terms are of the same order of magnitude if the Fourier number \(\textit{Fo} = \frac{L^2}{\alpha_l \, t_f}\) is of order one (in that case there is a “1” in front of both terms). When will the Fourier number be of order one? — good question!

The time scale \(t_f\) at which \(Fo = 1\) is exactly what we call the “thermal diffusion time scale”, e.g. \[ t_f = t_{\text{diff}} = \frac{L^2}{\alpha_l} \]

STUDENT TODO — Task 3: Compute the thermal diffusion time scale

We consider a characteristic length of \(L=1\,\mathrm{m}\). Compute the thermal diffusion time scale in seconds and assign it to the variable thermal_diffusion_timescale below.


thermal_diffusion_timescale = 1 # YOUR ASSIGNMENT HERE

print("Thermal diffusion timescale (in days):", thermal_diffusion_timescale / 3600 / 24)
Thermal diffusion timescale (in days): 1.1574074074074073e-05
Verification

The cell below checks the computed thermal diffusion timescale.

def test_diffusion_timescale():
    assert np.isclose(7000000.0, thermal_diffusion_timescale, atol=1e-10), \
        f"Expected thermal diffusion timescale to be 7000000.0, but got {thermal_diffusion_timescale}"
    
    print_rank0("TEST PASSED: thermal diffusion timescale is correct.")

test_diffusion_timescale()
TEST PASSED: thermal diffusion timescale is correct.

6. Visualization: temperature profiles over time

With that, we can plot the analytical temperature profile for \(t=\{0.1\,t_{\text{diff}}, 0.5\,t_{\text{diff}}, \,t_{\text{diff}}, 2\,t_{\text{diff}}\}\)

import matplotlib.pyplot as plt

t_diff = thermal_diffusion_timescale 
x = np.linspace(0, 1, 100)
plot_times = np.array([0.1, 0.5, 1, 2]) * t_diff

for time_snapshot in plot_times:
    temperature_profile = stefan_problem.temperature_analytical(T_liquid, x, time_snapshot)
    interface_location = stefan_problem.interface_location(time_snapshot)
    plt.axvline(interface_location, color='gray', linestyle='--', alpha = 0.6)
    plt.plot(x, temperature_profile, label=f't={time_snapshot / t_diff:.2f} $t_d$')

plt.xlabel('Position (m)')
plt.ylabel('Temperature (K)')
plt.title('Temperature Profile Over Time')
plt.legend()
plt.grid()

7. What if we increase the Stefan number?

Next, we investigate the effect of varying the Stefan number. For this problem, the Stefan number is defined as

\[ \textit{Ste} = \frac{c_{p,l}\,\left(T_w - T_m\right)}{L} \quad , \]

which is a ratio of a characteristic amount of sensible heat \(c_{p,l}\,\left(T_w - T_m\right)\) to latent heat \(L\) of the phase change.

Note that the sign is sometimes flipped or the Stefan number is defined as the inverse. Here we chose the sign to get a positive Stefan number for the case \(T_w > T_m\).

We saw earlier that it is the only parameter that determines the value of \(\lambda\), but how does the actual solution change according to that?

STUDENT TODO — Task 4: Modify the Stefan number

Based on the material parameters of our initial setup, modify the latent heat such that the resulting setup has a Stefan number 10 times lower (liquid_parameters_low_Ste) and 10 times higher (liquid_parameters_high_Ste) than the original setup.

from copy import deepcopy
liquid_parameters_low_Ste = deepcopy(liquid_parameters)
liquid_parameters_high_Ste = deepcopy(liquid_parameters)
# YOUR CODE HERE
# YOUR CODE HERE

The Stefan number has changed so we need to compute new \(\lambda\) similarity parameters for the problems with low and high Stefan number:

stefan_problem_low_ste = SinglePhaseStefanProblem(
    T_wall=300,  # K
    T_m=273,  # K
    liquid_material_params=liquid_parameters_low_Ste
)
# Recompute the similarity parameter lambda with the new latent heat
lambda_low_ste = find_similarity_parameter_lambda(stefan_problem_low_ste)
stefan_problem_low_ste.lam = lambda_low_ste


stefan_problem_high_ste = SinglePhaseStefanProblem(
    T_wall=300,  # K
    T_m=273,  # K
    liquid_material_params=liquid_parameters_high_Ste
)
lambda_high_ste = find_similarity_parameter_lambda(stefan_problem_high_ste)
stefan_problem_high_ste.lam = lambda_high_ste
Verification

The cell below checks that the modified setups indeed have a Stefan number 10 times lower and 10 times higher than the original setup.

def test_stefan_numbers():
    original_Ste = stefan_problem.Ste
    assert np.isclose(original_Ste / 10, stefan_problem_low_ste.Ste, rtol=1e-5), \
        f"Expected low Stefan number to be {original_Ste / 10}, but got {stefan_problem_low_ste.Ste}"
    assert np.isclose(original_Ste*10, stefan_problem_high_ste.Ste, rtol=1e-5), \
        f"Expected high Stefan number to be {original_Ste * 10}, but got {stefan_problem_high_ste.Ste}"
        
    print_rank0("TEST PASSED: Stefan numbers for low and high latent heat are correct.")

test_stefan_numbers()
TEST PASSED: Stefan numbers for low and high latent heat are correct.

Now we can plot the analytical temperature profiles at a given time for the three Stefan numbers:

import matplotlib.pyplot as plt

t_diff = thermal_diffusion_timescale 
x = np.linspace(0, 3, 200)
plot_time = 0.2 * t_diff

stefan_problems = [
    stefan_problem_low_ste,
    stefan_problem,
    stefan_problem_high_ste
]

for problem in stefan_problems:
    temperature_profile = problem.temperature_analytical(T_liquid, x, plot_time)
    interface_location = problem.interface_location(plot_time)
    plt.axvline(interface_location, color='gray', linestyle='--', alpha = 0.6)
    plt.plot(x, temperature_profile, label=f'Ste={problem.Ste:.2f}')

plt.xlabel('Position (m)')
plt.ylabel('Temperature (K)')
plt.title('Temperature Profile for different Stefan numbers')
plt.legend()
plt.grid()

P2 — The full two-phase Stefan problem

Until now we have only considered “single-phase” Stefan problem in the sense that nothing was happening in the solid phase. Due to the choice of initial and far-field condition, the heat flux in the solid phase was always zero and we only had to care about the liquid phase heat equation.

In the final part of the exercise, we show that the same analytical approach generalizes to the “full” Stefan problem, where temperature in both phases is allowed to vary.

Here, we assume a liquid phase that is initially at temperature \(T_0 = T_{\text{initial}} = 278\,K\). It is then suddenly cooled from the left side (\(x=0\)) to \(T_w = 268\,K\) below the freezing point (\(T_m = 273\,K\)). The right side has a homogeneous Neumann boundary condition.

In this scenario, we have to solve the heat equation in both phases and hence our expressions get more complicated.

What we want to show in the following is how we can build on the solution steps shown above and use a similar strategy to find the similarity parameter \(\lambda\) and the analytical solution for the temperature profile.

1. Material parameters

First, we need to define the solid phase parameters:

class SolidMaterialParameters():
    """ 
    A small dataclass to hold the material parameters for the solid phase.
    """
    
    specific_heat_capacity = 1877  # J / (kg * K)
    density = 1000  # kg / m³
    thermal_conductivity = 2.5  # W/(m * K)
    latent_heat = 333700  # J / kg
    thermal_diffusivity = thermal_conductivity / (density * specific_heat_capacity)

solid_parameters = SolidMaterialParameters()

2. The TwoPhaseStefanProblem class

The full TwoPhaseStefanProblem stores the material parameters for both phases and has a function to compute the interface location and temperature solution for our setup.

For the full Stefan problem in the setup discussed here, the interface location is based on the similarity solution

\[ X_m(t) = 2\,\sqrt{\alpha_s \, t}\,\lambda \]

Note that here the thermal diffusivity \(\alpha_s\) of the solid phase is used instead of the liquid phase diffusivity \(\alpha_l\).

class TwoPhaseStefanProblem():
    """Implements the analytical solution to the 1D Stefan problem.

    In this realization, the temperature in the solid phase is not assumed to be
    constant. For simplicity we fix the wall temperature and the initial temperature to
    be constant.

    Follow www.doi.org/10.1016/0017-9310(81)90062-4 and
    www.doi.org/10.1080/10407790.2023.2211231
    """

    def __init__(
        self,
        T_wall,
        T_m,
        T_initial,
        liquid_material_params: LiquidMaterialParameters,
        solid_material_params: SolidMaterialParameters,
    ):
        """
        :param T_wall: Liquid / Wall temperature (temperature at :math:`x=0`)
        :param T_m: Melting temperature
        :param T_initial: Initial temperature :math:`T(x,t=0)` (assumed constant throughout domain for simplicity)
        :param liquid_material_params: Liquid material parameters
        :param solid_material_params: Solid material parameters
        """
        self.cp_l = liquid_material_params.specific_heat_capacity
        self.L = liquid_material_params.latent_heat
        self.k_l = liquid_material_params.thermal_conductivity
        self.rho_l = liquid_material_params.density

        self.cp_s = solid_material_params.specific_heat_capacity
        self.k_s = solid_material_params.thermal_conductivity
        self.rho_s = solid_material_params.density

        self.T_wall = T_wall
        r"""Liquid / Wall temperature (temperature at :math:`x=0`)"""
        self.Tm = T_m
        r"""Melting temperature (equals temperature at :math:`x=\infty`)"""
        self.T_initial = T_initial
        r"""Initial temperature :math:`T(x,t=0)` (considered constant here for
        simplicity)"""
        
        self.alpha_l = self.k_l / (self.rho_l * self.cp_l)
        self.alpha_s = self.k_s / (self.rho_s * self.cp_s)

        self.lam = None
        """Similarity parameter lambda"""

    def interface_location(self, t) -> float:
        r"""Compute phase-change interface (PCI) location :math:`x_{\text{PCI}}(t)` at
        time `t`.

        :param t: time
        :return: PCI location at time `t`
        """
        return float(2 * self.lam * np.sqrt(self.alpha_s * t))

    def temperature_analytical(self, T_liquid, T_solid,  x, t):
        """ 
        As before the temperature at position `x` and time `t` based on the analytical solution
        of the Stefan problem.  

        This time, we also need the temperature of the solid phase, which is not constant!

        :param T_liquid: Temperature of the liquid phase, given as a Python function of `x` and `t` 
                         that returns the liquid temperature at the given position and time.
        :param T_solid: Temperature of the solid phase, given as a Python function of `x` and `t`
                         that returns the solid temperature at the given position and time.
        :param x: position in the domain (in dimensions)
        :param t: time (in dimensions)
        """
        interfaceLoc = self.interface_location(t)

        # Support both scalar and numpy array x
        x = np.asarray(x)
        interfaceLoc = np.asarray(interfaceLoc)
        # For each x, if x < interfaceLoc, use T_liquid(x, t), else use self.T_m
        mask = x < interfaceLoc
        result = np.where(mask, T_solid(self, x, t), T_liquid(self, x, t))
        return result

We initialize the problem with our setup parameters:

stefan_problem_full = TwoPhaseStefanProblem(
    T_wall=268,  # K
    T_m=273,  # K
    T_initial=278,  # K
    liquid_material_params=liquid_parameters,
    solid_material_params=solid_parameters
)

3. Finding the similarity parameter \(\lambda\) for the full Stefan problem

We can find the similarity parameter \(\lambda\) for the full Stefan problem in a similar way as before.

The equation for the similarity parameter \(\lambda\) is much more complicated than before, so we provide it as a preimplemented function eqn_lambda.

STUDENT TODO — Task 5: Find the similarity parameter for the full Stefan problem

Complete the find_similarity_parameter_lambda_full_stefan() function below by solving eqn_lambda numerically. Return the resulting float value of \(\lambda\).

from scipy.optimize import fsolve
import numpy as np
from scipy.special import erf, erfc
from cmath import pi 

def find_similarity_parameter_lambda_full_stefan(stefan_problem: TwoPhaseStefanProblem):
    """
    Find the similarity parameter lambda by solving the nonlinear lambda equation for the full Stefan problem
    """

    rootpi = np.sqrt(pi)
    root_alpha_l = np.sqrt(stefan_problem.alpha_l)
    root_alpha_s = np.sqrt(stefan_problem.alpha_s)

    L = stefan_problem.L
    cp_s = stefan_problem.cp_s
    k_l = stefan_problem.k_l
    k_s = stefan_problem.k_s
    T_wall = stefan_problem.T_wall
    T_m = stefan_problem.Tm
    T_initial = stefan_problem.T_initial
    alpha_s = stefan_problem.alpha_s
    alpha_l = stefan_problem.alpha_l

    eqn_lambda = (
            lambda lam: -((lam * L * rootpi) / (cp_s * (T_m - T_wall)))
            + np.exp(-(lam**2)) / erf(lam)
            - (k_l / k_s)
            * (
                root_alpha_s
                * (T_initial - T_m)
                * np.exp(-(alpha_s / alpha_l) * lam**2)
            )
            / (root_alpha_l * (T_m - T_wall) * erfc(lam * np.sqrt(alpha_s / alpha_l)))
        )
    
    lambda_sim = 0 # YOUR CODE HERE

    return lambda_sim

our_lambda = find_similarity_parameter_lambda_full_stefan(stefan_problem_full)
stefan_problem_full.lam = our_lambda
Verification

The cell below checks the computed similarity parameter \(\lambda\) against the expected value for the full Stefan problem.

def test_lambda_full_stefan_solution():
    assert np.isclose(0.11020095191626643, stefan_problem_full.lam, rtol=1e-5), \
        f"Expected lambda to be 0.11020095191626643, but got {stefan_problem_full.lam}"
    
    print_rank0("TEST PASSED: similarity parameter is correct for the full Stefan problem.")

test_lambda_full_stefan_solution()
TEST PASSED: similarity parameter is correct for the full Stefan problem.

4. Analytical temperature profiles in both phases

What remains is to implement the temperature solution in the two phases. This time, we do not get away with a single expression for the temperature profile, but have to distinguish between the two phases.

In the solid phase (\(x < X_m(t)\)), the temperature profile is given by

\[ T_{\text{solid}}(x, t) = \frac{T_m - T_\text{wall}}{\text{erf}(\lambda)} \cdot \text{erf}\left(\frac{x}{2\sqrt{\alpha_s t}}\right) + T_\text{wall} \quad . \]

The temperature profile in the liquid phase (\(x > X_m(t)\)) is given by

\[ T_{\text{liquid}}(x, t) = T_\text{initial} - \frac{T_\text{initial} - T_m}{\text{erfc}\left(\lambda \sqrt{\frac{\alpha_s}{\alpha_l}}\right)} \cdot \text{erfc}\left(\frac{x}{2\sqrt{\alpha_l t}}\right) \quad . \]

Note: the functions T_solid and T_liquid below redefine the single-phase T_liquid from P1 — from here on, T_liquid refers to the two-phase solution. If we want to rerun the P1 cells, we need to rerun the notebook from the top.

STUDENT TODO — Task 6: Two-phase analytical temperature profiles

Implement the two functions T_solid(stefan_problem, x, t) and T_liquid(stefan_problem, x, t) to return the analytical temperature profile in the solid and liquid phase respectively.

Note: the functions should work given sample points \(x\) as a numpy array and return the temperature profile as a numpy array of the same shape.

from scipy.special import erf, erfc

def T_solid(stefan_problem: TwoPhaseStefanProblem, x, t):
        alpha_s = stefan_problem.alpha_s
        T_wall = stefan_problem.T_wall
        T_m = stefan_problem.Tm
        lam = stefan_problem.lam

        result = x # YOUR CODE HERE
        return result

def T_liquid(stefan_problem: TwoPhaseStefanProblem, x, t):
        alpha_l = stefan_problem.alpha_l
        alpha_s = stefan_problem.alpha_s
        T_m = stefan_problem.Tm
        T_initial = stefan_problem.T_initial
        lam = stefan_problem.lam

        result = x # YOUR CODE HERE
        return result
Verification

The cells below evaluate T_solid and T_liquid at sample points and compare against the expected values.

def test_T_solid_analytical():
    """
    Test the temperature function for the solid phase.
    """

    t_eval = thermal_diffusion_timescale
    interface_location = stefan_problem_full.interface_location(t_eval)
    x_eval = np.linspace(0, interface_location, 3) 
    
    check_result = np.array([stefan_problem_full.T_wall,  270.50758629, stefan_problem_full.Tm])

    assert np.allclose(
        T_solid(stefan_problem_full, x_eval, t_eval), check_result, rtol=1e-5
    ), f"Expected {check_result}, but got {T_solid(stefan_problem_full, x_eval, t_eval)}"

    print_rank0("TEST PASSED: temperature function for the solid phase seems correct.")

def test_T_liquid_analytical():
    """
    Test the temperature function for the liquid phase
    """

    t_eval = thermal_diffusion_timescale
    interface_location = stefan_problem_full.interface_location(t_eval)
    x_eval = np.linspace(interface_location, 3*interface_location, 3) 

    check_result = np.array([stefan_problem_full.Tm,  275.30962335, 276.79050018])

    assert np.allclose(
        T_liquid(stefan_problem_full, x_eval, t_eval), check_result, rtol=1e-5
    ), f"Expected {check_result}, but got {T_liquid(stefan_problem_full, x_eval, t_eval)}"

    print_rank0("TEST PASSED: temperature function for the liquid phase seems correct.")

test_T_solid_analytical()
test_T_liquid_analytical()
TEST PASSED: temperature function for the solid phase seems correct.
TEST PASSED: temperature function for the liquid phase seems correct.

5. Visualization: the full analytical solution

The full analytical solution is then given by combining the two phases, \[ T(x, t) = \begin{cases} T_{\text{solid}}(x, t), & x < X_m(t) \\ T_{\text{liquid}}(x, t), & x \geq X_m(t) \end{cases} \quad , \] which is implemented in the temperature_analytical function of the TwoPhaseStefanProblem class.

In the following, we plot the analytical solution for the temperature profile at different times:

import matplotlib.pyplot as plt

t_diff = thermal_diffusion_timescale 
x = np.linspace(0, 1, 100)
plot_times = np.array([0.01, 0.1, 0.5, 1]) * t_diff

for time_snapshot in plot_times:
    temperature_profile = stefan_problem_full.temperature_analytical(T_liquid, T_solid, x, time_snapshot)
    interface_location = stefan_problem_full.interface_location(time_snapshot)
    plt.axvline(interface_location, color='gray', linestyle='--', alpha = 0.6)
    plt.plot(x, temperature_profile, label=f't={time_snapshot / t_diff:.2f} $t_d$')

plt.xlabel('Position (m)')
plt.ylabel('Temperature (K)')
plt.title('Temperature Profile Over Time')
plt.legend()
plt.grid()

Prepare for your presentation

During the presentation, be prepared to answer the following questions/tasks. We will ask 1-2 questions per person.

For questions that specifically ask for equations, please be ready to write them down by hand without additional notes

  • What was the exercise about? Give a short summary and name the tasks that you had to solve.
  • Sketch the setup for the Stefan problem (single phase or full two-phase). Be prepared to sketch the initial condition and the solution (qualitatively) at a given time \(t > 0\).
  • Explain the meaning of the similarity parameter \(\lambda\). How can you compute \(X_m(t)\) from \(\lambda\)?
  • What is the thermal diffusion time scale \(t_{\text{diff}}\) and how can you compute it from the characteristic scales of the problem?
  • How is the Stefan number defined in this exercise? Explain its interpretation. Which quantities does it compare?
  • How does the PCI (phase change interface) location \(X_m(t)\) for a fixed sample time \(T\) vary if the Stefan number is increased / decreased?
  • How does the single-phase Stefan problem differ from the full two-phase Stefan problem?