@@ -15,30 +15,31 @@ const end_drying_callback = ContinuousCallback(end_cond, terminate!, save_positi
1515# -------------------------------------------
1616# Incorporate the nonlinear algebraic part in a DAE formulation.
1717# This has the advantage that, afterward, temperatures can be cheaply interpolated by builtin solutions
18+ """
19+ $(SIGNATURES)
1820
19- @doc raw """
20- lyo_1d_dae!(du, u, params, t)
21+ With the Pikal model, compute the mass flow at `u=[hf, Tf]`, time `t`, and conditions `po`.
22+ Ideally `po` should be a `ParamObjPikal`.
2123
22- Internal implementation of the Pikal model.
23- See [`lyo_1d_dae_f`](@ref) for the wrapped version, which is more fully documented .
24+ This gets used inside the DAE solve, but for consistency is provided separately so to avoid
25+ the risk of making unit mistakes .
2426"""
25- function lyo_1d_dae! (du, u, params , t)
27+ function calc_md_Q ( u, po , t)
2628
27- if params isa ParamObj
29+ if po isa ParamObj
2830 (; Rp, hf0, csolid, ρsolution,
29- Kshf, Av, Ap, pch, Tsh) = params
31+ Kshf, Av, Ap, pch, Tsh) = po
3032 else
31- Rp, hf0, csolid, ρsolution = params [1 ]
32- Kshf, Av, Ap = params [2 ]
33- pch, Tsh = params [3 ]
33+ Rp, hf0, csolid, ρsolution = po [1 ]
34+ Kshf, Av, Ap = po [2 ]
35+ pch, Tsh = po [3 ]
3436 end
35-
37+
3638 td = t* u " hr" # Dimensional time
3739 hf = u[1 ]* u " cm"
3840 Tf = u[2 ]* u " K"
39- if Tf < 0.0 u " K"
40- du .= NaN
41- return nothing
41+ if Tf < 0.0 u " K" # Unphysical temperature needs an escape hatch
42+ return NaN * u " kg/s" , NaN * u " W"
4243 end
4344
4445 hd = hf0 - hf
@@ -47,6 +48,31 @@ function lyo_1d_dae!(du, u, params, t)
4748 Tsub = Tf - Qshf/ k_ice/ Ap* hf
4849 delta_p = calc_psub (Tsub)- pch (td)
4950 dmdt = - Ap* (delta_p)/ Rp (hd)
51+ return dmdt, Qshf
52+ end
53+
54+ @doc raw """
55+ lyo_1d_dae!(du, u, params, t)
56+
57+ Internal implementation of the Pikal model.
58+ See [`lyo_1d_dae_f`](@ref) for the wrapped version, which is more fully documented.
59+ """
60+ function lyo_1d_dae! (du, u, params, t)
61+
62+ # Need a handful of parameters in this function.
63+ if params isa ParamObj
64+ (; csolid, ρsolution, Ap) = params
65+ else
66+ csolid, ρsolution = params[1 ][3 : 4 ]
67+ Ap = params[2 ][3 ]
68+ end
69+ # This logic is carried out in a separate function,
70+ # so that it can be reused after the fact for computing mass flow.
71+ dmdt, Qshf = calc_md_Q (u, params, t)
72+ if isnan (dmdt)
73+ du .= NaN
74+ return nothing
75+ end
5076 Qsub = uconvert (u " W" , dmdt* ΔHsub)
5177
5278 dhf_dt = min (0.0 u " cm/hr" , dmdt/ (ρsolution- csolid)/ Ap |> u " cm/hr" ) # Cap dhf_dt at 0: no desublimation
@@ -214,7 +240,6 @@ function dae_Rp!(du, u, p, tn)
214240
215241 Tf = Tf_interp (t)
216242
217-
218243 Q = Kshf (pch (t))* Av* (Tsh (t) - Tf)
219244 Tsub = Tf - Q/ Ap/ LyoPronto. k_ice * (hf0- hd)
220245 md = Q/ LyoPronto. ΔH
0 commit comments