Struct SolidificationMelting 

pub struct SolidificationMelting { /* private fields */ }

Solidification and melting of a pure or eutectic material.

The enthalpy-porosity method

The difficulty with melting is that the solid/liquid interface moves, and tracking it explicitly on a fixed mesh is painful. The enthalpy-porosity method avoids tracking it at all. Every cell carries a liquid fraction α₁ between 0 and 1, and two things follow from it:

• Latent heat enters the energy equation as a source proportional to ∂(ρα₁)/∂t. A cell that is melting absorbs energy without rising in temperature, which is exactly what latent heat means.
• The solid is immobilised by treating the mushy zone as a porous medium: a Darcy drag -Cu(1-α₁)²/(α₁³+q) is applied to the momentum equation. In a fully liquid cell α₁ = 1 and the drag vanishes; as α₁ falls the drag grows without bound, pinning the velocity to zero. No moving boundary is ever needed.

The Carman-Kozeny form of that drag is not arbitrary — it is the permeability of a porous bed of fixed solid fraction, which is what a mushy zone physically is.

Buoyancy is included in the Boussinesq form, ρ_ref g β ΔT, since natural convection in the melt is usually what drives the process.

State

Unlike most of this crate, the model is stateful: it stores the liquid fraction per cell and advances it each time it is applied. Upstream hides this behind C++ mutable; here the apply methods take &mut self, which makes the state visible in the signature rather than hidden.

Not covered

Upstream also handles a thermophysical-model lookup for Cp and mesh topology changes. This port takes Cp as a coefficient (upstream’s CpRef mode) and does not track topology changes, since the mesh here is fixed at construction.

Implementations

impl SolidificationMelting

pub fn new( name: impl Into<String>, velocity_name: impl Into<String>, energy_name: impl Into<String>, energy_is_temperature: bool, selection: CellSelection, coefficients: SolidificationMeltingCoefficients, gravity: Vector3, n_cells: usize, ) -> Self

Build the model.

energy_is_temperature selects which form the energy equation is solved in: true for temperature, false for enthalpy. Upstream makes the same distinction by inspecting the solved field’s dimensions, and it matters — the latent-heat source is divided by Cp in the temperature form and not in the enthalpy form. Getting it wrong scales the latent heat by Cp, typically a factor of several hundred.

n_cells is the mesh cell count, needed to size the internal state.

pub fn name(&self) -> &str

The model’s name.

pub fn velocity_name(&self) -> &str

The velocity field this adds a momentum sink to.

pub fn energy_name(&self) -> &str

The energy field this adds latent heat to.

pub fn selection(&self) -> &CellSelection

The cells acted on.

pub fn liquid_fraction(&self) -> &[f64]

The current liquid fraction, indexed by position within the selection (not by mesh cell index).

pub fn effective_liquidus(&self, liquid_fraction: f64) -> f64

The effective liquidus temperature [K] at a given liquid fraction.

max(T_sol, T_sol + (T_liq - T_sol)(α₁ - α₁ₑ)/(1 - α₁ₑ)) — upstream’s expression. The max matters: below the eutectic fraction the linear term would fall below the solidus, and clamping there is what keeps the eutectic plateau flat.

pub fn update(&mut self, temperature: &VolScalarField)

Advance the liquid fraction from the current temperature field, at most once per timestep.

Upstream’s update(). Called automatically by the apply methods; public because a caller driving the model directly may want it.

The update is explicit and under-relaxed:

α₁ ← clamp(α₁ + relax · Cp · (T - T_liq_eff(α₁)) / L, 0, 1)

then ΔT is recomputed against the updated fraction — which is what upstream does, and is not the same as using the old one.

The once-per-timestep guard

A solver applies this model to two equations per timestep — momentum and energy — and each apply needs a current liquid fraction. Letting the update run on both would advance the under-relaxed explicit iteration twice per step, so the effective relaxation would be a function of how many equations happen to reference the model. Upstream guards against that with curTimeIndex_; this port uses a flag cleared by advance_time. Calls after the first in a step return without doing anything.

pub fn advance_time(&mut self)

Roll the liquid-fraction history forward and re-arm update.

Call once per completed timestep. It does two things, both required:

• The latent-heat source is proportional to ∂(ρα₁)/∂t, so it needs the previous step’s fraction; without this call that rate stays measured against an ever-staler state and the latent heat is silently wrong.
• It clears the once-per-timestep guard documented on update. Omitting the call therefore freezes the liquid fraction entirely, rather than merely staling it.

pub fn darcy_coefficient(&self, liquid_fraction: f64) -> f64

The Darcy drag coefficient at a given liquid fraction [kg/(m³·s)].

-Cu (1-α₁)² / (α₁³ + q). Negative, since it is a sink.

pub fn add_momentum_source( &mut self, temperature: &VolScalarField, eqn: &mut FvVectorMatrix, )

Add the momentum sink and buoyancy to a velocity equation.

Upstream’s addSup(U, eqn). The Darcy drag is the implicit coefficient and goes onto the diagonal (stabilising); the Boussinesq buoyancy is explicit and goes into the source. Both are scaled by the cell volume.

Why this is not a literal transcription of upstream

Upstream writes

Sp[celli] += Vc*S;                      // S  = -Cu(1-a)^2/(a^3+q)
Su[celli] += Vc*Sb;                     // Sb = rhoRef*g*beta*deltaT

directly into eqn.diag() and eqn.source() — but that eqn is not the solver’s equation. It is a separate fvMatrix that fvModels::sourceTerm builds and that the solver then combines with solve(UEqn == fvModels.source(rho, U)), where operator==(A, B) expands to A - B. Every coefficient an addSup writes is therefore negated before it reaches the system actually solved.

Transcribing the two lines literally into a matrix that is the solved system inverts both terms: the Darcy drag becomes a source of momentum inside the solid instead of a sink, and it lands on the diagonal with the wrong sign, destroying diagonal dominance exactly where the coefficient is largest. This port therefore places the terms per the convention documented on SourceContribution — source += V·explicit, diag -= V·implicit — which reproduces upstream’s solved equation rather than its intermediate matrix.

The resulting explicit buoyancy is -ρ_ref β ΔT g, which is the standard Boussinesq force: a cell hotter than its effective liquidus (ΔT > 0) is pushed against gravity.

pub fn add_energy_source( &mut self, rho: &VolScalarField, temperature: &VolScalarField, dt: f64, eqn: &mut FvMatrix, )

Add the latent-heat source to an energy equation.

Upstream’s addSup(rho, he, eqn) via apply(), which writes eqn -= L·∂(ρα₁)/∂t (or eqn -= (L/Cp)·∂(ρα₁)/∂t in the temperature form) into the intermediate fvModels matrix. Following that through the solve(hEqn == fvModels.source(...)) negation described on add_momentum_source, the term reaching the solved system’s right-hand side is

-L · ∂(ρα₁)/∂t, or -(L/Cp) · ∂(ρα₁)/∂t when the equation is solved in temperature rather than enthalpy.

The time derivative is discretised as first-order Euler against the previous step’s liquid fraction, so advance_time must be called once per completed timestep.

The sign is a sink: melting absorbs energy. That is why the contribution is subtracted from the equation, and it is the term that holds a melting cell at its melting point instead of letting it heat through.

Trait Implementations

impl Clone for SolidificationMelting

fn clone(&self) -> Self

Returns a duplicate of the value.

fn clone_from(&mut self, source: &Self)

Performs copy-assignment from source.

impl Debug for SolidificationMelting

fn fmt(&self, f: &mut Formatter<'_>) -> Result

Formats the value using the given formatter.

impl PartialEq for SolidificationMelting

fn eq(&self, other: &Self) -> bool

Equality operator ==.

fn ne(&self, other: &Rhs) -> bool

Inequality operator !=.

impl StructuralPartialEq for SolidificationMelting