Simulate freezing

We need a few imports:

import matplotlib.pyplot as plt
from matplotlib import rc as matplotlibrc
from ruamel.yaml import YAML
yaml = YAML()

from lyopronto import freezing, plot_styling

Then, we provide all the necessary simulation parameters.

# Set up the simulation settings
# This needs to be a dict with string keys, which can be expressed in YAML as well
sim = yaml.load("""
tool: Freezing Calculator
Kv_known: Y
Rp_known: Y
Variable_Pch: N
Variable_Tsh: N
""")
# Or, equivalently:
sim = {
    'tool': 'Freezing Calculator',
    'Kv_known': 'Y',
    'Rp_known': 'Y',
    'Variable_Pch': 'N',
    'Variable_Tsh': 'N'
}


# Vial and fill properties
vial = {
    # Av = Vial area in cm^2
    'Av': 3.80,
    # Ap = Product Area in cm^2
    'Ap': 3.14,
    # Vfill = Fill volume in mL
    'Vfill': 2.0
}

# Product properties for freezing
product = {
    # cSolid = Fractional concentration of solute in the frozen solution
    'cSolid': 0.0,
    # Initial product temperature
    'Tpr0': 15.8,  # in deg C
    # Freezing temperature
    'Tf': -1.54,  # in deg C
    # Nucleation temperature
    'Tn': -5.84,  # in deg C
    # Critical product temperature
    # At least 2 to 3 deg C below collapse or glass transition temperature
    'T_pr_crit': -5  # in deg C
}

# Shelf temperature profile
Tshelf = {
    # init = Initial shelf temperature in C
    'init': 15.0,
    # setpt = Shelf temperature set points in C
    'setpt': [-40.0],
    # dt_setpt = Time for which shelf temperature set points are held in min
    'dt_setpt': [180.0],
    # ramp_rate = Shelf temperature ramping rate in C/min
    'ramp_rate': 1.0
}

# Time step
dt = 0.01    # hr

# Equipment capability parameters (not used for freezing, but included for completeness)
eq_cap = {
    'a': -0.182,  # kg/hr
    'b': 11.7     # kg/hr/Torr
}

# Number of vials (not used for freezing, but included for completeness)
nVial = 398

# Freezing heat transfer coefficient
h_freezing = 38.0  # W/m^2/K

Now, we are ready to actually run the simulation, which is lyopronto.freezing.freeze.

Info

Here is the docstring for freezing.freeze:

Simulates the primary drying process for a vial.

Parameters:

Name	Type	Description	Default
vial	dict	Vial properties, including 'Vfill', 'Ap', and 'Av'..	required
product	dict	Product properties, including 'cSolid', initial temperature 'Tpr0', freezing temperature 'Tf', and nucleation temperature 'Tn'.	required
h_freezing	float	Heat transfer coefficient during freezing [W/m²/K].	required
Tshelf	dict	Shelf temperature set points and time (see docs).	required
dt	float	Time step for the simulation [hr]. Used only as a sampling rate for output.	required

Returns:

Type	Description
	np.ndarray: An array containing the time, shelf temperature, and product temperature at each time step.

Source code in lyopronto/freezing.py

def freeze(vial,product,h_freezing,Tshelf,dt):
    """Simulates the primary drying process for a vial.

    Args:
        vial (dict): Vial properties, including 'Vfill', 'Ap', and 'Av'..
        product (dict): Product properties, including 'cSolid', initial temperature 'Tpr0', freezing temperature 'Tf', and nucleation temperature 'Tn'.
        h_freezing (float): Heat transfer coefficient during freezing [W/m²/K].
        Tshelf (dict): Shelf temperature set points and time (see docs).
        dt (float): Time step for the simulation [hr]. Used only as a sampling rate for output.

    Returns:
        np.ndarray: An array containing the time, shelf temperature, and product temperature at each time step.
    """

    ##################  Initialization ################

    # Initial fill height
    Lpr0 = functions.Lpr0_FUN(vial['Vfill'],vial['Ap'],product['cSolid'])   # [cm]

    # Frozen product volume
    V_frozen = Lpr0*vial['Ap']    # [mL]

    # Initialization of time
    iStep = 0      # Time iteration number
    t = 0.0    # Time [hr]

    # Initial shelf temperature
    Tsh = Tshelf['init']        # [degC]

    # Shelf temperature and time triggers, ramping rates
    Tshr = RampInterpolator(Tshelf)
    Tsh_tr = Tshr.values
    t_tr = Tshr.times
    r = np.array([[0.0]])    # [degC/min]
    for i,T in enumerate(Tsh_tr[:-1]):
        if Tsh_tr[i+1]>T:
            r = np.append(r,Tshelf['ramp_rate'])    # [degC/min]
        elif Tsh_tr[i+1]<T:
            r = np.append(r,-Tshelf['ramp_rate'])    # [degC/min]
        else:
            r = np.append(r,0.0)    # [degC/min]

    # Initial product temperature
    Tpr = product['Tpr0']    # [degC]
    Tpr0 = Tpr
    i_prev = 1    

    ######################################################

    freezing_output_saved = np.array([[t, Tsh, Tpr]])

    ################ Cooling ######################

    while(Tpr>product['Tn']): # Till the product reaches the nucleation temperature

        iStep = iStep + 1 # Time iteration number
        t = iStep*dt # [hr]

        if np.all(t_tr<t):
            warn("Total time exceeded. Freezing incomplete, no nucleation occurred")    # Shelf temperature set point time exceeded, freezing not done
            return freezing_output_saved
        else:
            i = np.argmax(t_tr>t) # Get first index where time trigger exceeds current time
            if not(i == i_prev):
                Tpr0 = Tpr
                i_prev = i
            # Evaluate shelf temperature at current time point
            Tsh = Tshr(t)
            # Product temperature
            Tpr = functions.lumped_cap_Tpr_sol(t-t_tr[i-1],Tpr0,vial['Vfill'],h_freezing,vial['Av'],Tsh,Tsh_tr[i-1],r[i])    # [degC]

        # Update record as functions of the cycle time
            freezing_output_saved = np.append(freezing_output_saved, [[t, Tsh, Tpr]],axis=0)    

    ######################################################

    ################ Nucleation ######################

    freezing_output_saved = np.append(freezing_output_saved, [[t, Tsh, product['Tn']]],axis=0)

    ######################################################

    ################ Crystallization ######################

    tn = t    # Nucleation onset time [hr]
    dt_crystallization = functions.crystallization_time_FUN(vial['Vfill'],h_freezing,vial['Av'],product['Tf'],product['Tn'],Tshr, tn)    # Crystallization time [hr]
    ts = tn + dt_crystallization    # Solidification onset time [hr]

    while(t<ts):

        if np.all(t_tr<t):
            warn("Total time exceeded. Freezing incomplete, nucleated but not fully crystallized")    # Shelf temperature set point time exceeded, freezing not done
            return freezing_output_saved
        else:
            i = np.argmax(t_tr>t) # Get first index where time trigger exceeds current time
            if not(i == i_prev):
                i_prev = i
            # Evaluate shelf temperature at current time point 
            Tsh = Tshr(t)    # [degC]
            # Product temperature stays at freezing temperature
            Tpr = product['Tf']    # [degC]

        # Update record as functions of the cycle time
            freezing_output_saved = np.append(freezing_output_saved, [[t, Tsh, Tpr]],axis=0)

        iStep = iStep + 1 # Time iteration number
        t = iStep*dt # [hr]    

    ######################################################

    ################ Solidification ######################

    t_last = ts
    Tpr0 = Tpr    # [degC]
    Tsh0 = Tsh
    while(t<t_tr[-1]):

        i = np.argmax(t_tr>t) # Get first index where time trigger exceeds current time
        if not(i == i_prev):
            i_prev = i
            t_last = t_tr[i-1]
            Tpr0 = Tpr
            Tsh0 = Tsh

        # Evaluate shelf temperature at current time point 
        Tsh = Tshr(t)    # [degC]

        # Product temperature
        Tpr = functions.lumped_cap_Tpr_ice(t-t_last,Tpr0,V_frozen,h_freezing,vial['Av'],Tsh,Tsh0,r[i])
        # Update record as functions of the cycle time
        freezing_output_saved = np.append(freezing_output_saved, [[t, Tsh, Tpr]],axis=0)

        iStep = iStep + 1 # Time iteration number
        t = iStep*dt # [hr]

    ######################################################

    return freezing_output_saved    

output_table = freezing.freeze(vial, product, h_freezing, Tshelf, dt)

Now, let's plot the results.

# Figure dimensions
figwidth = 30
figheight = 20
# Line width
lineWidth = 5
# Font family
plt.rcParams['font.family'] = 'Arial'

# Plot product and shelf temperature vs time
fig, ax = plt.subplots(figsize=(figwidth, figheight))

ax.plot(
    output_table[:, 0],  # Time [hr]
    output_table[:, 2],  # Product Temperature [°C]
    linewidth=lineWidth,
    label='Product Temperature'
)

ax.plot(
    output_table[:, 0],  # Time [hr]
    output_table[:, 1],  # Shelf Temperature [°C]
    linewidth=lineWidth,
    label='Shelf Temperature',
    linestyle='--'
)

plot_styling.axis_style_temperature(ax)

plt.legend(fontsize=40,loc='best')
ll,ul = ax.get_ylim()
ax.set_ylim([ll,ul+5.0])
plt.tight_layout()