Sensitivity of Numerical Predictions to the Permeability Coefficient in Simulations of Melting and Solidification Using the Enthalpy-Porosity Method

The high degree of uncertainty and conflicting literature data on the value of the permeability coefficient (also known as the mushy zone constant), which aims to dampen fluid velocities in the mushy zone and suppress them in solid regions, is a critical drawback when using the fixed-grid enthalpy-porosity technique for modelling non-isothermal phase-change processes. In the present study, the sensitivity of numerical predictions to the value of this coefficient was scrutinised. Using finite-volume based numerical simulations of isothermal and non-isothermal melting and solidification problems, the causes of increased sensitivity were identified. It was found that depending on the mushy-zone thickness and the velocity field, the solid-liquid interface morphology and the rate of phase-change are sensitive to the permeability coefficient. It is demonstrated that numerical predictions of an isothermal phase-change problem are independent of the permeability coefficient for sufficiently fine meshes. It is also shown that sensitivity to the choice of permeability coefficient can be assessed by means of an appropriately defined P\'eclet number.


Introduction
Numerical simulations of melting and solidification processes are critical to develop our understanding of phase transformations that occur in various technologies such as additive manufacturing, thermal energy storage, anti-icing and materials processing. It is however challenging due to the involvement of the moving boundary problem [1] and the wide range of length and time scales in transport phenomena during solid-liquid phase transformations [2]. Different numerical techniques have been developed to resolve the transport phenomena and the solid-liquid phase transition at different scales, which have been reviewed in [3,4,5,6].
Techniques for numerically modelling solid-liquid phase transitions at the continuum level have generally been divided into transformed-grid and fixed-grid approaches [4]. Detailed information on the derivation and implementation of these approaches can be found in the literature [7,8,9]. The focus of the present work is on the fixed-grid approach in which latent heat effects and fluid flow near the liquid-solid interface are taken into account through the inclusion of thermal energy and momentum source terms in that region [3]. One of the advantages of the fixed-grid approach over the transformed-grid approach is its robustness in treating changes in the interface topology. However, interface smearing is an inherent disadvantage, which may be diminished by applying local grid refinement near the solid-liquid interface [10,11,12].
Flow velocities in the liquid phase, particularly in the so-called mushy zone close to the solid interface, where phase-change takes place over a melting temperature range [13], can have a significant influence on the local heat transfer [14]; it is therefore crucial to predict fluid flow in the mushy zone with a sufficient level of accuracy. Of the various approaches that have been employed for this purpose [15,16,17], the porosity approach is the most common. In this approach, it is assumed that the fluid flow in the mushy zone is analogous to that in a porous medium. Consequently, Darcy's law is assumed to govern the flow and a corresponding sink term ( S d ) is added to the momentum equation, where V is the fluid velocity vector and µ dynamic viscosity. Assuming isotropy for the solid-liquid morphology and using the Blake-Kozeny equation [18], the permeability (K) can be defined as a function of the liquid fraction (f L ) as follows: with C the so-called permeability coefficient (or the mushy zone constant). The Blake-Kozeny equation is valid for liquid fractions lower than 0.7 [19,20]; however, it is typically used for the entire range of liquid fractions. The value of the permeability coefficient C and its influence on numerical simulations of melting and solidification is generally uncertain. The values reported for the permeability coefficient C in the literature generally range between 10 3 and 10 15 kg s −1 m −3 [21,22]; however, values between 10 4 and 10 8 kg s −1 m −3 are often applied. The proposed values for C do not seem to have a one-to-one relation to the type of material being studied. For instance, values between O(10 4 ) [23] and O(10 8 ) [24] have been proposed for stainless steel, and values between O(10 4 ) [25] and O(10 15 ) [26] have been proposed for gallium. Thus, the values proposed for C in the literature seem to depend markedly on the process parameters and the associated boundary conditions, and thus lack generality. In other words, tuning the value of C is essential for every set of boundary conditions and material properties, which requires considerable trial and error evaluation. Previous studies on the influence of the permeability coefficient in simulations of melting and solidification have mainly concentrated on finding a value for C to diminish discrepancies between numerical and experimental data for a specific problem (see, for instance, [27,28,29,30,31,32,33]). Additionally, they have often focused on the phasechange materials for thermal energy storage applications for which material properties and operating conditions differ from those for materials processing applications (e.g., casting, welding and additive manufacturing). It is necessary to know under which circumstance and to what extent the numerical predictions of melting and solidification are sensitive to the value of the permeability coefficient to avoid the need for excessive computation. However, there is as yet no general guideline for assessing the appropriate value of C.
The degree of sensitivity of the numerical predictions to the value of the permeability coefficient C appears to be diverse. Several studies reported that the chosen value of C affect the numerical results predicted by the enthalpy-porosity method for isothermal phase-change problems [34,35,36], in which phase transformation occurs at the melting temperature of the material. This effect is not physically realistic for isothermal phase transformations and should be considered as a numerical artefact [37]. Pan et al. [28] studied the melting of calcium chloride hexahydrate (CaCl 2 H 12 O 6 ) in a vertical cylinder and stated that the value of C depends on the temperature difference that drives melting and should be tuned to obtain a reasonable agreement with experimental data. Fadl and Eames [21] reported that the numerical predictions are less sensitive to the value of C in regions where heat transfer is dominated by conduction. Assessing the degree of sensitivity of the numerical prediction to the value of C for different materials and boundary conditions is an open question, despite its significance has been emphasised in the literature [31,38,22,21].
Realising the degree of sensitivity of the numerical predictions to the value of the permeability coefficient in phase-change simulations with a particular interest in melting and solidification during welding and additive manufacturing was the motivation for the present study. Our aim was to quantitatively assess the influence of the permeability coefficient (C) on the results of the enthalpyporosity method in predicting heat and fluid flow and the position of the solid-liquid interface during solid-liquid phase transformations. The sensitivity of the results to the permeability coefficient was analysed for both isothermal and non-isothermal melting and solidification problems, and possible roots of errors in the simulation of solidification and melting processes were highlighted. Our study quantified the influence of the permeability coefficient on the numerical predictions as a function of numerical simulation parameters and physical process parameters and elucidated the limitations of the enthalpy-porosity method.

Problem Description
Isothermal and non-isothermal phase transformations in the two-dimensional rectangular enclosure shown in figure 1 were studied. The length of the enclosure (W) is twice its height (D = 0.1 m). The enclosure is initially filled with solid material at a temperature (T i ) below the melting temperature T m (for isothermal phase-change) or solidus temperature T s (for non-isothermal phase-change). The left and the right solid walls are isothermal walls and the upper and the lower walls are adiabatic. At the starting time of the simulation, the left wall temperature is suddenly raised from the initial temperature (T i ) to the hot wall temperature (T h ), whereas the cold wall is kept at a temperature T c = T i . The thermophysical properties of our artificial materials are presented in table 1, which represent a wide range of metallic and non-metallic phase-change materials. The liquid phase is assumed to be incompressible and Newtonian, with constant dynamic viscosity µ. All other thermophysical material properties are assumed to be the same for both the solid and liquid phases and temperature independent. Figure 1: The schematic of the system under consideration. From t = 0, the left vertical wall is kept at an isothermal temperature T h , and the right vertical wall is kept at an isothermal temperature Tc, whereas the horizontal walls are adiabatic.

Property
Value Unit

Mathematical Formulation
Since the Rayleigh number for the problems considered here is less than 10 6 , the fluid flow is assumed to be laminar. Thermal buoyancy effects are taken into account using the Boussinesq assumption [39]. Utilising the dimensionless variables Fo = tα/D 2 , [10], the unsteady momentum and energy equations are cast in conservative dimensionless form, respectively, as follows: 1 Pr where the Prandtl number (Pr) represents the ratio of momentum diffusivity to thermal diffusivity (α = k/ (ρc p )) in molten materials and is defined as To evaluate the relative importance of buoyancy to viscous forces acting on the molten materials, the Grashof number (Gr) and Rayleigh number (Ra) can be defined as The Stefan number (Ste) is the ratio of sensible to latent heat and is defined as In the above, ρ is density, µ dynamic viscosity, p the static pressure, k thermal conductivity, c p specific heat capacity at constant pressure, f L the local liquid fraction, L m the latent heat of melting or solidification, β the thermal expansion coefficient, g the gravitational acceleration, t time and V the fluid velocity vector.ŷ is the unit vector in the y-axis direction. T melt is the melting-temperature (for isothermal phase-change) or solidus temperature (for non-isothermal phase-change).
For the sake of simplicity, the local liquid fraction was considered to be a function of temperature only, which is a reasonable assumption for the cases where under-cooling is not significant [40]. Different relationships have been proposed for temperature-dependence of the liquid fraction, depending on the materials and the nature of the micro-segregation [41,42]. Based on the method suggested by Voller and Swaminathan [41], the relationship between the local liquid fraction and the temperature is defined to be linear as follows: where T l and T s are liquidus and solidus temperatures, respectively. In the case of isothermal phasetransformation, a step change in the liquid fraction occurs at the melting-temperature according to the method given by Voller and Prakash [16]. Additionally, the convective part of the enthalpy source term in the energy equation (i.e., ∇ V * ) takes the value zero for the isothermal phase transformation due to the step change in the latent heat and a zero velocity at the solid-liquid interface [16].
To deal with the fluid velocity in the mushy zone, a sink term based on the Blake-Kozeny equation (i.e., equation (1) and equation (2)) was introduced into the momentum equation [16].
where C is the permeability coefficient and is a small constant, here chosen to be equal to 10 −3 , to avoid division by zero. The sink term is zero in the liquid region (f L = 1) while its limiting value for f L = 0 should be large enough to dominate the other terms in the momentum equation to suppress the fluid velocities in the solid region. The value of the permeability coefficient C can be expected to depend on the morphology of the mushy zone. No consensus was found on the value of the permeability coefficient in the literature [43,44,45,24,46,47]. The effects of this parameter on numerical results are therefore reported here.
The following dimensionless parameters are utilised to construct a framework for analysing our results. The ratio θ of melting-temperature range (∆T m = T l − T s ) to the temperature difference between the hot and cold wall ( The volumetric fraction of liquid inside the computational domain at a time instant t is evaluated as We use the relative difference ∆φ(t) = |φ 2 (t) − φ 1 (t)| /φ 1 (t) between two solutions obtained using different permeability coefficients C 2 and C 1 to quantify the sensitivity of a solution to the chosen value of the permeability coefficient C. In addition, a dimensionless parameter (Γ) is employed to quantify the difference between the solid-liquid interface morphology predicted using different permeability coefficients C 1 and C 2 as follows: where x f is the position of the liquidus line.

Numerical Procedure
Numerical predictions obtained from an open-source (OpenFOAM) and a commercial (ANSYS Fluent) solver were compared for melting and solidification simulations, and the results were found to be rather identical (see figure 2). ANSYS Fluent (Release 18.1) was selected to carry out the calculations, and simulations were performed in parallel on eight cores of an Intel Xeon E5-2630 processor (2.20 GHz). The computational domain was discretised on a uniform mesh with quadrilateral grid cells. After performing a grid independence test (results presented in sections 5.2 and 5.3.1), a base grid with a cell size ∆x i /D = ∆y i /D = 4 × 10 −3 was chosen. To capture a sharp solid-liquid interface and to reduce errors associated with high-gradient regions, especially when the thickness of the mushy zone is smaller than the base grid size (i.e., (∆T m /∆T w ) W ≤ ∆x i ), a dynamic solution-adaptive mesh refinement was applied using the "gradient approach" [48] based on the liquid fraction gradients. Four levels of mesh refinement ∆x N = ∆x i /2 N , N = 1, · · · , 4 were applied every ten time-steps. A fixed time-step size of ∆Fo = k∆t/ ρc p ∆x 2 i = 10 −4 was selected, corresponding to a Courant number Co = | V |∆t/∆x less than 0.25.
The conservation equations were discretised using the finite-volume method. The centraldifferencing scheme was utilised for the discretisation of the convection and diffusion terms both with second-order accuracy. The Pressure-Implicit with Splitting of Operators (PISO) scheme [49] was used for pressure-velocity coupling. Additionally, the PRESTO (PREssure STaggering Option) scheme [50] was used for the pressure interpolation. The time derivative was discretised with a secondorder implicit scheme. Convergence requires that scaled residuals of the continuity, momentum and the energy equations fall below 10 −10 , 10 −12 and 10 −14 , respectively, and that the relative change in the volumetric fraction of liquid φ from one iteration to the next is less than 10 −10 .

Model Verification
The reliability of the present numerical model was verified against available experimental, theoretical and numerical data for various phase-change benchmark problems. The transient solidification of a semi-infinite slab of Al-4.5%Cu alloy was considered as a one-dimensional problem. The slab is initially in the liquid phase with a uniform temperature of 969 K above the liquidus temperature T l = 919 K. The phase-change process starts at t = 0 s by suddenly changing the temperature of one side of the slab to 573 K, below the solidus T s = 821 K. The results of the present model on a uniform mesh with cell size ∆x = 0.01 m are compared with semi-analytical [51] and numerical results [41] in figure 2a. Our present results are in better agreement with the semi-analytical solution from [51] than the numerical results presented in [41], with the maximum absolute difference between present results and the semi-analytical solution being 0.5%. Melting with convection of pure tin in a square enclosure was considered as a two-dimensional benchmark problem. Details of the problem can be found in [11]. The results obtained from the present model on a uniform mesh with quadrilateral grid cells, ∆x i /W = 4 × 10 −4 for isothermal phase-change including convection are compared with numerical results presented by Hannoun et al. [11] in figure 2b. The maximum difference between the results predicted by the present model and the reference case [11] is within 1%.
Numerical predictions of solid-liquid interface morphologies in a rectangular enclosure subject to natural convection was compared with the experimental observations of Kumar et al. [52] for isothermal melting of lead. Figure 3 indicates a reasonable agreement between the results obtained from the present numerical simulations and the reference data. The deviations between the numerical and experimental data are attributed to both the simplifying assumptions made in the numerical model and the uncertainties associated with the experiments such as the uncertainty in determining the position of solid-liquid interface, which in this case is reported to be about 8.7% [52]. The results obtained from the present model are also validated with experimental observations in section 6 for a laser spot melting process (see figure 14).

Grid Size and Sensitivity to the Permeability Coefficient for Isothermal Phase Change
Sensitivity of the solution to the grid size and the value of the permeability coefficient was studied for isothermal phase-change of gallium melting at T m = 302.78 K, in a side-heated rectangular enclosure with T i = T c = 301.30 K, and T h = 311.0 K, which was experimentally investigated by Gau and Viskanta [53]. This benchmark case has often been considered for validation of phasechange simulations in the literature [4]. Detailed information regarding the benchmark case is available in [53,54]. Figure 4 shows the influence of the permeability coefficient on the predicted liquid fraction, using uniform fixed meshes with 42 × 32, 105 × 80, 210 × 160 and 420 × 320 computational grid cells and C = 10 4 and 10 8 kg s −1 m −3 . When a coarse mesh is utilised, resulting in a non-physical mushy zone of significant thickness, the solution appears to be very sensitive to the value of the permeability coefficient. This sensitivity is attributed to the enhanced heat transfer from the hot fluid to the solidliquid interface due to higher fluid velocities in the mushy zone predicted with a smaller permeability coefficient. Reducing the grid cell size, and thus reducing the thickness of the non-physical mushy zone, decreases the total volume of the grid cell in which the permeability coefficient affects the numerical predictions. Consequently, the sensitivity of the solution to the value of the permeability coefficient decreases with grid refinement. A finite amount of time is required for the mass in a computational cell to absorb heat and melt. A change in the value of C therefore affects the convective heat transfer to the cells located at the solid-liquid interface that can lead to a change in the predicted interface morphology and the rate of melting during the transient phase. This effect reduces by refining the grid size adjacent to the solid-liquid interface [10]. In figure 5, the sensitivity to the value of C is quantified as a function of grid cell size by looking at the parameters ∆φ and Γ (defined below equation (12) and in equation (13), representing the relative difference in the predicted liquid fractions, and the relative difference between the solid-liquid interface morphologies) for C 1 = 1 × 10 8 kg s −1 m −3 and C 2 = 1 × 10 4 kg s −1 m −3 at t = 150 s. Reducing the cell size decreases the influence of the permeability coefficient on the results, with Γ scaling approximately linearly with cell size, and ∆φ scaling approximately quadratically with cell size. Thus, although in principle the chosen value of C should be irrelevant for isothermal phase change, sufficient grid refinement is needed to obtain solutions which are indeed insensitive to the value of C. In addition, sufficient grid resolution in the liquid zone was found to be required to predict the multicellular convection [55,56,14,57] early in the phase-change process, which leads to the formation of a wavy solid-liquid interface. According to equation (10), the magnitude of the sink-term in solid regions (f L = 0.0) should be large compared to the viscous term (∇ 2 V * ) to suppress fluid velocities (i.e., C D 2 / (µ ) 1). The influence of the sink-term magnitude on predicted liquid fraction φ and the ratio of the maximum velocity magnitude in solid regions to that in liquid regions (i.e., | V solid | Max / | V liquid | Max ) is shown in figure 6. Here, D is chosen to be the size of the heated wall. For magnitudes of C D 2 / (µ ) roughly larger than 2.5 × 10 7 , the predicted liquid fraction is independent of the sink-term. For smaller values of the sink-term, even though the velocity magnitude in the solid region is orders of magnitude smaller than that in the liquid region, the energy transfer to the solid material is affected, which leads to a change in the predicted liquid fraction. The fine-grid results presented in figure 4 represent a mathematically converged solution of the governing equations for the given set of boundary conditions. However, the grid independent results of this particular case do not closely match with the experimental results reported in [53]. Obtaining a better agreement between numerical predictions and experimental observations by using a coarser grid (i.e., grid-dependent results), or by tuning the mathematical model has been reported in the literature [34] that is coincidental, but such ad hoc tuning lacks generality [54,11,58].  (13)), circles in blue), and the relative difference between the liquid fractions (∆φ (defined below equation (12)

Grid Sensitivity
Having seen the importance of sufficient grid resolution for isothermal phase-change, we considered the non-isothermal phase-change problem defined in section 2, examining the impact of grid refinement on the sensitivity to the permeability coefficient C. For ∆T m = 50 K, figure 7 shows the sensitivity parameter ∆φ, for short time (Fo = 0.12) and steady-state (Fo = 9.0) obtained with C 1 = 1 × 10 8 kg s −1 m −3 and C 2 = 1 × 10 4 kg s −1 m −3 . For comparison, the same problem was also solved assuming isothermal phase change (∆T m = 0 K). For isothermal phase change, the results are insensitive to the chosen value of C in the steady-state, and exhibit a small sensitivity to the chosen value of C during the transient phase before reaching the steady-state that becomes negligible with reducing grid size. For non-isothermal phase change, the sensitivity to the chosen value of C is much stronger and does not vanish with decreasing grid size, even in the steady-state, as is to be expected from the fact that a physically realistic mushy zone is now present, in which flow and convective heat transfer are sensitive to the value of C. However, the sensitivity to C becomes grid independent for base grid sizes below ∆x i /D = 4 × 10 −3 . Consequently, this base grid size is used in the next sections. In addition, a four-level dynamic solution-adaptive mesh refinement (i.e., ∆x N = ∆x i /2 N , N = 1, · · · , 4) is applied to further enhance the accuracy with which we capture the solid-liquid interface.

Influence of the Permeability Coefficient on Predicted Results
For the non-isothermal phase change problem defined in section 2, the influence of the permeability coefficient on the predicted steady-state liquidus-line position (f L = 1) is shown in figure 8 for a fixed wall temperature difference (∆T w ), a fixed melting temperature range (∆T m ) and varying Rayleigh number (i.e., varying fluid velocities). In the virtual absence of flow, for Ra = 1, the results are insensitive to the value of the permeability coefficient. In this case, heat conduction dominates the total energy transfer around and in the mushy zone. By increasing the value of Ra, i.e. increasing fluid flow velocities and increasing convective heat transfer, the results become more sensitive to the chosen value of C. This can be understood from the fact that the value of C, through the momentum sink term, only influences the convective terms and therefore the results become more sensitive to C when convection plays an important role in total heat transfer in the mushy zone.  Figure 9 shows the effect of the permeability coefficient for a fixed melting temperature range (∆T m ) with different wall-temperature differences (∆T w ). Here, a higher value of ∆T w (i.e., smaller θ) leads to a smaller mushy zone thickness, and consequently a reduced sensitivity to the permeability coefficient. Similarly, for a fixed (∆T w ), reducing (∆T m ) decreases the mushy zone thickness and as a result lower the sensitivity to the permeability coefficient C, as shown in figure 10. In summary, the results show that non-isothermal phase-change simulations are less affected by the chosen value of the permeability coefficient C when the mushy zone has a smaller thickness (i.e., smaller θ), as this leads to the conductive heat transfer through the mushy zone being large compared to the convective heat transfer. A change in the thermal diffusivity of the material can therefore affect the sensitivity of the numerical predictions to the value of the permeability coefficient. Figure 11 indicates the sensitivity of the results to the chosen value of the permeability coefficient as a function of thermal diffusivity of the material for a fixed melting temperature range (∆T m ) and wall-temperature difference (∆T w ). It is seen that sensitivity to the permeability coefficient decreases with increasing thermal diffusivity of the material, which can be attributed to the enhancement of heat conduction contribution to the total heat transfer.   The influence of the permeability coefficient on the time evolution of the melt pool shape is presented in figure 12. The results are indeed independent of the permeability coefficient for an isothermal phase-change. A similar conclusion has also been drawn after monitoring the time variations of the liquid fraction obtained from numerical simulations [59,52]. However, for nonisothermal phase change with a thick mushy zone (θ = 1/15) and strong flow (Ra = 10 6 ), the results are very sensitive to the chosen value of the permeability coefficient, resulting is different pool shapes and rates of phase-change. Less sensitivity to the value of C is found for a thin mushy zone (θ = 1/30). The numerical predictions are also less sensitive to the permeability coefficient at early time instances, when the fluid flow is characterised by low velocities.

Discussion
The results reveal that, depending on the temperature gradient, velocity field and thermophysical properties of the phase-change material, which in turn determine the mushy zone thickness, numerical predictions of phase-change problems can show sensitivity to the value of the permeability coefficient. A general guideline is presented here that allows prediction and evaluation of the influence of the permeability coefficient in phase-change simulations. This concept involves both the heat transfer mechanism and the mushy zone thickness.
For phase-change problems without fluid flow, the mushy zone thickness (m t ) can be estimated as where L is a characteristic length scale. The fluid flow can alter the mushy zone thickness when convective heat transfer in the mushy zone is of significance. This can be expressed through a Péclet number, which expresses the ratio between the rate of heat advection and heat diffusion.
To identify the regions where there is a significant heat transfer enhancement or reduction due to convection, and consequently where the results are sensitive to the permeability coefficient, Pe * is defined as follows: Non-zero values of Pe * adjacent to the solid-liquid interface indicate increased sensitivity to the permeability coefficient. For isothermal phase-change, Pe * is zero and predictions are independent of the permeability coefficient.
Sensitivity of the numerical predictions to the value of permeability coefficient C has been appraised for the problem defined in section 2, using three different melting temperature ranges ∆T m and a fixed ∆T w = 750 K, and the results are shown in figure 13. Higher values of Pe * (≈ O(10)) along the melting front for the case with ∆T m = 55 K (figure 13c) compared to the case with ∆T m = 10 K (figure 13b) indicate that the numerical predictions are more sensitive to the permeability coefficient, which is consistent with the results presented in figures 10 and 12. When the values of Pe * along the interface are small (i.e., Pe * 1), the results appear to be insensitive to the permeability coefficient, while, for Pe * 1, the results are sensitive to C. The general applicability of the proposed Pe * criterion is also examined for a simulation of a laser spot melting process, which was carried out experimentally by Pitschender et al. [60]. A steel plate containing 20 ppm of sulphur was heated using a stationary laser beam with a power of 5200 W and a top-hat radius of 1.4 mm. The surface absorptivity is set to 0.13 [60]. Material properties of the steel plate are assumed to be constant and temperature independent, except for the surface tension of the molten material, and can be found in [60,61,62]. Variations of surface tension with temperature are modelled using the expression proposed by Sahoo et al. [63]. Melting of the material and associated heat and fluid flow are numerically predicted using permeability coefficients of C = 10 6 , 10 8 and 10 10 kg s −1 m −3 . To study the influence of the mushy-zone thickness on the sensitivity to C through Pe * , the melting-temperature range (∆T m ) of the material is changed artificially to 40 K and 200 K. Figure 14 shows the position of melting front (f L = 1) at t = 5 s, as well as the value of Pe * along it when using different permeability coefficients and melting temperature ranges. The laser-beam diameter is chosen here as the characteristic length scale L to calculate Pe * . Higher values of Pe * are found for larger ∆T m , predicting more sensitivity to the permeability coefficient as is indeed observed when comparing figure 14a and figure 14b. Due to a steep increase in the momentum sink term with liquid fraction for large permeability coefficients, all predictions are found to converge to an identical solution for very large C. However, a very large value of the permeability can lead to numerical instabilities.

Conclusions
A systematic study was performed to scrutinise the influence of the permeability coefficient on the numerical predictions of isothermal and non-isothermal phase-change simulations using the enthalpyporosity method.
For isothermal phase-change problems, reducing the cell size diminishes the influence of the permeability coefficient on the results, which become independent of the permeability coefficient for fine enough meshes. A grid independent solution of an isothermal phase-change problem is independent of the permeability coefficient. However, not every numerical result that is independent of the permeability coefficient is grid independent.
Numerical predictions of non-isothermal phase-change problems are inherently dependent on the permeability coefficient. The sensitivity of the numerical predictions to the permeability coefficient increases with increased mushy zone thickness and increased fluid flow velocities perpendicular to the solid-liquid interface. A method is proposed to predict and evaluate the influence of the permeability coefficient on numerical predictions, and verified for two-dimensional phase-change problems including laser spot melting. Large values of Pe * 1 adjacent to the solid-liquid interface indicate a strong sensitivity and Pe * 1 indicates insensitivity to the permeability coefficient C.