Construct design space with known Kv and Rp

We need a few imports:

from ruamel.yaml import YAML
yaml = YAML()

from lyopronto import design_space, generate_visualizations

# Set up the simulation settings
# This needs to be a dict with string keys, which can be expressed in YAML as well
sim = {
    'tool':'Design Space Generator',
    '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': -15 # 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 evaluation points in Torr
    'setpt': [0.05, 0.1, 0.15],
}

Tshelf = {
    # init = Intial shelf temperature in C
    'init': -35.0,
    # ramp_rate = Shelf temperature ramping rate in C/min
    'ramp_rate': 1.0,
    # setpt = Shelf temperature set points in C
    'setpt': [-15, 0, 30, 90],
}

# Equipment capability (linear fit of sublimation rate vs chamber pressure)
eq_cap = {
    'a': -0.182,  # intercept [kg/hr]
    'b': 11.7     # slope [kg/hr/Torr]
}

# Number of vials
nVial = 2000

# Time step
dt = 0.01    # hr

Now, we are ready to construct the design space, which is lyopronto.design_space.dry.

Info

Here is the docstring for design_space.dry:

Compute quantities necessary for constructing a graphical design space.

Parameters:

Name	Type	Description	Default
vial	dict	Vial properties, including 'Vfill', 'Ap', and 'Av'..	required
product	dict	Product properties, including 'cSolid', 'T_pr_crit', and Rp parameters 'R0', 'A1', and 'A2'.	required
ht	dict	Heat transfer properties, including 'KC', 'KP', and 'KD'.	required
Pchamber	dict	Chamber pressure set points and time (see docs).	required
Tshelf	dict	Shelf temperature set points and time (see docs).	required
dt	float	Fixed time step for output [hours]	required
eq_cap	dict	Equipment capability line, with 'a' slope and 'b' intercept	required
nVial	int	Number of vials in the load, used for equipment capability calculation	required

Returns:

Type	Description
tuple[ndarray, ndarray, ndarray]	A tuple containing: - table of results for shelf isotherms - table of results for product isotherms - table of results for equipment capability curve

The first two returns have 5 rows corresponding to

• Maximum product temperature [degC]
• Primary drying time [hr]
• Average sublimation flux [kg/hr/m^2]
• Maximum/minimum sublimation flux [kg/hr/m^2]
• Sublimation flux at the end of primary drying [kg/hr/m^2]

The third return has 3 rows corresponding to the first three of that list.

With nT setpoints in Tshelf['setpt'] and nP setpoints in Pchamber['setpt'], the returned arrays have the following shapes: - Shelf isotherms: (5, nT, nP) array - Product isotherms: (5, 2) array (for the lowest and highest Pchamber setpoints) - Equipment capability curve: (3, nP) array

Source code in lyopronto/design_space.py

def dry(vial,product,ht,Pchamber,Tshelf,dt,eq_cap,nVial):
    """Compute quantities necessary for constructing a graphical design space. 

    Args:
        vial (dict): Vial properties, including 'Vfill', 'Ap', and 'Av'..
        product (dict): Product properties, including 'cSolid', 'T_pr_crit', and Rp parameters 'R0', 'A1', and 'A2'.
        ht (dict): Heat transfer properties, including 'KC', 'KP', and 'KD'.
        Pchamber (dict): Chamber pressure set points and time (see docs).
        Tshelf (dict): Shelf temperature set points and time (see docs).
        dt (float): Fixed time step for output [hours]
        eq_cap (dict): Equipment capability line, with 'a' slope and 'b' intercept
        nVial (int): Number of vials in the load, used for equipment capability calculation

    Returns:
        (tuple[ndarray, ndarray, ndarray]): A tuple containing:
            - table of results for shelf isotherms
            - table of results for product isotherms
            - table of results for equipment capability curve

    The first two returns have 5 rows corresponding to:
        - Maximum product temperature [degC]
        - Primary drying time [hr]
        - Average sublimation flux [kg/hr/m^2]
        - Maximum/minimum sublimation flux [kg/hr/m^2]
        - Sublimation flux at the end of primary drying [kg/hr/m^2]
    The third return has 3 rows corresponding to the first three of that list.

    With nT setpoints in Tshelf['setpt'] and nP setpoints in Pchamber['setpt'], the returned arrays have the following shapes:
        - Shelf isotherms: (5, nT, nP) array
        - Product isotherms: (5, 2) array (for the lowest and highest Pchamber setpoints)
        - Equipment capability curve: (3, nP) array

    """

    T_max = np.zeros([np.size(Tshelf['setpt']),np.size(Pchamber['setpt'])])
    drying_time = np.zeros([np.size(Tshelf['setpt']),np.size(Pchamber['setpt'])])
    sub_flux_avg = np.zeros([np.size(Tshelf['setpt']),np.size(Pchamber['setpt'])])
    sub_flux_max = np.zeros([np.size(Tshelf['setpt']),np.size(Pchamber['setpt'])])
    sub_flux_end = np.zeros([np.size(Tshelf['setpt']),np.size(Pchamber['setpt'])])

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

    ############  Shelf temperature isotherms ##########

    for i_Tsh,Tsh_setpt in enumerate(Tshelf['setpt']):

        for i_Pch,Pch in enumerate(Pchamber['setpt']):

            # Check for feasibility
            if functions.Vapor_pressure(Tsh_setpt) < Pch:
                # TODO: decide about how to gracefully exit
                # For now, just set outputs to NaN and continue
                warn(f"At Tshelf={Tsh_setpt} and Pch={Pch}, sublimation is not feasible (vapor pressure < chamber pressure).")
                T_max[i_Tsh,i_Pch] = np.nan
                drying_time[i_Tsh,i_Pch] = np.nan
                sub_flux_avg[i_Tsh,i_Pch] = np.nan
                sub_flux_max[i_Tsh,i_Pch] = np.nan
                sub_flux_end[i_Tsh,i_Pch] = np.nan
                continue

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

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

            # Initialization of cake length
            Lck = 0.0    # Cake length [cm]

            # Initial shelf temperature
            Tsh = Tshelf['init']        # [degC]
            # Time at which shelf temperature reaches set point [hr]
            t_setpt = abs(Tsh_setpt-Tshelf['init'])/Tshelf['ramp_rate']/constant.hr_To_min

            # Intial product temperature
            T0=Tsh   # [degC]

            # Vial heat transfer coefficient [cal/s/K/cm^2]
            Kv = functions.Kv_FUN(ht['KC'],ht['KP'],ht['KD'],Pch) 

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

            ################ Primary drying ######################

            while(Lck<=Lpr0): # Dry the entire frozen product

                Rp = functions.Rp_FUN(Lck,product['R0'],product['A1'],product['A2'])  # Product resistance [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 [degC]
                # Use previous value as guess for next iteration
                T0 = Tsub
                dmdt = functions.sub_rate(vial['Ap'],Rp,Tsub,Pch)   # Total sublimation rate [kg/hr]
                if dmdt<0:
                    warn(f"At t={t}hr, shelf temperature Tsh={Tsh} is too low for sublimation.")
                    dmdt = 0.0
                Tbot = functions.T_bot_FUN(Tsub,Lpr0,Lck,Pch,Rp)    # Vial bottom temperature [degC]

                # Sublimated ice length
                dL = (dmdt*constant.kg_To_g)*dt/(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]

            # Update record as functions of the cycle time
                if (iStep==0):
                    output_saved = np.array([[t, float(Tbot), dmdt/(vial['Ap']*constant.cm_To_m**2)]])
                else:
                    output_saved = np.append(output_saved, [[t, float(Tbot), dmdt/(vial['Ap']*constant.cm_To_m**2)]],axis=0)

                # Advance counters
                Lck_prev = Lck # Previous cake length [cm]
                Lck = Lck + dL # Cake length [cm]
                if (Lck_prev < Lpr0) and (Lck > Lpr0):
                    Lck = Lpr0    # Final cake length [cm]
                    dL = Lck - Lck_prev   # Cake length dried [cm]
                    t = iStep*dt + dL/((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)) # [hr]
                else:
                    t = (iStep+1) * dt # Time [hr]

                # Shelf temperature
                if t<t_setpt:
                    # Ramp till set point is reached
                    if Tshelf['init'] < Tsh_setpt:
                        Tsh = Tsh + Tshelf['ramp_rate']*constant.hr_To_min*dt
                    else:
                        Tsh = Tsh - Tshelf['ramp_rate']*constant.hr_To_min*dt
                else:
                    Tsh = Tsh_setpt    # Maintain at set point
                iStep = iStep + 1 # Time iteration number

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

            T_max[i_Tsh,i_Pch] = np.max(output_saved[:,1])    # Maximum product temperature [degC]

            drying_time[i_Tsh,i_Pch] = t    # Total drying time [hr]
            # TODO: consider whether to make this error rather than return NaN
            if T_max[i_Tsh,i_Pch] > 0: # exceeds melting temperature, not a physically valid solution
                warn(f"At Tsh={Tsh} and Pch={Pch}, computed temperatures of {T_max[i_Tsh,i_Pch]} exceed melting point of ice: check inputs.")
                sub_flux_avg[i_Tsh,i_Pch] = np.nan
                sub_flux_max[i_Tsh,i_Pch] = np.nan
                sub_flux_end[i_Tsh,i_Pch] = np.nan
                continue
            if output_saved.shape[0] <= 2:
                warn(f"At Tsh={Tsh} and Pch={Pch}, drying completed in single timestep: check inputs.")
                sub_flux_avg[i_Tsh,i_Pch] = np.nan
                sub_flux_max[i_Tsh,i_Pch] = np.nan
                sub_flux_end[i_Tsh,i_Pch] = np.nan
                continue
            del_t = output_saved[1:,0]-output_saved[:-1,0]
            del_t = np.append(del_t,del_t[-1])
            sub_flux_avg[i_Tsh,i_Pch] = np.sum(output_saved[:,2]*del_t)/np.sum(del_t)    # Average sublimation flux [kg/hr/m^2]
            sub_flux_max[i_Tsh,i_Pch] = np.max(output_saved[:,2])    # Maximum sublimation flux [kg/hr/m^2]
            sub_flux_end[i_Tsh,i_Pch] = output_saved[-1,2]    # Sublimation flux at end of primary drying [kg/hr/m^2]

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

    drying_time_pr = np.zeros([2])
    sub_flux_avg_pr = np.zeros([2])
    sub_flux_min_pr = np.zeros([2])
    sub_flux_end_pr = np.zeros([2])

    ############  Product temperature isotherms ##########

    for j,Pch in enumerate([Pchamber['setpt'][0],Pchamber['setpt'][-1]]):

        # Check for feasibility
        if functions.Vapor_pressure(product['T_pr_crit']) <= Pch:
            # TODO: decide about how to gracefully exit
            # For now, just set outputs to NaN and continue
            warn(f"With maximum T of Tcrit={product['T_pr_crit']} and Pch={Pch}, sublimation is not feasible (vapor pressure <= chamber pressure).")
            drying_time_pr[j] = np.nan
            sub_flux_avg_pr[j] = np.nan
            sub_flux_min_pr[j] = np.nan
            sub_flux_end_pr[j] = np.nan
            continue


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

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

        # Initialization of cake length
        Lck = 0.0    # Cake length [cm]

        # Vial heat transfer coefficient [cal/s/K/cm^2]
        Kv = functions.Kv_FUN(ht['KC'],ht['KP'],ht['KD'],Pch) 

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

        ################ Primary drying ######################

        while(Lck<=Lpr0): # Dry the entire frozen product

            Rp = functions.Rp_FUN(Lck,product['R0'],product['A1'],product['A2'])  # Product resistance [cm^2-hr-Torr/g]

            Tsub = fsolve(functions.T_sub_fromTpr, product['T_pr_crit'], args = (product['T_pr_crit'],Lpr0,Lck,Pch,Rp))[0] # Sublimation front temperature [degC]
            dmdt = functions.sub_rate(vial['Ap'],Rp,Tsub,Pch)   # Total sublimation rate [kg/hr]

            # Sublimated ice length
            dL = (dmdt*constant.kg_To_g)*dt/(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]

            # Update record as functions of the cycle time
            if (iStep==0):
                output_saved = np.array([[t, dmdt/(vial['Ap']*constant.cm_To_m**2)]])
            else:
                output_saved = np.append(output_saved, [[t, dmdt/(vial['Ap']*constant.cm_To_m**2)]],axis=0)

            # Advance counters
            Lck_prev = Lck # Previous cake length [cm]
            Lck = Lck + dL # Cake length [cm]
            if (Lck_prev < Lpr0) and (Lck > Lpr0):
                Lck = Lpr0    # Final cake length [cm]
                dL = Lck - Lck_prev   # Cake length dried [cm]
                t = iStep*dt + dL/((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)) # [hr]
            else:
                t = (iStep+1) * dt # Time [hr]
            iStep = iStep + 1 # Time iteration number

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

        drying_time_pr[j] = t    # Total drying time [hr]
        # TODO: consider whether this should error rather than return NaN
        if output_saved.shape[0] <= 2:
            warn(f"At Pch={Pch} and critical temp {product['T_pr_crit']}, drying completed in single timestep: check inputs.")
            sub_flux_avg_pr[j] = np.nan
            sub_flux_min_pr[j] = np.nan
            sub_flux_end_pr[j] = np.nan
            continue
        del_t = output_saved[1:,0]-output_saved[:-1,0]
        del_t = np.append(del_t,del_t[-1])
        sub_flux_avg_pr[j] = np.sum(output_saved[:,1]*del_t)/np.sum(del_t)    # Average sublimation flux [kg/hr/m^2]
        sub_flux_min_pr[j] = np.min(output_saved[:,1])    # Minimum sublimation flux [kg/hr/m^2]
        sub_flux_end_pr[j] = output_saved[-1,1]    # Sublimation flux at end of primary drying [kg/hr/m^2]

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

    T_max_eq_cap = np.zeros([np.size(Pchamber['setpt'])])

    ############  Equipment Capability ##########

    dmdt_eq_cap = eq_cap['a'] + eq_cap['b']*np.array(Pchamber['setpt'])    # Sublimation rate [kg/hr]
    if np.any(dmdt_eq_cap < 0):
        warn("Equipment capability sublimation rate is negative for some chamber pressures; setting to nan.")
        # dmdt_eq_cap = np.maximum(dmdt_eq_cap, 0.0)
        dmdt_eq_cap[dmdt_eq_cap <=0.0] = np.nan
    sub_flux_eq_cap = dmdt_eq_cap/nVial/(vial['Ap']*constant.cm_To_m**2)    # Sublimation flux [kg/hr/m^2]

    drying_time_eq_cap = Lpr0/((dmdt_eq_cap/nVial*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))    # Drying time [hr]

    Lck = np.linspace(0,Lpr0,100)    # Cake length [cm]
    Rp = functions.Rp_FUN(Lck,product['R0'],product['A1'],product['A2'])    # Product resistance [cm^2-hr-Torr/g]
    for k,Pch in enumerate(Pchamber['setpt']):
        T_max_eq_cap[k] = functions.Tbot_max_eq_cap(Pch,dmdt_eq_cap[k],Lpr0,Lck,Rp,vial['Ap'])        # Maximum product temperature [degC]

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

    return np.array([T_max,drying_time,sub_flux_avg,sub_flux_max,sub_flux_end]), \
        np.array([np.array([product['T_pr_crit'],product['T_pr_crit']]),drying_time_pr,sub_flux_avg_pr,sub_flux_min_pr,sub_flux_end_pr]), \
        np.array([T_max_eq_cap,drying_time_eq_cap,sub_flux_eq_cap])

# Run the design space simulation
ds_shelf, ds_pr, ds_eq_cap = design_space.dry(vial, product, ht, Pchamber, Tshelf, dt, eq_cap, nVial)

# 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.

sim_setup = {
    'sim': sim,
    'vial': vial,
    'product': product,
    'ht': ht,
    'Pchamber': Pchamber,
    'Tshelf': Tshelf,
    'eq_cap': eq_cap,
    'nVial': nVial,
    '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()

plots = generate_visualizations((ds_shelf, ds_pr, ds_eq_cap), sim_setup, "", save_figures=False) # Don't save to disk here, but you should in general