2.1.1 Mathematical and numerical background

The heat transfer equation for the 1D model follows that given by the BFFTM (Mooney et al., 2012) but without the advection term because no translation of the sample takes place. This equation is expressed as:

where t is the time, ρ is the density of the cylindrical bar, c_p is the specific heat, T is the temperature, k is the thermal conductivity, h_I is the interfacial heat transfer coefficient, C is the circumference of the cylindrical bar, A is the cross-sectional area, T_Wall is the constant temperature of the crucible, and E is the latent heat term. An expansion of the latent heat term gives

where L is the specific latent heat and f_S is the solid fraction.

Figure 2Discretisation of the cylindrical rod with disc shaped control volumes.

Download

The relationship between the solid fraction and temperature, T, is given as follows:

In the FVM, the geometry of the cylindrical rod is discretised with a finite number of the control volume (CVs) discs having a width of Δx. Figure 2 depicts the details of the disc shaped CV, labelled as “i” and its adjacent CVs are labelled as “i−1” and “i+1”. The solution to Eq. (1) can be implemented by considering the energy balance for CV “i” at time step “n”:

where $T_{i,n+\mathrm{1}}$ and T_i,n, are the temperatures at the next time step “n+1” and current time step “n” of CV “i”, r is the radius of the bar, f_i,n is the solid fraction, and $q_{\mathrm{rad},i,n}$ is the radial heat loss from the side of the disc-shape CV, as indicated by Fig. 2. In this study, this is simplified as equivalent to the interfacial heat between the bar and the crucible, which can be expressed as:

The term $q_{i-\mathrm{1},n}$ in Eq. (4) is the interfacial heat flux of thermal diffusion from the previous CV “i−1” to CV “i”, and $q_{i+\mathrm{1},n}$ is the interfacial heat flux from CV “i” to the next CV “i+1”, which can be expressed as the following finite difference equations:

where $T_{i-\mathrm{1},n}$, T_i,n and $T_{i+\mathrm{1},n}$ are the temperatures of CVs “i−1”, “i”, and “i+1” respectively at time step “n”. Inserting Eqs. (5), (6), and (7) into Eq. (4) and rearranging gives Eq. (8), which can be used to compute the temperature distribution of CV “i” at the next time step “n+1”:

where T_L is the liquidus temperature and T_S is the solidus temperature. This approach is known as an explicit scheme.

Latent heat “L” is only considered in the mushy zone, i.e., when the temperature of CV is greater than the solidus temperature and lower than the liquidus temperature ($T_{\mathrm{L}}<T_{i,n}<T_{\mathrm{S}}$) as indicated by Eq. (3).

Figure 3a and b shows a flowchart for the main code algorithm and subroutine which has been used in the production of the non-commercial code.

Figure 3(a) Flow chart describing the relevant steps taken within the non-commercial code. (b) Flow chart describing the relevant steps taken within subroutine of the non-commercial code.

Download

2.2.1 Modelling phase change in COMSOL Multiphysics^®

Phase change in COMSOL Multiphysics^® was performed using an apparent heat capacity method (COMSOL Multiphysics ^® v.5.4, 2018). The inbuilt COMSOL^® phase change material sub-node was used to specify the properties of the material to change from phase 1 to phase 2. In our case phase 1 was designated as solid and phase 2 was designated as liquid. In COMSOL^®, the transition between the phases is assumed to occur smoothly over a temperature interval ΔT. COMSOL^® accounts for the release of latent heat over the temperature interval between $T_{\mathrm{pc}}-\mathrm{\Delta }T/$2 and $T_{\mathrm{pc}}+\mathrm{\Delta }T/\mathrm{2}$ where T_pc is the temperature which initiates the phase change in material. During this interval, the phenomenon of phase change within the material is assigned to thermal models using a smoothed function, θ, to represent the solid fraction. Equations (9) and (10) represents the fraction of phase prior to transition.

Material properties such as density, ρ, and enthalpy, H, were updated according to the phase change taking place in that temperature interval. Equations (11) and (12) were used to accommodate this phenomenon in the COMSOL Multiphysics^® model:

where, ph1 and ph2 indicate phase 1 and phase 2 of the material. Input parameters for the transition are given in Table 1. The COMSOL^® phase change input parameters are related to the materials liquidus and solidus temperatures, T_L and T_S respectively.

Table 1Parameters and values used within the COMSOL Multiphysics^® subnode.

Download Print Version | Download XLSX

2.2.2 The heat input

The 3D model accounted for a time-dependent power input to the left-hand boundary using a piecewise linear function. The function was set to 0 W at 0 s with full ramping up of the power from the heat source after 1 s. Thereafter, the power maintained its maximum value, 400 W, for the next 34 s. After 35 s of process time, the power ramped down to 0 W over 5 s. The simulation ran for 100 s (140 s of process time). A general stationary inward heat flux, q₀, has been considered in the present model. This adds to the total flux across selected boundaries.

The inward heat flux, q₀, may be calculated from the power input as:

where, q₀ is the general inward heat flux, P is the power of heat source from the piecewise linear function described earlier and A is the cross-sectional area of the cylindrical rod. A peak q₀ value as 5.093×10⁶ W m⁻² relates to a maximum power input of 400 W over a 10 mm diameter cross-section.