Simulate primary drying with known Kv and Rp

Since this documentation example is a Jupyter notebook, inside the LyoPRONTO file structure, it needs to be directed to the LyoPRONTO code, which means adding ../../ to sys.path. However, this step will in general not be necessary.

# import sys
# sys.path.append('../../')

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 calc_knownRp, 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: Primary Drying Calculator
Kv_known: Y
Rp_known: Y
Variable_Pch: N
Variable_Tsh: N
""")
# Or, equivalently:
sim = {
    'tool':'Primary Drying 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
product = {
    # cSolid = Fractional concentration of solute in the frozen solution
    'cSolid': 0.05,
    # Product Resistance Parameters
    'R0': 1.4, # cm^2-hr-Torr/g
    'A1': 16.0, # cm-hr-Torr/g
    'A2': 0.0, # 1/cm
    # Critical product temperature
    # At least 2 to 3 deg C below collapse or glass transition temperature
    'T_pr_crit': - 5 # in deg C
}

# Vial Heat Transfer Parameters
ht = {
    'KC': 2.75e-4, # cal/s-cm^2-K
    'KP': 8.93e-4, # cal/s-cm^2-K-Torr
    'KD': 0.46     # 1/Torr
}

# Chamber Pressure
Pchamber = {
    # setpt = Chamber pressure set points in Torr
    'setpt': [0.15],
    # dt_setpt = Time for which chamber pressure set points are held in min
    'dt_setpt': [1800.0],
    # ramp_rate = Chamber pressure ramping rate in Torr/min
    'ramp_rate': 0.5
}

Tshelf = {
    # init = Intial shelf temperature in C
    'init': -35.0,
    # setpt = Shelf temperature set points in C
    'setpt': [20.0],
    # dt_setpt = Time for which shelf temperature set points are held in min
    'dt_setpt': [1800.0],
    # ramp_rate = Shelf temperature ramping rate in C/min
    'ramp_rate': 1.0
}

# Time step
dt = 0.01    # hr

Now, we are ready to actually run the simulation, which is lyopronto.calc_knownRp.dry.

Info

Here is the docstring for calc_knownRp.dry:

Simulate the primary drying process for known condiditions and parameters.

Parameters:

Name	Type	Description	Default
vial	dict	see master simulation inputs	required
product	dict	see master simulation inputs	required
ht	dict	see master simulation inputs	required
Pchamber	dict	see master simulation inputs	required
Tshelf	dict	see master simulation inputs	required
dt	float	Fixed time step for output [hours]	required

Returns:

Name	Type	Description
output_table	ndarray	Simulation output table with columns for: 0. Time [hr], 1. Sublimation front temperature [°C], 2. Vial bottom temperature [°C], 3. Shelf temperature [°C], 4. Chamber pressure [mTorr], 5. Sublimation flux [kg/hr/m²], 6. Drying percent [%]

Source code in lyopronto/calc_knownRp.py

def dry(vial,product,ht,Pchamber,Tshelf,dt):
    """Simulate the primary drying process for known condiditions and parameters.

    Args:
        vial (dict): see master simulation inputs
        product (dict): see master simulation inputs
        ht (dict): see master simulation inputs
        Pchamber (dict): see master simulation inputs
        Tshelf (dict): see master simulation inputs
        dt (float): Fixed time step for output [hours]

    Returns:
        output_table (ndarray): Simulation output table with columns for:
            0. Time [hr],
            1. Sublimation front temperature [°C],
            2. Vial bottom temperature [°C],
            3. Shelf temperature [°C],
            4. Chamber pressure [mTorr],
            5. Sublimation flux [kg/hr/m²],
            6. Drying percent [%]
    """

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

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

    # Time-dependent functions for Pchamber and Tshelf, take time in hours
    Pch_t = functions.RampInterpolator(Pchamber)
    Tsh_t = functions.RampInterpolator(Tshelf)

    # Get maximum simulation time based on shelf and chamber setpoints
    # This may not really be necessary, but is part of legacy behavior
    # Could remove in a future release
    max_t = max(Pch_t.max_time(), Tsh_t.max_time())   # hr, add buffer

    if Pch_t.max_setpt() > functions.Vapor_pressure(Tsh_t.max_setpt()):
        warn("Chamber pressure setpoint exceeds vapor pressure at shelf temperature " +\
        "setpoint(s). Drying cannot proceed.")
        return np.array([[0.0, Tsh_t(0), Tsh_t(0), Tsh_t(0), Pch_t(0), 0.0, 0.0]])

    inputs = (vial, product, ht, Pch_t, Tsh_t, dt, Lpr0)

    Lck0 = [0.0]
    T0 = Tsh_t(0)

    ################ Set up dynamic equation ######################
    # This function is defined here because it uses local variables, rather than
    # taking them as arguments.
    def calc_dLdt(t, u):
        # Time in hours
        Lck = u[0] # cm
        Tsh = Tsh_t(t)
        Pch = Pch_t(t)
        Kv = functions.Kv_FUN(ht['KC'],ht['KP'],ht['KD'],Pch)  # Vial heat transfer coefficient in cal/s/K/cm^2
        Rp = functions.Rp_FUN(Lck,product['R0'],product['A1'],product['A2'])  # Product resistance in cm^2-hr-Torr/g
        Tsub = fsolve(functions.T_sub_solver_FUN, T0, args = (Pch,vial['Av'],vial['Ap'],Kv,Lpr0,Lck,Rp,Tsh))[0] # Sublimation front temperature array in degC
        dmdt = functions.sub_rate(vial['Ap'],Rp,Tsub,Pch)   # Total sublimation rate array in kg/hr
        if dmdt<0:
            # print("Shelf temperature is too low for sublimation.")
            dmdt = 0.0
            dLdt = 0
            return [dLdt]
        # Tbot = functions.T_bot_FUN(Tsub,Lpr0,Lck,Pch,Rp)    # Vial bottom temperature array in degC

        dLdt = (dmdt*constant.kg_To_g)/(1-product['cSolid']*constant.rho_solution/constant.rho_solute)/(vial['Ap']*constant.rho_ice)*(1-product['cSolid']*(constant.rho_solution-constant.rho_ice)/constant.rho_solute) # cm/hr
        return [dLdt]

    ### ------ Condition for ending simulation: completed drying
    def finish(t, L):
        return Lpr0 - L[0]
    finish.terminal = True


    # ------- Solve the equations
    sol = solve_ivp(calc_dLdt, (0, max_t), Lck0, events=finish, 
                    vectorized=False, dense_output=True, method="BDF")
    if sol.t[-1] == max_t:# and Lpr0 > sol.y[0, -1]:
        warn("Maximum simulation time (specified by Pchamber and Tshelf) reached before drying completion.")

    output = functions.fill_output(sol, inputs)

    return output    

output_table = calc_knownRp.dry(vial,product,ht,Pchamber,Tshelf,dt)

It's a good idea, particularly when you are exploring interactively, to write the simulation input and output to disk together so that as you do later analysis, you have a record of what you did. Uncommenting the following code will do so.

Previously the inputs were recorded as a two-column CSV, which records info but isn't very readable (and frankly not very useful to a machine either). YAML provides a somewhat more sensible format, and fortunately a Python dictionary can be readily represented in YAML.

sim_setup = {
    'sim': sim,
    'vial': vial,
    'product': product,
    'ht': ht,
    'Pchamber': Pchamber,
    'Tshelf': Tshelf,
    'dt': dt,
}

# # Write input data to disk as YAML
# import time
# current_time = time.strftime("%Y%m%d_%H%M%S")
# try:
#     yamlfile = open('lyopronto_input_'+current_time+'.yaml', 'w')
#     yaml.dump(sim_setup, yamlfile)
# finally:
#     yamlfile.close()

# # Write simulation data to disk as CSV
# try:
#     csvfile = open('lyopronto_output_'+current_time+'.csv', 'w')
#     writer = csv.writer(csvfile)
#     writer.writerow(['Time [hr]','Sublimation Temperature [C]','Vial Bottom Temperature [C]', 'Shelf Temperature [C]','Chamber Pressure [mTorr]','Sublimation Flux [kg/hr/m^2]','Percent Dried'])
#     for i in range(0,len(output_table)):
#         writer.writerow(output_table[i])
# finally:
#     csvfile.close()

Finally, it's a good idea to plot everything. Here, again, you could save plots to disk by uncommenting the lines below.

Some default styling of the axes is carried out with lyopronto.plot_styling.axis_style_*, like turning on minor ticks and setting tick labels to be bold. However, this is entirely up to aesthetic taste and leaving it out does not affect the basic logic of the plots.

matplotlibrc('text.latex', preamble=r'\usepackage{color}')
matplotlibrc('text',usetex=False)
plt.rcParams['font.family'] = 'Arial'

figwidth = 30
figheight = 20
lineWidth = 5
markerSize = 20

fig = plt.figure(0,figsize=(figwidth,figheight))
ax1 = fig.add_subplot(1,1,1)
ax2 = ax1.twinx()
ax1.plot(output_table[:,0],output_table[:,4],'-o',color='b',markevery=5,linewidth=lineWidth, markersize=markerSize, label = "Chamber Pressure")
ax2.plot(output_table[:,0],output_table[:,5],'-',color=[0,0.7,0.3],linewidth=lineWidth, label = "Sublimation Flux")

plot_styling.axis_style_pressure(ax1)
plot_styling.axis_style_subflux(ax2)

plt.tight_layout()
# figure_name = 'lyopronto_pressure_subflux_'+current_time+'.pdf'
# plt.savefig(figure_name)
# plt.close()

fig = plt.figure(0,figsize=(figwidth,figheight))
ax = fig.add_subplot(1,1,1)
plot_styling.axis_style_percdried(ax)
ax.plot(output_table[:,0],output_table[:,-1],'-k',linewidth=lineWidth, label = "Percent Dried")
plt.tight_layout()
# figure_name = 'lyopronto_percentdried_'+current_time+'.pdf'
# plt.savefig(figure_name)
# plt.close()

fig = plt.figure(0,figsize=(figwidth,figheight))
ax = fig.add_subplot(1,1,1)
plot_styling.axis_style_temperature(ax)
ax.plot(output_table[:,0],output_table[:,1],'-b',linewidth=lineWidth, label = "Sublimation Front Temperature")
ax.plot(output_table[:,0],output_table[:,2],'-r',linewidth=lineWidth, label = "Maximum Product Temperature")
ax.plot(output_table[:,0],output_table[:,3],'-o',color='k',markevery=5,linewidth=lineWidth, markersize=markerSize, label = "Shelf Temperature")
plt.legend(fontsize=40,loc='best')
ll,ul = ax.get_ylim()
ax.set_ylim([ll,ul+5.0])
plt.tight_layout()
# figure_name = 'lyopronto_temperatures_'+current_time+'.pdf'
# plt.savefig(figure_name)
# plt.close()