Imports

LyoPronto is this package. It reexports several other packages, so after using LyoPronto, you have effectively also done using Unitful and a few others.

using LyoPronto

These are other packages that I use in the test suite, but you can use others in their place. TypedTables provides a lightweight table structure, not as broadly flexible as a DataFrame but great for our needs

using TypedTables, CSV

Plots is a frontend for several plotting packages, and its companion package StatsPlots has a very nice macro I like.

using Plots
using StatsPlots: @df
using LaTeXStrings

For dealing with parameter structs and making copies, Accessors provides the @set and @reset macros

using Accessors

For helping make some nice plot labels

using Latexify: latexify, set_default
set_default(labelformat=:square) # Sets a Latexify default

Example Process Data

# Data start at 8th row of CSV file.
# This needs to point to the right file, which for documentation is kinda wonky
procdata_raw = CSV.read(joinpath(@__DIR__, "..", "..", "example", "2024-06-04-10_MFD_AH.csv"), Table, header=7)
t = uconvert.(u"hr", procdata_raw.CycleTime .- procdata_raw.CycleTime[1])
# At midnight, timestamps revert to zero, so catch that case
for i in eachindex(t)[begin+1:end]
    if t[i] < t[i-1]
        t[i:end] .+= 24u"hr"
    end
end
# Some of the dispatches don't like if time is not a float
t = float.(t)

# Rename the columns we will use, and add units
procdata = map(procdata_raw) do row
    # In the anonymous `do` function, `row` is a row of the table.
    # Return a new row as a NamedTuple
    (pirani = row.VacPirani * u"mTorr",
     cm = row.VacCPM * u"mTorr",
     T1 = row.TP1 * u"°C",
     T2 = row.TP2 * u"°C",
     T3 = row.TP4 * u"°C", # Quirk of this experimental run: TP3 slot was empty
     Tsh = row.ShelfSetPT * u"°C",
     phase = row.Phase # identify whether freezing, primary drying, or secondary
    )
end
procdata = Table(procdata, (;t)) # Append time to table

# Count time from the beginning of experiment
pd_data = filter(row->row.phase == 4, procdata)
tstart_pd = pd_data.t[1]
pd_data.t .-= pd_data.t[1]

t_end = identify_pd_end(pd_data.t, pd_data.pirani, Val(:onoff))

(18.66015837432897 hr, 20.604561173169635 hr)

Example models

Conventional lyophilization

# Vial geometry
# Ran with a 10mL vial, not strictly a 10R but with similar dimensions
Ap, Av = π.*get_vial_radii("10R") .^ 2

# Experimental conditions
T_shelf_0 = -40.0u"°C" |> u"K" # initial shelf temperature, in Kelvin for math reasons
T_shelf_final = -10.0u"°C" |> u"K" # final shelf temperature
ramp_rate = 0.5 *u"K/minute" # ramp rate
# Set points, followed, by ramp rate, followed by hold times if there are multiple ramps
Tsh = RampedVariable([T_shelf_0, T_shelf_final], ramp_rate)
# Single set point with no ramps
pch = RampedVariable(100u"mTorr")

# Formulation parameters
csolid = 0.05u"g/mL" # g solute / mL solution
ρsolution = 1u"g/mL" # g/mL total solution density
# Previously fitted values for Rp
R0 = 0.93u"cm^2*Torr*hr/g"
A1 = 21.1u"cm*Torr*hr/g"
A2 = 1.2u"1/cm"
Rp = RpFormFit(R0, A1, A2)
# Fit value for heat transfer coeff
Kshf = ConstPhysProp(13.9u"W/m^2/K")

# Fill
Vfill = 3u"mL"
hf0 = Vfill / Ap

po = ParamObjPikal([
    (Rp, hf0, csolid, ρsolution),
    (Kshf, Av, Ap),
    (pch, Tsh)
]);

prob = ODEProblem(po)
sol_conv = solve(prob, Rodas3());

System setpoints and conditions with `RampedVariable`s

Because it is very common to have a gradual ramp in temperature at the start of drying (and in general any time the set point changes), LyoPronto provides a tool for concisely describing the set point over time.

For a constant set point (with no ramps), provide a single value:

Tsh1 = RampedVariable(-15u"°C")

RampedVariable(-15 °C)

For a ramp from freezing temperature followed by a single hold (as is common in primary drying), provide the initial temperature, target temperature, and ramp rate:

# Convert from Celsius to Kelvin so that algebra can happen on backend
Tsh2 = RampedVariable([-40.0, -10]u"°C" .|> u"K", 1.0u"K/minute")

RampedVariable(Unitful.Quantity{Float64, 𝚯, Unitful.FreeUnits{(K,), 𝚯, nothing}}[233.14999999999998 K, 263.15 K], 1.0 K minute^-1)

For multiple ramps, provide set points, then ramp rates, then hold times in between ramps. In general, supply one less ramp than set points, and one less hold time than ramps.

Tsh3 = RampedVariable([-40, -20, 0]u"°C" .|> u"K", [2//3, 0.5]u"K/minute", [1u"hr"])

RampedVariable(Unitful.Quantity{Rational{Int64}, 𝚯, Unitful.FreeUnits{(K,), 𝚯, nothing}}[4663//20 K, 5063//20 K, 5463//20 K], Unitful.Quantity{Float64, 𝚯 𝐓^-1, Unitful.FreeUnits{(K, minute^-1), 𝚯 𝐓^-1, nothing}}[0.6666666666666666 K minute^-1, 0.5 K minute^-1], Unitful.Quantity{Int64, 𝐓, Unitful.FreeUnits{(hr,), 𝐓, nothing}}[1 hr])

A plot recipe is provided for all of these RampedVariables.

plot(xunit=u"hr", xlimit=(-0.1u"hr", 2.5u"hr"))
plot!(Tsh1, lw=3, label="No ramp")
plot!(Tsh2, lw=3, label="1 ramp")
plot!(Tsh3, lw=3, label="2 ramps")

Estimating Rp over time

The standard (Pikal) model for primary drying in lyophilization consists of, fundamentally, one ODE (change in frozen layer height) with a nonlinear algebraic constraint (pseudosteady heat and mass transfer, coupled at sublimation front). Since we have a system with one ODE and one nonlinear algebraic constraint, there is only one degree of freedom at a given time point. So with temperature measurements over time, we can compute corresponding mass flow over time or mass transfer resistance over length.

In LyoPronto, this functionality is implemented with the [calc_hRp_T](@ref) function, (think "compute $Rp(h_d)$ from $T_f(t)$").

To do so, we need to know about the experimental conditions; for that purpose, we pass a ParamObjPikal containing that information. To deal with the actual temperature series, we use a PrimaryDryFit object, which allows us to encode the way that, at some point, each temperature series deviates from the regular pseudosteady behavior governed by this model. (Strictly speaking, there are a variety of phenomena involved, but for here it is enough to say that at some point in time each temperature series experiences a sharp rise that is not described by the model.)

# The PrimaryDryFit:
fitdat_all = @df pd_data PrimaryDryFit(:t, (:T1[:t .< 15u"hr"],
                                    :T2[:t .< 13u"hr"],
                                    :T3[:t .< 16u"hr"]),)
plot(fitdat_all, nmarks=30, sampmarks=true)

Note that in this plot, T1 rises after 13 hours–I have deliberately included that to show what this will do in $R_p(h_d)$ space. Now, with the ParamObjPikal and PrimaryDryFit defined, we can calculate $R_p(h_d)$:

# Compute just for the first temperature series
calc_hRp_T(po, fitdat_all, i=1)
# Compute for all temperature series
hRps = [calc_hRp_T(po, fitdat_all; i) for i in 1:3]

3-element Vector{Tuple{Vector{Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(cm,), 𝐋, nothing}}}, Vector{Unitful.Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(g^-1, hr, cm^2, Torr), 𝐋 𝐓^-1, nothing}}}}}:
 ([0.0 cm, 0.00014725402256887766 cm, 0.0003092334473946433 cm, 0.0004859382744772783 cm, 0.0006773685038168201 cm, 0.0008816834601311172 cm, 0.0010988831434202123 cm, 0.0013308082289661702 cm, 0.0015756180414869302 cm, 0.0018333125809824396 cm  …  0.5641596072806395 cm, 0.5646234574517315 cm, 0.5650873076228234 cm, 0.5655511577939154 cm, 0.5660150079650074 cm, 0.5664788581360993 cm, 0.5669408676319092 cm, 0.5674010364524369 cm, 0.5678612052729646 cm, 0.5683213740934923 cm], [0.03210321556210788 hr cm^2 Torr g^-1, 0.12976891333942547 hr cm^2 Torr g^-1, 0.2112163499029719 hr cm^2 Torr g^-1, 0.2803447613231318 hr cm^2 Torr g^-1, 0.3398981320792421 hr cm^2 Torr g^-1, 0.5029892602064789 hr cm^2 Torr g^-1, 0.5430972484216493 hr cm^2 Torr g^-1, 0.5788441560484192 hr cm^2 Torr g^-1, 0.710497119553236 hr cm^2 Torr g^-1, 0.7352151600028615 hr cm^2 Torr g^-1  …  18.186756553943212 hr cm^2 Torr g^-1, 18.18762481078925 hr cm^2 Torr g^-1, 18.18849309915616 hr cm^2 Torr g^-1, 18.189361419045067 hr cm^2 Torr g^-1, 18.19022977045682 hr cm^2 Torr g^-1, 18.191098153392684 hr cm^2 Torr g^-1, 18.559143271472514 hr cm^2 Torr g^-1, 18.560012740436783 hr cm^2 Torr g^-1, 18.56088224043994 hr cm^2 Torr g^-1, 18.561751771482996 hr cm^2 Torr g^-1])
 ([0.0 cm, 0.0003386842519084204 cm, 0.0006902532307915796 cm, 0.0010528662613674415 cm, 0.0014283640189180404 cm, 0.0018185871787255673 cm, 0.0022216950655078722 cm, 0.0026358470039827987 cm, 0.003062883669432546 cm, 0.0035046457371391347 cm  …  0.5495309325098241 cm, 0.5501935756113842 cm, 0.5508562187129442 cm, 0.5515188618145043 cm, 0.5521815049160643 cm, 0.5528441480176244 cm, 0.5535067911191844 cm, 0.5541694342207445 cm, 0.5548431213739972 cm, 0.5554947204238645 cm], [0.0238139179112965 hr cm^2 Torr g^-1, 0.06798708683261871 hr cm^2 Torr g^-1, 0.16994261670529975 hr cm^2 Torr g^-1, 0.26773929838920146 hr cm^2 Torr g^-1, 0.3002635588185551 hr cm^2 Torr g^-1, 0.3307322812679421 hr cm^2 Torr g^-1, 0.4178681646279042 hr cm^2 Torr g^-1, 0.5022272473728764 hr cm^2 Torr g^-1, 0.5251117226031348 hr cm^2 Torr g^-1, 0.5468277482324843 hr cm^2 Torr g^-1  …  6.208043139712921 hr cm^2 Torr g^-1, 6.208789881749693 hr cm^2 Torr g^-1, 6.209536681739804 hr cm^2 Torr g^-1, 6.210283539687356 hr cm^2 Torr g^-1, 6.211030455596457 hr cm^2 Torr g^-1, 6.211777429471237 hr cm^2 Torr g^-1, 6.212524461315712 hr cm^2 Torr g^-1, 6.213271551134015 hr cm^2 Torr g^-1, 6.214031151884703 hr cm^2 Torr g^-1, 6.214765904708436 hr cm^2 Torr g^-1])
 ([0.0 cm, 0.00020063360575009632 cm, 0.0004159926137570574 cm, 0.0006442363487388188 cm, 0.0008853648106953309 cm, 0.0011393779996266458 cm, 0.0014062759155327087 cm, 0.0016860585584135774 cm, 0.0019787259282691917 cm, 0.0022824373498175027 cm  …  0.6149462925477891 cm, 0.6155353086380645 cm, 0.6161243247283399 cm, 0.6167133408186154 cm, 0.6173023569088908 cm, 0.6178913729991662 cm, 0.6184803890894416 cm, 0.6190694051797171 cm, 0.6196565805947104 cm, 0.6202419153344216 cm], [0.048170392546770024 hr cm^2 Torr g^-1, 0.018311855031507795 hr cm^2 Torr g^-1, 0.1837673932734037 hr cm^2 Torr g^-1, 0.2396936575284525 hr cm^2 Torr g^-1, 0.3820893215547277 hr cm^2 Torr g^-1, 0.42303350869457057 hr cm^2 Torr g^-1, 0.5483730700648033 hr cm^2 Torr g^-1, 0.5787104335807935 hr cm^2 Torr g^-1, 0.691059088155828 hr cm^2 Torr g^-1, 0.7978673964595586 hr cm^2 Torr g^-1  …  9.362312824148438 hr cm^2 Torr g^-1, 9.363121085137653 hr cm^2 Torr g^-1, 9.363929394811558 hr cm^2 Torr g^-1, 9.364737753172786 hr cm^2 Torr g^-1, 9.36554616022402 hr cm^2 Torr g^-1, 9.36635461596797 hr cm^2 Torr g^-1, 9.367163120407229 hr cm^2 Torr g^-1, 9.367971673544519 hr cm^2 Torr g^-1, 9.555029596335256 hr cm^2 Torr g^-1, 9.555840762732194 hr cm^2 Torr g^-1])

This returned a vector tuples, with a vector each for h_d(t) and R_p(t) for each set of temperature measurements. To plot this against a temperature fit as in the other tutorial here, we can do the following:

pl = plot(xlabel="h_d", ylabel="R_p", unitformat=latexify, xunit=u"cm", yunit=u"cm^2*Torr*hr/g")
for (i, hRp) in enumerate(hRps)
    plot!(hRp[1], hRp[2], label="T$i")
end
# For comparison, plot Rp as computed from fit in the other example
l = range(0u"cm", hf0, length=100)
plot!(l, Rp.(l), label="Direct fit to \$T(t)\$")