Solving steady-state heat equation

solve_T(u, Tf, dom::Domain, params)

Compute 2D axisymmetric T profile, returning ghost cell values, for given state u and Tf.

Neumann boundary conditions on all rectangular boundaries; Dirichlet on zero-level set.

This implementation uses second-order finite differences, with linear extrapolation into Ω⁻. Coefficients are all hard-coded here, unfortunately. (For more on extrapolation, see Gibou et al., 2002, "Second-Order-Accurate ... Poisson ... ") Neumann boundaries use a ghost point & BC to define ghost cell, then use same stencil as normal. Coefficients computed in gfm_extrap.ipynb, using Sympy. (For higher order, see Gibou and Fedkiw, 2005, "A fourth order accurate discretization ... Laplace ... ")

solve_p(u, T, dom::Domain, params[, p0]; maxit=20, reltol=1e-6) where G<:AbstractArray

Iteratively compute the pressure profile for system state u, Tf, and T.

There is a weak nonlinearity in the system, since the mass transfer coefficient b depends partially on pressure. To treat this, use a guessed pressure p0 (which, if not provided, is set everywhere to chamber pressure) to compute b, then perform a linear solve for p using solve_p_given_b. At this point, recompute b, then recompute p.

In preliminary testing, this usually converges within 5 or 10 iterations. Usually if it doesn't converge, it is because temperatures are outside the expected range, yielding crazy sublimation pressures. (Occasionally it means I incorrectly wrote a finite difference somewhere.)

solve_p_given_b(ϕ, b, Tf, dom::Domain, params)

Compute 2D axisymmetric pseudosteady pressure profile for given values of level set function ϕ, temperature T, and transport coefficient b.

b is a dusty-gas transport coefficient for the pressure, which can vary spatially. Homogeneous Neumann boundary conditions at r=0, r=R, z=0; Dirichlet on zero-level set (p=p_sub), Robin at top (dp/dz = (pch-p)/Rp0). params should have fields:

• Rp0 : zero-thickness resistance offset, often written R0 in lyo literature
• pch : chamber (or vial) pressure at top surface

This implementation uses second-order finite differences, with linear extrapolation into Ω⁻. (For details, see (Gibou and Fedkiw, 2005).) Coefficients are all hard-coded here, unfortunately. Neumann boundaries use a ghost point & BC to define ghost cell, then use same stencil as normal. Coefficients computed in gfm_extrap.ipynb, using Sympy.

function compute_Qice(ϕ, T, p, dom::Domain, params)

Compute the total heat input into frozen & dried domains. Also passes glass-ice heat as separate return.

See p. 112-113 from paperlike notes. At pseudosteady conditions, all heat getting added to the system goes to the frozen domain, so we don't actually need to treat different areas of the surface. In addition, the sublimation flux is simply evaluated at the top cake surface.

compute_topmassflux(ϕ, T, p, dom::Domain, params)

Compute total mass flow through top of the cake (that is, mass flux integrated across top surface).

compute_Qice_noflow(u, T, dom::Domain, params)

Compute the total heat input into frozen & dried domains, assuming mass flow is zero. Also passes glass-ice heat as separate return.

compute_Qice_nodry(u, T, dom::Domain, params)

Compute the total heat input into frozen domain from volumetric, shelf, and glass. Contrast with compute_Qice_noflow and compute_Qice, which include heat to dried domain.

compute_frontvel_mass(ϕ, Tf, T, p, dom::Domain, params; debug=false)

Generate an empty velocity field and compute velocity on Γ⁺ (i.e. cells on Γ with ϕ>0).

compute_frontvel_heat(u, T, p, dom::Domain, params; debug=false)

Generate an empty velocity field and compute velocity on Γ⁺ (i.e. cells on Γ with ϕ>0).

compute_frontvel_fixedspeed(v0, ϕ, dom::Domain)

Compute speed v0 times normal vector on Γ⁺, positive side of interface.

compute_icesh_area(ϕ, dom::Domain)

Compute area of where ice meets bottom surface of vial.

compute_icegl_area(ϕ, dom::Domain)

Compute area of where ice meets radial outer surface of vial.