GitHub

Running a Simulation with LevelSetSublimation

Main Simulation Method

LevelSetSublimation.sim_from_dict — Function
sim_from_dict(config; tf=1e5, verbose=false)

Given a simulation configuration config, run a simulation.

Maximum simulation time (not CPU time, but simulation time) is specified by tf, which is in seconds; 1e5 is 270 hours. (No matter how much ice is left at that point, simulation will end.) verbose=true will put out some info messages about simulation progress, i.e. at each reinitialization.

All passed parameters should have Unitful unit annotations, so they can be in whatever units are convenient; inside this function, they will be converted to SI marks then have units stripped off (for easier numerical implementation). The results will be in SI units, so times in seconds, temperature in Kelvin, pressure in Pa, etc. Employ Unitful to do unit conversions after simulation if necessary.

config must have the following fields:

The following fields have default values and are therefore optional:

  • simgridsize, a tuple giving number of grid points to use for simulation. Defaults to (51, 51).
  • Tf0, an initial ice temperature. Defaults to Tsh(0) if not provided
  • Tvw0, an initial glass temperature. Defaults to Tf0.
  • time_integ: defaults to Val(:exp_newton); other options are Val(:dae), Val(:dae_then_exp) and Val(:ode_implicit). Among these three, slightly different problem formulations are used (explicit ODE with internal Newton solve, DAE, and implicit ODE).
  • init_prof, a Symbol indicating a starting profile (from make_ϕ0)

The three problem formulations have different advantages; in my testing, the DAE formulation and implicit formulation tend to run faster and more closely reflect the problem structure, but they run into instabilities when the sublimation front peels away from the wall, whereas the explicit ODE formulation can jump over that point.

source

Setting Up for Simulation

LevelSetSublimation.PhysicalProperties — Type

An immutable type for representing a slate of physical properties and constants needed.

A default constructor is provided which uses all defaults (borosilicate glass, water, ice), but thanks to Parameters.jl you can call with just values that need to be changed, e.g.

PhysicalProperties() # all defaults
PhysicalProperties(ρ_vw = 1u"kg/m^3", cp_vw = 1u"J/kg/K", εpp_vw = 1e-1) # adjust vial wall properties

All properties stored here are:

  • ρ_vw: Density of vial wall (defaults to borosilicate glass). Default: LevelSetSublimation.ρ_gl

  • cp_vw: Heat capacity of vial wall (defaults to borosilicate glass). Default: LevelSetSublimation.cp_gl

  • εpp_vw: Dielectric loss coefficient of vial wall (defaults to borosilicate glass). Default: LevelSetSublimation.εpp_gl

  • ρf: Density of frozen material (taken as water ice). Default: ρ_ice

  • Cpf: Heat capacity of frozen material (taken as water ice). Default: cp_ice

  • kf: Thermal conductivity of frozen material. Default: LevelSetSublimation.kf

  • Mw: Molecular weight of sublimating species (defaults to water). Default: 0.018 * u"kg/mol"

  • μ: Gas phase viscosity (in rarefied regime) of sublimating species (defaults to water). Default: LevelSetSublimation.μ

  • ΔH: Heat of sublimation of sublimating species (defaults to water); give as positive number. Default: LevelSetSublimation.ΔH

  • εppf: Dielectric loss coefficient of frozen layer (defaults to water ice). Default: LyoPronto.εppf

  • εpp_d: Dielectric loss coefficient of dry layer (defaults to 0). Default: 0.0

  • R: Universal gas constant. Kept as a constant so that it gets nondimensionalized with the rest. Default: u"R"

  • ϵ0: Vacuum permittivity. Kept as a constant so that it gets nondimensionalized with the rest. Default: u"ϵ0" |> u"F/m"

source
LevelSetSublimation.TimeConstantProperties — Type
TimeConstantProperties(ϵ, l, κ, Rp0, kd, Kvwf, m_v, A_v, B_d, B_f, B_vw)

A struct for holding physical properties which are likely to change from case to case.

No default constructor is provided by intention–all of these parameters should be at least considered before running a simulation.

  • ϵ

  • l

  • κ

  • Rp0

  • kd

  • Kvwf

  • m_v

  • A_v

  • B_d

  • B_f

  • B_vw

source
LevelSetSublimation.TimeVaryingProperties — Type
TimeVaryingProperties(f_RF, P_per_vial, Tsh, pch, Kshf)

A struct for holding controlled parameters that may change over time.

To get the value of all those parameters at a given time, call the struct with a time, and it will evaluate each field at that time and provide a TimeVaryingPropertiesSnapshot.

No default constructor is provided by intention–all of these parameters should be at least considered before running a simulation.

Each should be passed as a callable, which returns the value of the parameter as a function of time, with the exception of Kshf which is a function of pressure. See the RampedVariable and RpFormFit types from LyoPronto, which are intended to make this more convenient.

  • f_RF

  • P_per_vial

  • Tsh

  • pch

  • Kshf

source

Here is a sample simulation setup, in a case where nothing complicated is happening.

# ---- Properties which do not change in time

# Mass transfer
ϵ = 0.95 # 90% porosity
κ = 0.0u"m^2" # ~size^2 of a pore
Rp0 = 1.4u"cm^2*hr*Torr/g"
A1 = 16u"cm*hr*Torr/g"
Tguess = 260u"K"
l = sqrt(base_props.R*Tguess/base_props.Mw) / A1
# Heat transfer
kd = LSS.k_sucrose * (1-ϵ)
m_v = LyoPronto.get_vial_mass(vialsize)
A_v = π*LyoPronto.get_vial_radii(vialsize)[2]^2
# Microwave
B_d = 0.0u"Ω/m^2"
tcprops = TimeConstantProperties(ϵ, l_bulk, κ, Rp0, kd, Kvwf, m_v, A_v, B_d, Bf, Bvw)

# Properties which may change in time
f_RF = RampedVariable(8.0u"GHz")
pch = RampedVariable(100.0u"mTorr")
Tsh = RampedVariable(uconvert.(u"K", [-40.0, 10]*u"°C"), 1u"K/minute")
P_per_vial = RampedVariable(10u"W"/17)
# Heat transfer coefficient as function of pressure
KC = 2.75e-4u"cal/s/K/cm^2"
KP = 8.93e-4u"cal/s/K/cm^2/Torr"
KD = 0.46u"1/Torr"
Kshf = RpFormFit(KC, KP, KD)
tvprops = TimeVaryingProperties(f_RF, P_per_vial, Tsh, pch, Kshf)

paramsd = LSS.base_props, tcprops, tvprops

vialsize = "6R"
fillvol = 5u"mL"
simgridsize = (41, 31)

config = Dict{Symbol, Any}()
@pack! config = paramsd, vialsize, fillvol, simgridsize

Run a simulation

res = sim_from_dict(config)

Visualize and work with simulation results

The Guts (not exported)

Heat-transfer-only simulation: OUTDATED

LevelSetSublimation.sim_heatonly — Function
sim_heatonly(config; tf=1e5, verbose=false)

Given a simulation configuration config, run a simulation.

This simulation is stripped-down: no mass transfer, no variation in ice & glass temperature

Maximum simulation time is specified by tf. verbose=true will put out some info messages about simulation progress, i.e. at each reinitialization.

CURRENTLY OUT OF DATE with parameter input structures.

source
LevelSetSublimation.dudt_heatonly — Function
dudt_heatonly(u, dom::Domain, params)
dudt_heatonly(u, config)

Compute the time derivative of u with given parameters.

Wraps a call on dudt_heatonly!, for convenience in debugging and elsewhere that efficiency is less important

source
LevelSetSublimation.dudt_heatonly! — Function
dudt_heatonly!(du, u, p, t)

Evaluate local time rate of change for u and put results in du.

This function leaves Tf and Tvw untouched, since there isn't a way to govern their dynamics without mass transfer.

u and du are both structured according to

This is a right-hand-side for ∂ₜϕ = -v⋅∇ϕ, where v = (vr, vz) is evaluated by computing and extrapolating front velocity using compute_frontvel_mass.

source