All Docstrings

PrimaryDryFit: a type for storing experimental data and indicating how it should be fit.

Provided constructors:

PrimaryDryFit(t, Tfs, Tvw, t_end)
PrimaryDryFit(t, Tfs) = PrimaryDryFit(t, Tfs, missing, missing)
PrimaryDryFit(t, Tfs, Tvw) = PrimaryDryFit(t, Tfs, Tvw, missing)
PrimaryDryFit(t, Tfs, t_end::Union{Unitful.Time}, Tuple{Tt, Tt}}) where Tt  = PrimaryDryFit(t, Tfs, missing, t_end)

The use of this struct is determined in large part by the implementation of LyoPronto.obj_expT. If a given field is not available, set it to missing and things should basically work. At least t and Tfs are expected to always be provided.

In the end, Tfs and Tvws are each stored as a tuple of vectors, but the constructor tries to be flexible about allowing a single vector to be passed in place of a tuple of vectors.

Tf_iend and Tvw_iend default to [length(Tf) for Tf in Tfs] and [length(Tvw) for Tvw in Tvws], respectively, with one value for each temperature series; they are used to dictate if a given temperature series should be truncated sooner than the full length in the fitting procedure. This implies that all the temperature series are valid initially at the same time points, then stop having measured values after a different number of measurements. If a single value is given for Tvw, then it is taken to be an endpoint, and Tvw_iend will be missing.

t_end indicates an end of drying, particularly if taken from other measurements (e.g. from Pirani-CM convergence). If set to missing, it is ignored in the objective function. If set to a tuple of two times, then in the objective function any time in that window is not penalized; outside that window, squared error takes over, as for the single time case.

Principal Cases:

• Conventional: provide only t, Tfs
• Conventional with Pirani ending: provide t, Tfs, t_end
• RF with measured vial wall: provide t, Tfs, Tvws,
• RF, matching model Tvw to experimental Tf[end] without measured vial wall: provide t, Tfs, Tvw

A convenience type for computing temperatures, pressures, etc. with multiple setpoints in sequence, and linear interpolation according to a fixed ramp rate between set points

Three main constructors are available: For a non-varying value, call with one argument:

RampedVariable(constant_setpt)

For one ramp from initial value to set point with indefinite hold, call with two arguments:

RampedVariable(setpts, ramprate)

And for multiple setpoints, call with three arguments:

RampedVariable(setpts, ramprates, holds)

With three arguments, setpts, ramprates, and holds should all be vectors, with lengths N+1, N, N-1 respectively.

The resulting RampedVariable rv = RampedVariable(...) can be called as rv(x) at any (dimensionally consistent) value of x, and will return the value at that time point along the ramp process.

A plot recipe is also provided for this type, e.g. plot(rv; tmax=10u"hr") where tmax indicates where to stop drawing the last setpoint hold.

A convenience type for dealing with the common functional form given to Rp and Kv.

An object Rp = RpFormFit(A, B, C) can be called as Rp(x), which simply computes A + B*x/(1 + C*x). Likewise, Kv = RpFormFit(Kc, Kp, Kd) can be called as Kv(p) to get Kc + Kp*p/(1 + Kd*p).

Be careful to pass dimensionally consistent values.

KBB_transform_basic(Kvwfg, Bfg, Bvwg)

Construct a typical transform for fitting Kvwf, Bf, and Bvw (as for a microwave cycle).

KRp_transform_basic(Kshfg, R0g, A1g, A2g)

Construct a typical transform for fitting both Kshf and Rp.

K_transform_basic(Kshfg)

Construct a typical transform for fitting Kshf (a.k.a. Kv).

Rp_transform_basic(R0g, A1g, A2g)

Construct a typical transform for fitting Rp.

err_expT(errs, sol, pdfit; tweight, verbose)

Evaluate the error between model solution sol to experimental data in pdfit.

In contrast to obj_expT(), this function makes an array of all the errors, which would be squared and summed to produce an objective function.

• errs is a vector of length num_errs(pdfit), which this function fills with the errors.

-sol is a solution to an appropriate model; see gen_sol_pd for a helper function.

• pdfit is an instance of PrimaryDryFit, which contains some information about what to compare.
• tweight = 1 gives the weighting (in K^2/hr^2) of the total drying time in the objective, as compared to the temperature error.

Each time series, plus the end time, is given equal weight by dividing by its length; error is given in K (but ustripped).

Note that if pdfit has vial wall temperatures (i.e. ismissing(pdfit.Tvws) == false), the third-index variable in sol is assumed to be temperature, as is true for solutions with ParamObjRF.

If there are multiple series of Tf in pdfit, error is computed for each separately; likewise for Tvw.

gen_nsol_pd(fitlog, tr, pos, fitdats; badprms, kwargs...)

Generate multiple solutions at once.

If the transformation tr makes something (e.g. NamedTuple) with properties separate and shared, then one each of separate is combined with shared, then they are matched up with each element of pos and fitdats. If pos is a single object, it is repeated. Further, if tr also has a field sep_inds, then those indices are used to map separate to the sets of pos and fitdats. This is useful for when e.g. 2 sets of separate parameters are to be applied across 5 different experiments.

gen_sol_pd(fitlog, tr, po, fitdat; badprms, kwargs...)

gen_sol_pd(fitlog, tr, po; saveat, badprms, kwargs...)

Simulate primary drying, given a vector of parameter guesses, a mapping tr from fitlog to named coefficients, and other parameters in po.

The equations used are determined by the type of po, which (with the magic of dispatch) is used to set up an ODE system.

tr should be a TransformTuple object, from TransformVariables, which maps e.g. a vector of 3 real numbers to a NamedTuple with R0, A1, A2 as keys and appropriate Unitful dimensions on the values. This small function runs

fitprm = transform(tr, fitlog)
new_params = setproperties(po, fitprm)
prob = ODEProblem(new_params; tspan=(0.0, 1000.0))
sol = solve(prob, Rodas4(autodiff=AutoForwardDiff(chunksize=2)); saveat, kwargs...)

which is wrapped to avoid code duplication.

So, to choose which parameters to fitting, all that is necessary is to provide an appropriate transform tr and add a method of setproperties for the desired parameters. Therefore this function can be used for both K-Rp fitting, or just Rp, or just a subset of the 3 Rp coefficients.

If given, fitdat is used to set saveat for the ODE solution.

Other kwargs are passed directly (as is) to the ODE solve call.

nls_pd(errs, fitlog, tpf; tweight, verbose)

Calculate the sum of squared error (objective function) for fitting parameters to primary drying data. This directly calls gen_sol_pd, then obj_expT, so see those docstrings.

obj_expT(sol, pdfit; tweight, verbose, Tvw_weight)

Evaluate an objective function which compares model solution computed by sol to experimental data in pdfit.

• sol is a solution to an appropriate model; see gen_sol_pd for a helper.
• pdfit is an instance of PrimaryDryFit, which contains some information about what to compare.
• tweight = 1 gives the weighting (in K^2/hr^2) of the total drying time in the objective, as compared to the temperature error.
• Tvw_weight = 1 gives the weighting of Tvw in the objective, as compared to Tf.

Note that if pdfit has vial wall temperatures (i.e. ismissing(pdfit.Tvws) == false), the third-index variable in sol is assumed to be temperature, as is true for the lumped capacitance model (see ParamObjRF`.

If there are multiple series of Tf in pdfit, squared error is computed for each separately then summed; likewise for Tvw.

I've considered writing several methods and dispatching on pdfit somehow, which would be cool and might individually be easier to read. But control flow might be harder to document and explain, and this should work just fine.

obj_pd(fitlog, tpf; tweight, badprms, verbose)

Calculate the sum of squared error (objective function) for fitting parameters to primary drying data. This directly calls gen_sol_pd, then obj_expT, so see those docstrings.

exptfplot(time, T1, [T2, ...]; labsuffix=", exp.")
exptfplot!(time, T1, [T2, ...]; labsuffix=", exp.")

Plot recipe for one or more experimentally measured product temperatures, all at same times. This recipe adds one series for each passed temperature series, with labels defaulting to "T_{fi}"*labsuffix.

exptfplot(time, T1, [T2, ...]; labsuffix=", exp.")
exptfplot!(time, T1, [T2, ...]; labsuffix=", exp.")

Plot recipe for one or more experimentally measured product temperatures, all at same times. This recipe adds one series for each passed temperature series, with labels defaulting to "T_{fi}"*labsuffix.

exptvwplot(time, T1, [T2, ...]; trim)
exptvwplot!(time, T1, [T2, ...]; trim)

Plot recipe for a set of experimentally measured vial wall temperatures. This recipe adds one series for each passed temperature series, with labels defaulting to "T_{vwi}"*labsuffix. trim is an integer, indicating how many points to skip at a time, so that the dotted line looks dotted even with noisy data.

exptvwplot(time, T1, [T2, ...]; trim)
exptvwplot!(time, T1, [T2, ...]; trim)

Plot recipe for a set of experimentally measured vial wall temperatures. This recipe adds one series for each passed temperature series, with labels defaulting to "T_{vwi}"*labsuffix. trim is an integer, indicating how many points to skip at a time, so that the dotted line looks dotted even with noisy data.

modconvtplot(sols; sampmarks=false, labsuffix = ", model")
modconvtplot!(sols; sampmarks=false, labsuffix = ", model")

Plot recipe for one or multiple solutions to the Pikal model, e.g. the output of gen_sol_pd. This adds a series to the plot for each passed solution, with labels defaulting to "T_{fi}"*labsuffix.

If sampmarks is true, the solution will be interpolated to evenly spaced time points and markers will be added to some of those points.

modconvtplot(sols; sampmarks=false, labsuffix = ", model")
modconvtplot!(sols; sampmarks=false, labsuffix = ", model")

Plot recipe for one or multiple solutions to the Pikal model, e.g. the output of gen_sol_pd. This adds a series to the plot for each passed solution, with labels defaulting to "T_{fi}"*labsuffix.

If sampmarks is true, the solution will be interpolated to evenly spaced time points and markers will be added to some of those points.

modrftplot(sol, labsuffix=", model", sampmarks=false, trimend=0)
modrftplot!(sol, labsuffix=", model", sampmarks=false, trimend=0)

Plot recipe for one solution to the lumped capacitance model.

This adds two series to the plot, with labels defaulting to ["T_f" "T_{vw}"] .* labsuffix. The optional argument trimend controls how many time points to trim from the end (which is helpful if temperature shoots up as mf -> 0).

If sampmarks is true, the solution will be interpolated to evenly spaced time points and markers will be added to some of those points.

Since this is a recipe, any Plots.jl keyword arguments can be passed to modify the plot.

modrftplot(sol, labsuffix=", model", sampmarks=false, trimend=0)
modrftplot!(sol, labsuffix=", model", sampmarks=false, trimend=0)

Plot recipe for one solution to the lumped capacitance model.

This adds two series to the plot, with labels defaulting to ["T_f" "T_{vw}"] .* labsuffix. The optional argument trimend controls how many time points to trim from the end (which is helpful if temperature shoots up as mf -> 0).

If sampmarks is true, the solution will be interpolated to evenly spaced time points and markers will be added to some of those points.

Since this is a recipe, any Plots.jl keyword arguments can be passed to modify the plot.

qrf_integrate(sol, RF_params)

Compute the integral over time of each heat transfer mode in the lumped capacitance model.

Returns as a Dict{String, Quantity{...}}, with string keys Qsub, Qshf, Qvwf, QRFf, QRFvw.

tendplot(t_end)
tendplot!(t_end)
tendplot(t_end1, t_end2)
tendplot!(t_end1, t_end2)

Plot recipe that adds a labeled vertical line to the plot at time t_end. A default label and styling are applied, but these can be modified by keyword arguments as usual If two time points are passed, a light shading is applied between instead of a vertical line.

tendplot(t_end)
tendplot!(t_end)
tendplot(t_end1, t_end2)
tendplot!(t_end1, t_end2)

Plot recipe that adds a labeled vertical line to the plot at time t_end. A default label and styling are applied, but these can be modified by keyword arguments as usual If two time points are passed, a light shading is applied between instead of a vertical line.

get_vial_mass(vialsize::String)

Return vial mass for given ISO vial size.

Uses a table provided by Schott, stored internally in a CSV.

get_vial_radii(vialsize::String)

Return inner and outer radius for passed ISO vial size.

Uses a table provided by Schott, stored internally in a CSV.

get_vial_shape(vialsize::String)

Return a Dict{Symbol, } with a slew of vial dimensions, useful for drawing the shape of the vial with make_outlines.

get_vial_thickness(vialsize::String)

Return vial wall thickness for given ISO vial size.

Uses a table provided by Schott, stored internally in a CSV.

make_outlines(dims, Vfill)

Return a sequence of points (ready to be made into Plots.Shapes for the vial and fill volume, with Unitful dimensions, for given vial dimensions.

This is a convenience function for making figures illustrating fill depth.

A callback for use in simulating either the Pikal or RF model.

Terminates the time integration when end_cond evaluates to true.

struct ParamObjPikal{__T_Rp, __T_hf0, __T_csolid, __T_ρsolution, __T_Kshf, __T_Av, __T_Ap, __T_pch, __T_Tsh} <: LyoPronto.ParamObj

The ParamObjPikal type is a container for the parameters used in the Pikal model.

struct ParamObjRF{__T_Rp, __T_hf0, __T_csolid, __T_ρsolution, __T_Kshf, __T_Av, __T_Ap, __T_pch, __T_Tsh, __T_P_per_vial, __T_mf0, __T_cpf, __T_mv, __T_cpv, __T_Arad, __T_f_RF, __T_eppf, __T_eppvw, __T_Kvwf, __T_Bf, __T_Bvw, __T_alpha} <: LyoPronto.ParamObj

The ParamObjRF type is a container for the parameters used in the RF model.

end_cond(u, t, integ)

Compute the end condition for primary drying (that mf or hf approaches zero).

lumped_cap_rf!(du, u, params, tn)
lumped_cap_rf!(du, u, params, tn, qret)

Compute the right-hand-side function for the ODEs making up the lumped-capacitance microwave-assisted model.

The optional argument qret defaults to Val(false); if set to Val(true), the function returns [Q_sub, Q_shf, Q_vwf, Q_RF_f, Q_RF_vw, Q_shw] with Q_... as Unitful quantities in watts. The extra results are helpful in investigating the significance of the various heat transfer modes, but are not necessary in the ODE integration.

du refers to [dmf/dt, dTf/dt, dTvw/dt], with u = [mf, Tf, Tvw]. u is taken without units but assumed to have the units of [g, K, K] (which is internally added). tn is assumed to be in hours (internally added), so dudt is returned with assumed units [g/hr, K/hr, K/hr] to be consistent.

It is recommended to use the ParamObjRF type to hold the parameters, since it allows some more convenient access to the parameters, but they can be given in the form of a tuple-of-tuples:

params = (   
    (Rp, hf0, cSolid, ρsolution),
    (Kshf_f, Av, Ap),
    (pch, Tsh, P_per_vial),
    (mf0, cpf, mv, cpv, Arad),
    (f_RF, eppf, eppvw),
    (Kvwf, Bf, Bvw, alpha),
)

This tuple-of-tuples structure is also used in an extra constructor for the ParamObjRF type.

These should all be Unitful quantities with appropriate dimensions, with some exceptions which are callables returning quantities. See RpFormFit and RampedVariable for convenience types that can help with these cases.

• Rp(x) with x a length returns mass transfer resistance (as a Unitful quantity)

• Kshf_f(p) with p a pressure returns heat transfer coefficient (as a Unitful quantity).

• Tsh(t), pch(t), P_per_vial(t) return shelf temperature, chamber pressure, and microwave power respectively at time t.

• Arad and alpha are used only in prior versions of the model, and can be left out.

This is my updated version of the model. LC3: Q_shw evaluated with Kshf; shape factor included; α=0

lyo_1d_dae!(du, u, params, t)

Internal implementation of the Pikal model. See lyo_1d_dae_f for the wrapped version, which is more fully documented.

lyo_1d_dae_f = ODEFunction(lyo_1d_dae!, mass_matrix=lyo_1d_mm)

Compute the right hand side function for the Pikal model.

The DAE system which is the Pikal model (1 ODE, one nonlinear algebraic equation for pseudosteady conditions) is here treated as a constant-mass-matrix implicit ODE system. The implementation is in lyo_1d_dae!

The initial conditions u0 = [h_f, Tf] should be unitless, but are internally assigned to be in [cm, K]. The unitless time is taken to be in hours, so derivatives are given in unitless [cm/hr, K/hr].

params is a ParamObjPikal, which can be constructed with the following form (helping with readability):

params = ParamObjPikal((
    (Rp, hf0, csolid, ρsolution),
    (Kshf, Av, Ap),
    (pch, Tsh) ,
))

where those listed following are callables returning Quantitys, and the rest are Quantitys. See RpFormFit and RampedVariable for convenience types that can help with the callables.

• Rp(x) with x a length returns mass transfer resistance (as a Unitful quantity)
• Kshf(p) with p a pressure returns heat transfer coefficient (as a Unitful quantity).
• Tsh(t), pch(t) return shelf temperature and chamber pressure respectively at time t.

Molecular weight of water

Glass heat capacity Rough estimate from Bansal and Doremus 1985, "Handbook of Glass Properties"

Ice heat capacity IAPWS for ice: use triple point value for simplicity

Permittivity of free space

Glass thermal conductivity From Bansal and Doremus 1985, rough estimate based an a sodium borosilicate glass

Thermal conductivity of ice From Slack, 1980 200K : .032 W/cmK 250K : .024 W/cmK 273K : .0214 W/cmK

Thermal conductivity of sucrose, Caster grade powder From McCarthy and Fabre, 1989 book chapter

Density of ice at triple point.

Heat of sublimation of ice Approximate: coresponds to 254K or 225K From Feistel and Wagner, 2006

Dielectric loss coefficient of borosilicate glass, per Schott's testing

Water vapor viscosity in dilute limit

From Hellmann and Vogel, 2015 250K: 8.054 μPas 260K: 8.383 μPas 270K: 8.714 μPa*s

Glass density SCHOTT measurement, posted online

Density of sucrose, Caster grade powder From McCarthy and Fabre, 1989 book chapter

Stefan-Boltzmann Constant

Using the same Arrhenius fit as calc_psub, compute the sublimation temperature at a given pressure.

calc_psub(T) 
calc_psub(T::Q) where Q<:Quantity

Compute pressure (in Pascals) of sublimation at temperature T in Kelvin.

This is essentially an Arrhenius fit, where we compute: psub = pref * exp(-ΔHsub*Mw/R/T)

Single-purpose module for computing the dielectric loss coefficient of ice, as a function of temperature and frequency. This module exports a function ϵppf(T, f) for doing so. (Set as a separate module just to keep namespaces clean.) This code is a nearly-direct implementation of the correlation, Eqs. 3-6 presented in: Takeshi Matsuoka, Shuji Fujita, Shinji Mae; Effect of temperature on dielectric properties of ice in the range 5–39 GHz. J. Appl. Phys. 15 November 1996; 80 (10): 5884–5890. https://doi.org/10.1063/1.363582

ϵppf(T, f)

Compute the dielectric loss of ice as a function of temperature and frequency. Expects temperature in Kelvin and frequency in Hz, with Unitful units.