Model for a Freezing Food Slab

Initializing live version
Download to Desktop

Requires a Wolfram Notebook System

Interact on desktop, mobile and cloud with the free Wolfram Player or other Wolfram Language products.

This project numerically solves for the moving boundary (MB) displacement during freezing (thawing) of a one-dimensional moist-food slab. The analysis uses the pseudo-steady state assumption (PSS), convective cooling at the exposed surface, and a temperature-dependent thermal conductivity. We first compute the surface temperature assuming a MB (s) vector using the iterative Nest function. By knowing the surface temperature, the energy differential equation can be integrated over time, which satisfies the assumed s vector. This integration takes into account the thermal properties of the slab and the cooling conditions as input parameters. Subsequently, the surface heat flux, temperature profile, and the influence of perturbations in the normalized thermal conductivity parameter (Κ) on the MB location are calculated.This time-dependent perturbation effect is graphically presented as a dimensionless sensitivity σ. Our predictions are compared against the well-known Plank freezing solution. Noteworthy observations include the influence of the thermal conductivity parameter (Κ), the dimensionless latent heat number (λ), the Biot number (Bi), and the dimensionless cooling-medium temperature on the MB position and its sensitivity, which is driven by Κ. This Demonstration thus allows the estimation of the extent to which deviations in thermal conductivity affect the MB location during freezing (thawing).

Contributed by: Victor M. Chavarria (March 2020)
Open content licensed under CC BY-NC-SA


The moving boundary (MB) is the idealized infinitesimally-thin interface separating the unfrozen region from the freezing region. To locate it, and assuming the pseudo steady state (PSS) condition, requires solving the heat conduction differential Eq. (1) below that describes the freezing process of a one-dimensional slab [1]:

, (1)

where is the dimensionless slab temperature defined as:


and are the (initial) freezing points of pure water and of the food material, respectively,

is the dimensionless freezing front (length) growing from to ,

is the normalized distance (spanning from the convective surface to the slab center),

and the dimensionless latent heat number defined λ = Δ​H/(C pf·(T0-Ti));

where is the heat of fusion of the moist food material adjusted by its unfreezable water content and is the volumetric specific heat capacity of the fully frozen slab [1].

The convective boundary condition (BC) is:

, (2)

where is the convective-surface slab temperature at . We adopt Schwartzberg's [1] temperature-dependent thermal conductivity model given by:

, (3)

where is the thermal conductivity parameter given by ,

and and are the unfrozen and fully frozen state thermal conductivities, respectively. governs the temperature response of thermal conductivity and influences the freezing process in a nonlinear fashion.

At the moving interface, the BC gives: .

Other relevant key assumptions include the full release of the latent heat of fusion, the material temperature at time 0 and the phase change all take place at the initial freezing point of the material, and no sensible heat is released.

Equations (1) and (2) together with (3) were numerically solved, first using Mathematica's iterative Nest function for in terms of a discretized List. We then back-calculate the time steps from the known increments and s List applying the simple forward Euler method. The slab temperature profile was solved applying the Nest function to Eq. (1), subject to the known temperature conditions at the surface and at the interface (s). As the MB penetrates the slab thickness, the (time-dependent) temperature profile develops in the MB penetration region only, while the unfrozen slab region is kept at the initial food freezing point (Ti).

Plank's solution serves as a simple but well known reference to the freezing slab problem, Eq. (4) below, even though it assumes constant thermal properties [2] while the other assumptions are identical to our assumptions. This solution still provides a simple but useful approximations to this phase change problem:

s(τ) = + , (4)

where St is the Stefan number defined as St = .

The MB sensitivity to K perturbations is calculated with small but finite changes (=3.34%) in the thermal conductivity parameter Κ and is given by Eq. (5):

= s(τ)/s(τ))/(δ , (5)

This sensitivity quantity represents a time-dependent response of the MB position to errors or perturbations in Κ; and thus provides a quantifiable measure of the influence of this relevant material thermal property on the interface movement. A sensitivity of σ = ±0.25 at a given point in time indicates that a ±10% perturbation or error in Κ, for instance, generates a ±2.5% variation in the MB position.

Results are plotted against dimensionless time , where ; here is the slab half thickness and is dimensional time.

The interesting effects that can be controlled as parameters in the plots are:

1. thermal conductivity parameter

2. Biot number , where is the convective heat transfer coefficient,

3. dimensionless latent heat number

4. temperature driving force , based on , the medium-cooling temperature

5. slab temperature profile , which follows a close to but not perfectly linear pattern.


[1] H. G. Schwartzberg, "The Prediction of Freezing and Thawing Temperatures vs. Time through the Use of Effective Heat Capacity Equations," in Proceedings of the Joint Meeting of Commissions C1 and C2 of the International Institute of Refrigeration, 1977 pp. 311–317.

[2] V. M. Chavarria, "Modeling the Influence of Temperature-Dependent Thermal Properties on the Freezing Front," Journal of Food Research, 8(6), 2019 pp. 129–146. doi:10.5539/jfr.v8n6p129.


Feedback (field required)
Email (field required) Name
Occupation Organization
Note: Your message & contact information may be shared with the author of any specific Demonstration for which you give feedback.